function [transf,varargout] = mri2scs(MRI,varargin);
%MRI2SCS - Compute the transform to move form the MRI coordinate system (in mm) to the SCS 
% function [transf,varargout] = mri2scs(MRI,varargin);
% function [transf] = mri2scs(MRI);
%  MRI - a proper BrainStorm MRI structure (i.e. from any subjectimage file, with fiducial points 
%        and SCS system properly defined)
%  transf - a structure specifying the transform that is applied to the MRI mm coordinates
%  Definition of the transform is the following:
%  Xscs = transf.R Xmri + transf.T ; 
%  (Xmri in mm)
%  transf.Origin is the location of the SCS origin in MRI coordinates
%  
%  function [transf,scsCoord] = mri2scs(MRI,mriCoord);
%  mriCoord - a 3xN matrix of point coordinates in the MRI system (in mm)
%  scsCoord - a 3xN matric of corresponding point coordinates in the SCS system (in mm)

%<autobegin> ---------------------- 08-Jun-2004 15:11:31 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Utility - Numeric
%
% At Check-in: $Author: Mosher $  $Revision: 12 $  $Date: 6/08/04 1:38p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 28-May-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 08-Jun-2004 15:11:31 -----------------------

  
 
% /---Script Authors-------------------------------------\
% |  *** Sylvain Baillet, Ph.D.                          |
% |  Cognitive Neuroscience & Brain Imaging Laboratory   |
% |  CNRS UPR640 - LENA                                  | 
% |  Hopital de la Salpetriere, Paris, France            |
% |  sylvain.baillet@chups.jussieu.fr                    |
% |                                                      |
% |  *** Alexei Ossadtchi, Ph.D.                         |
% |  USC - SIPI                                          |
% |                                                      |
% \------------------------------------------------------/
%  
% Script History ---------------------------------------------------------------
% SB  15-Mar-2004 (BrainStorm v.1 update)
% Script History ---------------------------------------------------------------



% Memo
% FidNames0{1} = 'Nasion';
% FidNames0{2} = 'LEFT Preauricular';
% FidNames0{3} = 'RIGHT Preauricular';

if ~isfield(MRI,'SCS')
    errordlg({'MRI SCS was not defined or subjectimage file is from another version of BrainStorm',...
            'Please re-define the subject coordinate system on this MRI'},'MRI file error')
    return
end

% Fiducial coordinates in mm
% Nasion
NAS = MRI.SCS.mmCubeFiducial(:,1);%MRI.SCS.LookUpFiducials(1));
% MRI.SCS.Label{LookUpFiducials(1)} = 'Nasion';
% LPA 
LPA = MRI.SCS.mmCubeFiducial(:,2);%MRI.SCS.LookUpFiducials(2));
% MRI.SCS.Label{LookUpFiducials(2)} = 'LEFT Preauricular';
% RPA
RPA = MRI.SCS.mmCubeFiducial(:,3);%MRI.SCS.LookUpFiducials(3));
% MRI.SCS.Label{LookUpFiducials(3)} = 'RIGHT Preauricular';

% Move to SCS 
% Definition:
% Xscs = R Xmri + T ; Xmri in mm

switch(MRI.SCS.System)
    
case 'CTF'
    transf.Origin = .5*(LPA+RPA); % in mm
	%     nLR = cross(NAS-RPA,LPA-RPA);  % Vector normal to the NLR plane, pointing upwards
	%     nLR = nLR/norm(nLR);
    
    vx = (NAS-transf.Origin)/norm((NAS-transf.Origin));
	%     vy = inv([vx,nLR,cross(vx,nLR)])*[0 0 1]'; % y direction
	%     vy = vy/norm(vy);
    
	vz = cross(vx, LPA-RPA); vz = vz/norm(vz);  % Vector normal to the NLR plane, pointing upwards
	vy = cross(vz,vx); vy = vy/norm(vy);
	
	transf.R = inv([vx,vy,vz]); % Rotation
	transf.T = - transf.R * transf.Origin;
   
    %rot = [0 1 0;1 0 0;0 0 1]*inv([vx,vy,cross(vx,vy)]);
    %rot = [0 1 0;-1 0 0;0 0 1]* rot;
    
case 'NEUROMAG' % Warning - need to be re-checked !

    
    LR = LPA - RPA; % Left\Right Vector
    LR = LR/norm(LR);  
    LRp = [0 -1 0;1 0 0;0 0 0]* LR; % A vector orthogonal to LR that will define a plane perpendicular to the NLR plane
    LRp = LRp/norm(LRp);
    % note: The LR line is defined according to the intersection between 2 orthogonal planes.
    
    OR = dot(RPA-NAS,RPA-LPA)/dot(RPA-LPA,RPA-LPA) * (RPA-LPA);
    Origin = RPA - OR;
    vx = OR/norm(OR);
    vy = (NAS-Origin)/norm(NAS - Origin);
    vz = cross(vx,vy);


    transf.R = inv([vx,vy,vz]); % Rotation
    transf.T = - transf.R * Origin; % Translation
    transf.Origin = Origin;
end


if nargin == 2 % User has specified point coordinates to transform
    if nargout == 2 % ... and wants to get them back in SCS 

        varargout{1} = [transf.R, transf.T] * [varargin{1};ones(1,size(varargin{1},2))]; % Compact form of RX+T 
        %varargout{1} = transf.R * varargin{1} + repmat(transf.T,1,size(varargin{1},2)); % Compact form of RX+T 
        
    end
end

