diff --git a/config/coil_config.json b/config/coil_config.json index 3d75e088..89e20e76 100644 --- a/config/coil_config.json +++ b/config/coil_config.json @@ -1,15 +1,16 @@ { "name": "Prisma_fit", - "coef_channel_minmax": [ - [123100100, 123265000], - [-2300, 2300], - [-2300, 2300], - [-2300, 2300], - [-4959.01, 4959.01], - [-3551.29, 3551.29], - [-3503.299, 3503.299], - [-3551.29, 3551.29], - [-3487.302, 3487.302] - ], - "coef_sum_max": null + "coef_channel_minmax": + { + "0": [[123100100, 123265000]], + "1": [[-2300, 2300], + [-2300, 2300], + [-2300, 2300]], + "2": [[-4959.01, 4959.01], + [-3551.29, 3551.29], + [-3503.299, 3503.299], + [-3551.29, 3551.29], + [-3487.302, 3487.302]] + }, + "coef_sum_max": null } diff --git a/config/custom_coil_config.json b/config/custom_coil_config.json new file mode 100644 index 00000000..289dcc60 --- /dev/null +++ b/config/custom_coil_config.json @@ -0,0 +1,23 @@ +{ + "name": "custom", + "coef_channel_minmax": [ + [ + -2.5, + 2.5 + ], + [ + -2.5, + 2.5 + ], + [ + -2.5, + 2.5 + ], + [ + -2.5, + 2.5 + ] + ], + "coef_sum_max": 10, + "Units": "A" +} diff --git a/shimmingtoolbox/cli/b0shim.py b/shimmingtoolbox/cli/b0shim.py index 84ebab99..5eae1160 100644 --- a/shimmingtoolbox/cli/b0shim.py +++ b/shimmingtoolbox/cli/b0shim.py @@ -31,6 +31,7 @@ CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) logging.basicConfig(level=logging.INFO) logger = logging.getLogger(__name__) +AVAILABLE_ORDERS = [-1, 0, 1, 2] @click.group(context_settings=CONTEXT_SETTINGS, @@ -54,9 +55,11 @@ def b0shim_cli(): help="Anatomical image to apply the correction onto.") @click.option('--mask', 'fname_mask_anat', type=click.Path(exists=True), required=False, help="Mask defining the spatial region to shim.") -@click.option('--scanner-coil-order', type=click.Choice(['-1', '0', '1', '2']), default='-1', show_default=True, - help="Maximum order of the shim system. Note that specifying 1 will return " - "orders 0 and 1. The 0th order is the f0 frequency.") +@click.option('--scanner-coil-order', 'scanner_coil_order', type=click.STRING, default='-1', show_default=True, + help="Spherical harmonics orders to be used in optimization. " + f"Available orders: {AVAILABLE_ORDERS}. " + "Orders should be writen with a coma separating the values. (i.e. 0,1,2)" + "The 0th order is the f0 frequency.") @click.option('--scanner-coil-constraints', 'fname_sph_constr', type=click.Path(), default="", help=f"Constraints for the scanner coil. Example file located: {__dir_config_scanner_constraints__}") @click.option('--slices', type=click.Choice(['interleaved', 'sequential', 'volume', 'auto']), required=False, @@ -70,7 +73,7 @@ def b0shim_cli(): @click.option('--optimizer-method', 'method', type=click.Choice(['least_squares', 'pseudo_inverse', 'quad_prog']), required=False, default='quad_prog', show_default=True, help="Method used by the optimizer. LS, and QP will respect the constraints," - "PS will not respect the constraints") + "PS will not respect the constraints") @click.option('--regularization-factor', 'reg_factor', type=click.FLOAT, required=False, default=0.0, show_default=True, help="Regularization factor for the current when optimizing. A higher coefficient will penalize higher " "current values while 0 provides no regularization. Not relevant for 'pseudo-inverse' " @@ -136,12 +139,11 @@ def dynamic(fname_fmap, fname_anat, fname_mask_anat, method, opt_criteria, slice Example of use: st_b0shim dynamic --coil coil1.nii coil1_config.json --coil coil2.nii coil2_config.json --fmap fmap.nii --anat anat.nii --mask mask.nii --optimizer-method least_squares """ + + scanner_coil_order = parse_orders(scanner_coil_order) # Set logger level set_all_loggers(verbose) - # Input scanner_coil_order can be a string - scanner_coil_order = int(scanner_coil_order) - # Load the fieldmap nii_fmap_orig = nib.load(fname_fmap) @@ -236,7 +238,7 @@ def dynamic(fname_fmap, fname_anat, fname_mask_anat, method, opt_criteria, slice if output_value_format != 'delta': raise ValueError(f"Unsupported output value format: {output_value_format} for output file format: " f"{o_format_sph}") - if scanner_coil_order != 1: + if not (scanner_coil_order == [0, 1] or scanner_coil_order == [1]): raise ValueError(f"Unsupported scanner coil order: {scanner_coil_order} for output file format: " f"{o_format_sph}") if json_fm_data.get('Manufacturer') != 'Siemens': @@ -295,9 +297,10 @@ def dynamic(fname_fmap, fname_anat, fname_mask_anat, method, opt_criteria, slice if o_format_sph == 'gradient': logger.debug("Converting Siemens scanner coil from Shim CS (LAI) to Gradient CS") # First convert to RAS + orders = tuple([order for order in scanner_coil_order if order != 0]) for i_shim in range(coefs.shape[0]): # Convert coefficient - coefs_coil[i_shim, 1:] = shim_to_phys_cs(coefs_coil[i_shim, 1:], manufacturer) + coefs_coil[i_shim, 1:] = shim_to_phys_cs(coefs_coil[i_shim, 1:], manufacturer, orders) # Convert coef of 1st order sph harmonics to Gradient coord system coefs_freq, coefs_phase, coefs_slice = phys_to_gradient_cs(coefs_coil[:, 1], @@ -346,7 +349,8 @@ def dynamic(fname_fmap, fname_anat, fname_mask_anat, method, opt_criteria, slice # Select the coefficients for a coil coefs_coil = copy.deepcopy(coefs[:, start_channel:end_channel]) # Plot a figure of the coefficients - _plot_coefs(coil, list_slices, coefs_coil, path_output, i_coil, bounds=coil.coef_channel_minmax) + _plot_coefs(coil, list_slices, coefs_coil, path_output, i_coil, + bounds=[bound for bounds in coil.coef_channel_minmax.values() for bound in bounds]) logger.info(f"Finished plotting figure(s)") @@ -463,13 +467,23 @@ def _save_to_text_file_static(coil, coefs, list_slices, path_output, o_format, o @click.command(context_settings=CONTEXT_SETTINGS) -@click.option('--coil', 'coils', nargs=2, multiple=True, type=(click.Path(exists=True), click.Path(exists=True)), +@click.option('--coil', 'coils_static', nargs=2, multiple=True, type=(click.Path(exists=True), click.Path(exists=True)), help="Pair of filenames containing the coil profiles followed by the filename to the constraints " "e.g. --coil a.nii cons.json. If you have more than one coil, use this option more than once. " "The coil profiles and the fieldmaps (--fmap) must have matching units (if fmap is in Hz, the coil " "profiles must be in Hz/unit_shim). If you only want to shim using the scanner's gradient/shim " "coils, use the `--scanner-coil-order` option. For an example of a constraint file, " f"see: {__dir_config_scanner_constraints__}") +@click.option('--coil-riro', 'coils_riro', nargs=2, multiple=True, + type=(click.Path(exists=True), click.Path(exists=True)), required=False, + help="Pair of filenames containing the coil profiles followed by the filename to the constraints " + "e.g. --coil a.nii cons.json. If you have more than one coil, use this option more than once. " + "The coil profiles and the fieldmaps (--fmap) must have matching units (if fmap is in Hz, the coil " + "profiles must be in Hz/unit_shim). If this option is used, these coil profiles will be used for " + "the RIRO optimization, otherwise, the coils from the --coil options will be used." + "If you only want to shim using the scanner's gradient/shim " + "coils, use the `--scanner-coil-order` option. For an example of a constraint file, " + f"see: {__dir_config_scanner_constraints__}") @click.option('--fmap', 'fname_fmap', required=True, type=click.Path(exists=True), help="Timeseries of B0 fieldmap.") @click.option('--anat', 'fname_anat', type=click.Path(exists=True), required=True, @@ -481,9 +495,18 @@ def _save_to_text_file_static(coil, coefs, list_slices, path_output, o_format, o @click.option('--mask-riro', 'fname_mask_anat_riro', type=click.Path(exists=True), required=False, help="Mask defining the time varying (i.e. RIRO, Respiration-Induced Resonance Offset) " "region to shim.") -@click.option('--scanner-coil-order', type=click.Choice(['-1', '0', '1', '2']), default='-1', show_default=True, - help="Maximum order of the shim system. Note that specifying 1 will return " - "orders 0 and 1. The 0th order is the f0 frequency.") +@click.option('--scanner-coil-order', 'scanner_coil_order_static', type=click.STRING, default='-1', show_default=True, + help="Spherical harmonics orders to be used in static optimization. " + f"Available orders: {AVAILABLE_ORDERS}. " + "Orders should be writen with a coma separating the values. (i.e. 0,1,2)" + "The 0th order is the f0 frequency.") +@click.option('--scanner-coil-order-riro', 'scanner_coil_order_riro', type=click.STRING, default=None, + show_default=True, + help="Spherical harmonics orders to be used in RIRO optimization. If not set, the same orders as " + "--scanner-coil-order will be used for RIRO" + f"Available orders: {AVAILABLE_ORDERS}. " + "Orders should be writen with a coma separating the values. (i.e. 0,1,2)" + "The 0th order is the f0 frequency.") @click.option('--scanner-coil-constraints', 'fname_sph_constr', type=click.Path(), default="", help=f"Constraints for the scanner coil. Example file located: {__dir_config_scanner_constraints__}") @click.option('--slices', type=click.Choice(['interleaved', 'sequential', 'volume', 'auto']), required=False, @@ -550,9 +573,9 @@ def _save_to_text_file_static(coil, coefs, list_slices, path_output, o_format, o @click.option('-v', '--verbose', type=click.Choice(['info', 'debug']), default='info', help="Be more verbose") @timeit def realtime_dynamic(fname_fmap, fname_anat, fname_mask_anat_static, fname_mask_anat_riro, fname_resp, method, - opt_criteria, slices, slice_factor, coils, dilation_kernel_size, scanner_coil_order, - fname_sph_constr, fatsat, path_output, o_format_coil, o_format_sph, output_value_format, - reg_factor, verbose): + opt_criteria, slices, slice_factor, coils_static, coils_riro, dilation_kernel_size, + scanner_coil_order_static, scanner_coil_order_riro, fname_sph_constr, fatsat, path_output, + o_format_coil, o_format_sph, output_value_format, reg_factor, verbose): """ Realtime shim by fitting a fieldmap to a pressure monitoring unit. Use the option --optimizer-method to change the shimming algorithm used to optimize. Use the options --slices and --slice-factor to change the shimming order/size of the slices. @@ -560,9 +583,12 @@ def realtime_dynamic(fname_fmap, fname_anat, fname_mask_anat_static, fname_mask_ Example of use: st_b0shim realtime-dynamic --coil coil1.nii coil1_config.json --coil coil2.nii coil2_config.json --fmap fmap.nii --anat anat.nii --mask-static mask.nii --resp trace.resp --optimizer-method least_squares """ + # Set coils and scanner order for riro if none were indicated + if scanner_coil_order_riro is None: + scanner_coil_order_riro = scanner_coil_order_static - # Input can be a string - scanner_coil_order = int(scanner_coil_order) + scanner_coil_order_static = parse_orders(scanner_coil_order_static) + scanner_coil_order_riro = parse_orders(scanner_coil_order_riro) # Set logger level set_all_loggers(verbose) @@ -637,8 +663,9 @@ def realtime_dynamic(fname_fmap, fname_anat, fname_mask_anat_static, fname_mask_ if output_value_format == 'absolute': raise ValueError(f"Unsupported output value format: {output_value_format} for output file format: " f"{o_format_sph}") - if scanner_coil_order != 1: - raise ValueError(f"Unsupported scanner coil order: {scanner_coil_order} for output file format: " + if not (scanner_coil_order_static == [0, 1] or scanner_coil_order_static == [1]) or \ + not (scanner_coil_order_riro == [0, 1] or scanner_coil_order_riro == [1]): + raise ValueError(f"Unsupported scanner coil order: {scanner_coil_order_static} for output file format: " f"{o_format_sph}") if json_fm_data['Manufacturer'] != 'Siemens': raise ValueError(f"Unsupported manufacturer: {json_fm_data['manufacturer']} for output file format: " @@ -649,8 +676,12 @@ def realtime_dynamic(fname_fmap, fname_anat, fname_mask_anat_static, fname_mask_ options = {'scanner_shim': scanner_shim_settings.shim_settings} # Load the coils - list_coils = _load_coils(coils, scanner_coil_order, fname_sph_constr, nii_fmap, options['scanner_shim'], - json_fm_data['Manufacturer'], json_fm_data['ManufacturersModelName']) + list_coils_static = _load_coils(coils_static, scanner_coil_order_static, fname_sph_constr, nii_fmap, + options['scanner_shim'], json_fm_data['Manufacturer'], + json_fm_data['ManufacturersModelName']) + list_coils_riro = _load_coils(coils_riro, scanner_coil_order_riro, fname_sph_constr, nii_fmap, + options['scanner_shim'], json_fm_data['Manufacturer'], + json_fm_data['ManufacturersModelName']) if logger.level <= getattr(logging, 'DEBUG'): # Save inputs @@ -670,7 +701,7 @@ def realtime_dynamic(fname_fmap, fname_anat, fname_mask_anat_static, fname_mask_ # 1 ) Create the real time pmu sequencer object sequencer = RealTimeSequencer(nii_fmap_orig, json_fm_data, nii_anat, nii_mask_anat_static, nii_mask_anat_riro, - list_slices, pmu, list_coils, + list_slices, pmu, list_coils_static, list_coils_riro, method=method, opt_criteria=opt_criteria, mask_dilation_kernel='sphere', @@ -685,74 +716,164 @@ def realtime_dynamic(fname_fmap, fname_anat, fname_mask_anat_static, fname_mask_ # Load output options options['fatsat'] = _get_fatsat_option(json_anat_data, fatsat) - list_fname_output = [] - end_channel = 0 - for i_coil, coil in enumerate(list_coils): + # Get common coils between static and riro // Comparison based on coil name + coil_static_only = [coil for coil in list_coils_static if coil not in list_coils_riro] + coil_riro_only = [coil for coil in list_coils_riro if coil not in list_coils_static] + list_coils_common = [coil for coil in list_coils_static if coil in list_coils_riro] + # Create a list of all coils used in optimization + all_coils = list_coils_common + coil_static_only + coil_riro_only + + index = 0 + coil_indexes_static = {} + for coil in list_coils_static: + if type(coil) == Coil: + coil_indexes_static[coil.name] = [index, index + len(coil.coef_channel_minmax['coil'])] + index += len(coil.coef_channel_minmax['coil']) + else: + coil_indexes_static[coil.name] = {} + for key in coil.coef_channel_minmax: + coil_indexes_static[coil.name][key] = [index, index + len(coil.coef_channel_minmax[key])] + index += len(coil.coef_channel_minmax[key]) + + index = 0 + coil_indexes_riro = {} + for coil in list_coils_riro: + if type(coil) == Coil: + coil_indexes_riro[coil.name] = [index, index + len(coil.coef_channel_minmax['coil'])] + index += len(coil.coef_channel_minmax['coil']) + else: + coil_indexes_riro[coil.name] = {} + for key in coil.coef_channel_minmax: + coil_indexes_riro[coil.name][key] = [index, index + len(coil.coef_channel_minmax[key])] + index += len(coil.coef_channel_minmax[key]) + list_fname_output = [] + for i_coil, coil in enumerate(all_coils): # Figure out the start and end channels for a coil to be able to select it from the coefs - n_channels = coil.dim[3] - start_channel = end_channel - end_channel = start_channel + n_channels - - # Select the coefficients for a coil - coefs_coil_static = copy.deepcopy(coefs_static[:, start_channel:end_channel]) - coefs_coil_riro = copy.deepcopy(coefs_riro[:, start_channel:end_channel]) # If it's a scanner if type(coil) == ScannerCoil: - manufacturer = json_anat_data['Manufacturer'] - # If outputting in the gradient CS, it must be the 1st order and must be in the delta CS and Siemens - # The check has already been done earlier in the program to avoid processing and throw an error afterwards. - # Therefore, we can only check for the o_format_sph. - if o_format_sph == 'gradient': - logger.debug("Converting scanner coil from Shim CS to Gradient CS") + if coil in list_coils_common: + keys = [str(order) for order in AVAILABLE_ORDERS + if (order != -1 and (str(order) in coil_indexes_riro[coil.name] + or str(order) in coil_indexes_static[coil.name]))] + elif coil in coil_static_only: + keys = [str(order) for order in AVAILABLE_ORDERS + if (order != -1 and str(order) in coil_indexes_static[coil.name])] + elif coil in coil_riro_only: + keys = [str(order) for order in AVAILABLE_ORDERS + if (order != -1 and str(order) in coil_indexes_riro[coil.name])] + + for key in keys: + if coil in list_coils_riro: + if key in coil_indexes_riro[coil.name]: + coefs_coil_riro = copy.deepcopy( + coefs_riro[:, coil_indexes_riro[coil.name][key][0]:coil_indexes_riro[coil.name][key][1]]) + else: + coefs_coil_riro = np.zeros_like(coefs_static[:, coil_indexes_static[coil.name][key][0]: + coil_indexes_static[coil.name][key][1]]) + else: + coefs_coil_riro = np.zeros_like( + coefs_static[:, coil_indexes_static[coil.name][key][0]:coil_indexes_static[coil.name][key][1]]) - # First convert coefficients from Shim CS to RAS - for i_shim in range(coefs_coil_static.shape[0]): - # Convert coefficient - coefs_coil_static[i_shim, 1:] = shim_to_phys_cs(coefs_coil_static[i_shim, 1:], manufacturer) - coefs_coil_riro[i_shim, 1:] = shim_to_phys_cs(coefs_coil_riro[i_shim, 1:], manufacturer) - - # RAS to gradient - coefs_st_freq, coefs_st_phase, coefs_st_slice = phys_to_gradient_cs( - coefs_coil_static[:, 1], - coefs_coil_static[:, 2], - coefs_coil_static[:, 3], - fname_anat) - coefs_coil_static[:, 1] = coefs_st_freq - coefs_coil_static[:, 2] = coefs_st_phase - coefs_coil_static[:, 3] = coefs_st_slice - - coefs_riro_freq, coefs_riro_phase, coefs_riro_slice = phys_to_gradient_cs( - coefs_coil_riro[:, 1], - coefs_coil_riro[:, 2], - coefs_coil_riro[:, 3], - fname_anat) - coefs_coil_riro[:, 1] = coefs_riro_freq - coefs_coil_riro[:, 2] = coefs_riro_phase - coefs_coil_riro[:, 3] = coefs_riro_slice + if coil in list_coils_static: + if key in coil_indexes_static[coil.name]: + coefs_coil_static = copy.deepcopy(coefs_static[:, coil_indexes_static[coil.name][key][0]: + coil_indexes_static[coil.name][key][1]]) + else: + coefs_coil_static = np.zeros_like(coefs_coil_riro) + else: + coefs_coil_static = np.zeros_like(coefs_coil_riro) + + manufacturer = json_anat_data['Manufacturer'] + # If outputting in the gradient CS, it must be the 1st order and must be in the delta CS and Siemens + # The check has already been done earlier in the program to avoid processing and throw an error + # afterwards. + # Therefore, we can only check for the o_format_sph. + if o_format_sph == 'gradient': + if key == '0': + save_coefs_static = coefs_coil_static + if coefs_coil_riro is not None: + save_coefs_riro = coefs_coil_riro + else: + save_coefs_riro = None + has0 = True + continue + elif key == '1' and has0: + save_coefs_static = np.concatenate((save_coefs_static, coefs_coil_static), axis=1) + if save_coefs_riro is not None: + save_coefs_riro = np.concatenate((save_coefs_riro, coefs_coil_riro), axis=1) + elif coefs_coil_riro is not None: + save_coefs_riro = coefs_coil_riro + else: + raise ValueError("Orders do not match gradient") + coefs_coil_static = save_coefs_static + coefs_coil_riro = save_coefs_riro + logger.debug("Converting scanner coil from Shim CS to Gradient CS") + orders_static = tuple([order for order in scanner_coil_order_static if order != 0]) + orders_riro = tuple([order for order in scanner_coil_order_riro if order != 0]) + # First convert coefficients from Shim CS to RAS + for i_shim in range(coefs_coil_static.shape[0]): + # Convert coefficient + coefs_coil_static[i_shim, 1:] = shim_to_phys_cs(coefs_coil_static[i_shim, 1:], manufacturer, + orders_static) + coefs_coil_riro[i_shim, 1:] = shim_to_phys_cs(coefs_coil_riro[i_shim, 1:], manufacturer, + orders_riro) + + # RAS to gradient + coefs_st_freq, coefs_st_phase, coefs_st_slice = phys_to_gradient_cs( + coefs_coil_static[:, 1], + coefs_coil_static[:, 2], + coefs_coil_static[:, 3], + fname_anat) + coefs_coil_static[:, 1] = coefs_st_freq + coefs_coil_static[:, 2] = coefs_st_phase + coefs_coil_static[:, 3] = coefs_st_slice + + coefs_riro_freq, coefs_riro_phase, coefs_riro_slice = phys_to_gradient_cs( + coefs_coil_riro[:, 1], + coefs_coil_riro[:, 2], + coefs_coil_riro[:, 3], + fname_anat) + coefs_coil_riro[:, 1] = coefs_riro_freq + coefs_coil_riro[:, 2] = coefs_riro_phase + coefs_coil_riro[:, 3] = coefs_riro_slice - else: + else: - # If the output format is absolute, add the initial coefs - if output_value_format == 'absolute': - initial_coefs = scanner_shim_settings.concatenate_shim_settings(scanner_coil_order) - for i_channel in range(n_channels): - # abs_coef = delta + initial - coefs_coil_static[:, i_channel] = coefs_coil_static[:, i_channel] + initial_coefs[i_channel] - # riro does not change + # If the output format is absolute, add the initial coefs + if output_value_format == 'absolute' and coefs_coil_static is not None: + initial_coefs = scanner_shim_settings.concatenate_shim_settings(scanner_coil_order_static) + for i_channel in range(coefs_coil_static.shape[-1]): + # abs_coef = delta + initial + coefs_coil_static[:, i_channel] = coefs_coil_static[:, i_channel] + initial_coefs[i_channel] + # riro does not change - list_fname_output += _save_to_text_file_rt(coil, coefs_coil_static, coefs_coil_riro, mean_p, - list_slices, path_output, o_format_sph, options, i_coil, - default_st_coefs=initial_coefs) - continue + list_fname_output += _save_to_text_file_rt(coil, coefs_coil_static, coefs_coil_riro, mean_p, + list_slices, path_output, o_format_sph, options, + i_coil, int(key) ** 2, + default_st_coefs=initial_coefs) + continue - list_fname_output += _save_to_text_file_rt(coil, coefs_coil_static, coefs_coil_riro, mean_p, list_slices, - path_output, o_format_sph, options, i_coil) + list_fname_output += _save_to_text_file_rt(coil, coefs_coil_static, coefs_coil_riro, mean_p, + list_slices, + path_output, o_format_sph, options, i_coil, int(key) ** 2) else: # Custom coil + if coil in list_coils_riro: + coefs_coil_riro = copy.deepcopy( + coefs_riro[:, coil_indexes_riro[coil.name][0]:coil_indexes_riro[coil.name][1]]) + else: + coefs_coil_riro = np.zeros_like( + coefs_static[:, coil_indexes_static[coil.name][0]:coil_indexes_static[coil.name][1]]) + if coil in list_coils_static: + coefs_coil_static = copy.deepcopy( + coefs_static[:, coil_indexes_static[coil.name][0]:coil_indexes_static[coil.name][1]]) + else: + coefs_coil_static = np.zeros_like(coefs_coil_riro) + list_fname_output += _save_to_text_file_rt(coil, coefs_coil_static, coefs_coil_riro, mean_p, list_slices, - path_output, o_format_coil, options, i_coil) + path_output, o_format_coil, options, i_coil, 0) logger.info(f"Coil txt file(s) are here:\n{os.linesep.join(list_fname_output)}") logger.info(f"Plotting figure(s)") @@ -760,36 +881,43 @@ def realtime_dynamic(fname_fmap, fname_anat, fname_mask_anat_static, fname_mask_ logger.info(f"Plotting Currents") # Plot the coefs after outputting the currents to the text file end_channel = 0 - for i_coil, coil in enumerate(list_coils): - # Figure out the start and end channels for a coil to be able to select it from the coefs - n_channels = coil.dim[3] - start_channel = end_channel - end_channel = start_channel + n_channels + for i_coil, coil in enumerate(all_coils): + # Figure out the start and end channels for a coil to be able to select it from the coefs if type(coil) != ScannerCoil: - # Select the coefficients for a coil - coefs_coil_static = copy.deepcopy(coefs_static[:, start_channel:end_channel]) - coefs_coil_riro = copy.deepcopy(coefs_riro[:, start_channel:end_channel]) + if coil in list_coils_riro: + coefs_coil_riro = copy.deepcopy( + coefs_riro[:, coil_indexes_riro[coil.name][0]:coil_indexes_riro[coil.name][1]]) + else: + coefs_coil_riro = None + if coil in list_coils_static: + coefs_coil_static = copy.deepcopy( + coefs_static[:, coil_indexes_static[coil.name][0]:coil_indexes_static[coil.name][1]]) + else: + coefs_coil_static = np.zeros_like(coefs_coil_riro) # Plot a figure of the coefficients _plot_coefs(coil, list_slices, coefs_coil_static, path_output, i_coil, coefs_coil_riro, pres_probe_max=pmu.max - mean_p, pres_probe_min=pmu.min - mean_p, - bounds=coil.coef_channel_minmax) + bounds=[bound for bounds in coil.coef_channel_minmax.values() for bound in bounds]) logger.info(f"Finished plotting figure(s)") def _save_to_text_file_rt(coil, currents_static, currents_riro, mean_p, list_slices, path_output, o_format, - options, coil_number, default_st_coefs=None): + options, coil_number, channel_start, default_st_coefs=None): """o_format can either be 'chronological-ch', 'chronological-coil', 'gradient'""" list_fname_output = [] - n_channels = coil.dim[3] - + if currents_riro is not None: + n_channels = currents_riro.shape[-1] + else: + n_channels = currents_static.shape[-1] # Write a file for each channel for i_channel in range(n_channels): if o_format == 'chronological-ch': - fname_output = os.path.join(path_output, f"coefs_coil{coil_number}_ch{i_channel}_{coil.name}.txt") + fname_output = os.path.join(path_output, + f"coefs_coil{coil_number}_ch{channel_start + i_channel}_{coil.name}.txt") with open(fname_output, 'w', encoding='utf-8') as f: # Each row will have 3 coef representing the static, riro and mean_p in chronological order for i_shim in range(len(list_slices)): @@ -801,19 +929,25 @@ def _save_to_text_file_rt(coil, currents_static, currents_riro, mean_p, list_sli else: # Output initial coefs (absolute) f.write(f"{default_st_coefs[i_channel]:.1f}, {0:.1f}, {mean_p:.4f},\n") - f.write(f"{currents_static[i_shim, i_channel]:.6f}, ") - f.write(f"{currents_riro[i_shim, i_channel]:.12f}, ") + if currents_static is not None: + f.write(f"{currents_static[i_shim, i_channel]:.6f}, ") + if currents_riro is not None: + f.write(f"{currents_riro[i_shim, i_channel]:.12f}, ") f.write(f"{mean_p:.4f},\n") elif o_format == 'slicewise-ch': - fname_output = os.path.join(path_output, f"coefs_coil{coil_number}_ch{i_channel}_{coil.name}.txt") + fname_output = os.path.join(path_output, + f"coefs_coil{coil_number}_ch{channel_start + i_channel}_{coil.name}.txt") with open(fname_output, 'w', encoding='utf-8') as f: # Each row will have one coef representing the static, riro and mean_p in slicewise order n_slices = np.sum([len(a_tuple) for a_tuple in list_slices]) for i_slice in range(n_slices): i_shim = [list_slices.index(i) for i in list_slices if i_slice in i][0] - f.write(f"{currents_static[i_shim, i_channel]:.6f}, ") - f.write(f"{currents_riro[i_shim, i_channel]:.12f}, ") + + if currents_static is not None: + f.write(f"{currents_static[i_shim, i_channel]:.6f}, ") + if currents_riro is not None: + f.write(f"{currents_riro[i_shim, i_channel]:.12f}, ") f.write(f"{mean_p:.4f},\n") else: # o_format == 'gradient': @@ -835,18 +969,22 @@ def _save_to_text_file_rt(coil, currents_static, currents_riro, mean_p, list_sli if i_channel == 0: # f0, Output is in Hz - f.write(f"corr_vec[0][{i_slice}]= " - f"{currents_static[i_shim, i_channel]:.6f}\n") - f.write(f"corr_vec[1][{i_slice}]= " - f"{currents_riro[i_shim, i_channel]:.12f}\n") + if currents_static is not None: + f.write(f"corr_vec[0][{i_slice}]= " + f"{currents_static[i_shim, i_channel]:.6f}\n") + if currents_riro is not None: + f.write(f"corr_vec[1][{i_slice}]= " + f"{currents_riro[i_shim, i_channel]:.12f}\n") f.write(f"corr_vec[2][{i_slice}]= {mean_p:.3f}\n") else: # For Gx, Gy, Gz: Divide by 1000 for mT/m - f.write(f"corr_vec[0][{i_slice}]= " - f"{currents_static[i_shim, i_channel] / 1000:.6f}\n") - f.write(f"corr_vec[1][{i_slice}]= " - f"{currents_riro[i_shim, i_channel] / 1000:.12f}\n") + if currents_static is not None: + f.write(f"corr_vec[0][{i_slice}]= " + f"{currents_static[i_shim, i_channel] / 1000:.6f}\n") + if currents_riro is not None: + f.write(f"corr_vec[1][{i_slice}]= " + f"{currents_riro[i_shim, i_channel] / 1000:.12f}\n") f.write(f"corr_vec[2][{i_slice}]= {mean_p:.3f}\n") list_fname_output.append(os.path.abspath(fname_output)) @@ -854,16 +992,29 @@ def _save_to_text_file_rt(coil, currents_static, currents_riro, mean_p, list_sli return list_fname_output -def _load_coils(coils, order, fname_constraints, nii_fmap, scanner_shim_settings, manufacturer, +def parse_orders(orders: str): + orders = orders.split(',') + try: + orders = [int(order) for order in orders] + orders.sort() + if any(order not in AVAILABLE_ORDERS for order in orders): + raise ValueError(f'Orders must be selected from: {AVAILABLE_ORDERS}') + return orders + except ValueError: + raise ValueError(f"Invalid orders: {orders}\n Orders must be integers ") + + +def _load_coils(coils, orders, fname_constraints, nii_fmap, scanner_shim_settings, manufacturer, manufacturers_model_name): + # ! Modify description if everything works """ Loads the Coil objects from filenames Args: coils (list): List of tuples(fname_nii, fname_json) of coil profiles and constraints - order (int): Order of the scanner coils (0 or 1 or 2) + orders (list): Orders of the scanner coils (0 or 1 or 2) fname_constraints (str): Filename of the constraints of the scanner coils nii_fmap (nib.Nifti1Image): Nibabel object of the fieldmap - scanner_shim_settings (dict): Dictionary containing the shim settings of the scanner ('f0', 'order1', 'order2') + scanner_shim_settings (dict): Dictionary containing the shim settings of the scanner ('0', '1', '2') manufacturer (str): Name of the MRI manufacturer manufacturers_model_name (str): Name of the scanner @@ -879,19 +1030,26 @@ def _load_coils(coils, order, fname_constraints, nii_fmap, scanner_shim_settings constraints = json.load(json_file) list_coils.append(Coil(nii_coil_profiles.get_fdata(), nii_coil_profiles.affine, constraints)) + if len(list_coils) != len(set(list_coils)): + raise ValueError("Coils must be unique. Make sure different coils have different names.") + # Create the spherical harmonic coil profiles of the scanner - if 0 <= order <= 2: + if -1 not in orders: if os.path.isfile(fname_constraints): with open(fname_constraints) as json_file: sph_contraints = json.load(json_file) + orders_to_delete = [] + for key in sph_contraints['coef_channel_minmax']: + if key not in str(orders): + orders_to_delete.append(key) + for key in orders_to_delete: + del sph_contraints['coef_channel_minmax'][key] else: - sph_contraints = get_scanner_constraints(manufacturers_model_name) - - sph_contraints_calc = calculate_scanner_constraints(sph_contraints, scanner_shim_settings, order, manufacturer) + sph_contraints = get_scanner_constraints(manufacturers_model_name, orders) - # Create a ScannerCoil object - scanner_coil = ScannerCoil(nii_fmap.shape[:3], nii_fmap.affine, sph_contraints_calc, order, + sph_contraints_calc = calculate_scanner_constraints(sph_contraints, scanner_shim_settings, orders) + scanner_coil = ScannerCoil(nii_fmap.shape[:3], nii_fmap.affine, sph_contraints_calc, orders, manufacturer=manufacturer) list_coils.append(scanner_coil) @@ -902,70 +1060,68 @@ def _load_coils(coils, order, fname_constraints, nii_fmap, scanner_shim_settings return list_coils -def calculate_scanner_constraints(constraints, scanner_shim_settings, order, manufacturer): +def calculate_scanner_constraints(constraints: dict, scanner_shim_settings, orders): + # ! Modify description if everything works """ Calculate the constraints that should be used for the scanner by considering the current shim settings and the absolute bounds Args: constraints (dict): Constraints of the scanner coils - scanner_shim_settings (dict): Dictionary containing the shim settings of the scanner ('f0', 'order1', 'order2') - order (int): Order of the scanner coils (0 or 1 or 2) + scanner_shim_settings (dict): Dictionary containing the shim settings of the scanner ('0', '1', '2') + orders (list): Order of the scanner coils (0 or 1 or 2) manufacturer (str): Name of the MRI manufacturer Returns: dict: Updated constraints of the scanner """ - def _initial_in_bounds(coefs, bounds): + + def _initial_in_bounds(coefs: dict, bounds: dict): """Makes sure the initial values are within the bounds of the constraints""" - if len(coefs) != len(bounds): + if coefs.keys() != bounds.keys(): + raise RuntimeError("The scanner coil's orders is not the same length as the initial orders") + if any(len(coefs[key]) != len(bounds[key]) for key in bounds): raise RuntimeError("The scanner coil's bounds is not the same length as the initial bounds") - for i_bound in range(len(bounds)): - if bounds[i_bound][0] is None and bounds[i_bound][1] is None: - continue - if bounds[i_bound][1] is not None: - if not coefs[i_bound] <= bounds[i_bound][1]: - logger.warning(f"Initial scanner coefs are outside the bounds allowed in the constraints: " - f"{bounds[i_bound]}, initial: {coefs[i_bound]}") - if bounds[i_bound][0] is not None: - if not bounds[i_bound][0] <= coefs[i_bound]: - logger.warning(f"Initial scanner coefs are outside the bounds allowed in the constraints: " - f"{bounds[i_bound]}, initial: {coefs[i_bound]}") + + for key in coefs: + for (bound, coef) in zip(bounds[key], coefs[key]): + if bound[0] is None and bound[1] is None: + continue + if bound[1] is not None: + if not coef <= bound[1]: + logger.warning(f"Initial scanner coefs are outside the bounds allowed in the constraints: " + f"{bound}, initial: {coef}") + if bound[0] is not None: + if not bound[0] <= coef: + logger.warning(f"Initial scanner coefs are outside the bounds allowed in the constraints: " + f"{bound}, initial: {coef}") # Set the initial coefficients to 0 - if order == 0: - # Order 0 has 1 coefficient (f0) - initial_coefs = [0] - elif order == 1: - # Order 1 has 3 more coefficients than order 0 (f0, x, y, z) - initial_coefs = [0] * 4 - elif order == 2: - # Order 2 has 5 more coefficients than order 1 (f0, X, Y, Z, Z2, ZX, ZY, X2-Y2, XY) - initial_coefs = [0] * 9 - else: + initial_coefs = {} + for order in orders: + initial_coefs[str(order)] = [0] * (order * 2 + 1) + if initial_coefs == {}: initial_coefs = None # Restrict the constraints to the provided order - constraints['coef_channel_minmax'] = restrict_sph_constraints(constraints['coef_channel_minmax'], order) + constraints['coef_channel_minmax'] = restrict_sph_constraints(constraints['coef_channel_minmax'], orders) # If the scanner coefficients are valid, update the initial coefficients if scanner_shim_settings['has_valid_settings']: - - if scanner_shim_settings['f0'] is not None and order >= 0: - initial_coefs[0] = scanner_shim_settings['f0'] - if scanner_shim_settings['order1'] is not None and order >= 1: - initial_coefs[1:4] = scanner_shim_settings['order1'] - if scanner_shim_settings['order2'] is not None and order >= 2: - initial_coefs[4:9] = scanner_shim_settings['order2'] + if scanner_shim_settings['0'] is not None and 0 in orders: + initial_coefs['0'] = np.array([scanner_shim_settings['0']]) + if scanner_shim_settings['1'] is not None and 1 in orders: + initial_coefs['1'] = scanner_shim_settings['1'] + if scanner_shim_settings['2'] is not None and 2 in orders: + initial_coefs['2'] = scanner_shim_settings['2'] # Make sure the initial coefficients are within the specified bounds _initial_in_bounds(initial_coefs, constraints['coef_channel_minmax']) # Update the bounds to what they should be by taking into account that the fieldmap was acquired using some # shimming - constraints['coef_channel_minmax'] = new_bounds_from_currents(np.array([initial_coefs]), + constraints['coef_channel_minmax'] = new_bounds_from_currents(initial_coefs, constraints['coef_channel_minmax'] - )[0] - + ) return constraints diff --git a/shimmingtoolbox/cli/create_coil_profiles.py b/shimmingtoolbox/cli/create_coil_profiles.py index 7f4eb4ab..d5c858c4 100644 --- a/shimmingtoolbox/cli/create_coil_profiles.py +++ b/shimmingtoolbox/cli/create_coil_profiles.py @@ -375,7 +375,7 @@ def from_cad(fname_txt, fname_fmap, offset, dims_to_flip, software, coil_name, m # Create the coil profiles json file if max_current_sum is None: max_current_sum = nb_channels - coef_channel_minmax = [[min_current, max_current]] * nb_channels + coef_channel_minmax = {"coil": [[min_current, max_current]] * nb_channels} config_coil = { 'name': coil_name, 'coef_channel_minmax': coef_channel_minmax, diff --git a/shimmingtoolbox/cli/download_data.py b/shimmingtoolbox/cli/download_data.py index d9ff0ac4..58e08a90 100755 --- a/shimmingtoolbox/cli/download_data.py +++ b/shimmingtoolbox/cli/download_data.py @@ -15,7 +15,7 @@ URL_DICT: Dict[str, Tuple[List[str], str]] = { "testing_data": ( - ["https://github.com/shimming-toolbox/data-testing/archive/r20230809.zip"], + ["https://github.com/shimming-toolbox/data-testing/archive/r20230929.zip"], "Light-weighted dataset for testing purpose.", ), "prelude": ( @@ -27,8 +27,8 @@ "B0 dynamic shimming dataset.", ), "data_create_coil_profiles": ( - ["https://github.com/shimming-toolbox/data-create-coil-profiles/releases/download/r20230815/dicoms.zip"], - ["https://github.com/shimming-toolbox/data-create-coil-profiles/archive/refs/tags/r20230815.zip"], + ["https://github.com/shimming-toolbox/data-create-coil-profiles/releases/download/r20230929/dicoms.zip"], + ["https://github.com/shimming-toolbox/data-create-coil-profiles/archive/refs/tags/r20230929.zip"], "B0 coil profile dataset.", ), "data_b1_shimming": ( diff --git a/shimmingtoolbox/coils/coil.py b/shimmingtoolbox/coils/coil.py index 39f6c693..522c22a2 100644 --- a/shimmingtoolbox/coils/coil.py +++ b/shimmingtoolbox/coils/coil.py @@ -1,5 +1,6 @@ #!/usr/bin/python3 # -*- coding: utf-8 -*- +import copy import logging import numpy as np @@ -30,7 +31,7 @@ class Coil(object): profile. This transformation relates to the physical coordinates of the scanner (qform). required_constraints (list): List containing the required keys for ``constraints`` coef_sum_max (float): Contains the maximum value for the sum of the coefficients - coef_channel_minmax (list): List of ``(min, max)`` pairs for each coil channels. (None, None) is + coef_channel_minmax (dict): Dict of ``(min, max)`` pairs for each coil channels. (None, None) is used to specify no bounds. name (str): Name of the coil. """ @@ -85,27 +86,27 @@ def profile(self, profile): self.dim = profile.shape self._profile = profile - def load_constraints(self, constraints): + def load_constraints(self, constraints: dict): """Loads the constraints named in required_constraints as attribute to this class""" - # global `required_constraints` for key_name in required_constraints: if key_name in constraints: - if key_name == "coef_channel_minmax": - if len(constraints["coef_channel_minmax"]) != self.dim[3]: + if sum([len(constraints["coef_channel_minmax"][key]) for key in + constraints["coef_channel_minmax"]]) != self.dim[3]: raise ValueError(f"length of 'coef_channel_max' must be the same as the number of channels: " - f"{self.dim[3]}") - - for i_channel in range(self.dim[3]): - if constraints["coef_channel_minmax"][i_channel] is None: - constraints["coef_channel_minmax"][i_channel] = (-np.inf, np.inf) - if constraints["coef_channel_minmax"][i_channel][0] is None: - constraints["coef_channel_minmax"][i_channel] = \ - (-np.inf, constraints["coef_channel_minmax"][i_channel][1]) - if constraints["coef_channel_minmax"][i_channel][1] is None: - constraints["coef_channel_minmax"][i_channel] = \ - (constraints["coef_channel_minmax"][i_channel][0], np.inf) + f"{self.dim[3]} {sum([len(constraints['coef_channel_minmax'][key]) for key in constraints['coef_channel_minmax']])}") + + for key in constraints["coef_channel_minmax"]: + for i in range(len(constraints["coef_channel_minmax"][key])): + if constraints["coef_channel_minmax"][key][i] is None: + constraints["coef_channel_minmax"][key][i] = (-np.inf, np.inf) + if constraints["coef_channel_minmax"][key][i][0] is None: + constraints["coef_channel_minmax"][key][i] = \ + (-np.inf, constraints["coef_channel_minmax"][key][i][1]) + if constraints["coef_channel_minmax"][key][i][1] is None: + constraints["coef_channel_minmax"][key][i] = \ + (constraints["coef_channel_minmax"][key][i][0], np.inf) if key_name == "coef_sum_max": if constraints["coef_sum_max"] is None: @@ -115,12 +116,19 @@ def load_constraints(self, constraints): else: raise KeyError(f"Missing required constraint: {key_name}") + def __hash__(self): + return hash(self.name) + + def __eq__(self, __value: object) -> bool: + return self.name == __value.name + class ScannerCoil(Coil): """Coil class for scanner coils as they require extra arguments""" - def __init__(self, dim_volume, affine, constraints, order, manufacturer=""): - self.order = order + def __init__(self, dim_volume, affine, constraints, orders, manufacturer=""): + + self.orders = orders manufacturer = manufacturer.upper() if manufacturer in shim_cs: @@ -134,45 +142,57 @@ def __init__(self, dim_volume, affine, constraints, order, manufacturer=""): # Create the spherical harmonics with the correct order, dim and affine sph_coil_profile = self._create_coil_profile(dim_volume, manufacturer) # Restricts the constraints to the specified order - constraints['coef_channel_minmax'] = restrict_sph_constraints(constraints['coef_channel_minmax'], self.order) + constraints['coef_channel_minmax'] = restrict_sph_constraints(constraints['coef_channel_minmax'], self.orders) super().__init__(sph_coil_profile, affine, constraints) def _create_coil_profile(self, dim, manufacturer=None): # Define profile for Tx (constant volume) - profile_order_0 = -np.ones(dim) - - # Create spherical harmonics coil profiles - if self.order == 0: + if 0 in self.orders: + profile_order_0 = -np.ones(dim) + else: + profile_order_0 = None + # define the coil profiles + if self.orders == [0]: # f0 --> [1] sph_coil_profile = profile_order_0[..., np.newaxis] else: # f0, orders + temp_orders = [order for order in self.orders if order != 0] mesh1, mesh2, mesh3 = generate_meshgrid(dim, self.affine) + if manufacturer == 'SIEMENS': - profile_orders = siemens_basis(mesh1, mesh2, mesh3, orders=tuple(range(1, self.order + 1)), + profile_orders = siemens_basis(mesh1, mesh2, mesh3, orders=tuple(temp_orders), shim_cs=self.coord_system) elif manufacturer == 'GE': - profile_orders = ge_basis(mesh1, mesh2, mesh3, orders=tuple(range(1, self.order + 1)), + profile_orders = ge_basis(mesh1, mesh2, mesh3, orders=tuple(temp_orders), shim_cs=self.coord_system) else: logger.warning(f"{manufacturer} manufacturer not implemented. Outputting in Hz, uT/m, uT/m^2 for order " f"0, 1 and 2 respectively") - profile_orders = siemens_basis(mesh1, mesh2, mesh3, orders=tuple(range(1, self.order + 1)), + profile_orders = siemens_basis(mesh1, mesh2, mesh3, orders=tuple(temp_orders), shim_cs=self.coord_system) - - sph_coil_profile = np.concatenate((profile_order_0[..., np.newaxis], profile_orders), axis=3) + if profile_order_0 is None: + sph_coil_profile = profile_orders + else: + sph_coil_profile = np.concatenate((profile_order_0[..., np.newaxis], profile_orders), axis=3) return sph_coil_profile + def __hash__(self): + return hash(self.name) + + def __eq__(self, __value: object) -> bool: + return super().__eq__(__value) -def get_scanner_constraints(manufacturers_model_name, order=2): + +def get_scanner_constraints(manufacturers_model_name, orders): """ Returns the scanner spherical harmonics constraints depending on the manufacturer's model name and required order Args: manufacturers_model_name (str): Name of the scanner - order (int): Maximum order of the shim system + orders (list): List of all orders of the shim system to be used Returns: dict: The constraints including the scanner name, bounds and the maximum sum of currents. @@ -181,80 +201,82 @@ def get_scanner_constraints(manufacturers_model_name, order=2): if manufacturers_model_name == "Prisma_fit": constraints = { "name": "Prisma_fit", - "coef_channel_minmax": [], + "coef_channel_minmax": {"0": [], "1": [], "2": []}, "coef_sum_max": None } - if order >= 0: - constraints["coef_channel_minmax"].append([123100100, 123265000]) - if order >= 1: + if 0 in orders: + constraints["coef_channel_minmax"]["0"].append([123100100, 123265000]) + if 1 in orders: for _ in range(3): - constraints["coef_channel_minmax"].append([-2300, 2300]) - if order >= 2: - constraints["coef_channel_minmax"].extend([[-4959.01, 4959.01], - [-3551.29, 3551.29], - [-3503.299, 3503.299], - [-3551.29, 3551.29], - [-3487.302, 3487.302]]) + constraints["coef_channel_minmax"]["1"].append([-2300, 2300]) + if 2 in orders: + constraints["coef_channel_minmax"]["2"].extend([[-4959.01, 4959.01], + [-3551.29, 3551.29], + [-3503.299, 3503.299], + [-3551.29, 3551.29], + [-3487.302, 3487.302]]) elif manufacturers_model_name == "Investigational_Device_7T": constraints = { "name": "Investigational_Device_7T", - "coef_channel_minmax": [], + "coef_channel_minmax": {"0": [], "1": [], "2": []}, "coef_sum_max": None } - if order >= 0: + if 0 in orders: pass # todo: f0 min and max is wrong - constraints["coef_channel_minmax"].append([None, None]) - if order >= 1: + constraints["coef_channel_minmax"]["0"].append([None, None]) + if 1 in orders: for _ in range(3): - constraints["coef_channel_minmax"].append([-5000, 5000]) - if order >= 2: - constraints["coef_channel_minmax"].extend([[-1839.63, 1839.63], - [-791.84, 791.84], - [-791.84, 791.84], - [-615.87, 615.87], - [-615.87, 615.87]]) + constraints["coef_channel_minmax"]["1"].append([-5000, 5000]) + if 2 in orders: + constraints["coef_channel_minmax"]["2"].extend([[-1839.63, 1839.63], + [-791.84, 791.84], + [-791.84, 791.84], + [-615.87, 615.87], + [-615.87, 615.87]]) else: logger.warning(f"Scanner: {manufacturers_model_name} constraints not yet implemented, constraints might not be " "respected.") constraints = { "name": "Unknown", + "coef_channel_minmax": {"0": [], "1": [], "2": []}, "coef_sum_max": None } - if order == 0: - constraints["coef_channel_minmax"] = [[None, None]] - elif order == 1: - constraints["coef_channel_minmax"] = [[None, None] for _ in range(4)] - elif order == 2: - constraints["coef_channel_minmax"] = [[None, None] for _ in range(9)] + if 0 in orders: + constraints["coef_channel_minmax"]["0"] = [[None, None]] + if 1 in orders: + constraints["coef_channel_minmax"]["1"] = [[None, None] for _ in range(3)] + if 2 in orders: + constraints["coef_channel_minmax"]["2"] = [[None, None] for _ in range(5)] return constraints -def restrict_sph_constraints(bounds, order): +def restrict_sph_constraints(bounds: dict, orders): + # ! Modify description if everything works """ Select bounds according to the order specified Args: - bounds (list): 2D list (n_channels, 2) containing the min and max currents for multiple spherical harmonics + bounds (dict): Dictionary containing the min and max currents for multiple spherical harmonics orders - order (int): Maximum order of spherical harmonics + orders (list): Lsit of all spherical harmonics orders to be used Returns: - list: 2D list with the bounds of order 0 to the specified order + dict: Dictionary with the bounds of all specified orders """ - - if order == 0: + minmax_out = {} + if 0 in orders: # f0 --> [1] - minmax_out = bounds[:1] - elif order == 1: + minmax_out["0"] = bounds["0"] + if 1 in orders: # f0, ch1, ch2, ch3 -- > [4] - minmax_out = bounds[:4] - elif order == 2: + minmax_out["1"] = bounds["1"] + if 2 in orders: # f0, ch1, ch2, ch3, ch4, ch5, ch6, ch7, ch8 -- > [9] - minmax_out = bounds[:9] - else: + minmax_out["2"] = bounds["2"] + if minmax_out == {}: raise NotImplementedError("Order must be between 0 and 2") return minmax_out diff --git a/shimmingtoolbox/coils/spher_harm_basis.py b/shimmingtoolbox/coils/spher_harm_basis.py index 4f277ef8..e35c128f 100644 --- a/shimmingtoolbox/coils/spher_harm_basis.py +++ b/shimmingtoolbox/coils/spher_harm_basis.py @@ -9,7 +9,7 @@ logger = logging.getLogger(__name__) GYROMAGNETIC_RATIO = 42.5774785178325552 # [MHz/T] - +ORDER_INDEXES = {1: 3, 2: 5} def siemens_basis(x, y, z, orders=(1, 2), shim_cs='LAI'): """ @@ -52,14 +52,20 @@ def siemens_basis(x, y, z, orders=(1, 2), shim_cs='LAI'): _check_basis_inputs(x, y, z, orders) # Create spherical harmonics from first to second order - all_orders = np.array(range(1, 3)) + all_orders = np.array([order for order in orders if order != 0]) spher_harm = scaled_spher_harm(x, y, z, all_orders, shim_cs=shim_cs) # Reorder according to siemens convention: X, Y, Z, Z2, ZX, ZY, X2-Y2, XY - reordered_spher = _reorder_to_siemens(spher_harm) + reordered_spher = _reorder_to_siemens(spher_harm, orders) # Select order - range_per_order = {1: list(range(3)), 2: list(range(3, 8))} + range_per_order = {} + index = 0 + for order in orders: + range_per_order[order] = list(range(index, index + ORDER_INDEXES[order])) + index += ORDER_INDEXES[order] + + # range_per_order = {1: list(range(3)), 2: list(range(3, 8))} length_dim3 = np.sum([len(values) for key, values in range_per_order.items() if key in orders]) output = np.zeros(reordered_spher[..., 0].shape + (length_dim3,), dtype=reordered_spher.dtype) start_index = 0 @@ -72,7 +78,7 @@ def siemens_basis(x, y, z, orders=(1, 2), shim_cs='LAI'): return output -def _reorder_to_siemens(spher_harm): +def _reorder_to_siemens(spher_harm, orders): """ Reorder 1st - 2nd order coefficients along the last dim. From 1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics), to @@ -81,15 +87,28 @@ def _reorder_to_siemens(spher_harm): Args: spher_harm (numpy.ndarray): Coefficients with the last dimensions representing the different order channels. ``spher_harm.shape[-1]`` must equal 8. + orders (tuple): Spherical harmonics orders to use Returns: numpy.ndarray: Coefficients ordered following Siemens convention """ + if orders == (1, 2): + if spher_harm.shape[-1] != 8: + raise RuntimeError("Input arrays should have 4th dimension's shape equal to 8") + + reordered = spher_harm[..., [2, 0, 1, 5, 6, 4, 7, 3]] - if spher_harm.shape[-1] != 8: - raise RuntimeError("Input arrays should have 4th dimension's shape equal to 8") + if orders == (2,): + if spher_harm.shape[-1] != 5: + raise RuntimeError("Input arrays should have 4th dimension's shape equal to 5") - reordered = spher_harm[..., [2, 0, 1, 5, 6, 4, 7, 3]] + reordered = spher_harm[..., [2, 3, 1, 4, 0]] + + if orders == (1,): + if spher_harm.shape[-1] != 3: + raise RuntimeError(f"Input arrays should have 4th dimension's shape equal to 3 {spher_harm.shape}") + + reordered = spher_harm[..., [2, 0, 1]] return reordered @@ -132,7 +151,7 @@ def ge_basis(x, y, z, orders=(1, 2), shim_cs='LPI'): _check_basis_inputs(x, y, z, orders) # Create spherical harmonics from first to second order - all_orders = np.array(range(1, 3)) + all_orders = np.array([order for order in orders if order != 0]) spher_harm = scaled_spher_harm(x, y, z, all_orders, shim_cs=shim_cs) # The following matrix (8 x 5) refers to the following: @@ -163,31 +182,41 @@ def ge_basis(x, y, z, orders=(1, 2), shim_cs='LPI'): # [0.85228, 2.3916, -0.10486, 0.48776, -305.75]]) # Reorder according to GE convention: x, y, z, xy, zy, zx, X2-Y2, z2 - reordered_spher = _reorder_to_ge(spher_harm) + reordered_spher = _reorder_to_ge(spher_harm, orders) # Scale # Hz/cm2/A, -> uT/m2/A = order2_to_order2 * 1e6 * (100 ** 2) / (GYROMAGNETIC_RATIO * 1e6) # = order2_to_order2 * (100 ** 2) / GYROMAGNETIC_RATIO orders_to_order2_uT = order2_to_order2 * (100 ** 2) / GYROMAGNETIC_RATIO + # Order 2 scaled = np.zeros_like(reordered_spher) - for i_channel in range(reordered_spher.shape[3]): - if i_channel in [0, 1, 2]: - - # Rescale to unit-shims that are G/cm - # They are currently in uT/m - # 1G = 1e-4T, 1T = 1e4G - # uT/m --> G/cm = reordered_spher * (1/1e6) * 1e4 * 100 = reordered_sphere - scaled[..., i_channel] = reordered_spher[..., i_channel] - else: + index = 0 + if 1 in orders: + + # Rescale to unit-shims that are G/cm + # They are currently in uT/m + # 1G = 1e-4T, 1T = 1e4G + # uT/m --> G/cm = reordered_spher * (1/1e6) * 1e4 * 100 = reordered_sphere + scaled[..., 0:ORDER_INDEXES[1]] = reordered_spher[..., 0:ORDER_INDEXES[1]] + index += ORDER_INDEXES[1] + + if 2 in orders: + for i_channel in range(index, index + ORDER_INDEXES[2]): # Since reordered_spher contains the values of 1uT/m^2 in Hz/mm^2. We simply multiply by the amount of # uT/m^2 / A # This gives us a value in Hz/mm^2 / A which we need to modify to Hz/mm^2 / mA - scaled[..., i_channel] = np.matmul(reordered_spher[..., 3:], orders_to_order2_uT[i_channel - 3, :]) / 1000 + scaled[..., i_channel] = np.matmul(reordered_spher[..., index:index + ORDER_INDEXES[2]], + orders_to_order2_uT[i_channel - index, :]) / 1000 # Output according to the specified order - range_per_order = {1: list(range(3)), 2: list(range(3, 8))} + range_per_order = {} + index = 0 + for order in orders: + range_per_order[order] = list(range(index, index + ORDER_INDEXES[order])) + index += ORDER_INDEXES[order] + length_dim3 = np.sum([len(values) for key, values in range_per_order.items() if key in orders]) output = np.zeros(scaled[..., 0].shape + (length_dim3,), dtype=scaled.dtype) start_index = 0 @@ -200,7 +229,7 @@ def ge_basis(x, y, z, orders=(1, 2), shim_cs='LPI'): return output -def _reorder_to_ge(spher_harm): +def _reorder_to_ge(spher_harm, orders): """ Reorder 1st - 2nd order coefficients along the last dim. From 1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics), to @@ -214,15 +243,25 @@ def _reorder_to_ge(spher_harm): numpy.ndarray: Coefficients ordered following GE convention """ - if spher_harm.shape[-1] != 8: - raise RuntimeError("Input arrays should have 4th dimension's shape equal to 8") + if orders == (1, 2): + if spher_harm.shape[-1] != 8: + raise RuntimeError("Input arrays should have 4th dimension's shape equal to 8") + reordered = spher_harm[..., [2, 0, 1, 3, 4, 6, 7, 5]] + + if orders == (2,): + if spher_harm.shape[-1] != 5: + raise RuntimeError("Input arrays should have 4th dimension's shape equal to 5") + reordered = spher_harm[..., [0, 1, 3, 4, 2]] - reordered = spher_harm[..., [2, 0, 1, 3, 4, 6, 7, 5]] + if orders == (1,): + if spher_harm.shape[-1] != 3: + raise RuntimeError(f"Input arrays should have 4th dimension's shape equal to 3 {spher_harm.shape}") + reordered = spher_harm[..., [2, 0, 1]] return reordered -def scaled_spher_harm(x, y, z, orders=(1, 2), shim_cs='ras'): +def scaled_spher_harm(x, y, z, all_orders, shim_cs='ras'): """ The function first wraps ``shimmingtoolbox.coils.spherical_harmonics`` to generate 1st and 2nd order spherical harmonic ``basis`` fields at the grid positions given by arrays ``X,Y,Z``. It is then: @@ -246,22 +285,25 @@ def scaled_spher_harm(x, y, z, orders=(1, 2), shim_cs='ras'): numpy.ndarray: 4d basis set of spherical harmonics scaled """ # Check inputs - _check_basis_inputs(x, y, z, orders) + _check_basis_inputs(x, y, z, all_orders) # Create spherical harmonics from first to second order - all_orders = np.array([1, 2]) spher_harm = spherical_harmonics(all_orders, x, y, z) # 1. Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 (output by shimmingtoolbox.coils.spherical_harmonics.spherical_harmonics) - spher_harm_cs = get_flip_matrix(shim_cs) * spher_harm + # spher_harm_cs = _reorder_to_siemens(spher_harm, tuple(all_orders)) + spher_harm_cs = get_flip_matrix(shim_cs, all_orders) * spher_harm # scale according to # - 1 micro-T/m for *X,Y,Z* gradients (= 0.042576 Hz/mm) # - 1 micro-T/m^2 for 2nd order terms (= 0.000042576 Hz/mm^2) scaling_factors = _get_scaling_factors() scaled = np.zeros_like(spher_harm_cs) - for i_channel in range(0, spher_harm_cs.shape[3]): - scaled[:, :, :, i_channel] = scaling_factors[i_channel] * spher_harm_cs[:, :, :, i_channel] + index = 0 + for order in all_orders: + scaled[:, :, :, index:index + ORDER_INDEXES[order]] = np.array(scaling_factors[order]).squeeze() \ + * spher_harm_cs[:, :, :, index:index + ORDER_INDEXES[order]] + index += ORDER_INDEXES[order] # 1 uT/m, 1 uT/m2 in Hz/mm, Hz/mm2 return scaled @@ -309,12 +351,13 @@ def _get_scaling_factors(): r = [1, 1, 1, np.sqrt(2), np.sqrt(2), 1, np.sqrt(2), 1] # scaling: - orders = [1, 1, 1, 2, 2, 2, 2, 2] - - for i_ch in range(0, n_channels): - field = sh[:, :, :, i_ch] - scaling_factors[i_ch] = GYROMAGNETIC_RATIO * ((r[i_ch] * 0.001) ** orders[i_ch]) / field[i_ref[i_ch]] - + orders = [1, 2] + scaling_factors = {} + index = 0 + for order in orders: + scaling_factors[order] = [GYROMAGNETIC_RATIO * ((r[i_ch] * 0.001) ** order) + / sh[:, :, :, i_ch][i_ref[i_ch]] for i_ch in range(index, index + ORDER_INDEXES[order])] + index += ORDER_INDEXES[order] return scaling_factors @@ -330,7 +373,7 @@ def _check_basis_inputs(x, y, z, orders): raise NotImplementedError("Spherical harmonics not implemented for order 3 and up") -def get_flip_matrix(shim_cs='ras', manufacturer=None, xyz=False): +def get_flip_matrix(shim_cs='ras', orders=(1,2), manufacturer=None): """ Return a matrix to flip the spherical harmonics basis set from ras to the desired coordinate system. @@ -360,30 +403,31 @@ def get_flip_matrix(shim_cs='ras', manufacturer=None, xyz=False): if shim_cs[2] == 'I': xyz_cs[2] = -1 + temp_list_out = [] + if 1 in orders: + # Y, Z, X + out_matrix = temp_list_out.extend([xyz_cs[1], xyz_cs[2], xyz_cs[0]]) + if 2 in orders: + # XY, ZY, ZX, X2 - Y2 + out_matrix = temp_list_out.extend([xyz_cs[0] * xyz_cs[1], + xyz_cs[2] * xyz_cs[1], 1, + xyz_cs[2] * xyz_cs[0], 1]) # Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 - out_matrix = np.array([xyz_cs[1], xyz_cs[2], - xyz_cs[0], xyz_cs[0] * xyz_cs[1], - xyz_cs[2] * xyz_cs[1], 1, - xyz_cs[2] * xyz_cs[0], 1]) + out_matrix = np.array(temp_list_out) if manufacturer is not None: manufacturer = manufacturer.upper() - if manufacturer == 'SIEMENS': - out_matrix = _reorder_to_siemens(out_matrix) + out_matrix = _reorder_to_siemens(out_matrix, orders) elif manufacturer == 'GE': - out_matrix = _reorder_to_ge(out_matrix) + out_matrix = _reorder_to_ge(out_matrix, orders) elif manufacturer == 'PHILIPS': logger.warning("Philips shim CS not implemented yet") else: # Do not reorder if the manufacturer is not specified pass - if xyz: - # X, Y, Z - return out_matrix[:3] - else: - # None: Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 - # GE: x, y, z, xy, zy, zx, X2 - Y2, z2 - # Siemens: X, Y, Z, Z2, ZX, ZY, X2 - Y2, XY - return out_matrix + # None: Y, Z, X, XY, ZY, Z2, ZX, X2 - Y2 + # GE: x, y, z, xy, zy, zx, X2 - Y2, z2 + # Siemens: X, Y, Z, Z2, ZX, ZY, X2 - Y2, XY + return out_matrix diff --git a/shimmingtoolbox/optimizer/basic_optimizer.py b/shimmingtoolbox/optimizer/basic_optimizer.py index 517e1ebf..13c3a4a1 100644 --- a/shimmingtoolbox/optimizer/basic_optimizer.py +++ b/shimmingtoolbox/optimizer/basic_optimizer.py @@ -81,7 +81,7 @@ def set_merged_bounds(self, merged_bounds): merged_bounds: Concatenated coil profile bounds """ if len(self.merged_bounds) != len(merged_bounds): - raise ValueError(f"Size of merged bounds: must match the number of total channel: {len(self.merged_bounds)}") + raise ValueError(f"Size of merged bounds: must match the number of total channel: {len(self.merged_bounds)} not {len(merged_bounds)}") self.merged_bounds = merged_bounds def optimize(self, mask): @@ -177,8 +177,9 @@ def merge_bounds(self): bounds = [] for coil in self.coils: # Concat coils and bounds - for a_bound in coil.coef_channel_minmax: - bounds.append(a_bound) + for key in coil.coef_channel_minmax: + for a_bound in coil.coef_channel_minmax[key]: + bounds.append(a_bound) return bounds diff --git a/shimmingtoolbox/shim/sequencer.py b/shimmingtoolbox/shim/sequencer.py index 2b2c00cd..68cc17b8 100644 --- a/shimmingtoolbox/shim/sequencer.py +++ b/shimmingtoolbox/shim/sequencer.py @@ -798,9 +798,9 @@ class RealTimeSequencer(Sequencer): acq_timestamps (np.ndarray) : 1D array that contains the acquisitions timestamps """ - def __init__(self, nii_fieldmap, json_fmap, nii_anat, nii_static_mask, nii_riro_mask, slices, pmu: PmuResp, coils, - method='least_squares', opt_criteria='mse', mask_dilation_kernel='sphere', mask_dilation_kernel_size=3, - reg_factor=0, path_output=None): + def __init__(self, nii_fieldmap, json_fmap, nii_anat, nii_static_mask, nii_riro_mask, slices, pmu: PmuResp, coils_static, + coils_riro, method='least_squares', opt_criteria='mse', mask_dilation_kernel='sphere', + mask_dilation_kernel_size=3, reg_factor=0, path_output=None): """ Initialization of the RealTimeSequencer class @@ -843,7 +843,8 @@ def __init__(self, nii_fieldmap, json_fmap, nii_anat, nii_static_mask, nii_riro_ super().__init__(slices, mask_dilation_kernel, mask_dilation_kernel_size, reg_factor, path_output) self.json_fmap = json_fmap self.pmu = pmu - self.coils = coils + self.coils_static = coils_static + self.coils_riro = coils_riro self.method = method self.bounds = None @@ -1062,7 +1063,8 @@ def shim(self): # RIRO optimization # Use the currents to define a list of new coil bounds for the riro optimization - self.bounds = new_bounds_from_currents(coef_static, self.optimizer.merged_bounds) + self.bounds = new_bounds_from_currents_static_to_riro(coef_static, self.optimizer.merged_bounds, self.coils_static, + self.coils_riro) logger.info("Realtime optimization") coef_riro = self.optimize_riro(riro_mask_resampled) @@ -1089,10 +1091,10 @@ def select_optimizer(self, unshimmed, affine, pmu: PmuResp = None): # global supported_optimizers if self.method in supported_optimizers: if self.method == 'least_squares': - self.optimizer = supported_optimizers[self.method](self.coils, unshimmed, affine, self.opt_criteria, + self.optimizer = supported_optimizers[self.method](self.coils_static, unshimmed, affine, self.opt_criteria, reg_factor=self.reg_factor) elif self.method == 'quad_prog': - self.optimizer = supported_optimizers[self.method](self.coils, unshimmed, affine, + self.optimizer = supported_optimizers[self.method](self.coils_static, unshimmed, affine, reg_factor=self.reg_factor) elif self.method == 'least_squares_rt': @@ -1101,7 +1103,7 @@ def select_optimizer(self, unshimmed, affine, pmu: PmuResp = None): raise ValueError(f"pmu parameter is required if using the optimization method: {self.method}") # Add pmu to the realtime optimizer(s) - self.optimizer_riro = supported_optimizers[self.method](self.coils, unshimmed, affine, + self.optimizer_riro = supported_optimizers[self.method](self.coils_riro, unshimmed, affine, self.opt_criteria, pmu, reg_factor=self.reg_factor) elif self.method == 'quad_prog_rt': @@ -1110,14 +1112,14 @@ def select_optimizer(self, unshimmed, affine, pmu: PmuResp = None): raise ValueError(f"pmu parameter is required if using the optimization method: {self.method}") # Add pmu to the realtime optimizer(s) - self.optimizer_riro = supported_optimizers[self.method](self.coils, unshimmed, affine, pmu, + self.optimizer_riro = supported_optimizers[self.method](self.coils_riro, unshimmed, affine, pmu, reg_factor=self.reg_factor) else: if pmu is None: - self.optimizer = supported_optimizers[self.method](self.coils, unshimmed, affine) + self.optimizer = supported_optimizers[self.method](self.coils_static, unshimmed, affine) else: - self.optimizer_riro = supported_optimizers[self.method](self.coils, unshimmed, affine) + self.optimizer_riro = supported_optimizers[self.method](self.coils_riro, unshimmed, affine) else: raise KeyError(f"Method: {self.method} is not part of the supported optimizers") @@ -1223,7 +1225,7 @@ def eval(self, coef_static, coef_riro, mean_p, pressure_rms): cval=0).get_fdata()), 0, 1) for i_shim in range(len(self.slices)): # Calculate static correction - correction_static = self.optimizer_riro.merged_coils @ coef_static[i_shim] + correction_static = self.optimizer.merged_coils @ coef_static[i_shim] # Calculate the riro coil profiles riro_profile = self.optimizer_riro.merged_coils @ coef_riro[i_shim] @@ -1700,33 +1702,103 @@ def plot_full_time_std(self, unshimmed, masked_shim_static_riro, mask_fmap_cs, m fig.savefig(fname_figure, bbox_inches='tight') -def new_bounds_from_currents(currents, old_bounds): +def new_bounds_from_currents(currents:dict, old_bounds:dict): + """ + Uses the currents to determine the appropriate bounds for the next optimization. It assumes that + "old_coef + next_bound < old_bound". + + Args: + currents (dict): Dictionary with n_shims as keys each with a list of n_channels values. + old_bounds (dict): Dictionary with orders as keys containing (min, max) containing the merged bounds of the previous + optimization. + Returns: + dict: Modified bounds (same shape as old_bounds) + """ + new_bounds = {} + for key in old_bounds: + new_bounds[key] = [] + for i, bound in enumerate(old_bounds[key]): + if bound == [None, None]: + new_bounds[key].append(bound) + elif bound[0] is None: + new_bounds[key].append([None, bound[1] - currents[key][i]]) + elif bound[1] is None: + new_bounds[key].append([bound[0] - currents[key][i], None]) + else: + new_bounds[key].append([bound[0] - currents[key][i], bound[1] - currents[key][i]]) + return new_bounds + + +def new_bounds_from_currents_static_to_riro(currents, old_bounds, coils_static=[], coils_riro=[]): + #! Modify description if everything works """ Uses the currents to determine the appropriate bounds for the next optimization. It assumes that "old_coef + next_bound < old_bound". Args: - currents (np.ndarray): 2D array (n_shims x n_channels). - old_bounds (list): 2D list (n_channels, 2) containing (min, max) containing the merged bounds of the previous + currents (np.ndarray): 2D array (n_shims x n_channels). Direct output from :func:`_optimize`. + old_bounds (list): 1d list (n_channels) of tuples (min, max) containing the merged bounds of the previous optimization. Returns: list: 2d list (n_shim_groups x n_channels) of bounds (min, max) corresponding to each shim group and channel. """ + currents_riro = np.empty((currents.shape[0], 0)) + old_bounds_riro = [] + static_coil_names = [c.name for c in coils_static] + + index = 0 + coil_indexes = {} + for coil in coils_static: + if type(coil) == Coil: + coil_indexes[coil.name] = [index, index + len(coil.coef_channel_minmax['coil'])] + index += len(coil.coef_channel_minmax['coil']) + else: + coil_indexes[coil.name] = {} + for key in coil.coef_channel_minmax: + coil_indexes[coil.name][key] = [index, index + len(coil.coef_channel_minmax[key])] + index += len(coil.coef_channel_minmax[key]) + + for i, coil in enumerate(coils_riro): + if coil.name in static_coil_names: + if type(coil) == Coil: + currents_riro = np.append(currents_riro, + currents[:, coil_indexes[coil.name][0]:coil_indexes[coil.name][1]], + axis=1) + old_bounds_riro.extend(old_bounds[coil_indexes[coil.name][0]:coil_indexes[coil.name][1]]) + else: + for order in coil.coef_channel_minmax: + if order in coils_static[static_coil_names.index(coil.name)].coef_channel_minmax.keys(): + currents_riro = np.append(currents_riro, + currents[:, coil_indexes[coil.name][order][0]:coil_indexes[coil.name][order][1]], + axis = 1) + old_bounds_riro.extend(old_bounds[coil_indexes[coil.name][order][0]:coil_indexes[coil.name][order][1]]) + else: + currents_riro = np.append(currents_riro, + np.zeros((currents.shape[0], len(coil.coef_channel_minmax[order]))), + axis=1) + old_bounds_riro.extend(coil.coef_channel_minmax[order]) + + else: + if type(coil) == Coil: + currents_riro = np.append(currents_riro, + np.zeros((currents.shape[0], len(coil.coef_channel_minmax['coil']))), + axis=1) + old_bounds_riro.extend(coil.coef_channel_minmax['coil']) + + else: + for order in coil.coef_channel_minmax: + currents_riro = np.append(currents_riro, + np.zeros((currents.shape[0], len(coil.coef_channel_minmax[order]))), + axis=1) + old_bounds_riro.extend(coil.coef_channel_minmax[order]) + new_bounds = [] - for i_shim in range(currents.shape[0]): + for i_shim in range(currents_riro.shape[0]): shim_bound = [] - for i_channel in range(len(old_bounds)): - if old_bounds[i_channel] == [None, None]: - a_bound = old_bounds[i_channel] - elif old_bounds[i_channel][0] is None: - a_bound = [None, old_bounds[i_channel][1] - currents[i_shim, i_channel]] - elif old_bounds[i_channel][1] is None: - a_bound = [old_bounds[i_channel][0] - currents[i_shim, i_channel], None] - else: - a_bound = [old_bounds[i_channel][0] - currents[i_shim, i_channel], - old_bounds[i_channel][1] - currents[i_shim, i_channel]] - shim_bound.append(a_bound) + for i_channel in range(len(old_bounds_riro)): + a_bound = old_bounds_riro[i_channel] - currents_riro[i_shim, i_channel] + shim_bound.append(tuple(a_bound)) new_bounds.append(shim_bound) return new_bounds diff --git a/shimmingtoolbox/shim/shim_utils.py b/shimmingtoolbox/shim/shim_utils.py index f7c358d5..adfa5123 100644 --- a/shimmingtoolbox/shim/shim_utils.py +++ b/shimmingtoolbox/shim/shim_utils.py @@ -156,7 +156,7 @@ def calculate_metric_within_mask(array, mask, metric='mean', axis=None): return output -def phys_to_shim_cs(coefs, manufacturer): +def phys_to_shim_cs(coefs, manufacturer, orders): """Convert a list of coefficients from RAS to the Shim Coordinate System Args: @@ -165,6 +165,7 @@ def phys_to_shim_cs(coefs, manufacturer): more coefficients, they are of higher order and must correspond to the implementation of the manufacturer. i.e. Siemens: *X, Y, Z, Z2, ZX, ZY, X2-Y2, XY* manufacturer (str): Name of the manufacturer + orders (tuple): Tuple containing the spherical harmonic orders Returns: np.ndarray: Coefficients in the shim coordinate system of the manufacturer @@ -175,10 +176,10 @@ def phys_to_shim_cs(coefs, manufacturer): flip_mat = np.ones(len(coefs)) # Order 1 if len(coefs) == 3: - flip_mat[:3] = get_flip_matrix(shim_cs[manufacturer], manufacturer=manufacturer, xyz=True) + flip_mat[:3] = get_flip_matrix(shim_cs[manufacturer], orders, manufacturer=manufacturer) # Order 2 elif len(coefs) >= 8: - flip_mat[:8] = get_flip_matrix(shim_cs[manufacturer], manufacturer=manufacturer) + flip_mat[:8] = get_flip_matrix(shim_cs[manufacturer], orders, manufacturer=manufacturer) else: logger.warning("Order not supported") @@ -190,7 +191,7 @@ def phys_to_shim_cs(coefs, manufacturer): return coefs -def shim_to_phys_cs(coefs, manufacturer): +def shim_to_phys_cs(coefs, manufacturer, orders): """ Convert coefficients from the shim coordinate system to the physical RAS coordinate system Args: @@ -199,12 +200,15 @@ def shim_to_phys_cs(coefs, manufacturer): more coefficients, they are of higher order and must correspond to the implementation of the manufacturer. Siemens: *X, Y, Z, Z2, ZX, ZY, X2-Y2, XY* manufacturer (str): Name of the manufacturer + orders (tuple): Tuple containing the spherical harmonic orders Returns: np.ndarray: Coefficients in the physical RAS coordinate system + """ + # It's sign flips so the same function can be used for shimCS <--> phys RAS - coefs = phys_to_shim_cs(coefs, manufacturer) + coefs = phys_to_shim_cs(coefs, manufacturer, orders) return coefs @@ -222,26 +226,26 @@ def get_scanner_shim_settings(bids_json_dict): """ scanner_shim = { - 'f0': None, - 'order1': None, - 'order2': None, - 'order3': None, + '0': None, + '1': None, + '2': None, + '3': None, 'has_valid_settings': False } # get_imaging_frequency if bids_json_dict.get('ImagingFrequency'): - scanner_shim['f0'] = int(bids_json_dict.get('ImagingFrequency') * 1e6) + scanner_shim['0'] = int(bids_json_dict.get('ImagingFrequency') * 1e6) # get_shim_orders if bids_json_dict.get('ShimSetting'): n_shim_values = len(bids_json_dict.get('ShimSetting')) if n_shim_values == 3: - scanner_shim['order1'] = bids_json_dict.get('ShimSetting') + scanner_shim['1'] = bids_json_dict.get('ShimSetting') scanner_shim['has_valid_settings'] = True elif n_shim_values == 8: - scanner_shim['order2'] = bids_json_dict.get('ShimSetting')[3:] - scanner_shim['order1'] = bids_json_dict.get('ShimSetting')[:3] + scanner_shim['2'] = bids_json_dict.get('ShimSetting')[3:] + scanner_shim['1'] = bids_json_dict.get('ShimSetting')[:3] scanner_shim['has_valid_settings'] = True else: logger.warning(f"ShimSetting tag has an unsupported number of values: {n_shim_values}") @@ -269,28 +273,28 @@ def convert_to_mp(manufacturers_model_name, shim_settings): scanner_shim_mp = shim_settings if manufacturers_model_name == "Prisma_fit": - if shim_settings.get('order1'): + if shim_settings.get('1'): # One can use the Siemens commandline AdjValidate tool to get all the values below: max_current_mp_order1 = np.array([2300] * 3) max_current_dcm_order1 = np.array([14436, 14265, 14045]) - order1_mp = np.array(shim_settings['order1']) * max_current_mp_order1 / max_current_dcm_order1 + order1_mp = np.array(shim_settings['1']) * max_current_mp_order1 / max_current_dcm_order1 if np.any(np.abs(order1_mp) > max_current_mp_order1): scanner_shim_mp['has_valid_settings'] = False raise ValueError("Multipole values exceed known system limits.") else: - scanner_shim_mp['order1'] = order1_mp + scanner_shim_mp['1'] = order1_mp - if shim_settings.get('order2'): + if shim_settings.get('2'): max_current_mp_order2 = np.array([4959.01, 3551.29, 3503.299, 3551.29, 3487.302]) max_current_dcm_order2 = np.array([9998, 9998, 9998, 9998, 9998]) - order2_mp = np.array(shim_settings['order2']) * max_current_mp_order2 / max_current_dcm_order2 + order2_mp = np.array(shim_settings['2']) * max_current_mp_order2 / max_current_dcm_order2 if np.any(np.abs(order2_mp) > max_current_mp_order2): scanner_shim_mp['has_valid_settings'] = False raise ValueError("Multipole values exceed known system limits.") else: - scanner_shim_mp['order2'] = order2_mp + scanner_shim_mp['2'] = order2_mp else: logger.debug(f"Manufacturer model {manufacturers_model_name} not implemented, could not convert shim settings") @@ -304,24 +308,24 @@ def __init__(self, bids_json_dict): shim_settings_dac = get_scanner_shim_settings(bids_json_dict) self.shim_settings = convert_to_mp(bids_json_dict.get('ManufacturersModelName'), shim_settings_dac) - def concatenate_shim_settings(self, order=2): + def concatenate_shim_settings(self, orders=[2]): coefs = [] if not self.shim_settings['has_valid_settings']: logger.warning("Invalid Shim Settings") return coefs - if self.shim_settings.get('f0') is not None and order >= 0: + if self.shim_settings.get('0') is not None and any(order>=0 for order in orders): # Concatenate 2 lists - coefs = [self.shim_settings.get('f0')] + coefs = [self.shim_settings.get('0')] else: coefs = [0] - for i_order in range(1, order + 1): - if self.shim_settings.get(f'order{i_order}') is not None: + for order in orders: + if self.shim_settings.get(order) is not None: # Concatenate 2 lists - coefs.extend(self.shim_settings.get(f'order{i_order}')) + coefs.extend(self.shim_settings.get(order)) else: - n_coefs = (i_order + 1) * 2 + n_coefs = (order + 1) * 2 coefs.extend([0] * n_coefs) return coefs diff --git a/test/cli/test_cli_b0shim.py b/test/cli/test_cli_b0shim.py index 38873d18..fb30cb4b 100644 --- a/test/cli/test_cli_b0shim.py +++ b/test/cli/test_cli_b0shim.py @@ -10,6 +10,7 @@ import nibabel as nib import numpy as np import json +import csv from shimmingtoolbox.cli.b0shim import define_slices_cli from shimmingtoolbox.cli.b0shim import b0shim_cli @@ -88,13 +89,19 @@ def test_cli_dynamic_sph(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data) '--fmap', fname_fmap, '--anat', fname_anat, '--mask', fname_mask, - '--scanner-coil-order', '2', + '--scanner-coil-order', '1,2', '--regularization-factor', '0.1', '--output', tmp], catch_exceptions=False) assert res.exit_code == 0 assert os.path.isfile(os.path.join(tmp, "coefs_coil0_Prisma_fit.txt")) + with open(os.path.join(tmp, "coefs_coil0_Prisma_fit.txt"), 'r') as file: + lines = file.readlines() + line = lines[8].strip().split(',') + values = [float(val) for val in line if val.strip()] + + assert values == [0.002985, -14.587414, -57.016499, -2.745062, -0.401786, -3.580623, 0.668977, -0.105560] def test_cli_dynamic_no_mask(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): """Test cli with scanner coil profiles of order 1 with default constraints""" @@ -113,7 +120,7 @@ def test_cli_dynamic_no_mask(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_d res = runner.invoke(b0shim_cli, ['dynamic', '--fmap', fname_fmap, '--anat', fname_anat, - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--output', tmp], catch_exceptions=False) @@ -311,7 +318,7 @@ def test_cli_dynamic_format_chronological_ch(self, nii_fmap, nii_anat, nii_mask, '--fmap', fname_fmap, '--anat', fname_anat, '--mask', fname_mask, - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--slice-factor', '2', '--output-file-format-scanner', 'chronological-ch', '--output', tmp], @@ -344,7 +351,7 @@ def test_cli_dynamic_format_slicewise_ch(self, nii_fmap, nii_anat, nii_mask, fm_ '--fmap', fname_fmap, '--anat', fname_anat, '--mask', fname_mask, - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--slice-factor', '2', '--output-file-format-scanner', 'slicewise-ch', '--output', tmp], @@ -376,7 +383,7 @@ def test_cli_dynamic_format_gradient(self, nii_fmap, nii_anat, nii_mask, fm_data '--fmap', fname_fmap, '--anat', fname_anat, '--mask', fname_mask, - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--slice-factor', '2', '--output-file-format-scanner', 'gradient', '--output', tmp], @@ -586,7 +593,7 @@ def test_cli_rt_sph(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): '--mask-riro', fname_mask, '--resp', fname_resp, '--slice-factor', '2', - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--output', tmp, '-v', 'debug'], catch_exceptions=False) @@ -596,6 +603,11 @@ def test_cli_rt_sph(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch1_Prisma_fit.txt")) assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch2_Prisma_fit.txt")) assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch3_Prisma_fit.txt")) + with open(os.path.join(tmp, "coefs_coil0_ch0_Prisma_fit.txt"), 'r') as file: + lines = file.readlines() + line = lines[5].strip().split(',') + values = [float(val) for val in line if val.strip()] + assert values == [11.007908, -0.014577058094, 1311.6784] def test_cli_rt_sph_order_0(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp: @@ -630,6 +642,44 @@ def test_cli_rt_sph_order_0(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_da assert res.exit_code == 0 assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch0_Prisma_fit.txt")) + def test_cli_rt_sph_order_02(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): + with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp: + # Save the inputs to the new directory + fname_fmap = os.path.join(tmp, 'fmap.nii.gz') + fname_fm_json = os.path.join(tmp, 'fmap.json') + fname_mask = os.path.join(tmp, 'mask.nii.gz') + fname_anat = os.path.join(tmp, 'anat.nii.gz') + fname_anat_json = os.path.join(tmp, 'anat.json') + _save_inputs(nii_fmap=nii_fmap, fname_fmap=fname_fmap, + nii_anat=nii_anat, fname_anat=fname_anat, + nii_mask=nii_mask, fname_mask=fname_mask, + fm_data=fm_data, fname_fm_json=fname_fm_json, + anat_data=anat_data, fname_anat_json=fname_anat_json) + + # Input pmu fname + fname_resp = os.path.join(__dir_testing__, 'ds_b0', 'derivatives', 'sub-realtime', + 'sub-realtime_PMUresp_signal.resp') + + runner = CliRunner() + res = runner.invoke(b0shim_cli, ['realtime-dynamic', + '--fmap', fname_fmap, + '--anat', fname_anat, + '--mask-static', fname_mask, + '--mask-riro', fname_mask, + '--resp', fname_resp, + '--slice-factor', '2', + '--scanner-coil-order', '0,2', + '--output', tmp], + catch_exceptions=False) + + assert res.exit_code == 0 + assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch0_Prisma_fit.txt")) + assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch4_Prisma_fit.txt")) + assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch5_Prisma_fit.txt")) + assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch6_Prisma_fit.txt")) + assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch7_Prisma_fit.txt")) + assert os.path.isfile(os.path.join(tmp, "coefs_coil0_ch8_Prisma_fit.txt")) + def test_cli_rt_custom_coil(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): with tempfile.TemporaryDirectory(prefix='st_' + pathlib.Path(__file__).stem) as tmp: # Save the inputs to the new directory @@ -661,6 +711,7 @@ def test_cli_rt_custom_coil(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_da runner = CliRunner() res = runner.invoke(b0shim_cli, ['realtime-dynamic', '--coil', fname_dummy_coil, fname_constraints, + '--coil-riro', fname_dummy_coil, fname_constraints, '--fmap', fname_fmap, '--anat', fname_anat, '--mask-static', fname_mask, @@ -700,7 +751,7 @@ def test_cli_rt_debug(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): '--mask-riro', fname_mask, '--resp', fname_resp, '--slice-factor', '2', - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '-v', 'debug', '--output', tmp], catch_exceptions=False) @@ -739,7 +790,7 @@ def test_cli_rt_no_mask(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data): '--anat', fname_anat, '--resp', fname_resp, '--slice-factor', '2', - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--output', tmp], catch_exceptions=False) @@ -775,7 +826,7 @@ def test_cli_rt_chronological_ch(self, nii_fmap, nii_anat, nii_mask, fm_data, an '--mask-riro', fname_mask, '--resp', fname_resp, '--slice-factor', '2', - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--output-file-format-scanner', 'chronological-ch', '--output', tmp], catch_exceptions=False) @@ -812,7 +863,7 @@ def test_cli_rt_gradient(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data) '--mask-riro', fname_mask, '--resp', fname_resp, '--slice-factor', '2', - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--output-file-format-scanner', 'gradient', '--output', tmp], catch_exceptions=False) @@ -849,7 +900,7 @@ def test_cli_rt_absolute(self, nii_fmap, nii_anat, nii_mask, fm_data, anat_data) '--mask-riro', fname_mask, '--resp', fname_resp, '--slice-factor', '2', - '--scanner-coil-order', '2', + '--scanner-coil-order', '0,1,2', '--output-value-format', 'absolute', '--output', tmp], catch_exceptions=False) @@ -886,7 +937,7 @@ def test_cli_rt_pseudo_inverse(self, nii_fmap, nii_anat, nii_mask, fm_data, anat '--mask-riro', fname_mask, '--resp', fname_resp, '--slice-factor', '2', - '--scanner-coil-order', '1', + '--scanner-coil-order', '0,1', '--optimizer-method', 'pseudo_inverse', '--output', tmp], catch_exceptions=False) @@ -1066,8 +1117,8 @@ def _create_dummy_coil(nii_fmap): # Dummy constraints dummy_coil_constraints = { "name": "Dummy_coil", - "coef_channel_minmax": [(-6000, 6000), (-2405, 2194), (-1120, 3479), (-754, 3845), (-4252.2, 5665.8), - (-3692, 3409), (-3417, 3588), (-3568, 3534), (-3483, 3490)], + "coef_channel_minmax": {"coil": [(-6000, 6000), (-2405, 2194), (-1120, 3479), (-754, 3845), (-4252.2, 5665.8), + (-3692, 3409), (-3417, 3588), (-3568, 3534), (-3483, 3490)]}, "coef_sum_max": None } diff --git a/test/cli/test_cli_create_coil_profiles.py b/test/cli/test_cli_create_coil_profiles.py index 1dda1ceb..386271d0 100644 --- a/test/cli/test_cli_create_coil_profiles.py +++ b/test/cli/test_cli_create_coil_profiles.py @@ -181,8 +181,7 @@ def test_create_coil_profiles_from_cad(): '--offset', '0', '-111', '-47', '-o', fname_output], catch_exceptions=False) - print(fname_txt) - print(res.output) + assert res.exit_code == 0 fname_cp = os.path.join(fname_output, coil_name + '_coil_profiles.nii.gz') @@ -197,4 +196,5 @@ def test_create_coil_profiles_from_cad(): with open(fname_config, 'rb') as f: config_test = json.load(f) + assert are_jsons_equal(config_test, ref_config) diff --git a/test/shim/test_sequencer.py b/test/shim/test_sequencer.py index fdf9968c..dbd673b6 100644 --- a/test/shim/test_sequencer.py +++ b/test/shim/test_sequencer.py @@ -66,9 +66,9 @@ def create_unshimmed_affine(): def create_constraints(max_coef, min_coef, sum_coef, n_channels=8): # Set up bounds for output currents - bounds = [] + bounds = {"coil": []} for _ in range(n_channels): - bounds.append((min_coef, max_coef)) + bounds['coil'].append((min_coef, max_coef)) constraints = { "name": "test", @@ -414,7 +414,7 @@ def test_shim_realtime_pmu_sequencer_fake_data(self, nii_fieldmap, json_data, ni # Find optimal currents sequencer_realtime_test = RealTimeSequencer(nii_fieldmap, json_data, nii_anat, nii_mask_static, nii_mask_riro, - slices, pmu, [coil], + slices, pmu, [coil], [coil], mask_dilation_kernel='sphere') output = sequencer_realtime_test.shim() currents_static, currents_riro, mean_p, p_rms = output @@ -496,7 +496,7 @@ def test_shim_sequencer_rt_larger_coil(self, nii_fieldmap, json_data, nii_anat, # Find optimal currents output = RealTimeSequencer(nii_fieldmap, json_data, nii_anat, nii_mask_static, nii_mask_riro, - slices, pmu, [new_coil]).shim() + slices, pmu, [new_coil], [new_coil]).shim() currents_static, currents_riro, mean_p, p_rms = output print(f"\nSlices: {slices}" @@ -511,7 +511,7 @@ def test_shim_sequencer_rt_kernel_line(self, nii_fieldmap, json_data, nii_anat, nii_mask_riro, slices, pmu, coil): # Optimize output = RealTimeSequencer(nii_fieldmap, json_data, nii_anat, nii_mask_static, nii_mask_riro, - slices, pmu, [coil], mask_dilation_kernel='line').shim() + slices, pmu, [coil], [coil], mask_dilation_kernel='line').shim() assert output[0].shape == (20, 3) @@ -522,7 +522,7 @@ def test_shim_sequencer_rt_wrong_fmap_dim(self, nii_fieldmap, json_data, nii_ana header=nii_fieldmap.header) with pytest.raises(ValueError, match="Fieldmap must be 4d"): RealTimeSequencer(nii_wrong_fmap, json_data, nii_anat, nii_mask_static, nii_mask_riro, - slices, pmu, [coil]).shim() + slices, pmu, [coil], [coil]).shim() def test_shim_sequencer_rt_wrong_anat_dim(self, nii_fieldmap, json_data, nii_anat, nii_mask_static, nii_mask_riro, slices, pmu, coil): @@ -530,7 +530,7 @@ def test_shim_sequencer_rt_wrong_anat_dim(self, nii_fieldmap, json_data, nii_ana nii_wrong_anat = nib.Nifti1Image(nii_anat.get_fdata()[..., 0], nii_anat.affine, header=nii_anat.header) with pytest.raises(ValueError, match="Anatomical image must be in 3d"): RealTimeSequencer(nii_fieldmap, json_data, nii_wrong_anat, nii_mask_static, nii_mask_riro, - slices, pmu, [coil]).shim() + slices, pmu, [coil], [coil]).shim() def test_shim_sequencer_rt_diff_mask_shape_static(self, nii_fieldmap, json_data, nii_anat, nii_mask_static, nii_mask_riro, slices, pmu, coil): @@ -538,7 +538,7 @@ def test_shim_sequencer_rt_diff_mask_shape_static(self, nii_fieldmap, json_data, nii_diff_mask = nib.Nifti1Image(nii_mask_static.get_fdata()[5:, ...], nii_mask_static.affine, header=nii_mask_static.header) output = RealTimeSequencer(nii_fieldmap, json_data, nii_anat, nii_diff_mask, nii_mask_riro, - slices, pmu, [coil]).shim() + slices, pmu, [coil], [coil]).shim() assert output[0].shape == (20, 3) def test_shim_sequencer_rt_diff_mask_shape_riro(self, nii_fieldmap, json_data, nii_anat, nii_mask_static, @@ -548,7 +548,7 @@ def test_shim_sequencer_rt_diff_mask_shape_riro(self, nii_fieldmap, json_data, n header=nii_mask_riro.header) output = RealTimeSequencer(nii_fieldmap, json_data, nii_anat, nii_mask_static, nii_diff_mask, - slices, pmu, [coil]).shim() + slices, pmu, [coil], [coil]).shim() assert output[0].shape == (20, 3) def test_shim_sequencer_rt_diff_mask_affine(self, nii_fieldmap, json_data, nii_anat, nii_mask_static, @@ -558,7 +558,7 @@ def test_shim_sequencer_rt_diff_mask_affine(self, nii_fieldmap, json_data, nii_a diff_affine[0, 0] = 2 nii_diff_mask = nib.Nifti1Image(nii_mask_static.get_fdata(), diff_affine, header=nii_mask_static.header) output = RealTimeSequencer(nii_fieldmap, json_data, nii_anat, nii_diff_mask, nii_mask_riro, - slices, pmu, [coil]).shim() + slices, pmu, [coil], [coil]).shim() assert output[0].shape == (20, 3) @@ -604,7 +604,7 @@ def test_shim_realtime_pmu_sequencer_rt_zshim_data(): # Find optimal currents output = RealTimeSequencer(nii_fieldmap, json_data, nii_anat, nii_mask_static, nii_mask_riro, slices, pmu, - [coil], method='least_squares').shim() + [coil], [coil], method='least_squares').shim() currents_static, currents_riro, mean_p, p_rms = output # Scale according to rms diff --git a/test/shim/test_shim_utils.py b/test/shim/test_shim_utils.py index 3bd8c28e..718c3469 100644 --- a/test/shim/test_shim_utils.py +++ b/test/shim/test_shim_utils.py @@ -18,17 +18,17 @@ def test_convert_to_mp_unknown_scanner(caplog): def test_convert_to_mp_outside_bounds(): - dac_units = {'order1': [20000, 14265, 14045], 'order2': [9998, 9998, 9998, 9998, 9998]} + dac_units = {'1': [20000, 14265, 14045], '2': [9998, 9998, 9998, 9998, 9998]} with pytest.raises(ValueError, match="Multipole values exceed known system limits."): convert_to_mp('Prisma_fit', dac_units) def test_phys_to_shim_cs(): - out = phys_to_shim_cs(np.array([1, 1, 1]), 'Siemens') + out = phys_to_shim_cs(np.array([1, 1, 1]), 'Siemens', orders=(1,)) assert np.all(out == [-1, 1, -1]) def test_shim_to_phys_cs(): - out = shim_to_phys_cs(np.array([1, 1, 1]), 'Siemens') + out = shim_to_phys_cs(np.array([1, 1, 1]), 'Siemens', orders=(1,)) assert np.all(out == [-1, 1, -1]) diff --git a/test/test_coil.py b/test/test_coil.py index 58eccd52..c406f2f8 100644 --- a/test/test_coil.py +++ b/test/test_coil.py @@ -10,7 +10,6 @@ def test_coil_siemens_basis(): - grid_x, grid_y, grid_z = np.meshgrid(np.array(range(-1, 2)), np.array(range(-1, 2)), np.array(range(-1, 2)), indexing='ij') profiles = siemens_basis(grid_x, grid_y, grid_z) @@ -18,7 +17,7 @@ def test_coil_siemens_basis(): constraints = { "name": "Siemens Basis", "coef_sum_max": 40, - "coef_channel_minmax": [(-2, 2), (-2, 2), (-2, 2), (-2, 2), (-3, 3), (-3, 3), (-3, 3), (-3, 3)], + "coef_channel_minmax": {"coil": [(-2, 2), (-2, 2), (-2, 2), (-2, 2), (-3, 3), (-3, 3), (-3, 3), (-3, 3)]}, } a_coil = Coil(profiles, np.eye(4), constraints) @@ -36,7 +35,7 @@ def test_coil_custom_coil(): def test_create_scanner_coil_order0(): sph_contraints = json.load(open(__dir_config_scanner_constraints__)) - scanner_coil = ScannerCoil((4, 5, 6), np.eye(4), sph_contraints, 0) + scanner_coil = ScannerCoil((4, 5, 6), np.eye(4), sph_contraints, [0]) assert scanner_coil.profile[0, 0, 0, 0] == -1.0 @@ -44,13 +43,13 @@ def test_create_scanner_coil_order0(): def test_create_scanner_coil_order1(): sph_contraints = json.load(open(__dir_config_scanner_constraints__)) - scanner_coil = ScannerCoil((4, 5, 6), np.eye(4), sph_contraints, 1) + scanner_coil = ScannerCoil((4, 5, 6), np.eye(4), sph_contraints, [0, 1]) assert scanner_coil.profile[0, 0, 0, 0] == -1.0 def test_create_scanner_coil_order2(): sph_contraints = json.load(open(__dir_config_scanner_constraints__)) - scanner_coil = ScannerCoil((4, 5, 6), np.eye(4), sph_contraints, 2) + scanner_coil = ScannerCoil((4, 5, 6), np.eye(4), sph_contraints, [0, 1, 2]) assert scanner_coil.profile[0, 0, 0, 0] == -1.0 diff --git a/test/test_spher_harm_basis.py b/test/test_spher_harm_basis.py index 4ed97bbe..08cd2819 100644 --- a/test/test_spher_harm_basis.py +++ b/test/test_spher_harm_basis.py @@ -83,15 +83,15 @@ def test_ge_basis(x, y, z): class TestGetFlipMatrix: def test_flip_cs(self): - out = get_flip_matrix('RAS', xyz=True) + out = get_flip_matrix('RAS', orders=(1,)) assert np.all(out == [1, 1, 1]) def test_flip_cs_lpi(self): - out = get_flip_matrix('LPI', xyz=True) + out = get_flip_matrix('LPI', orders=(1,)) assert np.all(out == [-1, -1, -1]) def test_flip_cs_order2(self): - out = get_flip_matrix('LAI', xyz=False) + out = get_flip_matrix('LAI') assert np.all(out == [1, -1, -1, -1, -1, 1, 1, 1]) def test_flip_cs_len4(self): @@ -103,9 +103,9 @@ def test_flip_cs_lap(self): get_flip_matrix('LAP') def test_flip_siemens(self): - out = get_flip_matrix('LAI', xyz=False, manufacturer='Siemens') + out = get_flip_matrix('LAI', manufacturer='Siemens') assert np.all(out == [-1, 1, -1, 1, 1, -1, 1, -1]) def test_flip_ge(self): - out = get_flip_matrix('LAI', xyz=False, manufacturer='GE') + out = get_flip_matrix('LAI', manufacturer='GE') assert np.all(out == [-1, 1, -1, -1, -1, 1, 1, 1])