search

cpp-lln-lab / bidspm

an SPM centric BIDS app
https://bidspm.readthedocs.io/en/latest/general_information.html
GNU General Public License v3.0
6 stars 14 forks source link

[BUG] Anatomical preprocessing only happens for the first session in querry #1184

Open JacMatu opened 9 months ago

JacMatu commented 9 months ago

Is there an existing issue for this?

Operating system

Operating system version

SPM 12 version

Platform

Platform version

bidspm version

v3.0.0

bidspm branch / commit number

No response

Expected Behavior

Using preprocessing with 'anat_only' = true with data from multiple sessions (varrying across subjects) to end up in preprocessing done for each session instead of just the first one from query. 

parfor i = 1:numel(subject_label)
       bidspm(bids_dir, output_dir, 'subject', ...
              'action', 'preprocess', ...
              'participant_label', subject_label(i), ...
              'space', {'individual', 'IXI549Space'}, ...
              'anat_only', true, ...
              'options', opt, ...
              'dry_run', false, ...
              'verbosity', 2, ...
              'skip_validation', true);
end

When I query for specific session preproc happens for it (e.g. ses-02 or 03, but when I grab 02, 03 and 04 it only happens for 02.)

Anything else?

matlabbatch (job) is built for one session only.

matlabbatch{1}.cfg_basicio.cfg_named_file.name = 'Anatomical';
matlabbatch{1}.cfg_basicio.cfg_named_file.files = {{'/Volumes/Slim_Reaper/Projects/High_res_7T_Liege_myelin/analysis_high-res_BLAM_anat-freesurfer-laynii/outputs/test_preproc_then_average/derivatives/bidspm-preproc/sub-SC01/ses-01/anat/sub-SC01_ses-01_acq-r0p75_UNIT1.nii'}};
matlabbatch{2}.spm.spatial.preproc.channel.write = [false true];
matlabbatch{2}.spm.spatial.preproc.channel.vols(1) = cfg_dep('Named File Selector: Anatomical(1) - Files', substruct('.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','files', '{}',{1}));
matlabbatch{2}.spm.spatial.preproc.channel.biasreg = 0.001;
matlabbatch{2}.spm.spatial.preproc.channel.biasfwhm = 60;
matlabbatch{2}.spm.spatial.preproc.tissue(1).tpm = {'/Users/jacekmatuszewski/Documents/spm12/tpm/TPM.nii,1'};
matlabbatch{2}.spm.spatial.preproc.tissue(1).ngaus = 1;
matlabbatch{2}.spm.spatial.preproc.tissue(1).native = [1 0];
matlabbatch{2}.spm.spatial.preproc.tissue(1).warped = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(2).tpm = {'/Users/jacekmatuszewski/Documents/spm12/tpm/TPM.nii,2'};
matlabbatch{2}.spm.spatial.preproc.tissue(2).ngaus = 1;
matlabbatch{2}.spm.spatial.preproc.tissue(2).native = [1 0];
matlabbatch{2}.spm.spatial.preproc.tissue(2).warped = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(3).tpm = {'/Users/jacekmatuszewski/Documents/spm12/tpm/TPM.nii,3'};
matlabbatch{2}.spm.spatial.preproc.tissue(3).ngaus = 2;
matlabbatch{2}.spm.spatial.preproc.tissue(3).native = [1 0];
matlabbatch{2}.spm.spatial.preproc.tissue(3).warped = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(4).tpm = {'/Users/jacekmatuszewski/Documents/spm12/tpm/TPM.nii,4'};
matlabbatch{2}.spm.spatial.preproc.tissue(4).ngaus = 3;
matlabbatch{2}.spm.spatial.preproc.tissue(4).native = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(4).warped = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(5).tpm = {'/Users/jacekmatuszewski/Documents/spm12/tpm/TPM.nii,5'};
matlabbatch{2}.spm.spatial.preproc.tissue(5).ngaus = 4;
matlabbatch{2}.spm.spatial.preproc.tissue(5).native = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(5).warped = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(6).tpm = {'/Users/jacekmatuszewski/Documents/spm12/tpm/TPM.nii,6'};
matlabbatch{2}.spm.spatial.preproc.tissue(6).ngaus = 2;
matlabbatch{2}.spm.spatial.preproc.tissue(6).native = [0 0];
matlabbatch{2}.spm.spatial.preproc.tissue(6).warped = [0 0];
matlabbatch{2}.spm.spatial.preproc.warp.mrf = 1;
matlabbatch{2}.spm.spatial.preproc.warp.cleanup = 1;
matlabbatch{2}.spm.spatial.preproc.warp.reg = [0 0.001 0.5 0.05 0.2];
matlabbatch{2}.spm.spatial.preproc.warp.affreg = 'mni';
matlabbatch{2}.spm.spatial.preproc.warp.fwhm = 0;
matlabbatch{2}.spm.spatial.preproc.warp.samp = 3;
matlabbatch{2}.spm.spatial.preproc.warp.write = [1 1];
matlabbatch{3}.spm.util.imcalc.input(1) = cfg_dep('Segment: Bias Corrected (1)', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','channel', '()',{1}, '.','biascorr', '()',{':'}));
matlabbatch{3}.spm.util.imcalc.input(2) = cfg_dep('Segment: c1 Images', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{1}, '.','c', '()',{':'}));
matlabbatch{3}.spm.util.imcalc.input(3) = cfg_dep('Segment: c2 Images', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{2}, '.','c', '()',{':'}));
matlabbatch{3}.spm.util.imcalc.input(4) = cfg_dep('Segment: c3 Images', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{3}, '.','c', '()',{':'}));
matlabbatch{3}.spm.util.imcalc.output = 'sub-SC01_ses-01_acq-r0p75_space-individual_desc-skullstripped_UNIT1.nii';
matlabbatch{3}.spm.util.imcalc.outdir = {'/Volumes/Slim_Reaper/Projects/High_res_7T_Liege_myelin/analysis_high-res_BLAM_anat-freesurfer-laynii/outputs/test_preproc_then_average/derivatives/bidspm-preproc/sub-SC01/ses-01/anat'};
matlabbatch{3}.spm.util.imcalc.expression = 'i1.*((i2+i3+i4)>0.750000)';
matlabbatch{3}.spm.util.imcalc.options.dtype = 16;
matlabbatch{4}.spm.util.imcalc.input(1) = cfg_dep('Segment: Bias Corrected (1)', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','channel', '()',{1}, '.','biascorr', '()',{':'}));
matlabbatch{4}.spm.util.imcalc.input(2) = cfg_dep('Segment: c1 Images', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{1}, '.','c', '()',{':'}));
matlabbatch{4}.spm.util.imcalc.input(3) = cfg_dep('Segment: c2 Images', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{2}, '.','c', '()',{':'}));
matlabbatch{4}.spm.util.imcalc.input(4) = cfg_dep('Segment: c3 Images', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{3}, '.','c', '()',{':'}));
matlabbatch{4}.spm.util.imcalc.output = 'sub-SC01_ses-01_acq-r0p75_space-individual_desc-brain_mask.nii';
matlabbatch{4}.spm.util.imcalc.outdir = {'/Volumes/Slim_Reaper/Projects/High_res_7T_Liege_myelin/analysis_high-res_BLAM_anat-freesurfer-laynii/outputs/test_preproc_then_average/derivatives/bidspm-preproc/sub-SC01/ses-01/anat'};
matlabbatch{4}.spm.util.imcalc.expression = '(i2+i3+i4)>0.750000';
matlabbatch{4}.spm.util.imcalc.options.dtype = 16;
matlabbatch{5}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{5}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Segment: bias corrected image', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','channel', '()',{1}, '.','biascorr', '()',{':'}));
matlabbatch{5}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
matlabbatch{6}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{6}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Segment: grey matter image', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{1}, '.','c', '()',{':'}));
matlabbatch{7}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{7}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Segment: white matter image', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{2}, '.','c', '()',{':'}));
matlabbatch{8}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{8}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Segment: csf matter image', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','tiss', '()',{3}, '.','c', '()',{':'}));
matlabbatch{9}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{9}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Image Calculator: skullstripped anatomical', substruct('.','val', '{}',{3}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','files'));
matlabbatch{9}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
matlabbatch{10}.spm.spatial.normalise.write.subj.def(1) = cfg_dep('Segment: Forward Deformations', substruct('.','val', '{}',{2}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','fordef', '()',{':'}));
matlabbatch{10}.spm.spatial.normalise.write.subj.resample(1) = cfg_dep('Image Calculator: skullstrip mask', substruct('.','val', '{}',{4}, '.','val', '{}',{1}, '.','val', '{}',{1}), substruct('.','files'));
matlabbatch{10}.spm.spatial.normalise.write.woptions.vox = [1 1 1];
matlabbatch{10}.spm.spatial.normalise.write.woptions.interp = 0;
github-actions[bot] commented 9 months ago

Thank you for your issue. Give us a little time to review it.

PS. You might want to check the FAQ if you haven't done so already.

This is an automated reply, generated by FAQtory

Remi-Gau commented 9 months ago

OK I can indeed reproduce the bug.

The question is what is the expected anat only preprocessing pipeline for several sessions of anat.

average anat across sessions then preprocess

Assuming one T1w per session (to keep things simple), we coregister all T1w to that of the first session, then we average them, then we preprocess that single image.

preprocess sessions independently

Assuming one T1w per session (to keep things simple), we coregister all T1w to that of the first session, then we preprocess that each image from each session.

The second one is easier to do for now. The question is which one do you need?

JacMatu commented 9 months ago

I defo need the second one, Marco pointed out that the bias corrections should be done independently for each session since the noise distribution might be session-specific, so even in case of averaging an image across multiple sessions it might be better to first preprocess them individually and then average the corrected products.

If you will include the realignment I will change my scripts to just average them. One thing that I didn't test is I switched the degree of interpolation to max (7th degree) for more precision (& more computational time) - which might be better for high-res studies.

I might also be nice to include in the anatomical preproc only batch, since there are typically fewer images to realign (in opposition to bazillion frames of func acquisitions).

In realign and unwarp batch it was the field (instead of default 3) matlabbatch{1}.spm.spatial.realignunwarp.uwroptions.rinterp = 7;

What do you think?

Remi-Gau commented 9 months ago

One thing that I didn't test is I switched the degree of interpolation to max (7th degree) for more precision (& more computational time) - which might be better for high-res studies.

In realign and unwarp batch it was the field (instead of default 3) matlabbatch{1}.spm.spatial.realignunwarp.uwroptions.rinterp = 7;

OK for this I would let users modify it by adding an spm_my_default.m file in the matlab path

Here is the one used by BIDSPM: https://github.com/cpp-lln-lab/bidspm/blob/main/src/defaults/spm_my_defaults.m

This allows you to edit the default values by SPM in the batches without having to hard code it in your script.

Remi-Gau commented 9 months ago

I might also be nice to include in the anatomical preproc only batch, since there are typically fewer images to realign (in opposition to bazillion frames of func acquisitions).

I may end up doing this if I can refactor it out of the current preprocessing.

Remi-Gau commented 9 months ago

I defo need the second one, Marco pointed out that the bias corrections should be done independently for each session since the noise distribution might be session-specific, so even in case of averaging an image across multiple sessions it might be better to first preprocess them individually and then average the corrected products.

Fair enough, will start with this then.