Hello everyone, I am currently experimenting with a denoising strategy that utilizes the first 5 aCompCor components from WM and the first 5 from CSF, as in
https://doi.org/10.1016/j.neuroimage.2014.03.028
I have a few questions regarding the application of this method and would greatly appreciate any insights. Here’s an outline of what I am currently doing:
fmri_img = nilearn.image.load_img(fmri_filename)
fmri_data = fmri_img.get_fdata()
wm_signals = fmri_data[wm_mask]
wm_signals_cleaned = nilearn.signal.clean(wm_signals.T, detrend=True, standardize=‘zscore_sample’).T
wm_pca = sklearn.decomposition.PCA(n_components=5).fit_transform(wm_signals_cleaned.T)
masker = nilearn.maskers.NiftiLabelsMasker(labels_img=atlas_img, detrend=True, standardize=‘zscore_sample’, standardize_confounds=“zscore_sample”, low_pass=0.1, high_pass=0.01, t_r=2.0, ensure_finite=False, extrapolate=True)
BOLD = masker.fit_transform(fmri_img, confounds=confounds)
The confounds include motion regressors and the first 5 aCompCor components from each mask. As I understand it, the signal before PCA should not be averaged, but depending on the source, it can be detrended and filtered. Should these steps be applied to all voxels before masking, or only after, for example, applying the CSF mask? How do detrending and filtering interact when applied again (?) during signal cleaning (in my case, using a masker)? I noticed, for instance, that nipype uses linear detrending and cosine term for high-pass filtering by default, whereas in my masker, I set Butterworth filtering. Although I observed only marginal differences during experimentation, I remain curious about the nuances.
Additionally, if anyone would be kind enough to review my mask erosion procedure, I would be very grateful:
fmri_img = nilearn.image.load_img(fmri_filename)
wm_prob_img = nib.load(wm_prob_filename)
csf_prob_img = nib.load(csf_prob_filename)
wm_data = (wm_prob_img.get_fdata() > 0.9).astype(np.float32)
csf_data = (csf_prob_img.get_fdata() > 0.9).astype(np.float32)
cubic_kernel = np.ones((3, 3, 3), dtype=bool)
def process_mask(data, max_iterations, mask_type):
eroded_mask = data.copy()
last_valid_mask = None
for i in range(max_iterations):
eroded_mask = scipy.ndimage.binary_erosion(eroded_mask, structure=cubic_kernel).astype(np.float32)
resampled_mask = resample_to_img(nib.Nifti1Image(eroded_mask, wm_prob_img.affine), fmri_img, interpolation='nearest')
final_mask = np.isclose(resampled_mask.get_fdata(), 1)
final_voxel_count = np.sum(final_mask)
if final_voxel_count >= 5:
last_valid_mask = final_mask
else:
break
return last_valid_mask
wm_mask = process_mask(wm_data, max_iterations=5, mask_type="WM")
csf_mask = process_mask(csf_data, max_iterations=2, mask_type="CSF")
Thank you for all the help.