Skip to content

Commit 3e0d751

Browse files
etpetersonetpeterson
committed
Merge branch 'main' into testing/add_ETP_and_add_configuration
2 parents 17f3033 + 6d6207b commit 3e0d751

File tree

14 files changed

+677
-1
lines changed

14 files changed

+677
-1
lines changed

‎doc/code_contributions_record.csv

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,4 +16,5 @@ IVIM,Fitting,Linear fit,Linear fit for D with extrapolation for f. Supports unit
1616
IVIM,Fitting,sIVIM fit,NLLS of the simplified IVIM model (sIVIM). Supports units in mm2/s and µm2/ms,IAR_LundUniversity,TF2.4_IVIM-MRI_CodeCollection/src/original/IAR_LundUniversity/ivim_fit_method_modified_sivim.py,Modified by Ivan A. Rashid,Lund University,IvimModelsIVIM,tba,tbd,
1717
IVIM,Fitting,Segmented NLLS fitting,MATLAB code,OJ_GU,TF2.4_IVIM-MRI_CodeCollection/src/original/OJ_GU/,Oscar Jalnefjord,University of Gothenburg,IVIM_seg,https://doi.org/10.1007/s10334-018-0697-5,tbd,
1818
IVIM,Fitting,Bayesian,MATLAB code,OJ_GU,TF2.4_IVIM-MRI_CodeCollection/src/original/OJ_GU/,Oscar Jalnefjord,University of Gothenburg,IVIM_bayes,https://doi.org/10.1002/mrm.26783,tbd,
19+
IVIM,Fitting,Segmented NLLS fitting,Specifically tailored algorithm for NLLS segmented fitting,OJ_GU,TF2.4_IVIM-MRI_CodeCollection/src/original/OJ_GU/,Oscar Jalnefjord,University of Gothenburg,seg,https://doi.org/10.1007/s10334-018-0697-5,tbd,
1920
IVIM,Fitting,Linear fit,Linear fit for D and D* and f. Intended to be extremely fast but not always accurate,ETP_SRI,TF2.4_IVIM-MRI_CodeCollection/src/original/ETP_SRI/LinearFitting.py,Eric Peterson,SRI International,,,tbd,

