%mt = wfu_findobject2target(object,target)
%__________________________________________________________________________
%Returns the affine transformation matrix between 2 image files
%Syntax:  mt = wfu_findobject2target(object,target)
%It uses the wfu_image2mni function to compute the respective transformation
%JAM  12-16-04
%__________________________________________________________________________
% The `.mat' files.
%
% This simply contains a 4x4 affine transformation matrix in a variable `M'.
% These files are normally generated by the `realignment' and
% `coregistration' modules.  What these matrixes contain is a mapping from
% the voxel coordinates (x0,y0,z0) (where the first voxel is at coordinate
% (1,1,1)), to coordinates in millimeters (x1,y1,z1).  By default, the
% the new coordinate system is derived from the `origin' and `vox' fields
% of the image header.
%  
% x1 = M(1,1)*x0 + M(1,2)*y0 + M(1,3)*z0 + M(1,4)
% y1 = M(2,1)*x0 + M(2,2)*y0 + M(2,3)*z0 + M(2,4)
% z1 = M(3,1)*x0 + M(3,2)*y0 + M(3,3)*z0 + M(3,4)
%
% Assuming that image1 has a transformation matrix M1, and image2 has a
% transformation matrix M2, the mapping from image1 to image2 is: M2\M1
% (ie. from the coordinate system of image1 into millimeters, followed
% by a mapping from millimeters into the space of image2).
%
% These `.mat' files allow several realignment or coregistration steps to be
% combined into a single operation (without the necessity of resampling the
% images several times).  The `.mat' files are also used by the spatial
% normalisation module.
%__________________________________________________________________________

%---------------------------------------------------
%Here is the logic for the sn3d transformations
%---------------------------------------------------
%object2target = target2tal \ object2tal
%magnet_ax2MNIv = Affine;
%ax2magnet = MF;
%epi2magnet = Mepi;
%epi2ax= MF\Mepi;
%ax2axMNIv = Affine;
%ax2MNIw = MG/Affine;
%epi2MNIw = (MG/Affine) * (MF\Mepi)

%Affine = native_display -> MNI display
%MG	= MNI_display	 -> full tal
%MF	= native_display -> magnet
%sag->ax = M_ax\M_sag
%MNI_origin = inv(MG) * acpc' = MG\acpc'
%full_tal->native_display = Affine/MG
%native_display->full_tal = MG/Affine = inv(Affine/MG)
%native_display->full_tal = MNI_display->full_tal) / (native_display->MNI_display)
%object2fulltal = (MG/Affine) * (MF\M0)
%sag2fulltal	= (native_ax_display ->full_tal) * (sag2ax)

%If normalized object2tal = Mo
%If not normalized:
%	[Dims,Affine,MF,MG,Tr] = load_params(sn3d);
%	object2tal = (MG/Affine) * (MF\Mo) ;

%If normalized target2tal = Mt
%If not normalized:
%	[Dims,Affine,MF,MG,Tr] = load_params(sn3d);
%	target2tal = (MG/Affine) * (MF\Mt) ;

%_______________________________________________
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%_______________________________________________

function [MT,object2mni,target2mni] = wfu_findobject2target(object, target,compute)
%-------------------------------------------------------------------
%Computes the affine transformation for object space to target space
%-------------------------------------------------------------------
object2mni = wfu_image2mni(object);
target2mni = wfu_image2mni(target);
MT = target2mni \ object2mni;

if exist('compute')
	[pth,fn,ext] = fileparts(object);
	flags.mean=0;
	flags.hold = 0;
	flags.which=1;
	flags.mask=0;
	P = {target,object};
	wfu_spm_reslice(P,flags,MT);
	%---------------------------------------------------
	%If you want to write a matfile instead
	%you will need to copy the .img/.hdr to r.img/r.hdr
	%---------------------------------------------------
	%matfile = fullfile(pth,['r' fn '.mat']);
	%V = spm_vol(target);
	%[normalized,mat0,M0,spm99_norm] = wfu_isnormalized(target,V);
	%M = M0 * MT;
	%mat = mat0 * MT;
	%save(matfile,'M','mat');

end
return

%_______________________________________________
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%_______________________________________________


function wfu_sn_compute(object,target)
Vt = spm_vol(target);
Vo = spm_vol(object);
MG = Vt.mat;
MF = Vo.mat;
Transform = [];
d1 = Vt.dim(1:3);
d2 = [8 8 8];
[thdr,tvox,torigin,toff] = wfu_header_info(target);
d3 = tvox;
d4 = torigin;
d5 = Vo.dim(1:3);
[ohdr,ovox,oorigin,ooff] = wfu_header_info(object);
d6 = ovox;
Dims = [d1;d2;d3;d4;d5;d6];
mgc = 960209;
Affine = target2mni/object2mni;
[pth,fn,ext] = fileparts(object);
sn3dfile = fullfile(pth,'object_sn3d.mat');
save(sn3dfile,'mgc','MG','MF','Dims','Transform','Affine');
VG = Vt;
VF = Vo;
flags.smosrc = 8;
flags.smoref = 0;
flags.regtype = 'mni'
flags.weight = '';
flags.cutoff = 25;
flags.nits = 16;
flags.reg = 1;
flags.wtsrc = 0;
Tr = [];
snfile = fullfile(pth,'object_sn.mat');
save(snfile,'VG','VF','flags','Tr','Affine');
%spm_write_sn(object,sn3dfile,flags);
return




function [hdr,vox,origin,off] = wfu_header_info(imagename)
[pth,nam,ext] = fileparts(imagename);
hfname        = fullfile(pth,[nam '.hdr']);
hdr=[];
vox = [];
origin = [];
off = [];
if exist(hfname) == 2,
	% Read as much information as is stored in the ANALYZE header.
	% Assume transverse slices.
	hdr = spm_read_hdr(hfname);
	if isempty(hdr),
		error(sprintf('Can''t read header for "%s"\n', imagename));
	end
	if any(hdr.hist.origin(1:3)),
		origin = hdr.hist.origin(1:3);
	else,
		origin = (hdr.dime.dim(2:4)+1)/2;
	end;
	vox    = hdr.dime.pixdim(2:4);
	if all(vox == 0), vox = [1 1 1]; end;
	off    = -vox.*origin;
end
return



