1+ import numpy as np
2+ import numpy .typing as npt
3+
4+ def seg (Y : npt .NDArray [np .float64 ], b : npt .NDArray [np .float64 ], bthr : float = 200 , verbose : bool = False ) -> dict :
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 : npt .NDArray [np .float64 ]) -> npt .NDArray [np .float64 ]:
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 : npt .NDArray [np .float64 ], D : npt .NDArray [np .float64 ]) -> npt .NDArray [np .float64 ]:
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 : npt .NDArray [np .float64 ], b : npt .NDArray [np .float64 ], lim : list = [0 , 3e-3 ], validate : bool = True , verbose : bool = False ) -> tuple [npt .NDArray [np .float64 ], npt .NDArray [np .float64 ]]:
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 _get_S0 (Y : npt .NDArray [np .float64 ], b : npt .NDArray [np .float64 ]) -> npt .NDArray [np .float64 ]:
74+ """ Return the signal values at b = 0."""
75+ return np .mean (Y [:, b == 0 ], axis = 1 )
76+
77+
78+ def _f_from_intercept (A : npt .NDArray [np .float64 ], S0 : npt .NDArray [np .float64 ]) -> npt .NDArray [np .float64 ]:
79+ """ Calculate f from S(b=0) and extrapolated y axis intercept A."""
80+ f = 1 - A / S0
81+ f [f < 0 ] = 0
82+ return f
83+
84+ def _optimizeD (Y : npt .NDArray [np .float64 ], b : npt .NDArray [np .float64 ], lim : list , optlim : float = 1e-6 , disp_prog : bool = False ) -> tuple [npt .NDArray [np .float64 ], npt .NDArray [np .float64 ]]:
85+ """ Specfically tailored function for finding the least squares solution of monoexponenital fit. """
86+
87+ n = Y .shape [0 ]
88+ D = np .zeros (n )
89+ yb = Y * np .tile (b , (n , 1 )) # Precalculate for speed.
90+
91+ ##############################################
92+ # Check if a minimum is within the interval. #
93+ ##############################################
94+ # Check that all diff < 0 for Dlow.
95+ Dlow = lim [0 ] * np .ones (n )
96+ difflow ,_ = _Ddiff (Y , yb , b , Dlow )
97+ low_check = difflow < 0 # difflow must be < 0 if the optimum is within the interval.
98+
99+ # Check that all diff > 0 for Dhigh
100+ Dhigh = lim [1 ] * np .ones (n )
101+ diffhigh ,_ = _Ddiff (Y , yb , b , Dhigh )
102+ high_check = diffhigh > 0 # diffhigh must be > 0 if the optimum is within the interval.
103+
104+ # Set parameter value with optimum out of bounds.
105+ D [~ low_check ] = lim [0 ] # difflow > 0 means that the mimimum has been passed .
106+ D [~ high_check ] = lim [1 ] # diffhigh < 0 means that the minium is beyond the interval.
107+
108+ # Only the voxels with a possible minimum should be estimated.
109+ mask = low_check & high_check
110+ if disp_prog :
111+ print (f'Discarding { np .count_nonzero (~ mask )} )
112+
113+ # Allocate all variables.
114+ D_lin = np .zeros (n )
115+ diff_lin = np .zeros (n )
116+ D_mid = np .zeros (n )
117+ diff_mid = np .zeros (n )
118+ ratio_lin = np .zeros (n )
119+ ratio_mid = np .zeros (n )
120+
121+ ##########################################################
122+ # Iterative method for finding the point where diff = 0. #
123+ ##########################################################
124+ k = 0
125+ while np .any (mask ): # Continue if there are voxels left to optimize.
126+ # Assume diff is linear within the search interval [Dlow Dhigh].
127+ D_lin [mask ] = Dlow [mask ] - difflow [mask ] * (Dhigh [mask ]- Dlow [mask ]) / (diffhigh [mask ]- difflow [mask ])
128+ # Calculate diff in the point of intersection given by the previous expression.
129+ diff_lin [mask ], ratio_lin [mask ] = _Ddiff (Y [mask , :], yb [mask , :], b , D_lin [mask ])
130+
131+ # As a potential speed up, the mean of Dlow and Dhigh is also calculated.
132+ D_mid [mask ] = (Dlow [mask ]+ Dhigh [mask ]) / 2
133+ diff_mid [mask ], ratio_mid [mask ] = _Ddiff (Y [mask , :], yb [mask , :], b , D_mid [mask ])
134+
135+ # If diff < 0, then the point of intersection or mean is used as the
136+ # new Dlow. Only voxels with diff < 0 are updated at this step. Linear
137+ # interpolation or the mean is used depending of which method that
138+ # gives the smallest diff.
139+ updatelow_lin = (diff_lin < 0 ) & ((diff_mid > 0 ) | ((D_lin > D_mid ) & (diff_mid < 0 )))
140+ updatelow_mid = (diff_mid < 0 ) & ((diff_lin > 0 ) | ((D_mid > D_lin ) & (diff_lin < 0 )))
141+ Dlow [updatelow_lin ] = D_lin [updatelow_lin ]
142+ Dlow [updatelow_mid ] = D_mid [updatelow_mid ]
143+
144+ # If diff > 0, then the point of intersection or mean is used as the
145+ # new Dhigh. Only voxels with diff > 0 are updated at this step.
146+ # Linear interpolation or the mean is used depending of which method
147+ # that gives the smallest diff.
148+ updatehigh_lin = (diff_lin > 0 ) & ((diff_mid < 0 ) | ((D_lin < D_mid ) & (diff_mid > 0 )))
149+ updatehigh_mid = (diff_mid > 0 ) & ((diff_lin < 0 ) | ((D_mid < D_lin ) & (diff_lin > 0 )))
150+ Dhigh [updatehigh_lin ] = D_lin [updatehigh_lin ]
151+ Dhigh [updatehigh_mid ] = D_mid [updatehigh_mid ]
152+
153+ # Update the mask to exclude voxels that fulfills the optimization
154+ # limit from the mask.
155+ opt_lin = np .abs (1 - ratio_lin ) < optlim
156+ opt_mid = np .abs (1 - ratio_mid ) < optlim
157+
158+ D [opt_lin ] = D_lin [opt_lin ]
159+ D [opt_mid ] = D_mid [opt_mid ]
160+ # Not optimal if both D_lin and D_mean fulfills the optimization limit,
161+ # but has a small impact on the result as long as optlim is small.
162+
163+ # Update the mask.
164+ mask = mask & (~ (opt_lin | opt_mid ))
165+
166+ # Calculate diff for the new bounds.
167+ if np .any (mask ):
168+ difflow [mask ],_ = _Ddiff (Y [mask , :], yb [mask , :], b , Dlow [mask ])
169+ diffhigh [mask ],_ = _Ddiff (Y [mask , :], yb [mask , :], b , Dhigh [mask ])
170+
171+ k += 1
172+ if disp_prog :
173+ print (f'Iteration { k } { np .count_nonzero (mask )} )
174+
175+ A = np .sum (Y * np .exp (- np .outer (D , b )), axis = 1 ) / np .sum (np .exp (- 2 * np .outer (b , D )), axis = 0 )
176+
177+ return D , A
178+
179+
180+ def _Ddiff (Y : npt .NDArray [np .float64 ], yb : npt .NDArray [np .float64 ], b : npt .NDArray [np .float64 ], D : npt .NDArray [np .float64 ]):
181+ """
182+ Return the difference between q1 = e^(-2*b*D)*yb*e^(-b*D) and
183+ q2 = Y*e^(-b*D)*b*e^(-2*b*D) summed over b as well as the ratio q1/q2
184+ summed over b, setting divisions by zero as infinite.
185+ """
186+ q1 = np .sum (np .exp (- 2 * np .outer (b , D )), axis = 0 ) * np .sum (yb * np .exp (- np .outer (D , b )), axis = 1 )
187+ 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 )
188+ diff = q1 - q2
189+ ratio = np .full (q1 .shape , np .inf )
190+ ratio [q2 != 0 ] = q1 [q2 != 0 ] / q2 [q2 != 0 ]
191+ return diff , ratio
192+
193+ if b .ndim != 1 :
194+ raise ValueError ('b must a 1D array' )
195+
196+ if Y .ndim != 2 :
197+ if Y .ndim != 1 :
198+ raise ValueError ('Y must be a 2D array' )
199+ else :
200+ Y = Y [np .newaxis ,:]
201+
202+ if b .size != Y .shape [1 ]:
203+ raise ValueError ('Number of b-values must match the second dimension of Y.' )
204+
205+ bmask = b >= bthr
206+
207+ D , A = _monoexp (Y [:, bmask ], b [bmask ], verbose = verbose )
208+ Ysub = (Y - A [:, np .newaxis ]* monoexp_model (b , D )) # Remove signal related to diffusion
209+
210+ Dstar , Astar = _monoexp (Ysub , b , lim = [3e-3 , 0.1 ], validate = False , verbose = verbose )
211+ S0 = A + Astar
212+
213+ f = _f_from_intercept (A , S0 )
214+
215+ pars = {'D' : D , 'f' : f , 'Dstar' : Dstar , 'S0' : S0 }
216+ return pars
0 commit comments