-
Notifications
You must be signed in to change notification settings - Fork 11
/
coil.py
262 lines (213 loc) · 10.7 KB
/
coil.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
#!/usr/bin/python3
# -*- coding: utf-8 -*-
import logging
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.coordinates import generate_meshgrid
logger = logging.getLogger(__name__)
required_constraints = [
"name",
"coef_channel_minmax",
"coef_sum_max"
]
class Coil(object):
"""
Coil profile object that stores coil profiles and there constraints
Attributes:
dim (Tuple[int]): Dimension along specific axis. dim: 0,1,2 are spatial axes, while dim: 3 corresponds to the
coil channel.
profile (np.ndarray): (dim1, dim2, dim3, channels) 4d array of N 3d coil profiles
affine (np.ndarray): 4x4 array containing the affine transformation associated with the NIfTI file of the coil
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
used to specify no bounds.
name (str): Name of the coil.
"""
def __init__(self, profile, affine, constraints):
""" Initialize Coil
Args:
profile (np.ndarray): Coil profile (dim1, dim2, dim3, channels) 4d array of N 3d coil profiles
affine (np.ndarray): 4x4 array containing the qform affine transformation for the coil profiles
constraints (dict): dict containing the constraints for the coil profiles. Required keys:
* name (str): Name of the coil.
* coef_sum_max (float): Contains the maximum value for the sum of the coefficients. None is used to
specify no bounds
* coef_channel_max (list): List of ``(min, max)`` pairs for each coil channels. (None, None) is
used to specify no bounds.
Examples:
::
# Example of constraints
constraints = {
'name': "dummy coil",
'coef_sum_max': 10,
# 8 channel coil
'coef_channel_minmax': [(-2, 2), (-2, 2), (-2, 2), (-2, 2), (-3, 3), (-3, 3), (-3, 3), (-3, 3)],
}
"""
self.dim = (np.nan,) * 4
self.profile = profile
self.required_constraints = required_constraints
if affine.shape != (4, 4):
raise ValueError("Shape of affine matrix should be 4x4")
self.affine = affine
self.name = ""
self.coef_channel_minmax = self.coef_sum_max = -1
self.load_constraints(constraints)
@property
def profile(self):
return self._profile
@profile.setter
def profile(self, profile):
if profile.ndim != 4:
raise ValueError(f"Coil profile has {profile.ndim} dimensions, expected 4 (dim1, dim2, dim3, channel)")
self.dim = profile.shape
self._profile = profile
def load_constraints(self, constraints):
"""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]:
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)
if key_name == "coef_sum_max":
if constraints["coef_sum_max"] is None:
constraints["coef_sum_max"] = np.inf
setattr(self, key_name, constraints[key_name])
else:
raise KeyError(f"Missing required constraint: {key_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
manufacturer = manufacturer.upper()
if manufacturer in SHIM_CS:
self.coord_system = SHIM_CS[manufacturer.upper()]
else:
logger.warning(f"Unknown manufacturer {manufacturer}, assuming RAS")
self.coord_system = 'RAS'
self.affine = affine
# 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)
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:
# f0 --> [1]
sph_coil_profile = profile_order_0[..., np.newaxis]
else:
# f0, orders
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)),
shim_cs=self.coord_system)
elif manufacturer == 'GE':
profile_orders = ge_basis(mesh1, mesh2, mesh3, orders=tuple(range(1, self.order + 1)),
shim_cs=self.coord_system)
elif manufacturer == 'PHILIPS':
profile_orders = philips_basis(mesh1, mesh2, mesh3, orders=tuple(range(1, self.order + 1)),
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)),
shim_cs=self.coord_system)
sph_coil_profile = np.concatenate((profile_order_0[..., np.newaxis], profile_orders), axis=3)
return sph_coil_profile
def get_scanner_constraints(manufacturers_model_name, order=2):
""" 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
Returns:
dict: The constraints including the scanner name, bounds and the maximum sum of currents.
"""
if manufacturers_model_name == "Prisma_fit":
constraints = {
"name": "Prisma_fit",
"coef_channel_minmax": [],
"coef_sum_max": None
}
if order >= 0:
constraints["coef_channel_minmax"].append([123100100, 123265000])
if order >= 1:
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]])
elif manufacturers_model_name == "Investigational_Device_7T":
constraints = {
"name": "Investigational_Device_7T",
"coef_channel_minmax": [],
"coef_sum_max": None
}
if order >= 0:
pass
# todo: f0 min and max is wrong
constraints["coef_channel_minmax"].append([None, None])
if order >= 1:
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]])
else:
logger.warning(f"Scanner: {manufacturers_model_name} constraints not yet implemented, constraints might not be "
"respected.")
constraints = {
"name": "Unknown",
"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)]
return constraints
def restrict_sph_constraints(bounds, order):
""" 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
orders
order (int): Maximum order of spherical harmonics
Returns:
list: 2D list with the bounds of order 0 to the specified order
"""
if order == 0:
# f0 --> [1]
minmax_out = bounds[:1]
elif order == 1:
# f0, ch1, ch2, ch3 -- > [4]
minmax_out = bounds[:4]
elif order == 2:
# f0, ch1, ch2, ch3, ch4, ch5, ch6, ch7, ch8 -- > [9]
minmax_out = bounds[:9]
else:
raise NotImplementedError("Order must be between 0 and 2")
return minmax_out