function [PSCBrain_AllWindow, Header] = infodr_dyPerAF_4D(AllVolume, outfile_PerAF, infile_Mask, Parameter)
%-----------------------------------------------------------
%   Copyright(c) 2015
%	Center for Cognition and Brain Disorders, Hangzhou Normal University, Hangzhou 310015, China
%	Written by JIA Xi-Ze 201412
%	http://www.restfmri.net/
% 	Mail to Authors: jxz.rest@gmail.com, jiaxize@foxmail.com


WindowSize = Parameter.WindowSize;
WindowStep = Parameter.WindowStep;
WindowType = Parameter.WindowType;
AResultFilename = outfile_PerAF;
AMaskFilename = infile_Mask;

[PSCBrain_AllWindow, Header]  = dy_PerAF(WindowSize, WindowStep, WindowType, ...
    AllVolume, AMaskFilename, AResultFilename);
end

function [PSCBrain_AllWindow, Header] = dy_PerAF(WindowSize, WindowStep, WindowType, AllVolume, AMaskFilename, AResultFilename)
% Calculate dynamic regional homogeneity (i.e. ReHo)
% FORMAT     [ReHoBrain_AllWindow, Header] = y_reho_Window(WindowSize, WindowStep, WindowType, AllVolume, NVoxel, AMaskFilename, AResultFilename, IsNeedDetrend, Band, TR, TemporalMask, ScrubbingMethod, ScrubbingTiming, Header, CUTNUMBER)
% Input:
%   WindowSize      -   the size of the sliding window
%   WindowStep      -   the step size
%   WindowType      -   the type of window (e.g., hamming)
% 	AllVolume		-	4D data matrix (DimX*DimY*DimZ*DimTimePoints) or the directory of 3D image data file or the filename of one 4D data file
%   NVoxel              The number of the voxel for a given cluster during calculating the KCC (e.g. 27, 19, or 7); Recommand: NVoxel=27;
% 	AMaskFilename		the mask file name, I only compute the point within the mask
%	AResultFilename		the output filename
% Output:
%	PSCBrain_AllWindow       -   The PerAF results of the windows
%   Header          -   The NIfTI Header
%	AResultFilename	the filename of PerAF result
theElapsedTime = cputime;
%  Read the functional images 
% -------------------------------------------------------------------------
fprintf('\n\t Computing PerAF....\twait...');

[volume_data, VoxelSize, theImgFileList, Header] = rp_to4d(AllVolume);

[nDim1, nDim2, nDim3, nDimTimePoints] = size(volume_data);
BrainSize = [nDim1 nDim2 nDim3];
VoxelSize = sqrt(sum(Header.mat(1:3,1:3).^2));
fprintf('\nLoad mask "%s".\n', AMaskFilename);

if ischar(AMaskFilename)
    if ~isempty(AMaskFilename)
        [mask, MaskVox, MaskHead] = rp_readfile(AMaskFilename);
        if ~all(size(mask)==[nDim1 nDim2 nDim3])
            error('The size of Mask (%dx%dx%d) doesn''t match the required size (%dx%dx%d).\n', size(mask), [nDim1 nDim2 nDim3]);
        end
        mask = double(logical(mask));
    else
        mask = ones(nDim1, nDim2, nDim3);
    end
end

volume_data_tmp = reshape(volume_data, [], size(volume_data, 4));

nWindow = fix((nDimTimePoints - WindowSize)/WindowStep) + 1; 
PSCBrain_AllWindow = zeros([nDim1 nDim2 nDim3 nWindow]);

if ischar(WindowType)
    eval(['WindowType = ', WindowType, '(WindowSize);'])
end
WindowMultiplier = repmat(WindowType(:), 1, size(volume_data_tmp, 1));
WindowMultiplier = WindowMultiplier';
for iWindow = 1:nWindow
    fprintf('\nProcessing window %g of total %g windows\n', iWindow, nWindow);
    AllVolumeWindow = volume_data_tmp(:, (iWindow-1)*WindowStep+1:(iWindow-1)*WindowStep+WindowSize, :);
    AllVolumeWindow = AllVolumeWindow.*WindowMultiplier;
    indata_meantp = mean(AllVolumeWindow, 2);
    indata_meantp = repmat(indata_meantp, 1, size(AllVolumeWindow, 2));
    PSC_Brain4d = (AllVolumeWindow-indata_meantp)./indata_meantp;
    PSC_Brain4d = reshape(PSC_Brain4d, size(volume_data, 1), size(volume_data, 2), size(volume_data, 3), size(AllVolumeWindow, 2));
    PSC_Brain4d  = abs(PSC_Brain4d);
    [ndim1,ndim2,ndim3,ndim4] = size(PSC_Brain4d);
    volume_data_2d = reshape(PSC_Brain4d, ndim1*ndim2*ndim3, ndim4);
    matrix_meantp = mean(volume_data_2d, 2);
    matrix_3d = reshape(matrix_meantp, ndim1, ndim2, ndim3);
    PSC_Brain = matrix_3d.*mask;
    PSC_Brain(find(isnan(PSC_Brain))) = 0;
    PSCBrain_AllWindow(:, :, :, iWindow) = PSC_Brain;  
end

Header.pinfo = [1;0;0];
Header.dt    = [16, 0];

rp_Write4DNIfTI(PSCBrain_AllWindow, Header, AResultFilename);
theElapsedTime = cputime - theElapsedTime;
fprintf('\n\tPerAF computation over, elapsed time: %g seconds\n', theElapsedTime);
end