Merge 291638a into 4013cae

setup.py

@@ -29,7 +29,8 @@
             "st_dump_env_info=shimmingtoolbox.cli.check_env:dump_env_info",
             "st_image=shimmingtoolbox.cli.image:image_cli",
             "st_maths=shimmingtoolbox.cli.maths:maths_cli",
-            "st_b0shim=shimmingtoolbox.cli.b0shim:b0shim_cli"
+            "st_b0shim=shimmingtoolbox.cli.b0shim:b0shim_cli",
+            "st_create_coil_profiles=shimmingtoolbox.cli.create_coil_profiles:create_coil_profiles_cli"
         ]
     },
     packages=find_packages(exclude=["docs"]),

shimmingtoolbox/cli/create_coil_profiles.py

@@ -0,0 +1,223 @@
+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+
+import click
+import json
+import os
+import logging
+import nibabel as nib
+import numpy as np
+
+from shimmingtoolbox.coils.create_coil_profiles import create_coil_profiles
+from shimmingtoolbox.cli.prepare_fieldmap import prepare_fieldmap_uncli
+from shimmingtoolbox.utils import create_output_dir, save_nii_json, set_all_loggers
+from shimmingtoolbox.masking.threshold import threshold as mask_threshold
+
+CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help'])
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+
+@click.command(
+    context_settings=CONTEXT_SETTINGS,
+)
+@click.option('-i', '--input', 'fname_json', type=click.Path(exists=True), required=True,
+              help="Input filename of json config file")
+@click.option('--unwrapper', type=click.Choice(['prelude']), default='prelude', show_default=True,
+              help="Algorithm for unwrapping")
+@click.option('--threshold', type=float, help="Threshold for masking.")
+@click.option('--autoscale-phase', 'autoscale', type=click.BOOL, default=True, show_default=True,
+              help="Tells whether to auto rescale phase inputs according to manufacturer standards. If you have non "
+                   "standard data, it would be preferable to set this option to False and input your phase data from "
+                   "-pi to pi to avoid unwanted rescaling")
+@click.option('--gaussian-filter', 'gaussian_filter', type=bool, show_default=True, help="Gaussian filter for B0 maps")
+@click.option('--sigma', type=float, default=1, help="Standard deviation of gaussian filter. Used for: gaussian_filter")
+@click.option('-o', '--output', 'fname_output', type=click.Path(), required=False,
+              default=os.path.join(os.path.curdir, 'coil_profiles.nii.gz'),
+              help="Output path filename of coil profile nifti file. Supported types : '.nii', '.nii.gz'")
+@click.option('-v', '--verbose', type=click.Choice(['info', 'debug']), default='info', help="Be more verbose")
+def create_coil_profiles_cli(fname_json, autoscale, unwrapper, threshold, gaussian_filter, sigma, fname_output,
+                             verbose):
+    """Create b0 coil profiles from acquisitions defined in the input json file"""
+
+    # Set logger level
+    set_all_loggers(verbose)
+
+    # Config file is set up:
+    # phase[i_channel][min_max][i_echo]
+    # mag[i_channel][min_max][i_echo]
+    #     "setup_currents": [
+    #     [-0.5, 0.5],
+    #     [-0.5, 0.5],
+    #     [-0.5, 0.5],
+    #     [-0.5, 0.5],
+    #     [-0.5, 0.5],
+    #     [-0.5, 0.5],
+    #     [-0.5, 0.5]
+    #   ],
+    #   "name": "greg_coil",
+    #   "n_channels": 7,
+    #   "units": "A",
+    #   "coef_channel_minmax": [i_channel][min_max]
+    #   "coef_sum_max": null
+
+    # Get directory
+    path_output = os.path.dirname(os.path.abspath(fname_output))
+    create_output_dir(path_output)
+
+    # Open json config file
+    with open(fname_json) as json_file:
+        json_data = json.load(json_file)
+
+    # Init variables
+    phases = json_data["phase"]
+    mags = json_data["mag"]
+    list_setup_currents = json_data["setup_currents"]
+    n_channels = json_data["n_channels"]
+
+    dead_channels = []
+    for i_channel in range(n_channels):
+        if not phases[i_channel][0]:
+            dead_channels.append(i_channel)
+            continue
+
+    # Create a mask containing the threshold of all channels and currents
+    fname_mask = os.path.join(path_output, 'mask.nii.gz')
+    fname_mag = mags[0][0][0]
+    nii_mag = nib.load(fname_mag)
+    mask = np.full_like(nii_mag.get_fdata(), True, bool)
+
+    for i_channel in range(n_channels):
+
+        # If channel is dead, dont include it in the mask
+        if i_channel in dead_channels:
+            continue
+
+        n_currents = len(phases[i_channel])
+        channel_mask = np.full_like(nii_mag.get_fdata(), True, bool)
+        for i_current in range(n_currents):
+            n_echoes = len(phases[i_channel][i_current])
+            # Calculate the average mag image
+            current_mean = np.zeros_like(nii_mag.get_fdata())
+            for i_echo in range(n_echoes):
+                fname_mag = mags[i_channel][i_current][i_echo]
+                mag = nib.load(fname_mag).get_fdata()
+                current_mean += mag
+
+            # Threshold mask for i_channel, i_current and all echoes
+            current_mean /= n_echoes
+            tmp_mask = mask_threshold(current_mean, threshold)
+
+            # And mask for a i_channel but all currents
+            channel_mask = np.logical_and(tmp_mask, channel_mask)
+
+        # And mask for all channels
+        mask = np.logical_and(channel_mask, mask)
+
+    # Mask contains the region where all channels get enough signal
+    nii_mask = nib.Nifti1Image(mask.astype(int), nii_mag.affine, header=nii_mag.header)
+    nib.save(nii_mask, fname_mask)
+
+    # In 4d
+    nii_mask_4d = nib.Nifti1Image(np.repeat(mask.astype(int)[..., np.newaxis], n_currents, -1), nii_mag.affine,
+                                  header=nii_mag.header)
+    fname_mask_4d = os.path.join(path_output, f"mask_4d.nii.gz")
+    nib.save(nii_mask_4d, fname_mask_4d)
+
+    if not dead_channels:
+        logger.warning(f"Channels: {dead_channels} do not have phase data. They will be set to 0.")
+
+    fnames_fmap = []
+    index_channel = -1
+    # For each channel
+    for i_channel in range(n_channels):
+        if i_channel in dead_channels:
+            continue
+
+        # Keeps track of the index since there could be dead channels
+        index_channel += 1
+
+        n_currents = len(phases[i_channel])
+        # Method 1: Unwrap each current individually
+        # for each current
+        fnames_fmap.append([])
+        for i_current in range(n_currents):
+            phase = phases[i_channel][i_current]
+
+            # Calculate fieldmap and save to a file
+            fname_fmap = os.path.join(path_output, f"channel{i_channel}_{i_current}_fieldmap.nii.gz")
+            prepare_fieldmap_uncli(phase, fname_mag, unwrapper, fname_fmap, autoscale,
+                                   fname_mask=fname_mask,
+                                   gaussian_filter=gaussian_filter,
+                                   sigma=sigma)
+
+            fnames_fmap[index_channel].append(fname_fmap)
+
+    # Remove dead channels from the list of currents
+    for i_channel in dead_channels:
+        list_setup_currents.pop(i_channel)
+
+    # Create coil profiles
+    profiles = create_coil_profiles(fnames_fmap, list_currents=list_setup_currents)
+
+    # Add dead channels as 0s
+    for i_dead_channel in dead_channels:
+        profiles = np.insert(profiles, i_dead_channel, np.zeros(profiles.shape[:3]), axis=3)
+
+    # If not debug, remove junk output
+    if not logger.level <= getattr(logging, 'DEBUG'):
+        os.remove(fname_mask)
+        # For each channel
+        for list_fnames in fnames_fmap:
+            # For each fieldmap
+            for fname_nifti in list_fnames:
+                fname_json_m = fname_nifti.rsplit('.nii', 1)[0] + '.json'
+                # Delete nifti
+                # os.remove(fname_nifti)
+                # Delete json
+                # os.remove(fname_json_m)
+
+    # Use header and json info from first file in the list
+    fname_json_phase = phases[0][0][0].rsplit('.nii', 1)[0] + '.json'
+    with open(fname_json_phase) as json_file:
+        json_data_phase = json.load(json_file)
+    nii_phase = nib.load(phases[0][0][0])
+    nii_profiles = nib.Nifti1Image(profiles, nii_phase.affine, header=nii_phase.header)
+
+    # Save nii and json
+    save_nii_json(nii_profiles, json_data_phase, fname_output)
+    logger.info(f"\n\n Filename of the coil profiles is: {fname_output}")
+
+    # Create coil config file
+    coil_name = json_data['name']
+
+    config_coil = {
+        'name': coil_name,
+        'coef_channel_minmax': json_data['coef_channel_minmax'],
+        'coef_sum_max': json_data['coef_sum_max']
+    }
+
+    # write json
+    fname_coil_config = os.path.join(path_output, coil_name + '_config.json')
+    with open(fname_coil_config, mode='w') as f:
+        json.dump(config_coil, f, indent=4)
+
+    logger.info(f"Filename of the coil config file is: {fname_coil_config}")
+
+
+def _concat_and_save_nii(list_fnames_nii, fname_output):
+    res = []
+    for _, fname in enumerate(list_fnames_nii):
+        nii = nib.load(fname)
+        nii.get_fdata()
+        res.append(nii.get_fdata())
+
+    fname_json = fname.split('.nii')[0] + '.json'
+    # Read from json file
+    with open(fname_json) as json_file:
+        json_data = json.load(json_file)
+
+    array_4d = np.moveaxis(np.array(res), 0, 3)
+    nii_4d = nib.Nifti1Image(array_4d, nii.affine, header=nii.header)
+    save_nii_json(nii_4d, json_data, fname_output)

