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