‎phantoms/brain/README.md

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
# Code for generating ivim phantoms of the brain
2+
3+
Two sets of ground truth data are used for simulation:
4+
- Simulation in the diffusive regime (IVIM parameters D, f and D*): reference values from Ryghög et al. 2014 and b-values from Federau et al. 2012
5+
- Simulation in the ballistic regime (IVIM parameters D, f and vd): reference values and sequence parameters from Ahlgren et al. 2016
6+
7+
The segmentation used by the simulations is the ICBM 2009a nonlinear symmetric 3T atlas (https://nist.mni.mcgill.ca/icbm-152-nonlinear-atlases-2009/), the same as in e.g. ASLDRO (https://asldro.readthedocs.io/en/stable/).

‎phantoms/brain/ground_truth/ballistic.bval

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
0 10 20 30 40 50 60 70 80 90 100 120 140 160 180 200 10 20 30 40 50 60 70 80 90 100 120 140 160 180 200

‎phantoms/brain/ground_truth/ballistic.cval

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.471 0.666 0.816 0.942 1.054 1.154 1.247 1.333 1.414 1.490 1.632 1.763 1.885 1.999 2.107

‎phantoms/brain/ground_truth/ballistic_calc_cval.py

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
import os
2+
import numpy as np
3+
4+
# Data from Ahlgren et al 2016 NMRinBiomed
5+
Delta = 7.5e-3 # 7.5 ms
6+
delta = 7.3e-3 # 7.3 ms
7+
8+
# For DDE sequence we have
9+
# b = 2y^2G^2d^2(D-d/3), c = 2yGdD => c = sqrt(b 2D^2/(D-d/3))
10+
bval_file = os.path.join(os.path.dirname(__file__),'ballistic.bval')
11+
b = np.loadtxt(bval_file)
12+
c = np.sqrt(b*2*Delta**2/(Delta-delta/3))
13+
c[1:(c.size-1)//2+1] = 0 # flow compensated => c = 0
14+
cval_file = bval_file.replace('bval','cval')
15+
np.savetxt(cval_file,c,fmt='%.3f',newline=' ')

‎phantoms/brain/ground_truth/ballistic_groundtruth.json

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
{
2+
"D":[1.2e-3,0.98e-3,3e-3],
3+
"f":[0.0243,0.0164,0],
4+
"vd":[1.71,1.73,0],
5+
"Db":1.75e-3
6+
}

‎phantoms/brain/ground_truth/diffusive.bval

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
0 10 20 30 40 80 110 140 170 200 300 400 500 600 700 800 900

‎phantoms/brain/ground_truth/diffusive_groundtruth.json

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
{
2+
"D":[0.81e-3,0.86e-3,3e-3],
3+
"f":[0.044,0.033,0],
4+
"Dstar":[84e-3,76e-3,0]
5+
}

‎phantoms/brain/sim_brain_phantom.py

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
import os
2+
import shutil
3+
import json
4+
import numpy as np
5+
import nibabel as nib
6+
from scipy.ndimage import zoom
7+
8+
DIFFUSIVE_REGIME = 'diffusive'
9+
BALLISTIC_REGIME = 'ballistic'
10+
11+
folder = os.path.dirname(__file__)
12+
13+
###########################
14+
# Simulation parameters #
15+
regime = DIFFUSIVE_REGIME
16+
snr = 200
17+
resolution = [3,3,3]
18+
# #
19+
###########################
20+
21+
22+
# Ground truth
23+
nii = nib.load(os.path.join(folder,'ground_truth','hrgt_icbm_2009a_nls_3t.nii.gz'))
24+
segmentation = np.squeeze(nii.get_fdata()[...,-1])
25+
26+
with open(os.path.join(folder,'ground_truth',regime+'_groundtruth.json'), 'r') as f:
27+
ivim_pars = json.load(f)
28+
S0 = 1
29+
30+
# Sequence parameters
31+
bval_file = os.path.join(folder,'ground_truth',regime+'.bval')
32+
b = np.loadtxt(bval_file)
33+
if regime == BALLISTIC_REGIME:
34+
cval_file = bval_file.replace('bval','cval')
35+
c = np.loadtxt(cval_file)
36+
37+
# Calculate signal
38+
S = np.zeros(list(np.shape(segmentation))+[b.size])
39+
40+
if regime == BALLISTIC_REGIME:
41+
Db = ivim_pars["Db"]
42+
for i,(D,f,vd) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["vd"])):
43+
S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Db-c**2*vd**2))
44+
else:
45+
for i,(D,f,Dstar) in enumerate(zip(ivim_pars["D"],ivim_pars["f"],ivim_pars["Dstar"])):
46+
S[segmentation==i+1,:] = S0*((1-f)*np.exp(-b*D)+f*np.exp(-b*Dstar))
47+
48+
# Resample to suitable resolution
49+
im = zoom(S,np.append(np.diag(nii.affine)[:3]/np.array(resolution),1),order=1)
50+
sz = im.shape
51+
52+
# Add noise
53+
im_noise = np.abs(im + S0/snr*(np.random.randn(sz[0],sz[1],sz[2],sz[3])+1j*np.random.randn(sz[0],sz[1],sz[2],sz[3])))
54+
55+
# Save as image and sequence parameters
56+
nii_out = nib.Nifti1Image(im_noise,np.eye(4))
57+
base_name = os.path.join(folder,'data','{}_snr{}'.format(regime,snr))
58+
nib.save(nii_out,base_name+'.nii.gz')
59+
shutil.copyfile(bval_file,base_name+'.bval')
60+
if regime == BALLISTIC_REGIME:
61+
shutil.copyfile(cval_file,base_name+'.cval')

‎src/original/OJ_GU/ivim_seg.py

