Skip to content

Commit 3360ba8

Browse files
acsenrafilhoacsenrafilho
authored
Merge pull request #23 from LOAMRI/develop
Develop
2 parents cd1404d + 24a0d58 commit 3360ba8

23 files changed

+1407
-374
lines changed

‎asltk/asldata.py‎

Lines changed: 20 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -82,13 +82,21 @@ def __init__(
8282
if kwargs.get('dw_values'):
8383
self._parameters['dw'] = kwargs.get('dw_values')
8484

85-
def set_image(self, full_path: str, spec: str):
85+
def set_image(self, image, spec: str):
8686
"""Insert a image necessary to define de ASL data processing.
8787
8888
The `spec` parameters specifies what is the type of image to be used in
8989
ASL processing step. Choose one of the options: `m0` for the M0 volume,
9090
`pcasl` for the pCASL data.
9191
92+
Note:
93+
The image can be a full path to the image file or a numpy array.
94+
In case the image parameter is a path, then the method will load
95+
the image file directly and associate it with the `spec` parameter.
96+
However, if the image is a numpy array, then the method will
97+
pass it to the ASLData object image data regarding the `spec`
98+
parameter as well.
99+
92100
Examples:
93101
>>> data = ASLData()
94102
>>> path_m0 = './tests/files/m0.nii.gz' # M0 file with shape (5,35,35)
@@ -97,13 +105,19 @@ def set_image(self, full_path: str, spec: str):
97105
(5, 35, 35)
98106
99107
Args:
100-
full_path (str): The full path with filename to be loaded
108+
image (str): The image to be used.
101109
spec (str): The type of image being used in the ASL processing.
102110
"""
103-
if spec == 'm0':
104-
self._m0_image = load_image(full_path)
105-
elif spec == 'pcasl':
106-
self._asl_image = load_image(full_path)
111+
if isinstance(image, str) and os.path.exists(image):
112+
if spec == 'm0':
113+
self._m0_image = load_image(image)
114+
elif spec == 'pcasl':
115+
self._asl_image = load_image(image)
116+
elif isinstance(image, np.ndarray):
117+
if spec == 'm0':
118+
self._m0_image = image
119+
elif spec == 'pcasl':
120+
self._asl_image = image
107121

108122
def get_ld(self):
109123
"""Obtain the LD array values"""

‎asltk/models/__init__.py‎

Whitespace-only changes.

‎asltk/models/signal_dynamic.py‎

Lines changed: 336 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,336 @@
1+
import math
2+
3+
import numpy as np
4+
5+
6+
def asl_model_buxton(
7+
tau: list,
8+
w: list,
9+
m0: float,
10+
cbf: float,
11+
att: float,
12+
lambda_value: float = 0.98,
13+
t1b: float = 1650.0,
14+
alpha: float = 0.85,
15+
):
16+
"""Buxton model to calculate the ASL magnetization values.
17+
18+
It is assumed that the LD and PLD values are coherent with the ASl Buxton
19+
model, i.e. the both has the same array size.
20+
21+
The calculations is given assuming a voxel value. Hence, all the `tau`,
22+
`w`, `cbf` and `att` values must representas a voxel in the image.
23+
24+
Note:
25+
The CBF value is the original scale, without assuming the normalized
26+
CBF value. See more details at the CBFMapping class documentation.
27+
28+
Args:
29+
tau (list): LD values
30+
w (list): PLD values
31+
m0 (float): The M0 magnetization value
32+
cbf (float): The CBF value, not been assumed as normalized.
33+
att (float): The ATT value
34+
lambda_value (float, optional): The blood-brain partition coefficient (0 to 1.0). Defaults to 0.98.
35+
t1b (float, optional): The T1 relaxation value of the blood. Defaults to 1650.0.
36+
alpha (float, optional): The labeling efficiency. Defaults to 0.85.
37+
38+
Returns:
39+
(numpy.ndarray): A numpy array with the magnetization values calculated
40+
"""
41+
tau = tau.tolist() if isinstance(tau, np.ndarray) else tau
42+
w = w.tolist() if isinstance(w, np.ndarray) else w
43+
44+
if not (isinstance(tau, list) ^ isinstance(tau, tuple)):
45+
raise ValueError('tau parameter must be a list or tuple of values.')
46+
47+
if not isinstance(w, list) ^ isinstance(w, tuple):
48+
raise ValueError('w parameter must be a list or tuple of values.')
49+
50+
for v in tau:
51+
if not isinstance(v, float) ^ isinstance(v, int):
52+
raise ValueError('tau list must contain float or int values')
53+
54+
for v in w:
55+
if not isinstance(v, float) ^ isinstance(v, int):
56+
raise ValueError('w list must contain float or int values')
57+
58+
# if len(tau) != len(w):
59+
# raise SyntaxError("tau and w parameters must be at the same size.")
60+
61+
t = np.add(tau, w).tolist()
62+
63+
t1bp = 1 / ((1 / t1b) + (cbf / lambda_value))
64+
m_values = np.zeros(len(tau))
65+
66+
for i in range(0, len(tau)):
67+
try:
68+
if t[i] < att:
69+
m_values[i] = 0.0
70+
elif (att <= t[i]) and (t[i] < tau[i] + att):
71+
q = 1 - math.exp(-(t[i] - att) / t1bp)
72+
m_values[i] = (
73+
2.0 * m0 * cbf * t1bp * alpha * q * math.exp(-att / t1b)
74+
)
75+
else:
76+
q = 1 - math.exp(-tau[i] / t1bp)
77+
m_values[i] = (
78+
2.0
79+
* m0
80+
* cbf
81+
* t1bp
82+
* alpha
83+
* q
84+
* math.exp(-att / t1b)
85+
* math.exp(-(t[i] - tau[i] - att) / t1bp)
86+
)
87+
except OverflowError: # pragma: no cover
88+
m_values[i] = 0.0
89+
90+
return m_values
91+
92+
93+
def asl_model_multi_te(
94+
tau: list,
95+
w: list,
96+
te: list,
97+
m0: float,
98+
cbf: float,
99+
att: float,
100+
t2b: float = 165.0,
101+
t2csf: float = 75.0,
102+
tblcsf: float = 1400.0,
103+
alpha: float = 0.85,
104+
t1b: float = 1650.0,
105+
t1csf: float = 1400.0,
106+
):
107+
"""Multi Time of Echos (TE) ASL model to calculate the T1 relaxation time for
108+
blood and Grey Matter exchange.
109+
110+
This model is directly used on the MultiTE_ASLMapping class.
111+
112+
Reference: Ultra-long-TE arterial spin labeling reveals rapid and
113+
brain-wide blood-to-CSF water transport in humans, NeuroImage,
114+
doi: 10.1016/j.neuroimage.2021.118755
115+
116+
Args:
117+
tau (list): The LD values
118+
w (list): The PLD values
119+
te (list): The TE values
120+
m0 (float): The M0 voxel value
121+
cbf (float): The CBF voxel value
122+
att (float): The ATT voxel value
123+
t2b (float, optional): The T2 relaxation value for blood. Defaults to 165.0.
124+
t2csf (float, optional): The T2 relaxation value for CSF. Defaults to 75.0.
125+
tblcsf (float, optional): The T1 relaxation value between blood and CSF. Defaults to 1400.0.
126+
alpha (float, optional): The pulse labeling efficiency. Defaults to 0.85.
127+
t1b (float, optional): The T1 relaxation value for blood. Defaults to 1650.0.
128+
t1csf (float, optional): The T1 relaxation value for CSF. Defaults to 1400.0.
129+
130+
Returns:
131+
(nd.ndarray): The magnetization values for T1-Blood-GM
132+
"""
133+
t1bp = 1 / ((1 / t1b) + (1 / tblcsf))
134+
t1csfp = 1 / ((1 / t1csf) + (1 / tblcsf))
135+
136+
t2bp = 1 / ((1 / t2b) + (1 / tblcsf))
137+
t2csfp = 1 / ((1 / t2csf) + (1 / tblcsf))
138+
139+
t = np.add(tau, w).tolist()
140+
141+
mag_total = np.zeros(len(tau))
142+
143+
for i in range(0, len(tau)):
144+
try:
145+
if t[i] < att:
146+
S1b = 0.0
147+
S1csf = 0.0
148+
if te[i] < (att - t[i]):
149+
Sb = 0
150+
Scsf = 0
151+
elif (att - t[i]) <= te[i] and te[i] < (att + tau[i] - t[i]):
152+
Sb = (
153+
2
154+
* alpha
155+
* m0
156+
* cbf
157+
* t2bp
158+
* math.exp(-att / t1b)
159+
* math.exp(-te[i] / t2b)
160+
* (1 - math.exp(-(te[i] - att + t[i]) / t2bp))
161+
) #% measured signal = S2
162+
Scsf = (
163+
2
164+
* alpha
165+
* m0
166+
* cbf
167+
* math.exp(-att / t1b)
168+
* math.exp(-te[i] / t2b)
169+
* (
170+
t2csf
171+
* (1 - math.exp(-(te[i] - att + t[i]) / t2csf))
172+
- t2csfp
173+
* (1 - math.exp(-(te[i] - att + t[i]) / t2csfp))
174+
)
175+
)
176+
else: #% att + tau - t <= te
177+
Sb = (
178+
2
179+
* alpha
180+
* m0
181+
* cbf
182+
* t2bp
183+
* math.exp(-att / t1b)
184+
* math.exp(-te[i] / t2b)
185+
* math.exp(-(te[i] - att + t[i]) / t2bp)
186+
* (math.exp(tau[i] / t2bp) - 1)
187+
)
188+
Scsf = (
189+
2
190+
* alpha
191+
* m0
192+
* cbf
193+
* math.exp(-att / t1b)
194+
* math.exp(-te[i] / t2b)
195+
* (
196+
t2csf
197+
* math.exp(-(te[i] - att + t[i]) / t2csf)
198+
* (math.exp(tau[i] / t2csf) - 1)
199+
- t2csfp
200+
* math.exp(-(te[i] - att + t[i]) / t2csfp)
201+
* (math.exp(tau[i] / t2csfp) - 1)
202+
)
203+
)
204+
elif (att <= t[i]) and (t[i] < (att + tau[i])):
205+
S1b = (
206+
2
207+
* alpha
208+
* m0
209+
* cbf
210+
* t1bp
211+
* math.exp(-att / t1b)
212+
* (1 - math.exp(-(t[i] - att) / t1bp))
213+
)
214+
S1csf = (
215+
2
216+
* alpha
217+
* m0
218+
* cbf
219+
* math.exp(-att / t1b)
220+
* (
221+
t1csf * (1 - math.exp(-(t[i] - att) / t1csf))
222+
- t1csfp * (1 - math.exp(-(t[i] - att) / t1csfp))
223+
)
224+
)
225+
if te[i] < (att + tau[i] - t[i]):
226+
Sb = S1b * math.exp(
227+
-te[i] / t2bp
228+
) + 2 * alpha * m0 * cbf * t2bp * math.exp(
229+
-att / t1b
230+
) * math.exp(
231+
-te[i] / t2b
232+
) * (
233+
1 - math.exp(-te[i] / t2bp)
234+
)
235+
Scsf = (
236+
S1b
237+
* (1 - math.exp(-te[i] / tblcsf))
238+
* math.exp(-te[i] / t2csf)
239+
+ S1csf * math.exp(-te[i] / t2csf)
240+
+ 2
241+
* alpha
242+
* m0
243+
* cbf
244+
* math.exp(-att / t1b)
245+
* math.exp(-te[i] / t2b)
246+
* (
247+
t2csf * (1 - math.exp(-te[i] / t2csf))
248+
- t2csfp * (1 - math.exp(-te[i] / t2csfp))
249+
)
250+
)
251+
else: # att + tau - t <= te
252+
Sb = S1b * math.exp(
253+
-te[i] / t2bp
254+
) + 2 * alpha * m0 * cbf * t2bp * math.exp(
255+
-att / t1b
256+
) * math.exp(
257+
-te[i] / t2b
258+
) * math.exp(
259+
-te[i] / t2bp
260+
) * (
261+
math.exp((att + tau[i] - t[i]) / t2bp) - 1
262+
)
263+
Scsf = (
264+
S1b
265+
* (1 - math.exp(-te[i] / tblcsf))
266+
* math.exp(-te[i] / t2csf)
267+
+ S1csf * math.exp(-te[i] / t2csf)
268+
+ 2
269+
* alpha
270+
* m0
271+
* cbf
272+
* math.exp(-att / t1b)
273+
* math.exp(-te[i] / t2b)
274+
* (
275+
t2csf
276+
* math.exp(-te[i] / t2csf)
277+
* (math.exp((att + tau[i] - t[i]) / t2csf) - 1)
278+
- t2csfp
279+
* math.exp(-te[i] / t2csfp)
280+
* (math.exp((att + tau[i] - t[i]) / t2csfp) - 1)
281+
)
282+
)
283+
else: # att+tau < t
284+
S1b = (
285+
2
286+
* alpha
287+
* m0
288+
* cbf
289+
* t1bp
290+
* math.exp(-att / t1b)
291+
* math.exp(-(t[i] - att) / t1bp)
292+
* (math.exp(tau[i] / t1bp) - 1)
293+
)
294+
S1csf = (
295+
2
296+
* alpha
297+
* m0
298+
* cbf
299+
* math.exp(-att / t1b)
300+
* (
301+
t1csf
302+
* math.exp(-(t[i] - att) / t1csf)
303+
* (math.exp(tau[i] / t1csf) - 1)
304+
- t1csfp
305+
* math.exp(-(t[i] - att) / t1csfp)
306+
* (math.exp(tau[i] / t1csfp) - 1)
307+
)
308+
)
309+
310+
Sb = S1b * math.exp(-te[i] / t2bp)
311+
Scsf = S1b * (1 - math.exp(-te[i] / tblcsf)) * math.exp(
312+
-te[i] / t2csf
313+
) + S1csf * math.exp(-te[i] / t2csf)
314+
except (OverflowError, RuntimeError): # pragma: no cover
315+
Sb = 0.0
316+
Scsf = 0.0
317+
318+
mag_total[i] = Sb + Scsf
319+
320+
return mag_total
321+
322+
323+
def asl_model_multi_dw(
324+
b_values: list, A1: list, D1: float, A2: list, D2: float
325+
):
326+
mag_total = np.zeros(len(b_values))
327+
328+
for i in range(0, len(b_values)):
329+
try:
330+
mag_total[i] = A1 * math.exp(-b_values[i] * D1) + A2 * math.exp(
331+
-b_values[i] * D2
332+
)
333+
except (OverflowError, RuntimeError): # pragma: no cover
334+
mag_total[i] = 0.0
335+
336+
return mag_total

‎asltk/reconstruction/cbf_mapping.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,8 @@
88

99
from asltk.asldata import ASLData
1010
from asltk.aux_methods import _check_mask_values
11+
from asltk.models.signal_dynamic import asl_model_buxton
1112
from asltk.mri_parameters import MRIParameters
12-
from asltk.utils import asl_model_buxton
1313

1414
# Global variables to assist multi cpu threading
1515
cbf_map = None

0 commit comments

Comments
 (0)