% This is the main function for MAGEE method.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Before running this program, please enter the required parameters,
%filenames and input image and output image directories bellow
 tic
 DimSPD=1;   %Dimension of response
 num_Subjects=5;   %Number of subjects
 num_cov=3;   %Number of covariates for each dimension of response
 MaxTimePoints=1;    %Maximum number of time points among all subjects
 No_Row=1;   %Number of rows of constraint matrix
 num_iteration=6;   %Number of iterations, at least 5
 if num_iteration<5
     
     exit
 
 
 ExampleDir='D:\PROJECTS\SSPM\MAGEE_example';
 InputFileFolderDir=fullfile(ExampleDir,'InputFileFolder');
 TimePointsFname='TimePointFilename.txt';
 DesignMatrixFname='DesignMatrixName.txt';
 MaskFname='Mask.hdr';
 ConsMatrixFname='linearConstraintMatrixName.txt';
 InputImageFn='InputImageFilename.txt';
 
 InputImageFilesDirectory=fullfile(ExampleDir,'InputImages');
 
 
 %This is the directory containing files and folder above. 
 OutputFilesDirectory=fullfile(ExampleDir,'OutputImages');  % The output images directory.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %You don not need to chage anything bellow.  
 
 TimePointsFname=fullfile(InputFileFolderDir,TimePointsFname);
 DesignMatrixFname=fullfile(InputFileFolderDir,DesignMatrixFname);
 MaskFname=fullfile(InputFileFolderDir,MaskFname);
 ConsMatrixFname=fullfile(InputFileFolderDir,ConsMatrixFname);
 InputImageFn=fullfile(InputFileFolderDir,InputImageFn);
 
 
fprintf('Dimension of response:                    %d\n',DimSPD);
fprintf('Number of sujects:                        %d\n',num_Subjects);
fprintf('Number of covariates:                     %d\n',num_cov);
fprintf('Maximum number of time points:            %d\n',MaxTimePoints); 
fprintf('Number of rows of constraint matrix:      %d\n',No_Row);
fprintf('Number of iterations:                     %d\n\n',num_iteration);
fprintf('Filename containing time points:\n%s\n',TimePointsFname);
fprintf('Filename containing design matrix:\n%s\n',DesignMatrixFname);
fprintf('Mask image name: \n%s\n',MaskFname);
fprintf('Filename containing constraint matrix:\n%s\n',ConsMatrixFname);
fprintf('Filename containing input images:\n%s\n\n',InputImageFn);
fprintf('Input image file directory:\n%s\n',InputImageFilesDirectory);
fprintf('Output image file directory:\n%s\n\n',OutputFilesDirectory);

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% transform mask image into dat file
  Mask_im=load_nii(MaskFname);
  Mask_matrix_1=Mask_im.img;
  [dimX,dimY,dimZ]=size(Mask_matrix_1);
  Mask_matrix=zeros(dimY,dimX,dimZ,'uint8');
  for kk=1:dimZ
     Mask_matrix(:,:,kk)=(Mask_matrix_1(:,:,kk))';
  end
  voxels=dimX*dimY*dimZ;
  pathstr = fileparts(MaskFname);  
  MaskMatrixName=fullfile(pathstr,'mask.dat');
  fid=fopen(MaskMatrixName,'w');
  fwrite(fid,Mask_matrix,'float32');
  fclose(fid);
  
  fprintf('The mask image size : [%d, %d, %d]\n',dimX,dimY,dimZ);
  fprintf('Number of voxels: %d\n', voxels); 
  fprintf('MaskMatrixName: %s\n', MaskMatrixName); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%% transform input images into dat files
 
 InputImageMatrixFilesFullName=fullfile(InputImageFilesDirectory,'InputImageMatrixFilesname.txt');
 fprintf('InputImageMatrixFilesFullName:%s\n',InputImageMatrixFilesFullName);
 fid=fopen(InputImageMatrixFilesFullName,'wt');
 fid1=fopen(InputImageFn,'rt');
 tline = fgetl(fid1);
 ii=1;
 index=0;
while ischar(tline)
   str=strcat(InputImageFilesDirectory,'\',tline,'.nii');
   nii_ID=exist(str,'file');
   if nii_ID~=2   
      str=strcat(InputImageFilesDirectory,'\',tline,'.hdr');
   end
   
   im=load_nii(str);
   matr_1=im.img;
   [dimxr,dimyr,dimzr]=size(matr_1);
   fprintf('The %dth real image data size : [%d, %d, %d]\n',ii,dimxr,dimyr,dimzr);
   matr=zeros(dimyr,dimxr,dimzr,'uint8');
   
   if dimX~=dimxr ||dimY~=dimyr || dimZ~=dimzr
      fprintf('The %dth real image data size does not match those of the mask image!\n',ii);
      index=1;
      break;
   else       
      for kk=1:dimZ
         matr(:,:,kk)=(matr_1(:,:,kk))';
      end
   end
   
   str1=strcat(InputImageFilesDirectory,'/',tline,'.dat');
   fid2=fopen(str1,'w');
   fwrite(fid2,matr,'float32'); 
   fclose(fid2);
   fprintf(fid,'%s\n',str1);
   tline = fgetl(fid1);   
   ii=ii+1;   
end
fclose(fid1);
fclose(fid);
if index==0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%mex MAGEE.cpp  predef.cpp newdelete.cpp %%compile mex file
% MAGEE main program
fprintf('\nData processing, please wait...\n');

  MAGEE(DimSPD,num_Subjects,num_cov,MaxTimePoints,TimePointsFname,DesignMatrixFname,...
      MaskMatrixName,InputImageMatrixFilesFullName,No_Row,ConsMatrixFname,OutputFilesDirectory,dimX,dimY,dimZ,num_iteration);

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%% transform output data into imge file
 fprintf('\nTransforming output data into imges,please wait...\n\n');
 OutputImageMatrixFilesFullName=fullfile(OutputFilesDirectory,'OutputFilenames.txt'); 
 fid=fopen(OutputImageMatrixFilesFullName,'rt');
 Matr_temp=zeros(dimX,dimY,dimZ,'single');
 tline = fgetl(fid);
 mm=1;
 while ischar(tline)
   fidinside=fopen(tline,'r');
   for kk=1:dimZ
       for ii=1:dimX
          Matr_temp(ii,:,kk)=fread(fidinside,dimY,'float32');
       end
   end
   fclose(fidinside);
   %fprintf('%dth output images is being created...\n',mm); 
   [Path,Name]=fileparts(tline);
   imName=fullfile(Path,Name);  
   imName=strcat(imName,'.nii');        
   IM=make_nii(Matr_temp);
   save_nii(IM,imName);
   tline = fgetl(fid);
   mm=mm+1;
 end
 fclose(fid);    

end
fprintf('++++++Your job is finished!!++++++++\n');

toc
