Source code for CPAC.seg_preproc.utils

# Except for functions ``_erode``, ``erosion`` and ``mask_erosion``:
# Copyright (C) 2012 - 2024 C-PAC Developers

# This file is part of C-PAC.

# C-PAC is free software: you can redistribute it and/or modify it under
# the terms of the GNU Lesser General Public License as published by the
# Free Software Foundation, either version 3 of the License, or (at your
# option) any later version.

# C-PAC is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
# License for more details.

# You should have received a copy of the GNU Lesser General Public
# License along with C-PAC. If not, see <https://www.gnu.org/licenses/>.

# STATEMENT OF CHANGES:
#     Functions in this file are derived from sources licensed under the Apache-2.0 terms,
#     and these functions have been changed.

# CHANGES:
#     * Split ``_tpm2roi`` (https://github.com/nipreps/niworkflows/blob/805cdd0/niworkflows/interfaces/utils.py#L904-L955) into ``_erode`` (common logic), ``erosion`` ("erode_out") and ``mask_erosion`` ("erode_in")
#     * Relies on mask already being generated rather than generating the mask in any of these functions
#     * Cast ``orig_vol`` as float in target volume ratio condition
#     * Modified docstrings to reflect local changes
#     * Updated style to match C-PAC codebase

# ORIGINAL WORK'S ATTRIBUTION NOTICE:
#    Copyright (c) 2016, the CRN developers team.
#    All rights reserved.

#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are met:

#    * Redistributions of source code must retain the above copyright notice, this
#      list of conditions and the following disclaimer.

#    * Redistributions in binary form must reproduce the above copyright notice,
#      this list of conditions and the following disclaimer in the documentation
#      and/or other materials provided with the distribution.

#   * Neither the name of niworkflows nor the names of its
#      contributors may be used to endorse or promote products derived from
#      this software without specific prior written permission.

#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
#    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
#    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
#    DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
#    FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
#    DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
#    SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
#    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
#    OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at

#        http://www.apache.org/licenses/LICENSE-2.0

#    Unless required by applicable law or agreed to in writing, software
#    distributed under the License is distributed on an "AS IS" BASIS,
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#    See the License for the specific language governing permissions and
#    limitations under the License.

# Modifications copyright (C) 2019 - 2024  C-PAC Developers
import os

import numpy as np
import nibabel as nib
import scipy.ndimage as nd


