function [ sync_signal, rotation_matrix ] = brain_sync( reference_signal, moving_signal, mean_0_norm_1 )
%% BRAIN_SYNC: AA Joshi, M Chong, RM Leahy, 2017
%
% Syncronize the moving signal to the reference signal, using BrainSync.
%
% Input:
% reference_signal: Time series as the reference of BrainSync (vertices
% x time)
% moving_signal: Time series (individual, e.g.) to be rotated to sync
% with reference (vertices x time)
% mean_0_norm_1: Flag to indicate if the signal is already
% pre-processed to be zero mean and unit norm. (binary, default 0)
% Output:
% sync_signal: Syncronized time series, with respect to reference
% signal (vertices x time)
% rotation_matrix: The rotation matrix that syncronize the two time
% series (time x time)
%
%
% Please cite the following publication:
% AA Joshi, M Chong, RM Leahy, BrainSync: An Orthogonal Transformation for Synchronization of fMRI Data Across Subjects,, Proc. MICCAI 2017, in press.
%
%% Pre-Processing
if ~(exist('mean_0_norm_1', 'var'))
mean_0_norm_1 = 0 ;
end
if ~(mean_0_norm_1)
[reference_signal, reference_mean, reference_std] = brain_sync_pre(reference_signal) ;
[moving_signal, ~, ~] = brain_sync_pre(moving_signal) ;
end

%% BrainSync
cross_corr = moving_signal' * reference_signal ;
[U, ~, V] = svd(cross_corr) ;
rotation_matrix = V * U' ;
sync_signal = moving_signal * rotation_matrix' ;

%%
% Map back if necessary
if ~(mean_0_norm_1)
sync_signal = brain_sync_back(sync_signal, reference_mean, reference_std) ;
end

end

%% Auxillary functions
% To normalize signal to zero mean and unit norm
function [normed_signal, mean_vector, std_vector] = brain_sync_pre(pre_signal)
ones_vector = ones(1, size(pre_signal, 2)) ;
pre_signal(isnan(pre_signal)) = 0 ;
mean_vector = mean(pre_signal, 2) ;
normed_signal = pre_signal - mean_vector * ones_vector ;
std_vector = std(normed_signal, 1, 2) ;
std_vector(std_vector == 0) = 1 ;
normed_signal = normed_signal ./ (std_vector * ones_vector) ;
end

% Map the signal back to match the reference signal norm and mean
function [back_signal] = brain_sync_back(normed_signal, mean_vector, std_vector)
ones_vector = ones(1, size(normed_signal, 2)) ;
back_signal = normed_signal .* (std_vector * ones_vector) ;
back_signal = back_signal + mean_vector * ones_vector ;
end