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 )} )
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 )} )
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