Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Split load_and_verify_mti function #933

Merged
merged 4 commits into from
Mar 7, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
117 changes: 82 additions & 35 deletions scilpy/io/mti.py
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,8 @@ def add_common_args_mti(p):
'smoothing, in number of voxels. [%(default)s]')


def verifications_and_loading_mti(args, parser, input_maps_lists,
extended_dir, affine, contrast_names):
def load_and_verify_mti(args, parser, input_maps_lists, extended_dir, affine,
contrast_names):
"""
Common verifications and loading for both MT and ihMT scripts.

Expand Down Expand Up @@ -125,44 +125,14 @@ def verifications_and_loading_mti(args, parser, input_maps_lists,

# Set TR and FlipAngle parameters. Required with --in_mtoff_t1, in which
# case one of --in_aqc_parameters or --in_jsons is set.
rep_times = None
flip_angles = None
if args.in_acq_parameters:
flip_angles = np.asarray(args.in_acq_parameters[:2]) * np.pi / 180.
rep_times = np.asarray(args.in_acq_parameters[2:]) * 1000
if rep_times[0] > 10000 or rep_times[1] > 10000:
logging.warning('Given repetition times do not seem to be given '
'in seconds. MTsat results might be affected.')
elif args.in_jsons:
rep_times = []
flip_angles = []
for curr_json in args.in_jsons:
acq_parameter = get_acq_parameters(curr_json,
['RepetitionTime', 'FlipAngle'])
if acq_parameter[0] > 10:
logging.warning('Repetition time found in {} does not seem to '
'be given in seconds. MTsat and ihMTsat '
'results might be affected.'.format(curr_json))
rep_times.append(acq_parameter[0] * 1000) # convert to ms.
flip_angles.append(np.deg2rad(acq_parameter[1]))
rep_times, flip_angles = _parse_acquisition_parameters(args)

# Fix issue from the presence of invalid value and division by zero
np.seterr(divide='ignore', invalid='ignore')

# Load B1 image
B1_map = None
if args.in_B1_map and args.in_mtoff_t1:
B1_img = nib.load(args.in_B1_map)
B1_map = B1_img.get_fdata(dtype=np.float32)
B1_map = adjust_B1_map_intensities(B1_map, nominal=args.B1_nominal)
B1_map = smooth_B1_map(B1_map, wdims=args.B1_smooth_dims)
if args.B1_correction_method == 'model_based':
# Apply the B1 map to the flip angles for model-based correction
flip_angles[0] *= B1_map
flip_angles[1] *= B1_map
if args.extended:
nib.save(nib.Nifti1Image(B1_map, affine),
os.path.join(extended_dir, "B1_map.nii.gz"))
B1_map, flip_angles = _prepare_B1_map(args, flip_angles, extended_dir,
affine)

# Define contrasts maps names
if args.filtering:
Expand Down Expand Up @@ -193,3 +163,80 @@ def verifications_and_loading_mti(args, parser, input_maps_lists,
contrast_names[idx] + '.nii.gz'))

return single_echo, flip_angles, rep_times, B1_map, contrast_maps


def _parse_acquisition_parameters(args):
"""
Parse the acquisition parameters from MTI, either from json files or
directly inputed parameters.

Parameters
----------
args: Namespace

Returns
-------
flip_angles: list[float]
The flip angles, in radian
rep_times: list[float]
The rep times, in ms.
"""
rep_times = None
flip_angles = None
if args.in_acq_parameters:
flip_angles = np.asarray(args.in_acq_parameters[:2]) * np.pi / 180.
rep_times = np.asarray(args.in_acq_parameters[2:]) * 1000
if rep_times[0] > 10000 or rep_times[1] > 10000:
logging.warning('Given repetition times do not seem to be given '
'in seconds. MTsat results might be affected.')
elif args.in_jsons:
rep_times = []
flip_angles = []
for curr_json in args.in_jsons:
acq_parameter = get_acq_parameters(curr_json,
['RepetitionTime', 'FlipAngle'])
if acq_parameter[0] > 10:
logging.warning('Repetition time found in {} does not seem to '
'be given in seconds. MTsat and ihMTsat '
'results might be affected.'.format(curr_json))
rep_times.append(acq_parameter[0] * 1000) # convert to ms.
flip_angles.append(np.deg2rad(acq_parameter[1]))
return rep_times, flip_angles


def _prepare_B1_map(args, flip_angles, extended_dir, affine):
"""
Prepare the B1 map for MTI B1+ inhomogeneity correction. Flip angles might
also be affected.

