Skip to content

Commit

Permalink
Merge 7ad62db into 6170cf3
Browse files Browse the repository at this point in the history
  • Loading branch information
@po09i
po09i committed Jan 24, 2024
2 parents 6170cf3 + 7ad62db commit 60797d1
Show file tree
Hide file tree
Showing 8 changed files with 441 additions and 426 deletions.
23 changes: 13 additions & 10 deletions shimmingtoolbox/cli/b0shim.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +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]
AVAILABLE_ORDERS = [-1, 0, 1, 2, 3]


@click.group(context_settings=CONTEXT_SETTINGS,
Expand Down Expand Up @@ -139,11 +139,12 @@ 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)

# Parse scanner_coil_order
scanner_coil_order = parse_orders(scanner_coil_order)

# Load the fieldmap
nii_fmap_orig = nib.load(fname_fmap)

Expand Down Expand Up @@ -1006,7 +1007,6 @@ def parse_orders(orders: str):

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:
Expand Down Expand Up @@ -1046,9 +1046,9 @@ def _load_coils(coils, orders, fname_constraints, nii_fmap, scanner_shim_setting
for key in orders_to_delete:
del sph_contraints['coef_channel_minmax'][key]
else:
sph_contraints = get_scanner_constraints(manufacturers_model_name, orders)
sph_contraints = get_scanner_constraints(manufacturers_model_name, orders, manufacturer)

sph_contraints_calc = calculate_scanner_constraints(sph_contraints, scanner_shim_settings, orders)
sph_contraints_calc = calculate_scanner_constraints(sph_contraints, scanner_shim_settings, orders, manufacturer)
scanner_coil = ScannerCoil(nii_fmap.shape[:3], nii_fmap.affine, sph_contraints_calc, orders,
manufacturer=manufacturer)
list_coils.append(scanner_coil)
Expand All @@ -1060,8 +1060,7 @@ def _load_coils(coils, orders, fname_constraints, nii_fmap, scanner_shim_setting
return list_coils


def calculate_scanner_constraints(constraints: dict, scanner_shim_settings, orders):
# ! Modify description if everything works
def calculate_scanner_constraints(constraints: dict, scanner_shim_settings, orders, manufacturer):
""" Calculate the constraints that should be used for the scanner by considering the current shim settings and the
absolute bounds
Expand All @@ -1070,7 +1069,6 @@ def calculate_scanner_constraints(constraints: dict, scanner_shim_settings, orde
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
"""
Expand Down Expand Up @@ -1098,7 +1096,10 @@ def _initial_in_bounds(coefs: dict, bounds: dict):
# Set the initial coefficients to 0
initial_coefs = {}
for order in orders:
initial_coefs[str(order)] = [0] * (order * 2 + 1)
if order == 3 and manufacturer == 'Siemens':
initial_coefs['3'] = [0] * 4
else:
initial_coefs[str(order)] = [0] * (order * 2 + 1)
if initial_coefs == {}:
initial_coefs = None

Expand All @@ -1113,6 +1114,8 @@ def _initial_in_bounds(coefs: dict, bounds: dict):
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']
if scanner_shim_settings['3'] is not None and 3 in orders:
initial_coefs['3'] = scanner_shim_settings['3']

# Make sure the initial coefficients are within the specified bounds
_initial_in_bounds(initial_coefs, constraints['coef_channel_minmax'])
Expand Down
107 changes: 45 additions & 62 deletions shimmingtoolbox/coils/coil.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
import numpy as np
from typing import Tuple

from shimmingtoolbox.coils.spher_harm_basis import siemens_basis, ge_basis, philips_basis, SHIM_CS
from shimmingtoolbox.coils.spher_harm_basis import siemens_basis, ge_basis, philips_basis, SHIM_CS, channels_per_order
from shimmingtoolbox.coils.coordinates import generate_meshgrid

logger = logging.getLogger(__name__)
Expand Down Expand Up @@ -99,13 +99,13 @@ def load_constraints(self, constraints: dict):
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)
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])
[-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)
[constraints["coef_channel_minmax"][key][i][0], np.inf]

if key_name == "coef_sum_max":
if constraints["coef_sum_max"] is None:
Expand Down Expand Up @@ -142,42 +142,26 @@ def __init__(self, dim_volume, affine, constraints, orders, manufacturer=""):
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.orders)

super().__init__(sph_coil_profile, affine, constraints)