shimmingtoolbox/cli/prepare_fieldmap.py

@@ -9,7 +9,7 @@
 
 from shimmingtoolbox.load_nifti import read_nii
 from shimmingtoolbox.prepare_fieldmap import prepare_fieldmap
-from shimmingtoolbox.utils import create_fname_from_path, set_all_loggers, create_output_dir
+from shimmingtoolbox.utils import create_fname_from_path, set_all_loggers, create_output_dir, save_nii_json
 
 logging.basicConfig(level=logging.INFO)
 logger = logging.getLogger(__name__)
@@ -56,13 +56,68 @@ def prepare_fieldmap_cli(phase, fname_mag, unwrapper, fname_output, autoscale, f
     # Set logger level
     set_all_loggers(verbose)
 
-    # Make sure output filename is valid
+    prepare_fieldmap_uncli(phase, fname_mag, unwrapper, fname_output, autoscale, fname_mask, threshold, fname_save_mask,
+                           gaussian_filter, sigma)
+
+
+def prepare_fieldmap_uncli(phase, fname_mag, unwrapper='prelude',
+                           fname_output=os.path.join(os.curdir, FILE_OUTPUT_DEFAULT), autoscale=True,
+                           fname_mask=None, threshold=0.05, fname_save_mask=None, gaussian_filter=False, sigma=1):
+    """ Prepare fieldmap cli without the click decorators. This allows this function to be imported and called from
+    other python modules.
+
+    Args:
+        phase (list): Input path of phase nifti file(s), in ascending order: echo1, echo2, etc.
+        fname_mag (str): Input path of mag nifti file
+        unwrapper (str): Algorithm for unwrapping. Supported unwrapper: 'prelude'.
+        fname_output (str): Output filename for the fieldmap, supported types : '.nii', '.nii.gz'
+        autoscale (bool): Tells whether to auto rescale phase inputs according to manufacturer standards. If you have
+                          non siemens data not automatically converted from dcm2niix, you should set this to False and
+                          input phase data from -pi to pi.
+        fname_mask (str): Input path for a mask. Mask must be the same shape as the array of each PHASE input.
+                          Used for PRELUDE
+        threshold (float): Threshold for masking. Used for: PRELUDE
+        fname_save_mask (str): Filename of the mask calculated by the unwrapper
+        gaussian_filter (bool): Gaussian filter for B0 map
+        sigma (float): Standard deviation of gaussian filter. Used for: gaussian_filter
+    """
+    # Return fieldmap and json file
+    nii_fieldmap, json_fieldmap = prepare_fieldmap_cli_inputs(phase, fname_mag, unwrapper, autoscale, fname_mask,
+                                                              threshold, fname_save_mask, gaussian_filter, sigma)
+
+    # Create filename for the output if it,s a path
     fname_output_v2 = create_fname_from_path(fname_output, FILE_OUTPUT_DEFAULT)
-    if fname_output_v2[-4:] != '.nii' and fname_output_v2[-7:] != '.nii.gz':
-        raise ValueError("Output filename must have one of the following extensions: '.nii', '.nii.gz'")
 
-    # Prepare the output
-    create_output_dir(fname_output_v2, is_file=True)
+    # Save fieldmap and json file to their respective filenames
+    save_nii_json(nii_fieldmap, json_fieldmap, fname_output_v2)
+
+    # Log output file
+    logger.info(f"Filename of the fieldmap is: {fname_output_v2}")
+
+
+def prepare_fieldmap_cli_inputs(phase, fname_mag, unwrapper, autoscale, fname_mask, threshold, fname_save_mask,
+                                gaussian_filter, sigma):
+    """Prepare fieldmap using click inputs
+
+    Args:
+        phase (list): Input path of phase nifti file(s), in ascending order: echo1, echo2, etc.
+        fname_mag (str): Input path of mag nifti file
+        unwrapper (str): Algorithm for unwrapping. Supported unwrapper: 'prelude'.
+        autoscale (bool): Tells whether to auto rescale phase inputs according to manufacturer standards. If you have
+                          non siemens data not automatically converted from dcm2niix, you should set this to False and
+                          input phase data from -pi to pi.
+        fname_mask (str): Input path for a mask. Mask must be the same shape as the array of each PHASE input.
+                          Used for PRELUDE
+        threshold (float): Threshold for masking. Used for: PRELUDE
+        gaussian_filter (bool): Gaussian filter for B0 map
+        sigma (float): Standard deviation of gaussian filter. Used for: gaussian_filter
+
+    Returns:
+        (tuple): tuple containing:
+
+            * nib.Nifti1Image: Nibabel object containing the fieldmap in hz.
+            * dict: Dictionary containing the json sidecar associated with the nibabel object fieldmap.
+    """
 
     # Save mask
     if fname_save_mask is not None:
@@ -106,21 +161,16 @@ def prepare_fieldmap_cli(phase, fname_mag, unwrapper, fname_output, autoscale, f
 
     # Save fieldmap
     nii_fieldmap = nib.Nifti1Image(fieldmap_hz, affine, header=nii_phase.header)
-    nib.save(nii_fieldmap, fname_output_v2)
 
     # Save fieldmap json
     json_fieldmap = json_phase
     if len(phase) > 1:
         for i_echo in range(len(echo_times)):
             json_fieldmap[f'EchoTime{i_echo + 1}'] = echo_times[i_echo]
-    fname_json = fname_output_v2.rsplit('.nii', 1)[0] + '.json'
-    with open(fname_json, 'w') as outfile:
-        json.dump(json_fieldmap, outfile, indent=2)
-
     # save mask json
     if fname_save_mask is not None:
         fname_mask_json = fname_save_mask.rsplit('.nii', 1)[0] + '.json'
         with open(fname_mask_json, 'w') as outfile:
             json.dump(json_mag, outfile, indent=2)
 
-    logger.info(f"Filename of the fieldmap is: {fname_output_v2}")
+    return nii_fieldmap, json_fieldmap