[docs] def check_if_file_is_empty(in_file): """Raise exception if regressor fie is empty. Parameters ---------- in_file : nii file (string) regressor file Returns ------- in_file : string return same file """ import numpy as np import nibabel as nib nii = nib.load(in_file) data = nii.get_fdata() if data.size == 0 or np.all(data == 0) or np.all(data == np.nan): msg = f"File {in_file} is empty. Use a lower threshold or turn off regressors." raise ValueError(msg) return in_file
def _erode(roi_mask, erosion_mm, erosion_prop): """Perform in-common erosion steps. Parameters ---------- mask : str Path to mask erosion_mm : int, float or None Kernel width in mm erosion_prop : float or None Target volume ratio, 0 < erosion_prop < 1 Returns ------- mask_img : nibabel image Original mask image erode : bool erosion_mm or erosion_prop mask_data : numpy array eroded mask data """ # This functionality is adapted from poldracklab/niworkflows: # https://github.com/nipreps/niworkflows/blob/805cdd0/niworkflows/interfaces/utils.py#L904-L955 # https://fmriprep.readthedocs.io/ # https://poldracklab.stanford.edu/ # We are temporarily maintaining our own copy for more granular control. mask_img = nib.load(roi_mask) mask_data = mask_img.get_fdata() erode = (erosion_mm is not None and erosion_mm > 0) or ( erosion_prop is not None and 0 < erosion_prop < 1 ) if erode: if erosion_mm: iter_n = max(int(erosion_mm / max(mask_img.header.get_zooms())), 1) mask_data = nd.binary_erosion(mask_data, iterations=iter_n) else: orig_vol = np.sum(mask_data > 0) while np.sum(mask_data > 0) / (orig_vol * 1.0) > erosion_prop: mask_data = nd.binary_erosion(mask_data, iterations=1) return mask_img, erode, mask_data
[docs] def pick_wm_prob_0(probability_maps): """Returns the csf probability map from the list of segmented probability maps. Parameters ---------- probability_maps : list (string) List of Probability Maps Returns ------- file : string Path to segment_prob_0.nii.gz is returned """ if isinstance(probability_maps, list): if len(probability_maps) == 1: probability_maps = probability_maps[0] for filename in probability_maps: if filename.endswith("prob_0.nii.gz"): return filename return None
[docs] def pick_wm_prob_1(probability_maps): """Returns the gray matter probability map from the list of segmented probability maps. Parameters ---------- probability_maps : list (string) List of Probability Maps Returns ------- file : string Path to segment_prob_1.nii.gz is returned """ if isinstance(probability_maps, list): if len(probability_maps) == 1: probability_maps = probability_maps[0] for filename in probability_maps: if filename.endswith("prob_1.nii.gz"): return filename return None
[docs] def pick_wm_prob_2(probability_maps): """Returns the white matter probability map from the list of segmented probability maps. Parameters ---------- probability_maps : list (string) List of Probability Maps Returns ------- file : string Path to segment_prob_2.nii.gz is returned """ if isinstance(probability_maps, list): if len(probability_maps) == 1: probability_maps = probability_maps[0] for filename in probability_maps: if filename.endswith("prob_2.nii.gz"): return filename return None
[docs] def pick_wm_class_0(tissue_class_files): """Returns the csf tissue class file from the list of segmented tissue class files. Parameters ---------- tissue_class_files : list (string) List of tissue class files Returns ------- file : string Path to segment_seg_0.nii.gz is returned """ if isinstance(tissue_class_files, list): if len(tissue_class_files) == 1: tissue_class_files = tissue_class_files[0] for filename in tissue_class_files: if filename.endswith("seg_0.nii.gz"): return filename return None
[docs] def pick_wm_class_1(tissue_class_files): """Returns the gray matter tissue class file from the list of segmented tissue class files. Parameters ---------- tissue_class_files : list (string) List of tissue class files Returns ------- file : string Path to segment_seg_1.nii.gz is returned """ if isinstance(tissue_class_files, list): if len(tissue_class_files) == 1: tissue_class_files = tissue_class_files[0] for filename in tissue_class_files: if filename.endswith("seg_1.nii.gz"): return filename return None
[docs] def pick_wm_class_2(tissue_class_files): """Returns the white matter tissue class file from the list of segmented tissue class files. Parameters ---------- tissue_class_files : list (string) List of tissue class files Returns ------- file : string Path to segment_seg_2.nii.gz is returned """ if isinstance(tissue_class_files, list): if len(tissue_class_files) == 1: tissue_class_files = tissue_class_files[0] for filename in tissue_class_files: if filename.endswith("seg_2.nii.gz"): return filename return None
[docs] def mask_erosion( roi_mask=None, skullstrip_mask=None, mask_erosion_mm=None, mask_erosion_prop=None ): """Returns eroded segment mask and skull-stripped brain mask. Parameters ---------- roi_mask : string Path to binarized segment mask skullstrip_mask : string Path to skull-stripped brain mask mask_erosion_prop : float Target volume ratio for skull-stripped brain mask Returns ------- output_roi_mask : string Path to eroded segment mask eroded_skullstrip_mask : string Path to eroded skull-stripped brain mask """ # This functionality is adapted from poldracklab/niworkflows: # https://github.com/nipreps/niworkflows/blob/805cdd0/niworkflows/interfaces/utils.py#L916-L935 # https://fmriprep.readthedocs.io/ # https://poldracklab.stanford.edu/ # We are temporarily maintaining our own copy for more granular control. roi_mask_img = nib.load(roi_mask) roi_mask_data = roi_mask_img.get_fdata() skullstrip_mask_img, erode_in, skullstrip_mask_data = _erode( skullstrip_mask, mask_erosion_mm, mask_erosion_prop ) if erode_in: # pylint: disable=invalid-unary-operand-type roi_mask_data[~skullstrip_mask_data] = 0 hdr = roi_mask_img.header output_roi_mask_img = nib.Nifti1Image( roi_mask_data, header=hdr, affine=roi_mask_img.affine ) output_roi_mask = os.path.join(os.getcwd(), "segment_tissue_eroded_mask.nii.gz") output_roi_mask_img.to_filename(output_roi_mask) hdr = skullstrip_mask_img.header output_skullstrip_mask_img = nib.Nifti1Image( skullstrip_mask_data, header=hdr, affine=skullstrip_mask_img.affine ) eroded_skullstrip_mask = os.path.join(os.getcwd(), "eroded_skullstrip_mask.nii.gz") output_skullstrip_mask_img.to_filename(eroded_skullstrip_mask) return output_roi_mask, eroded_skullstrip_mask
[docs] def erosion(roi_mask=None, erosion_mm=None, erosion_prop=None): """Returns eroded tissue segment mask. Parameters ---------- roi_mask : string Path to binarized segment (ROI) mask erosion_prop : float Target volume ratio for erosion segment mask Returns ------- eroded_roi_mask : string Path to eroded segment mask """ # This functionality is adapted from poldracklab/niworkflows: # https://github.com/nipreps/niworkflows/blob/805cdd0/niworkflows/interfaces/utils.py#L937-L954 # https://fmriprep.readthedocs.io/ # https://poldracklab.stanford.edu/ # We are temporarily maintaining our own copy for more granular control. roi_mask_img, _, roi_mask_data = _erode(roi_mask, erosion_mm, erosion_prop) hdr = roi_mask_img.header output_img = nib.Nifti1Image(roi_mask_data, header=hdr, affine=roi_mask_img.affine) eroded_roi_mask = os.path.join(os.getcwd(), "segment_tissue_mask.nii.gz") output_img.to_filename(eroded_roi_mask) return eroded_roi_mask
[docs] def hardcoded_antsJointLabelFusion( anatomical_brain, anatomical_brain_mask, template_brain_list, template_segmentation_list, ): """Run antsJointLabelFusion.sh. Parameters ---------- anatomical_brain : string (nifti file) Target image to be labeled. anatomical_brain_mask : string (nifti file) Target mask image template_brain_list : list Atlas to be warped to target image. template_segmentation_list : list Labels corresponding to atlas. Returns ------- multiatlas_Intensity : string (nifti file) multiatlas_Labels : string (nifti file) """ import os import subprocess cmd = ["${ANTSPATH}${ANTSPATH:+/}antsJointLabelFusion.sh"] cmd.append( f" -d 3 -o ants_multiatlas_ -t {anatomical_brain} -x {anatomical_brain_mask} -y b -c 0" ) if not len(template_brain_list) == len(template_segmentation_list): err_msg = ( "\n\n[!] C-PAC says: Please check ANTs Prior-based " "Segmentation setting. For performing ANTs Prior-based " "segmentation method the number of specified " "segmentations should be identical to the number of atlas " "image sets.\n\n" ) raise Exception(err_msg) else: for index in range(len(template_brain_list)): cmd.append( f" -g {template_brain_list[index]} -l {template_segmentation_list[index]}" ) # write out the actual command-line entry for testing/validation later command_file = os.path.join(os.getcwd(), "command.txt") with open(command_file, "wt") as f: f.write(" ".join(cmd)) str = "" bash_cmd = str.join(cmd) try: subprocess.check_output(bash_cmd, shell=True) # pylint: disable=unused-variable except Exception as e: # pylint: disable=broad-except,invalid-name # pylint: disable=raise-missing-from msg = ( "[!] antsJointLabel segmentation method did not " "complete successfully.\n\nError " "details:\n{0}\n{1}\n".format(e, getattr(e, "output", "")) ) raise Exception(msg) multiatlas_Intensity = None multiatlas_Labels = None files = [f for f in os.listdir(".") if os.path.isfile(f)] for f in files: if "Intensity" in f: multiatlas_Intensity = os.getcwd() + "/" + f if "Labels" in f: multiatlas_Labels = os.getcwd() + "/" + f if not multiatlas_Labels: msg = ( "\n\n[!] No multiatlas labels file found. " "antsJointLabelFusion may not have completed " "successfully.\n\n" ) raise Exception(msg) return multiatlas_Intensity, multiatlas_Labels
def pick_tissue_from_labels_file( multiatlas_Labels, csf_label=[4, 14, 15, 24, 43], gm_label=[3, 42], wm_label=[2, 41] ): """Pick tissue mask from multiatlas labels file based off of FreeSurferColorLUT https://surfer.nmr.mgh.harvard.edu/fswiki/FsTutorial/AnatomicalROI/FreeSurferColorLUT or user provided label value. Parameters ---------- multiatlas_Labels : string (nifti file) csf_label : list a list of integer label values corresponding to CSF in multiatlas file gm_label : list a list of integer label value corresponding to Gray Matter in multiatlas file wm_label : list a list of integer label value corresponding to White Matter in multiatlas file Returns ------- csf_mask : string (nifti file) gm_mask : string (nifti file) wm_mask : string (nifti file) """ # pylint: disable=import-outside-toplevel,redefined-outer-name,reimported import os import numpy as np import nibabel as nib img = nib.load(multiatlas_Labels) data = img.get_fdata() # pick tissue mask from multiatlas labels file # based off of FreeSurferColorLUT or user provided label values # hard-coded csf/gm/wm label values are based off of FreeSurferColorLUT # FreeSurfer Ventricle Labels: # Left-Lateral-Ventricle 4, 3rd-Ventricle 14, 4th-Ventricle 15, Right-Lateral-Ventricle 43 csf = np.zeros(np.size(data)) csf[np.where(np.in1d(data, np.array(csf_label)))] = 1 csf = csf.reshape(data.shape) gm = np.zeros(np.size(data)) gm[np.where(np.in1d(data, np.array(gm_label)))] = 1 gm = gm.reshape(data.shape) wm = np.zeros(np.size(data)) wm[np.where(np.in1d(data, np.array(wm_label)))] = 1 wm = wm.reshape(data.shape) save_img_csf = nib.Nifti1Image(csf, header=img.header, affine=img.affine) save_img_gm = nib.Nifti1Image(gm, header=img.header, affine=img.affine) save_img_wm = nib.Nifti1Image(wm, header=img.header, affine=img.affine) save_img_csf.to_filename("csf_mask.nii.gz") save_img_gm.to_filename("gm_mask.nii.gz") save_img_wm.to_filename("wm_mask.nii.gz") csf_mask = os.path.join(os.getcwd(), "csf_mask.nii.gz") gm_mask = os.path.join(os.getcwd(), "gm_mask.nii.gz") wm_mask = os.path.join(os.getcwd(), "wm_mask.nii.gz") return csf_mask, gm_mask, wm_mask