From e281458284d56bc339d7df345248fcf795b9ca11 Mon Sep 17 00:00:00 2001 From: Gab-D-G Date: Thu, 28 Sep 2023 14:01:14 -0400 Subject: [PATCH 1/2] After smoothing, the image is divided by the smoothed mask to correct for edge effects, and the mask is applied to remove values outside the brain. Closes #320 --- .../confound_correction.py | 3 ++- rabies/confound_correction_pkg/utils.py | 18 +++++++++++++++++- 2 files changed, 19 insertions(+), 2 deletions(-) diff --git a/rabies/confound_correction_pkg/confound_correction.py b/rabies/confound_correction_pkg/confound_correction.py index 0dd704c..bf25942 100755 --- a/rabies/confound_correction_pkg/confound_correction.py +++ b/rabies/confound_correction_pkg/confound_correction.py @@ -490,7 +490,8 @@ def _run_interface(self, runtime): ''' import nibabel as nb affine = nb.load(bold_file).affine[:3,:3] # still not sure how to match nibabel's affine reliably - timeseries_img = smooth_image(timeseries_img, affine, cr_opts.smoothing_filter) + mask_img = sitk.ReadImage(brain_mask_file, sitk.sitkFloat32) + timeseries_img = smooth_image(timeseries_img, affine, cr_opts.smoothing_filter, mask_img) cleaned_path = cr_out+'/'+filename_split[0]+'_cleaned.nii.gz' sitk.WriteImage(timeseries_img, cleaned_path) diff --git a/rabies/confound_correction_pkg/utils.py b/rabies/confound_correction_pkg/utils.py index 4566c85..f0dabc3 100644 --- a/rabies/confound_correction_pkg/utils.py +++ b/rabies/confound_correction_pkg/utils.py @@ -447,7 +447,7 @@ def phase_randomized_regressors(confounds_array, frame_mask, TR): return randomized_confounds_array -def smooth_image(img, affine, fwhm): +def smooth_image(img, affine, fwhm, mask_img): # apply nilearn's Gaussian smoothing on a SITK image from nilearn.image.image import _smooth_array from rabies.utils import copyInfo_4DImage, copyInfo_3DImage @@ -466,12 +466,28 @@ def smooth_image(img, affine, fwhm): array_3d = sitk.GetArrayFromImage(img) arr = array_3d.transpose(2,1,0) # re-orient the array to match the affine smoothed_arr = _smooth_array(arr, affine, fwhm=fwhm, ensure_finite=True, copy=True) + + # smoothing creates leakage around mask boundaries + # correct for edge effects by dividing by the smoothed mask like FSL https://johnmuschelli.com/fslr/reference/fslsmooth.html + mask_arr = sitk.GetArrayFromImage(mask_img).transpose(2,1,0) # re-orient the array to match the affine + smoothed_mask = _smooth_array(mask_arr, affine, fwhm=fwhm, ensure_finite=True, copy=True) + smoothed_mask = smoothed_mask.transpose(2,1,0) # recover SITK orientation + + # recover SITK img if dim==4: smoothed_arr = smoothed_arr.transpose(3,2,1,0) + smoothed_arr /= smoothed_mask # correct for edge effect + smoothed_arr *= sitk.GetArrayFromImage(mask_img) # re-apply mask + smoothed_arr[np.isnan(smoothed_arr)] = 0 + smoothed_img = copyInfo_4DImage(sitk.GetImageFromArray( smoothed_arr, isVector=False), img, img) elif dim==3: smoothed_arr = smoothed_arr.transpose(2,1,0) + smoothed_arr /= smoothed_mask # correct for edge effect + smoothed_arr *= sitk.GetArrayFromImage(mask_img) # re-apply mask + smoothed_arr[np.isnan(smoothed_arr)] = 0 + smoothed_img = copyInfo_3DImage(sitk.GetImageFromArray( smoothed_arr, isVector=False), img) From a9b06dee7be807b37ea8f9ddaf8926f851dcc81b Mon Sep 17 00:00:00 2001 From: Gab-D-G Date: Thu, 28 Sep 2023 14:33:46 -0400 Subject: [PATCH 2/2] smoothing takes the mask as input --- rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py b/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py index 46718d6..46fb1df 100644 --- a/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py +++ b/rabies/analysis_pkg/diagnosis_pkg/analysis_QC.py @@ -57,10 +57,11 @@ def get_maps(prior, prior_list, corr_variable, mask_file, smoothing=False, non_p maps.append(corr_map) if smoothing: + mask_img = sitk.ReadImage(mask_file) import nibabel as nb affine = nb.load(mask_file).affine[:3,:3] for i in range(len(maps)): - maps[i] = sitk.GetArrayFromImage(smooth_image(recover_3D(mask_file,maps[i]), affine, 0.3))[volume_indices] + maps[i] = sitk.GetArrayFromImage(smooth_image(recover_3D(mask_file,maps[i]), affine, 0.3, mask_img))[volume_indices] return maps