Lines changed: 211 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,211 @@
1+
import numpy as np
2+
import numpy.typing as npt
3+
4+
def seg(Y, b, bthr = 200, verbose = False):
5+
"""
6+
Segmented fitting of the IVIM model.
7+
8+
Arguments:
9+
Y: v x b matrix with data
10+
b: vector of size b with b-values
11+
bthr: (optional) threshold b-value from which signal is included in first fitting step
12+
verbose: (optional) if True, diagnostics during fitting is printet to terminal
13+
"""
14+
15+
def valid_signal(Y):
16+
"""
17+
Return a mask representing all rows in Y with valid values (not non-positive, NaN or infinite).
18+
19+
Arguments:
20+
Y -- v x b matrix with data
21+
22+
Output:
23+
mask -- vector of size v indicating valid rows in Y
24+
"""
25+
26+
mask = ~np.any((Y<=0) | np.isnan(Y) | np.isinf(Y), axis=1)
27+
return mask
28+
29+
def monoexp_model(b, D):
30+
"""
31+
Return the monoexponential e^(-b*D).
32+
33+
Arguments:
34+
b: vector of b-values [s/mm2]
35+
D: ND array of diffusion coefficients [mm2/s]
36+
37+
Output:
38+
S: (N+1)D array of signal values
39+
"""
40+
b = np.atleast_1d(b)
41+
D = np.atleast_1d(D)
42+
S = np.exp(-np.outer(D, b))
43+
return np.reshape(S, list(D.shape) + [b.size]) # reshape as np.outer flattens D is ndim > 1
44+
45+
46+
def _monoexp(Y, b, lim = [0, 3e-3], validate = True, verbose = False):
47+
""" Estimate D and A (y axis intercept) from a monoexponential model. """
48+
49+
if validate:
50+
mask = valid_signal(Y) # Avoids signal with obvious errors (non-positive, nan, inf)
51+
else:
52+
mask = np.full(Y.shape[0], True)
53+
54+
D = np.full(mask.shape, np.nan)
55+
A = np.full(mask.shape, np.nan)
56+
57+
if b.size == 2:
58+
if b[1] == b[0]:
59+
raise ZeroDivisionError("Two b-values can't be equal or both zero.")
60+
D[mask] = (np.log(Y[mask, 0]) - np.log(Y[mask, 1])) / (b[1] - b[0])
61+
62+
D[mask & (D<lim[0])] = lim[0]
63+
D[mask & (D>lim[1])] = lim[1]
64+
65+
A[mask] = Y[mask, 0] * np.exp(b[0]*D[mask])
66+
elif b.size > 2:
67+
D[mask], A[mask] = _optimizeD(Y[mask, :], b, lim, disp_prog = verbose)
68+
else:
69+
raise ValueError('Too few b-values.')
70+
71+
return D, A
72+
73+
def _f_from_intercept(A, S0):
74+
""" Calculate f from S(b=0) and extrapolated y axis intercept A."""
75+
f = 1 - A/S0
76+
f[f<0] = 0
77+
return f
78+
79+
def _optimizeD(Y, b, lim, optlim = 1e-6, disp_prog = False):
80+
""" Specfically tailored function for finding the least squares solution of monoexponenital fit. """
81+
82+
n = Y.shape[0]
83+
D = np.zeros(n)
84+
yb = Y * np.tile(b, (n, 1)) # Precalculate for speed.
85+
86+
##############################################
87+
# Check if a minimum is within the interval. #
88+
##############################################
89+
# Check that all diff < 0 for Dlow.
90+
Dlow = lim[0] * np.ones(n)
91+
difflow,_ = _Ddiff(Y, yb, b, Dlow)
92+
low_check = difflow < 0 # difflow must be < 0 if the optimum is within the interval.
93+
94+
# Check that all diff > 0 for Dhigh
95+
Dhigh = lim[1] * np.ones(n)
96+
diffhigh,_ = _Ddiff(Y, yb, b, Dhigh)
97+
high_check = diffhigh > 0 # diffhigh must be > 0 if the optimum is within the interval.
98+
99+
# Set parameter value with optimum out of bounds.
100+
D[~low_check] = lim[0] # difflow > 0 means that the mimimum has been passed .
101+
D[~high_check] = lim[1] # diffhigh < 0 means that the minium is beyond the interval.
102+
103+
# Only the voxels with a possible minimum should be estimated.
104+
mask = low_check & high_check
105+
if disp_prog:
106+
print(f'Discarding {np.count_nonzero(~mask)} voxels due to parameters out of bounds.')
107+
108+
# Allocate all variables.
109+
D_lin = np.zeros(n)
110+
diff_lin = np.zeros(n)
111+
D_mid = np.zeros(n)
112+
diff_mid = np.zeros(n)
113+
ratio_lin = np.zeros(n)
114+
ratio_mid = np.zeros(n)
115+
116+
##########################################################
117+
# Iterative method for finding the point where diff = 0. #
118+
##########################################################
119+
k = 0
120+
while np.any(mask): # Continue if there are voxels left to optimize.
121+
# Assume diff is linear within the search interval [Dlow Dhigh].
122+
D_lin[mask] = Dlow[mask] - difflow[mask] * (Dhigh[mask]-Dlow[mask]) / (diffhigh[mask]-difflow[mask])
123+
# Calculate diff in the point of intersection given by the previous expression.
124+
diff_lin[mask], ratio_lin[mask] = _Ddiff(Y[mask, :], yb[mask, :], b, D_lin[mask])
125+
126+
# As a potential speed up, the mean of Dlow and Dhigh is also calculated.
127+
D_mid[mask] = (Dlow[mask]+Dhigh[mask]) / 2
128+
diff_mid[mask], ratio_mid[mask] = _Ddiff(Y[mask, :], yb[mask, :], b, D_mid[mask])
129+
130+
# If diff < 0, then the point of intersection or mean is used as the
131+
# new Dlow. Only voxels with diff < 0 are updated at this step. Linear
132+
# interpolation or the mean is used depending of which method that
133+
# gives the smallest diff.
134+
updatelow_lin = (diff_lin<0) & ((diff_mid>0) | ((D_lin>D_mid) & (diff_mid<0)))
135+
updatelow_mid = (diff_mid<0) & ((diff_lin>0) | ((D_mid>D_lin) & (diff_lin<0)))
136+
Dlow[updatelow_lin] = D_lin[updatelow_lin]
137+
Dlow[updatelow_mid] = D_mid[updatelow_mid]
138+
139+
# If diff > 0, then the point of intersection or mean is used as the
140+
# new Dhigh. Only voxels with diff > 0 are updated at this step.
141+
# Linear interpolation or the mean is used depending of which method
142+
# that gives the smallest diff.
143+
updatehigh_lin = (diff_lin>0) & ((diff_mid<0) | ((D_lin<D_mid) & (diff_mid>0)))
144+
updatehigh_mid = (diff_mid>0) & ((diff_lin<0) | ((D_mid<D_lin) & (diff_lin>0)))
145+
Dhigh[updatehigh_lin] = D_lin[updatehigh_lin]
146+
Dhigh[updatehigh_mid] = D_mid[updatehigh_mid]
147+
148+
# Update the mask to exclude voxels that fulfills the optimization
149+
# limit from the mask.
150+
opt_lin = np.abs(1-ratio_lin) < optlim
151+
opt_mid = np.abs(1-ratio_mid) < optlim
152+
153+
D[opt_lin] = D_lin[opt_lin]
154+
D[opt_mid] = D_mid[opt_mid]
155+
# Not optimal if both D_lin and D_mean fulfills the optimization limit,
156+
# but has a small impact on the result as long as optlim is small.
157+
158+
# Update the mask.
159+
mask = mask & (~(opt_lin|opt_mid))
160+
161+
# Calculate diff for the new bounds.
162+
if np.any(mask):
163+
difflow[mask],_ = _Ddiff(Y[mask, :], yb[mask, :], b, Dlow[mask])
164+
diffhigh[mask],_ = _Ddiff(Y[mask, :], yb[mask, :], b, Dhigh[mask])
165+
166+
k += 1
167+
if disp_prog:
168+
print(f'Iteration {k}: {np.count_nonzero(mask)} voxels left.')
169+
170+
A = np.sum(Y*np.exp(-np.outer(D, b)), axis=1) / np.sum(np.exp(-2*np.outer(b, D)), axis=0)
171+
172+
return D, A
173+
174+
175+
def _Ddiff(Y, yb, b, D):
176+
"""
177+
Return the difference between q1 = e^(-2*b*D)*yb*e^(-b*D) and
178+
q2 = Y*e^(-b*D)*b*e^(-2*b*D) summed over b as well as the ratio q1/q2
179+
summed over b, setting divisions by zero as infinite.
180+
"""
181+
q1 = np.sum(np.exp(-2*np.outer(b, D)), axis=0) * np.sum(yb*np.exp(-np.outer(D, b)), axis=1)
182+
q2 = np.sum(Y*np.exp(-np.outer(D, b)), axis=1) * np.sum(b[:, np.newaxis]*np.exp(-2*np.outer(b, D)), axis=0)
183+
diff = q1 - q2
184+
ratio = np.full(q1.shape, np.inf)
185+
ratio[q2!=0] = q1[q2!=0] / q2[q2!=0]
186+
return diff, ratio
187+
188+
if b.ndim != 1:
189+
raise ValueError('b must a 1D array')
190+
191+
if Y.ndim != 2:
192+
if Y.ndim != 1:
193+
raise ValueError('Y must be a 2D array')
194+
else:
195+
Y = Y[np.newaxis,:]
196+
197+
if b.size != Y.shape[1]:
198+
raise ValueError('Number of b-values must match the second dimension of Y.')
199+
200+
bmask = b >= bthr
201+
202+
D, A = _monoexp(Y[:, bmask], b[bmask], verbose=verbose)
203+
Ysub = (Y - A[:, np.newaxis]*monoexp_model(b, D)) # Remove signal related to diffusion
204+
205+
Dstar, Astar = _monoexp(Ysub, b, lim=[3e-3, 0.1], validate = False, verbose = verbose)
206+
S0 = A + Astar
207+
208+
f = _f_from_intercept(A, S0)
209+
210+
pars = {'D': D, 'f': f, 'Dstar': Dstar, 'S0': S0}
211+
return pars

0 commit comments

Comments
 (0)