%File: ASLtbx_normalize_nseg8.m
% spatially normalizing brain images using the spm8 new segment deformation
% field.
% deffile - the deformation field generated during new segmentation. If
% ASLtbx_spm8nseg is used for calling spm8 new segment, the deformation
% field will be automatically created.
% img2benorm - the images to be normalized
% interp - 1 means trilinear interpolation, 4 means the 4-th order b-spline
% etc
% opath - output path
% res2mm - 0 will yield images with a resolution of 1.5x1.5x1.5
% and image matrix of 121x145x121, 1 (default) will have a resolution of 2x2x2 and
% 91x109x91
% Ze Wang. All rights reserved. Upenn
function ASLtbx_normalize_nseg8(deffile, img2benorm, interp, res2mm, opath)
% deformation field and mat from the segment output
if nargin<3, interp=4;  end
if nargin<4, res2mm=1;  end  % res2mm will write the image using 2mm isotropic voxel
if nargin<5, opath='';  end
% deffile='/big2/ZEWANG/boldrest/SMK_MED/VSMK010_WK1/VSMK010_WK1_t1/y_VSMK010_WK1_t1_struct.nii';
P      = [repmat(deffile,3,1), [',1,1';',1,2';',1,3']];
V      = spm_vol(P);
Def1    = cell(3,1);
Def1{1} = spm_load_float(V(1));
Def1{2} = spm_load_float(V(2));
Def1{3} = spm_load_float(V(3));
mat1    = V(1).mat;
if res2mm
    % getting an identity mat from the widely used MNI template with 91x109x91
    % matrix and 2x2x2 voxel size
    SPM_path=spm('Dir');
    targetimg=[SPM_path '/canonical/single_subj_T1.nii,1'];
    % Get an identity transform based on an image volume.
    N   = nifti(targetimg);
    d   = [size(N.dat),1];
    d   = d(1:3);
    mat = N.mat;
    Def = cell(3,1);
    [y1,y2,y3] = ndgrid(1:d(1),1:d(2),1:d(3));
    Def{1} = single(y1*mat(1,1) + y2*mat(1,2) + y3*mat(1,3) + mat(1,4));
    Def{2} = single(y1*mat(2,1) + y2*mat(2,2) + y3*mat(2,3) + mat(2,4));
    Def{3} = single(y1*mat(3,1) + y2*mat(3,2) + y3*mat(3,3) + mat(3,4));
    M    = inv(mat1);
    clear d;
    for j=1:size(Def{1},3)
        d0    = {double(Def{1}(:,:,j)), double(Def{2}(:,:,j)),double(Def{3}(:,:,j))};
        d{1}  = M(1,1)*d0{1}+M(1,2)*d0{2}+M(1,3)*d0{3}+M(1,4);
        d{2}  = M(2,1)*d0{1}+M(2,2)*d0{2}+M(2,3)*d0{3}+M(2,4);
        d{3}  = M(3,1)*d0{1}+M(3,2)*d0{2}+M(3,3)*d0{3}+M(3,4);
        Def{1}(:,:,j) = single(spm_sample_vol(Def1{1},d{:},[1,NaN]));
        Def{2}(:,:,j) = single(spm_sample_vol(Def1{2},d{:},[1,NaN]));
        Def{3}(:,:,j) = single(spm_sample_vol(Def1{3},d{:},[1,NaN]));

    end;
else
    Def=Def1;
    mat=mat1;
end
% img2benorm='/big2/ZEWANG/boldrest/SMK_MED/VSMK010_WK1/VSMK010_WK1_t1/t1_noskull.nii,1';
% opath='';
% interp=4;  % trilinear, 4 mean 4-th order b-spline
apply_def(Def,mat, img2benorm,opath, interp);

function ofnames = apply_def(Def,mat,fnames,odir,intrp)
% Warp an image or series of images according to a deformation field

intrp = [intrp*[1 1 1], 0 0 0];
ofnames = cell(size(fnames,1),1);

for i=1:size(fnames,1),
    V = spm_vol(fnames(i,:));
    M = inv(V.mat);
    [pth,nam,ext,num] = spm_fileparts(deblank(fnames(i,:)));
    if isempty(odir)
        % use same path as source image
        opth = pth;
    else
        % use prespecified path
        opth = odir;
    end
    ofnames{i} = fullfile(opth,['w',nam,ext]);
    Vo = struct('fname',ofnames{i},...
                'dim',[size(Def{1},1) size(Def{1},2) size(Def{1},3)],...
                'dt',V.dt,...
                'pinfo',V.pinfo,...
                'mat',mat,...
                'n',V.n,...
                'descrip',V.descrip);
    ofnames{i} = [ofnames{i} num];
    C  = spm_bsplinc(V,intrp);
    Vo = spm_create_vol(Vo);
    for j=1:size(Def{1},3)
        d0    = {double(Def{1}(:,:,j)), double(Def{2}(:,:,j)),double(Def{3}(:,:,j))};
        d{1}  = M(1,1)*d0{1}+M(1,2)*d0{2}+M(1,3)*d0{3}+M(1,4);
        d{2}  = M(2,1)*d0{1}+M(2,2)*d0{2}+M(2,3)*d0{3}+M(2,4);
        d{3}  = M(3,1)*d0{1}+M(3,2)*d0{2}+M(3,3)*d0{3}+M(3,4);
        dat   = spm_bsplins(C,d{:},intrp);
        Vo    = spm_write_plane(Vo,dat,j);
    end;
end;
return;