function data = mcaoAssess(mainDir, key, model)
% mainDir is the directory which the input image together with template,
% target, brainmask image belong to.
% key is a 25*3 matrix saving the 3D coordinates of the 25 SIFT keypoints
% model is the trained SVM for prediction.

% First step: change voxel size and orientation
P = spm_select('List',mainDir,'^image.*\.img$'); % select all the images with file names starting from 'image'
V = spm_vol(P);
l = size(V,1);
%matrix = [9 0 0 -576; 0 0 14.4 -428.8; 0 9 0 -695; 0 0 0 1];
% matrix = [9 0 0 -576; 0 0 14.4 -388.8; 0 9 0 -585; 0 0 0 1]; % this is the matrix in the header file that specifies voxel size and orientation
%parfor i = 1:l
%    V0 = V(i);
%    image = spm_read_vols(V0);
%    V0.mat = matrix;
%    spm_write_vol(V0,image);
%end

% Second step: realign each image to the reference image 'BC0011.img', and then
% apply the same transformation matrix so that each image is roughly
% aligned to the template.
% clear V
% P = char('BC0011.img',P);
% spm_realign(P,struct('quality',0.9,'fwhm',5,'rtm',0,'sep',4,'interp',2,'wrap',[0 0 0],'weight','')); % realign to the first image 'BC0011.img'
% P = spm_select('List',mainDir,'^image.*\.img$');
% V = spm_vol(P);
% for i = 1:length
%     V0 = V(i);
%     image = spm_read_vols(V0);
%     matrix = V0.mat; % read the transformation matrix
%     matrix(1,4) = matrix(1,4);
%     matrix(2,4) = matrix(2,4)-40;  % Since all images are already realigned to the first one...
%     matrix(3,4) = matrix(3,4)-110;  % they have the same transformation parameters to the PET template
%     V0.mat = matrix; % write the new transformation matrix
%     spm_write_vol(V0,image);
% end

