Merge pull request #321 from CoBrALab/fix_smoothing

After smoothing, the image is divided by the smoothed mask to correct…

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
 

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)

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)