Parameters
----------
args: Namespace
flip_angles: list[float]
The flip angles, in radian
extended_dir: str
The folder for extended savings (with option args.extended).
affine: np.ndarray
A reference affine to save files.

Returns
-------
B1_map: np.ndarray
The loaded map, with adjusted intensities, smoothed, corrected.
flip_angles: list[float]
The modified flip angles, in radian
"""
B1_map = None
if args.in_B1_map and args.in_mtoff_t1:
B1_img = nib.load(args.in_B1_map)
B1_map = B1_img.get_fdata(dtype=np.float32)
B1_map = adjust_B1_map_intensities(B1_map, nominal=args.B1_nominal)
B1_map = smooth_B1_map(B1_map, wdims=args.B1_smooth_dims)
if args.B1_correction_method == 'model_based':
# Apply the B1 map to the flip angles for model-based correction
flip_angles[0] *= B1_map
flip_angles[1] *= B1_map
if args.extended:
nib.save(nib.Nifti1Image(B1_map, affine),
os.path.join(extended_dir, "B1_map.nii.gz"))
return B1_map, flip_angles
8 changes: 5 additions & 3 deletions scripts/scil_mti_maps_MT.py
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@
import nibabel as nib
import numpy as np

from scilpy.io.mti import add_common_args_mti, verifications_and_loading_mti
from scilpy.io.mti import add_common_args_mti, load_and_verify_mti
from scilpy.io.utils import (add_overwrite_arg,
assert_inputs_exist, add_verbose_arg,
assert_output_dirs_exist_and_empty)
Expand Down Expand Up @@ -138,7 +138,9 @@ def _build_arg_parser():
'calculation of MTsat. \nAcquisition '
'parameters should also be set with this image.')

# Other MTI arguments are gathered here.
add_common_args_mti(p)

add_verbose_arg(p)
add_overwrite_arg(p)

Expand Down Expand Up @@ -186,8 +188,8 @@ def main():

# Other checks, loading, saving contrast_maps.
single_echo, flip_angles, rep_times, B1_map, contrast_maps = \
verifications_and_loading_mti(args, parser, input_maps_lists,
extended_dir, affine, contrast_names)
load_and_verify_mti(args, parser, input_maps_lists, extended_dir,
affine, contrast_names)

# Compute MTR
if 'positive' in contrast_names_og and 'negative' in contrast_names_og:
Expand Down
8 changes: 5 additions & 3 deletions scripts/scil_mti_maps_ihMT.py
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,7 @@
import nibabel as nib
import numpy as np

from scilpy.io.mti import add_common_args_mti, verifications_and_loading_mti
from scilpy.io.mti import add_common_args_mti, load_and_verify_mti
from scilpy.io.utils import (add_overwrite_arg,
assert_inputs_exist, add_verbose_arg,
assert_output_dirs_exist_and_empty)
Expand Down Expand Up @@ -157,7 +157,9 @@ def _build_arg_parser():
'calculation of MTsat and ihMTsat. \nAcquisition '
'parameters should also be set with this image.')

# Other MTI arguments are gathered here.
add_common_args_mti(p)

add_verbose_arg(p)
add_overwrite_arg(p)

Expand Down Expand Up @@ -203,8 +205,8 @@ def main():

# Other checks, loading, saving contrast_maps.
single_echo, flip_angles, rep_times, B1_map, contrast_maps = \
verifications_and_loading_mti(args, parser, input_maps_lists,
extended_dir, affine, contrast_names)
load_and_verify_mti(args, parser, input_maps_lists, extended_dir,
affine, contrast_names)

# Compute ratio maps
MTR, ihMTR = compute_ratio_map((contrast_maps[2] + contrast_maps[3]) / 2,
Expand Down