Merge 31455dd into 6170cf3

setup.py

@@ -56,6 +56,7 @@
         "raven",
         "joblib",
         "quadprog",
+        "spec2nii",
     ],
     extras_require={
         'docs': ["sphinx>=1.7", "sphinx_rtd_theme>=1.2.2", "sphinx-click"],

shimmingtoolbox/cli/mask.py

@@ -8,6 +8,7 @@
 
 from shimmingtoolbox.masking.shapes import shape_square, shape_cube, shape_sphere
 import shimmingtoolbox.masking.threshold
+from shimmingtoolbox.masking.mask_mrs import mask_mrs
 from shimmingtoolbox.utils import run_subprocess, create_output_dir, set_all_loggers
 
 CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
@@ -303,6 +304,41 @@ def sct(fname_input, fname_output, contrast, centerline, file_centerline, brain,
     return fname_output
 
 
+@mask_cli.command(context_settings=CONTEXT_SETTINGS,
+                  help="Create a mask to shim single voxel MRS. "
+                       "Voxel position and size can be directly given or these info can be read "
+                       "from the twix raw-data. "
+                       "The mask is stored by default under the name 'mask_mrs.nii.gz' in the output "
+                       "folder. Return an output nifti file to be used as a mask for MRS shimming.")
+@click.option('-i', '--input', 'fname_input', type=click.Path(), required=True,
+               help="Input path of the fieldmap to be shimmed.")
+@click.option('-r', '--raw_data', type=click.Path(),
+              help="Input path of the of the twix raw-data (supported extention .dat)")
+@click.option('-o', '--output', type=click.Path(), default=os.path.join(os.curdir, 'mask_mrs.nii.gz'),
+              show_default=True, help="Name of the output mask. Supported extensions are .nii or .nii.gz. (default: "
+                   "(os.curdir, 'mask_mrs.nii.gz'))")
+@click.option('-c', '--center', nargs=3, type=click.FLOAT, help="Voxel's center position in mm of the x, y and z of "
+              "the scanner's coordinate")
+@click.option('-s', '--size', nargs=3, type=click.FLOAT, help="Voxel size in mm of the x, y and z of the scanner's "
+              "coordinate")
+@click.option('--verbose', type=click.Choice(['info', 'debug']), default='info', help="Be more verbose")
+def mrs(fname_input, output, raw_data, center, size, verbose):
+
+    # Set all loggers
+    set_all_loggers(verbose)
+
+    # Prepare the output
+    create_output_dir(output, is_file=True)
+
+    nii = nib.load(fname_input)
+    output_mask = mask_mrs(fname_input, raw_data, center, size) # creation of the MRS mask
+    output_mask = output_mask.astype(np.int32)
+    nii_img = nib.Nifti1Image(output_mask, nii.affine, header=nii.header)
+    nib.save(nii_img, output)
+    click.echo(f"The filename for the output mask is: {os.path.abspath(output)}")
+    return output
+
+
 # def _get_centerline(fname_process, fname_output, method='optic', contrast='t2', centerline_algo='bspline',
 #                     centerline_smooth='30', verbose='1'):
 #     """ Wrapper to sct_get_centerline. Allows to get the centerline of the spinal cord and outputs a nifti file

shimmingtoolbox/masking/mask_mrs.py

@@ -0,0 +1,83 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+"""
+Creating MRS mask API
+"""
+import numpy as np
+import numpy.linalg as npl
+import nibabel as nib
+import logging
+from shimmingtoolbox.utils import splitext
+from shimmingtoolbox.utils import run_subprocess
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+
+def mask_mrs(fname_input, raw_data, center, size):
+    """
+    Create a mask to shim single voxel MRS
+
+    Args:
+        fname_input (str): Input path of the fieldmap to be shimmed (supported extention .nii and .nii.gz)
+        raw_data (str): Input path of the of the twix raw-data (supported extention .dat)
+        center (float): voxel's center position in mm of the x, y and z of the scanner's coordinate
+        size (float): voxel size in mm of the x, y and z of the scanner's coordinate
+    Returns:
+        numpy.ndarray: Cubic mask with same dimensions as MRS voxel.
+    """
+
+    if fname_input is None:
+        raise TypeError(f"The fname_input is cannot be None. See: st_mask mrs-mask -h")
+
+    if raw_data is None:
+        logger.info(" The raw_data is not provided, creating the mask with the given voxel position and size info")
+        if center is None or size is None:
+             raise TypeError('The raw_data is not provided; the voxel position and size are required to proceed')
+        else:
+            mrs_voxel_size = size
+            logger.debug('center:', center)
+            scanner_coordinate = np.array(center + (1,))
+            logger.debug(f"Scanner position is: {scanner_coordinate}")
+
+    else:
+        logger.info("Reading the twix raw_data")
+        run_subprocess(['spec2nii', 'twix', raw_data, '-e', 'image'])
+        name_nii, ext = splitext(raw_data)
+        nii = nib.load(name_nii + '.nii.gz')
+        header_twix = nii.header
+        affine = nii.affine
+        position_sag = header_twix['qoffset_x']
+        position_cor = header_twix['qoffset_y']
+        position_tra = header_twix['qoffset_z']
+        logger.debug(f"twix header: {header_twix}")
+        logger.debug(f"affine: {affine}")
+        scanner_coordinate = np.array([position_sag, position_cor, position_tra, 1])
+        logger.debug(f"Scanner position is: {scanner_coordinate}")
+        mrs_voxel_size = header_twix['pixdim'][1:4]
+
+    logger.debug('mrs_voxel_size is:', mrs_voxel_size)
+    fmap_nii = nib.load(fname_input)
+    fmap_array = fmap_nii.get_fdata()
+    fmap_affine = fmap_nii.affine
+    fmap_header = fmap_nii.header
+    logger.debug('reference_affine:', fmap_affine)
+    logger.debug('reference_affine shape:', np.shape(fmap_affine))
+    voxel_position = npl.inv(fmap_affine).dot(scanner_coordinate)
+    voxel_position = np.round(voxel_position)
+    logger.debug('voxel_position:', voxel_position)
+
+    # The given coordinate (i, j, k) is the voxel's center position.
+    i, j, k = map(int, voxel_position[:3])
+    fmap_resolution = fmap_header['pixdim'][1:4]
+
+    # The distance from the center of the MRS voxel to it's edges is calculated based on number of fieldmap voxels.
+    half_voxel_distances = np.ceil(np.divide(mrs_voxel_size, 2 * fmap_resolution)).astype(int)
+    sd = half_voxel_distances
+
+    # create a zero mask with the same size as the input fieldmap to be shimmed.
+    mask = np.zeros(fmap_array.shape)
+
+    # Change the MRS voxel position to have 1 value.
+    mask[i - sd[0]:i + sd[0], j - sd[1]:j + sd[1], k - sd[2]:k + sd[2]] = 1
+    return mask