% Third & fourth step: coregister to the template and reslice by
% 'target.img'
clear V
clear matlabbatch
matlabbatch{1}.spm.spatial.coreg.write.ref = {fullfile(mainDir,'target.img')};
matlabbatch{1}.spm.spatial.coreg.write.source = {fullfile(mainDir,'BC0011.img')};
matlabbatch{1}.spm.spatial.coreg.write.roptions.interp = 4;
matlabbatch{1}.spm.spatial.coreg.write.roptions.wrap = [0 0 0];
matlabbatch{1}.spm.spatial.coreg.write.roptions.mask = 0;
matlabbatch{1}.spm.spatial.coreg.write.roptions.prefix = 'r';
V = spm_vol(P);
VG = spm_vol('template_PET.img');
for i = 1:l
     V0 = V(i);
     x = spm_coreg(VG,V0,struct('cost_fun','mi'));
     Vx = spm_matrix(x(:)');  % transformation matrix for the coregistration.
     image = spm_read_vols(V0);
     V0.mat = Vx\V0.mat;
     spm_write_vol(V0,image);
     matlabbatch{1}.spm.spatial.coreg.write.source = {fullfile(mainDir,V0.fname)};
     spm_jobman('run',matlabbatch);
end

% Fifth & sixth step: replace the right injured hemisphere by the left half
% brain, then normalize to the template.
clear matlabbatch
matlabbatch{1}.spm.spatial.normalise.estwrite.subj.vol = {fullfile(mainDir,'synthetic.img')};
matlabbatch{1}.spm.spatial.normalise.estwrite.subj.resample = {fullfile(mainDir,'rBC0011.img')};
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.biasreg = 0.0001;
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.biasfwhm = 60;
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.tpm = {fullfile(mainDir,'template_PET.img')};
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.affreg = 'mni';
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.reg = [0 0.001 0.5 0.05 0.2];
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.fwhm = 0;
matlabbatch{1}.spm.spatial.normalise.estwrite.eoptions.samp = 3;
matlabbatch{1}.spm.spatial.normalise.estwrite.woptions.bb = [-78 -154 -118
                                                             80 60 6];
matlabbatch{1}.spm.spatial.normalise.estwrite.woptions.vox = [2 2 2];
matlabbatch{1}.spm.spatial.normalise.estwrite.woptions.interp = 4;
matlabbatch{1}.spm.spatial.normalise.estwrite.woptions.prefix = 'w';

P = spm_select('List',mainDir,'^r.*\.img$'); % read the images with prefix 'r' from the last step.
V = spm_vol(P);
l = size(V,1);
for i = 1:l
    V0 = V(i);
    fname = V0.fname;
    image = spm_read_vols(V0);
    imageLeft = image(71:end,:,:);
    image = [flipud(imageLeft);imageLeft];
    V0.fname = 'synthetic.img';
    spm_write_vol(V0,image);
    matlabbatch{1}.spm.spatial.normalise.estwrite.subj.vol = {fullfile(mainDir,'synthetic.img')};
    matlabbatch{1}.spm.spatial.normalise.estwrite.subj.resample = {fullfile(mainDir,fname)};
    spm_jobman('run',matlabbatch);
end

% Seventh step: remove the extracranial part by a mask image.
% This step is optional for users. However, they can still use the brainmask in
% further analyses to remove the extracranial tissues.
% P = spm_select('List',mainDir,'^wr.*\.img$'); % read the images with prefix 'wr' from the last step.
% V = spm_vol(P);
% l = size(V,1);
% Vmask = spm_vol('brainmask.img');
% mask = spm_read_vols(Vmask);
% maskloc = mask >= 0.8;
% parfor i = 1:l
%     V0 = V(i);
%     image = spm_read_vols(V0);
%     image = image.*maskloc;
%     V0.fname = strcat('m',V0.fname);
%     spm_write_vol(V0,image);
% end
% end

% Last step: Smooth images.
clear matlabbatch
matlabbatch{1}.spm.spatial.smooth.data = {fullfile(mainDir,'test.img')};
matlabbatch{1}.spm.spatial.smooth.fwhm = [2 2 4];
matlabbatch{1}.spm.spatial.smooth.dtype = 0;
matlabbatch{1}.spm.spatial.smooth.im = 0;
matlabbatch{1}.spm.spatial.smooth.prefix = 's';
P = spm_select('List',mainDir,'^wr.*\.img$'); % read the images with prefix 'wr' from normalization step. Change to 'mwr' if brainmask is used in the previous step.
V = spm_vol(P);
l = size(V,1);
for i = 1:l
    V0 = V(i);
    fname = V0.fname;
    matlabbatch{1}.spm.spatial.smooth.data = {fullfile(mainDir,fname)};
    spm_jobman('run',matlabbatch);
end
% all the final images will have file names starting with 'smwr...'
% End of preprocessing

% extract the SIFT descriptors
data = dir(fullfile(mainDir,'swr*.img')); % change to 'smwr' if brainmask is used.
nd = length(data);
nk = size(key,1); % in our paper we choose nk=25, the number of keypoints
desc = zeros(nd*nk,768); % since the 3D SIFT algorithm we use will extract a descriptor as a 768 dimensional vector
% key(:,1) = 41-key(:,1); % the 3D SIFT toolbox will flip the images left to right
for i = 1:nd
    datapath = fullfile(mainDir,data(i).name);
    [im, units] = imRead3D(datapath);
    for j = 1:nk
        keys = keypoint3D(key(j,:));
        [keys, ~] = orientation3D(keys, im, units);
        [desc((i-1)*nk+j,:), ~] = extractSift3D(keys, im, units); % refer to the 3D SIFT script and usage
    end
end
% the result desc is a matrix with 25N times 768 dimension
% where N is the number of images.
% End of 3D SIFT

% SVM classification
descriptor = zeros(nd,nk*768); % concatenate all the descriptors of the same image together and form a predictor
parfor i = 1:nd
    features = desc(nk*(i-1)+1:nk*i,:);
    features = features';
    features = features(:);
    descriptor(i,:) = features';
end % rows of descriptor represent predictors/features of each image
label = predict(model,descriptor);
% the i-th entry of label represent the classification result for the i-th
% image. '1' indicates healthy and '-1' indicates mcao.
data = rmfield(data,'folder');
data = rmfield(data,'date');
data = rmfield(data,'bytes');
data = rmfield(data,'isdir');
data = rmfield(data,'datenum');
for i = 1:nd
    data(i).name = data(i).name(10:end);
    data(i).results = label(i);
end
end
