function [GIS,SSF,VAF]=rp_PCA_Image(path1,path2,maskfile,output_path, NumPC)
% This program is used to do SSMPCA on voxel-wise metrics, such as ALFF, ReHo, DC
%
% INPUTS:
%     path1: the folder contains all the metric maps from one group, such as Group A
%     path2: the folder contains all the metric maps from the other group, such as Group B
%     maskfile: the filename of the mask
%     output_path: the directory to save the results, such as GIS, SSF,VAF matrixes and Z-standard GIS map
%     NumPC: the number of principal component (PC) to save as images
%
% OUTPUTS:
%     GIS: Group invariant subprofile, namley networks or PC
%     SSF: Subject scaling factor, namely the projection of  subject data on each PC or the load of each subject for each PC
%     VAF:Variance accounting for, nameley the percentage of the variance for each PC
%
%  EXAMPLES:
%     path1=uigetdir(pwd,'select the folder containing Group A');
%     path2=uigetdir(pwd,'select the folder containing Group B');
%     [FileName,PathName]=uigetfile(pwd,'select the mask file');
%     maskfile=fullfile(PathName,FileName);
%     output_path=uigetdir(pwd,'select the folder to save the results');
%     [GIS,SSF,VAF]=SSMPCA_YLX(path1,path2,maskfile,output_path,5);
%
% Auther: Li-Xia Yuan
% Date: 2020.04.26
% Some hints have been added  SUN Jia-Wei 200506

%% ===============1. read  data===========================
%read the mask
header=spm_vol(maskfile);
brainmask=spm_read_vols(header);%voxel size: 3*3*3
Ind_Mask=find(brainmask>0);
Dim=size(brainmask);

file1=dir([path1,filesep,'*.nii']);
file2=dir([path2,filesep,'*.nii']);
sub=length(file1)+length(file2);
Data=zeros([Dim,sub]);
fprintf('\n\tImage Files in Group 1: \n');% Added by SUN Jia-Wei 200506
% read data from Group A
for i=1:length(file1)
    header=spm_vol([path1,filesep,file1(i).name]);
    ALFF=spm_read_vols(header);
    Data(:,:,:,i)=ALFF;
    fprintf('\t%s\n',file1(i).name)% Added by SUN Jia-Wei 200506
end
% read data from Group B
fprintf('\n\tImage Files in Group 2: \n');% Added by SUN Jia-Wei 200506
for i=1:length(file2)
    header=spm_vol([path2,filesep,file2(i).name]);
    ALFF=spm_read_vols(header);
    Data(:,:,:,length(file1)+i)=ALFF;
    fprintf('\t%s\n',file2(i).name)% Added by SUN Jia-Wei 200506
end

fprintf('\n\tPCA Calculating...\n');% Added by SUN Jia-Wei 200506
%% ===============3. SSM-PCA preprocessing ==================

% mask the data
Dim=size(Data);
Data=reshape(Data,prod(Dim(1:3)),[] );
Data_mask=Data(Ind_Mask,:);

%demean the data
Data_mask=bsxfun(@minus, Data_mask, mean(Data_mask,1));
Data_mask=bsxfun(@minus, Data_mask, mean(Data_mask,2));

%% ===============4. PCA: by SVD===============
[GIS,S,V]=svd(Data_mask,'econ');
SSF=V*S;
S=diag(S);
VAF=S.^2./sum(S.^2)*100;

%% ===============5. standardize GIS===============
zGIS=bsxfun(@minus,GIS,mean(GIS,1));
zGIS=bsxfun(@rdivide,zGIS, std(GIS,0,1)); %rdivide the equivalent A./B divides each entry of A by the corresponding entry of B.

%% ===============Save the results===============
save([output_path,filesep,'SSMPCA_Results.mat'], 'VAF','SSF','zGIS');
% write the zGIS into nii files
for i=1:NumPC;
    rp_waitbar(i/NumPC, 'PCA Calculating...', 'PCA','Child');% Added by SUN Jia-Wei 200506
    Pattern=zeros(prod(Dim(1:3)),1);
    Pattern(Ind_Mask)=zGIS(:,i);
    Pattern=reshape(Pattern,Dim(1:3));
    header.fname=[output_path,filesep,'zGIS',num2str(i),'.nii'];
    spm_write_vol(header,Pattern);
end
rp_waitbar;
fprintf('\n\tPCA Calculation finished.\n');% Added by SUN Jia-Wei 200506
end
