Skip to content

Commit

Permalink
Split coils and scaner order for optimization + alow specific scaner …
Browse files Browse the repository at this point in the history 
…orders for each (#484)

* First step: split coils and order for optimization

* remove unnecesary prints

* first try at specific orders for rt only

* start of dictionary

* temp commit to compare

* corrections

* Refactoring to dictionary is working

* uniform dictinary keys

* Modified descriptions

* parsed scanner-coil-order

* different coil orders between static and riro implemented

* Dynamic split opt implementation

* test for doc tests

* coil profile modifications to fit with dict

* Tests pass

* fix sequencers test

* more text fixes

* Modification of coeficients output

* coefs fix

* Update shimmingtoolbox/cli/b0shim.py

Co-authored-by: Alex Dastous <47249340+po09i@users.noreply.github.com>

* review corrections

* Test fix + unique coil verification

* GE implementation

* add tests

* Test coil doubles fix

* Test coil doubles fix dynamic

* Review and fix error for different coils than the Prisma

* Add example custom coil config file

* Remove comments

* Resolve comment

---------

Co-authored-by: Alex Dastous <47249340+po09i@users.noreply.github.com>
Co-authored-by: Alexandre D'Astous <po09i@hotmail.com>
  • Loading branch information
@4rnaudB @po09i
3 people committed Jan 9, 2024
1 parent d86c825 commit 072aaa8
Show file tree
Hide file tree
Showing 16 changed files with 759 additions and 386 deletions.
25 changes: 13 additions & 12 deletions config/coil_config.json
View file
@@ -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
}
23 changes: 23 additions & 0 deletions config/custom_coil_config.json
View file
@@ -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"
}
458 changes: 307 additions & 151 deletions shimmingtoolbox/cli/b0shim.py
View file

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion shimmingtoolbox/cli/create_coil_profiles.py
View file
Expand Up @@ -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,
Expand Down
6 changes: 3 additions & 3 deletions shimmingtoolbox/cli/download_data.py
View file
Expand Up @@ -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": (
Expand All @@ -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": (
Expand Down
162 changes: 92 additions & 70 deletions shimmingtoolbox/coils/coil.py
View file
@@ -1,5 +1,6 @@
#!/usr/bin/python3
# -*- coding: utf-8 -*-
import copy
import logging

import numpy as np
Expand Down Expand Up @@ -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.
"""
Expand Down Expand Up @@ -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:
Expand All @@ -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:
Expand All @@ -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.
Expand All @@ -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

0 comments on commit 072aaa8

Please sign in to comment.