diff --git a/docs/source/api_reference/api.rst b/docs/source/api_reference/api.rst index 48c144f07..d051026df 100644 --- a/docs/source/api_reference/api.rst +++ b/docs/source/api_reference/api.rst @@ -56,9 +56,9 @@ Shim .. Error while importing sklearn and skimage .. automodule:: shimmingtoolbox.shim.realtime_shim :members: - -.. automodule:: shimmingtoolbox.optimizer.sequential - :members: +.. Error while importing sklearn (Linear Regression) + .. automodule:: shimmingtoolbox.shim.sequencer + :members: Optimizer --------- diff --git a/shimmingtoolbox/optimizer/basic_optimizer.py b/shimmingtoolbox/optimizer/basic_optimizer.py index bd9de0023..b4114d1dd 100644 --- a/shimmingtoolbox/optimizer/basic_optimizer.py +++ b/shimmingtoolbox/optimizer/basic_optimizer.py @@ -43,20 +43,44 @@ def __init__(self, coils: ListCoil, unshimmed, affine): logging.basicConfig(filename='test_optimizer.log', filemode='w', level=logging.DEBUG) self.coils = coils + self.unshimmed = np.array([]) + self.unshimmed_affine = [] + self.merged_coils = [] + self.merged_bounds = [] + self.set_unshimmed(unshimmed, affine) + def set_unshimmed(self, unshimmed, affine): + """ + Set the unshimmed array to a new array. Resamples coil profiles accordingly. + + Args: + unshimmed (numpy.ndarray): 3d array of unshimmed volume + affine: (numpy.ndarray): 4x4 array containing the qform affine transformation for the unshimmed array + """ # Check dimensions of unshimmed map if unshimmed.ndim != 3: raise ValueError(f"Unshimmed profile has {unshimmed.ndim} dimensions, expected 3 (dim1, dim2, dim3)") - self.unshimmed = unshimmed # Check dimensions of affine if affine.shape != (4, 4): raise ValueError("Shape of affine matrix should be 4x4") + + # Define coil profiles if unshimmed or affine is different than previously + if (self.unshimmed.shape != unshimmed.shape) or not np.all(self.unshimmed_affine == affine): + self.merged_coils, self.merged_bounds = self.merge_coils(unshimmed, affine) + + self.unshimmed = unshimmed self.unshimmed_affine = affine - # Define coil profiles - merged_coils, merged_bounds = self.merge_coils(unshimmed, affine) - self.merged_coils = merged_coils + def set_merged_bounds(self, merged_bounds): + """ + Changes the default bounds set in the coil profile + + Args: + 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)}") self.merged_bounds = merged_bounds def optimize(self, mask): @@ -99,7 +123,6 @@ def merge_coils(self, unshimmed, affine): """ coil_profiles_list = [] - bounds = [] # Define the nibabel unshimmed array nii_unshimmed = nib.Nifti1Image(unshimmed, affine) @@ -109,15 +132,48 @@ def merge_coils(self, unshimmed, affine): # Resample a coil on the unshimmed image resampled_coil = resample_from_to(nii_coil, nii_unshimmed).get_fdata() + coil_profiles_list.append(resampled_coil) + coil_profiles = np.concatenate(coil_profiles_list, axis=3) + + bounds = self.merge_bounds() + + return coil_profiles, bounds + + def merge_bounds(self): + """ + Merge the coil profile bounds into a single array. + + Returns: + list: list of bounds corresponding to each merged coils + """ + + bounds = [] + for coil in self.coils: # Concat coils and bounds - coil_profiles_list.append(resampled_coil) for a_bound in coil.coef_channel_minmax: bounds.append(a_bound) - coil_profiles = np.concatenate(coil_profiles_list, axis=3) + return bounds - return coil_profiles, bounds + def initial_guess_mean_bounds(self): + """ + Calculates the initial guess from the bounds, sets it to the mean of the bounds + + Returns: + np.ndarray: 1d array (n_channels) of coefficient representing the initial guess + + """ + current_0 = [] + for bounds in self.merged_bounds: + avg = np.mean(bounds) + + if np.isnan(avg): + current_0.append(0) + else: + current_0.append(avg) + + return np.array(current_0) def _check_sizing(self, mask): """ diff --git a/shimmingtoolbox/optimizer/lsq_optimizer.py b/shimmingtoolbox/optimizer/lsq_optimizer.py index 63906e081..4bb9112e7 100644 --- a/shimmingtoolbox/optimizer/lsq_optimizer.py +++ b/shimmingtoolbox/optimizer/lsq_optimizer.py @@ -68,7 +68,7 @@ def _apply_sum_constraint(inputs, indexes, coef_sum_max): start_index = end_index # Set up output currents - currents_0 = np.zeros(n_channels) + currents_0 = self.initial_guess_mean_bounds() # Optimize currents_sp = opt.minimize(self._residuals, currents_0, diff --git a/shimmingtoolbox/optimizer/sequential.py b/shimmingtoolbox/optimizer/sequential.py deleted file mode 100644 index 055500d35..000000000 --- a/shimmingtoolbox/optimizer/sequential.py +++ /dev/null @@ -1,89 +0,0 @@ -#!/usr/bin/python3 -# -*- coding: utf-8 -*- - -import numpy as np -from typing import List - -from shimmingtoolbox.optimizer.lsq_optimizer import LsqOptimizer -from shimmingtoolbox.optimizer.basic_optimizer import Optimizer -from shimmingtoolbox.coils.coil import Coil - -ListCoil = List[Coil] - -supported_optimizers = { - 'least_squares': LsqOptimizer, - 'pseudo_inverse': Optimizer -} - - -def sequential_zslice(unshimmed, affine, coils: ListCoil, mask, z_slices, method='least_squares'): - """ - Performs shimming slice by slice using one of the supported optimizers - - Args: - unshimmed (numpy.ndarray): 3D B0 map - affine (np.ndarray): 4x4 array containing the affine transformation for the unshimmed array - coils (ListCoil): List of Coils containing the coil profiles - mask (numpy.ndarray): 3D mask used for the optimizer (only consider voxels with non-zero values). - z_slices (numpy.ndarray): 1D array containing z slices to shim - method (str): Supported optimizer: 'least_squares', 'pseudo_inverse' - Returns: - numpy.ndarray: Coefficients to shim (channels x z_slices.size) - """ - - # Select and initialize the optimizer - optimizer = select_optimizer(method, unshimmed, affine, coils) - - # Optimize slice by slice - currents = optimize_slicewise(optimizer, mask, z_slices) - - return currents - - -def select_optimizer(method, unshimmed, affine, coils: ListCoil): - """ - Select and initialize the optimizer - - Args: - method (str): Supported optimizer: 'least_squares', 'pseudo_inverse' - unshimmed (numpy.ndarray): 3D B0 map - affine (np.ndarray): 4x4 array containing the affine transformation for the unshimmed array - coils (ListCoil): List of Coils containing the coil profiles - - Returns: - Optimizer: Initialized Optimizer object - """ - - # global supported_optimizers - if method in supported_optimizers: - optimizer = supported_optimizers[method](coils, unshimmed, affine) - else: - raise KeyError(f"Method: {method} is not part of the supported optimizers") - - return optimizer - - -def optimize_slicewise(optimizer: Optimizer, mask, z_slices): - """ - Shim slicewise in the specified ROI - - Args: - optimizer (Optimizer): Initialized Optimizer object - mask (numpy.ndarray): 3D mask used for the optimizer (only consider voxels with non-zero values). - z_slices (numpy.ndarray): 1D array containing z slices to shim - - Returns: - numpy.ndarray: Coefficients to shim (channels x z_slices.size) - """ - # Count number of channels - n_channels = optimizer.merged_coils.shape[3] - - z_slices.reshape(z_slices.size) - currents = np.zeros((z_slices.size, n_channels)) - for i in range(z_slices.size): - sliced_mask = np.full_like(mask, fill_value=False) - z = z_slices[i] - sliced_mask[:, :, z:z+1] = mask[:, :, z:z+1] - currents[i, :] = optimizer.optimize(sliced_mask) - - return currents diff --git a/shimmingtoolbox/pmu.py b/shimmingtoolbox/pmu.py index 514e36282..f35e1999e 100644 --- a/shimmingtoolbox/pmu.py +++ b/shimmingtoolbox/pmu.py @@ -37,6 +37,8 @@ def __init__(self, fname_pmu): self.stop_time_mdh = attributes['stop_time_mdh'] self.start_time_mpcu = attributes['start_time_mpcu'] self.stop_time_mpcu = attributes['stop_time_mpcu'] + self.max = attributes['max'] + self.min = attributes['min'] def read_resp(self, fname_pmu): """ @@ -114,7 +116,9 @@ def read_resp(self, fname_pmu): 'start_time_mdh': start_time_mdh, 'stop_time_mdh': stop_time_mdh, 'start_time_mpcu': start_time_mpcu, - 'stop_time_mpcu': stop_time_mpcu + 'stop_time_mpcu': stop_time_mpcu, + 'max': 4095, + 'min': 0 } return attributes diff --git a/shimmingtoolbox/shim/sequencer.py b/shimmingtoolbox/shim/sequencer.py new file mode 100644 index 000000000..9cdb27b64 --- /dev/null +++ b/shimmingtoolbox/shim/sequencer.py @@ -0,0 +1,234 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- + +import numpy as np +from typing import List +from sklearn.linear_model import LinearRegression + +from shimmingtoolbox.optimizer.lsq_optimizer import LsqOptimizer +from shimmingtoolbox.optimizer.basic_optimizer import Optimizer +from shimmingtoolbox.coils.coil import Coil +from shimmingtoolbox.load_nifti import get_acquisition_times +from shimmingtoolbox.pmu import PmuResp + +ListCoil = List[Coil] + +supported_optimizers = { + 'least_squares': LsqOptimizer, + 'pseudo_inverse': Optimizer +} + + +def shim_sequencer(unshimmed, affine, coils: ListCoil, mask, slices, method='least_squares'): + """ + Performs shimming slice by slice using one of the supported optimizers + + Args: + unshimmed (numpy.ndarray): 3D B0 map + affine (np.ndarray): 4x4 array containing the affine transformation for the unshimmed array + coils (ListCoil): List of Coils containing the coil profiles + mask (numpy.ndarray): 3D mask used for the optimizer (only consider voxels with non-zero values). + slices (list): 1D array containing tuples of z slices to shim + method (str): Supported optimizer: 'least_squares', 'pseudo_inverse' + + Returns: + numpy.ndarray: Coefficients to shim (len(slices) x channels) + """ + + # Select and initialize the optimizer + optimizer = select_optimizer(method, unshimmed, affine, coils) + + # Optimize slice by slice + currents = optimize(optimizer, mask, slices) + + return currents + + +def shim_realtime_pmu_sequencer(nii_fieldmap, json_fmap, pmu: PmuResp, coils: ListCoil, static_mask, riro_mask, slices, + opt_method='least_squares'): + """ + Performs realtime shimming using one of the supported optimizers and an external respiratory trace. + + Args: + nii_fieldmap (nibabel.Nifti1Image): Nibabel object containing fieldmap data in 4d where the 4th dimension is the + timeseries. + json_fmap (dict): dict of the json sidecar corresponding to the fieldmap data (Used to find the acquisition + timestamps). + pmu (PmuResp): Filename of the file of the respiratory trace. + coils (ListCoil): List of Coils containing the coil profiles. The coil profiles and the fieldmaps must have + matching units (if fmap is in Hz, the coil profiles must be in hz/unit_shim). + static_mask (numpy.ndarray): 3D mask used for the optimizer to shim the region for the static component. + riro_mask (numpy.ndarray): 3D mask used for the optimizer to shim the region for the riro component. + slices (list): 1D array containing tuples of z slices to shim. + opt_method (str): Supported optimizer: 'least_squares', 'pseudo_inverse'. + + Returns: + (tuple): tuple containing: + + * numpy.ndarray: Static coefficients to shim (len(slices) x channels) [Hz] + * numpy.ndarray: Static coefficients to shim (len(slices) x channels) [Hz/unit_pressure] + * float: Mean pressure of the respiratory trace. + * float: Root mean squared of the pressure. This is provided to compare results between scans, multiply the + riro coefficients by rms of the pressure to do so. + """ + + # Make sure fieldmap has the appropriate dimensions + fieldmap = nii_fieldmap.get_fdata() + affine = nii_fieldmap.affine + if fieldmap.ndim != 4: + raise RuntimeError("Fieldmap must be 4d (x, y, z, t)") + + # Fetch PMU timing + acq_timestamps = get_acquisition_times(nii_fieldmap, json_fmap) + # TODO: deal with saturation + # fit PMU and fieldmap values + acq_pressures = pmu.interp_resp_trace(acq_timestamps) + # # DEBUG: + # acq_pressures = np.array([3000, 2000, 1000, 2000]) + + # regularization --> static, riro + # field(i_vox) = riro(i_vox) * (acq_pressures - mean_p) + static(i_vox) + mean_p = np.mean(acq_pressures) + pressure_rms = np.sqrt(np.mean((acq_pressures - mean_p) ** 2)) + reg = LinearRegression().fit(acq_pressures.reshape(-1, 1) - mean_p, fieldmap.reshape(-1, fieldmap.shape[-1]).T) + + # static/riro contains a 3d matrix of static/riro coefficients + static = reg.intercept_.reshape(fieldmap.shape[:-1]) + riro = reg.coef_.reshape(fieldmap.shape[:-1]) # [unit_shim/unit_pressure], ex: [Hz/unit_pressure] + + # Static shim + optimizer = select_optimizer(opt_method, static, affine, coils) + currents_static = optimize(optimizer, static_mask, slices) + + # Use the currents to define a list of new bounds for the riro optimization + bounds = new_bounds_from_currents(currents_static, optimizer.merged_bounds) + + # Riro shim + # We multiply by the max offset of the siemens pmu [max - min = 4095] so that the bounds take effect on the maximum + # value that the pressure probe can acquire. The equation "riro(i_vox) * (acq_pressures - mean_p)" becomes + # "riro(i_vox) * max_offset" which is the maximum shim we will have. We solve for that to make sure the coils can + # support it. The units of riro * max_offset are: [unit_shim], ex: [Hz] + max_offset = max((pmu.max - pmu.min) - mean_p, mean_p) + + # Set the riro map to shim + optimizer.set_unshimmed(riro * max_offset, affine) + currents_max_riro = optimize(optimizer, riro_mask, slices, shimwise_bounds=bounds) + # Once the currents are solved, we divide by max_offset to return to units of + # [unit_shim/unit_pressure], ex: [Hz/unit_pressure] + currents_riro = currents_max_riro / max_offset + + # Multiplying by the RMS of the pressure allows to make abstraction of the tightness of the bellow + # between scans. This allows to compare results between scans. + # currents_riro_rms = currents_riro * pressure_rms + # [unit_shim/unit_pressure] * rms_pressure, ex: [Hz/unit_pressure] * rms_pressure + + return currents_static, currents_riro, mean_p, pressure_rms + + +def new_bounds_from_currents(currents, old_bounds): + """ + Uses the currents to determine the appropriate bound for a next optimization. It assumes that + "old_current + next_current < old_bound". + + Args: + currents: 2D array (n_shims x n_channels). Direct output from ``optimize``. + old_bounds: 1d list (n_channels) of the merged bounds of the previous optimization. + + Returns: + list: 2d list (n_shims x n_channels) of bounds (min, max) corresponding to each shim and channel. + """ + + new_bounds = [] + for i_shim in range(currents.shape[0]): + shim_bound = [] + for i_channel in range(len(old_bounds)): + a_bound = old_bounds[i_channel] - currents[i_shim, i_channel] + shim_bound.append(tuple(a_bound)) + new_bounds.append(shim_bound) + + return new_bounds + + +def select_optimizer(method, unshimmed, affine, coils: ListCoil): + """ + Select and initialize the optimizer + + Args: + method (str): Supported optimizer: 'least_squares', 'pseudo_inverse' + unshimmed (numpy.ndarray): 3D B0 map + affine (np.ndarray): 4x4 array containing the affine transformation for the unshimmed array + coils (ListCoil): List of Coils containing the coil profiles + + Returns: + Optimizer: Initialized Optimizer object + """ + + # global supported_optimizers + if method in supported_optimizers: + optimizer = supported_optimizers[method](coils, unshimmed, affine) + else: + raise KeyError(f"Method: {method} is not part of the supported optimizers") + + return optimizer + + +def optimize(optimizer: Optimizer, mask, slices, shimwise_bounds=None): + """ + Optimize in the specified ROI according to a specified Optimizer and specified slices + + Args: + optimizer (Optimizer): Initialized Optimizer object + mask (numpy.ndarray): 3D mask used for the optimizer (only consider voxels with non-zero values). + slices (list): 1D array containing tuples of z slices to shim + shimwise_bounds (list): list (n shims) of list (n channels) of tuples (min_coef, max_coef) + + Returns: + numpy.ndarray: Coefficients to shim (len(slices) x channels) + """ + # Count number of channels + n_channels = optimizer.merged_coils.shape[3] + n_shims = len(slices) + currents = np.zeros((n_shims, n_channels)) + for i in range(n_shims): + sliced_mask = np.full_like(mask, fill_value=False) + sliced_mask[:, :, slices[i]] = mask[:, :, slices[i]] + if shimwise_bounds is not None: + optimizer.set_merged_bounds(shimwise_bounds[i]) + currents[i, :] = optimizer.optimize(sliced_mask) + + return currents + + +def define_slices(n_slices: int, factor: int): + """ + Define the slices to shim according to the output covention. + + Args: + n_slices (int): Number of total slices + factor (int): Number of slices per shim + + Returns: + list: 1D list containing tuples of z slices to shim + + Examples: + + :: + + slices = define_slices(10, 2) + print(slices) # [(0, 5), (1, 6), (2, 7), (3, 8), (4, 9)] + + """ + if n_slices <= 0: + return [tuple()] + + slices = [] + n_shims = n_slices // factor + leftover = n_slices % factor + + for i_shim in range(n_shims): + slices.append(tuple(range(i_shim, n_shims * factor, n_shims))) + + if leftover != 0: + slices.append(tuple(range(n_shims * factor, n_slices))) + + return slices diff --git a/test/shim/test_sequencer.py b/test/shim/test_sequencer.py new file mode 100644 index 000000000..ed516a7f0 --- /dev/null +++ b/test/shim/test_sequencer.py @@ -0,0 +1,477 @@ +#!/usr/bin/python3 +# -*- coding: utf-8 -*- + +import numpy as np +import pytest + +from shimmingtoolbox.shim.sequencer import shim_sequencer +from shimmingtoolbox.coils.siemens_basis import siemens_basis +from shimmingtoolbox.simulate.numerical_model import NumericalModel +from shimmingtoolbox.masking.shapes import shapes +from shimmingtoolbox.coils.coil import Coil +from shimmingtoolbox.optimizer.basic_optimizer import Optimizer +from shimmingtoolbox.coils.coordinates import generate_meshgrid + +import os +import nibabel as nib +import json +from matplotlib.figure import Figure + +from shimmingtoolbox.shim.sequencer import shim_realtime_pmu_sequencer +from shimmingtoolbox import __dir_testing__ +from shimmingtoolbox.pmu import PmuResp +from shimmingtoolbox.load_nifti import get_acquisition_times +from shimmingtoolbox.shim.sequencer import define_slices + +DEBUG = False + + +def create_unshimmed(): + # Set up 2-dimensional unshimmed fieldmaps + num_vox = 100 + model_obj = NumericalModel('shepp-logan', num_vox=num_vox) + model_obj.generate_deltaB0('linear', [0.025, 2]) + tr = 0.025 # in s + te = [0.004, 0.008] # in s + model_obj.simulate_measurement(tr, te) + phase_meas1 = model_obj.get_phase() + phase_e1 = phase_meas1[:, :, 0, 0] + phase_e2 = phase_meas1[:, :, 0, 1] + b0_map = ((phase_e2 - phase_e1) / (te[1] - te[0])) / (2 * np.pi) + + # Construct a 3-dimensional synthetic field map by stacking different z-slices along the 3rd dimension. Each + # slice is subjected to a manipulation of model_obj across slices (e.g. rotation, squared) in order to test + # various shim configurations. + unshimmed = np.zeros([num_vox, num_vox, nz]) + for i_n in range(nz // 3): + unshimmed[:, :, 3 * i_n] = b0_map + unshimmed[:, :, (3 * i_n) + 1] = (np.rot90(unshimmed[:, :, 0]) + unshimmed[:, :, 0]) / 2 + unshimmed[:, :, (3 * i_n) + 2] = unshimmed[:, :, 0] ** 2 + + return unshimmed + + +def create_unshimmed_affine(): + # return np.array([[0., 0., 3., 1], + # [-2.91667008, 0., 0., 2], + # [0., 2.91667008, 0., 3], + # [0., 0., 0., 1.]]) + return np.eye(4) + + +def create_constraints(max_coef, min_coef, sum_coef, n_channels=8): + # Set up bounds for output currents + bounds = [] + for _ in range(n_channels): + bounds.append((min_coef, max_coef)) + + constraints = { + "coef_sum_max": sum_coef, + "coef_channel_minmax": bounds + } + return constraints + + +def create_coil(n_x, n_y, n_z, constraints, coil_affine): + # Set up spherical harmonics coil profile + mesh_x, mesh_y, mesh_z = generate_meshgrid((n_x, n_y, n_z), coil_affine) + profiles = siemens_basis(mesh_x, mesh_y, mesh_z) + + # Define coil1 + coil = Coil(profiles[..., :8], coil_affine, constraints) + return coil + + +def create_mask(ref): + # Set up mask + mask = shapes(ref, 'cube', len_dim1=40, len_dim2=40, len_dim3=nz) + return mask + + +nz = 3 # Must be multiple of 3 +a_unshimmed = create_unshimmed() +unshimmed_affine = create_unshimmed_affine() +coil_affine = unshimmed_affine * 2 +coil_affine[3, 3] = 1 +coil1 = create_coil(100, 100, nz, create_constraints(1000, -1000, 2000), unshimmed_affine) +affine = coil_affine * 0.75 +affine[3, 3] = 1 +coil2 = create_coil(150, 120, nz + 10, create_constraints(500, -500, 1500), affine) + + +@pytest.mark.parametrize( + "unshimmed,un_affine,sph_coil,sph_coil2,mask", [( + a_unshimmed, + unshimmed_affine, + coil1, + coil2, + create_mask(a_unshimmed) + )] +) +class TestSequencer(object): + def test_shim_sequencer_lsq(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + # Optimize + z_slices = [] + for i in range(nz): + z_slices.append((i,)) + currents = shim_sequencer(unshimmed, un_affine, [sph_coil], mask, z_slices, method='least_squares') + + assert_results([sph_coil], unshimmed, un_affine, currents, mask, z_slices) + + def test_shim_sequencer_pseudo(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + # Optimize + z_slices = [] + for i in range(nz): + z_slices.append((i,)) + currents = shim_sequencer(unshimmed, un_affine, [sph_coil], mask, z_slices, method='pseudo_inverse') + + assert_results([sph_coil], unshimmed, un_affine, currents, mask, z_slices) + + def test_shim_sequencer_2_coils_lsq(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + + # Optimize + z_slices = [] + for i in range(nz): + z_slices.append((i,)) + currents = shim_sequencer(unshimmed, un_affine, [sph_coil, sph_coil2], mask, z_slices, method='least_squares') + + assert_results([sph_coil, sph_coil2], unshimmed, un_affine, currents, mask, z_slices) + + def test_shim_sequencer_coefs_are_none(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + + coil = create_coil(5, 5, nz, create_constraints(None, None, None), affine) + + # Optimize + z_slices = [] + for i in range(nz): + z_slices.append((i,)) + currents = shim_sequencer(unshimmed, un_affine, [coil], mask, z_slices, method='least_squares') + + assert_results([coil], unshimmed, un_affine, currents, mask, z_slices) + + def test_shim_sequencer_slab_slices(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + """Test for slices arranged as a slab""" + # Optimize + z_slices = [(0, 1, 2)] + currents = shim_sequencer(unshimmed, un_affine, [sph_coil], mask, z_slices, method='least_squares') + + assert_results([sph_coil], unshimmed, un_affine, currents, mask, z_slices) + + def test_shim_sequencer_dynamic_slices(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + """Test for slices arranged for dynamic shimming""" + # Optimize + z_slices = [(0,), (1,), (2,)] + currents = shim_sequencer(unshimmed, un_affine, [sph_coil], mask, z_slices, method='least_squares') + + assert_results([sph_coil], unshimmed, un_affine, currents, mask, z_slices) + + def test_shim_sequencer_multi_slices(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + """Test for slices arranged for multi slice""" + # Optimize + z_slices = [(0, 2), (1,)] + currents = shim_sequencer(unshimmed, un_affine, [sph_coil], mask, z_slices, method='least_squares') + + assert_results([sph_coil], unshimmed, un_affine, currents, mask, z_slices) + + # def test_speed_huge_matrix(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): + # # Create 1 huge coil which essentially is siemens basis concatenated 4 times + # coils = [sph_coil, sph_coil, sph_coil, sph_coil] + # + # coil_profiles_list = [] + # + # for coil in coils: + # # Concat coils and bounds + # coil_profiles_list.append(coil.profile) + # + # coil_profiles = np.concatenate(coil_profiles_list, axis=3) + # constraints = create_constraints(1000, 2000, 32) + # + # huge_coil = Coil(coil_profiles, sph_coil.affine, constraints) + # + # z_slices = np.array(range(unshimmed.shape[2])) + # currents = sequential_zslice(unshimmed, un_affine, [huge_coil], mask, z_slices) + # + # # Calculate theoretical shimmed map + # opt = Optimizer([huge_coil], unshimmed, un_affine) + # shimmed = unshimmed + np.sum(currents * opt.merged_coils, axis=3, keepdims=False) + # + # for i_slice in z_slices: + # sum_shimmed = np.sum(np.abs(mask[:, :, i_slice] * shimmed[:, :, i_slice])) + # sum_unshimmed = np.sum(np.abs(mask[:, :, i_slice] * unshimmed[:, :, i_slice])) + # # print(f"\nshimmed: {sum_shimmed}, unshimmed: {sum_unshimmed}, current: \n{currents[i_slice, :]}") + # assert sum_shimmed <= sum_unshimmed + + +def assert_results(coil, unshimmed, un_affine, currents, mask, z_slices): + # Calculate theoretical shimmed map + opt = Optimizer(coil, unshimmed, un_affine) + + if DEBUG: + # Save fieldmap + fname_fieldmap_2 = os.path.join(os.curdir, 'fig_fieldmap.nii.gz') + nii_fieldmap = nib.Nifti1Image(unshimmed, un_affine) + nib.save(nii_fieldmap, fname_fieldmap_2) + + # Save coil profiles as nifti + fname_coil = os.path.join(os.curdir, 'fig_coil_orig.nii.gz') + nii_coil = nib.Nifti1Image(coil[0].profile, coil[0].affine) + nib.save(nii_coil, fname_coil) + + # save resampled coil profiles + fname_coil_res = os.path.join(os.curdir, 'fig_coil_resampled.nii.gz') + nii_coil = nib.Nifti1Image(opt.merged_coils, opt.unshimmed_affine) + nib.save(nii_coil, fname_coil_res) + + correction_per_channel = np.zeros_like(opt.merged_coils) + shimmed = np.zeros_like(unshimmed) + for i_shim in range(len(z_slices)): + correction_per_channel[..., z_slices[i_shim], :] = (currents[i_shim] * opt.merged_coils)[:, :, z_slices[i_shim], :] + correction = np.sum(correction_per_channel[..., z_slices[i_shim], :], axis=3, keepdims=False) + shimmed[..., z_slices[i_shim]] = unshimmed[..., z_slices[i_shim]] + correction + + sum_shimmed = np.sum(np.abs(mask[:, :, z_slices[i_shim]] * shimmed[:, :, z_slices[i_shim]])) + sum_unshimmed = np.sum(np.abs(mask[:, :, z_slices[i_shim]] * unshimmed[:, :, z_slices[i_shim]])) + print(f"\nshimmed: {sum_shimmed}, unshimmed: {sum_unshimmed}, current: \n{currents[i_shim, :]}") + # assert sum_shimmed < sum_unshimmed + + if DEBUG: + # save resampled coil profiles + fname_correction = os.path.join(os.curdir, 'fig_correction.nii.gz') + nii_correction = nib.Nifti1Image(correction_per_channel, opt.unshimmed_affine) + nib.save(nii_correction, fname_correction) + + +def test_realtime_sequencer(): + # Fieldmap + fname_fieldmap = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + 'sub-example_fieldmap.nii.gz') + nii_fieldmap = nib.load(fname_fieldmap) + unshimmed = nii_fieldmap.get_fdata() + # Create affine closer to isocenter + new_affine = np.eye(4) + new_affine[:3, :3] = nii_fieldmap.affine[:3, :3] + nii_fieldmap = nib.Nifti1Image(unshimmed, new_affine, header=nii_fieldmap.header) + unshimmed = nii_fieldmap.get_fdata() + + # TODO: add linear riro, riro could require 0th order shim, is this done through frequency adjust? No, we should do + # frequency adjust, for now, single slice can be 0th order shim using through slice "gradient" + # fake[..., 0] contains the original linear fieldmap. This repeats the linear fieldmap over the 3rd dim and scale + # down + # Dont forget to change the pmu trace + # fake = create_unshimmed() + # # fake_temp = np.zeros([100, 100, 3, 4]) + # # for i_temp in range(4): + # # fake_temp[..., i_temp] = fake + # fake_temp = np.zeros([100, 100, 3, 4]) + # lin = np.repeat(fake[:, :, 0, np.newaxis], 3, axis=2) / 10 + # fake_temp[..., 0] = fake + lin + # fake_temp[..., 1] = fake + # fake_temp[..., 2] = fake - lin + # fake_temp[..., 3] = fake + # nii_fieldmap = nib.Nifti1Image(fake_temp, create_unshimmed_affine(), header=nii_fieldmap.header) + # unshimmed = nii_fieldmap.get_fdata() + + # Set up mask + # static + nx, ny, nz, _ = nii_fieldmap.shape + static_mask = shapes(unshimmed[..., 0], 'cube', + center_dim1=int(nx / 2) - 6, + center_dim2=int(ny / 2) - 5, + len_dim1=5, len_dim2=10, len_dim3=nz) + + # Riro + riro_mask = shapes(unshimmed[..., 0], 'cube', + center_dim1=int(nx / 2) - 6, + center_dim2=int(ny / 2) - 5, + len_dim1=5, len_dim2=10, len_dim3=nz) + + # Pmu + fname_resp = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'PMUresp_signal.resp') + pmu = PmuResp(fname_resp) + + # Path for json file + fname_json = os.path.join(__dir_testing__, 'realtime_zshimming_data', 'nifti', 'sub-example', 'fmap', + 'sub-example_magnitude1.json') + with open(fname_json) as json_file: + json_data = json.load(json_file) + + # Create Coil + coil_affine = nii_fieldmap.affine + coil = create_coil(150, 150, nz + 10, create_constraints(np.inf, -np.inf, np.inf), coil_affine) + + slices = define_slices(unshimmed.shape[2], 1) + + currents_static, currents_riro, mean_p, p_rms = shim_realtime_pmu_sequencer(nii_fieldmap, json_data, pmu, [coil], + static_mask, riro_mask, slices, + opt_method='least_squares') + # currents_static, currents_riro, mean_p, p_rms = shim_realtime_pmu_sequencer(nii_fieldmap, json_data, pmu, [coil], + # static_mask, riro_mask, slices, + # opt_method='pseudo_inverse') + currents_riro_rms = currents_riro * p_rms + + print(f"\nSlices: {slices}" + f"\nFieldmap affine:\n{nii_fieldmap.affine}\n" + f"Coil affine:\n{coil_affine}\n" + f"Static currents:\n{currents_static}\n" + f"Riro currents * p_rms:\n{currents_riro_rms}\n") + + # Calculate theoretical shimmed map + # Calc pressure + acq_timestamps = get_acquisition_times(nii_fieldmap, json_data) + acq_pressures = pmu.interp_resp_trace(acq_timestamps) + # # DEBUG + # acq_pressures = [3000, 2000, 1000, 2000] + + # shim + opt = Optimizer([coil], unshimmed[..., 0], nii_fieldmap.affine) + shimmed_static_riro = np.zeros_like(unshimmed) + shimmed_static = np.zeros_like(unshimmed) + shimmed_riro = np.zeros_like(unshimmed) + masked_shim_static_riro = np.zeros_like(unshimmed) + masked_shim_static = np.zeros_like(unshimmed) + masked_shim_riro = np.zeros_like(unshimmed) + masked_unshimmed = np.zeros_like(unshimmed) + shim_trace_static_riro = [] + shim_trace_static = [] + shim_trace_riro = [] + unshimmed_trace = [] + for i_shim in range(len(slices)): + # Calculate static correction + correction_static = np.sum(currents_static[i_shim] * + opt.merged_coils, axis=3, keepdims=False)[..., slices[i_shim]] + + riro_profile = np.sum(currents_riro[i_shim] * opt.merged_coils, axis=3, keepdims=False)[..., slices[i_shim]] + for i_t in range(nii_fieldmap.shape[3]): + correction_riro = riro_profile * (acq_pressures[i_t] - mean_p) + shimmed_static[..., slices[i_shim], i_t] = unshimmed[..., slices[i_shim], i_t] + correction_static + shimmed_static_riro[..., slices[i_shim], i_t] = shimmed_static[..., slices[i_shim], i_t] + correction_riro + shimmed_riro[..., slices[i_shim], i_t] = unshimmed[..., slices[i_shim], i_t] + correction_riro + + sum_shimmed_static = np.sum(np.abs(static_mask[:, :, slices[i_shim]] * shimmed_static[:, :, slices[i_shim], i_t])) + sum_shimmed_static_riro = np.sum(np.abs(riro_mask[:, :, slices[i_shim]] * shimmed_static_riro[:, :, slices[i_shim], i_t])) + sum_shimmed_riro = np.sum(np.abs(riro_mask[:, :, slices[i_shim]] * shimmed_riro[:, :, slices[i_shim], i_t])) + sum_unshimmed = np.sum(np.abs(riro_mask[:, :, slices[i_shim]] * unshimmed[:, :, slices[i_shim], i_t])) + print(f"\ni_shim: {i_shim}, t: {i_t}" + f"\nshimmed: {sum_shimmed_static_riro}, unshimmed: {sum_unshimmed}, " + f"Static currents:\n{currents_static[i_shim]}\n" + f"Riro currents:\n{currents_riro[i_shim] * (acq_pressures[i_t] - mean_p)}\n") + + shim_trace_static.append(sum_shimmed_static) + shim_trace_static_riro.append(sum_shimmed_static_riro) + shim_trace_riro.append(sum_shimmed_riro) + unshimmed_trace.append(sum_unshimmed) + + assert sum_shimmed_static_riro < sum_unshimmed + + if DEBUG: + n_t = nii_fieldmap.shape[3] + for i_t in range(n_t): + # static mask and riro are the same + masked_shim_static[..., i_t] = static_mask * shimmed_static[..., i_t] + masked_shim_static_riro[..., i_t] = riro_mask * shimmed_static_riro[..., i_t] + masked_shim_riro[..., i_t] = shimmed_riro[..., i_t] + masked_unshimmed[..., i_t] = riro_mask * unshimmed[..., i_t] + + # Plot Static and RIRO + i_t = 0 + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(2, 3, 1) + im = ax.imshow(np.rot90(masked_shim_static_riro[..., 0, i_t])) + fig.colorbar(im) + ax.set_title("masked_shim static + riro") + ax = fig.add_subplot(2, 3, 2) + im = ax.imshow(np.rot90(masked_shim_static[..., 0, i_t])) + fig.colorbar(im) + ax.set_title("masked_shim static") + ax = fig.add_subplot(2, 3, 3) + im = ax.imshow(np.rot90(masked_unshimmed[..., 0, i_t])) + fig.colorbar(im) + ax.set_title("masked_unshimmed") + + ax = fig.add_subplot(2, 3, 4) + im = ax.imshow(np.rot90(shimmed_static_riro[..., 0, i_t])) + fig.colorbar(im) + ax.set_title("shim static + riro") + ax = fig.add_subplot(2, 3, 5) + im = ax.imshow(np.rot90(shimmed_static[..., 0, i_t])) + fig.colorbar(im) + ax.set_title("shim static") + ax = fig.add_subplot(2, 3, 6) + im = ax.imshow(np.rot90(unshimmed[..., 0, i_t])) + fig.colorbar(im) + ax.set_title("unshimmed") + fname_figure = os.path.join(os.curdir, 'fig_realtime_masked_shimmed_vs_unshimmed.png') + fig.savefig(fname_figure) + + # plot shimmed and unshimmed trace + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(111) + ax.plot(shim_trace_static_riro, label='shimmed static + riro') + ax.plot(shim_trace_static, label='shimmed static') + ax.plot(shim_trace_riro, label='shimmed_riro') + ax.plot(unshimmed_trace, label='unshimmed') + ax.legend() + ax.set_ylim(0, max(unshimmed_trace)) + ax.set_title("Unshimmed vs shimmed values") + fname_figure = os.path.join(os.curdir, 'fig_trace_shimmed_vs_unshimmed.png') + fig.savefig(fname_figure) + + fig = Figure(figsize=(10, 10)) + ax = fig.add_subplot(111) + ax.plot(acq_pressures, label='pressures') + ax.legend() + ax.set_ylim(0, 4095) + ax.set_title("Pressures vs time points") + fname_figure = os.path.join(os.curdir, 'fig_trace_pressures.png') + fig.savefig(fname_figure) + + # Save fieldmap + fname_fieldmap_2 = os.path.join(os.curdir, 'fig_fieldmap.nii.gz') + nib.save(nii_fieldmap, fname_fieldmap_2) + + # Save coil profiles as nifti + fname_coil = os.path.join(os.curdir, 'fig_coil_orig.nii.gz') + nii_coil = nib.Nifti1Image(coil.profile, coil.affine) + nib.save(nii_coil, fname_coil) + + # save resampled coil profiles + fname_coil_res = os.path.join(os.curdir, 'fig_coil_resampled.nii.gz') + nii_coil = nib.Nifti1Image(opt.merged_coils, opt.unshimmed_affine) + nib.save(nii_coil, fname_coil_res) + + # metric to isolate temporal and static component + # Temporal: Compute the STD across time pixelwise, and then compute the mean across pixels. + # Static: Compute the MEAN across time pixelwise, and then compute the STD across pixels. + + rep_mask = np.repeat(static_mask[..., np.newaxis], nii_fieldmap.shape[3], 3) + ma_unshimmed = np.ma.array(unshimmed, mask=rep_mask == False) + ma_shim_static = np.ma.array(shimmed_static, mask=rep_mask == False) + ma_shim_static_riro = np.ma.array(shimmed_static_riro, mask=rep_mask == False) + ma_shim_riro = np.ma.array(shimmed_riro, mask=rep_mask == False) + + # Temporal + temp_shim_static = np.ma.mean(np.ma.std(ma_shim_static, 3)) + temp_shim_static_riro = np.ma.mean(np.ma.std(ma_shim_static_riro, 3)) + temp_shim_riro = np.ma.mean(np.ma.std(ma_shim_riro, 3)) + temp_unshimmed = np.ma.mean(np.ma.std(ma_unshimmed, 3)) + + # Static + static_shim_static = np.ma.std(np.ma.mean(ma_shim_static, 3)) + static_shim_static_riro = np.ma.std(np.ma.mean(ma_shim_static_riro, 3)) + static_shim_riro = np.ma.std(np.ma.mean(ma_shim_riro, 3)) + static_unshimmed = np.ma.std(np.ma.mean(ma_unshimmed, 3)) + print(f"\nTemporal: Compute the STD across time pixelwise, and then compute the mean across pixels." + f"\ntemp_shim_static: {temp_shim_static}" + f"\ntemp_shim_static_riro: {temp_shim_static_riro}" + f"\ntemp_shim_riro: {temp_shim_riro}" + f"\ntemp_unshimmed: {temp_unshimmed}" + f"\nStatic: Compute the MEAN across time pixelwise, and then compute the STD across pixels." + f"\nstatic_shim_static: {static_shim_static}" + f"\nstatic_shim_static_riro: {static_shim_static_riro}" + f"\nstatic_shim_riro: {static_shim_riro}" + f"\nstatic_unshimmed: {static_unshimmed}") + + # Todo: + # Use same affine for coil and for defining siemens basis + # + # Look at input of tests if the coil profiles are at the same place as the unshimmed map diff --git a/test/test_sequential_zslice.py b/test/test_sequential_zslice.py deleted file mode 100644 index c31f26a9c..000000000 --- a/test/test_sequential_zslice.py +++ /dev/null @@ -1,170 +0,0 @@ -#!/usr/bin/python3 -# -*- coding: utf-8 -*- - -import numpy as np -import pytest - -from shimmingtoolbox.optimizer.sequential import sequential_zslice -from shimmingtoolbox.coils.siemens_basis import siemens_basis -from shimmingtoolbox.simulate.numerical_model import NumericalModel -from shimmingtoolbox.masking.shapes import shapes -from shimmingtoolbox.coils.coil import Coil -from shimmingtoolbox.optimizer.basic_optimizer import Optimizer -from shimmingtoolbox.coils.coordinates import generate_meshgrid - - -def create_unshimmed(): - # Set up 2-dimensional unshimmed fieldmaps - num_vox = 100 - model_obj = NumericalModel('shepp-logan', num_vox=num_vox) - model_obj.generate_deltaB0('linear', [0.0, 20]) - tr = 0.025 # in s - te = [0.004, 0.008] # in s - model_obj.simulate_measurement(tr, te) - phase_meas1 = model_obj.get_phase() - phase_e1 = phase_meas1[:, :, 0, 0] - phase_e2 = phase_meas1[:, :, 0, 1] - b0_map = ((phase_e2 - phase_e1) / (te[1] - te[0])) / (2 * np.pi) - - # Construct a 3-dimensional synthetic field map by stacking different z-slices along the 3rd dimension. Each - # slice is subjected to a manipulation of model_obj across slices (e.g. rotation, squared) in order to test - # various shim configurations. - unshimmed = np.zeros([num_vox, num_vox, nz]) - for i_n in range(nz // 3): - unshimmed[:, :, 3 * i_n] = b0_map - unshimmed[:, :, (3 * i_n) + 1] = (np.rot90(unshimmed[:, :, 0]) + unshimmed[:, :, 0]) / 2 - unshimmed[:, :, (3 * i_n) + 2] = unshimmed[:, :, 0] ** 2 - - return unshimmed - - -def create_unshimmed_affine(): - return np.array([[0., 0., 3., 1], - [-2.91667008, 0., 0., 2], - [0., 2.91667008, 0., 3], - [0., 0., 0., 1.]]) - - -def create_constraints(max_coef, min_coef, sum_coef, n_channels=8): - # Set up bounds for output currents - bounds = [] - for _ in range(n_channels): - bounds.append((min_coef, max_coef)) - - constraints = { - "coef_sum_max": sum_coef, - "coef_channel_minmax": bounds - } - return constraints - - -def create_coil(n_x, n_y, n_z, constraints, coil_affine): - # Set up spherical harmonics coil profile - mesh_x, mesh_y, mesh_z = generate_meshgrid((n_x, n_y, n_z), coil_affine) - profiles = siemens_basis(mesh_x, mesh_y, mesh_z) - - # Define coil1 - coil = Coil(profiles, coil_affine, constraints) - return coil - - -def create_mask(ref): - # Set up mask - mask = shapes(ref, 'cube', len_dim1=40, len_dim2=40, len_dim3=nz) - return mask - - -nz = 3 # Must be multiple of 3 -a_unshimmed = create_unshimmed() -affine = np.eye(4) * 2 -affine[3, 3] = 1 -coil1 = create_coil(150, 150, nz + 10, create_constraints(1000, -1000, 2000), affine) -affine = np.eye(4) * 0.75 -affine[3, 3] = 1 -coil2 = create_coil(150, 120, nz + 10, create_constraints(500, -500, 1500), affine) - - -@pytest.mark.parametrize( - "unshimmed,un_affine,sph_coil,sph_coil2,mask", [( - a_unshimmed, - create_unshimmed_affine(), - coil1, - coil2, - create_mask(a_unshimmed) - ) - ] -) -class TestSequentialZSlice(object): - def test_zslice_lsq(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): - # Optimize - z_slices = np.array(range(unshimmed.shape[2])) - currents = sequential_zslice(unshimmed, un_affine, [sph_coil], mask, z_slices, method='least_squares') - - assert_results([sph_coil], unshimmed, un_affine, currents, mask, z_slices) - - def test_zslice_pseudo(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): - # Optimize - z_slices = np.array(range(unshimmed.shape[2])) - currents = sequential_zslice(unshimmed, un_affine, [sph_coil], mask, z_slices, - method='pseudo_inverse') - - assert_results([sph_coil], unshimmed, un_affine, currents, mask, z_slices) - - def test_zslice_2_coils_lsq(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): - - # Optimize - z_slices = np.array(range(unshimmed.shape[2])) - currents = sequential_zslice(unshimmed, un_affine, [sph_coil, sph_coil2], mask, - z_slices, method='least_squares') - - assert_results([sph_coil, sph_coil2], unshimmed, un_affine, currents, mask, z_slices) - - def test_coefs_are_none(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): - - coil = create_coil(5, 5, nz, create_constraints(None, None, None), affine) - - # Optimize - z_slices = np.array(range(unshimmed.shape[2])) - currents = sequential_zslice(unshimmed, un_affine, [coil], mask, z_slices, method='least_squares') - - assert_results([coil], unshimmed, un_affine, currents, mask, z_slices) - - # def test_speed_huge_matrix(self, unshimmed, un_affine, sph_coil, sph_coil2, mask): - # # Create 1 huge coil which essentially is siemens basis concatenated 4 times - # coils = [sph_coil, sph_coil, sph_coil, sph_coil] - # - # coil_profiles_list = [] - # - # for coil in coils: - # # Concat coils and bounds - # coil_profiles_list.append(coil.profile) - # - # coil_profiles = np.concatenate(coil_profiles_list, axis=3) - # constraints = create_constraints(1000, 2000, 32) - # - # huge_coil = Coil(coil_profiles, sph_coil.affine, constraints) - # - # z_slices = np.array(range(unshimmed.shape[2])) - # currents = sequential_zslice(unshimmed, un_affine, [huge_coil], mask, z_slices) - # - # # Calculate theoretical shimmed map - # opt = Optimizer([huge_coil], unshimmed, un_affine) - # shimmed = unshimmed + np.sum(currents * opt.merged_coils, axis=3, keepdims=False) - # - # for i_slice in z_slices: - # sum_shimmed = np.sum(np.abs(mask[:, :, i_slice] * shimmed[:, :, i_slice])) - # sum_unshimmed = np.sum(np.abs(mask[:, :, i_slice] * unshimmed[:, :, i_slice])) - # # print(f"\nshimmed: {sum_shimmed}, unshimmed: {sum_unshimmed}, current: \n{currents[i_slice, :]}") - # assert sum_shimmed <= sum_unshimmed - - -def assert_results(coil, unshimmed, un_affine, currents, mask, z_slices): - # Calculate theoretical shimmed map - opt = Optimizer(coil, unshimmed, un_affine) - shimmed = unshimmed + np.sum(currents * opt.merged_coils, axis=3, keepdims=False) - - for i_slice in z_slices: - sum_shimmed = np.sum(np.abs(mask[:, :, i_slice] * shimmed[:, :, i_slice])) - sum_unshimmed = np.sum(np.abs(mask[:, :, i_slice] * unshimmed[:, :, i_slice])) - print(f"\nshimmed: {sum_shimmed}, unshimmed: {sum_unshimmed}, current: \n{currents[i_slice, :]}") - assert sum_shimmed < sum_unshimmed