From 56a86b31e10d4e2fc30733f5c81c6991c98be233 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 1 Jun 2022 16:08:58 +0900 Subject: [PATCH 01/29] add working files for merge --- ...nhighpassandvariancenormalize_riken copy.m | 207 ++++++++++++++++++ ...highpassandvariancenormalize_riken_merge.m | 133 +++++++++++ 2 files changed, 340 insertions(+) create mode 100755 ICAFIX/functionhighpassandvariancenormalize_riken copy.m create mode 100755 ICAFIX/functionhighpassandvariancenormalize_riken_merge.m diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken copy.m b/ICAFIX/functionhighpassandvariancenormalize_riken copy.m new file mode 100755 index 0000000..a8e8212 --- /dev/null +++ b/ICAFIX/functionhighpassandvariancenormalize_riken copy.m @@ -0,0 +1,207 @@ +function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) + + + +%UNTITLED4 Summary of this function goes here +% Detailed explanation goes here +% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser +fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) + + +regstring = ''; +dovol = 1; +%%if length(varargin) > 0 && ~strcmp(varargin{1}, '') +%% dovol = 0;%regname is only used for a new surface registration, so will never require redoing volume +%% regstring = varargin{1};%this has the underscore on the front already +if length(varargin) == 1 && ~isempty(varargin{1}) + dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume + regstring = varargin{1}; %this has the underscore on the front already + ndhpvol=2 + ndhpcifti=3 + ndvol=3 + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end +elseif length(varargin) >= 3 + if ~isempty(varagin{4}) + dovol = 0 + regstring = varargin{4} + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end + end + ndhpvol=varargin{1} + ndhpcifti=varargin{2} + ndvol=varargin{3} + if ~isint(ndhpvol) + error('%s: NDHPVOL should be an integer', mfilename); + end + if ~isint(ndhpcifti) + error('%s: NDHPCIFTI should be an integer', mfilename); + end + if ~isint(ndvol) + error('%s: NDVOL should be an integer', mfilename); + end +end + +if dovol > 0 +cts=single(read_avw([fmri '.nii.gz'])); +ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4); +cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); +end + +if hp>=0 + confounds=load([fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf.par']); + confounds=confounds(:,1:6); + confounds=functionnormalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); + confounds=functionnormalise([confounds confounds.*confounds]); + + BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); +else + if dovol > 0 + cts=demean(cts')'; + end +end + +if hp==0 + if dovol > 0 + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + confounds=detrend(confounds); + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); + + cts=detrend(cts')'; + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + end + + BO.cdata=detrend(BO.cdata')'; + ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); +end +if hp>0 + if dovol > 0 + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + call_fsl(sprintf(['fslmaths ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); + if exist([fmri '_hp' num2str(hp) '.nii.gz'], 'file') == 0 % added by Takuya Hayashi + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); + call_fsl(['fslmaths ' fmri '_hp' num2str(hp) '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri '_hp' num2str(hp) '.nii.gz']); + end + cts=single(read_avw([fmri '_hp' num2str(hp) '.nii.gz'])); + cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + end + + BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); + save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); + call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); + grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; + call_fsl('rm Atlas.nii.gz'); + ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); + +end + +%Compute VN +if hp>=0 + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end + + OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping +else + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,ndvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end +end + +%Save VN +if hp>=0 + if dovol > 0 + save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),[fmri '_hp' num2str(hp) '_vn.nii.gz'],'f',[1 1 1 1]); + call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fmri '_hp' num2str(hp) '_vn.nii.gz -d']); + end + + VN=BO; + VN.cdata=OutBO.noise_unst_std; + disp(['saving ' fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii']); + ciftisavereset(VN,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii'],WBC); +end + +%Apply VN and Save HP_VN TCS +if dovol > 0 + cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); + if hp>=0 + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '_vnts.nii.gz'],'f',[1 1 1 1]); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '_vnts.nii.gz -d']); + else + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_vnts.nii.gz'],'f',[1 1 1 1]); + call_fsl(['fslmaths ' fmri '.nii.gz -sub ' fmri '.nii.gz -add ' fmri '_vnts.nii.gz ' fmri '_vnts.nii.gz']); + %call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_vnf.nii.gz -d']); + end +end + +if hp>=0 + BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); + ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dtseries.nii'],WBC); +end + +%Echo Dims +%TSC: add the regstring to the output filename to avoid overwriting +if hp>=0 + if dovol > 0 + dlmwrite([fmri regstring '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); + else + dlmwrite([fmri regstring '_dims.txt'],[OutBO.calcDim],'\t'); + end +else + %TSC: this mode never gets called with a regstring + dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); +end + +dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); +end +%% ---------------------------------------------- + +%% Polynomial detrending function +function Y = detrendpoly(X,p); + +% X: Input data (column major order) +% p: Order of polynomial to remove +% Y: Detrended output + +% Need to define a function to accomplish this, because MATLAB's own DETREND +% is only capable of removing a *linear* trend (i.e., "p=1" only) + + % Check data, must be in column order + [m, n] = size(X); + if (m == 1) + X = X'; + r=n; + else + r=m; + end + + if (~isscalar(p) || p < 0 || p ~= fix(p)) + error('order of polynomial (p) must be a non-negative integer'); + end + + % 5/1/2019 -- Construct the "Vandermonde matrix" (V) scaled to a maximum of 1, for better numerical properties. + % Note that Octave's DETREND function supports arbitrary polynomial orders, but computes V by taking powers + % of [1:r] (rather than [1:r]/r), which is not numerically robust as p increases. + + V = ([1 : r]'/r * ones (1, p + 1)) .^ (ones (r, 1) * [0 : p]); % "Vandermonde" design matrix + + % Cast design matrix to 'single' if the input is also 'single' (which CIFTI will be) + if strcmp(class(X),'single') + V = single(V); + end + + % Use mldivide ('\') as the linear solver, as it has the nice property of generating a warning + % if the solution is rank deficient (in MATLAB at least; Octave doesn't appear to generate a similar warning). + % [In contrast, PINV does NOT generate a warning if the singular values are less than its internal tolerance]. + % Note that even with the scaling of the Vandermonde matrix to a maximum of 1, rank deficiency starts + % becoming a problem at p=6 for data of class 'single' and 1200 time points. + % Rather than explicitly restricting the allowed order here, we'll code a restriction into the calling scripts. + + Y = X - V * (V \ X); % Remove polynomial fit + +end + diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m new file mode 100755 index 0000000..3480e82 --- /dev/null +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -0,0 +1,133 @@ +function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC,ndhpvol,ndhpcifti,ndvol,varargin) + +%UNTITLED4 Summary of this function goes here +% Detailed explanation goes here +% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser +fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) + + +regstring = ''; +dovol = 1; +if length(varargin) > 0 && ~strcmp(varargin{1}, '') + dovol = 0;%regname is only used for a new surface registration, so will never require redoing volume + regstring = varargin{1};%this has the underscore on the front already +end + +if dovol > 0 +cts=single(read_avw([fmri '.nii.gz'])); +ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4); +cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); +end + +if hp>=0 + confounds=load([fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf.par']); + confounds=confounds(:,1:6); + confounds=functionnormalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); + confounds=functionnormalise([confounds confounds.*confounds]); + + BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); +else + if dovol > 0 + cts=demean(cts')'; + end +end + +if hp==0 + if dovol > 0 + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + confounds=detrend(confounds); + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); + + cts=detrend(cts')'; + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + end + + BO.cdata=detrend(BO.cdata')'; + ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); +end +if hp>0 + if dovol > 0 + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + call_fsl(sprintf(['fslmaths ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); + if exist([fmri '_hp' num2str(hp) '.nii.gz'], 'file') == 0 % added by Takuya Hayashi + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); + call_fsl(['fslmaths ' fmri '_hp' num2str(hp) '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri '_hp' num2str(hp) '.nii.gz']); + end + cts=single(read_avw([fmri '_hp' num2str(hp) '.nii.gz'])); + cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + end + + BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); + save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); + call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); + grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; + call_fsl('rm Atlas.nii.gz'); + ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); + +end + +%Compute VN +if hp>=0 + if dovol > 0 + %Outcts=icaDim(cts,0,1,2,2); %Volume fits two distributions to deal with processing interpolation and multiband + Outcts=icaDim(cts,0,1,2,ndhpvol); %Volume fits two distributions to deal with processing interpolation and multiband + end + %OutBO=icaDim(BO.cdata,0,1,2,3); %CIFTI fits three because of the effects of volume to CIFTI mapping and regularization + OutBO=icaDim(BO.cdata,0,1,2,ndhpcifti); %CIFTI fits three because of the effects of volume to CIFTI mapping and regularization +else + if dovol > 0 + %Outcts=icaDim(cts,0,1,2,3); %Volume fits three distributions to deal with processing interpolation and multiband and concatination effects + Outcts=icaDim(cts,0,1,2,ndvol); % More repeated WF + end +end + +%Save VN +if hp>=0 + if dovol > 0 + save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),[fmri '_hp' num2str(hp) '_vn.nii.gz'],'f',[1 1 1 1]); + call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fmri '_hp' num2str(hp) '_vn.nii.gz -d']); + end + + VN=BO; + VN.cdata=OutBO.noise_unst_std; + disp(['saving ' fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii']); + ciftisavereset(VN,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii'],WBC); +end + +%Apply VN and Save HP_VN TCS +if dovol > 0 + cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); + if hp>=0 + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '_vnts.nii.gz'],'f',[1 1 1 1]); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '_vnts.nii.gz -d']); + else + save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_vnts.nii.gz'],'f',[1 1 1 1]); + call_fsl(['fslmaths ' fmri '.nii.gz -sub ' fmri '.nii.gz -add ' fmri '_vnts.nii.gz ' fmri '_vnts.nii.gz']); + %call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_vnf.nii.gz -d']); + end +end + +if hp>=0 + BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); + ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dtseries.nii'],WBC); +end + +%Echo Dims +%TSC: add the regstring to the output filename to avoid overwriting +if hp>=0 + if dovol > 0 + dlmwrite([fmri regstring '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); + else + dlmwrite([fmri regstring '_dims.txt'],[OutBO.calcDim],'\t'); + end +else + %TSC: this mode never gets called with a regstring + dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); +end + +dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); + +end + From 8c64bd488dccb925b87e6ebe6ba4cc00aea583ba Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 1 Jun 2022 16:10:15 +0900 Subject: [PATCH 02/29] add new params and condition --- ...highpassandvariancenormalize_riken_merge.m | 94 +++++++++++++++++-- 1 file changed, 84 insertions(+), 10 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m index 3480e82..a8e8212 100755 --- a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -1,4 +1,6 @@ -function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC,ndhpvol,ndhpcifti,ndvol,varargin) +function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) + + %UNTITLED4 Summary of this function goes here % Detailed explanation goes here @@ -8,9 +10,38 @@ regstring = ''; dovol = 1; -if length(varargin) > 0 && ~strcmp(varargin{1}, '') - dovol = 0;%regname is only used for a new surface registration, so will never require redoing volume - regstring = varargin{1};%this has the underscore on the front already +%%if length(varargin) > 0 && ~strcmp(varargin{1}, '') +%% dovol = 0;%regname is only used for a new surface registration, so will never require redoing volume +%% regstring = varargin{1};%this has the underscore on the front already +if length(varargin) == 1 && ~isempty(varargin{1}) + dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume + regstring = varargin{1}; %this has the underscore on the front already + ndhpvol=2 + ndhpcifti=3 + ndvol=3 + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end +elseif length(varargin) >= 3 + if ~isempty(varagin{4}) + dovol = 0 + regstring = varargin{4} + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end + end + ndhpvol=varargin{1} + ndhpcifti=varargin{2} + ndvol=varargin{3} + if ~isint(ndhpvol) + error('%s: NDHPVOL should be an integer', mfilename); + end + if ~isint(ndhpcifti) + error('%s: NDHPCIFTI should be an integer', mfilename); + end + if ~isint(ndvol) + error('%s: NDVOL should be an integer', mfilename); + end end if dovol > 0 @@ -71,15 +102,13 @@ %Compute VN if hp>=0 if dovol > 0 - %Outcts=icaDim(cts,0,1,2,2); %Volume fits two distributions to deal with processing interpolation and multiband - Outcts=icaDim(cts,0,1,2,ndhpvol); %Volume fits two distributions to deal with processing interpolation and multiband + Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform end - %OutBO=icaDim(BO.cdata,0,1,2,3); %CIFTI fits three because of the effects of volume to CIFTI mapping and regularization - OutBO=icaDim(BO.cdata,0,1,2,ndhpcifti); %CIFTI fits three because of the effects of volume to CIFTI mapping and regularization + + OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping else if dovol > 0 - %Outcts=icaDim(cts,0,1,2,3); %Volume fits three distributions to deal with processing interpolation and multiband and concatination effects - Outcts=icaDim(cts,0,1,2,ndvol); % More repeated WF + Outcts=icaDim(cts,0,1,-1,ndvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform end end @@ -128,6 +157,51 @@ end dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); +end +%% ---------------------------------------------- + +%% Polynomial detrending function +function Y = detrendpoly(X,p); + +% X: Input data (column major order) +% p: Order of polynomial to remove +% Y: Detrended output + +% Need to define a function to accomplish this, because MATLAB's own DETREND +% is only capable of removing a *linear* trend (i.e., "p=1" only) + + % Check data, must be in column order + [m, n] = size(X); + if (m == 1) + X = X'; + r=n; + else + r=m; + end + + if (~isscalar(p) || p < 0 || p ~= fix(p)) + error('order of polynomial (p) must be a non-negative integer'); + end + + % 5/1/2019 -- Construct the "Vandermonde matrix" (V) scaled to a maximum of 1, for better numerical properties. + % Note that Octave's DETREND function supports arbitrary polynomial orders, but computes V by taking powers + % of [1:r] (rather than [1:r]/r), which is not numerically robust as p increases. + + V = ([1 : r]'/r * ones (1, p + 1)) .^ (ones (r, 1) * [0 : p]); % "Vandermonde" design matrix + + % Cast design matrix to 'single' if the input is also 'single' (which CIFTI will be) + if strcmp(class(X),'single') + V = single(V); + end + + % Use mldivide ('\') as the linear solver, as it has the nice property of generating a warning + % if the solution is rank deficient (in MATLAB at least; Octave doesn't appear to generate a similar warning). + % [In contrast, PINV does NOT generate a warning if the singular values are less than its internal tolerance]. + % Note that even with the scaling of the Vandermonde matrix to a maximum of 1, rank deficiency starts + % becoming a problem at p=6 for data of class 'single' and 1200 time points. + % Rather than explicitly restricting the allowed order here, we'll code a restriction into the calling scripts. + + Y = X - V * (V \ X); % Remove polynomial fit end From 68b4edee66895f023b2bcf52f349995b2a3301b7 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 1 Jun 2022 16:11:00 +0900 Subject: [PATCH 03/29] Delete an unnecessary file --- ...nhighpassandvariancenormalize_riken copy.m | 207 ------------------ 1 file changed, 207 deletions(-) delete mode 100755 ICAFIX/functionhighpassandvariancenormalize_riken copy.m diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken copy.m b/ICAFIX/functionhighpassandvariancenormalize_riken copy.m deleted file mode 100755 index a8e8212..0000000 --- a/ICAFIX/functionhighpassandvariancenormalize_riken copy.m +++ /dev/null @@ -1,207 +0,0 @@ -function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) - - - -%UNTITLED4 Summary of this function goes here -% Detailed explanation goes here -% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser -fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) - - -regstring = ''; -dovol = 1; -%%if length(varargin) > 0 && ~strcmp(varargin{1}, '') -%% dovol = 0;%regname is only used for a new surface registration, so will never require redoing volume -%% regstring = varargin{1};%this has the underscore on the front already -if length(varargin) == 1 && ~isempty(varargin{1}) - dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume - regstring = varargin{1}; %this has the underscore on the front already - ndhpvol=2 - ndhpcifti=3 - ndvol=3 - if ~ischar(regstring) - error('%s: REGSTRING should be a string', mfilename); - end -elseif length(varargin) >= 3 - if ~isempty(varagin{4}) - dovol = 0 - regstring = varargin{4} - if ~ischar(regstring) - error('%s: REGSTRING should be a string', mfilename); - end - end - ndhpvol=varargin{1} - ndhpcifti=varargin{2} - ndvol=varargin{3} - if ~isint(ndhpvol) - error('%s: NDHPVOL should be an integer', mfilename); - end - if ~isint(ndhpcifti) - error('%s: NDHPCIFTI should be an integer', mfilename); - end - if ~isint(ndvol) - error('%s: NDVOL should be an integer', mfilename); - end -end - -if dovol > 0 -cts=single(read_avw([fmri '.nii.gz'])); -ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4); -cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); -end - -if hp>=0 - confounds=load([fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf.par']); - confounds=confounds(:,1:6); - confounds=functionnormalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); - confounds=functionnormalise([confounds confounds.*confounds]); - - BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); -else - if dovol > 0 - cts=demean(cts')'; - end -end - -if hp==0 - if dovol > 0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - confounds=detrend(confounds); - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); - - cts=detrend(cts')'; - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); - end - - BO.cdata=detrend(BO.cdata')'; - ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); -end -if hp>0 - if dovol > 0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - call_fsl(sprintf(['fslmaths ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); - if exist([fmri '_hp' num2str(hp) '.nii.gz'], 'file') == 0 % added by Takuya Hayashi - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslmaths ' fmri '_hp' num2str(hp) '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri '_hp' num2str(hp) '.nii.gz']); - end - cts=single(read_avw([fmri '_hp' num2str(hp) '.nii.gz'])); - cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); - end - - BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); - save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); - call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); - grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; - call_fsl('rm Atlas.nii.gz'); - ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); - -end - -%Compute VN -if hp>=0 - if dovol > 0 - Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform - end - - OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping -else - if dovol > 0 - Outcts=icaDim(cts,0,1,-1,ndvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform - end -end - -%Save VN -if hp>=0 - if dovol > 0 - save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),[fmri '_hp' num2str(hp) '_vn.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fmri '_hp' num2str(hp) '_vn.nii.gz -d']); - end - - VN=BO; - VN.cdata=OutBO.noise_unst_std; - disp(['saving ' fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii']); - ciftisavereset(VN,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii'],WBC); -end - -%Apply VN and Save HP_VN TCS -if dovol > 0 - cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); - if hp>=0 - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '_vnts.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '_vnts.nii.gz -d']); - else - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_vnts.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslmaths ' fmri '.nii.gz -sub ' fmri '.nii.gz -add ' fmri '_vnts.nii.gz ' fmri '_vnts.nii.gz']); - %call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_vnf.nii.gz -d']); - end -end - -if hp>=0 - BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); - ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dtseries.nii'],WBC); -end - -%Echo Dims -%TSC: add the regstring to the output filename to avoid overwriting -if hp>=0 - if dovol > 0 - dlmwrite([fmri regstring '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); - else - dlmwrite([fmri regstring '_dims.txt'],[OutBO.calcDim],'\t'); - end -else - %TSC: this mode never gets called with a regstring - dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); -end - -dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); -end -%% ---------------------------------------------- - -%% Polynomial detrending function -function Y = detrendpoly(X,p); - -% X: Input data (column major order) -% p: Order of polynomial to remove -% Y: Detrended output - -% Need to define a function to accomplish this, because MATLAB's own DETREND -% is only capable of removing a *linear* trend (i.e., "p=1" only) - - % Check data, must be in column order - [m, n] = size(X); - if (m == 1) - X = X'; - r=n; - else - r=m; - end - - if (~isscalar(p) || p < 0 || p ~= fix(p)) - error('order of polynomial (p) must be a non-negative integer'); - end - - % 5/1/2019 -- Construct the "Vandermonde matrix" (V) scaled to a maximum of 1, for better numerical properties. - % Note that Octave's DETREND function supports arbitrary polynomial orders, but computes V by taking powers - % of [1:r] (rather than [1:r]/r), which is not numerically robust as p increases. - - V = ([1 : r]'/r * ones (1, p + 1)) .^ (ones (r, 1) * [0 : p]); % "Vandermonde" design matrix - - % Cast design matrix to 'single' if the input is also 'single' (which CIFTI will be) - if strcmp(class(X),'single') - V = single(V); - end - - % Use mldivide ('\') as the linear solver, as it has the nice property of generating a warning - % if the solution is rank deficient (in MATLAB at least; Octave doesn't appear to generate a similar warning). - % [In contrast, PINV does NOT generate a warning if the singular values are less than its internal tolerance]. - % Note that even with the scaling of the Vandermonde matrix to a maximum of 1, rank deficiency starts - % becoming a problem at p=6 for data of class 'single' and 1200 time points. - % Rather than explicitly restricting the allowed order here, we'll code a restriction into the calling scripts. - - Y = X - V * (V \ X); % Remove polynomial fit - -end - From a3c9601f58031372407653241b0b30ee77cce79c Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 8 Jun 2022 14:39:51 +0900 Subject: [PATCH 04/29] add newest version file --- ICAFIX/functionhighpassandvariancenormalize.m | 326 +++++++++++++++--- 1 file changed, 269 insertions(+), 57 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize.m b/ICAFIX/functionhighpassandvariancenormalize.m index 3d0b733..c8a64dd 100644 --- a/ICAFIX/functionhighpassandvariancenormalize.m +++ b/ICAFIX/functionhighpassandvariancenormalize.m @@ -1,94 +1,306 @@ -function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC) -%UNTITLED4 Summary of this function goes here -% Detailed explanation goes here +function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) +% +% FUNCTIONHIGHPASSANDVARIANCENORMALIZE(TR,HP,FMRI,WBC,[REGSTRING]) +% +% This function: +% (1) Detrends / high-pass filters motion confounds, NIFTI (volume) and CIFTI files +% (2) Calls icaDim.m to estimate the dimensionalities of structured and unstructured noise, +% and compute a spatial variance normalization map. +% It is written to support 'hcp_fix_multi_run', and is not intended as a general (stand-alone) +% function supporting detrending and variance normalization. +% +% TR: repetition time between frames, in sec +% HP: determines what high-pass filtering to apply to the motion confounds and data +% HP<0: No highpass applied (demeaning only). No 'hp' string will be added as part of the file names. +% HP>0: The full-width (2*sigma), in sec, to apply using 'fslmaths -bptf'. +% HP='pd#': If HP is a string, in which the first two characters are 'pd', followed by an integer value, +% then polynomial detrending is applied, with the order specified by the integer value +% embedded at the end of the string. +% The output files will include the string '_hppd#' +% HP=0: Gets interpreted as a linear (1st order) detrend +% This is consistent with the convention in FIX of using HP=0 to specify a linear detrend. +% The output file names will include the string '_hp0'. +% FMRI: The base string of the fmri time series (no extensions) +% WBC: wb_command (full path) +% [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. + +% Note: HP='pd0' would be interpreted as a true 0th order detrend, which is +% the same as demeaning. Mathematically, this is the same as the HP<0 condition, +% but the output files will be named differently (i.e., include '_hppd0'), and additional +% output files will be written relative to the HP<0 condition. -cts=single(read_avw([fmri '.nii.gz'])); -ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4); -cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); +% Authors: M. Glasser and M. Harms +CIFTIMatlabReaderWriter=getenv('FSL_FIX_CIFTIRW'); +if (~isdeployed) + % addpath does not work and should not be used in compiled MATLAB + fprintf('Adding %s to MATLAB path\n', CIFTIMatlabReaderWriter); + addpath(CIFTIMatlabReaderWriter); +end + +%% Defaults +dovol = 1; +regstring = ''; +pdflag = false; % Polynomial detrend +pdstring = 'pd'; % Expected string at start of HP variable to request a "polynomial detrend" +hpstring = ''; + +%% Parse varargin +if length(varargin) >= 1 && ~isempty(varargin{1}) + dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume + regstring = varargin{1}; %this has the underscore on the front already + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end +end + +% Check whether polynomial detrend is being requested for the high-pass filtering. +if ischar(hp) + hp = lower(hp); + if strncmp(hp,pdstring,numel(pdstring)) + pdflag = true; + hp = str2double(hp(numel(pdstring)+1:end)); % hp is now a numeric representing the order of the polynomial detrend + elseif ~isempty(str2double(hp)) % Allow for hp to be provided as a string that contains purely numeric elements + hp = str2double(hp); + else error('%s: Invalid specification for the high-pass filter', mfilename); + end +end + +%% Allow for compiled matlab (hp as a string is already handled above) +if (isdeployed) + TR = str2double(TR); +end + +%% Argument checking +if hp==0 && ~pdflag + warning('%s: hp=0 will be interpreted as polynomial detrending of order 1', mfilename); + pdflag=true; + hp=1; + % But in this case, preserve use of "_hp0" in the filename to indicate a linear detrend (for backwards compatibility in file naming) + hpstring = '_hp0'; +end +if pdflag + if (~isscalar(hp) || hp < 0 || hp ~= fix(hp)) + error('%s: Invalid specification for the order of the polynomial detrending', mfilename); + end +end + +%% Finish setting up hpstring to use in file names +if ~pdflag + pdstring = ''; +end +if isempty(hpstring) && hp>=0 + hpstring = ['_hp' pdstring num2str(hp)]; +end + +%% Load the motion confounds, and the CIFTI (if hp>=0) (don't need either if hp<0) if hp>=0 - confounds=load([fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf.par']); + confounds=load([fmri hpstring '.ica/mc/prefiltered_func_data_mcf.par']); confounds=confounds(:,1:6); - confounds=functionnormalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); - confounds=functionnormalise([confounds confounds.*confounds]); + %% normalise function is in HCPPIPEDIR/global/matlab/normalise.m + confounds=normalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); + confounds=normalise([confounds confounds.*confounds]); - BO=ciftiopen([fmri '_Atlas.dtseries.nii'],WBC); -else - cts=demean(cts')'; + BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); end -if hp==0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - confounds=detrend(confounds); - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); - - cts=detrend(cts')'; - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); +%% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI +if pdflag % polynomial detrend case + if dovol > 0 + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + confounds=detrendpoly(confounds,hp); + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + % Note: Use 'range' to identify non-zero voxels (which is very memory efficient) + % rather than 'std' (which requires additional memory equal to the size of the input) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + % Polynomial detrend + cts=detrendpoly(cts',hp)'; + + % Write out result, restoring to original size + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]); + clear ctsfull; + % Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); + end - BO.cdata=detrend(BO.cdata')'; - ciftisave(BO,[fmri '_Atlas_hp' num2str(hp) '.dtseries.nii'],WBC); -end -if hp>0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - call_fsl(sprintf(['fslmaths ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); - - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslmaths ' fmri '_hp' num2str(hp) '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri '_hp' num2str(hp) '.nii.gz']); - cts=single(read_avw([fmri '_hp' num2str(hp) '.nii.gz'])); - cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + BO.cdata=detrendpoly(BO.cdata',hp)'; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + +elseif hp>0 % "fslmaths -bptf" based filtering + if dovol > 0 + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + call_fsl(sprintf(['fslmaths ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); + + % bptf filtering; no masking here, so this is probably memory inefficient + call_fsl(['fslmaths ' fmri '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); + + % Load in result and reduce to just the non-zero voxels (for memory efficiency going forward) + ctsfull=single(read_avw([fmri hpstring '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + end BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); - grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; - call_fsl('rm Atlas.nii.gz'); - ciftisave(BO,[fmri '_Atlas_hp' num2str(hp) '.dtseries.nii'],WBC); + grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + delete('Atlas.nii.gz'); + +elseif hp<0 % If no hp filtering, still need to at least demean the volumetric time series + if dovol > 0 + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + cts=demean(cts')'; + end + end -%Compute VN + +% NOTE: To interpret the logic of the following "Compute, Save, Apply" sections need to know that this function expects that: +% (1) individual run time series will be passed in with hp requested +% (2) concatenated (multi-run) time series will be passed in with hp<0 +%% [This is a bit of hack: would be cleaner if we disentangled the operation of this function from the specific needs of 'hcp_fix_multi_run']. + +%% Compute VN map if hp>=0 - Outcts=icaDim(cts,0,1,2,2); %Volume fits two distributions to deal with processing interpolation and multiband - OutBO=icaDim(BO.cdata,0,1,2,3); %CIFTI fits three because of the effects of volume to CIFTI mapping and regularization + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,2); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end + + OutBO=icaDim(BO.cdata,0,1,-1,3); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping else - Outcts=icaDim(cts,0,1,2,3); %Volume fits three distributions to deal with processing interpolation and multiband and concatination effects + if dovol > 0 + Outcts=icaDim(cts,0,1,-1,2); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + end end -%Save VN +%% Save VN map +% (Only need these for the individual runs, which, per above, are expected to be passed in with hp requested). if hp>=0 - save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),[fmri '_vn.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fmri '_vn.nii.gz -d']); -end + if dovol > 0 + fname=[fmri hpstring '_vn.nii.gz']; + vnfull=zeros(ctsX*ctsY*ctsZ,1, 'single'); + vnfull(ctsmask)=Outcts.noise_unst_std; + + save_avw(reshape(vnfull,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]); + clear vnfull; + call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']); + end -if hp>=0 VN=BO; VN.cdata=OutBO.noise_unst_std; - ciftisavereset(VN,[fmri '_Atlas_vn.dscalar.nii'],WBC); + fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii']; + disp(['saving ' fname]); + ciftisavereset(VN,fname,WBC); end -%Apply VN and Save HP_VN TCS -cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); -if hp>=0 - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '_vn.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '_vn.nii.gz -d']); -else - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_vnf.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslmaths ' fmri '.nii.gz -sub ' fmri '.nii.gz -add ' fmri '_vnf.nii.gz ' fmri '_vnf.nii.gz']); - %call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_vnf.nii.gz -d']); +%% Apply VN and Save HP_VN TCS +% Here, NIFTI VN'ed TCS gets saved regardless of whether hp>=0 (i.e., regardless of whether individual or concatenated run) +% But CIFTI version of TCS only saved if hp>0 +if dovol > 0 + cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); + % Use '_vnts' (volume normalized time series) as the suffix for the volumetric VN'ed TCS + fname=[fmri hpstring '_vnts.nii.gz']; + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]); + clear ctsfull; + % N.B. Version of 'fslcpgeom' in FSL 6.0.0 requires a patch because it doesn't copy both the qform and sform faithfully + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']); end - +% For CIFTI, we can use the extension to distinguish between VN maps (.dscalar) and VN'ed time series (.dtseries) if hp>=0 BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); - ciftisave(BO,[fmri '_Atlas_hp' num2str(hp) '_vn.dtseries.nii'],WBC); + ciftisave(BO,[fmri '_Atlas' regstring hpstring '_vn.dtseries.nii'],WBC); end %Echo Dims +%TSC: add the regstring to the output filename to avoid overwriting if hp>=0 - dlmwrite([fmri '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); + if dovol > 0 + dlmwrite([fmri regstring '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); + else + dlmwrite([fmri regstring '_dims.txt'],[OutBO.calcDim],'\t'); + end else + %TSC: this mode never gets called with a regstring dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); end end + +%% ---------------------------------------------- + +%% Polynomial detrending function +function Y = detrendpoly(X,p); + +% X: Input data (column major order) +% p: Order of polynomial to remove +% Y: Detrended output + +% Need to define a function to accomplish this, because MATLAB's own DETREND +% is only capable of removing a *linear* trend (i.e., "p=1" only) + + % Check data, must be in column order + [m, n] = size(X); + if (m == 1) + X = X'; + r=n; + else + r=m; + end + + if (~isscalar(p) || p < 0 || p ~= fix(p)) + error('order of polynomial (p) must be a non-negative integer'); + end + + % 5/1/2019 -- Construct the "Vandermonde matrix" (V) scaled to a maximum of 1, for better numerical properties. + % Note that Octave's DETREND function supports arbitrary polynomial orders, but computes V by taking powers + % of [1:r] (rather than [1:r]/r), which is not numerically robust as p increases. + + V = ([1 : r]'/r * ones (1, p + 1)) .^ (ones (r, 1) * [0 : p]); % "Vandermonde" design matrix + + % Cast design matrix to 'single' if the input is also 'single' (which CIFTI will be) + if strcmp(class(X),'single') + V = single(V); + end + + % Use mldivide ('\') as the linear solver, as it has the nice property of generating a warning + % if the solution is rank deficient (in MATLAB at least; Octave doesn't appear to generate a similar warning). + % [In contrast, PINV does NOT generate a warning if the singular values are less than its internal tolerance]. + % Note that even with the scaling of the Vandermonde matrix to a maximum of 1, rank deficiency starts + % becoming a problem at p=6 for data of class 'single' and 1200 time points. + % Rather than explicitly restricting the allowed order here, we'll code a restriction into the calling scripts. + + Y = X - V * (V \ X); % Remove polynomial fit + +end + From eab493ec9042b3f7baa0de746abbc8a9f41cd738 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 8 Jun 2022 14:54:39 +0900 Subject: [PATCH 05/29] modified to match HCPv4.3.0 --- ...highpassandvariancenormalize_riken_merge.m | 185 +++++++++++++++--- 1 file changed, 154 insertions(+), 31 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m index a8e8212..022b4b8 100755 --- a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -1,6 +1,45 @@ -function [ output_args ] = functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) +function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) +% +% FUNCTIONHIGHPASSANDVARIANCENORMALIZE(TR,HP,FMRI,WBC,[REGSTRING]) +% +% This function: +% (1) Detrends / high-pass filters motion confounds, NIFTI (volume) and CIFTI files +% (2) Calls icaDim.m to estimate the dimensionalities of structured and unstructured noise, +% and compute a spatial variance normalization map. +% It is written to support 'hcp_fix_multi_run', and is not intended as a general (stand-alone) +% function supporting detrending and variance normalization. +% +% TR: repetition time between frames, in sec +% HP: determines what high-pass filtering to apply to the motion confounds and data +% HP<0: No highpass applied (demeaning only). No 'hp' string will be added as part of the file names. +% HP>0: The full-width (2*sigma), in sec, to apply using 'fslmaths -bptf'. +% HP='pd#': If HP is a string, in which the first two characters are 'pd', followed by an integer value, +% then polynomial detrending is applied, with the order specified by the integer value +% embedded at the end of the string. +% The output files will include the string '_hppd#' +% HP=0: Gets interpreted as a linear (1st order) detrend +% This is consistent with the convention in FIX of using HP=0 to specify a linear detrend. +% The output file names will include the string '_hp0'. +% FMRI: The base string of the fmri time series (no extensions) +% WBC: wb_command (full path) +% [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. + +% Note: HP='pd0' would be interpreted as a true 0th order detrend, which is +% the same as demeaning. Mathematically, this is the same as the HP<0 condition, +% but the output files will be named differently (i.e., include '_hppd0'), and additional +% output files will be written relative to the HP<0 condition. + +% Authors: M. Glasser and M. Harms +CIFTIMatlabReaderWriter=getenv('FSL_FIX_CIFTIRW'); +if (~isdeployed) + % addpath does not work and should not be used in compiled MATLAB + fprintf('Adding %s to MATLAB path\n', CIFTIMatlabReaderWriter); + addpath(CIFTIMatlabReaderWriter); +end +%% Defaults +dovol = 1; %UNTITLED4 Summary of this function goes here % Detailed explanation goes here @@ -9,10 +48,11 @@ regstring = ''; -dovol = 1; -%%if length(varargin) > 0 && ~strcmp(varargin{1}, '') -%% dovol = 0;%regname is only used for a new surface registration, so will never require redoing volume -%% regstring = varargin{1};%this has the underscore on the front already +pdflag = false; % Polynomial detrend +pdstring = 'pd'; % Expected string at start of HP variable to request a "polynomial detrend" +hpstring = ''; + +%% Parse varargin if length(varargin) == 1 && ~isempty(varargin{1}) dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume regstring = varargin{1}; %this has the underscore on the front already @@ -50,11 +90,52 @@ cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); end +% Check whether polynomial detrend is being requested for the high-pass filtering. +if ischar(hp) + hp = lower(hp); + if strncmp(hp,pdstring,numel(pdstring)) + pdflag = true; + hp = str2double(hp(numel(pdstring)+1:end)); % hp is now a numeric representing the order of the polynomial detrend + elseif ~isempty(str2double(hp)) % Allow for hp to be provided as a string that contains purely numeric elements + hp = str2double(hp); + else error('%s: Invalid specification for the high-pass filter', mfilename); + end +end + +%% Allow for compiled matlab (hp as a string is already handled above) +if (isdeployed) + TR = str2double(TR); +end + +%% Argument checking +if hp==0 && ~pdflag + warning('%s: hp=0 will be interpreted as polynomial detrending of order 1', mfilename); + pdflag=true; + hp=1; + % But in this case, preserve use of "_hp0" in the filename to indicate a linear detrend (for backwards compatibility in file naming) + hpstring = '_hp0'; +end +if pdflag + if (~isscalar(hp) || hp < 0 || hp ~= fix(hp)) + error('%s: Invalid specification for the order of the polynomial detrending', mfilename); + end +end + +%% Finish setting up hpstring to use in file names +if ~pdflag + pdstring = ''; +end +if isempty(hpstring) && hp>=0 + hpstring = ['_hp' pdstring num2str(hp)]; +end + +%% Load the motion confounds, and the CIFTI (if hp>=0) (don't need either if hp<0) if hp>=0 - confounds=load([fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf.par']); + confounds=load([fmri hpstring '.ica/mc/prefiltered_func_data_mcf.par']); confounds=confounds(:,1:6); - confounds=functionnormalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); - confounds=functionnormalise([confounds confounds.*confounds]); + %% normalise function is in HCPPIPEDIR/global/matlab/normalise.m + confounds=normalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); + confounds=normalise([confounds confounds.*confounds]); BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); else @@ -63,39 +144,81 @@ end end -if hp==0 +%% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI +if pdflag % polynomial detrend case if dovol > 0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - confounds=detrend(confounds); - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); - - cts=detrend(cts')'; - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + confounds=detrendpoly(confounds,hp); + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],'f',[1 1 1 TR]); + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + % Note: Use 'range' to identify non-zero voxels (which is very memory efficient) + % rather than 'std' (which requires additional memory equal to the size of the input) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + % Polynomial detrend + cts=detrendpoly(cts',hp)'; + + % Write out result, restoring to original size + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),[fmri hpstring '.nii.gz'],'f',[1 1 1 TR]); + clear ctsfull; + % Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); end - BO.cdata=detrend(BO.cdata')'; - ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); -end -if hp>0 + BO.cdata=detrendpoly(BO.cdata',hp)'; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + +elseif hp>0 % "fslmaths -bptf" based filtering if dovol > 0 - save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); - call_fsl(sprintf(['fslmaths ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri '_hp' num2str(hp) '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); - if exist([fmri '_hp' num2str(hp) '.nii.gz'], 'file') == 0 % added by Takuya Hayashi - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '.nii.gz'],'f',[1 1 1 TR]); - call_fsl(['fslmaths ' fmri '_hp' num2str(hp) '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri '_hp' num2str(hp) '.nii.gz']); - end - cts=single(read_avw([fmri '_hp' num2str(hp) '.nii.gz'])); - cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '.nii.gz']); + % Save and filter confounds, as a NIFTI, as expected by FIX + save_avw(reshape(confounds',size(confounds,2),1,1,size(confounds,1)),[fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf'],'f',[1 1 1 TR]); + call_fsl(sprintf(['fslmaths ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf -bptf %f -1 ' fmri hpstring '.ica/mc/prefiltered_func_data_mcf_conf_hp'],0.5*hp/TR)); + + % bptf filtering; no masking here, so this is probably memory inefficient + call_fsl(['fslmaths ' fmri '.nii.gz -bptf ' num2str(0.5*hp/TR) ' -1 ' fmri hpstring '.nii.gz']); + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); + + % Load in result and reduce to just the non-zero voxels (for memory efficiency going forward) + ctsfull=single(read_avw([fmri hpstring '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; end BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; - call_fsl('rm Atlas.nii.gz'); - ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '.dtseries.nii'],WBC); + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + delete('Atlas.nii.gz'); + +elseif hp<0 % If no hp filtering, still need to at least demean the volumetric time series + if dovol > 0 + + % Load volume time series and reduce to just the non-zero voxels (for memory efficiency) + ctsfull=single(read_avw([fmri '.nii.gz'])); + ctsX=size(ctsfull,1); ctsY=size(ctsfull,2); ctsZ=size(ctsfull,3); ctsT=size(ctsfull,4); + ctsfull=reshape(ctsfull,ctsX*ctsY*ctsZ,ctsT); + ctsmask=range(ctsfull, 2) > 0; + fprintf('Non-empty voxels: %d (= %.2f%% of %d)\n', sum(ctsmask), 100*sum(ctsmask)/size(ctsfull,1), size(ctsfull,1)); + cts=ctsfull(ctsmask,:); + clear ctsfull; + + cts=demean(cts')'; + end end From 3d10102ccb7c776a09f81a32da62491d05b4f6b1 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 8 Jun 2022 14:59:48 +0900 Subject: [PATCH 06/29] modified to match HCPv4.3.0 --- ...highpassandvariancenormalize_riken_merge.m | 37 ++++++++++++++----- 1 file changed, 27 insertions(+), 10 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m index 022b4b8..be22a54 100755 --- a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -222,7 +222,13 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) end -%Compute VN + +% NOTE: To interpret the logic of the following "Compute, Save, Apply" sections need to know that this function expects that: +% (1) individual run time series will be passed in with hp requested +% (2) concatenated (multi-run) time series will be passed in with hp<0 +%% [This is a bit of hack: would be cleaner if we disentangled the operation of this function from the specific needs of 'hcp_fix_multi_run']. + +%% Compute VN map if hp>=0 if dovol > 0 Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform @@ -235,20 +241,29 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) end end -%Save VN +%% Save VN map +% (Only need these for the individual runs, which, per above, are expected to be passed in with hp requested). if hp>=0 if dovol > 0 - save_avw(reshape(Outcts.noise_unst_std,ctsX,ctsY,ctsZ,1),[fmri '_hp' num2str(hp) '_vn.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fmri '_hp' num2str(hp) '_vn.nii.gz -d']); - end + fname=[fmri hpstring '_vn.nii.gz']; + vnfull=zeros(ctsX*ctsY*ctsZ,1, 'single'); + vnfull(ctsmask)=Outcts.noise_unst_std; + + save_avw(reshape(vnfull,ctsX,ctsY,ctsZ,1),fname,'f',[1 1 1 1]); + clear vnfull; + call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']); + end VN=BO; VN.cdata=OutBO.noise_unst_std; - disp(['saving ' fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii']); - ciftisavereset(VN,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dscalar.nii'],WBC); + fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii']; + disp(['saving ' fname]); + ciftisavereset(VN,fname,WBC); end -%Apply VN and Save HP_VN TCS +%% Apply VN and Save HP_VN TCS +% Here, NIFTI VN'ed TCS gets saved regardless of whether hp>=0 (i.e., regardless of whether individual or concatenated run) +% But CIFTI version of TCS only saved if hp>0 if dovol > 0 cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); if hp>=0 @@ -260,10 +275,10 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) %call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_vnf.nii.gz -d']); end end - +% For CIFTI, we can use the extension to distinguish between VN maps (.dscalar) and VN'ed time series (.dtseries) if hp>=0 BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); - ciftisave(BO,[fmri '_Atlas' regstring '_hp' num2str(hp) '_vn.dtseries.nii'],WBC); + ciftisave(BO,[fmri '_Atlas' regstring hpstring '_vn.dtseries.nii'],WBC); end %Echo Dims @@ -281,6 +296,8 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); end + + %% ---------------------------------------------- %% Polynomial detrending function From 8d682e50acfda4352c87ff54fc9b80e1625189d2 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 8 Jun 2022 15:01:51 +0900 Subject: [PATCH 07/29] Fixed varargin to accept new params(ndhpvol etc) --- ICAFIX/functionhighpassandvariancenormalize.m | 31 ++++++++++++++++--- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize.m b/ICAFIX/functionhighpassandvariancenormalize.m index c8a64dd..ff19311 100644 --- a/ICAFIX/functionhighpassandvariancenormalize.m +++ b/ICAFIX/functionhighpassandvariancenormalize.m @@ -46,12 +46,35 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) hpstring = ''; %% Parse varargin -if length(varargin) >= 1 && ~isempty(varargin{1}) +if length(varargin) == 1 && ~isempty(varargin{1}) dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume regstring = varargin{1}; %this has the underscore on the front already + ndhpvol=2 + ndhpcifti=3 + ndvol=3 if ~ischar(regstring) error('%s: REGSTRING should be a string', mfilename); end +elseif length(varargin) >= 3 + if ~isempty(varagin{4}) + dovol = 0 + regstring = varargin{4} + if ~ischar(regstring) + error('%s: REGSTRING should be a string', mfilename); + end + end + ndhpvol=varargin{1} + ndhpcifti=varargin{2} + ndvol=varargin{3} + if ~isint(ndhpvol) + error('%s: NDHPVOL should be an integer', mfilename); + end + if ~isint(ndhpcifti) + error('%s: NDHPCIFTI should be an integer', mfilename); + end + if ~isint(ndvol) + error('%s: NDVOL should be an integer', mfilename); + end end % Check whether polynomial detrend is being requested for the high-pass filtering. @@ -191,13 +214,13 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) %% Compute VN map if hp>=0 if dovol > 0 - Outcts=icaDim(cts,0,1,-1,2); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform end - OutBO=icaDim(BO.cdata,0,1,-1,3); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping + OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping else if dovol > 0 - Outcts=icaDim(cts,0,1,-1,2); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + Outcts=icaDim(cts,0,1,-1,ndvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform end end From 00014a25e75e77df9153d6ac38ca9a9a124d8a96 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 8 Jun 2022 15:32:39 +0900 Subject: [PATCH 08/29] mod make concatenated folder to the same as HCP --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 15666b9..1068c70 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -566,17 +566,21 @@ for fmri in $fmris ; do done #AlreadyHP="-1" #Don't run highpass on concatenated data +## --------------------------------------------------------------------------- +## Concatenate the individual runs and create necessary files +## --------------------------------------------------------------------------- #Make Concatenated Folder ConcatFolder=`dirname ${ConcatName}` log_Msg "Making concatenated folder: ${ConcatFolder}" if [ ! -e ${ConcatFolder} ] ; then - mkdir -p ${ConcatFolder} + mkdir ${ConcatFolder} else if [ $OW = 1 ] ; then - rm -r ${ConcatFolder} + /bin/rm -r ${ConcatFolder} + log_Warn "Previous ${ConcatFolder} directory removed" fi - mkdir -p ${ConcatFolder} + mkdir ${ConcatFolder} fi ConcatNameNoExt=$($FSLDIR/bin/remove_ext $ConcatName) # No extension, but still includes the directory path From 219254b777d80842d7e37149ea00013fabecfa36 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 8 Jun 2022 15:35:07 +0900 Subject: [PATCH 09/29] Modified merge volumes to the same as HCP --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 33 ++++++++++++++++++---------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 1068c70..185c60e 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -585,20 +585,31 @@ fi ConcatNameNoExt=$($FSLDIR/bin/remove_ext $ConcatName) # No extension, but still includes the directory path +# Merge volumes from the individual runs if [ -e `remove_ext ${ConcatName}`_fmris.txt ] ; then rm `remove_ext ${ConcatName}`_fmris.txt;fi # Add by TH May 2019 echo ${FileListSTRING} >> `remove_ext ${ConcatName}`_fmris.txt # Add by TH May 2019 -fslmerge -tr `remove_ext ${ConcatName}`_demean ${NIFTIvolMergeSTRING} $tr -fslmerge -tr `remove_ext ${ConcatName}`_hp${hp}_vnts ${NIFTIvolhpVNMergeSTRING} $tr -fslmerge -t `remove_ext ${ConcatName}`_SBRef ${SBRefVolSTRING} -fslmerge -t `remove_ext ${ConcatName}`_mean ${MeanVolSTRING} -fslmerge -t `remove_ext ${ConcatName}`_hp${hp}_vn ${VNVolSTRING} -fslmaths `remove_ext ${ConcatName}`_SBRef -Tmean `remove_ext ${ConcatName}`_SBRef -fslmaths `remove_ext ${ConcatName}`_mean -Tmean `remove_ext ${ConcatName}`_mean -fslmaths `remove_ext ${ConcatName}`_hp${hp}_vn -Tmean `remove_ext ${ConcatName}`_hp${hp}_vn -fslmaths `remove_ext ${ConcatName}`_hp${hp}_vnts -mul `remove_ext ${ConcatName}`_hp${hp}_vn `remove_ext ${ConcatName}`_hp${hp} -fslmaths `remove_ext ${ConcatName}`_demean -add `remove_ext ${ConcatName}`_mean `remove_ext ${ConcatName}` -fslmaths `remove_ext ${ConcatName}`_SBRef -bin `remove_ext ${ConcatName}`_brain_mask # Inserted to create mask to be used in melodic for suppressing memory error - Takuya Hayashi +fslmerge -tr ${ConcatNameNoExt}_demean ${NIFTIvolMergeSTRING} $tr +fslmerge -tr ${ConcatNameNoExt}_hp${hp}_vnts ${NIFTIvolhpVNMergeSTRING} $tr +fslmerge -t ${ConcatNameNoExt}_SBRef ${SBRefVolSTRING} +fslmerge -t ${ConcatNameNoExt}_mean ${MeanVolSTRING} +fslmerge -t ${ConcatNameNoExt}_hp${hp}_vn ${VNVolSTRING} +# Average across runs +fslmaths ${ConcatNameNoExt}_SBRef -Tmean ${ConcatNameNoExt}_SBRef +fslmaths ${ConcatNameNoExt}_mean -Tmean ${ConcatNameNoExt}_mean # "Grand" mean across runs +fslmaths ${ConcatNameNoExt}_demean -add ${ConcatNameNoExt}_mean ${ConcatNameNoExt} + # Preceding line adds back in the "grand" mean + # Resulting file not used below, but want this concatenated version (without HP or VN) to exist +fslmaths ${ConcatNameNoExt}_hp${hp}_vn -Tmean ${ConcatNameNoExt}_hp${hp}_vn # Mean VN map across the individual runs +fslmaths ${ConcatNameNoExt}_hp${hp}_vnts -mul ${ConcatNameNoExt}_hp${hp}_vn ${ConcatNameNoExt}_hp${hp} + # Preceding line adds back in the "grand" mean + # Resulting file not used below, but want this concatenated version (without HP or VN) to exist +fslmaths ${ConcatNameNoExt}_hp${hp}_vn -Tmean ${ConcatNameNoExt}_hp${hp}_vn # Mean VN map across the individual runs +fslmaths ${ConcatNameNoExt}_hp${hp}_vnts -mul ${ConcatNameNoExt}_hp${hp}_vn ${ConcatNameNoExt}_hp${hp} + # Preceding line restores the mean VN map + # The resulting TCS becomes the input to a 2nd pass of VN, which then becomes the input to melodic +fslmaths ${ConcatNameNoExt}_SBRef -bin ${ConcatNameNoExt}_brain_mask + # Preceding line creates mask to be used in melodic for suppressing memory error - Takuya Hayashi # Same thing for the CIFTI From cead08a7e99fbda2b1e6b59a0e617dfcbc273d19 Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 23 Jun 2022 16:28:35 +0900 Subject: [PATCH 10/29] new args for functionhighpassandvariancenormalize --- ICAFIX/ReApplyFixMultiRunPipeline.sh | 2 +- .../ReApplyFixMultiRunPipeline_RIKEN_merge.sh | 31 ++----------------- 2 files changed, 3 insertions(+), 30 deletions(-) diff --git a/ICAFIX/ReApplyFixMultiRunPipeline.sh b/ICAFIX/ReApplyFixMultiRunPipeline.sh index e9e4277..6fe33d5 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline.sh @@ -715,7 +715,7 @@ main() fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${RegString}');" + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndcifti}', '${RegString}');" log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" diff --git a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh index 22e1181..0f06d54 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh @@ -715,8 +715,8 @@ fi fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${RegString}');" - + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndcifti}', '${RegString}');" + log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" @@ -1088,22 +1088,6 @@ fi Start=`echo "${Start} + ${NumTPS}" | bc -l` done - ## --------------------------------------------------------------------------- - ## Remove all the large time series files in ${ConcatFolder} - ## --------------------------------------------------------------------------- - - ## Deleting these files would save a lot of space. - ## But downstream scripts (e.g., RestingStateStats) assume they exist, and - ## if deleted they would therefore need to be re-created "on the fly" later - - # cd ${ConcatFolder} - # log_Msg "Removing large (concatenated) time series files from ${ConcatFolder}" - # $FSLDIR/bin/imrm ${concatfmri} - # $FSLDIR/bin/imrm ${concatfmri}_hp${hp} - # $FSLDIR/bin/imrm ${concatfmri}_hp${hp}_clean - # /bin/rm -f ${concatfmri}_Atlas${RegString}.dtseries.nii - # /bin/rm -f ${concatfmri}_Atlas${RegString}_hp${hp}.dtseries.nii - # /bin/rm -f ${concatfmri}_Atlas${RegString}_hp${hp}_clean.dtseries.nii cd ${DIR} log_Msg "Completed!" @@ -1113,14 +1097,6 @@ fi # "Global" processing - everything above here should be in a function # ------------------------------------------------------------------------------ -set -e # If any command exits with non-zero value, this script exits - -# Establish defaults -G_DEFAULT_REG_NAME="NONE" -G_DEFAULT_LOW_RES_MESH=32 -G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB -G_DEFAULT_MOTION_REGRESSION="FALSE" - # Set global variables g_script_name=$(basename "${0}") @@ -1159,9 +1135,6 @@ log_Check_Env_Var FSL_FIXDIR # Show tool versions show_tool_versions -# Establish default MATLAB run mode -G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB - # Establish default low res mesh for NHP if [ "$SPECIES" = "Macaque" ] ; then G_DEFAULT_LOW_RES_MESH=10 From 76591cdb43a87f3eadf80254c149c79edd38e0ee Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 23 Jun 2022 16:47:13 +0900 Subject: [PATCH 11/29] Modified def params to the same as HCP --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 31 ++++++++++++++-------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 185c60e..705c775 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -437,17 +437,21 @@ FileListSTRING="" # Added by TH May 2019 for fmri in $fmris ; do log_Msg "Top of loop through fmris: fmri: ${fmri}" - NIFTIvolMergeSTRING=`echo "${NIFTIvolMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_demean "` - NIFTIvolhpVNMergeSTRING=`echo "${NIFTIvolhpVNMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}_vnts "` - SBRefVolSTRING=`echo "${SBRefVolSTRING}$($FSLDIR/bin/remove_ext $fmri)_SBRef "` - MeanVolSTRING=`echo "${MeanVolSTRING}$($FSLDIR/bin/remove_ext $fmri)_mean "` - VNVolSTRING=`echo "${VNVolSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}_vn "` - CIFTIMergeSTRING=`echo "${CIFTIMergeSTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii"` - CIFTIhpVNMergeSTRING=`echo "${CIFTIhpVNMergeSTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_hp${hp}_vn.dtseries.nii"` - MeanCIFTISTRING=`echo "${MeanCIFTISTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii "` - VNCIFTISTRING=`echo "${VNCIFTISTRING} -cifti $($FSLDIR/bin/remove_ext $fmri)_Atlas_hp${hp}_vn.dscalar.nii "` - MovementNIFTIMergeSTRING=`echo "${MovementNIFTIMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf.nii.gz "` - MovementNIFTIhpMergeSTRING=`echo "${MovementNIFTIhpMergeSTRING}$($FSLDIR/bin/remove_ext $fmri)_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_hp.nii.gz "` + fmriNoExt=$($FSLDIR/bin/remove_ext $fmri) # $fmriNoExt still includes leading directory components + + # Create necessary strings for merging across runs + # N.B. Some of these files don't exist yet, and are about to get created + NIFTIvolMergeSTRING+="${fmriNoExt}_demean " + NIFTIvolhpVNMergeSTRING+="${fmriNoExt}_hp${hp}_vnts " #These are the individual run, VN'ed *time series* + SBRefVolSTRING+="${fmriNoExt}_SBRef " + MeanVolSTRING+="${fmriNoExt}_mean " + VNVolSTRING+="${fmriNoExt}_hp${hp}_vn " #These are the individual run, VN'ed NIFTI *maps* (created by functionhighpassandvariancenormalize) + CIFTIMergeSTRING+="-cifti ${fmriNoExt}_Atlas_demean.dtseries.nii " + CIFTIhpVNMergeSTRING+="-cifti ${fmriNoExt}_Atlas_hp${hp}_vn.dtseries.nii " + MeanCIFTISTRING+="-cifti ${fmriNoExt}_Atlas_mean.dscalar.nii " + VNCIFTISTRING+="-cifti ${fmriNoExt}_Atlas_hp${hp}_vn.dscalar.nii " #These are the individual run, VN'ed CIFTI *maps* (created by functionhighpassandvariancenormalize) + MovementNIFTIMergeSTRING+="${fmriNoExt}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf " + MovementNIFTIhpMergeSTRING+="${fmriNoExt}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_hp " cd `dirname $fmri` log_Msg "pwd: "$(pwd) @@ -601,17 +605,12 @@ fslmaths ${ConcatNameNoExt}_demean -add ${ConcatNameNoExt}_mean ${ConcatNameNoEx # Preceding line adds back in the "grand" mean # Resulting file not used below, but want this concatenated version (without HP or VN) to exist fslmaths ${ConcatNameNoExt}_hp${hp}_vn -Tmean ${ConcatNameNoExt}_hp${hp}_vn # Mean VN map across the individual runs -fslmaths ${ConcatNameNoExt}_hp${hp}_vnts -mul ${ConcatNameNoExt}_hp${hp}_vn ${ConcatNameNoExt}_hp${hp} - # Preceding line adds back in the "grand" mean - # Resulting file not used below, but want this concatenated version (without HP or VN) to exist -fslmaths ${ConcatNameNoExt}_hp${hp}_vn -Tmean ${ConcatNameNoExt}_hp${hp}_vn # Mean VN map across the individual runs fslmaths ${ConcatNameNoExt}_hp${hp}_vnts -mul ${ConcatNameNoExt}_hp${hp}_vn ${ConcatNameNoExt}_hp${hp} # Preceding line restores the mean VN map # The resulting TCS becomes the input to a 2nd pass of VN, which then becomes the input to melodic fslmaths ${ConcatNameNoExt}_SBRef -bin ${ConcatNameNoExt}_brain_mask # Preceding line creates mask to be used in melodic for suppressing memory error - Takuya Hayashi - # Same thing for the CIFTI ${FSL_FIX_WBC} -cifti-merge `remove_ext ${ConcatName}`_Atlas_demean.dtseries.nii ${CIFTIMergeSTRING} ${FSL_FIX_WBC} -cifti-average `remove_ext ${ConcatName}`_Atlas_mean.dscalar.nii ${MeanCIFTISTRING} From 532018655802807b4359a9340c74f2cbfee46c3e Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 29 Jun 2022 14:24:29 +0900 Subject: [PATCH 12/29] add new params --- ICAFIX/functionhighpassandvariancenormalize.m | 9 ++++ ...highpassandvariancenormalize_riken_merge.m | 41 ++++++++----------- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize.m b/ICAFIX/functionhighpassandvariancenormalize.m index ff19311..82d52e5 100644 --- a/ICAFIX/functionhighpassandvariancenormalize.m +++ b/ICAFIX/functionhighpassandvariancenormalize.m @@ -22,7 +22,9 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) % The output file names will include the string '_hp0'. % FMRI: The base string of the fmri time series (no extensions) % WBC: wb_command (full path) +% varargin: % [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. +% [ndhpvol ndhpcifti ndvol]: number of distributions highpassed volume, highpassed cifti and concatenated volume. OPTIONAL. % Note: HP='pd0' would be interpreted as a true 0th order detrend, which is % the same as demeaning. Mathematically, this is the same as the HP<0 condition, @@ -77,6 +79,12 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) end end +%UNTITLED4 Summary of this function goes here +% Detailed explanation goes here +% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser +fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) + + % Check whether polynomial detrend is being requested for the high-pass filtering. if ischar(hp) hp = lower(hp); @@ -277,6 +285,7 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); end +dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); end diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m index be22a54..899fa2d 100755 --- a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -22,8 +22,10 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) % The output file names will include the string '_hp0'. % FMRI: The base string of the fmri time series (no extensions) % WBC: wb_command (full path) +% varargin: % [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. - +% [ndhpvol ndhpcifti ndvol]: number of distributions highpassed volume, highpassed cifti and concatenated volume. OPTIONAL. + % Note: HP='pd0' would be interpreted as a true 0th order detrend, which is % the same as demeaning. Mathematically, this is the same as the HP<0 condition, % but the output files will be named differently (i.e., include '_hppd0'), and additional @@ -40,13 +42,6 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) %% Defaults dovol = 1; - -%UNTITLED4 Summary of this function goes here -% Detailed explanation goes here -% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser -fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) - - regstring = ''; pdflag = false; % Polynomial detrend pdstring = 'pd'; % Expected string at start of HP variable to request a "polynomial detrend" @@ -84,11 +79,11 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) end end -if dovol > 0 -cts=single(read_avw([fmri '.nii.gz'])); -ctsX=size(cts,1); ctsY=size(cts,2); ctsZ=size(cts,3); ctsT=size(cts,4); -cts=reshape(cts,ctsX*ctsY*ctsZ,ctsT); -end +%UNTITLED4 Summary of this function goes here +% Detailed explanation goes here +% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser +fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) + % Check whether polynomial detrend is being requested for the high-pass filtering. if ischar(hp) @@ -138,10 +133,6 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) confounds=normalise([confounds confounds.*confounds]); BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); -else - if dovol > 0 - cts=demean(cts')'; - end end %% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI @@ -266,14 +257,14 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) % But CIFTI version of TCS only saved if hp>0 if dovol > 0 cts=cts./repmat(Outcts.noise_unst_std,1,ctsT); - if hp>=0 - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_hp' num2str(hp) '_vnts.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_hp' num2str(hp) '_vnts.nii.gz -d']); - else - save_avw(reshape(cts,ctsX,ctsY,ctsZ,ctsT),[fmri '_vnts.nii.gz'],'f',[1 1 1 1]); - call_fsl(['fslmaths ' fmri '.nii.gz -sub ' fmri '.nii.gz -add ' fmri '_vnts.nii.gz ' fmri '_vnts.nii.gz']); - %call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri '_vnf.nii.gz -d']); - end + % Use '_vnts' (volume normalized time series) as the suffix for the volumetric VN'ed TCS + fname=[fmri hpstring '_vnts.nii.gz']; + ctsfull=zeros(ctsX*ctsY*ctsZ,ctsT, 'single'); + ctsfull(ctsmask,:)=cts; + save_avw(reshape(ctsfull,ctsX,ctsY,ctsZ,ctsT),fname,'f',[1 1 1 1]); + clear ctsfull; + % N.B. Version of 'fslcpgeom' in FSL 6.0.0 requires a patch because it doesn't copy both the qform and sform faithfully + call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']); end % For CIFTI, we can use the extension to distinguish between VN maps (.dscalar) and VN'ed time series (.dtseries) if hp>=0 From 46739894d4d799493cdde6e965d35b245675168c Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 29 Jun 2022 14:57:20 +0900 Subject: [PATCH 13/29] Modified to the same as HCP --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 165 +++++++++++++++++---------- 1 file changed, 104 insertions(+), 61 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 705c775..253cda5 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -454,25 +454,29 @@ for fmri in $fmris ; do MovementNIFTIhpMergeSTRING+="${fmriNoExt}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_hp " cd `dirname $fmri` - log_Msg "pwd: "$(pwd) - fmri=`basename $fmri` - fmri=`$FSLDIR/bin/imglob $fmri` - [ `imtest $fmri` != 1 ] && echo No valid 4D_FMRI input file specified && exit 1 + log_Debug_Msg "pwd: "$(pwd) + fmri=`basename $fmri` # After this, $fmri no longer includes the leading directory components + fmri=`$FSLDIR/bin/imglob $fmri` # After this, $fmri will no longer have an extension (if there was one initially) + if [ `$FSLDIR/bin/imtest $fmri` != 1 ]; then + log_Err_Abort "Invalid 4D_FMRI input file specified: ${fmri}" + fi fmri_orig=$fmri FileListSTRING="${FileListSTRING}${fmri}@" # Added by TH May 2019 - mkdir -p ${fmri}_hp$hp.ica/mc + # Get Movement_Regressors.txt into the format expected by functionhighpassandvariancenormalize.m + mkdir -p ${fmri}_hp${hp}.ica/mc if [ -f Movement_Regressors.txt ] ; then - log_Msg "About to create ${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf.par file" - cat Movement_Regressors.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > ${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf.par + log_Debug_Msg "About to create ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf.par file" + cat Movement_Regressors.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf.par else - log_Error "Movement_Regressors.txt not retrieved properly." - exit -1 + log_Err_Abort "Movement_Regressors.txt not retrieved properly." fi - tr=`$FSLDIR/bin/fslval $fmri pixdim4` - log_Msg "tr: $tr" - log_Msg "processing FMRI file $fmri with highpass $hp" + #Demean movement regressors + demeanMovementRegressors Movement_Regressors.txt Movement_Regressors_demean.txt + MovementTXTMergeSTRING+="$(pwd)/Movement_Regressors_demean.txt " + + if [ $OW = 1 ] ; then imrm ${fmri}_demean;fi if [ ! -e ${fmri}_demean.nii.gz ] ; then @@ -482,8 +486,7 @@ for fmri in $fmris ; do ${FSLDIR}/bin/fslmaths $fmri -sub ${fmri}_mean ${fmri}_demean fi - demeanMovementRegressors Movement_Regressors.txt Movement_Regressors_demean.txt - MovementTXTMergeSTRING=`echo "${MovementTXTMergeSTRING}$(pwd)/Movement_Regressors_demean.txt "` + if [ $OW = 1 ] ; then rm -rf $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ; fi if [ ! -e $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ] ; then @@ -491,6 +494,11 @@ for fmri in $fmris ; do ${FSL_FIX_WBC} -cifti-reduce $($FSLDIR/bin/remove_ext $fmri)_Atlas.dtseries.nii MEAN $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii ${FSL_FIX_WBC} -cifti-math "TCS - MEAN" $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii -var TCS $($FSLDIR/bin/remove_ext $fmri)_Atlas.dtseries.nii -var MEAN $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii -select 1 1 -repeat fi + + ## "1st pass" VN on the individual runs; high-pass gets done here as well + tr=`$FSLDIR/bin/fslval $fmri pixdim4` #No checking currently that TR is same across runs + log_Msg "tr: ${tr}" + log_Msg "processing FMRI file ${fmri} with highpass ${hp}" #Highpass Filter and Variance Normalize Data #if [ $hp -gt 0 ] ; then @@ -557,9 +565,12 @@ for fmri in $fmris ; do log_Msg "Dims was not properly estimated by IcaDim" fi - fslmaths $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf.nii.gz -Tmean $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf_mean.nii.gz - fslmaths $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf.nii.gz -sub $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf_mean.nii.gz $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf.nii.gz - $FSLDIR/bin/imrm $(pwd)/${fmri}_hp$hp.ica/mc/prefiltered_func_data_mcf_conf_mean.nii.gz + # Demean the movement regressors (in the 'fake-NIFTI' format returned by functionhighpassandvariancenormalize) + fslmaths ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf -Tmean ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_mean + fslmaths ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf -sub ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_mean ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf + $FSLDIR/bin/imrm ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_mean + + cd ${DIR} # Return to directory where script was launched fmri=${fmri}_hp$hp #cd ${fmri}.ica @@ -567,9 +578,9 @@ for fmri in $fmris ; do #$FSLDIR/bin/imln ../$fmri filtered_func_data #cd .. log_Msg "Bottom of loop through fmris: fmri: ${fmri}" -done -#AlreadyHP="-1" #Don't run highpass on concatenated data +done ### END LOOP (for fmri in $fmris; do) + ## --------------------------------------------------------------------------- ## Concatenate the individual runs and create necessary files ## --------------------------------------------------------------------------- @@ -612,13 +623,12 @@ fslmaths ${ConcatNameNoExt}_SBRef -bin ${ConcatNameNoExt}_brain_mask # Preceding line creates mask to be used in melodic for suppressing memory error - Takuya Hayashi # Same thing for the CIFTI -${FSL_FIX_WBC} -cifti-merge `remove_ext ${ConcatName}`_Atlas_demean.dtseries.nii ${CIFTIMergeSTRING} -${FSL_FIX_WBC} -cifti-average `remove_ext ${ConcatName}`_Atlas_mean.dscalar.nii ${MeanCIFTISTRING} -${FSL_FIX_WBC} -cifti-average `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dscalar.nii ${VNCIFTISTRING} -${FSL_FIX_WBC} -cifti-math "TCS + MEAN" `remove_ext ${ConcatName}`_Atlas.dtseries.nii -var TCS `remove_ext ${ConcatName}`_Atlas_demean.dtseries.nii -var MEAN `remove_ext ${ConcatName}`_Atlas_mean.dscalar.nii -select 1 1 -repeat -${FSL_FIX_WBC} -cifti-merge `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dtseries.nii ${CIFTIhpVNMergeSTRING} -${FSL_FIX_WBC} -cifti-math "TCS * VN" `remove_ext ${ConcatName}`_Atlas_hp${hp}.dtseries.nii -var TCS `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dtseries.nii -var VN `remove_ext ${ConcatName}`_Atlas_hp${hp}_vn.dscalar.nii -select 1 1 -repeat - +${Caret7_Command} -cifti-merge ${ConcatNameNoExt}_Atlas_demean.dtseries.nii ${CIFTIMergeSTRING} +${Caret7_Command} -cifti-average ${ConcatNameNoExt}_Atlas_mean.dscalar.nii ${MeanCIFTISTRING} +${Caret7_Command} -cifti-math "TCS + MEAN" ${ConcatNameNoExt}_Atlas.dtseries.nii -var TCS ${ConcatNameNoExt}_Atlas_demean.dtseries.nii -var MEAN ${ConcatNameNoExt}_Atlas_mean.dscalar.nii -select 1 1 -repeat +${Caret7_Command} -cifti-merge ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dtseries.nii ${CIFTIhpVNMergeSTRING} +${Caret7_Command} -cifti-average ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dscalar.nii ${VNCIFTISTRING} +${Caret7_Command} -cifti-math "TCS * VN" ${ConcatNameNoExt}_Atlas_hp${hp}.dtseries.nii -var TCS ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dtseries.nii -var VN ${ConcatNameNoExt}_Atlas_hp${hp}_vn.dscalar.nii -select 1 1 -repeat # At this point the concatenated VN'ed time series (both volume and CIFTI, following the "1st pass" VN) can be deleted log_Msg "Removing the concatenated VN'ed time series" @@ -652,19 +662,21 @@ done ## Prepare for melodic on concatenated file ## --------------------------------------------------------------------------- -#ConcatFolder=`dirname ${ConcatName}` -cd ${ConcatFolder} - -#Concat fMRI is intensity normalized (not variance normalized) and has no overall mean but is detrended and within run demeaned and variance normalized for both volume and CIFTI -#concatfmri=`basename $(remove_ext ${ConcatName})`_hp${hp} -#concat_fmri_orig=`basename $(remove_ext ${ConcatName})` +# concatfmri is intensity normalized and has no overall mean. +# It is also not variance normalized (the mean VN map *across runs* was "restored" above). +# But is detrended and within run demeaned and variance normalized for both volume and CIFTI +concatfmri=`basename ${ConcatNameNoExt}` # Directory path is now removed +concatfmrihp=${concatfmri}_hp${hp} + +#MPH: Need to be in original directory here, for the ${MovementNIFTI*STRING} variables to work +#if user did not use absolute paths for input arguments +mkdir -p ${ConcatFolder}/${concatfmrihp}.ica +log_Debug_Msg "About to put contents of ${MovementTXTMergeSTRING} in Movement_Regressors_demean.txt file" +cat ${MovementTXTMergeSTRING} > ${ConcatFolder}/Movement_Regressors_demean.txt +mkdir ${ConcatFolder}/${concatfmrihp}.ica/mc +fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp ${MovementNIFTIhpMergeSTRING} $tr +fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf ${MovementNIFTIMergeSTRING} $tr -mkdir -p ${concatfmri}.ica -log_Msg "About to put contents of ${MovementTXTMergeSTRING} in Movement_Regressors_demean.txt file" -cat ${MovementTXTMergeSTRING} > Movement_Regressors_demean.txt -mkdir -p ${concatfmri}.ica/mc -fslmerge -tr ${concatfmri}.ica/mc/prefiltered_func_data_mcf_conf_hp ${MovementNIFTIhpMergeSTRING} $tr -fslmerge -tr ${concatfmri}.ica/mc/prefiltered_func_data_mcf_conf ${MovementNIFTIMergeSTRING} $tr # Added the following line - needed for FIX feature extraction - Takuya Hayashi July 2018 cat Movement_Regressors_demean.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > ${concatfmri}.ica/mc/prefiltered_func_data_mcf.par @@ -672,7 +684,13 @@ cd `dirname $fmri` } RunICA () { -AlreadyHP="-1" #Don't run highpass on concatenated data +## "2nd pass" VN +## This is solely so that we can apply our own variance normalization as part of the spatial ICA, +## rather than using the VN algorithm internal to melodic +## (The resulting ${concatfmrihp}_vnts time series is NOT the one that gets cleaned by FIX). +## See comments above about regarding the use of call_matlab here. + +AlreadyHP="-1" #Don't run highpass on concatenated data (since hp was done on the individual runs) cd ${ConcatFolder} if [ "$tr" = "" ] ; then fmri=`echo ${fmris} | tr ' ' '\n' | tail -1` @@ -739,15 +757,13 @@ log_Msg "melodic_cmd: ${melodic_cmd}" ${melodic_cmd} return_code=$? log_Msg "melodic has been run: return_code = ${return_code}" -log_Msg "melodic has been run: Contents of ${concatfmri}.ica follow" -if [ "${DEBUG}" = "TRUE" ] ; then - ls -lRa `remove_ext ${concatfmri}`.ica +log_Debug_Msg "melodic has been run: Contents of ${concatfmrihp}.ica follow" +if [ ! -z "${log_debugOn}" ] ; then + ls -lRa ${concatfmrihp}.ica fi if [ "${return_code}" -ne "0" ] ; then - log_Error "melodic has returned a non-zero code" - log_Error "Exiting this script with -1 return value." - exit -1 + log_Err_Abort "melodic has returned a non-zero code" fi # At this point (following melodic), ${concatfmrihp}_vnts is no longer necessary @@ -762,13 +778,25 @@ $FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data.ica/concat_data ## Housekeeping related to files expected for FIX ## --------------------------------------------------------------------------- -cd `remove_ext ${concatfmri}`.ica +cd ${concatfmrihp}.ica $FSLDIR/bin/imln ../${concatfmri} filtered_func_data -$FSLDIR/bin/imln filtered_func_data.ica/mask mask +# This is the concated volume time series from the 1st pass VN, with requested +# hp-filtering applied and with the mean VN map multiplied back in +$FSLDIR/bin/imln ../${concatfmrihp} filtered_func_data + +# This is the concated CIFTI time series from the 1st pass VN, with requested +# hp-filtering applied and with the mean VN map multiplied back in +# Unlike 'hcp_fix' (single-run), here we symlink to the hp-filtered CIFTI and +# use "AlreadyHP=-1" to skip any additional filtering in fix_3_clean. +if [ -f ../${concatfmri}_Atlas_hp${hp}.dtseries.nii ] ; then + $FSLDIR/bin/imln ../${concatfmri}_Atlas_hp${hp}.dtseries.nii Atlas.dtseries.nii +fi -if [ `$FSLDIR/bin/imtest ../${concat_fmri_orig}_SBRef` = 1 ] ; then - $FSLDIR/bin/imln ../${concat_fmri_orig}_SBRef mean_func +# Other necessary files (N.B. Movement regressors already prep'ed above) +$FSLDIR/bin/imln filtered_func_data.ica/mask mask +if [ `$FSLDIR/bin/imtest ../${concatfmri}_SBRef` = 1 ] ; then + $FSLDIR/bin/imln ../${concatfmri}_SBRef mean_func else $FSLDIR/bin/imln filtered_func_data.ica/mean mean_func fi @@ -781,7 +809,7 @@ mkdir -p reg cd reg i_am_at=`pwd` -log_Msg "current folder ${i_am_at}" +log_Debug_Msg "current folder ${i_am_at}" $FSLDIR/bin/imln ../../../../T1w_restore_brain highres $FSLDIR/bin/imln ../../../../wmparc wmparc @@ -801,13 +829,20 @@ if [ `$FSLDIR/bin/imtest ../../../../T2w` = 1 ] ; then $FSLDIR/bin/flirt -in veins -ref example_func -applyxfm -init highres2example_func.mat -out veins_exf $FSLDIR/bin/fslmaths veins_exf -mas example_func veins_exf fi + +# Return to ${ConcatFolder} +# Do not use 'cd ${ConcatFolder}', because ${ConcatFolder} may not be an absolute path cd ../.. } RunFix () { +## --------------------------------------------------------------------------- +## Actually run FIX +## --------------------------------------------------------------------------- + AlreadyHP="-1" #Don't run highpass on concatenated data cd ${ConcatFolder} -log_Msg "running FIX" +log_Msg "Running FIX" # Changes to handle user specified training data file if [ "X${TrainingData}" != X ]; then @@ -827,24 +862,28 @@ if [ "X${TrainingData}" != X ]; then # finally, if the TrainingData file is not found, report an error and get out of here if [ ! -f "${TrainingData}" ]; then - log_Error "FIX training data not found: ${TrainingData}" - exit -1 + log_Err_Abort "FIX training data not found: ${TrainingData}" fi # added for fix_0c_reg_masks - Takuya Hayashi cp ${FSLDIR}/etc/flirtsch/ident.mat ${concatfmri}.ica/reg/standard2example_func.mat else - log_Msg "User has NOT specified a training data file" - log_Msg "using training data ${FSL_FIXDIR}/training_files/HCP_hp2000.RData" + # User has not specified a training data file + #TSC: so, let's look for it and give a useful error, rather than pretending they didn't want to do what they said + #TSC: if you WANT to use, say, hp150 with hp2000 training, it should be explicitly requested by the user + #TSC: exception: we will recommend polynomial detrending for performance, and hp 2000 is already trained, and largely indistinguishable + log_Msg "training data file not specified" training_hp="2000" TrainingData=${FSL_FIXDIR}/training_files/HCP_hp${training_hp}.RData fi +log_Msg "using training data file: ${TrainingData}" - +# set up fix command +# use array for whitespace safety, even if the rest of the script isn't if [[ ${doMotionRegression} == "TRUE" ]]; then if [ $SURFGEOM = 0 ] ; then - fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmri}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") + fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmrihp}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") else fix_cmd=("${FSL_FIXDIR}/fix_hcp" "${concatfmri}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") fi @@ -854,14 +893,17 @@ else # Update (11/8/2019): No longer true as of FIX 1.06.12, but since filtering of the CIFTI in MR-FIX occurs in # functionhighpassandvariancenormalize it doesn't matter, so we leave the code unchanged so that we don't # need to add a FIX version dependency to the script - fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmri}.ica" "${TrainingData}" "${FixThresh}") + fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmrihp}.ica" "${TrainingData}" "${FixThresh}") fi - log_Msg "fix_cmd: ${fix_cmd[*]}" +## MPH: The 'fix' script itself will continue to log to its own custom files +## Alert user to where those are +log_Msg "Check ${concatfmrihp}.ica/fix/logMatlab.txt for log output from feature extraction" +log_Msg "Check ${concatfmrihp}.ica/.fix_2b_predict.log for log output from component classification" +log_Msg "Check ${concatfmrihp}.ica/.fix.log for log output from cleanup stage" "${fix_cmd[@]}" - return_code=$? -if [ "${return_code}" -ne "0" ] ; then +if [ "${return_code}" -ne "0" ]; then log_Err_Abort "return_code from fix_cmd: ${return_code}" fi @@ -943,3 +985,4 @@ elif [ "$RunMode" = 3 ] ; then fi cd ${DIR} +log_Msg "Completed!" From a7574e13c276e51777ef5a824074d0134f8bbef3 Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 15:11:03 +0900 Subject: [PATCH 14/29] fix unnecessary change add template mask --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 85 +++++++--------------------- 1 file changed, 21 insertions(+), 64 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 253cda5..9cc2f18 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -79,17 +79,11 @@ e.g. hcp_fix_multi_run rfMRI_REST1_RL/rfMRI_REST1_RL.nii.gz@rfMRI_REST1_RL/rfMR (if launching the script from the '\${StudyFolder}/\${Subject}/MNINonLinear/Results' directory) Optional single character arguments: - -t : Training file (default is HCP_Style_Single_Multirun_Dedrift.RData for Human, NHPHCP_Macaque.RData for Macaque, - NHPHCP_Marmoset.RData for Marmoset) -d : Dimension for group-ICA (default is estimated by icaDIM) - -s : Species - 0: Human (default), 1: Macaque, 2: Marmoset - -r : Run mode - 0: Run all (PreICA+FIX, ICA+FIX, PostICA+FIX, default), 1: Run ICA,FIX and PostICA+FIX, - 2: Run FIX and PostICA+FIX, 3: Run PostICA+FIX) - -g : use surface geometry FIX -w ,, : Wishart filter Ndist (default is 2,3,3) - -T : FIX threshold 0 to 100 (%) (default is 10 in HUman, 40 in Macaque and Marmoset) -h : help message -n : Do not overwrite hp if already calculated + -b : Brain mask EOF } @@ -277,7 +271,6 @@ function interpret_as_bool() ############################################################# - # RData if [ $Species = 0 ] ; then if [ "$TrainingData" = "" ] ; then @@ -368,18 +361,14 @@ WF="2,3,3" OW="1" FixThresh="" SURFGEOM=0 +TemplateMask="" UserDim=NONE; while getopts t:d:s:r:w:nT:hg OPT do case "$OPT" in - "t" ) TrainingData="$OPTARG";; "d" ) UserDim="$OPTARG";; - "s" ) Species="$OPTARG";; - "r" ) RunMode="$OPTARG";; "w" ) WF="$OPTARG";; - "n" ) OW=0;; - "T" ) FixThresh="$OPTARG";; - "g" ) SURFGEOM=1;; + "b" ) TemplateMask="$OPTARG";; "h" ) usage_exit;; * ) usage_exit;; esac @@ -432,7 +421,6 @@ VNCIFTISTRING="" MovementNIFTIMergeSTRING="" MovementNIFTIhpMergeSTRING="" MovementTXTMergeSTRING="" -FileListSTRING="" # Added by TH May 2019 for fmri in $fmris ; do log_Msg "Top of loop through fmris: fmri: ${fmri}" @@ -461,7 +449,6 @@ for fmri in $fmris ; do log_Err_Abort "Invalid 4D_FMRI input file specified: ${fmri}" fi fmri_orig=$fmri - FileListSTRING="${FileListSTRING}${fmri}@" # Added by TH May 2019 # Get Movement_Regressors.txt into the format expected by functionhighpassandvariancenormalize.m mkdir -p ${fmri}_hp${hp}.ica/mc @@ -491,8 +478,9 @@ for fmri in $fmris ; do if [ $OW = 1 ] ; then rm -rf $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ; fi if [ ! -e $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ] ; then log_Msg "Calculating mean and demeaned fmri Atlas" - ${FSL_FIX_WBC} -cifti-reduce $($FSLDIR/bin/remove_ext $fmri)_Atlas.dtseries.nii MEAN $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii - ${FSL_FIX_WBC} -cifti-math "TCS - MEAN" $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii -var TCS $($FSLDIR/bin/remove_ext $fmri)_Atlas.dtseries.nii -var MEAN $($FSLDIR/bin/remove_ext $fmri)_Atlas_mean.dscalar.nii -select 1 1 -repeat + #Demean CIFTI + ${Caret7_Command} -cifti-reduce ${fmri}_Atlas.dtseries.nii MEAN ${fmri}_Atlas_mean.dscalar.nii + ${Caret7_Command} -cifti-math "TCS - MEAN" ${fmri}_Atlas_demean.dtseries.nii -var TCS ${fmri}_Atlas.dtseries.nii -var MEAN ${fmri}_Atlas_mean.dscalar.nii -select 1 1 -repeat fi ## "1st pass" VN on the individual runs; high-pass gets done here as well @@ -526,22 +514,7 @@ for fmri in $fmris ; do rm .fix.functionhighpassandvariancenormalize_riken.log fi - if [ $OW = 1 ] ; then - if [ -e ${fmri}_hp${hp}.nii.gz ] ; then - log_Msg "Found and removing ${fmri}_hp${hp}" - imrm ${fmri}_hp${hp} - fi - if [ -e ${fmri}_Atlas_hp${hp}.dtseries.nii ] ; then - log_Msg "Found and removing ${fmri}_Atlas_hp${hp}.dtseries.nii" - rm -rf ${fmri}_Atlas_hp${hp}.dtseries.nii - fi - fi - if [ -e ${fmri}_hp${hp}.nii.gz ] ; then - log_Msg "Found ${fmri}_hp${hp} and confirming it is already filtered" - v=`fslstats ${fmri}_hp${hp} -M` - if [ `echo "$v > 0.000001" | bc` -eq 1 ] || [ `echo "$v < -0.000001" | bc` -eq 1 ] ; then - log_Msg "Found ${fmri}_hp${hp} is not yet filtered" - imrm ${fmri}_hp${hp} + else log_Msg "Found ${fmri}_hp${hp} is already filtered" fi @@ -591,19 +564,14 @@ log_Msg "Making concatenated folder: ${ConcatFolder}" if [ ! -e ${ConcatFolder} ] ; then mkdir ${ConcatFolder} else - if [ $OW = 1 ] ; then /bin/rm -r ${ConcatFolder} log_Warn "Previous ${ConcatFolder} directory removed" - fi mkdir ${ConcatFolder} fi ConcatNameNoExt=$($FSLDIR/bin/remove_ext $ConcatName) # No extension, but still includes the directory path # Merge volumes from the individual runs -if [ -e `remove_ext ${ConcatName}`_fmris.txt ] ; then rm `remove_ext ${ConcatName}`_fmris.txt;fi # Add by TH May 2019 -echo ${FileListSTRING} >> `remove_ext ${ConcatName}`_fmris.txt # Add by TH May 2019 - fslmerge -tr ${ConcatNameNoExt}_demean ${NIFTIvolMergeSTRING} $tr fslmerge -tr ${ConcatNameNoExt}_hp${hp}_vnts ${NIFTIvolhpVNMergeSTRING} $tr fslmerge -t ${ConcatNameNoExt}_SBRef ${SBRefVolSTRING} @@ -677,9 +645,6 @@ mkdir ${ConcatFolder}/${concatfmrihp}.ica/mc fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp ${MovementNIFTIhpMergeSTRING} $tr fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf ${MovementNIFTIMergeSTRING} $tr -# Added the following line - needed for FIX feature extraction - Takuya Hayashi July 2018 -cat Movement_Regressors_demean.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > ${concatfmri}.ica/mc/prefiltered_func_data_mcf.par - cd `dirname $fmri` } @@ -704,17 +669,17 @@ if [ "${DEBUG}" = "TRUE" ] ; then ls -lRa `remove_ext ${concatfmri}`.ica fi -log_Msg "Running melodic located at: ${FSLDIR}/bin/melodic" -log_Msg "Beginning of melodic version log, help, and checksum" -if [ "${DEBUG}" = "TRUE" ] ; then - log_Msg "${FSLDIR}/bin/melodic --version" - ${FSLDIR}/bin/melodic --version - log_Msg "${FSLDIR}/bin/melodic --help" - ${FSLDIR}/bin/melodic --help - log_Msg "md5sum ${FSLDIR}/bin/melodic" - md5sum ${FSLDIR}/bin/melodic +log_Msg "Running MELODIC located at: $MELODIC" +log_Debug_Msg "Beginning of melodic version log, help, and checksum" +if [ ! -z "${log_debugOn}" ] ; then + log_Debug_Msg "$MELODIC --version" + $MELODIC --version + log_Debug_Msg "$MELODIC --help" + $MELODIC --help + log_Debug_Msg "md5sum $MELODIC" + md5sum $MELODIC fi -log_Msg "End of melodic version log, help, and checksum" +log_Debug_Msg "End of melodic version log, help, and checksum" ###Run on concatenated file -- @@ -818,12 +783,11 @@ $FSLDIR/bin/makerot --theta=0 > highres2example_func.mat if [ `$FSLDIR/bin/imtest ../../../../T2w` = 1 ] ; then $FSLDIR/bin/fslmaths ../../../../T1w -div ../../../../T2w veins -odt float # Added for NHP - if [ $SPECIES = Human ] ; then - $FSLDIR/bin/flirt -in ${FSL_FIXDIR}/mask_files/hcp_0.7mm_brain_mask -ref veins -out veinbrainmask -applyxfm + if [[ $TemplateMask != "" ]] ; then + $FSLDIR/bin/fslmaths $TemplateMask -bin veinbrainmask else - $FSLDIR/bin/fslmaths $TemplateMask -bin veinbrainmask + $FSLDIR/bin/flirt -in ../../../../brainmask_fs -ref veins -out veinbrainmask -applyxfm fi - #$FSLDIR/bin/flirt -in ${FSL_FIXDIR}/mask_files/hcp_0.7mm_brain_mask -ref veins -out veinbrainmask -applyxfm $FSLDIR/bin/fslmaths veinbrainmask -bin veinbrainmask $FSLDIR/bin/fslmaths veins -div `$FSLDIR/bin/fslstats veins -k veinbrainmask -P 50` -mul 2.18 -thr 10 -min 50 -div 50 veins $FSLDIR/bin/flirt -in veins -ref example_func -applyxfm -init highres2example_func.mat -out veins_exf @@ -835,7 +799,6 @@ fi cd ../.. } -RunFix () { ## --------------------------------------------------------------------------- ## Actually run FIX ## --------------------------------------------------------------------------- @@ -882,12 +845,7 @@ log_Msg "using training data file: ${TrainingData}" # set up fix command # use array for whitespace safety, even if the rest of the script isn't if [[ ${doMotionRegression} == "TRUE" ]]; then - if [ $SURFGEOM = 0 ] ; then fix_cmd=("${FSL_FIXDIR}/fix" "${concatfmrihp}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") - else - fix_cmd=("${FSL_FIXDIR}/fix_hcp" "${concatfmri}.ica" "${TrainingData}" "${FixThresh}" -m -h "${AlreadyHP}") - fi - else #-h is actually a subargument to -m, and will cause problems if specified without (or even not directly following) -m # Update (11/8/2019): No longer true as of FIX 1.06.12, but since filtering of the CIFTI in MR-FIX occurs in @@ -909,7 +867,7 @@ fi log_Msg "Done running FIX" -} + RunPostFix () { cd ${ConcatFolder} @@ -923,7 +881,6 @@ if [ -f ${concatfmri}.ica/Atlas_clean.dtseries.nii ] ; then log_Msg "Moving ${concatfmri}.ica/Atlas_clean.dtseries.nii to ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii" /bin/mv ${concatfmri}.ica/Atlas_clean.dtseries.nii ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii /bin/mv ${concatfmri}.ica/Atlas_clean_vn.dscalar.nii ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii - /bin/mv ${concatfmri}.ica/Atlas_unst.dtseries.nii ${concat_fmri_orig}_Atlas_hp${hp}_unst.dtseries.nii ${FSL_FIX_WBC} -cifti-math "TCS/VN" ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dtseries.nii -var TCS ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii -var VN ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii -select 1 1 -repeat fi From 36242edd5820ce1f473a1e131435a8b2c779d953 Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 15:13:27 +0900 Subject: [PATCH 15/29] rm OW option --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 16 +--------------- 1 file changed, 1 insertion(+), 15 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 9cc2f18..201cfe3 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -82,7 +82,6 @@ e.g. hcp_fix_multi_run rfMRI_REST1_RL/rfMRI_REST1_RL.nii.gz@rfMRI_REST1_RL/rfMR -d : Dimension for group-ICA (default is estimated by icaDIM) -w ,, : Wishart filter Ndist (default is 2,3,3) -h : help message - -n : Do not overwrite hp if already calculated -b : Brain mask EOF } @@ -358,7 +357,6 @@ shift;shift;shift;shift Species=0 RunMode=0 WF="2,3,3" -OW="1" FixThresh="" SURFGEOM=0 TemplateMask="" @@ -463,25 +461,13 @@ for fmri in $fmris ; do demeanMovementRegressors Movement_Regressors.txt Movement_Regressors_demean.txt MovementTXTMergeSTRING+="$(pwd)/Movement_Regressors_demean.txt " - - - if [ $OW = 1 ] ; then imrm ${fmri}_demean;fi - if [ ! -e ${fmri}_demean.nii.gz ] ; then - log_Msg "processing FMRI mean and demean" - #Demean volume, movement regressors, and data + #Demean volumes ${FSLDIR}/bin/fslmaths $fmri -Tmean ${fmri}_mean ${FSLDIR}/bin/fslmaths $fmri -sub ${fmri}_mean ${fmri}_demean - fi - - - if [ $OW = 1 ] ; then rm -rf $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ; fi - if [ ! -e $($FSLDIR/bin/remove_ext $fmri)_Atlas_demean.dtseries.nii ] ; then - log_Msg "Calculating mean and demeaned fmri Atlas" #Demean CIFTI ${Caret7_Command} -cifti-reduce ${fmri}_Atlas.dtseries.nii MEAN ${fmri}_Atlas_mean.dscalar.nii ${Caret7_Command} -cifti-math "TCS - MEAN" ${fmri}_Atlas_demean.dtseries.nii -var TCS ${fmri}_Atlas.dtseries.nii -var MEAN ${fmri}_Atlas_mean.dscalar.nii -select 1 1 -repeat - fi ## "1st pass" VN on the individual runs; high-pass gets done here as well tr=`$FSLDIR/bin/fslval $fmri pixdim4` #No checking currently that TR is same across runs From 43d8fd30a78acdb44428f4552a97b0d5f60e004b Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 15:15:24 +0900 Subject: [PATCH 16/29] rm RunMode --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 201cfe3..fc63536 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -355,7 +355,6 @@ fi # Inserted by Takuya Hayashi, July 2018 shift;shift;shift;shift Species=0 -RunMode=0 WF="2,3,3" FixThresh="" SURFGEOM=0 @@ -634,7 +633,6 @@ fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_co cd `dirname $fmri` } -RunICA () { ## "2nd pass" VN ## This is solely so that we can apply our own variance normalization as part of the spatial ICA, ## rather than using the VN algorithm internal to melodic @@ -783,7 +781,6 @@ fi # Return to ${ConcatFolder} # Do not use 'cd ${ConcatFolder}', because ${ConcatFolder} may not be an absolute path cd ../.. -} ## --------------------------------------------------------------------------- ## Actually run FIX @@ -855,7 +852,6 @@ log_Msg "Done running FIX" -RunPostFix () { cd ${ConcatFolder} log_Msg "Moving ${concatfmri}.ica/filtered_func_data_clean to ${concatfmri}_clean" @@ -915,17 +911,9 @@ done log_Msg "Finished" -} -if [ "$RunMode" = 0 ] ; then - RunPreFix;RunICA;RunFix;RunPostFix -elif [ "$RunMode" = 1 ] ; then - RunICA; RunFix;RunPostFix -elif [ "$RunMode" = 2 ] ; then - RunFix;RunPostFix -elif [ "$RunMode" = 3 ] ; then - RunPostFix -fi + + cd ${DIR} log_Msg "Completed!" From 26e1f243bbb8ef76f6aea556af439c3473050cdd Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 15:16:21 +0900 Subject: [PATCH 17/29] rm RData --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 27 --------------------------- 1 file changed, 27 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index fc63536..b0feeeb 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -270,33 +270,6 @@ function interpret_as_bool() ############################################################# -# RData -if [ $Species = 0 ] ; then - if [ "$TrainingData" = "" ] ; then - TrainingData=${FSL_FIXDIR}/training_files/HCP_Style_Single_Multirun_Dedrift.RData; - fi - export SPECIES=Human - if [ "$FixThresh" = "" ] ; then - FixThresh=10 - fi -elif [ $Species = 1 ] ; then - export SPECIES=Macaque - if [ "$TrainingData" = "" ] ; then - TrainingData=${FSL_FIXDIR}/training_files/NHPHCP_Macaque.RData - fi - if [ "$FixThresh" = "" ] ; then - FixThresh=40 - fi -elif [ $Species = 2 ] ; then - export SPECIES=Marmoset - if [ "$TrainingData" = "" ] ; then - TrainingData=${FSL_FIXDIR}/training_files/NHPHCP_Marmoset.RData - fi - if [ "$FixThresh" = "" ] ; then - FixThresh=40 - fi -fi - ## --------------------------------------------------------------------------- ## Check whether FSL version is recent enough ## --------------------------------------------------------------------------- From 72c4a1631a4cca214a5a848f4cc57039d91c72a7 Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 15:39:35 +0900 Subject: [PATCH 18/29] fix matlab run --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 97 ++++++++++++++-------------- 1 file changed, 47 insertions(+), 50 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index b0feeeb..8ad0654 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -362,7 +362,6 @@ log_Msg "doMotionRegression: ${doMotionRegression}" log_Msg "TrainingData: ${TrainingData}" log_Msg "FixThresh: ${FixThresh}" log_Msg "DeleteIntermediates: ${DeleteIntermediates}" -log_Msg "SPECIES: ${SPECIES}" log_Msg "WF Ndist: ${WF}" log_Msg "UserDim: ${UserDim}" log_Msg "CARET7DIR: ${CARET7DIR}" @@ -418,7 +417,6 @@ for fmri in $fmris ; do if [ `$FSLDIR/bin/imtest $fmri` != 1 ]; then log_Err_Abort "Invalid 4D_FMRI input file specified: ${fmri}" fi - fmri_orig=$fmri # Get Movement_Regressors.txt into the format expected by functionhighpassandvariancenormalize.m mkdir -p ${fmri}_hp${hp}.ica/mc @@ -446,55 +444,59 @@ for fmri in $fmris ; do log_Msg "tr: ${tr}" log_Msg "processing FMRI file ${fmri} with highpass ${hp}" - #Highpass Filter and Variance Normalize Data - #if [ $hp -gt 0 ] ; then - # log_Msg "running highpass" - # hptr=`echo "10 k $hp 2 / $tr / p" | dc -` - # log_Msg "hptr $hptr" - # ${FSLDIR}/bin/fslmaths $fmri -sub ${fmri}_mean -bptf $hptr -1 ${fmri}_hp$hp - #fi - - #log_Msg "functionmotionconfounds log file is to be named: .fix.functionmotionconfounds.log instead of .fix.log" - #cd ${fmri}_hp$hp.ica - #${FSL_FIXDIR}/call_matlab.sh -l .fix.functionmotionconfounds.log -f functionmotionconfounds $tr $hp - #cd .. - - #${FSL_FIX_WBC} -cifti-convert -to-nifti ${fmri}_Atlas.dtseries.nii ${fmri}_Atlas_FAKENIFTI.nii.gz - #${FSLDIR}/bin/fslmaths ${fmri}_Atlas_FAKENIFTI.nii.gz -bptf $hptr -1 ${fmri}_Atlas_hp$hp_FAKENIFTI.nii.gz - #${FSL_FIX_WBC} -cifti-convert -from-nifti ${fmri}_Atlas_hp$hp_FAKENIFTI.nii.gz ${fmri}_Atlas.dtseries.nii ${fmri}_Atlas_hp$hp.dtseries.nii - #$FSLDIR/bin/imrm ${fmri}_Atlas_FAKENIFTI ${fmri}_Atlas_hp$hp_FAKENIFTI - - #cd ${fmri}_hp$hp.ica + ## MPH: Log output to existing stdout/stderr, rather than a separate .fix.functionhighpassandvariancenormalize.log file #if [ -e .fix.functionhighpassandvariancenormalize.log ] ; then - # rm .fix.functionhighpassandvariancenormalize.log + # /bin/rm .fix.functionhighpassandvariancenormalize.log #fi - if [ -e .fix.functionhighpassandvariancenormalize_riken.log ] ; then - rm .fix.functionhighpassandvariancenormalize_riken.log - fi - + ## We will use the call_matlab.sh script available in ${FSL_FIXDIR}. + ## But the arguments need to be customized depending on the matlab mode, because: + ## (1) the -r argument does not support the use of compiled matlab, but is needed for interpreted matlab/octave + ## to allow addition of 'addpath' as a command, so that functionhighpassandvariancenormalize.m can be found within matlab/octave. + ## (2) for compiled matlab, we need to use the -c and -b flags to specify the location of the MCR and the compiled binary. + ## (We cannot use the defaults in ${FSL_FIXDIR}/settings.sh here, since functionhighpassandvariancenormalize.m is an HCP + ## specific script that is not part of the FIX distribution). + + # For INTERPRETED MODES, make sure that matlab/octave has access to the functions it needs. + # normalise.m (needed by functionhighpassandvariancenormalize.m) is in '${HCPPIPEDIR}/global/matlab' + # Note that the call_matlab.sh script itself adds ${FSL_FIXDIR} (all the fix-related functions) + # and ${FSL_MATLAB_PATH} (e.g., read_avw.m, save_avw.m) to the matlab path. + # Note that fix_3_clean.m *appends* ${FSL_FIX_CIFTIRW} (defined in ${FSL_FIXDIR}/settings.sh) + # to the matlab path (i.e., the ciftiopen.m, ciftisave.m functions). + # BUT the CIFTI I/O functions are also now included within '${HCPPIPEDIR}/global/matlab' as well. + # AND since 'addpath' *prepends* to the matlab path, the versions in '${HCPPIPEDIR}/global/matlab' will therefore take precedence! + ML_PATHS="addpath('${HCPPIPEDIR}/global/matlab'); addpath('${this_script_dir}/scripts');" + + case ${FSL_FIX_MATLAB_MODE} in + 0) + # Use compiled Matlab + if [ -z "${MATLAB_COMPILER_RUNTIME}" ]; then + log_Err_Abort "To use MATLAB run mode: ${FSL_FIX_MATLAB_MODE}, the MATLAB_COMPILER_RUNTIME environment variable must be set" else - log_Msg "Found ${fmri}_hp${hp} is already filtered" + log_Msg "MATLAB_COMPILER_RUNTIME: ${MATLAB_COMPILER_RUNTIME}" fi - fi - if [ -e ${fmri}_Atlas_hp${hp}.dtseries.nii ] ; then - log_Msg "Found and using ${fmri}_Atlas_hp${hp}.dtseries.nii" - fi - - log_Msg "${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken $tr $hp $fmri ${FSL_FIX_WBC} $ndhpvol $ndhpcifti $ndcifti" - #${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken $tr $hp $fmri ${FSL_FIX_WBC} - ${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken $tr $hp $fmri ${FSL_FIX_WBC} $ndhpvol $ndhpcifti $ndcifti - - #log_Msg "matlab -nodesktop -nosplash -r \"addpath('$FSL_FIXDIR'); addpath('${FSLDIR}/etc/matlab'); functionhighpassandvariancenormalize_riken($tr,$hp,'$fmri','$FSL_FIX_WBC',$ndhpvol,$ndhpcifti,$ndcifti);\"" - #matlab -nodesktop -nosplash -r "addpath('$FSL_FIXDIR'); addpath('${FSLDIR}/etc/matlab'); functionhighpassandvariancenormalize_riken($tr,$hp,'$fmri','$FSL_FIX_WBC',$ndhpvol,$ndhpcifti,$ndcifti);" >> .fix.functionhighpassandvariancenormalize_riken.log 2>&1 - - #cd .. - if [ "$(cat ${fmri}_dims.txt)" != "" ] ; then - log_Msg "Dims: $(cat ${fmri}_dims.txt)" - else - log_Msg "Dims was not properly estimated by IcaDim" - fi + CSUBDIR=${HCPPIPEDIR}/ICAFIX/scripts/Compiled_functionhighpassandvariancenormalize + matlab_cmd="${FSL_FIXDIR}/call_matlab.sh -c ${MATLAB_COMPILER_RUNTIME} -b ${CSUBDIR} -f functionhighpassandvariancenormalize ${tr} ${hp} ${fmri} ${Caret7_Command} ${ndhpvol} ${ndhpcifti} ${ndcifti}" + log_Msg "Run compiled MATLAB with command..." + log_Msg "${matlab_cmd}" + ${matlab_cmd} + ;; + 1 | 2) + # Use interpreted MATLAB or Octave + # ${hp} needs to be passed in as a string, to handle the hp=pd* case + matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}, '${ndhpcifti}', '${ndcifti});" + log_Msg "Run interpreted MATLAB/Octave with command..." + log_Msg "${FSL_FIXDIR}/call_matlab.sh -r ${matlab_cmd}" + ${FSL_FIXDIR}/call_matlab.sh -r "${matlab_cmd}" # Note: the semicolon within ${matlab_cmd} needs to be protected with quotation marks + ;; + *) + # Unsupported MATLAB run mode + log_Err_Abort "Unsupported MATLAB run mode value (FSL_FIX_MATLAB_MODE) in $FSL_FIXDIR/settings.sh: ${FSL_FIX_MATLAB_MODE}" + ;; + esac + + log_Msg "Dims: $(cat ${fmri}_dims.txt)" # Demean the movement regressors (in the 'fake-NIFTI' format returned by functionhighpassandvariancenormalize) fslmaths ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf -Tmean ${fmri}_hp${hp}.ica/mc/prefiltered_func_data_mcf_conf_mean @@ -503,11 +505,6 @@ for fmri in $fmris ; do cd ${DIR} # Return to directory where script was launched - fmri=${fmri}_hp$hp - #cd ${fmri}.ica - # Per https://github.com/Washington-University/Pipelines/issues/60, the following line doesn't appear to be necessary - #$FSLDIR/bin/imln ../$fmri filtered_func_data - #cd .. log_Msg "Bottom of loop through fmris: fmri: ${fmri}" done ### END LOOP (for fmri in $fmris; do) From 5ec56d775348a7de84f51ce8bb5e17c976a602af Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 15:48:19 +0900 Subject: [PATCH 19/29] Modified to the same as HCP --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index 8ad0654..e6fe882 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -495,7 +495,7 @@ for fmri in $fmris ; do log_Err_Abort "Unsupported MATLAB run mode value (FSL_FIX_MATLAB_MODE) in $FSL_FIXDIR/settings.sh: ${FSL_FIX_MATLAB_MODE}" ;; esac - + log_Msg "Dims: $(cat ${fmri}_dims.txt)" # Demean the movement regressors (in the 'fake-NIFTI' format returned by functionhighpassandvariancenormalize) @@ -589,6 +589,7 @@ done # It is also not variance normalized (the mean VN map *across runs* was "restored" above). # But is detrended and within run demeaned and variance normalized for both volume and CIFTI concatfmri=`basename ${ConcatNameNoExt}` # Directory path is now removed +concat_fmri_orig=`basename $(remove_ext ${ConcatName})` concatfmrihp=${concatfmri}_hp${hp} #MPH: Need to be in original directory here, for the ${MovementNIFTI*STRING} variables to work @@ -600,8 +601,8 @@ mkdir ${ConcatFolder}/${concatfmrihp}.ica/mc fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp ${MovementNIFTIhpMergeSTRING} $tr fslmerge -tr ${ConcatFolder}/${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf ${MovementNIFTIMergeSTRING} $tr -cd `dirname $fmri` -} +#MPH: Now we can switch to ${ConcatFolder} +cd ${ConcatFolder} ## "2nd pass" VN ## This is solely so that we can apply our own variance normalization as part of the spatial ICA, @@ -756,8 +757,6 @@ cd ../.. ## Actually run FIX ## --------------------------------------------------------------------------- -AlreadyHP="-1" #Don't run highpass on concatenated data -cd ${ConcatFolder} log_Msg "Running FIX" # Changes to handle user specified training data file @@ -774,7 +773,6 @@ if [ "X${TrainingData}" != X ]; then if [ ! -f "${TrainingData}" ]; then TrainingData=${FSL_FIXDIR}/training_files/${TrainingData} fi - log_Msg "User has specified a training data file: ${TrainingData}" # finally, if the TrainingData file is not found, report an error and get out of here if [ ! -f "${TrainingData}" ]; then @@ -790,7 +788,16 @@ else #TSC: if you WANT to use, say, hp150 with hp2000 training, it should be explicitly requested by the user #TSC: exception: we will recommend polynomial detrending for performance, and hp 2000 is already trained, and largely indistinguishable log_Msg "training data file not specified" - training_hp="2000" + training_hp="${hp}" + if [[ ! -f "${FSL_FIXDIR}/training_files/HCP_hp${hp}.RData" ]]; then + if [[ "${pdFlag}" == "TRUE" && "${hp}" != "pd0" ]]; then # "hp=pd0" would be a demeaning only, so doesn't qualify here + #hack: hp 2000 is close enough to a polynomial detrend, but the latter is far faster + #we know the hp 2000 training data exists, so use it if we don't see an explicit hppd# set + training_hp=2000 + else + log_Err_Abort "no standard training data found for specified high pass (${hp}), please specify training data manually or use a standard high pass setting" + fi + fi TrainingData=${FSL_FIXDIR}/training_files/HCP_hp${training_hp}.RData fi log_Msg "using training data file: ${TrainingData}" @@ -879,11 +886,6 @@ for fmri in $fmris ; do Start=`echo "${Start} + ${NumTPS}" | bc -l` done -log_Msg "Finished" - - - - cd ${DIR} log_Msg "Completed!" From 76ecf44d52d6199442602453c18e811a73adf02b Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 16:53:55 +0900 Subject: [PATCH 20/29] fix matlab run --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 107 +++++++++++++++------------ 1 file changed, 61 insertions(+), 46 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index e6fe882..de30e7b 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -477,11 +477,12 @@ for fmri in $fmris ; do fi CSUBDIR=${HCPPIPEDIR}/ICAFIX/scripts/Compiled_functionhighpassandvariancenormalize - matlab_cmd="${FSL_FIXDIR}/call_matlab.sh -c ${MATLAB_COMPILER_RUNTIME} -b ${CSUBDIR} -f functionhighpassandvariancenormalize ${tr} ${hp} ${fmri} ${Caret7_Command} ${ndhpvol} ${ndhpcifti} ${ndcifti}" + matlab_cmd="${FSL_FIXDIR}/call_matlab.sh -c ${MATLAB_COMPILER_RUNTIME} -b ${CSUBDIR} -f functionhighpassandvariancenormalize ${tr} ${hp} ${fmri} ${Caret7_Command}" log_Msg "Run compiled MATLAB with command..." log_Msg "${matlab_cmd}" ${matlab_cmd} ;; + 1 | 2) # Use interpreted MATLAB or Octave # ${hp} needs to be passed in as a string, to handle the hp=pd* case @@ -611,17 +612,61 @@ cd ${ConcatFolder} ## See comments above about regarding the use of call_matlab here. AlreadyHP="-1" #Don't run highpass on concatenated data (since hp was done on the individual runs) -cd ${ConcatFolder} -if [ "$tr" = "" ] ; then - fmri=`echo ${fmris} | tr ' ' '\n' | tail -1` - tr=`$FSLDIR/bin/fslval $fmri pixdim4` - log_Msg "tr: $tr" + +case ${FSL_FIX_MATLAB_MODE} in + 0) + # Use compiled Matlab + CSUBDIR=${HCPPIPEDIR}/ICAFIX/scripts/Compiled_functionhighpassandvariancenormalize + matlab_cmd="${FSL_FIXDIR}/call_matlab.sh -c ${MATLAB_COMPILER_RUNTIME} -b ${CSUBDIR} -f functionhighpassandvariancenormalize ${tr} ${AlreadyHP} ${concatfmrihp} ${Caret7_Command}" + log_Msg "Run compiled MATLAB with command..." + log_Msg "${matlab_cmd}" + ${matlab_cmd} + ;; + + 1 | 2) + # Use interpreted MATLAB or Octave + # ${hp} needs to be passed in as a string, to handle the hp=pd* case + matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}, '${ndhpcifti}', '${ndcifti});" + log_Msg "Run interpreted MATLAB/Octave with command..." + log_Msg "${FSL_FIXDIR}/call_matlab.sh -r ${matlab_cmd}" + ${FSL_FIXDIR}/call_matlab.sh -r "${matlab_cmd}" # Note: the semicolon within ${matlab_cmd} needs to be protected with quotation marks + ;; + *) + # Unsupported MATLAB run mode + log_Err_Abort "Unsupported MATLAB run mode value (FSL_FIX_MATLAB_MODE) in $FSL_FIXDIR/settings.sh: ${FSL_FIX_MATLAB_MODE}" + ;; +esac + +# Takuya Hayashi inserted, July 2018 +if [ -e ${concatfmri}_WF.txt ] ; then rm ${concatfmri}_WF.txt;fi # Add by TH Aug 2019 +echo $ndhpvol $ndhpcifti $ndcifti > ${concatfmri}_WF.txt + +if [ "$UserDim" != "NONE" ] ; then + Dim=$UserDim + log_Msg "Melodic use Dim=${Dim} user specified" +else + Dim=`cat ${concatfmrihp}_dims.txt` + + log_Msg "Dim: ${Dim}" fi -log_Msg "running MELODIC" -log_Msg "About to run melodic: Contents of ${concatfmri}.ica follow" -if [ "${DEBUG}" = "TRUE" ] ; then - ls -lRa `remove_ext ${concatfmri}`.ica +## --------------------------------------------------------------------------- +## Run melodic on concatenated file +## --------------------------------------------------------------------------- + +log_Debug_Msg "About to run melodic: Contents of ${concatfmrihp}.ica follow" +if [ ! -z "${log_debugOn}" ] ; then + ls -lRa ${concatfmrihp}.ica +fi + +#grab melodic from $PATH by default, don't hardcode it with respect to $FSLDIR +#we need to do "if which ..." because the script currently uses ERR trap +if which melodic &> /dev/null +then + MELODIC=$(which melodic 2> /dev/null) +else + #if it isn't even in $PATH, fall back on FSLDIR + MELODIC="${FSLDIR}/bin/melodic" fi log_Msg "Running MELODIC located at: $MELODIC" @@ -636,43 +681,14 @@ if [ ! -z "${log_debugOn}" ] ; then fi log_Debug_Msg "End of melodic version log, help, and checksum" -###Run on concatenated file -- - -# Takuya Hayashi inserted, July 2018 - #${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize.log -f functionhighpassandvariancenormalize $tr $AlreadyHP $concatfmri ${FSL_FIX_WBC} - -if [ -e ${concatfmri}_WF.txt ] ; then rm ${concatfmri}_WF.txt;fi # Add by TH Aug 2019 -echo $ndhpvol $ndhpcifti $ndcifti > ${concatfmri}_WF.txt - -log_Msg "${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken ${tr} ${AlreadyHP} ${concatfmri} ${FSL_FIX_WBC} ${ndhpvol} ${ndhpcifti} ${ndcifti}" - - ${FSL_FIXDIR}/call_matlab.sh -l .fix.functionhighpassandvariancenormalize_riken.log -f functionhighpassandvariancenormalize_riken ${tr} ${AlreadyHP} ${concatfmri} ${FSL_FIX_WBC} ${ndhpvol} ${ndhpcifti} ${ndcifti} - -if [ "$UserDim" != "NONE" ] ; then - Dim=$UserDim - log_Msg "Melodic use Dim=${Dim} user specified" -else - Dim=`cat ${concatfmri}_dims.txt` - log_Msg "Melodic use Dim=${Dim} estimated by IcaDim" +# Use the Wishart-based dim estimate instead of the melodic-based estimate (--dim=${Dim}) +# Also turn off melodic VN and feed in the input that has already been VN'ed instead (--vn flag, which turns OFF melodic's vn) +# Added "-m ${concatfmri}_brain_mask" to avoid memory issue - Takuya Hayahsi July 2018 +melodic_cmd="${MELODIC} -i ${concatfmrihp}_vnts -o ${concatfmrihp}.ica/filtered_func_data.ica --nobet --report --Oall --tr=${tr} --vn --dim=${Dim} -m ${concatfmri}_brain_mask" +if [ ! -z "${log_debugOn}" ] ; then + melodic_cmd="${melodic_cmd} --verbose --debug" fi -#Need to decide if using Wishart-based dim estimate or melodic-based estimate. Also need to decide if turning off melodic VN and feeding in ours --vn --dim=${Dim} - -# Original melodic command without logging of command -# ${FSLDIR}/bin/melodic -i $concatfmri -o ${concatfmri}.ica/filtered_func_data.ica -d -250 --nobet --report --Oall --tr=$tr -# -# melodic command with verbose output and debugging output, used when trying to determine why melodic was crashing -# melodic_cmd="${FSLDIR}/bin/melodic -i $concatfmri -o ${concatfmri}.ica/filtered_func_data.ica --nobet --report --Oall --tr=$tr --verbose --debug" - -#melodic_cmd="${FSLDIR}/bin/melodic -i $concatfmri -o ${concatfmri}.ica/filtered_func_data.ica --nobet --report --Oall --tr=$tr" - -MELODIC="${FSLDIR}/bin/melodic" -MELODIC=/media/myelin/brainmappers/Connectome_Project/MultiRunICA+FIXTest/melodic -MELODIC=/mnt/pub/devel/FSL-devel/MELODIC/melodic_20170712 - -# Added "-m `remove_ext ${ConcatName}`_brain_mask" to avoid memory issue - Takuya Hayahsi July 2018 -melodic_cmd="${MELODIC} -i ${concatfmri}_vnts -o ${concatfmri}.ica/filtered_func_data.ica --nobet --report --Oall --tr=$tr --vn --dim=${Dim} --verbose --debug -m `remove_ext ${ConcatName}`_brain_mask" - log_Msg "melodic_cmd: ${melodic_cmd}" ${melodic_cmd} return_code=$? @@ -700,7 +716,6 @@ $FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data.ica/concat_data cd ${concatfmrihp}.ica -$FSLDIR/bin/imln ../${concatfmri} filtered_func_data # This is the concated volume time series from the 1st pass VN, with requested # hp-filtering applied and with the mean VN map multiplied back in $FSLDIR/bin/imln ../${concatfmrihp} filtered_func_data From 358f6b395ac60138f9083af35d5dbed3b8622db6 Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 30 Jun 2022 17:18:46 +0900 Subject: [PATCH 21/29] fix PostFix --- ICAFIX/hcp_fix_multi_run_RIKEN_merge | 55 +++++++++++++++++++++++----- 1 file changed, 45 insertions(+), 10 deletions(-) diff --git a/ICAFIX/hcp_fix_multi_run_RIKEN_merge b/ICAFIX/hcp_fix_multi_run_RIKEN_merge index de30e7b..afdae20 100755 --- a/ICAFIX/hcp_fix_multi_run_RIKEN_merge +++ b/ICAFIX/hcp_fix_multi_run_RIKEN_merge @@ -842,20 +842,55 @@ fi log_Msg "Done running FIX" +## --------------------------------------------------------------------------- +## Rename some files (relative to the default names coded in fix_3_clean) +## --------------------------------------------------------------------------- +$FSLDIR/bin/immv ${concatfmrihp}.ica/filtered_func_data_clean ${concatfmrihp}_clean +if [ "$?" -ne "0" ]; then + log_Err_Abort "Something went wrong; ${concatfmrihp}.ica/filtered_func_data_clean wasn't created" +fi +if [ -f ${concatfmrihp}.ica/Atlas_clean.dtseries.nii ]; then + /bin/mv ${concatfmrihp}.ica/Atlas_clean.dtseries.nii ${concatfmri}_Atlas_hp${hp}_clean.dtseries.nii +else + log_Err_Abort "Something went wrong; ${concatfmrihp}.ica/Atlas_clean.dtseries.nii wasn't created" +fi -cd ${ConcatFolder} +# The variance normalization ("_vn") outputs of fix (fix_3_clean) require use of fix1.067 or later +# So check whether those files exist before moving/renaming them +if [ `$FSLDIR/bin/imtest ${concatfmrihp}.ica/filtered_func_data_clean_vn` = 1 ]; then + $FSLDIR/bin/immv ${concatfmrihp}.ica/filtered_func_data_clean_vn ${concatfmrihp}_clean_vn +fi +if [ -f ${concatfmrihp}.ica/Atlas_clean_vn.dscalar.nii ]; then + /bin/mv ${concatfmrihp}.ica/Atlas_clean_vn.dscalar.nii ${concatfmri}_Atlas_hp${hp}_clean_vn.dscalar.nii + if [ -f ${concatfmri}_Atlas_hp${hp}_clean.dtseries.nii ] ; then + ${FSL_FIX_WBC} -cifti-math "TCS/VN" ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dtseries.nii -var TCS ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii -var VN ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii -select 1 1 -repeat + fi +fi +log_Msg "Done renaming files" + +# Remove the 'fake-NIFTI' file created in fix_3_clean for high-pass filtering of the CIFTI (if it exists) +$FSLDIR/bin/imrm ${concatfmrihp}.ica/Atlas -log_Msg "Moving ${concatfmri}.ica/filtered_func_data_clean to ${concatfmri}_clean" -$FSLDIR/bin/immv ${concatfmri}.ica/filtered_func_data_clean ${concatfmri}_clean -$FSLDIR/bin/immv ${concatfmri}.ica/filtered_func_data_clean_vn ${concatfmri}_clean_vn +# Always delete things with too-generic names +$FSLDIR/bin/imrm ${concatfmrihp}.ica/filtered_func_data +rm -f ${concatfmrihp}.ica/Atlas.dtseries.nii + +# Optional deletion of highpass intermediates +# (hp<0 not supported in this script currently, so no need to condition on value of hp) +if [ "${DeleteIntermediates}" == "TRUE" ] +then + $FSLDIR/bin/imrm ${concatfmri} ${concatfmrihp} + rm -f ${concatfmri}_Atlas.dtseries.nii ${concatfmri}_Atlas_hp${hp}.dtseries.nii +fi -log_Msg "Checking for existence of ${concatfmri}.ica/Atlas_clean.dtseries.nii" -if [ -f ${concatfmri}.ica/Atlas_clean.dtseries.nii ] ; then - log_Msg "Moving ${concatfmri}.ica/Atlas_clean.dtseries.nii to ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii" - /bin/mv ${concatfmri}.ica/Atlas_clean.dtseries.nii ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii - /bin/mv ${concatfmri}.ica/Atlas_clean_vn.dscalar.nii ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii - ${FSL_FIX_WBC} -cifti-math "TCS/VN" ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dtseries.nii -var TCS ${concat_fmri_orig}_Atlas_hp${hp}_clean.dtseries.nii -var VN ${concat_fmri_orig}_Atlas_hp${hp}_clean_vn.dscalar.nii -select 1 1 -repeat +# Convert sICA+FIX cleaned movement regressors to text +if [ -f ${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp_clean.nii.gz ] ; then + fslmeants -i ${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp_clean.nii.gz -o Movement_Regressors_hp${hp}_clean.txt --showall + # Strip header lines included as part of '--showall' flag + nVols=$($FSLDIR/bin/fslnvols ${concatfmrihp}.ica/mc/prefiltered_func_data_mcf_conf_hp_clean.nii.gz) + # Execute tail in a subshell, so we can successfully overwrite file with same name + echo "$(tail -n ${nVols} Movement_Regressors_hp${hp}_clean.txt)" > Movement_Regressors_hp${hp}_clean.txt fi #Cleaned volume and CIFTI have no mean and are intensity normalized by multiplying in the overall variance normalization From 2c2d516be88baead28539d479c820850073bf80b Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 6 Jul 2022 14:17:01 +0900 Subject: [PATCH 22/29] rename param ndvol to ndconcatvol --- ICAFIX/ReApplyFixMultiRunPipeline.sh | 6 +++--- .../ReApplyFixMultiRunPipeline_RIKEN_merge.sh | 6 +++--- ...nhighpassandvariancenormalize_riken_merge.m | 18 +++++++++--------- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/ICAFIX/ReApplyFixMultiRunPipeline.sh b/ICAFIX/ReApplyFixMultiRunPipeline.sh index 6fe33d5..80fd8d1 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline.sh @@ -564,11 +564,11 @@ main() if [ -e ${ConcatNameNoExt}_hp${hp}_wf.txt ] ; then ndhpvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $1}'`" ndhpcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $2}'`" - ndcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" + ndconcatvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" elif [ -z "$p_WF" ] ; then ndhpvol="`echo $p_WF | cut -d ',' -f1`" ndhpcifti="`echo $p_WF | cut -d ',' -f2`" - ndcifti="`echo $p_WF | cut -d ',' -f3`" + ndconcatvol="`echo $p_WF | cut -d ',' -f3`" else log_Err_Abort "ERROR: cannot find Ndist used in MR-FIX. Please use option --wf to set Ndist" fi @@ -715,7 +715,7 @@ main() fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndcifti}', '${RegString}');" + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}', '${RegString}');" log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" diff --git a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh index 0f06d54..0e9ac9f 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh @@ -564,11 +564,11 @@ main() if [ -e ${ConcatNameNoExt}_hp${hp}_wf.txt ] ; then ndhpvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $1}'`" ndhpcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $2}'`" - ndcifti="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" + ndconcatvol="`cat ${ConcatNameNoExt}_hp${hp}_wf.txt | awk '{print $3}'`" elif [ -z "$p_WF" ] ; then ndhpvol="`echo $p_WF | cut -d ',' -f1`" ndhpcifti="`echo $p_WF | cut -d ',' -f2`" - ndcifti="`echo $p_WF | cut -d ',' -f3`" + ndconcatvol="`echo $p_WF | cut -d ',' -f3`" else log_Err_Abort "ERROR: cannot find Ndist used in MR-FIX. Please use option --wf to set Ndist" fi @@ -715,7 +715,7 @@ fi fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndcifti}', '${RegString}');" + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}', '${RegString}');" log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m index 899fa2d..638acba 100755 --- a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -24,7 +24,7 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) % WBC: wb_command (full path) % varargin: % [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. -% [ndhpvol ndhpcifti ndvol]: number of distributions highpassed volume, highpassed cifti and concatenated volume. OPTIONAL. +% [ndhpvol ndhpcifti ndconcatvol]: number of distributions highpassed volume, highpassed cifti and concatenated volume. OPTIONAL. % Note: HP='pd0' would be interpreted as a true 0th order detrend, which is % the same as demeaning. Mathematically, this is the same as the HP<0 condition, @@ -53,7 +53,7 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) regstring = varargin{1}; %this has the underscore on the front already ndhpvol=2 ndhpcifti=3 - ndvol=3 + ndconcatvol=3 if ~ischar(regstring) error('%s: REGSTRING should be a string', mfilename); end @@ -67,22 +67,22 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) end ndhpvol=varargin{1} ndhpcifti=varargin{2} - ndvol=varargin{3} + ndconcatvol=varargin{3} if ~isint(ndhpvol) error('%s: NDHPVOL should be an integer', mfilename); end if ~isint(ndhpcifti) error('%s: NDHPCIFTI should be an integer', mfilename); end - if ~isint(ndvol) - error('%s: NDVOL should be an integer', mfilename); + if ~isint(ndconcatvol) + error('%s: NDCONCATVOL should be an integer', mfilename); end end %UNTITLED4 Summary of this function goes here % Detailed explanation goes here -% Default value: ndhpvol=2;ndhpcifti=3; ndvol=3; - Takuya Hayashi, Matt Glasser -fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndvol) +% Default value: ndhpvol=2;ndhpcifti=3; ndconcatvol=3; - Takuya Hayashi, Matt Glasser +fprintf('hp=%d, ndist=%d,%d,%d',hp,ndhpvol,ndhpcifti,ndconcatvol) % Check whether polynomial detrend is being requested for the high-pass filtering. @@ -228,7 +228,7 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping else if dovol > 0 - Outcts=icaDim(cts,0,1,-1,ndvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform + Outcts=icaDim(cts,0,1,-1,ndconcatvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform end end @@ -285,7 +285,7 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) dlmwrite([fmri '_dims.txt'],[Outcts.calcDim],'\t'); end -dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndvol],'\t'); +dlmwrite([fmri '_wf.txt'],[ndhpvol ndhpcifti ndconcatvol],'\t'); end From 561382a937c1db8fc792b60aa098d5b4d230c4a4 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 6 Jul 2022 15:00:02 +0900 Subject: [PATCH 23/29] fix varargin, add docifti to varargin --- ...highpassandvariancenormalize_riken_merge.m | 39 +++++++++++++++---- 1 file changed, 31 insertions(+), 8 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m index 638acba..9ee0191 100755 --- a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -1,6 +1,6 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) % -% FUNCTIONHIGHPASSANDVARIANCENORMALIZE(TR,HP,FMRI,WBC,[REGSTRING]) +% FUNCTIONHIGHPASSANDVARIANCENORMALIZE(TR,HP,FMRI,WBC,[REGSTRING],[ndhpvol ndhpcifti ndconcatvol]) % % This function: % (1) Detrends / high-pass filters motion confounds, NIFTI (volume) and CIFTI files @@ -25,6 +25,8 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) % varargin: % [REGSTRING]: Additional registration-related string to add to output file names. OPTIONAL. % [ndhpvol ndhpcifti ndconcatvol]: number of distributions highpassed volume, highpassed cifti and concatenated volume. OPTIONAL. +% [docifti]: Calculate cifti. If you use this option, the [REGSTRING] and [ndhpvol ndhpcifti ndconcatvol] options must be specified. +% If "" is specified for the [REGSTRING] and [ndhpvol ndhpcifti ndconcatvol] options, the default values are used. % Note: HP='pd0' would be interpreted as a true 0th order detrend, which is % the same as demeaning. Mathematically, this is the same as the HP<0 condition, @@ -42,6 +44,7 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) %% Defaults dovol = 1; +docifti = 1; regstring = ''; pdflag = false; % Polynomial detrend pdstring = 'pd'; % Expected string at start of HP variable to request a "polynomial detrend" @@ -58,16 +61,33 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) error('%s: REGSTRING should be a string', mfilename); end elseif length(varargin) >= 3 - if ~isempty(varagin{4}) - dovol = 0 - regstring = varargin{4} + if ~isempty(varargin{4}) + if ~isempty(varargin{1}) + dovol = 0 + end + regstring = varargin{1} if ~ischar(regstring) error('%s: REGSTRING should be a string', mfilename); end + ndhpvol=varargin{2} + ndhpcifti=varargin{3} + ndconcatvol=varargin{4} + else + ndhpvol=varargin{1} + ndhpcifti=varargin{2} + ndconcatvol=varargin{3} + if isempty(ndvol) + ndvol=2 + end + if isempty(ndvol) + ndhpcifti=3 + end + if isempty(ndconcatvol) + ndvol=3 + end + if ~isempty(varargin{5}) + docifti=varargin{5} end - ndhpvol=varargin{1} - ndhpcifti=varargin{2} - ndconcatvol=varargin{3} if ~isint(ndhpvol) error('%s: NDHPVOL should be an integer', mfilename); end @@ -75,7 +95,10 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) error('%s: NDHPCIFTI should be an integer', mfilename); end if ~isint(ndconcatvol) - error('%s: NDCONCATVOL should be an integer', mfilename); + error('%s: NDCONCATVOL should be an integer', mfilename); + end + if ~isint(docifti) + error('%s: DOCIFTI should be an integer', mfilename); end end From c0d1ffd1be46e4089e5df55dfe8034c67f0cd9ad Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 6 Jul 2022 15:02:49 +0900 Subject: [PATCH 24/29] fix positional arg RegString --- ICAFIX/ReApplyFixMultiRunPipeline.sh | 2 +- ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ICAFIX/ReApplyFixMultiRunPipeline.sh b/ICAFIX/ReApplyFixMultiRunPipeline.sh index 80fd8d1..887fb9d 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline.sh @@ -715,7 +715,7 @@ main() fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}', '${RegString}');" + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}','${RegString}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}');" log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" diff --git a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh index 0e9ac9f..c3639fb 100755 --- a/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh +++ b/ICAFIX/ReApplyFixMultiRunPipeline_RIKEN_merge.sh @@ -715,7 +715,7 @@ fi fi # ${hp} needs to be passed in as a string, to handle the hp=pd* case - local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}', '${RegString}');" + local matlab_cmd="${ML_PATHS} functionhighpassandvariancenormalize(${tr}, '${hp}', '${fmri}', '${Caret7_Command}','${RegString}', '${ndhpvol}', '${ndhpcifti}', '${ndconcatvol}');" log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." log_Msg "${matlab_cmd}" From 97e4212ae936eda167715001df750ede454408eb Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 6 Jul 2022 15:36:11 +0900 Subject: [PATCH 25/29] Add docifti conditional branch --- ...highpassandvariancenormalize_riken_merge.m | 76 ++++++++++--------- 1 file changed, 42 insertions(+), 34 deletions(-) diff --git a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m index 9ee0191..3c0f78b 100755 --- a/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m +++ b/ICAFIX/functionhighpassandvariancenormalize_riken_merge.m @@ -54,39 +54,39 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) if length(varargin) == 1 && ~isempty(varargin{1}) dovol = 0; %regname is only used for a new surface registration, so will never require redoing volume regstring = varargin{1}; %this has the underscore on the front already - ndhpvol=2 - ndhpcifti=3 - ndconcatvol=3 + ndhpvol=2; + ndhpcifti=3; + ndconcatvol=3; if ~ischar(regstring) error('%s: REGSTRING should be a string', mfilename); end elseif length(varargin) >= 3 if ~isempty(varargin{4}) if ~isempty(varargin{1}) - dovol = 0 + dovol = 0; end - regstring = varargin{1} + regstring = varargin{1}; if ~ischar(regstring) error('%s: REGSTRING should be a string', mfilename); end - ndhpvol=varargin{2} - ndhpcifti=varargin{3} - ndconcatvol=varargin{4} + ndhpvol=varargin{2}; + ndhpcifti=varargin{3}; + ndconcatvol=varargin{4}; else - ndhpvol=varargin{1} - ndhpcifti=varargin{2} - ndconcatvol=varargin{3} + ndhpvol=varargin{1}; + ndhpcifti=varargin{2}; + ndconcatvol=varargin{3}; if isempty(ndvol) - ndvol=2 + ndvol=2; end if isempty(ndvol) - ndhpcifti=3 + ndhpcifti=3; end if isempty(ndconcatvol) - ndvol=3 + ndvol=3; end if ~isempty(varargin{5}) - docifti=varargin{5} + docifti=varargin{5}; end if ~isint(ndhpvol) error('%s: NDHPVOL should be an integer', mfilename); @@ -154,8 +154,9 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) %% normalise function is in HCPPIPEDIR/global/matlab/normalise.m confounds=normalise([confounds [zeros(1,size(confounds,2)); confounds(2:end,:)-confounds(1:end-1,:)] ]); confounds=normalise([confounds confounds.*confounds]); - - BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); + if docifti > 0 + BO=ciftiopen([fmri '_Atlas' regstring '.dtseries.nii'],WBC); + end end %% Apply hp filtering of the motion confounds, volume (if requested), and CIFTI @@ -188,9 +189,10 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) % Use -d flag in fslcpgeom (even if not technically necessary) to reduce possibility of mistakes when editing script call_fsl(['fslcpgeom ' fmri '.nii.gz ' fmri hpstring '.nii.gz -d']); end - - BO.cdata=detrendpoly(BO.cdata',hp)'; - ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + if docifti > 0 + BO.cdata=detrendpoly(BO.cdata',hp)'; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + end elseif hp>0 % "fslmaths -bptf" based filtering if dovol > 0 @@ -212,13 +214,15 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) clear ctsfull; end - BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); - save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); - call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); - grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; - ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); - delete('Atlas.nii.gz'); + if docifti > 0 + BOdimX=size(BO.cdata,1); BOdimZnew=ceil(BOdimX/100); BOdimT=size(BO.cdata,2); + save_avw(reshape([BO.cdata ; zeros(100*BOdimZnew-BOdimX,BOdimT)],10,10,BOdimZnew,BOdimT),'Atlas','f',[1 1 1 TR]); + call_fsl(sprintf('fslmaths Atlas -bptf %f -1 Atlas',0.5*hp/TR)); + grot=reshape(single(read_avw('Atlas')),100*BOdimZnew,BOdimT); BO.cdata=grot(1:BOdimX,:); clear grot; + ciftisave(BO,[fmri '_Atlas' regstring hpstring '.dtseries.nii'],WBC); + end + delete('Atlas.nii.gz'); elseif hp<0 % If no hp filtering, still need to at least demean the volumetric time series if dovol > 0 @@ -247,8 +251,9 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) if dovol > 0 Outcts=icaDim(cts,0,1,-1,ndhpvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform end - - OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping + if docifti > 0 + OutBO=icaDim(BO.cdata,0,1,-1,ndhpcifti); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, CIFTI fits three distributions to deal with volume to CIFTI mapping + end else if dovol > 0 Outcts=icaDim(cts,0,1,-1,ndconcatvol); %0=Don't detrend, 1=Initialize variance normalization at 1, -1=Converge with running dim average, Volume fits two distributions to deal with MNI transform @@ -267,12 +272,15 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) clear vnfull; call_fsl(['fslcpgeom ' fmri '_mean.nii.gz ' fname ' -d']); end - - VN=BO; - VN.cdata=OutBO.noise_unst_std; + if docifti > 0 + VN=BO; + VN.cdata=OutBO.noise_unst_std; + end fname=[fmri '_Atlas' regstring hpstring '_vn.dscalar.nii']; disp(['saving ' fname]); - ciftisavereset(VN,fname,WBC); + if docifti > 0 + ciftisavereset(VN,fname,WBC); + end end %% Apply VN and Save HP_VN TCS @@ -290,14 +298,14 @@ function functionhighpassandvariancenormalize(TR,hp,fmri,WBC,varargin) call_fsl(['fslcpgeom ' fmri '.nii.gz ' fname ' -d']); end % For CIFTI, we can use the extension to distinguish between VN maps (.dscalar) and VN'ed time series (.dtseries) -if hp>=0 +if hp>=0 && docifti BO.cdata=BO.cdata./repmat(OutBO.noise_unst_std,1,size(BO.cdata,2)); ciftisave(BO,[fmri '_Atlas' regstring hpstring '_vn.dtseries.nii'],WBC); end %Echo Dims %TSC: add the regstring to the output filename to avoid overwriting -if hp>=0 +if hp>=0 && docifti > 0 if dovol > 0 dlmwrite([fmri regstring '_dims.txt'],[Outcts.calcDim OutBO.calcDim],'\t'); else From 1528a2f61b53e4899a73f8e9427d2d74d39c3b5d Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 6 Jul 2022 16:08:25 +0900 Subject: [PATCH 26/29] add merge working file --- ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh | 647 +++++++++++++++++++++++ 1 file changed, 647 insertions(+) create mode 100755 ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh diff --git a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh new file mode 100755 index 0000000..e971f27 --- /dev/null +++ b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh @@ -0,0 +1,647 @@ +#!/bin/bash + +# +# # ReApplyFixPipeline.sh +# +# ## Copyright Notice +# +# Copyright (C) 2015-2017 The Human Connectome Project/Connectome Coordination Facility +# +# * Washington University in St. Louis +# * University of Minnesota +# * Oxford University +# +# ## Author(s) +# +# * Matthew F. Glasser, Department of Anatomy and Neurobiology, Washington University in St. Louis +# * Timothy B. Brown, Neuroinformatics Research Group, Washington University in St. Louis +# +# ## Product +# +# [Human Connectome Project][HCP] (HCP) Pipelines +# +# ## License +# +# See the [LICENSE](https://github.com/Washington-Univesity/Pipelines/blob/master/LICENSE.md) file +# +# +# [HCP]: http://www.humanconnectome.org +# + +# ------------------------------------------------------------------------------ +# Show usage information for this script +# ------------------------------------------------------------------------------ + +usage() +{ + local script_name + script_name=$(basename "${0}") + + cat < = user supplied value + + Note: The PARAMETERS can be specified positionally (i.e. without using the --param=value + form) by simply specifying all values on the command line in the order they are + listed below. + + e.g. ${script_name} /path/to/study/folder 100307 rfMRI_REST1_LR 2000 ... + + [--help] : show usage information and exit + --path= OR --study-folder= + --subject= + --fmri-name= (Do not include path, nifti extension, or the 'hp' string). + --high-pass= Number to represent the ${HighPass} variable used in ICA+FIX + [--reg-name= defaults to ${G_DEFAULT_REG_NAME}. (Use NONE for MSMSulc registration) + [--low-res-mesh=] defaults to ${G_DEFAULT_LOW_RES_MESH} + [--matlab-run-mode={0, 1, 2}] defaults to ${G_DEFAULT_MATLAB_RUN_MODE} + 0 = Use compiled MATLAB + 1 = Use interpreted MATLAB + 2 = Use interpreted Octave + [--motion-regression={TRUE, FALSE}] defaults to ${G_DEFAULT_MOTION_REGRESSION} + +EOF +exit 1; +} +[ "$1" = "" ] && usage +# ------------------------------------------------------------------------------ +# Get the command line options for this script. +# ------------------------------------------------------------------------------ + +get_options() +{ + local arguments=($@) + + # initialize global output variables + unset p_StudyFolder # ${1} + unset p_Subject # ${2} + unset p_fMRIName # ${3} + unset p_HighPass # ${4} + unset p_RegName # ${5} + unset p_LowResMesh # ${6} + unset p_MatlabRunMode # ${7} + unset p_MotionRegression # ${8} + + # set default values + p_RegName=${G_DEFAULT_REG_NAME} + p_LowResMesh=${G_DEFAULT_LOW_RES_MESH} + p_MatlabRunMode=${G_DEFAULT_MATLAB_RUN_MODE} + p_MotionRegression=${G_DEFAULT_MOTION_REGRESSION} + + # parse arguments + local num_args=${#arguments[@]} + local argument + local index=0 + + while [ "${index}" -lt "${num_args}" ]; do + argument=${arguments[index]} + + case ${argument} in + --help) + usage + exit 1 + ;; + --path=*) + p_StudyFolder=${argument#*=} + index=$(( index + 1 )) + ;; + --study-folder=*) + p_StudyFolder=${argument#*=} + index=$(( index + 1 )) + ;; + --subject=*) + p_Subject=${argument#*=} + index=$(( index + 1 )) + ;; + --fmri-name=*) + p_fMRIName=${argument#*=} + index=$(( index + 1 )) + ;; + --high-pass=*) + p_HighPass=${argument#*=} + index=$(( index + 1 )) + ;; + --reg-name=*) + p_RegName=${argument#*=} + index=$(( index + 1 )) + ;; + --low-res-mesh=*) + p_LowResMesh=${argument#*=} + index=$(( index + 1 )) + ;; + --matlab-run-mode=*) + p_MatlabRunMode=${argument#*=} + index=$(( index + 1 )) + ;; + --motion-regression=*) + p_MotionRegression=${argument#*=} + index=$(( index + 1 )) + ;; + *) + usage + log_Err_Abort "unrecognized option: ${argument}" + ;; + esac + done + + local error_count=0 + + # check required parameters + if [ -z "${p_StudyFolder}" ]; then + log_Err "Study Folder (--path= or --study-folder=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Study Folder: ${p_StudyFolder}" + fi + + if [ -z "${p_Subject}" ]; then + log_Err "Subject ID (--subject=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Subject ID: ${p_Subject}" + fi + + if [ -z "${p_fMRIName}" ]; then + log_Err "fMRI Name (--fmri-name=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "fMRI Name: ${p_fMRIName}" + fi + + if [ -z "${p_HighPass}" ]; then + log_Err "High Pass (--high-pass=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "High Pass: ${p_HighPass}" + fi + + if [ -z "${p_RegName}" ]; then + log_Err "Reg Name (--reg-name=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Reg Name: ${p_RegName}" + fi + + if [ -z "${p_LowResMesh}" ]; then + log_Err "Low Res Mesh (--low-res-mesh=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Low Res Mesh: ${p_LowResMesh}" + fi + + if [ -z "${p_MatlabRunMode}" ]; then + log_Err "MATLAB run mode value (--matlab-run-mode=) required" + error_count=$(( error_count + 1 )) + else + case ${p_MatlabRunMode} in + 0) + log_Msg "MATLAB Run Mode: ${p_MatlabRunMode} - Use compiled MATLAB" + ;; + 1) + log_Msg "MATLAB Run Mode: ${p_MatlabRunMode} - Use interpreted MATLAB" + ;; + 2) + log_Msg "MATLAB Run Mode: ${p_MatlabRunMode} - Use interpreted Octave" + ;; + *) + log_Err "MATLAB Run Mode value must be 0, 1, or 2" + error_count=$(( error_count + 1 )) + ;; + esac + fi + + if [ -z "${p_MotionRegression}" ]; then + log_Err "motion regression setting (--motion-regression=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Motion Regression: ${p_MotionRegression}" + fi + + if [ ${error_count} -gt 0 ]; then + log_Err_Abort "For usage information, use --help" + fi +} + +# ------------------------------------------------------------------------------ +# Show Tool Versions +# ------------------------------------------------------------------------------ + +show_tool_versions() +{ + # Show HCP pipelines version + log_Msg "Showing HCP Pipelines version" + cat ${HCPPIPEDIR}/version.txt + + # Show wb_command version + log_Msg "Showing Connectome Workbench (wb_command) version" + ${CARET7DIR}/wb_command -version + + # Show FSL version + log_Msg "Showing FSL version" + fsl_version_get fsl_ver + log_Msg "FSL version: ${fsl_ver}" +} + +# ------------------------------------------------------------------------------ +# Check for whether or not we have hand reclassification files +# ------------------------------------------------------------------------------ + +have_hand_reclassification() +{ + local StudyFolder="${1}" + local Subject="${2}" + local fMRIName="${3}" + local HighPass="${4}" + + if (( HighPass >= 0 )); then + [ -e "${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName}/${fMRIName}_hp${HighPass}.ica/HandNoise.txt" ] + else + [ -e "${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName}/${fMRIName}.ica/HandNoise.txt" ] + fi +} + +# ------------------------------------------------------------------------------ +# Main processing of script. +# ------------------------------------------------------------------------------ + +main() +{ + log_Msg "Starting main functionality" + + # Retrieve positional parameters + local StudyFolder="${1}" + local Subject="${2}" + local fMRIName="${3}" + local HighPass="${4}" + local RegName="${5}" + if [ -z "${5}" ]; then + RegName=${G_DEFAULT_REG_NAME} + else + RegName="${5}" + fi + + local LowResMesh + if [ -z "${6}" ]; then + LowResMesh=${G_DEFAULT_LOW_RES_MESH} + else + LowResMesh="${6}" + fi + + local MatlabRunMode + if [ -z "${7}" ]; then + MatlabRunMode=${G_DEFAULT_MATLAB_RUN_MODE} + else + MatlabRunMode="${7}" + fi + + local MotionRegression + if [ -z "${8}" ]; then + MotionRegression="${G_DEFAULT_MOTION_REGRESSION}" + else + MotionRegression="${8}" + fi + + # Turn MotionRegression into an appropriate numeric value for fix_3_clean + case $(echo ${MotionRegression} | tr '[:upper:]' '[:lower:]') in + ( true | yes | 1) + MotionRegression=1 + ;; + ( false | no | none | 0) + MotionRegression=0 + ;; + *) + log_Err_Abort "motion regression setting must be TRUE or FALSE" + ;; + esac + + # Log values retrieved from positional parameters + log_Msg "StudyFolder: ${StudyFolder}" + log_Msg "Subject: ${Subject}" + log_Msg "fMRIName: ${fMRIName}" + log_Msg "HighPass: ${HighPass}" + log_Msg "RegName: ${RegName}" + log_Msg "LowResMesh: ${LowResMesh}" + log_Msg "MatlabRunMode: ${MatlabRunMode}" + log_Msg "MotionRegression: ${MotionRegression}" + + # Naming Conventions and other variables + local Caret7_Command="${CARET7DIR}/wb_command" + log_Msg "Caret7_Command: ${Caret7_Command}" + + local RegString + if [ ${RegName} != "NONE" ] ; then + RegString="_${RegName}" + else + RegString="" + fi + +# if [ ! -z ${LowResMesh} ] && [ ${LowResMesh} != ${G_DEFAULT_LOW_RES_MESH} ]; then +# RegString="${RegString}.${LowResMesh}k" +# fi # This was tentatively commented out for NHP - Takuya Hayashi, perhaps G_DEFAULT_LOW_RES_MESH should be defined somewhere else for each species + + log_Msg "RegString: ${RegString}" + + # For interpreted modes, make sure that matlab/octave has access to the functions it needs. + # Since we are NOT using the ${FSL_FIXDIR}/call_matlab.sh script to invoke matlab (unlike 'hcp_fix') + # we need to explicitly add ${FSL_FIXDIR} (all the fix-related functions) + # and ${FSL_MATLAB_PATH} (e.g., read_avw.m, save_avw.m) to the matlab path. + # Several additional necessary environment variables (e.g., ${FSL_FIX_CIFTIRW} and ${FSL_FIX_WBC}) + # are set in FSL_FIXDIR/settings.sh, which is sourced below for interpreted modes. + # Note that the ciftiopen.m, ciftisave.m functions are added to the path through the ${FSL_FIX_WBC} + # environment variable within fix_3_clean.m itself. + export FSL_MATLAB_PATH="${FSLDIR}/etc/matlab" + local ML_PATHS="addpath('${FSL_FIXDIR}'); addpath('${FSL_MATLAB_PATH}');" + + #export FSL_FIX_CIFTIRW="${HCPPIPEDIR}/ReApplyFix/scripts" + #export FSL_FIX_WBC="${Caret7_Command}" + + # Make appropriate files if they don't exist + + local aggressive=0 + #local domot=1 + local hp=${HighPass} + local DoVol=0 + local fixlist=".fix" + + #local fixlist + + if (( hp >= 0 )); then + hpStr="_hp${hp}" + else + hpStr="" + fi + + # fMRIName is expected to NOT include path info, or a nifti extension; make sure that is indeed the case + # (although if someone includes the hp string as part fMRIName itself, the script will still break) + fMRIName=$(basename $($FSLDIR/bin/remove_ext $fMRIName)) + + if have_hand_reclassification ${StudyFolder} ${Subject} ${fMRIName} ${HighPass} ; then + fixlist="HandNoise.txt" + #TSC: if regname (which applies to the surface) isn't NONE, assume the hand classification was previously already applied to the volume data + if [[ "${RegName}" == "NONE" ]] + then + DoVol=1 + fi + fi + # WARNING: fix_3_clean doesn't actually do anything different based on the value of DoVol (its 5th argument). + # Rather, if a 5th argument is present, fix_3_clean does NOT apply cleanup to the volume, *regardless* of whether + # that 5th argument is 0 or 1 (or even a non-sensical string such as 'foo'). + # It is for that reason that the code below needs to use separate calls to fix_3_clean, with and without DoVol + # as an argument, rather than simply passing in the value of DoVol as set within this script. + # Not sure if/when this non-intuitive behavior of fix_3_clean will change, but this is accurate as of fix1.067 + log_Msg "Use fixlist=$fixlist" + + #local fmri_orig="${fMRIName}" + #local fmri=${fMRIName} + + DIR=$(pwd) + cd ${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName}/${fMRIName}${hpStr}.ica + + # Note: fix_3_clean does NOT filter the volume (NIFTI) data -- it assumes + # that any desired filtering has already been done outside of fix. + # So here, we need to symlink to the hp-filtered volume data. + $FSLDIR/bin/imln ../${fMRIName}${hpStr} filtered_func_data + + # However, hp-filtering of the CIFTI (dtseries) occurs within fix_3_clean. + # So here, we just create a symlink with the file name expected by + # fix_3_clean ("Atlas.dtseries.nii") to the non-filtered data. + if [ -f ../${fMRIName}_Atlas${RegString}.dtseries.nii ] ; then + log_Msg "FOUND FILE: ../${fMRIName}_Atlas${RegString}.dtseries.nii" + log_Msg "Performing imln" + + /bin/rm -f Atlas.dtseries.nii + $FSLDIR/bin/imln ../${fMRIName}_Atlas${RegString}.dtseries.nii Atlas.dtseries.nii + + log_Msg "START: Showing linked files" + ls -l ../${fMRIName}_Atlas${RegString}.dtseries.nii + ls -l Atlas.dtseries.nii + log_Msg "END: Showing linked files" + else + log_Warn "FILE NOT FOUND: ../${fMRIName}_Atlas${RegString}.dtseries.nii" + fi + + # Get Movement_Regressors.txt into the format expected by functionmotionconfounds.m + mkdir -p mc + if [ -f ../Movement_Regressors.txt ] ; then + log_Msg "Creating mc/prefiltered_func_data_mcf.par file" + cat ../Movement_Regressors.txt | awk '{ print $4 " " $5 " " $6 " " $1 " " $2 " " $3}' > mc/prefiltered_func_data_mcf.par + else + log_Err_Abort "Movement_Regressors.txt not retrieved properly." + fi + + ## --------------------------------------------------------------------------- + ## Run fix_3_clean + ## --------------------------------------------------------------------------- + + # MPH: We need to invoke fix_3_clean directly, rather than through 'fix -a ' because + # the latter does not provide access to the "DoVol" option within the former. + # (Also, 'fix -a' is hard-coded to use '.fix' as the list of noise components, although that + # could be worked around). + + export FSL_FIX_WBC="${Caret7_Command}" + # WARNING: fix_3_clean uses the environment variable FSL_FIX_WBC, but most previous + # versions of FSL_FIXDIR/settings.sh (v1.067 and earlier) have a hard-coded value for + # FSL_FIX_WBC, and don't check whether it is already defined in the environment. + # Thus, when settings.sh file gets sourced, there is a possibility that the version of + # wb_command is no longer the same as that specified by ${Caret7_Command}. So, after + # sourcing settings.sh below, we explicitly set FSL_FIX_WBC back to value of ${Caret7_Command}. + # (This may only be relevant for interpreted matlab/octave modes). + + log_Msg "Running fix_3_clean" + + case ${MatlabRunMode} in + + # See important WARNING above regarding why ${DoVol} is NOT included as an argument when DoVol=1 !! + + 0) + # Use Compiled MATLAB + + local matlab_exe="${FSL_FIXDIR}/compiled/$(uname -s)/$(uname -m)/run_fix_3_clean.sh" + + # Do NOT enclose string variables inside an additional single quote because all + # variables are already passed into the compiled binary as strings + local matlab_function_arguments=("${fixlist}" "${aggressive}" "${MotionRegression}" "${hp}") + if (( ! DoVol )); then + matlab_function_arguments+=("${DoVol}") + fi + + # fix_3_clean is part of the FIX distribution. + # If ${FSL_FIX_MCR} is already defined in the environment, use that for the MCR location. + # If not, the appropriate MCR version for use with fix_3_clean should be set in $FSL_FIXDIR/settings.sh. + if [ -z "${FSL_FIX_MCR}" ]; then + set +e + source ${FSL_FIXDIR}/settings.sh + set -e + export FSL_FIX_WBC="${Caret7_Command}" + # If FSL_FIX_MCR is still not defined after sourcing settings.sh, we have a problem + if [ -z "${FSL_FIX_MCR}" ]; then + log_Err_Abort "To use MATLAB run mode: ${MatlabRunMode}, the FSL_FIX_MCR environment variable must be set" + fi + fi + log_Msg "FSL_FIX_MCR: ${FSL_FIX_MCR}" + + local matlab_cmd=("${matlab_exe}" "${FSL_FIX_MCR}" "${matlab_function_arguments[@]}") + + # redirect tokens must be parsed by bash before doing variable expansion, and thus can't be inside a variable + # MPH: Going to let Compiled MATLAB use the existing stdout and stderr, rather than creating a separate log file + #local matlab_logfile=".reapplyfix.${fMRIName}${RegString}.fix_3_clean.matlab.log" + #log_Msg "Run MATLAB command: ${matlab_cmd[*]} >> ${matlab_logfile} 2>&1" + #"${matlab_cmd[@]}" >> "${matlab_logfile}" 2>&1 + log_Msg "Run compiled MATLAB: ${matlab_cmd[*]}" + "${matlab_cmd[@]}" + ;; + + 1 | 2) + # Use interpreted MATLAB or Octave + if [[ ${MatlabRunMode} == "1" ]]; then + local interpreter=(matlab -nojvm -nodisplay -nosplash) + else + local interpreter=(octave-cli -q --no-window-system) + fi + + if (( DoVol )); then + local matlab_cmd="${ML_PATHS} fix_3_clean('${fixlist}',${aggressive},${MotionRegression},${hp});" + else + local matlab_cmd="${ML_PATHS} fix_3_clean('${fixlist}',${aggressive},${MotionRegression},${hp},${DoVol});" + fi + + log_Msg "Run interpreted MATLAB/Octave (${interpreter[@]}) with command..." + log_Msg "${matlab_cmd}" + + # Use bash redirection ("here-string") to pass multiple commands into matlab + # (Necessary to protect the semicolons that separate matlab commands, which would otherwise + # get interpreted as separating different bash shell commands) + (set +e; source "${FSL_FIXDIR}/settings.sh"; set -e; export FSL_FIX_WBC="${Caret7_Command}"; "${interpreter[@]}" <<<"${matlab_cmd}") + ;; + + *) + # Unsupported MATLAB run mode + log_Err_Abort "Unsupported MATLAB run mode value: ${MatlabRunMode}" + ;; + esac + + log_Msg "Done running fix_3_clean" + + ## --------------------------------------------------------------------------- + ## Rename some files (relative to the default names coded in fix_3_clean) + ## --------------------------------------------------------------------------- + + # Remove any existing old versions of the cleaned data (normally they should be overwritten + # in the renaming that follows, but this ensures that any old versions don't linger) + cd ${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName} + local fmri=${fMRIName} + if (( hp >= 0 )); then + local fmrihp=${fmri}_hp${hp} + else + local fmrihp=${fmri} + fi + + /bin/rm -f ${fmri}_Atlas${RegString}${hpStr}_clean.dtseries.nii + /bin/rm -f ${fmri}_Atlas${RegString}${hpStr}_clean_vn.dscalar.nii + + if (( DoVol )); then + $FSLDIR/bin/imrm ${fmrihp}_clean + $FSLDIR/bin/imrm ${fmrihp}_clean_vn + fi + + # Rename some of the outputs from fix_3_clean. + # Note that the variance normalization ("_vn") outputs require use of fix1.067 or later + # So check whether those files exist before moving/renaming them + if [ -f ${fmrihp}.ica/Atlas_clean.dtseries.nii ] ; then + /bin/mv ${fmrihp}.ica/Atlas_clean.dtseries.nii ${fmri}_Atlas${RegString}${hpStr}_clean.dtseries.nii + else + log_Err_Abort "Something went wrong; ${fmrihp}.ica/Atlas_clean.dtseries.nii wasn't created" + fi + if [ -f ${fmrihp}.ica/Atlas_clean_vn.dscalar.nii ]; then + /bin/mv ${fmrihp}.ica/Atlas_clean_vn.dscalar.nii ${fmri}_Atlas${RegString}${hpStr}_clean_vn.dscalar.nii + fi + + if (( DoVol )); then + $FSLDIR/bin/immv ${fmrihp}.ica/filtered_func_data_clean ${fmrihp}_clean + if [ "$?" -ne "0" ]; then + log_Err_Abort "Something went wrong; ${fmrihp}.ica/filtered_func_data_clean wasn't created" + fi + if [ `$FSLDIR/bin/imtest ${fmrihp}.ica/filtered_func_data_clean_vn` = 1 ]; then + $FSLDIR/bin/immv ${fmrihp}.ica/filtered_func_data_clean_vn ${fmrihp}_clean_vn + fi + fi + log_Msg "Done renaming files" + + # Remove the 'fake-NIFTI' file created in fix_3_clean for high-pass filtering of the CIFTI (if it exists) + $FSLDIR/bin/imrm ${fmrihp}.ica/Atlas + + cd ${DIR} + + log_Msg "Completed!" + +} + +# ------------------------------------------------------------------------------ +# "Global" processing - everything above here should be in a function +# ------------------------------------------------------------------------------ + +set -e # If any command exits with non-zero value, this script exits + +# Establish defaults +G_DEFAULT_REG_NAME="NONE" +G_DEFAULT_LOW_RES_MESH=32 +G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB +G_DEFAULT_MOTION_REGRESSION="FALSE" + +# Set global variables +g_script_name=$(basename "${0}") + +# Allow script to return a Usage statement, before any other output +if [ "$#" = "0" ]; then + usage + exit 1 +fi + +# Verify that HCPPIPEDIR environment variable is set +if [ -z "${HCPPIPEDIR}" ]; then + echo "$(basename ${0}): ABORTING: HCPPIPEDIR environment variable must be set" + exit 1 +fi + +# Load function libraries +source "${HCPPIPEDIR}/global/scripts/log.shlib" # Logging related functions +source "${HCPPIPEDIR}/global/scripts/fsl_version.shlib" # Function for getting FSL version +log_Msg "HCPPIPEDIR: ${HCPPIPEDIR}" + +# Verify any other needed environment variables are set +log_Check_Env_Var CARET7DIR +log_Check_Env_Var FSLDIR +log_Check_Env_Var FSL_FIXDIR + +# Show tool versions +show_tool_versions + +# Determine whether named or positional parameters are used +if [[ ${1} == --* ]]; then + # Named parameters (e.g. --parameter-name=parameter-value) are used + log_Msg "Using named parameters" + + # Get command line options + get_options "$@" + + # Invoke main functionality + # ${1} ${2} ${3} ${4} ${5} ${6} ${7} ${8} + main "${p_StudyFolder}" "${p_Subject}" "${p_fMRIName}" "${p_HighPass}" "${p_RegName}" "${p_LowResMesh}" "${p_MatlabRunMode}" "${p_MotionRegression}" + +else + # Positional parameters are used + log_Msg "Using positional parameters" + main "$@" + +fi From f6f3de362acc89a1b32bd86f2d840d0b55bb82f1 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 6 Jul 2022 16:20:06 +0900 Subject: [PATCH 27/29] modify same to HCP v4.30 --- ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh | 184 +++++++++++++++-------- 1 file changed, 125 insertions(+), 59 deletions(-) diff --git a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh index e971f27..c9db521 100755 --- a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh +++ b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh @@ -1,5 +1,4 @@ #!/bin/bash - # # # ReApplyFixPipeline.sh # @@ -32,14 +31,11 @@ # Show usage information for this script # ------------------------------------------------------------------------------ -usage() +show_usage() { - local script_name - script_name=$(basename "${0}") - cat < = user supplied value @@ -56,13 +52,13 @@ PARAMETERs are [ ] = optional; < > = user supplied value form) by simply specifying all values on the command line in the order they are listed below. - e.g. ${script_name} /path/to/study/folder 100307 rfMRI_REST1_LR 2000 ... + e.g. ${g_script_name} /path/to/study/folder 100307 rfMRI_REST1_LR 2000 ... [--help] : show usage information and exit --path= OR --study-folder= --subject= --fmri-name= (Do not include path, nifti extension, or the 'hp' string). - --high-pass= Number to represent the ${HighPass} variable used in ICA+FIX + --high-pass= [--reg-name= defaults to ${G_DEFAULT_REG_NAME}. (Use NONE for MSMSulc registration) [--low-res-mesh=] defaults to ${G_DEFAULT_LOW_RES_MESH} [--matlab-run-mode={0, 1, 2}] defaults to ${G_DEFAULT_MATLAB_RUN_MODE} @@ -70,18 +66,26 @@ PARAMETERs are [ ] = optional; < > = user supplied value 1 = Use interpreted MATLAB 2 = Use interpreted Octave [--motion-regression={TRUE, FALSE}] defaults to ${G_DEFAULT_MOTION_REGRESSION} + [--delete-intermediates={TRUE, FALSE}] defaults to ${G_DEFAULT_DELETE_INTERMEDIATES} + If TRUE, deletes the high-pass filtered files. EOF -exit 1; } -[ "$1" = "" ] && usage + +# Establish defaults +G_DEFAULT_REG_NAME="NONE" +G_DEFAULT_LOW_RES_MESH=32 +G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB +G_DEFAULT_MOTION_REGRESSION="FALSE" +G_DEFAULT_DELETE_INTERMEDIATES="FALSE" + # ------------------------------------------------------------------------------ # Get the command line options for this script. # ------------------------------------------------------------------------------ get_options() { - local arguments=($@) + local arguments=("$@") # initialize global output variables unset p_StudyFolder # ${1} @@ -92,12 +96,14 @@ get_options() unset p_LowResMesh # ${6} unset p_MatlabRunMode # ${7} unset p_MotionRegression # ${8} + unset p_DeleteIntermediates # ${9} # set default values p_RegName=${G_DEFAULT_REG_NAME} p_LowResMesh=${G_DEFAULT_LOW_RES_MESH} p_MatlabRunMode=${G_DEFAULT_MATLAB_RUN_MODE} p_MotionRegression=${G_DEFAULT_MOTION_REGRESSION} + p_DeleteIntermediates=${G_DEFAULT_DELETE_INTERMEDIATES} # parse arguments local num_args=${#arguments[@]} @@ -109,8 +115,8 @@ get_options() case ${argument} in --help) - usage - exit 1 + show_usage + exit 0 ;; --path=*) p_StudyFolder=${argument#*=} @@ -148,8 +154,12 @@ get_options() p_MotionRegression=${argument#*=} index=$(( index + 1 )) ;; + --delete-intermediates=*) + p_DeleteIntermediates=${argument#*=} + index=$(( index + 1 )) + ;; *) - usage + show_usage log_Err_Abort "unrecognized option: ${argument}" ;; esac @@ -228,6 +238,13 @@ get_options() log_Msg "Motion Regression: ${p_MotionRegression}" fi + if [ -z "${p_DeleteIntermediates}" ]; then + log_Err "delete intermediates setting (--delete-intermediates=) required" + error_count=$(( error_count + 1 )) + else + log_Msg "Delete Intermediates: ${p_DeleteIntermediates}" + fi + if [ ${error_count} -gt 0 ]; then log_Err_Abort "For usage information, use --help" fi @@ -251,6 +268,13 @@ show_tool_versions() log_Msg "Showing FSL version" fsl_version_get fsl_ver log_Msg "FSL version: ${fsl_ver}" + + # Show specific FIX version, if available + if [ -f ${FSL_FIXDIR}/fixversion ]; then + fixversion=$(cat ${FSL_FIXDIR}/fixversion ) + log_Msg "FIX version: $fixversion" + fi + } # ------------------------------------------------------------------------------ @@ -271,6 +295,21 @@ have_hand_reclassification() fi } +function interpret_as_bool() +{ + case $(echo "$1" | tr '[:upper:]' '[:lower:]') in + (true | yes | 1) + echo 1 + ;; + (false | no | none | 0) + echo 0 + ;; + (*) + log_Err_Abort "error: '$1' is not valid for this argument, please use TRUE or FALSE" + ;; + esac +} + # ------------------------------------------------------------------------------ # Main processing of script. # ------------------------------------------------------------------------------ @@ -284,8 +323,9 @@ main() local Subject="${2}" local fMRIName="${3}" local HighPass="${4}" - local RegName="${5}" - if [ -z "${5}" ]; then + + local RegName + if [ -z "${5}" ]; then RegName=${G_DEFAULT_REG_NAME} else RegName="${5}" @@ -307,9 +347,9 @@ main() local MotionRegression if [ -z "${8}" ]; then - MotionRegression="${G_DEFAULT_MOTION_REGRESSION}" + MotionRegression=$(interpret_as_bool "${G_DEFAULT_MOTION_REGRESSION}") else - MotionRegression="${8}" + MotionRegression=$(interpret_as_bool "${8}") fi # Turn MotionRegression into an appropriate numeric value for fix_3_clean @@ -325,6 +365,13 @@ main() ;; esac + local DeleteIntermediates + if [ -z "${9}" ]; then + DeleteIntermediates=$(interpret_as_bool "${G_DEFAULT_DELETE_INTERMEDIATES}") + else + DeleteIntermediates=$(interpret_as_bool "${9}") + fi + # Log values retrieved from positional parameters log_Msg "StudyFolder: ${StudyFolder}" log_Msg "Subject: ${Subject}" @@ -334,6 +381,7 @@ main() log_Msg "LowResMesh: ${LowResMesh}" log_Msg "MatlabRunMode: ${MatlabRunMode}" log_Msg "MotionRegression: ${MotionRegression}" + log_Msg "DeleteIntermediates: ${DeleteIntermediates}" # Naming Conventions and other variables local Caret7_Command="${CARET7DIR}/wb_command" @@ -352,30 +400,26 @@ main() log_Msg "RegString: ${RegString}" - # For interpreted modes, make sure that matlab/octave has access to the functions it needs. + # For INTERPRETED MODES, make sure that matlab/octave has access to the functions it needs. # Since we are NOT using the ${FSL_FIXDIR}/call_matlab.sh script to invoke matlab (unlike 'hcp_fix') # we need to explicitly add ${FSL_FIXDIR} (all the fix-related functions) # and ${FSL_MATLAB_PATH} (e.g., read_avw.m, save_avw.m) to the matlab path. # Several additional necessary environment variables (e.g., ${FSL_FIX_CIFTIRW} and ${FSL_FIX_WBC}) - # are set in FSL_FIXDIR/settings.sh, which is sourced below for interpreted modes. - # Note that the ciftiopen.m, ciftisave.m functions are added to the path through the ${FSL_FIX_WBC} + # are set in ${FSL_FIXDIR}/settings.sh, which is sourced below for interpreted modes. + # Note that the ciftiopen.m, ciftisave.m functions are *appended* to the path through the ${FSL_FIX_CIFTIRW} # environment variable within fix_3_clean.m itself. + # Note that we are NOT adding '${HCPPIPEDIR}/global/matlab' to the path (which also contains the + # CIFTI I/O functions as well), so the versions in ${FSL_FIX_CIFTIRW} will be the ones that get used + # (assuming that the user hasn't further customized their matlab path, e.g., via a startup.m addition). export FSL_MATLAB_PATH="${FSLDIR}/etc/matlab" local ML_PATHS="addpath('${FSL_FIXDIR}'); addpath('${FSL_MATLAB_PATH}');" - #export FSL_FIX_CIFTIRW="${HCPPIPEDIR}/ReApplyFix/scripts" - #export FSL_FIX_WBC="${Caret7_Command}" - - # Make appropriate files if they don't exist - + # Some defaults local aggressive=0 - #local domot=1 local hp=${HighPass} local DoVol=0 local fixlist=".fix" - #local fixlist - if (( hp >= 0 )); then hpStr="_hp${hp}" else @@ -386,7 +430,9 @@ main() # (although if someone includes the hp string as part fMRIName itself, the script will still break) fMRIName=$(basename $($FSLDIR/bin/remove_ext $fMRIName)) - if have_hand_reclassification ${StudyFolder} ${Subject} ${fMRIName} ${HighPass} ; then + # If we have a hand classification and no regname, reapply fix to the volume as well + if have_hand_reclassification ${StudyFolder} ${Subject} ${fMRIName} ${hp} + then fixlist="HandNoise.txt" #TSC: if regname (which applies to the surface) isn't NONE, assume the hand classification was previously already applied to the volume data if [[ "${RegName}" == "NONE" ]] @@ -400,11 +446,12 @@ main() # It is for that reason that the code below needs to use separate calls to fix_3_clean, with and without DoVol # as an argument, rather than simply passing in the value of DoVol as set within this script. # Not sure if/when this non-intuitive behavior of fix_3_clean will change, but this is accurate as of fix1.067 + # UPDATE (11/8/2019): As of FIX 1.06.12, fix_3_clean interprets its 5th argument ("DoVol") in the usual boolean + # manner. However, since we already had a work-around to this problem, we will leave the code unchanged so that + # we don't need to add a FIX version dependency to the script. + log_Msg "Use fixlist=$fixlist" - #local fmri_orig="${fMRIName}" - #local fmri=${fMRIName} - DIR=$(pwd) cd ${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName}/${fMRIName}${hpStr}.ica @@ -413,7 +460,7 @@ main() # So here, we need to symlink to the hp-filtered volume data. $FSLDIR/bin/imln ../${fMRIName}${hpStr} filtered_func_data - # However, hp-filtering of the CIFTI (dtseries) occurs within fix_3_clean. + # However, hp-filtering of the *CIFTI* (dtseries) occurs within fix_3_clean. # So here, we just create a symlink with the file name expected by # fix_3_clean ("Atlas.dtseries.nii") to the non-filtered data. if [ -f ../${fMRIName}_Atlas${RegString}.dtseries.nii ] ; then @@ -480,9 +527,9 @@ main() # If ${FSL_FIX_MCR} is already defined in the environment, use that for the MCR location. # If not, the appropriate MCR version for use with fix_3_clean should be set in $FSL_FIXDIR/settings.sh. if [ -z "${FSL_FIX_MCR}" ]; then - set +e + debug_disable_trap source ${FSL_FIXDIR}/settings.sh - set -e + debug_enable_trap export FSL_FIX_WBC="${Caret7_Command}" # If FSL_FIX_MCR is still not defined after sourcing settings.sh, we have a problem if [ -z "${FSL_FIX_MCR}" ]; then @@ -522,7 +569,7 @@ main() # Use bash redirection ("here-string") to pass multiple commands into matlab # (Necessary to protect the semicolons that separate matlab commands, which would otherwise # get interpreted as separating different bash shell commands) - (set +e; source "${FSL_FIXDIR}/settings.sh"; set -e; export FSL_FIX_WBC="${Caret7_Command}"; "${interpreter[@]}" <<<"${matlab_cmd}") + (debug_disable_trap; source "${FSL_FIXDIR}/settings.sh"; debug_enable_trap; export FSL_FIX_WBC="${Caret7_Command}"; "${interpreter[@]}" <<<"${matlab_cmd}") ;; *) @@ -558,7 +605,7 @@ main() # Rename some of the outputs from fix_3_clean. # Note that the variance normalization ("_vn") outputs require use of fix1.067 or later # So check whether those files exist before moving/renaming them - if [ -f ${fmrihp}.ica/Atlas_clean.dtseries.nii ] ; then + if [ -f ${fmrihp}.ica/Atlas_clean.dtseries.nii ]; then /bin/mv ${fmrihp}.ica/Atlas_clean.dtseries.nii ${fmri}_Atlas${RegString}${hpStr}_clean.dtseries.nii else log_Err_Abort "Something went wrong; ${fmrihp}.ica/Atlas_clean.dtseries.nii wasn't created" @@ -581,45 +628,64 @@ main() # Remove the 'fake-NIFTI' file created in fix_3_clean for high-pass filtering of the CIFTI (if it exists) $FSLDIR/bin/imrm ${fmrihp}.ica/Atlas - cd ${DIR} + # Always delete things with too-generic names + $FSLDIR/bin/imrm ${fmrihp}.ica/filtered_func_data + rm -f ${fmrihp}.ica/Atlas.dtseries.nii + + # Optional deletion of highpass intermediates + if [ "${DeleteIntermediates}" == "1" ] ; then + if (( hp > 0 )); then # fix_3_clean only writes out the hp-filtered time series if hp > 0 + $FSLDIR/bin/imrm ${fmri}_hp${hp} # Explicitly use _hp${hp} here (rather than $hpStr as a safeguard against accidental deletion of the non-hp-filtered timeseries) + rm -f ${fmrihp}.ica/Atlas_hp_preclean.dtseries.nii + fi + else + #even if we don't delete it, don't leave this file with a hard to interpret name + if (( hp > 0 )); then + # 'OR' mv command with "true" to avoid returning an error code if file doesn't exist for some reason + mv -f ${fmrihp}.ica/Atlas_hp_preclean.dtseries.nii ${fmri}_Atlas_hp${hp}.dtseries.nii || true + fi + fi - log_Msg "Completed!" + cd ${DIR} # Return to directory where script was launched + log_Msg "Completed!" } # ------------------------------------------------------------------------------ # "Global" processing - everything above here should be in a function # ------------------------------------------------------------------------------ -set -e # If any command exits with non-zero value, this script exits - -# Establish defaults -G_DEFAULT_REG_NAME="NONE" -G_DEFAULT_LOW_RES_MESH=32 -G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB -G_DEFAULT_MOTION_REGRESSION="FALSE" - # Set global variables g_script_name=$(basename "${0}") # Allow script to return a Usage statement, before any other output if [ "$#" = "0" ]; then - usage + show_usage exit 1 fi # Verify that HCPPIPEDIR environment variable is set if [ -z "${HCPPIPEDIR}" ]; then - echo "$(basename ${0}): ABORTING: HCPPIPEDIR environment variable must be set" + echo "${g_script_name}: ABORTING: HCPPIPEDIR environment variable must be set" exit 1 fi # Load function libraries -source "${HCPPIPEDIR}/global/scripts/log.shlib" # Logging related functions -source "${HCPPIPEDIR}/global/scripts/fsl_version.shlib" # Function for getting FSL version -log_Msg "HCPPIPEDIR: ${HCPPIPEDIR}" +source "${HCPPIPEDIR}/global/scripts/debug.shlib" "$@" # Debugging functions; also sources log.shlib +source ${HCPPIPEDIR}/global/scripts/opts.shlib # Command line option functions +source "${HCPPIPEDIR}/global/scripts/fsl_version.shlib" # Functions for getting FSL version -# Verify any other needed environment variables are set +opts_ShowVersionIfRequested $@ + +if opts_CheckForHelpRequest $@; then + show_usage + exit 0 +fi + +${HCPPIPEDIR}/show_version + +# Verify required environment variables are set and log value +log_Check_Env_Var HCPPIPEDIR log_Check_Env_Var CARET7DIR log_Check_Env_Var FSLDIR log_Check_Env_Var FSL_FIXDIR @@ -627,7 +693,7 @@ log_Check_Env_Var FSL_FIXDIR # Show tool versions show_tool_versions -# Determine whether named or positional parameters are used +# Determine whether named or positional parameters are used and invoke the 'main' function if [[ ${1} == --* ]]; then # Named parameters (e.g. --parameter-name=parameter-value) are used log_Msg "Using named parameters" @@ -636,9 +702,9 @@ if [[ ${1} == --* ]]; then get_options "$@" # Invoke main functionality - # ${1} ${2} ${3} ${4} ${5} ${6} ${7} ${8} - main "${p_StudyFolder}" "${p_Subject}" "${p_fMRIName}" "${p_HighPass}" "${p_RegName}" "${p_LowResMesh}" "${p_MatlabRunMode}" "${p_MotionRegression}" - + # ${1} ${2} ${3} ${4} ${5} ${6} ${7} ${8} ${9} + main "${p_StudyFolder}" "${p_Subject}" "${p_fMRIName}" "${p_HighPass}" "${p_RegName}" "${p_LowResMesh}" "${p_MatlabRunMode}" "${p_MotionRegression}" "${p_DeleteIntermediates}" + else # Positional parameters are used log_Msg "Using positional parameters" From 0280ceb2f95776e54309d1df4fc11f046792ff00 Mon Sep 17 00:00:00 2001 From: TakJim Date: Wed, 6 Jul 2022 16:39:34 +0900 Subject: [PATCH 28/29] modify same to HCP v4.30 --- ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh | 61 ++++++++++++++++++------ 1 file changed, 46 insertions(+), 15 deletions(-) diff --git a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh index c9db521..280e993 100755 --- a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh +++ b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh @@ -352,19 +352,6 @@ main() MotionRegression=$(interpret_as_bool "${8}") fi - # Turn MotionRegression into an appropriate numeric value for fix_3_clean - case $(echo ${MotionRegression} | tr '[:upper:]' '[:lower:]') in - ( true | yes | 1) - MotionRegression=1 - ;; - ( false | no | none | 0) - MotionRegression=0 - ;; - *) - log_Err_Abort "motion regression setting must be TRUE or FALSE" - ;; - esac - local DeleteIntermediates if [ -z "${9}" ]; then DeleteIntermediates=$(interpret_as_bool "${G_DEFAULT_DELETE_INTERMEDIATES}") @@ -453,12 +440,56 @@ main() log_Msg "Use fixlist=$fixlist" DIR=$(pwd) - cd ${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName}/${fMRIName}${hpStr}.ica + cd ${StudyFolder}/${Subject}/MNINonLinear/Results/${fMRIName} # Note: fix_3_clean does NOT filter the volume (NIFTI) data -- it assumes # that any desired filtering has already been done outside of fix. # So here, we need to symlink to the hp-filtered volume data. - $FSLDIR/bin/imln ../${fMRIName}${hpStr} filtered_func_data + # HOWEVER, if missing, only need to generate the hp-filtered volume data if DoVol=1. + # Otherwise (if DoVol=0), the only role of filtered_func_data in fix_3_clean is to determine the TR. + # In that case, we will just symlink the *non-filtered* data to filtered_func_data + # (as a hack for fix_3_clean to determine the TR without a time-consuming filtering step + # on the volume) + + useNonFilteredAsFilteredFunc=0 + if (( $($FSLDIR/bin/imtest "${fMRIName}${hpStr}") )); then + log_Warn "Using existing $($FSLDIR/bin/imglob -extension ${fMRIName}${hpStr}) (not re-filtering)" + else # hp filtered volume file doesn't exist + if (( DoVol )); then # need to actually refilter the volume data + if (( hp > 0 )); then + tr=`$FSLDIR/bin/fslval ${fMRIName} pixdim4` + log_Msg "tr: ${tr}" + log_Msg "processing FMRI file ${fMRIName} with highpass ${hp}" + hptr=$(echo "scale = 10; $hp / (2 * $tr)" | bc -l) + + # Starting with FSL 5.0.7, 'fslmaths -bptf' no longer includes the temporal mean in its output. + # A work-around to this, which works with both the pre- and post-5.0.7 behavior is to compute + # the temporal mean, remove it, run -bptf, and then add the mean back in. + ${FSLDIR}/bin/fslmaths ${fMRIName} -Tmean ${fMRIName}${hpStr} + highpass_cmd="${FSLDIR}/bin/fslmaths ${fMRIName} -sub ${fMRIName}${hpStr} -bptf ${hptr} -1 -add ${fMRIName}${hpStr} ${fMRIName}${hpStr}" + log_Msg "highpass_cmd: ${highpass_cmd}" + ${highpass_cmd} + elif (( hp == 0 )); then + # Nothing in script currently detrends the volume if hp=0 is requested (which is the intended meaning of hp=0) + log_Err_Abort "hp = ${hp} not currently supported" + fi + else + useNonFilteredAsFilteredFunc=1 + fi + fi + + if [ ! -e ${fMRIName}${hpStr}.ica ]; then + log_Err_Abort "${fMRIName}${hpStr}.ica is expected to already exist, but does not" + fi + + cd ${fMRIName}${hpStr}.ica + + # Create symlink for filtered_func_data (per comments above) + if (( useNonFilteredAsFilteredFunc )); then + $FSLDIR/bin/imln ../${fMRIName} filtered_func_data + else + $FSLDIR/bin/imln ../${fMRIName}${hpStr} filtered_func_data + fi # However, hp-filtering of the *CIFTI* (dtseries) occurs within fix_3_clean. # So here, we just create a symlink with the file name expected by From 5737dfaf9c4e34d56e88456f7d70c63f480245b7 Mon Sep 17 00:00:00 2001 From: TakJim Date: Thu, 7 Jul 2022 09:21:52 +0900 Subject: [PATCH 29/29] Set G_DEFAULT_LOW_RES_MESH for each species --- ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh index 280e993..05264ee 100755 --- a/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh +++ b/ICAFIX/ReApplyFixPipeline_RIKEN_merge.sh @@ -74,7 +74,13 @@ EOF # Establish defaults G_DEFAULT_REG_NAME="NONE" -G_DEFAULT_LOW_RES_MESH=32 +if [[ $SPECIES =~ "Macaque ]] ; then + G_DEFAULT_LOW_RES_MESH=10 +elif [[ $SPECIES =~ "Marmoset" ]] ; then + G_DEFAULT_LOW_RES_MESH=4 +else + G_DEFAULT_LOW_RES_MESH=32 +fi G_DEFAULT_MATLAB_RUN_MODE=1 # Use interpreted MATLAB G_DEFAULT_MOTION_REGRESSION="FALSE" G_DEFAULT_DELETE_INTERMEDIATES="FALSE" @@ -381,9 +387,9 @@ main() RegString="" fi -# if [ ! -z ${LowResMesh} ] && [ ${LowResMesh} != ${G_DEFAULT_LOW_RES_MESH} ]; then -# RegString="${RegString}.${LowResMesh}k" -# fi # This was tentatively commented out for NHP - Takuya Hayashi, perhaps G_DEFAULT_LOW_RES_MESH should be defined somewhere else for each species + if [ ! -z ${LowResMesh} ] && [ ${LowResMesh} != ${G_DEFAULT_LOW_RES_MESH} ]; then + RegString="${RegString}.${LowResMesh}k" + fi log_Msg "RegString: ${RegString}"