def _create_coil_profile(self, dim, manufacturer=None):
# Define profile for Tx (constant volume)
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]
# Create spherical harmonics coil profiles
# f0, orders
mesh1, mesh2, mesh3 = generate_meshgrid(dim, self.affine)
if manufacturer == 'SIEMENS':
sph_coil_profile = siemens_basis(mesh1, mesh2, mesh3, orders=tuple(self.orders),
shim_cs=self.coord_system)
elif manufacturer == 'GE':
sph_coil_profile = ge_basis(mesh1, mesh2, mesh3, orders=tuple(self.orders),
shim_cs=self.coord_system)
elif manufacturer == 'PHILIPS':
sph_coil_profile = philips_basis(mesh1, mesh2, mesh3, orders=tuple(self.orders),
shim_cs=self.coord_system)
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(temp_orders),
shim_cs=self.coord_system)
elif manufacturer == 'GE':
profile_orders = ge_basis(mesh1, mesh2, mesh3, orders=tuple(temp_orders),
shim_cs=self.coord_system)
elif manufacturer == 'PHILIPS':
profile_orders = philips_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(temp_orders),
shim_cs=self.coord_system)
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)
logger.warning(f"{manufacturer} manufacturer not implemented. Outputting in Hz, uT/m, uT/m^2, uT/m^3 for "
"order 0, 1, 2 and 3 respectively")
sph_coil_profile = siemens_basis(mesh1, mesh2, mesh3, orders=tuple(self.orders),
shim_cs=self.coord_system)

return sph_coil_profile

Expand All @@ -188,24 +172,24 @@ def __eq__(self, __value: object) -> bool:
return super().__eq__(__value)


def get_scanner_constraints(manufacturers_model_name, orders):
def get_scanner_constraints(manufacturers_model_name, orders, manufacturer):
""" 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
orders (list): List of all orders of the shim system to be used
manufacturer (str): Manufacturer of the scanner
Returns:
dict: The constraints including the scanner name, bounds and the maximum sum of currents.
"""

constraints = {
"coef_channel_minmax": {"0": [], "1": [], "2": [], "3": []},
"coef_sum_max": None
}
if manufacturers_model_name == "Prisma_fit":
constraints = {
"name": "Prisma_fit",
"coef_channel_minmax": {"0": [], "1": [], "2": []},
"coef_sum_max": None
}
constraints["name"] = "Prisma_fit"
if 0 in orders:
constraints["coef_channel_minmax"]["0"].append([123100100, 123265000])
if 1 in orders:
Expand All @@ -217,16 +201,17 @@ def get_scanner_constraints(manufacturers_model_name, orders):
[-3503.299, 3503.299],
[-3551.29, 3551.29],
[-3487.302, 3487.302]])
if 3 in orders:
n_channels = channels_per_order(3, manufacturer)
logger.warning(f"3rd order not available on the {manufacturers_model_name}, using the typical {n_channels} "
"unbounded 3rd order shim terms")
constraints["coef_channel_minmax"]["3"] = [[None, None] for _ in range(n_channels)]

elif manufacturers_model_name == "Investigational_Device_7T":
constraints = {
"name": "Investigational_Device_7T",
"coef_channel_minmax": {"0": [], "1": [], "2": []},
"coef_sum_max": None
}
constraints["name"] = "Investigational_Device_7T"
if 0 in orders:
pass
# todo: f0 min and max is wrong
# todo: f0 min and max are wrong
constraints["coef_channel_minmax"]["0"].append([None, None])
if 1 in orders:
for _ in range(3):
Expand All @@ -237,27 +222,22 @@ def get_scanner_constraints(manufacturers_model_name, orders):
[-791.84, 791.84],
[-615.87, 615.87],
[-615.87, 615.87]])
if 3 in orders:
# Todo: Get 3rd order shim constraints
constraints["coef_channel_minmax"]["3"] = [[None, None] for _ in range(channels_per_order(3, manufacturer))]

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 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)]
constraints["name"] = "Unknown"
for order in orders:
constraints["coef_channel_minmax"][str(order)] = [[None, None] for _ in range(
channels_per_order(order, manufacturer))]

return constraints


def restrict_sph_constraints(bounds: dict, orders):
# ! Modify description if everything works
""" Select bounds according to the order specified
Args:
Expand All @@ -278,6 +258,9 @@ def restrict_sph_constraints(bounds: dict, orders):
if 2 in orders:
# f0, ch1, ch2, ch3, ch4, ch5, ch6, ch7, ch8 -- > [9]
minmax_out["2"] = bounds["2"]
if 3 in orders:
minmax_out["3"] = bounds["3"]

if minmax_out == {}:
raise NotImplementedError("Order must be between 0 and 2")

Expand Down
Loading

0 comments on commit 60797d1

Please sign in to comment.