|
| 1 | +"""Utilities for confounds manipulation.""" |
| 2 | + |
| 3 | + |
| 4 | +def mask2vf(in_file, zooms=None, out_file=None): |
| 5 | + """ |
| 6 | + Convert a binary mask on a volume fraction map. |
| 7 | +
|
| 8 | + The algorithm simply applies a Gaussian filter with the kernel size scaled |
| 9 | + by the zooms given as argument. |
| 10 | +
|
| 11 | + """ |
| 12 | + import numpy as np |
| 13 | + import nibabel as nb |
| 14 | + from scipy.ndimage import gaussian_filter |
| 15 | + |
| 16 | + img = nb.load(in_file) |
| 17 | + imgzooms = np.array(img.header.get_zooms()[:3], dtype=float) |
| 18 | + if zooms is None: |
| 19 | + zooms = imgzooms |
| 20 | + |
| 21 | + zooms = np.array(zooms, dtype=float) |
| 22 | + sigma = 0.5 * (zooms / imgzooms) |
| 23 | + |
| 24 | + data = gaussian_filter(img.get_fdata(dtype=np.float32), sigma=sigma) |
| 25 | + |
| 26 | + max_data = np.percentile(data[data > 0], 99) |
| 27 | + data = np.clip(data / max_data, a_min=0, a_max=1) |
| 28 | + |
| 29 | + if out_file is None: |
| 30 | + return data |
| 31 | + |
| 32 | + hdr = img.header.copy() |
| 33 | + hdr.set_data_dtype(np.float32) |
| 34 | + nb.Nifti1Image(data.astype(np.float32), img.affine, hdr).to_filename(out_file) |
| 35 | + return out_file |
| 36 | + |
| 37 | + |
| 38 | +def acompcor_masks(in_files, is_aseg=False, zooms=None): |
| 39 | + """ |
| 40 | + Generate aCompCor masks. |
| 41 | +
|
| 42 | + This function selects the CSF partial volume map from the input, |
| 43 | + and generates the WM and combined CSF+WM masks for aCompCor. |
| 44 | +
|
| 45 | + The implementation deviates from Behzadi et al. |
| 46 | + Their original implementation thresholded the CSF and the WM partial-volume |
| 47 | + masks at 0.99 (i.e., 99% of the voxel volume is filled with a particular tissue), |
| 48 | + and then binary eroded that 2 voxels: |
| 49 | +
|
| 50 | + > Anatomical data were segmented into gray matter, white matter, |
| 51 | + > and CSF partial volume maps using the FAST algorithm available |
| 52 | + > in the FSL software package (Smith et al., 2004). Tissue partial |
| 53 | + > volume maps were linearly interpolated to the resolution of the |
| 54 | + > functional data series using AFNI (Cox, 1996). In order to form |
| 55 | + > white matter ROIs, the white matter partial volume maps were |
| 56 | + > thresholded at a partial volume fraction of 0.99 and then eroded by |
| 57 | + > two voxels in each direction to further minimize partial voluming |
| 58 | + > with gray matter. CSF voxels were determined by first thresholding |
| 59 | + > the CSF partial volume maps at 0.99 and then applying a threedimensional |
| 60 | + > nearest neighbor criteria to minimize multiple tissue |
| 61 | + > partial voluming. Since CSF regions are typically small compared |
| 62 | + > to white matter regions mask, erosion was not applied. |
| 63 | +
|
| 64 | + This particular procedure is not generalizable to BOLD data with different voxel zooms |
| 65 | + as the mathematical morphology operations will be scaled by those. |
| 66 | + Also, from reading the excerpt above and the tCompCor description, I (@oesteban) |
| 67 | + believe that they always operated slice-wise given the large slice-thickness of |
| 68 | + their functional data. |
| 69 | +
|
| 70 | + Instead, *fMRIPrep*'s implementation deviates from Behzadi's implementation on two |
| 71 | + aspects: |
| 72 | +
|
| 73 | + * the masks are prepared in high-resolution, anatomical space and then |
| 74 | + projected into BOLD space; and, |
| 75 | + * instead of using binary erosion, a dilated GM map is generated -- thresholding |
| 76 | + the corresponding PV map at 0.05 (i.e., pixels containing at least 5% of GM tissue) |
| 77 | + and then subtracting that map from the CSF, WM and CSF+WM (combined) masks. |
| 78 | + This should be equivalent to eroding the masks, except that the erosion |
| 79 | + only happens at direct interfaces with GM. |
| 80 | +
|
| 81 | + When the probseg maps provene from FreeSurfer's ``recon-all`` (i.e., they are |
| 82 | + discrete), binary maps are *transformed* into some sort of partial volume maps |
| 83 | + by means of a Gaussian smoothing filter with sigma adjusted by the size of the |
| 84 | + BOLD data. |
| 85 | +
|
| 86 | + """ |
| 87 | + from pathlib import Path |
| 88 | + import numpy as np |
| 89 | + import nibabel as nb |
| 90 | + from scipy.ndimage import binary_dilation |
| 91 | + from skimage.morphology import ball |
| 92 | + |
| 93 | + csf_file = in_files[2] # BIDS labeling (CSF=2; last of list) |
| 94 | + # Load PV maps (fast) or segments (recon-all) |
| 95 | + gm_vf = nb.load(in_files[0]) |
| 96 | + wm_vf = nb.load(in_files[1]) |
| 97 | + csf_vf = nb.load(csf_file) |
| 98 | + |
| 99 | + # Prepare target zooms |
| 100 | + imgzooms = np.array(gm_vf.header.get_zooms()[:3], dtype=float) |
| 101 | + if zooms is None: |
| 102 | + zooms = imgzooms |
| 103 | + zooms = np.array(zooms, dtype=float) |
| 104 | + |
| 105 | + if not is_aseg: |
| 106 | + gm_data = gm_vf.get_fdata() > 0.05 |
| 107 | + wm_data = wm_vf.get_fdata() |
| 108 | + csf_data = csf_vf.get_fdata() |
| 109 | + else: |
| 110 | + csf_file = mask2vf( |
| 111 | + csf_file, |
| 112 | + zooms=zooms, |
| 113 | + out_file=str(Path("acompcor_csf.nii.gz").absolute()), |
| 114 | + ) |
| 115 | + csf_data = nb.load(csf_file).get_fdata() |
| 116 | + wm_data = mask2vf(in_files[1], zooms=zooms) |
| 117 | + |
| 118 | + # We do not have partial volume maps (recon-all route) |
| 119 | + gm_data = np.asanyarray(gm_vf.dataobj, np.uint8) > 0 |
| 120 | + |
| 121 | + # Dilate the GM mask |
| 122 | + gm_data = binary_dilation(gm_data, structure=ball(3)) |
| 123 | + |
| 124 | + # Output filenames |
| 125 | + wm_file = str(Path("acompcor_wm.nii.gz").absolute()) |
| 126 | + combined_file = str(Path("acompcor_wmcsf.nii.gz").absolute()) |
| 127 | + |
| 128 | + # Prepare WM mask |
| 129 | + wm_data[gm_data] = 0 # Make sure voxel does not contain GM |
| 130 | + nb.Nifti1Image(wm_data, gm_vf.affine, gm_vf.header).to_filename(wm_file) |
| 131 | + |
| 132 | + # Prepare combined CSF+WM mask |
| 133 | + comb_data = csf_data + wm_data |
| 134 | + comb_data[gm_data] = 0 # Make sure voxel does not contain GM |
| 135 | + nb.Nifti1Image(comb_data, gm_vf.affine, gm_vf.header).to_filename(combined_file) |
| 136 | + return [csf_file, wm_file, combined_file] |
0 commit comments