function outputNAME = nic_ImageResample(imageIN,imageOUT,imageREF,method)
% Resampling the input image (imageIN) to generate a new image (imageOUT),
% using head information from reference image (imageREF). i.e., the mask
% image will be rearranged and output with the same format as the reference
% image, which is usually an arbitrary image one wish to apply the mask to.
% The interpolation method can be set in [method].
% NOTE:imageIN,imageREF must be images with nifti format in the standard space
%
% FORMAT [outputNAME]=nic_ImageResample(imageIN,imageOUT,imageREF,method)
%
% Input:
% 1. imageIN    - full file path and name for the image to resample 
%                 [path\*.img]
% 2. imageOUT   - full file path and name for resampled image, [path\*.img]                 
% 3. imageREF   - full file path and name for the reference image 
%                 [path\*.img].The reference image is the image which has 
%                 the same orientation, voxel size, dimension, origin, and 
%                 bounding box as the image you wish to apply the mask to. 
%                 If already having the images to apply the mask to, you 
%                 can simply put one of them as the reference image.
% 4. method     - sets the interpolation method for the resampling (ref.from SPM)
%                 0          Zero-order hold (nearest neighbour).
%                 1          First-order hold (trilinear interpolation).
%                 2->127     Higher order Lagrange (polynomial) interpolation 
%                            using different holds (second-order upwards).
%                -127 - -1   Different orders of sinc interpolation.
%                 if don't set the method, default is 1 (e.g. trilinear 
%                 interpolation) (NOTE: suggest to set to 1.)
%
% Output:
% 1. outputNAME - full file path and name of the resampled image for output
%                 [path\*.img]
%_______________________________________________________________________
% Copyright (C) 2007 Neuroimage Computing Group, State Key Laboratory of
% Cognitive Neuroscience and Learning
% 
% Liang Wang, wanglbit@gmail.com
% Modified by Han Zhang, napoleon1982@gmail.com
% @(#)nic_ImageResample.m  ver 1.0, 07/11/27

if ~exist('method','var')
    method = 1;
end
spmPath=nic_getSPMPath('nifti');
oldDir = pwd;
addpath(spmPath);
cd(spmPath); 
[dataIN, headIN]   = nic_ReadNiftiImage(imageIN);
[dataREF, headREF] = nic_ReadNiftiImage(imageREF);
orientIN           = headIN.mat;
orientREF          = headREF.mat;
dimIN              = headIN.dim;
dataTEMP_0         = dataIN;
headTEMP           = headIN;

% compare input head info with reference head info, if they are not equal,
% change the input head in order to match the reference head. 
if sign(orientIN(1,1))~=sign(orientREF(1,1))
    dataTEMP_1 = flipdim(dataTEMP_0,1);%2
    headTEMP.mat(1,:) = -1*headTEMP.mat(1,:);
else
    dataTEMP_1 = dataTEMP_0;
    headTEMP.mat(1,:) = headTEMP.mat(1,:);
end

if sign(orientIN(2,2))~=sign(orientREF(2,2))
    dataTEMP_2 = flipdim(dataTEMP_1,2);%1
    headTEMP.mat(2,:) = -1*headTEMP.mat(2,:);
else
    dataTEMP_2 = dataTEMP_1;
    headTEMP.mat(2,:) = headTEMP.mat(2,:);
end

if sign(orientIN(3,3))~=sign(orientREF(3,3))
    dataTEMP_3 = flipdim(dataTEMP_2,3);
    headTEMP.mat(3,:) = -1*headTEMP.mat(3,:);
else
    dataTEMP_3 = dataTEMP_2;
    headTEMP.mat(3,:) = headTEMP.mat(3,:);
end

dataTEMP                      = dataTEMP_3;  
% comparing finished, get the rearranged data in dataTEMP

[pathstrIN,fNameIN,extnIN]    = fileparts(deblank(imageIN));
[pathstrOUT,fNameOUT,extnOUT] = fileparts(deblank(imageOUT));
headTEMP.fname                = strcat('x',[fNameIN extnOUT]);
headTEMP                      = nic_spm_create_vol(headTEMP);
% get the modified head infomation as that from reference image

flag1 = nic_WriteMaskImage(dataTEMP,headTEMP,fullfile(pathstrOUT,headTEMP.fname));
% generate the rearranged *.img and *.hdr files

VI       = nic_spm_vol(fullfile(pathstrOUT,headTEMP.fname));
VO       = VI;
VO.fname = [fNameOUT extnOUT];
VR       = nic_spm_vol(imageREF);

VO.mat      = VR.mat;
VO.dim(1:3) = VR.dim(1:3);
VO          = nic_spm_create_vol(VO);  % update the head info using imageREF

for slice_index = 1:VO.dim(3)
    M  = inv(nic_spm_matrix([0 0 -slice_index])*inv(VO.mat)*VI.mat);
    % VI.mat - affine mat of imageIN
    % VO.mat - affine mat of imageREF
    % M is slice selection indicator
    v(:,:,slice_index)  = nic_spm_slice_vol(VI,M,VO.dim(1:2),method);
end

flag2 = nic_WriteMaskImage(v,VO,fullfile(pathstrOUT , VO.fname));

outputNAME = imageOUT;

% delete unuseful files generated
if extnOUT == '.img'
    delete(fullfile(pathstrOUT , headTEMP.fname));
    delete(fullfile(pathstrOUT , ['x' fNameIN '.hdr']));
    if size(deblank(pwd),2)==size(deblank(pathstrOUT),2)
        if sum(pwd==pathstrOUT,2)==size(pwd,2)
        else
            delete(fullfile(pwd , [fNameOUT '.hdr']));
            delete(fullfile(pwd , ['x' fNameIN '.hdr']));
        end
    else
        delete(fullfile(pwd , [fNameOUT '.hdr']));
        delete(fullfile(pwd , ['x' fNameIN '.hdr']));
    end
elseif extnOUT == '.nii'
    delete(fullfile(pathstrOUT, headTEMP.fname));
     if size(deblank(pwd),2)==size(deblank(pathstrOUT),2)
        if sum(pwd==pathstrOUT,2)==size(pwd,2)
        else
            delete(fullfile(pwd , [fNameOUT '.nii']));
            delete(fullfile(pwd,['x' fNameIN '.nii']));
        end
     else
         delete(fullfile(pwd,[fNameOUT '.nii']));
         delete(fullfile(pwd ,['x' fNameIN '.nii']));
    end
end
cd(oldDir); 
