1#!/usr/bin/python 2""" Classes and functions for fitting the diffusion kurtosis model """ 3 4import numpy as np 5import scipy.optimize as opt 6import dipy.core.sphere as dps 7from dipy.reconst.dti import (TensorFit, mean_diffusivity, 8 from_lower_triangular, 9 lower_triangular, decompose_tensor, 10 MIN_POSITIVE_SIGNAL, nlls_fit_tensor, 11 restore_fit_tensor) 12from dipy.reconst.utils import dki_design_matrix as design_matrix 13from dipy.reconst.recspeed import local_maxima 14from dipy.reconst.base import ReconstModel 15from dipy.core.ndindex import ndindex 16from dipy.core.geometry import (sphere2cart, cart2sphere, 17 perpendicular_directions) 18from dipy.data import get_sphere, get_fnames 19from dipy.reconst.vec_val_sum import vec_val_vect 20from dipy.core.gradients import check_multi_b 21 22 23def _positive_evals(L1, L2, L3, er=2e-7): 24 """ Helper function that indentifies which voxels in a array have all 25 eigenvalues significantly larger than zero 26 27 Parameters 28 ---------- 29 L1 : ndarray 30 First independent variable of the integral. 31 L2 : ndarray 32 Second independent variable of the integral. 33 L3 : ndarray 34 Third independent variable of the integral. 35 er : float, optional 36 A eigenvalues is classified as larger than zero if it is larger than er 37 38 Returns 39 ------- 40 ind : boolean (n,) 41 Array that marks the voxels that have all eigenvalues are larger than 42 zero. 43 """ 44 ind = np.logical_and(L1 > er, np.logical_and(L2 > er, L3 > er)) 45 46 return ind 47 48 49def carlson_rf(x, y, z, errtol=3e-4): 50 r""" Computes the Carlson's incomplete elliptic integral of the first kind 51 defined as: 52 53 .. math:: 54 55 R_F = \frac{1}{2} \int_{0}^{\infty} \left [(t+x)(t+y)(t+z) \right ] 56 ^{-\frac{1}{2}}dt 57 58 Parameters 59 ---------- 60 x : ndarray 61 First independent variable of the integral. 62 y : ndarray 63 Second independent variable of the integral. 64 z : ndarray 65 Third independent variable of the integral. 66 errtol : float 67 Error tolerance. Integral is computed with relative error less in 68 magnitude than the defined value 69 70 Returns 71 ------- 72 RF : ndarray 73 Value of the incomplete first order elliptic integral 74 75 Notes 76 ------ 77 x, y, and z have to be nonnegative and at most one of them is zero. 78 79 References 80 ---------- 81 .. [1] Carlson, B.C., 1994. Numerical computation of real or complex 82 elliptic integrals. arXiv:math/9409227 [math.CA] 83 """ 84 xn = x.copy() 85 yn = y.copy() 86 zn = z.copy() 87 An = (xn + yn + zn) / 3.0 88 Q = (3. * errtol) ** (-1 / 6.) * \ 89 np.max(np.abs([An - xn, An - yn, An - zn]), axis=0) 90 # Convergence has to be done voxel by voxel 91 index = ndindex(x.shape) 92 for v in index: 93 n = 0 94 # Convergence condition 95 while 4.**(-n) * Q[v] > abs(An[v]): 96 xnroot = np.sqrt(xn[v]) 97 ynroot = np.sqrt(yn[v]) 98 znroot = np.sqrt(zn[v]) 99 lamda = xnroot * (ynroot + znroot) + ynroot * znroot 100 n = n + 1 101 xn[v] = (xn[v] + lamda) * 0.250 102 yn[v] = (yn[v] + lamda) * 0.250 103 zn[v] = (zn[v] + lamda) * 0.250 104 An[v] = (An[v] + lamda) * 0.250 105 106 # post convergence calculation 107 X = 1. - xn / An 108 Y = 1. - yn / An 109 Z = - X - Y 110 E2 = X * Y - Z * Z 111 E3 = X * Y * Z 112 RF = An**(-1 / 2.) * \ 113 (1 - E2 / 10. + E3 / 14. + (E2**2) / 24. - 3 / 44. * E2 * E3) 114 115 return RF 116 117 118def carlson_rd(x, y, z, errtol=1e-4): 119 r""" Computes the Carlson's incomplete elliptic integral of the second kind 120 defined as: 121 122 .. math:: 123 124 R_D = \frac{3}{2} \int_{0}^{\infty} (t+x)^{-\frac{1}{2}} 125 (t+y)^{-\frac{1}{2}}(t+z) ^{-\frac{3}{2}} 126 127 Parameters 128 ---------- 129 x : ndarray 130 First independent variable of the integral. 131 y : ndarray 132 Second independent variable of the integral. 133 z : ndarray 134 Third independent variable of the integral. 135 errtol : float 136 Error tolerance. Integral is computed with relative error less in 137 magnitude than the defined value 138 139 Returns 140 ------- 141 RD : ndarray 142 Value of the incomplete second order elliptic integral 143 144 Notes 145 ------ 146 x, y, and z have to be nonnegative and at most x or y is zero. 147 """ 148 xn = x.copy() 149 yn = y.copy() 150 zn = z.copy() 151 A0 = (xn + yn + 3. * zn) / 5.0 152 An = A0.copy() 153 Q = (errtol / 4.) ** (-1 / 6.) * \ 154 np.max(np.abs([An - xn, An - yn, An - zn]), axis=0) 155 sum_term = np.zeros(x.shape, dtype=x.dtype) 156 n = np.zeros(x.shape) 157 158 # Convergence has to be done voxel by voxel 159 index = ndindex(x.shape) 160 for v in index: 161 # Convergence condition 162 while 4.**(-n[v]) * Q[v] > abs(An[v]): 163 xnroot = np.sqrt(xn[v]) 164 ynroot = np.sqrt(yn[v]) 165 znroot = np.sqrt(zn[v]) 166 lamda = xnroot * (ynroot + znroot) + ynroot * znroot 167 sum_term[v] = sum_term[v] + \ 168 4.**(-n[v]) / (znroot * (zn[v] + lamda)) 169 n[v] = n[v] + 1 170 xn[v] = (xn[v] + lamda) * 0.250 171 yn[v] = (yn[v] + lamda) * 0.250 172 zn[v] = (zn[v] + lamda) * 0.250 173 An[v] = (An[v] + lamda) * 0.250 174 175 # post convergence calculation 176 X = (A0 - x) / (4.**(n) * An) 177 Y = (A0 - y) / (4.**(n) * An) 178 Z = - (X + Y) / 3. 179 E2 = X * Y - 6. * Z * Z 180 E3 = (3. * X * Y - 8. * Z * Z) * Z 181 E4 = 3. * (X * Y - Z * Z) * Z**2. 182 E5 = X * Y * Z**3. 183 RD = \ 184 4**(-n) * An**(-3 / 2.) * \ 185 (1 - 3 / 14. * E2 + 1 / 6. * E3 + 186 9 / 88. * (E2**2) - 3 / 22. * E4 - 9 / 52. * E2 * E3 + 187 3 / 26. * E5) + 3 * sum_term 188 189 return RD 190 191 192def _F1m(a, b, c): 193 r""" Helper function that computes function $F_1$ which is required to 194 compute the analytical solution of the Mean kurtosis. 195 196 Parameters 197 ---------- 198 a : ndarray 199 Array containing the values of parameter $\lambda_1$ of function $F_1$ 200 b : ndarray 201 Array containing the values of parameter $\lambda_2$ of function $F_1$ 202 c : ndarray 203 Array containing the values of parameter $\lambda_3$ of function $F_1$ 204 205 Returns 206 ------- 207 F1 : ndarray 208 Value of the function $F_1$ for all elements of the arrays a, b, and c 209 210 Notes 211 -------- 212 Function $F_1$ is defined as [1]_: 213 214 .. math:: 215 216 F_1(\lambda_1,\lambda_2,\lambda_3)= 217 \frac{(\lambda_1+\lambda_2+\lambda_3)^2} 218 {18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)} 219 [\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1} 220 R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ 221 \frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3- 222 \lambda_1\lambda_3} 223 {3\lambda_1 \sqrt{\lambda_2 \lambda_3}} 224 R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ] 225 226 References 227 ---------- 228 .. [1] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 229 Estimation of tensors and tensor-derived measures in diffusional 230 kurtosis imaging. Magn Reson Med. 65(3), 823-836 231 """ 232 # Eigenvalues are considered equal if they are not 2.5% different to each 233 # other. This value is adjusted according to the analysis reported in: 234 # http://gsoc2015dipydki.blogspot.co.uk/2015/08/rnh-post-13-start-wrapping-up-test.html 235 er = 2.5e-2 236 237 # Initialize F1 238 F1 = np.zeros(a.shape) 239 240 # Only computes F1 in voxels that have all eigenvalues larger than zero 241 cond0 = _positive_evals(a, b, c) 242 243 # Apply formula for non problematic plausible cases, i.e. a!=b and a!=c 244 cond1 = np.logical_and(cond0, np.logical_and(abs(a - b) >= a * er, 245 abs(a - c) >= a * er)) 246 if np.sum(cond1) != 0: 247 L1 = a[cond1] 248 L2 = b[cond1] 249 L3 = c[cond1] 250 RFm = carlson_rf(L1 / L2, L1 / L3, np.ones(len(L1))) 251 RDm = carlson_rd(L1 / L2, L1 / L3, np.ones(len(L1))) 252 F1[cond1] = ((L1 + L2 + L3) ** 2) / (18 * (L1 - L2) * (L1 - L3)) * \ 253 (np.sqrt(L2 * L3) / L1 * RFm + 254 (3 * L1**2 - L1 * L2 - L1 * L3 - L2 * L3) / 255 (3 * L1 * np.sqrt(L2 * L3)) * RDm - 1) 256 257 # Resolve possible singularity a==b 258 cond2 = np.logical_and(cond0, np.logical_and(abs(a - b) < a * er, 259 abs(a - c) > a * er)) 260 if np.sum(cond2) != 0: 261 L1 = (a[cond2] + b[cond2]) / 2. 262 L3 = c[cond2] 263 F1[cond2] = _F2m(L3, L1, L1) / 2. 264 265 # Resolve possible singularity a==c 266 cond3 = np.logical_and(cond0, np.logical_and(abs(a - c) < a * er, 267 abs(a - b) > a * er)) 268 if np.sum(cond3) != 0: 269 L1 = (a[cond3] + c[cond3]) / 2. 270 L2 = b[cond3] 271 F1[cond3] = _F2m(L2, L1, L1) / 2 272 273 # Resolve possible singularity a==b and a==c 274 cond4 = np.logical_and(cond0, np.logical_and(abs(a - c) < a * er, 275 abs(a - b) < a * er)) 276 if np.sum(cond4) != 0: 277 F1[cond4] = 1 / 5. 278 279 return F1 280 281 282def _F2m(a, b, c): 283 r""" Helper function that computes function $F_2$ which is required to 284 compute the analytical solution of the Mean kurtosis. 285 286 Parameters 287 ---------- 288 a : ndarray 289 Array containing the values of parameter $\lambda_1$ of function $F_2$ 290 b : ndarray 291 Array containing the values of parameter $\lambda_2$ of function $F_2$ 292 c : ndarray 293 Array containing the values of parameter $\lambda_3$ of function $F_2$ 294 295 Returns 296 ------- 297 F2 : ndarray 298 Value of the function $F_2$ for all elements of the arrays a, b, and c 299 300 Notes 301 -------- 302 Function $F_2$ is defined as [1]_: 303 304 .. math:: 305 306 F_2(\lambda_1,\lambda_2,\lambda_3)= 307 \frac{(\lambda_1+\lambda_2+\lambda_3)^2} 308 {3(\lambda_2-\lambda_3)^2} 309 [\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}} 310 R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ 311 \frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}} 312 R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2] 313 314 References 315 ---------- 316 .. [1] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 317 Estimation of tensors and tensor-derived measures in diffusional 318 kurtosis imaging. Magn Reson Med. 65(3), 823-836 """ 319 # Eigenvalues are considered equal if they are not 2.5% different to each 320 # other. This value is adjusted according to the analysis reported in: 321 # http://gsoc2015dipydki.blogspot.co.uk/2015/08/rnh-post-13-start-wrapping-up-test.html 322 er = 2.5e-2 323 324 # Initialize F2 325 F2 = np.zeros(a.shape) 326 327 # Only computes F2 in voxels that have all eigenvalues larger than zero 328 cond0 = _positive_evals(a, b, c) 329 330 # Apply formula for non problematic plausible cases, i.e. b!=c 331 cond1 = np.logical_and(cond0, (abs(b - c) > b * er)) 332 if np.sum(cond1) != 0: 333 L1 = a[cond1] 334 L2 = b[cond1] 335 L3 = c[cond1] 336 RF = carlson_rf(L1 / L2, L1 / L3, np.ones(len(L1))) 337 RD = carlson_rd(L1 / L2, L1 / L3, np.ones(len(L1))) 338 F2[cond1] = (((L1 + L2 + L3) ** 2) / (3. * (L2 - L3) ** 2)) * \ 339 (((L2 + L3) / (np.sqrt(L2 * L3))) * RF + 340 ((2. * L1 - L2 - L3) / (3. * np.sqrt(L2 * L3))) * RD - 2.) 341 342 # Resolve possible singularity b==c 343 cond2 = np.logical_and(cond0, np.logical_and(abs(b - c) < b * er, 344 abs(a - b) > b * er)) 345 if np.sum(cond2) != 0: 346 L1 = a[cond2] 347 L3 = (c[cond2] + b[cond2]) / 2. 348 349 # Compute alfa [1]_ 350 x = 1. - (L1 / L3) 351 alpha = np.zeros(len(L1)) 352 for i in range(len(x)): 353 if x[i] > 0: 354 alpha[i] = 1. / np.sqrt(x[i]) * np.arctanh(np.sqrt(x[i])) 355 else: 356 alpha[i] = 1. / np.sqrt(-x[i]) * np.arctan(np.sqrt(-x[i])) 357 358 F2[cond2] = \ 359 6. * ((L1 + 2. * L3)**2) / (144. * L3**2 * (L1 - L3)**2) * \ 360 (L3 * (L1 + 2. * L3) + L1 * (L1 - 4. * L3) * alpha) 361 362 # Resolve possible singularity a==b and a==c 363 cond3 = np.logical_and(cond0, np.logical_and(abs(b - c) < b * er, 364 abs(a - b) < b * er)) 365 if np.sum(cond3) != 0: 366 F2[cond3] = 6 / 15. 367 368 return F2 369 370 371def directional_diffusion(dt, V, min_diffusivity=0): 372 r""" Calculates the apparent diffusion coefficient (adc) in each direction 373 of a sphere for a single voxel [1]_. 374 375 Parameters 376 ---------- 377 dt : array (6,) 378 elements of the diffusion tensor of the voxel. 379 V : array (g, 3) 380 g directions of a Sphere in Cartesian coordinates 381 min_diffusivity : float (optional) 382 Because negative eigenvalues are not physical and small eigenvalues 383 cause quite a lot of noise in diffusion-based metrics, diffusivity 384 values smaller than `min_diffusivity` are replaced with 385 `min_diffusivity`. Default = 0 386 387 Returns 388 -------- 389 adc : ndarray (g,) 390 Apparent diffusion coefficient (adc) in all g directions of a sphere 391 for a single voxel. 392 393 References 394 ---------- 395 .. [1] Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). 396 Exploring the 3D geometry of the diffusion kurtosis tensor - 397 Impact on the development of robust tractography procedures and 398 novel biomarkers, NeuroImage 111: 85-99 399 """ 400 adc = \ 401 V[:, 0] * V[:, 0] * dt[0] + \ 402 2 * V[:, 0] * V[:, 1] * dt[1] + \ 403 V[:, 1] * V[:, 1] * dt[2] + \ 404 2 * V[:, 0] * V[:, 2] * dt[3] + \ 405 2 * V[:, 1] * V[:, 2] * dt[4] + \ 406 V[:, 2] * V[:, 2] * dt[5] 407 408 if min_diffusivity is not None: 409 adc = adc.clip(min=min_diffusivity) 410 return adc 411 412 413def directional_diffusion_variance(kt, V, min_kurtosis=-3/7): 414 r""" Calculates the apparent diffusion variance (adv) in each direction 415 of a sphere for a single voxel [1]_. 416 417 Parameters 418 ---------- 419 dt : array (6,) 420 elements of the diffusion tensor of the voxel. 421 kt : array (15,) 422 elements of the kurtosis tensor of the voxel. 423 V : array (g, 3) 424 g directions of a Sphere in Cartesian coordinates 425 min_kurtosis : float (optional) 426 Because high-amplitude negative values of kurtosis are not physicaly 427 and biologicaly pluasible, and these cause artefacts in 428 kurtosis-based measures, directional kurtosis values smaller than 429 `min_kurtosis` are replaced with `min_kurtosis`. Default = -3./7 430 (theoretical kurtosis limit for regions that consist of water confined 431 to spherical pores [2]_) 432 adc : ndarray(g,) (optional) 433 Apparent diffusion coefficient (adc) in all g directions of a sphere 434 for a single voxel. 435 adv : ndarray(g,) (optional) 436 Apparent diffusion variance coefficient (advc) in all g directions of 437 a sphere for a single voxel. 438 439 Returns 440 -------- 441 adv : ndarray (g,) 442 Apparent diffusion variance (adv) in all g directions of a sphere for 443 a single voxel. 444 445 References 446 ---------- 447 .. [1] Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). 448 Exploring the 3D geometry of the diffusion kurtosis tensor - 449 Impact on the development of robust tractography procedures and 450 novel biomarkers, NeuroImage 111: 85-99 451 """ 452 adv = \ 453 V[:, 0] * V[:, 0] * V[:, 0] * V[:, 0] * kt[0] + \ 454 V[:, 1] * V[:, 1] * V[:, 1] * V[:, 1] * kt[1] + \ 455 V[:, 2] * V[:, 2] * V[:, 2] * V[:, 2] * kt[2] + \ 456 4 * V[:, 0] * V[:, 0] * V[:, 0] * V[:, 1] * kt[3] + \ 457 4 * V[:, 0] * V[:, 0] * V[:, 0] * V[:, 2] * kt[4] + \ 458 4 * V[:, 0] * V[:, 1] * V[:, 1] * V[:, 1] * kt[5] + \ 459 4 * V[:, 1] * V[:, 1] * V[:, 1] * V[:, 2] * kt[6] + \ 460 4 * V[:, 0] * V[:, 2] * V[:, 2] * V[:, 2] * kt[7] + \ 461 4 * V[:, 1] * V[:, 2] * V[:, 2] * V[:, 2] * kt[8] + \ 462 6 * V[:, 0] * V[:, 0] * V[:, 1] * V[:, 1] * kt[9] + \ 463 6 * V[:, 0] * V[:, 0] * V[:, 2] * V[:, 2] * kt[10] + \ 464 6 * V[:, 1] * V[:, 1] * V[:, 2] * V[:, 2] * kt[11] + \ 465 12 * V[:, 0] * V[:, 0] * V[:, 1] * V[:, 2] * kt[12] + \ 466 12 * V[:, 0] * V[:, 1] * V[:, 1] * V[:, 2] * kt[13] + \ 467 12 * V[:, 0] * V[:, 1] * V[:, 2] * V[:, 2] * kt[14] 468 469 return adv 470 471 472def directional_kurtosis(dt, md, kt, V, min_diffusivity=0, min_kurtosis=-3/7, 473 adc=None, adv=None): 474 r""" Calculates the apparent kurtosis coefficient (akc) in each direction 475 of a sphere for a single voxel [1]_. 476 477 Parameters 478 ---------- 479 dt : array (6,) 480 elements of the diffusion tensor of the voxel. 481 md : float 482 mean diffusivity of the voxel 483 kt : array (15,) 484 elements of the kurtosis tensor of the voxel. 485 V : array (g, 3) 486 g directions of a Sphere in Cartesian coordinates 487 min_diffusivity : float (optional) 488 Because negative eigenvalues are not physical and small eigenvalues 489 cause quite a lot of noise in diffusion-based metrics, diffusivity 490 values smaller than `min_diffusivity` are replaced with 491 `min_diffusivity`. Default = 0 492 min_kurtosis : float (optional) 493 Because high-amplitude negative values of kurtosis are not physicaly 494 and biologicaly pluasible, and these cause artefacts in 495 kurtosis-based measures, directional kurtosis values smaller than 496 `min_kurtosis` are replaced with `min_kurtosis`. Default = -3./7 497 (theoretical kurtosis limit for regions that consist of water confined 498 to spherical pores [2]_) 499 adc : ndarray(g,) (optional) 500 Apparent diffusion coefficient (adc) in all g directions of a sphere 501 for a single voxel. 502 adv : ndarray(g,) (optional) 503 Apparent diffusion variance (advc) in all g directions of a sphere for 504 a single voxel. 505 506 Returns 507 -------- 508 akc : ndarray (g,) 509 Apparent kurtosis coefficient (AKC) in all g directions of a sphere for 510 a single voxel. 511 512 References 513 ---------- 514 .. [1] Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). 515 Exploring the 3D geometry of the diffusion kurtosis tensor - 516 Impact on the development of robust tractography procedures and 517 novel biomarkers, NeuroImage 111: 85-99 518 .. [2] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 519 Robust estimation from DW-MRI using homogeneous polynomials. 520 Proceedings of the 8th {IEEE} International Symposium on 521 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 522 doi: 10.1109/ISBI.2011.5872402 523 """ 524 if adc is None: 525 adc = directional_diffusion(dt, V, min_diffusivity=min_diffusivity) 526 if adv is None: 527 adv = directional_diffusion_variance(kt, V) 528 529 akc = adv * (md / adc) ** 2 530 531 if min_kurtosis is not None: 532 akc = akc.clip(min=min_kurtosis) 533 534 return akc 535 536 537def apparent_kurtosis_coef(dki_params, sphere, min_diffusivity=0, 538 min_kurtosis=-3./7): 539 r""" Calculates the apparent kurtosis coefficient (AKC) in each direction 540 of a sphere [1]_. 541 542 Parameters 543 ---------- 544 dki_params : ndarray (x, y, z, 27) or (n, 27) 545 All parameters estimated from the diffusion kurtosis model. 546 Parameters are ordered as follows: 547 1) Three diffusion tensor's eigenvalues 548 2) Three lines of the eigenvector matrix each containing the first, 549 second and third coordinates of the eigenvectors respectively 550 3) Fifteen elements of the kurtosis tensor 551 sphere : a Sphere class instance 552 The AKC will be calculated for each of the vertices in the sphere 553 min_diffusivity : float (optional) 554 Because negative eigenvalues are not physical and small eigenvalues 555 cause quite a lot of noise in diffusion-based metrics, diffusivity 556 values smaller than `min_diffusivity` are replaced with 557 `min_diffusivity`. Default = 0 558 min_kurtosis : float (optional) 559 Because high-amplitude negative values of kurtosis are not physicaly 560 and biologicaly pluasible, and these cause artefacts in 561 kurtosis-based measures, directional kurtosis values smaller than 562 `min_kurtosis` are replaced with `min_kurtosis`. Default = -3./7 563 (theoretical kurtosis limit for regions that consist of water confined 564 to spherical pores [2]_) 565 566 Returns 567 -------- 568 akc : ndarray (x, y, z, g) or (n, g) 569 Apparent kurtosis coefficient (AKC) for all g directions of a sphere. 570 571 Notes 572 ----- 573 For each sphere direction with coordinates $(n_{1}, n_{2}, n_{3})$, the 574 calculation of AKC is done using formula [1]_: 575 576 .. math :: 577 578 AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} 579 \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl} 580 581 where $W_{ijkl}$ are the elements of the kurtosis tensor, MD the mean 582 diffusivity and ADC the apparent diffusion coefficent computed as: 583 584 .. math :: 585 586 ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij} 587 588 where $D_{ij}$ are the elements of the diffusion tensor. 589 590 References 591 ---------- 592 .. [1] Neto Henriques R, Correia MM, Nunes RG, Ferreira HA (2015). 593 Exploring the 3D geometry of the diffusion kurtosis tensor - 594 Impact on the development of robust tractography procedures and 595 novel biomarkers, NeuroImage 111: 85-99 596 .. [2] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 597 Robust estimation from DW-MRI using homogeneous polynomials. 598 Proceedings of the 8th {IEEE} International Symposium on 599 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 600 doi: 10.1109/ISBI.2011.5872402 601 """ 602 # Flat parameters 603 outshape = dki_params.shape[:-1] 604 dki_params = dki_params.reshape((-1, dki_params.shape[-1])) 605 606 # Split data 607 evals, evecs, kt = split_dki_param(dki_params) 608 609 # Initialize AKC matrix 610 V = sphere.vertices 611 akc = np.zeros((len(kt), len(V))) 612 613 # select relevant voxels to process 614 rel_i = _positive_evals(evals[..., 0], evals[..., 1], evals[..., 2]) 615 kt = kt[rel_i] 616 evecs = evecs[rel_i] 617 evals = evals[rel_i] 618 akci = akc[rel_i] 619 620 # Compute MD and DT 621 md = mean_diffusivity(evals) 622 dt = lower_triangular(vec_val_vect(evecs, evals)) 623 624 # loop over all relevant voxels 625 for vox in range(len(kt)): 626 akci[vox] = directional_kurtosis(dt[vox], md[vox], kt[vox], V, 627 min_diffusivity=min_diffusivity, 628 min_kurtosis=min_kurtosis) 629 630 # reshape data according to input data 631 akc[rel_i] = akci 632 633 return akc.reshape((outshape + (len(V),))) 634 635 636def mean_kurtosis(dki_params, min_kurtosis=-3./7, max_kurtosis=3, 637 analytical=True): 638 r""" Computes mean Kurtosis (MK) from the kurtosis tensor. 639 640 Parameters 641 ---------- 642 dki_params : ndarray (x, y, z, 27) or (n, 27) 643 All parameters estimated from the diffusion kurtosis model. 644 Parameters are ordered as follows: 645 1) Three diffusion tensor's eigenvalues 646 2) Three lines of the eigenvector matrix each containing the first, 647 second and third coordinates of the eigenvector 648 3) Fifteen elements of the kurtosis tensor 649 min_kurtosis : float (optional) 650 To keep kurtosis values within a plausible biophysical range, mean 651 kurtosis values that are smaller than `min_kurtosis` are replaced with 652 `min_kurtosis`. Default = -3./7 (theoretical kurtosis limit for regions 653 that consist of water confined to spherical pores [4]_) 654 max_kurtosis : float (optional) 655 To keep kurtosis values within a plausible biophysical range, mean 656 kurtosis values that are larger than `max_kurtosis` are replaced with 657 `max_kurtosis`. Default = 10 658 analytical : bool (optional) 659 If True, MK is calculated using its analytical solution, otherwise an 660 exact numerical estimator is used (see Notes). Default is set to True 661 662 Returns 663 ------- 664 mk : array 665 Calculated MK. 666 667 Notes 668 -------- 669 The MK is defined as the average of directional kurtosis coefficients 670 across all spatial directions, which can be formulated by the following 671 surface integral[1]_: 672 673 .. math:: 674 675 MK \equiv \frac{1}{4\pi} \int d\Omega_\mathbf{n} K(\mathbf{n}) 676 677 This integral can be numerically solved by averaging directional 678 kurtosis values sampled for directions of a spherical t-design [2]_. 679 680 Alternatively, MK can be solved from the analytical solution derived by 681 Tabesh et al. [3]_. This solution is given by: 682 683 .. math:: 684 685 MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+ 686 F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+ 687 F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\ 688 F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+ 689 F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+ 690 F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122} 691 692 where $\hat{W}_{ijkl}$ are the components of the $W$ tensor in the 693 coordinates system defined by the eigenvectors of the diffusion tensor 694 $\mathbf{D}$ and 695 696 .. math:: 697 698 F_1(\lambda_1,\lambda_2,\lambda_3)= 699 \frac{(\lambda_1+\lambda_2+\lambda_3)^2} 700 {18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)} 701 [\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1} 702 R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ 703 \frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3- 704 \lambda_1\lambda_3} 705 {3\lambda_1 \sqrt{\lambda_2 \lambda_3}} 706 R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ] 707 708 F_2(\lambda_1,\lambda_2,\lambda_3)= 709 \frac{(\lambda_1+\lambda_2+\lambda_3)^2} 710 {3(\lambda_2-\lambda_3)^2} 711 [\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}} 712 R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ 713 \frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}} 714 R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2] 715 716 where $R_f$ and $R_d$ are the Carlson's elliptic integrals. 717 718 References 719 ---------- 720 .. [1] Jensen, J.H., Helpern, J.A., 2010. MRI quantification of 721 non-Gaussian water diffusion by kurtosis analysis. NMR in 722 Biomedicine 23(7): 698-710 723 .. [2] Hardin, R.H., Sloane, N.J.A., 1996. McLaren's Improved Snub Cube and 724 Other New Spherical Designs in Three Dimensions. Discrete and 725 Computational Geometry 15, 429-441. 726 .. [3] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 727 Estimation of tensors and tensor-derived measures in diffusional 728 kurtosis imaging. Magn Reson Med. 65(3), 823-836 729 .. [4] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 730 Robust estimation from DW-MRI using homogeneous polynomials. 731 Proceedings of the 8th {IEEE} International Symposium on 732 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 733 doi: 10.1109/ISBI.2011.5872402 734 """ 735 # Flat parameters. For numpy versions more recent than 1.6.0, this step 736 # isn't required 737 outshape = dki_params.shape[:-1] 738 dki_params = dki_params.reshape((-1, dki_params.shape[-1])) 739 740 if analytical: 741 # Split the model parameters to three variable containing the evals, 742 # evecs, and kurtosis elements 743 evals, evecs, kt = split_dki_param(dki_params) 744 745 # Rotate the kurtosis tensor from the standard Cartesian coordinate 746 # system to another coordinate system in which the 3 orthonormal 747 # eigenvectors of DT are the base coordinate 748 Wxxxx = Wrotate_element(kt, 0, 0, 0, 0, evecs) 749 Wyyyy = Wrotate_element(kt, 1, 1, 1, 1, evecs) 750 Wzzzz = Wrotate_element(kt, 2, 2, 2, 2, evecs) 751 Wxxyy = Wrotate_element(kt, 0, 0, 1, 1, evecs) 752 Wxxzz = Wrotate_element(kt, 0, 0, 2, 2, evecs) 753 Wyyzz = Wrotate_element(kt, 1, 1, 2, 2, evecs) 754 755 # Compute MK 756 MK = \ 757 _F1m(evals[..., 0], evals[..., 1], evals[..., 2]) * Wxxxx + \ 758 _F1m(evals[..., 1], evals[..., 0], evals[..., 2]) * Wyyyy + \ 759 _F1m(evals[..., 2], evals[..., 1], evals[..., 0]) * Wzzzz + \ 760 _F2m(evals[..., 0], evals[..., 1], evals[..., 2]) * Wyyzz + \ 761 _F2m(evals[..., 1], evals[..., 0], evals[..., 2]) * Wxxzz + \ 762 _F2m(evals[..., 2], evals[..., 1], evals[..., 0]) * Wxxyy 763 764 else: 765 # Numerical Solution using t-design of 45 directions 766 V = np.loadtxt(get_fnames("t-design")) 767 sph = dps.Sphere(xyz=V) 768 KV = apparent_kurtosis_coef(dki_params, sph, min_kurtosis=min_kurtosis) 769 MK = np.mean(KV, axis=-1) 770 771 if min_kurtosis is not None: 772 MK = MK.clip(min=min_kurtosis) 773 774 if max_kurtosis is not None: 775 MK = MK.clip(max=max_kurtosis) 776 777 return MK.reshape(outshape) 778 779 780def _G1m(a, b, c): 781 r""" Helper function that computes function $G_1$ which is required to 782 compute the analytical solution of the Radial kurtosis. 783 784 Parameters 785 ---------- 786 a : ndarray 787 Array containing the values of parameter $\lambda_1$ of function $G_1$ 788 b : ndarray 789 Array containing the values of parameter $\lambda_2$ of function $G_1$ 790 c : ndarray 791 Array containing the values of parameter $\lambda_3$ of function $G_1$ 792 793 Returns 794 ------- 795 G1 : ndarray 796 Value of the function $G_1$ for all elements of the arrays a, b, and c 797 798 Notes 799 -------- 800 Function $G_1$ is defined as [1]_: 801 802 .. math:: 803 804 G_1(\lambda_1,\lambda_2,\lambda_3)= 805 \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2- 806 \lambda_3)} \left (2\lambda_2 + 807 \frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}} 808 \right) 809 810 References 811 ---------- 812 .. [1] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 813 Estimation of tensors and tensor-derived measures in diffusional 814 kurtosis imaging. Magn Reson Med. 65(3), 823-836 815 """ 816 # Float error used to compare two floats, abs(l1 - l2) < er for l1 = l2 817 # Error is defined as five orders of magnitude larger than system's epslon 818 er = np.finfo(a.ravel()[0]).eps * 1e5 819 820 # Initialize G1 821 G1 = np.zeros(a.shape) 822 823 # Only computes G1 in voxels that have all eigenvalues larger than zero 824 cond0 = _positive_evals(a, b, c) 825 826 # Apply formula for non problematic plausible cases, i.e. b!=c 827 cond1 = np.logical_and(cond0, (abs(b - c) > er)) 828 if np.sum(cond1) != 0: 829 L1 = a[cond1] 830 L2 = b[cond1] 831 L3 = c[cond1] 832 G1[cond1] = \ 833 (L1 + L2 + L3)**2 / (18 * L2 * (L2 - L3)**2) * \ 834 (2. * L2 + (L3**2 - 3 * L2 * L3) / np.sqrt(L2 * L3)) 835 836 # Resolve possible singularity b==c 837 cond2 = np.logical_and(cond0, abs(b - c) < er) 838 if np.sum(cond2) != 0: 839 L1 = a[cond2] 840 L2 = b[cond2] 841 G1[cond2] = (L1 + 2. * L2)**2 / (24. * L2**2) 842 843 return G1 844 845 846def _G2m(a, b, c): 847 r""" Helper function that computes function $G_2$ which is required to 848 compute the analytical solution of the Radial kurtosis. 849 850 Parameters 851 ---------- 852 a : ndarray 853 Array containing the values of parameter $\lambda_1$ of function $G_2$ 854 b : ndarray 855 Array containing the values of parameter $\lambda_2$ of function $G_2$ 856 c : ndarray (n,) 857 Array containing the values of parameter $\lambda_3$ of function $G_2$ 858 859 Returns 860 ------- 861 G2 : ndarray 862 Value of the function $G_2$ for all elements of the arrays a, b, and c 863 864 Notes 865 -------- 866 Function $G_2$ is defined as [1]_: 867 868 .. math:: 869 870 G_2(\lambda_1,\lambda_2,\lambda_3)= 871 \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2} 872 \left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-2\right ) 873 874 References 875 ---------- 876 .. [1] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 877 Estimation of tensors and tensor-derived measures in diffusional 878 kurtosis imaging. Magn Reson Med. 65(3), 823-836 879 """ 880 # Float error used to compare two floats, abs(l1 - l2) < er for l1 = l2 881 # Error is defined as five order of magnitude larger than system's epsilon 882 er = np.finfo(a.ravel()[0]).eps * 1e5 883 884 # Initialize G2 885 G2 = np.zeros(a.shape) 886 887 # Only computes G2 in voxels that have all eigenvalues larger than zero 888 cond0 = _positive_evals(a, b, c) 889 890 # Apply formula for non problematic plausible cases, i.e. b!=c 891 cond1 = np.logical_and(cond0, (abs(b - c) > er)) 892 if np.sum(cond1) != 0: 893 L1 = a[cond1] 894 L2 = b[cond1] 895 L3 = c[cond1] 896 G2[cond1] = \ 897 (L1 + L2 + L3)**2 / (3 * (L2 - L3)**2) * \ 898 ((L2 + L3) / np.sqrt(L2 * L3) - 2) 899 900 # Resolve possible singularity b==c 901 cond2 = np.logical_and(cond0, abs(b - c) < er) 902 if np.sum(cond2) != 0: 903 L1 = a[cond2] 904 L2 = b[cond2] 905 G2[cond2] = (L1 + 2. * L2)**2 / (12. * L2**2) 906 907 return G2 908 909 910def radial_kurtosis(dki_params, min_kurtosis=-3./7, max_kurtosis=10, 911 analytical=True): 912 r""" Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]_, [2]_. 913 914 Parameters 915 ---------- 916 dki_params : ndarray (x, y, z, 27) or (n, 27) 917 All parameters estimated from the diffusion kurtosis model. 918 Parameters are ordered as follows: 919 1) Three diffusion tensor's eigenvalues 920 2) Three lines of the eigenvector matrix each containing the first, 921 second and third coordinates of the eigenvector 922 3) Fifteen elements of the kurtosis tensor 923 min_kurtosis : float (optional) 924 To keep kurtosis values within a plausible biophysical range, radial 925 kurtosis values that are smaller than `min_kurtosis` are replaced with 926 `min_kurtosis`. Default = -3./7 (theoretical kurtosis limit for regions 927 that consist of water confined to spherical pores [3]_) 928 max_kurtosis : float (optional) 929 To keep kurtosis values within a plausible biophysical range, radial 930 kurtosis values that are larger than `max_kurtosis` are replaced with 931 `max_kurtosis`. Default = 10 932 analytical : bool (optional) 933 If True, RK is calculated using its analytical solution, otherwise an 934 exact numerical estimator is used (see Notes). Default is set to True. 935 936 Returns 937 ------- 938 rk : array 939 Calculated RK. 940 941 Notes 942 ----- 943 RK is defined as the average of the directional kurtosis perpendicular 944 to the fiber's main direction e1 [1]_, [2]_: 945 946 .. math:: 947 948 RK \equiv \frac{1}{2\pi} \int d\Omega _\mathbf{\theta} K(\mathbf{\theta}) 949 \delta (\mathbf{\theta}\cdot \mathbf{e}_1) 950 951 This equation can be numerically computed by averaging apparent 952 directional kurtosis samples for directions perpendicular to e1. 953 954 Otherwise, RK can be calculated from its analytical solution [2]_: 955 956 .. math:: 957 958 K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} + 959 G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} + 960 G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233} 961 962 where: 963 964 .. math:: 965 966 G_1(\lambda_1,\lambda_2,\lambda_3)= 967 \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2- 968 \lambda_3)} \left (2\lambda_2 + 969 \frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}} 970 \right) 971 972 and 973 974 .. math:: 975 976 G_2(\lambda_1,\lambda_2,\lambda_3)= 977 \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2} 978 \left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}-2\right ) 979 980 References 981 ---------- 982 .. [1] Jensen, J.H., Helpern, J.A., 2010. MRI quantification of 983 non-Gaussian water diffusion by kurtosis analysis. NMR in 984 Biomedicine 23(7): 698-710 985 .. [2] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 986 Estimation of tensors and tensor-derived measures in diffusional 987 kurtosis imaging. Magn Reson Med. 65(3), 823-836 988 .. [3] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 989 Robust estimation from DW-MRI using homogeneous polynomials. 990 Proceedings of the 8th {IEEE} International Symposium on Biomedical 991 Imaging: From Nano to Macro, ISBI 2011, 262-265. 992 doi: 10.1109/ISBI.2011.5872402 993 """ 994 outshape = dki_params.shape[:-1] 995 dki_params = dki_params.reshape((-1, dki_params.shape[-1])) 996 997 # Split the model parameters to three variable containing the evals, 998 # evecs, and kurtosis elements 999 evals, evecs, kt = split_dki_param(dki_params) 1000 1001 if analytical: 1002 # Rotate the kurtosis tensor from the standard Cartesian coordinate 1003 # system to another coordinate system in which the 3 orthonormal 1004 # eigenvectors of DT are the base coordinate 1005 Wyyyy = Wrotate_element(kt, 1, 1, 1, 1, evecs) 1006 Wzzzz = Wrotate_element(kt, 2, 2, 2, 2, evecs) 1007 Wyyzz = Wrotate_element(kt, 1, 1, 2, 2, evecs) 1008 1009 # Compute RK 1010 RK = \ 1011 _G1m(evals[..., 0], evals[..., 1], evals[..., 2]) * Wyyyy + \ 1012 _G1m(evals[..., 0], evals[..., 2], evals[..., 1]) * Wzzzz + \ 1013 _G2m(evals[..., 0], evals[..., 1], evals[..., 2]) * Wyyzz 1014 1015 else: 1016 # Numerical Solution using 10 perpendicular directions samples 1017 npa = 10 1018 1019 # Initialize RK 1020 RK = np.zeros(kt.shape[:-1]) 1021 1022 # select relevant voxels to process 1023 rel_i = _positive_evals(evals[..., 0], evals[..., 1], evals[..., 2]) 1024 dki_params = dki_params[rel_i] 1025 evecs = evecs[rel_i] 1026 RKi = RK[rel_i] 1027 1028 # loop over all voxels 1029 KV = np.zeros((dki_params.shape[0], npa)) 1030 for vox in range(len(dki_params)): 1031 V = perpendicular_directions(np.array(evecs[vox, :, 0]), num=npa, 1032 half=True) 1033 sph = dps.Sphere(xyz=V) 1034 KV[vox, :] = apparent_kurtosis_coef(dki_params[vox], sph, 1035 min_kurtosis=min_kurtosis) 1036 RKi = np.mean(KV, axis=-1) 1037 1038 RK[rel_i] = RKi 1039 1040 if min_kurtosis is not None: 1041 RK = RK.clip(min=min_kurtosis) 1042 1043 if max_kurtosis is not None: 1044 RK = RK.clip(max=max_kurtosis) 1045 1046 return RK.reshape(outshape) 1047 1048 1049def axial_kurtosis(dki_params, min_kurtosis=-3./7, max_kurtosis=10, 1050 analytical=True): 1051 r""" Computes axial Kurtosis (AK) from the kurtosis tensor [1]_, [2]_. 1052 1053 Parameters 1054 ---------- 1055 dki_params : ndarray (x, y, z, 27) or (n, 27) 1056 All parameters estimated from the diffusion kurtosis model. 1057 Parameters are ordered as follows: 1058 1) Three diffusion tensor's eigenvalues 1059 2) Three lines of the eigenvector matrix each containing the first, 1060 second and third coordinates of the eigenvector 1061 3) Fifteen elements of the kurtosis tensor 1062 min_kurtosis : float (optional) 1063 To keep kurtosis values within a plausible biophysical range, axial 1064 kurtosis values that are smaller than `min_kurtosis` are replaced with 1065 `min_kurtosis`. Default = -3./7 (theoretical kurtosis limit for regions 1066 that consist of water confined to spherical pores [3]_) 1067 max_kurtosis : float (optional) 1068 To keep kurtosis values within a plausible biophysical range, axial 1069 kurtosis values that are larger than `max_kurtosis` are replaced with 1070 `max_kurtosis`. Default = 10 1071 analytical : bool (optional) 1072 If True, AK is calculated from rotated diffusion kurtosis tensor, 1073 otherwise it will be computed from the apparent diffusion kurtosis 1074 values along the principal axis of the diffusion tensor (see notes). 1075 Default is set to True. 1076 1077 Returns 1078 ------- 1079 ak : array 1080 Calculated AK. 1081 1082 Notes 1083 ----- 1084 AK is defined as the directional kurtosis parallel to the fiber's main 1085 direction e1 [1]_, [2]_. You can compute AK using to approaches: 1086 1087 1) AK is calculated from rotated diffusion kurtosis tensor [2]_, i.e.: 1088 1089 .. math:: 1090 AK = \hat{W}_{1111} 1091 \frac{(\lambda_{1}+\lambda_{2}+\lambda_{3})^2}{(9 \lambda_{1}^2)} 1092 1093 2) AK can be sampled from the principal axis of the diffusion tensor: 1094 1095 .. math:: 1096 AK = K(\mathbf{\mathbf{e}_1) 1097 1098 Although both approaches leads to an exact calculation of AK, the first 1099 approach will be referred to as the analytical method while the second 1100 approach will be referred to as the numerical method based on their analogy 1101 to the estimation strategies for MK and RK. 1102 1103 References 1104 ---------- 1105 .. [1] Jensen, J.H., Helpern, J.A., 2010. MRI quantification of 1106 non-Gaussian water diffusion by kurtosis analysis. NMR in 1107 Biomedicine 23(7): 698-710 1108 .. [2] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 1109 Estimation of tensors and tensor-derived measures in diffusional 1110 kurtosis imaging. Magn Reson Med. 65(3), 823-836 1111 .. [3] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 1112 Robust estimation from DW-MRI using homogeneous polynomials. 1113 Proceedings of the 8th {IEEE} International Symposium on 1114 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 1115 doi: 10.1109/ISBI.2011.5872402 1116 """ 1117 # Flat parameters 1118 outshape = dki_params.shape[:-1] 1119 dki_params = dki_params.reshape((-1, dki_params.shape[-1])) 1120 1121 # Split data 1122 evals, evecs, kt = split_dki_param(dki_params) 1123 1124 # Initialize AK 1125 AK = np.zeros(kt.shape[:-1]) 1126 1127 # select relevant voxels to process 1128 rel_i = _positive_evals(evals[..., 0], evals[..., 1], evals[..., 2]) 1129 kt = kt[rel_i] 1130 evecs = evecs[rel_i] 1131 evals = evals[rel_i] 1132 AKi = AK[rel_i] 1133 1134 # Compute mean diffusivity 1135 md = mean_diffusivity(evals) 1136 1137 if analytical: 1138 # Rotate the kurtosis tensor from the standard Cartesian coordinate 1139 # system to another coordinate system in which the 3 orthonormal 1140 # eigenvectors of DT are the base coordinate 1141 Wxxxx = Wrotate_element(kt, 0, 0, 0, 0, evecs) 1142 AKi = Wxxxx * (md ** 2) / (evals[..., 0] ** 2) 1143 1144 else: 1145 # Compute apparent directional kurtosis along evecs[0] 1146 dt = lower_triangular(vec_val_vect(evecs, evals)) 1147 for vox in range(len(kt)): 1148 AKi[vox] = directional_kurtosis(dt[vox], md[vox], kt[vox], 1149 np.array([evecs[vox, :, 0]])) 1150 1151 # reshape data according to input data 1152 AK[rel_i] = AKi 1153 1154 if min_kurtosis is not None: 1155 AK = AK.clip(min=min_kurtosis) 1156 1157 if max_kurtosis is not None: 1158 AK = AK.clip(max=max_kurtosis) 1159 1160 return AK.reshape(outshape) 1161 1162 1163def _kt_maximum_converge(ang, dt, md, kt): 1164 """ Helper function that computes the inverse of the directional kurtosis 1165 of a voxel along a given direction in polar coordinates. 1166 1167 Parameters 1168 ---------- 1169 ang : array (2,) 1170 array containing the two polar angles 1171 dt : array (6,) 1172 elements of the diffusion tensor of the voxel. 1173 md : float 1174 mean diffusivity of the voxel 1175 kt : array (15,) 1176 elements of the kurtosis tensor of the voxel. 1177 1178 Returns 1179 ------- 1180 neg_kt : float 1181 The inverse value of the apparent kurtosis for the given direction. 1182 1183 Notes 1184 ----- 1185 This function is used to refine the kurtosis maximum estimate 1186 1187 See also 1188 -------- 1189 dipy.reconst.dki.kurtosis_maximum 1190 """ 1191 n = np.array([sphere2cart(1, ang[0], ang[1])]) 1192 return -1. * directional_kurtosis(dt, md, kt, n) 1193 1194 1195def _voxel_kurtosis_maximum(dt, md, kt, sphere, gtol=1e-2): 1196 """ Computes the maximum value of a single voxel kurtosis tensor 1197 1198 Parameters 1199 ---------- 1200 dt : array (6,) 1201 elements of the diffusion tensor of the voxel. 1202 md : float 1203 mean diffusivity of the voxel 1204 kt : array (15,) 1205 elements of the kurtosis tensor of the voxel. 1206 sphere : Sphere class instance, optional 1207 The sphere providing sample directions for the initial search of the 1208 maximum value of kurtosis. 1209 gtol : float, optional 1210 This input is to refine kurtosis maximum under the precision of the 1211 directions sampled on the sphere class instance. The gradient of the 1212 convergence procedure must be less than gtol before successful 1213 termination. If gtol is None, fiber direction is directly taken from 1214 the initial sampled directions of the given sphere object 1215 1216 Returns 1217 ------- 1218 max_value : float 1219 kurtosis tensor maximum value 1220 max_dir : array (3,) 1221 Cartesian coordinates of the direction of the maximal kurtosis value 1222 """ 1223 # Estimation of maximum kurtosis candidates 1224 akc = directional_kurtosis(dt, md, kt, sphere.vertices) 1225 max_val, ind = local_maxima(akc, sphere.edges) 1226 n = len(max_val) 1227 1228 # case that none maximum was find (spherical or null kurtosis tensors) 1229 if n == 0: 1230 return np.mean(akc), np.zeros(3) 1231 1232 max_dir = sphere.vertices[ind] 1233 1234 # Select the maximum from the candidates 1235 max_value = max(max_val) 1236 max_direction = max_dir[np.argmax(max_val.argmax)] 1237 1238 # refine maximum direction 1239 if gtol is not None: 1240 for p in range(n): 1241 r, theta, phi = cart2sphere(max_dir[p, 0], max_dir[p, 1], 1242 max_dir[p, 2]) 1243 ang = np.array([theta, phi]) 1244 ang[:] = opt.fmin_bfgs(_kt_maximum_converge, ang, 1245 args=(dt, md, kt), disp=False, 1246 retall=False, gtol=gtol) 1247 k_dir = np.array([sphere2cart(1., ang[0], ang[1])]) 1248 k_val = directional_kurtosis(dt, md, kt, k_dir) 1249 if k_val > max_value: 1250 max_value = k_val 1251 max_direction = k_dir 1252 1253 return max_value, max_direction 1254 1255 1256def kurtosis_maximum(dki_params, sphere='repulsion100', gtol=1e-2, 1257 mask=None): 1258 """ Computes kurtosis maximum value 1259 1260 Parameters 1261 ---------- 1262 dki_params : ndarray (x, y, z, 27) or (n, 27) 1263 All parameters estimated from the diffusion kurtosis model. 1264 Parameters are ordered as follows: 1265 1) Three diffusion tensor's eingenvalues 1266 2) Three lines of the eigenvector matrix each containing the first, 1267 second and third coordinates of the eigenvector 1268 3) Fifteen elements of the kurtosis tensor 1269 sphere : Sphere class instance, optional 1270 The sphere providing sample directions for the initial search of the 1271 maximal value of kurtosis. 1272 gtol : float, optional 1273 This input is to refine kurtosis maximum under the precision of the 1274 directions sampled on the sphere class instance. The gradient of the 1275 convergence procedure must be less than gtol before successful 1276 termination. If gtol is None, fiber direction is directly taken from 1277 the initial sampled directions of the given sphere object 1278 mask : ndarray 1279 A boolean array used to mark the coordinates in the data that should be 1280 analyzed that has the shape dki_params.shape[:-1] 1281 1282 Returns 1283 -------- 1284 max_value : float 1285 kurtosis tensor maximum value 1286 max_dir : array (3,) 1287 Cartesian coordinates of the direction of the maximal kurtosis value 1288 """ 1289 shape = dki_params.shape[:-1] 1290 1291 # load gradient directions 1292 if not isinstance(sphere, dps.Sphere): 1293 sphere = get_sphere('repulsion100') 1294 1295 # select voxels where to find fiber directions 1296 if mask is None: 1297 mask = np.ones(shape, dtype='bool') 1298 else: 1299 if mask.shape != shape: 1300 raise ValueError("Mask is not the same shape as dki_params.") 1301 1302 evals, evecs, kt = split_dki_param(dki_params) 1303 1304 # select non-zero voxels 1305 pos_evals = _positive_evals(evals[..., 0], evals[..., 1], evals[..., 2]) 1306 mask = np.logical_and(mask, pos_evals) 1307 1308 kt_max = np.zeros(mask.shape) 1309 dt = lower_triangular(vec_val_vect(evecs, evals)) 1310 md = mean_diffusivity(evals) 1311 1312 for idx in ndindex(shape): 1313 if not mask[idx]: 1314 continue 1315 kt_max[idx], da = _voxel_kurtosis_maximum(dt[idx], md[idx], kt[idx], 1316 sphere, gtol=gtol) 1317 1318 return kt_max 1319 1320 1321def mean_kurtosis_tensor(dki_params, min_kurtosis=-3./7, max_kurtosis=10): 1322 r""" Computes mean of the kurtosis tensor (MKT) [1]_. 1323 1324 Parameters 1325 ---------- 1326 dki_params : ndarray (x, y, z, 27) or (n, 27) 1327 All parameters estimated from the diffusion kurtosis model. 1328 Parameters are ordered as follows: 1329 1) Three diffusion tensor's eigenvalues 1330 2) Three lines of the eigenvector matrix each containing the first, 1331 second and third coordinates of the eigenvector 1332 3) Fifteen elements of the kurtosis tensor 1333 min_kurtosis : float (optional) 1334 To keep kurtosis values within a plausible biophysical range, mean 1335 kurtosis values that are smaller than `min_kurtosis` are replaced with 1336 `min_kurtosis`. Default = -3./7 (theoretical kurtosis limit for regions 1337 that consist of water confined to spherical pores [2]_) 1338 max_kurtosis : float (optional) 1339 To keep kurtosis values within a plausible biophysical range, mean 1340 kurtosis values that are larger than `max_kurtosis` are replaced with 1341 `max_kurtosis`. Default = 10 1342 Returns 1343 ------- 1344 mkt : array 1345 Calculated mean kurtosis tensor. 1346 1347 Notes 1348 -------- 1349 The MKT is defined as [1]_: 1350 1351 .. math:: 1352 1353 MKT \equiv \frac{1}{4\pi} \int d 1354 \Omega_{\mathnbf{n}} n_i n_j n_k n_l W_{ijkl} 1355 1356 1357 which can be directly computed from the trace of the kurtosis tensor: 1358 1359 .. math:: 1360 1361 MKT = \frac{1}{5} Tr(\mathbf{W}) = \frac{1}{5} 1362 (W_{1111} + W_{2222} + W_{3333} + 2W_{1122} + 2W_{1133} + 2W_{2233}) 1363 1364 References 1365 ---------- 1366 .. [1] Hansen, B., Lund, T. E., Sangill, R., and Jespersen, S. N. (2013). 1367 Experimentally and computationally fast method for estimation of 1368 a mean kurtosis.Magnetic Resonance in Medicine69, 1754–1760.388 1369 doi:10.1002/mrm.24743 1370 .. [2] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 1371 Robust estimation from DW-MRI using homogeneous polynomials. 1372 Proceedings of the 8th {IEEE} International Symposium on 1373 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 1374 doi: 10.1109/ISBI.2011.5872402 1375 """ 1376 MKT = 1/5 * (dki_params[..., 12] + dki_params[..., 13] + 1377 dki_params[..., 14] + 2 * dki_params[..., 21] + 1378 2 * dki_params[..., 22] + 2 * dki_params[..., 23]) 1379 1380 if min_kurtosis is not None: 1381 MKT = MKT.clip(min=min_kurtosis) 1382 1383 if max_kurtosis is not None: 1384 MKT = MKT.clip(max=max_kurtosis) 1385 1386 return MKT 1387 1388 1389def kurtosis_fractional_anisotropy(dki_params): 1390 r""" Computes the anisotropy of the kurtosis tensor (KFA) [1]_. 1391 1392 Parameters 1393 ---------- 1394 dki_params : ndarray (x, y, z, 27) or (n, 27) 1395 All parameters estimated from the diffusion kurtosis model. 1396 Parameters are ordered as follows: 1397 1) Three diffusion tensor's eigenvalues 1398 2) Three lines of the eigenvector matrix each containing the first, 1399 second and third coordinates of the eigenvector 1400 3) Fifteen elements of the kurtosis tensor 1401 Returns 1402 ------- 1403 kfa : array 1404 Calculated mean kurtosis tensor. 1405 1406 Notes 1407 -------- 1408 The KFA is defined as [1]_: 1409 1410 .. math:: 1411 1412 KFA \equiv 1413 \frac{||\mathbf{W} - MKT \mathbf{I}^{(4)}||_F}{||\mathbf{W}||_F} 1414 1415 where $W$ is the kurtosis tensor, MKT the kurtosis tensor mean, $I^(4)$ is 1416 the fully symmetric rank 2 isotropic tensor and $||...||_F$ is the tensor's 1417 Frobenius norm [1]_. 1418 1419 References 1420 ---------- 1421 .. [1] Glenn, G. R., Helpern, J. A., Tabesh, A., and Jensen, J. H. (2015). 1422 Quantitative assessment of diffusional kurtosis anisotropy. 1423 NMR in Biomedicine 28, 448–459. doi:10.1002/nbm.3271 1424 """ 1425 Wxxxx = dki_params[..., 12] 1426 Wyyyy = dki_params[..., 13] 1427 Wzzzz = dki_params[..., 14] 1428 Wxxxy = dki_params[..., 15] 1429 Wxxxz = dki_params[..., 16] 1430 Wxyyy = dki_params[..., 17] 1431 Wyyyz = dki_params[..., 18] 1432 Wxzzz = dki_params[..., 19] 1433 Wyzzz = dki_params[..., 20] 1434 Wxxyy = dki_params[..., 21] 1435 Wxxzz = dki_params[..., 22] 1436 Wyyzz = dki_params[..., 23] 1437 Wxxyz = dki_params[..., 24] 1438 Wxyyz = dki_params[..., 25] 1439 Wxyzz = dki_params[..., 26] 1440 1441 W = 1.0/5.0 * (Wxxxx + Wyyyy + Wzzzz + 2*Wxxyy + 2*Wxxzz + 2*Wyyzz) 1442 1443 # Compute's equation numerator 1444 A = (Wxxxx - W) ** 2 + (Wyyyy - W) ** 2 + (Wzzzz - W) ** 2 + \ 1445 4 * Wxxxy ** 2 + 4 * Wxxxz ** 2 + 4 * Wxyyy ** 2 + 4 * Wyyyz ** 2 + \ 1446 4 * Wxzzz ** 2 + 4 * Wyzzz ** 2 + \ 1447 6 * (Wxxyy - W/3) ** 2 + 6 * (Wxxzz - W/3) ** 2 + \ 1448 6 * (Wyyzz - W/3) ** 2 + \ 1449 12 * Wxxyz ** 2 + 12 * Wxyyz ** 2 + 12 * Wxyzz ** 2 1450 1451 # Compute's equation denominator 1452 B = Wxxxx ** 2 + Wyyyy ** 2 + Wzzzz ** 2 + 4 * Wxxxy ** 2 + \ 1453 4 * Wxxxz ** 2 + 4 * Wxyyy ** 2 + 4 * Wyyyz ** 2 + 4 * Wxzzz ** 2 + \ 1454 4 * Wyzzz ** 2 + 6 * Wxxyy ** 2 + 6 * Wxxzz ** 2 + 6 * Wyyzz ** 2 + \ 1455 12 * Wxxyz ** 2 + 12 * Wxyyz ** 2 + 12 * Wxyzz ** 2 1456 1457 # Compute KFA 1458 KFA = np.zeros(A.shape) 1459 cond = B > 0 # Avoiding Singularity (if B = 0, KFA = 0) 1460 KFA[cond] = np.sqrt(A[cond]/B[cond]) 1461 1462 return KFA 1463 1464 1465def dki_prediction(dki_params, gtab, S0=1.): 1466 """ Predict a signal given diffusion kurtosis imaging parameters. 1467 1468 Parameters 1469 ---------- 1470 dki_params : ndarray (x, y, z, 27) or (n, 27) 1471 All parameters estimated from the diffusion kurtosis model. 1472 Parameters are ordered as follows: 1473 1) Three diffusion tensor's eigenvalues 1474 2) Three lines of the eigenvector matrix each containing the first, 1475 second and third coordinates of the eigenvector 1476 3) Fifteen elements of the kurtosis tensor 1477 gtab : a GradientTable class instance 1478 The gradient table for this prediction 1479 S0 : float or ndarray (optional) 1480 The non diffusion-weighted signal in every voxel, or across all 1481 voxels. Default: 1 1482 1483 Returns 1484 -------- 1485 S : (..., N) ndarray 1486 Simulated signal based on the DKI model: 1487 1488 .. math:: 1489 1490 S=S_{0}e^{-bD+\frac{1}{6}b^{2}D^{2}K} 1491 1492 """ 1493 evals, evecs, kt = split_dki_param(dki_params) 1494 1495 # Define DKI design matrix according to given gtab 1496 A = design_matrix(gtab) 1497 1498 # Flat parameters and initialize pred_sig 1499 fevals = evals.reshape((-1, evals.shape[-1])) 1500 fevecs = evecs.reshape((-1,) + evecs.shape[-2:]) 1501 fkt = kt.reshape((-1, kt.shape[-1])) 1502 pred_sig = np.zeros((len(fevals), len(gtab.bvals))) 1503 if isinstance(S0, np.ndarray): 1504 S0_vol = np.reshape(S0, (len(fevals))) 1505 else: 1506 S0_vol = S0 1507 # looping for all voxels 1508 for v in range(len(pred_sig)): 1509 DT = np.dot(np.dot(fevecs[v], np.diag(fevals[v])), fevecs[v].T) 1510 dt = lower_triangular(DT) 1511 MD = (dt[0] + dt[2] + dt[5]) / 3 1512 if isinstance(S0_vol, np.ndarray): 1513 this_S0 = S0_vol[v] 1514 else: 1515 this_S0 = S0_vol 1516 X = np.concatenate((dt, fkt[v] * MD * MD, 1517 np.array([-np.log(this_S0)])), 1518 axis=0) 1519 pred_sig[v] = np.exp(np.dot(A, X)) 1520 1521 # Reshape data according to the shape of dki_params 1522 pred_sig = pred_sig.reshape(dki_params.shape[:-1] + (pred_sig.shape[-1],)) 1523 1524 return pred_sig 1525 1526 1527class DiffusionKurtosisModel(ReconstModel): 1528 """ Class for the Diffusion Kurtosis Model 1529 """ 1530 1531 def __init__(self, gtab, fit_method="WLS", *args, **kwargs): 1532 """ Diffusion Kurtosis Tensor Model [1] 1533 1534 Parameters 1535 ---------- 1536 gtab : GradientTable class instance 1537 1538 fit_method : str or callable 1539 str can be one of the following: 1540 'OLS' or 'ULLS' for ordinary least squares 1541 dki.ols_fit_dki 1542 'WLS' or 'UWLLS' for weighted ordinary least squares 1543 dki.wls_fit_dki 1544 1545 callable has to have the signature: 1546 fit_method(design_matrix, data, *args, **kwargs) 1547 1548 args, kwargs : arguments and key-word arguments passed to the 1549 fit_method. See dki.ols_fit_dki, dki.wls_fit_dki for details 1550 1551 References 1552 ---------- 1553 .. [1] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 1554 Estimation of tensors and tensor-derived measures in diffusional 1555 kurtosis imaging. Magn Reson Med. 65(3), 823-836 1556 """ 1557 ReconstModel.__init__(self, gtab) 1558 1559 if not callable(fit_method): 1560 try: 1561 self.fit_method = common_fit_methods[fit_method] 1562 except KeyError: 1563 raise ValueError('"' + str(fit_method) + '" is not a known ' 1564 'fit method, the fit method should either be ' 1565 'a function or one of the common fit methods') 1566 1567 self.design_matrix = design_matrix(self.gtab) 1568 self.args = args 1569 self.kwargs = kwargs 1570 self.min_signal = self.kwargs.pop('min_signal', None) 1571 if self.min_signal is not None and self.min_signal <= 0: 1572 e_s = "The `min_signal` key-word argument needs to be strictly" 1573 e_s += " positive." 1574 raise ValueError(e_s) 1575 1576 # Check if at least three b-values are given 1577 enough_b = check_multi_b(self.gtab, 3, non_zero=False) 1578 if not enough_b: 1579 mes = "DKI requires at least 3 b-values (which can include b=0)" 1580 raise ValueError(mes) 1581 1582 def fit(self, data, mask=None): 1583 """ Fit method of the DKI model class 1584 1585 Parameters 1586 ---------- 1587 data : array 1588 The measured signal from one voxel. 1589 1590 mask : array 1591 A boolean array used to mark the coordinates in the data that 1592 should be analyzed that has the shape data.shape[-1] 1593 """ 1594 if mask is not None: 1595 # Check for valid shape of the mask 1596 if mask.shape != data.shape[:-1]: 1597 raise ValueError("Mask is not the same shape as data.") 1598 mask = np.array(mask, dtype=bool, copy=False) 1599 data_in_mask = np.reshape(data[mask], (-1, data.shape[-1])) 1600 1601 if self.min_signal is None: 1602 min_signal = MIN_POSITIVE_SIGNAL 1603 else: 1604 min_signal = self.min_signal 1605 1606 data_in_mask = np.maximum(data_in_mask, min_signal) 1607 params_in_mask = self.fit_method(self.design_matrix, data_in_mask, 1608 *self.args, **self.kwargs) 1609 1610 if mask is None: 1611 out_shape = data.shape[:-1] + (-1, ) 1612 dki_params = params_in_mask.reshape(out_shape) 1613 else: 1614 dki_params = np.zeros(data.shape[:-1] + (27,)) 1615 dki_params[mask, :] = params_in_mask 1616 1617 return DiffusionKurtosisFit(self, dki_params) 1618 1619 def predict(self, dki_params, S0=1.): 1620 """ Predict a signal for this DKI model class instance given 1621 parameters. 1622 1623 Parameters 1624 ---------- 1625 dki_params : ndarray (x, y, z, 27) or (n, 27) 1626 All parameters estimated from the diffusion kurtosis model. 1627 Parameters are ordered as follows: 1628 1) Three diffusion tensor's eigenvalues 1629 2) Three lines of the eigenvector matrix each containing the 1630 first, second and third coordinates of the eigenvector 1631 3) Fifteen elements of the kurtosis tensor 1632 S0 : float or ndarray (optional) 1633 The non diffusion-weighted signal in every voxel, or across all 1634 voxels. Default: 1 1635 """ 1636 return dki_prediction(dki_params, self.gtab, S0) 1637 1638 1639class DiffusionKurtosisFit(TensorFit): 1640 """ Class for fitting the Diffusion Kurtosis Model""" 1641 1642 def __init__(self, model, model_params): 1643 """ Initialize a DiffusionKurtosisFit class instance. 1644 1645 Since DKI is an extension of DTI, class instance is defined as subclass 1646 of the TensorFit from dti.py 1647 1648 Parameters 1649 ---------- 1650 model : DiffusionKurtosisModel Class instance 1651 Class instance containing the Diffusion Kurtosis Model for the fit 1652 model_params : ndarray (x, y, z, 27) or (n, 27) 1653 All parameters estimated from the diffusion kurtosis model. 1654 Parameters are ordered as follows: 1655 1) Three diffusion tensor's eigenvalues 1656 2) Three lines of the eigenvector matrix each containing the 1657 first, second and third coordinates of the eigenvector 1658 3) Fifteen elements of the kurtosis tensor 1659 """ 1660 TensorFit.__init__(self, model, model_params) 1661 1662 @property 1663 def kt(self): 1664 """ 1665 Returns the 15 independent elements of the kurtosis tensor as an array 1666 """ 1667 return self.model_params[..., 12:] 1668 1669 def akc(self, sphere): 1670 r""" Calculates the apparent kurtosis coefficient (AKC) in each 1671 direction on the sphere for each voxel in the data 1672 1673 Parameters 1674 ---------- 1675 sphere : Sphere class instance 1676 1677 Returns 1678 ------- 1679 akc : ndarray 1680 The estimates of the apparent kurtosis coefficient in every 1681 direction on the input sphere 1682 1683 Notes 1684 ----- 1685 For each sphere direction with coordinates $(n_{1}, n_{2}, n_{3})$, the 1686 calculation of AKC is done using formula: 1687 1688 .. math :: 1689 1690 AKC(n)=\frac{MD^{2}}{ADC(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} 1691 \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl} 1692 1693 where $W_{ijkl}$ are the elements of the kurtosis tensor, MD the mean 1694 diffusivity and ADC the apparent diffusion coefficent computed as: 1695 1696 .. math :: 1697 1698 ADC(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij} 1699 1700 where $D_{ij}$ are the elements of the diffusion tensor. 1701 """ 1702 return apparent_kurtosis_coef(self.model_params, sphere) 1703 1704 def mk(self, min_kurtosis=-3./7, max_kurtosis=10, analytical=True): 1705 r""" Computes mean Kurtosis (MK) from the kurtosis tensor. 1706 1707 Parameters 1708 ---------- 1709 min_kurtosis : float (optional) 1710 To keep kurtosis values within a plausible biophysical range, mean 1711 kurtosis values that are smaller than `min_kurtosis` are replaced 1712 with `min_kurtosis`. Default = -3./7 (theoretical kurtosis limit 1713 for regions that consist of water confined to spherical pores [4]_) 1714 max_kurtosis : float (optional) 1715 To keep kurtosis values within a plausible biophysical range, mean 1716 kurtosis values that are larger than `max_kurtosis` are replaced 1717 with `max_kurtosis`. Default = 10 1718 analytical : bool (optional) 1719 If True, MK is calculated using its analytical solution, otherwise 1720 an exact numerical estimator is used (see Notes). Default is set to 1721 True. 1722 1723 Returns 1724 ------- 1725 mk : array 1726 Calculated MK. 1727 1728 Notes 1729 -------- 1730 The MK is defined as the average of directional kurtosis coefficients 1731 across all spatial directions, which can be formulated by the following 1732 surface integral[1]_: 1733 1734 .. math:: 1735 1736 MK \equiv \frac{1}{4\pi} \int d\Omega_\mathbf{n} K(\mathbf{n}) 1737 1738 This integral can be numerically solved by averaging directional 1739 kurtosis values sampled for directions of a spherical t-design [2]_. 1740 1741 Alternatively, MK can be solved from the analytical solution derived by 1742 Tabesh et al. [3]_. This solution is given by: 1743 1744 .. math:: 1745 1746 MK=F_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{1111}+ 1747 F_1(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{2222}+ 1748 F_1(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{3333}+ \\ 1749 F_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233}+ 1750 F_2(\lambda_2,\lambda_1,\lambda_3)\hat{W}_{1133}+ 1751 F_2(\lambda_3,\lambda_2,\lambda_1)\hat{W}_{1122} 1752 1753 where $\hat{W}_{ijkl}$ are the components of the $W$ tensor in the 1754 coordinates system defined by the eigenvectors of the diffusion tensor 1755 $\mathbf{D}$ and 1756 1757 .. math:: 1758 1759 F_1(\lambda_1,\lambda_2,\lambda_3)= 1760 \frac{(\lambda_1+\lambda_2+\lambda_3)^2} 1761 {18(\lambda_1-\lambda_2)(\lambda_1-\lambda_3)} 1762 [\frac{\sqrt{\lambda_2\lambda_3}}{\lambda_1} 1763 R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ 1764 \frac{3\lambda_1^2-\lambda_1\lambda_2-\lambda_2\lambda_3- 1765 \lambda_1\lambda_3} 1766 {3\lambda_1 \sqrt{\lambda_2 \lambda_3}} 1767 R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-1 ] 1768 1769 F_2(\lambda_1,\lambda_2,\lambda_3)= 1770 \frac{(\lambda_1+\lambda_2+\lambda_3)^2} 1771 {3(\lambda_2-\lambda_3)^2} 1772 [\frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}} 1773 R_F(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)+\\ 1774 \frac{2\lambda_1-\lambda_2-\lambda_3}{3\sqrt{\lambda_2 \lambda_3}} 1775 R_D(\frac{\lambda_1}{\lambda_2},\frac{\lambda_1}{\lambda_3},1)-2] 1776 1777 where $R_f$ and $R_d$ are the Carlson's elliptic integrals. 1778 1779 References 1780 ---------- 1781 .. [1] Jensen, J.H., Helpern, J.A., 2010. MRI quantification of 1782 non-Gaussian water diffusion by kurtosis analysis. NMR in 1783 Biomedicine 23(7): 698-710 1784 .. [2] Hardin, R.H., Sloane, N.J.A., 1996. McLaren's Improved Snub Cube 1785 and Other New Spherical Designs in Three Dimensions. Discrete 1786 and Computational Geometry 15, 429-441. 1787 .. [3] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 1788 Estimation of tensors and tensor-derived measures in diffusional 1789 kurtosis imaging. Magn Reson Med. 65(3), 823-836 1790 .. [4] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 1791 Robust estimation from DW-MRI using homogeneous polynomials. 1792 Proceedings of the 8th {IEEE} International Symposium on 1793 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 1794 doi: 10.1109/ISBI.2011.5872402 1795 """ 1796 return mean_kurtosis(self.model_params, min_kurtosis, max_kurtosis, 1797 analytical) 1798 1799 def ak(self, min_kurtosis=-3./7, max_kurtosis=10, analytical=True): 1800 r""" 1801 Axial Kurtosis (AK) of a diffusion kurtosis tensor [1]_. 1802 1803 Parameters 1804 ---------- 1805 min_kurtosis : float (optional) 1806 To keep kurtosis values within a plausible biophysical range, axial 1807 kurtosis values that are smaller than `min_kurtosis` are replaced 1808 with -3./7 (theoretical kurtosis limit 1809 for regions that consist of water confined to spherical pores [2]_) 1810 max_kurtosis : float (optional) 1811 To keep kurtosis values within a plausible biophysical range, axial 1812 kurtosis values that are larger than `max_kurtosis` are replaced 1813 with `max_kurtosis`. Default = 10 1814 analytical : bool (optional) 1815 If True, AK is calculated from rotated diffusion kurtosis tensor, 1816 otherwise it will be computed from the apparent diffusion kurtosis 1817 values along the principal axis of the diffusion tensor 1818 (see notes). Default is set to True. 1819 1820 Returns 1821 ------- 1822 ak : array 1823 Calculated AK. 1824 1825 Notes 1826 ----- 1827 AK is defined as the directional kurtosis parallel to the fiber's main 1828 direction e1 [1]_, [2]_. You can compute AK using to approaches: 1829 1830 1) AK is calculated from rotated diffusion kurtosis tensor [2]_, i.e.: 1831 1832 .. math:: 1833 AK = \hat{W}_{1111} 1834 \frac{(\lambda_{1}+\lambda_{2}+\lambda_{3})^2}{(9 \lambda_{1}^2)} 1835 1836 2) AK can be sampled from the principal axis of the diffusion tensor: 1837 1838 .. math:: 1839 AK = K(\mathbf{\mathbf{e}_1) 1840 1841 Although both approaches leads to an exact calculation of AK, the 1842 first approach will be referred to as the analytical method while the 1843 second approach will be referred to as the numerical method based on 1844 their analogy to the estimation strategies for MK and RK. 1845 1846 References 1847 ---------- 1848 .. [1] Jensen, J.H., Helpern, J.A., 2010. MRI quantification of 1849 non-Gaussian water diffusion by kurtosis analysis. NMR in 1850 Biomedicine 23(7): 698-710 1851 .. [2] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 1852 Estimation of tensors and tensor-derived measures in diffusional 1853 kurtosis imaging. Magn Reson Med. 65(3), 823-836 1854 .. [3] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 1855 Robust estimation from DW-MRI using homogeneous polynomials. 1856 Proceedings of the 8th {IEEE} International Symposium on 1857 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 1858 doi: 10.1109/ISBI.2011.5872402 1859 """ 1860 return axial_kurtosis(self.model_params, min_kurtosis, max_kurtosis, 1861 analytical) 1862 1863 def rk(self, min_kurtosis=-3./7, max_kurtosis=10, analytical=True): 1864 r""" Radial Kurtosis (RK) of a diffusion kurtosis tensor [1]_. 1865 1866 Parameters 1867 ---------- 1868 min_kurtosis : float (optional) 1869 To keep kurtosis values within a plausible biophysical range, 1870 radial kurtosis values that are smaller than `min_kurtosis` are 1871 replaced with `min_kurtosis`. Default = -3./7 (theoretical kurtosis 1872 limit for regions that consist of water confined to spherical pores 1873 [3]_) 1874 max_kurtosis : float (optional) 1875 To keep kurtosis values within a plausible biophysical range, 1876 radial kurtosis values that are larger than `max_kurtosis` are 1877 replaced with `max_kurtosis`. Default = 10 1878 analytical : bool (optional) 1879 If True, RK is calculated using its analytical solution, otherwise 1880 an exact numerical estimator is used (see Notes). Default is set to 1881 True 1882 1883 Returns 1884 ------- 1885 rk : array 1886 Calculated RK. 1887 1888 Notes 1889 -------- 1890 RK is defined as the average of the directional kurtosis perpendicular 1891 to the fiber's main direction e1 [1]_, [2]_: 1892 1893 .. math:: 1894 1895 RK \equiv \frac{1}{2\pi} \int d\Omega _\mathbf{\theta} 1896 K(\mathbf{\theta}) \delta (\mathbf{\theta}\cdot \mathbf{e}_1) 1897 1898 This equation can be numerically computed by averaging apparent 1899 directional kurtosis samples for directions perpendicular to e1. 1900 1901 Otherwise, RK can be calculated from its analytical solution [2]_: 1902 1903 .. math:: 1904 1905 K_{\bot} = G_1(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2222} + 1906 G_1(\lambda_1,\lambda_3,\lambda_2)\hat{W}_{3333} + 1907 G_2(\lambda_1,\lambda_2,\lambda_3)\hat{W}_{2233} 1908 1909 where: 1910 1911 .. math:: 1912 1913 G_1(\lambda_1,\lambda_2,\lambda_3)= 1914 \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{18\lambda_2(\lambda_2- 1915 \lambda_3)} \left (2\lambda_2 + 1916 \frac{\lambda_3^2-3\lambda_2\lambda_3}{\sqrt{\lambda_2\lambda_3}} 1917 \right) 1918 1919 and 1920 1921 .. math:: 1922 1923 G_2(\lambda_1,\lambda_2,\lambda_3)= 1924 \frac{(\lambda_1+\lambda_2+\lambda_3)^2}{(\lambda_2-\lambda_3)^2} 1925 \left ( \frac{\lambda_2+\lambda_3}{\sqrt{\lambda_2\lambda_3}}- 1926 2\right ) 1927 1928 References 1929 ---------- 1930 .. [1] Jensen, J.H., Helpern, J.A., 2010. MRI quantification of 1931 non-Gaussian water diffusion by kurtosis analysis. NMR in 1932 Biomedicine 23(7): 698-710 1933 .. [2] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 1934 Estimation of tensors and tensor-derived measures in diffusional 1935 kurtosis imaging. Magn Reson Med. 65(3), 823-836 1936 .. [3] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 1937 Robust estimation from DW-MRI using homogeneous polynomials. 1938 Proceedings of the 8th {IEEE} International Symposium on 1939 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 1940 doi: 10.1109/ISBI.2011.5872402 1941 """ 1942 return radial_kurtosis(self.model_params, min_kurtosis, max_kurtosis, 1943 analytical) 1944 1945 def kmax(self, sphere='repulsion100', gtol=1e-5, mask=None): 1946 r""" Computes the maximum value of a single voxel kurtosis tensor 1947 1948 Parameters 1949 ---------- 1950 sphere : Sphere class instance, optional 1951 The sphere providing sample directions for the initial search of 1952 the maximum value of kurtosis. 1953 gtol : float, optional 1954 This input is to refine kurtosis maximum under the precision of the 1955 directions sampled on the sphere class instance. The gradient of 1956 the convergence procedure must be less than gtol before successful 1957 termination. If gtol is None, fiber direction is directly taken 1958 from the initial sampled directions of the given sphere object 1959 1960 Returns 1961 -------- 1962 max_value : float 1963 kurtosis tensor maximum value 1964 """ 1965 return kurtosis_maximum(self.model_params, sphere, gtol, mask) 1966 1967 def mkt(self, min_kurtosis=-3./7, max_kurtosis=10): 1968 r""" Computes mean of the kurtosis tensor (MKT) [1]_. 1969 1970 Parameters 1971 ---------- 1972 min_kurtosis : float (optional) 1973 To keep kurtosis values within a plausible biophysical range, mean 1974 kurtosis values that are smaller than `min_kurtosis` are replaced 1975 with `min_kurtosis`. Default = -3./7 (theoretical kurtosis limit 1976 for regions that consist of water confined to spherical pores [2]_) 1977 max_kurtosis : float (optional) 1978 To keep kurtosis values within a plausible biophysical range, mean 1979 kurtosis values that are larger than `max_kurtosis` are replaced 1980 with `max_kurtosis`. Default = 10 1981 1982 Returns 1983 ------- 1984 mkt : array 1985 Calculated mean kurtosis tensor. 1986 1987 Notes 1988 -------- 1989 The MKT is defined as [1]_: 1990 1991 .. math:: 1992 1993 MKT \equiv \frac{1}{4\pi} \int d 1994 \Omega_{\mathnbf{n}} n_i n_j n_k n_l W_{ijkl} 1995 1996 1997 which can be directly computed from the trace of the kurtosis tensor: 1998 1999 .. math:: 2000 2001 MKT = \frac{1}{5} Tr(\mathbf{W}) = \frac{1}{5} 2002 (W_{1111} + W_{2222} + W_{3333} + 2W_{1122} + 2W_{1133} + 2W_{2233}) 2003 2004 References 2005 ---------- 2006 .. [1] Hansen, B., Lund, T. E., Sangill, R., and Jespersen, S. N. 2013. 2007 Experimentally and computationally fast method for estimation 2008 of a mean kurtosis. Magnetic Resonance in Medicine69, 1754–1760. 2009 388. doi:10.1002/mrm.24743 2010 .. [2] Barmpoutis, A., & Zhuo, J., 2011. Diffusion kurtosis imaging: 2011 Robust estimation from DW-MRI using homogeneous polynomials. 2012 Proceedings of the 8th {IEEE} International Symposium on 2013 Biomedical Imaging: From Nano to Macro, ISBI 2011, 262-265. 2014 doi: 10.1109/ISBI.2011.5872402 2015 """ 2016 return mean_kurtosis_tensor(self.model_params, min_kurtosis, 2017 max_kurtosis) 2018 2019 @property 2020 def kfa(self): 2021 r""" Returns the kurtosis tensor (KFA) [1]_. 2022 2023 Notes 2024 -------- 2025 The KFA is defined as [1]_: 2026 2027 .. math:: 2028 2029 KFA \equiv 2030 \frac{||\mathbf{W} - MKT \mathbf{I}^{(4)}||_F}{||\mathbf{W}||_F} 2031 2032 where $W$ is the kurtosis tensor, MKT the kurtosis tensor mean, $I^(4)$ 2033 is the fully symmetric rank 2 isotropic tensor and $||...||_F$ is the 2034 tensor's Frobenius norm [1]_. 2035 2036 References 2037 ---------- 2038 .. [1] Glenn, G. R., Helpern, J. A., Tabesh, A., and Jensen, J. H. 2039 (2015). Quantitative assessment of diffusional kurtosis 2040 anisotropy. NMR in Biomedicine 28, 448–459. doi:10.1002/nbm.3271 2041 """ 2042 return kurtosis_fractional_anisotropy(self.model_params) 2043 2044 def predict(self, gtab, S0=1.): 2045 r""" Given a DKI model fit, predict the signal on the vertices of a 2046 gradient table 2047 2048 Parameters 2049 ---------- 2050 gtab : a GradientTable class instance 2051 The gradient table for this prediction 2052 2053 S0 : float or ndarray (optional) 2054 The non diffusion-weighted signal in every voxel, or across all 2055 voxels. Default: 1 2056 2057 Notes 2058 ----- 2059 The predicted signal is given by: 2060 2061 .. math:: 2062 2063 S(n,b)=S_{0}e^{-bD(n)+\frac{1}{6}b^{2}D(n)^{2}K(n)} 2064 2065 $\mathbf{D(n)}$ and $\mathbf{K(n)}$ can be computed from the DT and KT 2066 using the following equations: 2067 2068 .. math:: 2069 2070 D(n)=\sum_{i=1}^{3}\sum_{j=1}^{3}n_{i}n_{j}D_{ij} 2071 2072 and 2073 2074 .. math:: 2075 2076 K(n)=\frac{MD^{2}}{D(n)^{2}}\sum_{i=1}^{3}\sum_{j=1}^{3} 2077 \sum_{k=1}^{3}\sum_{l=1}^{3}n_{i}n_{j}n_{k}n_{l}W_{ijkl} 2078 2079 where $D_{ij}$ and $W_{ijkl}$ are the elements of the second-order DT 2080 and the fourth-order KT tensors, respectively, and $MD$ is the mean 2081 diffusivity. 2082 """ 2083 return dki_prediction(self.model_params, gtab, S0) 2084 2085 2086def _ols_iter(inv_design, sig, min_diffusivity): 2087 """ Helper function used by ols_fit_dki - Applies OLS fit of the diffusion 2088 kurtosis model to single voxel signals. 2089 2090 Parameters 2091 ---------- 2092 inv_design : array (g, 22) 2093 Inverse of the design matrix holding the covariants used to solve for 2094 the regression coefficients. 2095 sig : array (g,) 2096 Diffusion-weighted signal for a single voxel data. 2097 min_diffusivity : float 2098 Because negative eigenvalues are not physical and small eigenvalues, 2099 much smaller than the diffusion weighting, cause quite a lot of noise 2100 in metrics such as fa, diffusivity values smaller than 2101 `min_diffusivity` are replaced with `min_diffusivity`. 2102 2103 Returns 2104 ------- 2105 dki_params : array (27,) 2106 All parameters estimated from the diffusion kurtosis model. 2107 Parameters are ordered as follows: 2108 1) Three diffusion tensor's eigenvalues 2109 2) Three lines of the eigenvector matrix each containing the first, 2110 second and third coordinates of the eigenvector 2111 3) Fifteen elements of the kurtosis tensor 2112 """ 2113 # DKI ordinary linear least square solution 2114 log_s = np.log(sig) 2115 result = np.dot(inv_design, log_s) 2116 2117 # Extracting the diffusion tensor parameters from solution 2118 DT_elements = result[:6] 2119 evals, evecs = decompose_tensor(from_lower_triangular(DT_elements), 2120 min_diffusivity=min_diffusivity) 2121 2122 # Extracting kurtosis tensor parameters from solution 2123 MD_square = (evals.mean(0))**2 2124 KT_elements = result[6:21] / MD_square 2125 2126 # Write output 2127 dki_params = np.concatenate((evals, evecs[0], evecs[1], evecs[2], 2128 KT_elements), axis=0) 2129 2130 return dki_params 2131 2132 2133def ols_fit_dki(design_matrix, data): 2134 r""" Computes the diffusion and kurtosis tensors using an ordinary linear 2135 least squares (OLS) approach [1]_. 2136 2137 Parameters 2138 ---------- 2139 design_matrix : array (g, 22) 2140 Design matrix holding the covariants used to solve for the regression 2141 coefficients. 2142 data : array (N, g) 2143 Data or response variables holding the data. Note that the last 2144 dimension should contain the data. It makes no copies of data. 2145 2146 Returns 2147 ------- 2148 dki_params : array (N, 27) 2149 All parameters estimated from the diffusion kurtosis model. 2150 Parameters are ordered as follows: 2151 1) Three diffusion tensor's eigenvalues 2152 2) Three lines of the eigenvector matrix each containing the first, 2153 second and third coordinates of the eigenvector 2154 3) Fifteen elements of the kurtosis tensor 2155 2156 See Also 2157 -------- 2158 wls_fit_dki, nls_fit_dki 2159 2160 References 2161 ---------- 2162 [1] Lu, H., Jensen, J. H., Ramani, A., and Helpern, J. A. (2006). 2163 Three-dimensional characterization of non-gaussian water diffusion in 2164 humans using diffusion kurtosis imaging. NMR in Biomedicine 19, 2165 236–247. doi:10.1002/nbm.1020 2166 """ 2167 tol = 1e-6 2168 2169 # preparing data and initializing parameters 2170 data = np.asarray(data) 2171 data_flat = data.reshape((-1, data.shape[-1])) 2172 dki_params = np.empty((len(data_flat), 27)) 2173 2174 # inverting design matrix and defining minimum diffusion 2175 min_diffusivity = tol / -design_matrix.min() 2176 inv_design = np.linalg.pinv(design_matrix) 2177 2178 # looping OLS solution on all data voxels 2179 for vox in range(len(data_flat)): 2180 dki_params[vox] = _ols_iter(inv_design, data_flat[vox], 2181 min_diffusivity) 2182 2183 # Reshape data according to the input data shape 2184 dki_params = dki_params.reshape((data.shape[:-1]) + (27,)) 2185 2186 return dki_params 2187 2188 2189def _wls_iter(design_matrix, inv_design, sig, min_diffusivity): 2190 """ Helper function used by wls_fit_dki - Applies WLS fit of the diffusion 2191 kurtosis model to single voxel signals. 2192 2193 Parameters 2194 ---------- 2195 design_matrix : array (g, 22) 2196 Design matrix holding the covariants used to solve for the regression 2197 coefficients 2198 inv_design : array (g, 22) 2199 Inverse of the design matrix. 2200 sig : array (g, ) 2201 Diffusion-weighted signal for a single voxel data. 2202 min_diffusivity : float 2203 Because negative eigenvalues are not physical and small eigenvalues, 2204 much smaller than the diffusion weighting, cause quite a lot of noise 2205 in metrics such as fa, diffusivity values smaller than 2206 `min_diffusivity` are replaced with `min_diffusivity`. 2207 2208 Returns 2209 ------- 2210 dki_params : array (27, ) 2211 All parameters estimated from the diffusion kurtosis model. 2212 Parameters are ordered as follows: 2213 1) Three diffusion tensor's eigenvalues 2214 2) Three lines of the eigenvector matrix each containing the first, 2215 second and third coordinates of the eigenvector 2216 3) Fifteen elements of the kurtosis tensor 2217 """ 2218 A = design_matrix 2219 2220 # DKI ordinary linear least square solution (initial guess) 2221 log_s = np.log(sig) 2222 ols_result = np.dot(inv_design, log_s) 2223 2224 # Define weights as diag(yn**2) 2225 W = np.diag(np.exp(2 * np.dot(A, ols_result))) 2226 2227 # DKI weighted linear least square solution 2228 inv_AT_W_A = np.linalg.pinv(np.dot(np.dot(A.T, W), A)) 2229 AT_W_LS = np.dot(np.dot(A.T, W), log_s) 2230 wls_result = np.dot(inv_AT_W_A, AT_W_LS) 2231 2232 # Extracting the diffusion tensor parameters from solution 2233 DT_elements = wls_result[:6] 2234 evals, evecs = decompose_tensor(from_lower_triangular(DT_elements), 2235 min_diffusivity=min_diffusivity) 2236 2237 # Extracting kurtosis tensor parameters from solution 2238 MD_square = (evals.mean(0))**2 2239 KT_elements = wls_result[6:21] / MD_square 2240 2241 # Write output 2242 dki_params = np.concatenate((evals, evecs[0], evecs[1], evecs[2], 2243 KT_elements), axis=0) 2244 2245 return dki_params 2246 2247 2248def wls_fit_dki(design_matrix, data): 2249 r""" Computes the diffusion and kurtosis tensors using a weighted linear 2250 least squares (WLS) approach [1]_. 2251 2252 Parameters 2253 ---------- 2254 design_matrix : array (g, 22) 2255 Design matrix holding the covariants used to solve for the regression 2256 coefficients. 2257 data : array (N, g) 2258 Data or response variables holding the data. Note that the last 2259 dimension should contain the data. It makes no copies of data. 2260 2261 Returns 2262 ------- 2263 dki_params : array (N, 27) 2264 All parameters estimated from the diffusion kurtosis model for all N 2265 voxels. 2266 Parameters are ordered as follows: 2267 1) Three diffusion tensor's eigenvalues 2268 2) Three lines of the eigenvector matrix each containing the first 2269 second and third coordinates of the eigenvector 2270 3) Fifteen elements of the kurtosis tensor 2271 2272 References 2273 ---------- 2274 [1] Veraart, J., Sijbers, J., Sunaert, S., Leemans, A., Jeurissen, B., 2275 2013. Weighted linear least squares estimation of diffusion MRI 2276 parameters: Strengths, limitations, and pitfalls. Magn Reson Med 81, 2277 335-346. 2278 """ 2279 2280 tol = 1e-6 2281 2282 # preparing data and initializing parameters 2283 data = np.asarray(data) 2284 data_flat = data.reshape((-1, data.shape[-1])) 2285 dki_params = np.empty((len(data_flat), 27)) 2286 2287 # inverting design matrix and defining minimum diffusion 2288 min_diffusivity = tol / -design_matrix.min() 2289 inv_design = np.linalg.pinv(design_matrix) 2290 2291 # looping WLS solution on all data voxels 2292 for vox in range(len(data_flat)): 2293 dki_params[vox] = _wls_iter(design_matrix, inv_design, data_flat[vox], 2294 min_diffusivity) 2295 2296 # Reshape data according to the input data shape 2297 dki_params = dki_params.reshape((data.shape[:-1]) + (27,)) 2298 2299 return dki_params 2300 2301 2302def Wrotate(kt, Basis): 2303 r""" Rotate a kurtosis tensor from the standard Cartesian coordinate system 2304 to another coordinate system basis 2305 2306 Parameters 2307 ---------- 2308 kt : (15,) 2309 Vector with the 15 independent elements of the kurtosis tensor 2310 Basis : array (3, 3) 2311 Vectors of the basis column-wise oriented 2312 inds : array(m, 4) (optional) 2313 Array of vectors containing the four indexes of m specific elements of 2314 the rotated kurtosis tensor. If not specified all 15 elements of the 2315 rotated kurtosis tensor are computed. 2316 2317 Returns 2318 -------- 2319 Wrot : array (m,) or (15,) 2320 Vector with the m independent elements of the rotated kurtosis tensor. 2321 If 'indices' is not specified all 15 elements of the rotated kurtosis 2322 tensor are computed. 2323 2324 Notes 2325 ------ 2326 KT elements are assumed to be ordered as follows: 2327 2328 .. math:: 2329 2330 \begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz} 2331 & ... \\ 2332 & W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy} 2333 & ... \\ 2334 & W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz} 2335 & & )\end{matrix} 2336 2337 References 2338 ---------- 2339 [1] Hui ES, Cheung MM, Qi L, Wu EX, 2008. Towards better MR 2340 characterization of neural tissues using directional diffusion kurtosis 2341 analysis. Neuroimage 42(1): 122-34 2342 """ 2343 inds = np.array([[0, 0, 0, 0], [1, 1, 1, 1], [2, 2, 2, 2], 2344 [0, 0, 0, 1], [0, 0, 0, 2], [0, 1, 1, 1], 2345 [1, 1, 1, 2], [0, 2, 2, 2], [1, 2, 2, 2], 2346 [0, 0, 1, 1], [0, 0, 2, 2], [1, 1, 2, 2], 2347 [0, 0, 1, 2], [0, 1, 1, 2], [0, 1, 2, 2]]) 2348 2349 Wrot = np.zeros(kt.shape) 2350 2351 for e in range(len(inds)): 2352 Wrot[..., e] = Wrotate_element(kt, inds[e][0], inds[e][1], inds[e][2], 2353 inds[e][3], Basis) 2354 2355 return Wrot 2356 2357 2358# Defining keys to select a kurtosis tensor element with indexes (i, j, k, l) 2359# on a kt vector that contains only the 15 independent elements of the kurtosis 2360# tensor: Considering y defined by (i+1) * (j+1) * (k+1) * (l+1). Two elements 2361# of the full 4D kurtosis tensor are equal if y obtain from the indexes of 2362# these two element are equal. Therefore, the possible values of y (1, 16, 81, 2363# 2, 3, 8, 24 27, 54, 4, 9, 36, 6, 12, 18) are used to point each element of 2364# the kurtosis tensor on the format of a vector containing the 15 independent 2365# elements. 2366ind_ele = {1: 0, 16: 1, 81: 2, 2: 3, 3: 4, 8: 5, 24: 6, 27: 7, 54: 8, 4: 9, 2367 9: 10, 36: 11, 6: 12, 12: 13, 18: 14} 2368 2369 2370def Wrotate_element(kt, indi, indj, indk, indl, B): 2371 r""" Computes the the specified index element of a kurtosis tensor rotated 2372 to the coordinate system basis B. 2373 2374 Parameters 2375 ---------- 2376 kt : ndarray (x, y, z, 15) or (n, 15) 2377 Array containing the 15 independent elements of the kurtosis tensor 2378 indi : int 2379 Rotated kurtosis tensor element index i (0 for x, 1 for y, 2 for z) 2380 indj : int 2381 Rotated kurtosis tensor element index j (0 for x, 1 for y, 2 for z) 2382 indk : int 2383 Rotated kurtosis tensor element index k (0 for x, 1 for y, 2 for z) 2384 indl: int 2385 Rotated kurtosis tensor element index l (0 for x, 1 for y, 2 for z) 2386 B: array (x, y, z, 3, 3) or (n, 15) 2387 Vectors of the basis column-wise oriented 2388 2389 Returns 2390 ------- 2391 Wre : float 2392 rotated kurtosis tensor element of index ind_i, ind_j, ind_k, ind_l 2393 2394 Notes 2395 ------ 2396 It is assumed that initial kurtosis tensor elementes are defined on the 2397 Cartesian coordinate system. 2398 2399 References 2400 ---------- 2401 [1] Hui ES, Cheung MM, Qi L, Wu EX, 2008. Towards better MR 2402 characterization of neural tissues using directional diffusion kurtosis 2403 analysis. Neuroimage 42(1): 122-34 2404 """ 2405 Wre = 0 2406 2407 xyz = [0, 1, 2] 2408 for il in xyz: 2409 for jl in xyz: 2410 for kl in xyz: 2411 for ll in xyz: 2412 key = (il + 1) * (jl + 1) * (kl + 1) * (ll + 1) 2413 multiplyB = \ 2414 B[..., il, indi] * B[..., jl, indj] * \ 2415 B[..., kl, indk] * B[..., ll, indl] 2416 Wre = Wre + multiplyB * kt[..., ind_ele[key]] 2417 return Wre 2418 2419 2420def Wcons(k_elements): 2421 r""" Construct the full 4D kurtosis tensors from its 15 independent 2422 elements 2423 2424 Parameters 2425 ---------- 2426 k_elements : (15,) 2427 elements of the kurtosis tensor in the following order: 2428 2429 .. math:: 2430 2431 \begin{matrix} ( & W_{xxxx} & W_{yyyy} & W_{zzzz} & W_{xxxy} & W_{xxxz} 2432 & ... \\ 2433 & W_{xyyy} & W_{yyyz} & W_{xzzz} & W_{yzzz} & W_{xxyy} 2434 & ... \\ 2435 & W_{xxzz} & W_{yyzz} & W_{xxyz} & W_{xyyz} & W_{xyzz} 2436 & & )\end{matrix} 2437 2438 Returns 2439 ------- 2440 W : array(3, 3, 3, 3) 2441 Full 4D kurtosis tensor 2442 """ 2443 W = np.zeros((3, 3, 3, 3)) 2444 2445 xyz = [0, 1, 2] 2446 for ind_i in xyz: 2447 for ind_j in xyz: 2448 for ind_k in xyz: 2449 for ind_l in xyz: 2450 key = (ind_i + 1) * (ind_j + 1) * (ind_k + 1) * (ind_l + 1) 2451 W[ind_i][ind_j][ind_k][ind_l] = k_elements[ind_ele[key]] 2452 2453 return W 2454 2455 2456def split_dki_param(dki_params): 2457 r""" Extract the diffusion tensor eigenvalues, the diffusion tensor 2458 eigenvector matrix, and the 15 independent elements of the kurtosis tensor 2459 from the model parameters estimated from the DKI model 2460 2461 Parameters 2462 ---------- 2463 dki_params : ndarray (x, y, z, 27) or (n, 27) 2464 All parameters estimated from the diffusion kurtosis model. 2465 Parameters are ordered as follows: 2466 1) Three diffusion tensor's eigenvalues 2467 2) Three lines of the eigenvector matrix each containing the first, 2468 second and third coordinates of the eigenvector 2469 3) Fifteen elements of the kurtosis tensor 2470 2471 Returns 2472 -------- 2473 eigvals : array (x, y, z, 3) or (n, 3) 2474 Eigenvalues from eigen decomposition of the tensor. 2475 eigvecs : array (x, y, z, 3, 3) or (n, 3, 3) 2476 Associated eigenvectors from eigen decomposition of the tensor. 2477 Eigenvectors are columnar (e.g. eigvecs[:,j] is associated with 2478 eigvals[j]) 2479 kt : array (x, y, z, 15) or (n, 15) 2480 Fifteen elements of the kurtosis tensor 2481 """ 2482 evals = dki_params[..., :3] 2483 evecs = dki_params[..., 3:12].reshape(dki_params.shape[:-1] + (3, 3)) 2484 kt = dki_params[..., 12:] 2485 2486 return evals, evecs, kt 2487 2488 2489common_fit_methods = {'WLS': wls_fit_dki, 2490 'OLS': ols_fit_dki, 2491 'NLS': nlls_fit_tensor, 2492 'UWLLS': wls_fit_dki, 2493 'ULLS': ols_fit_dki, 2494 'WLLS': wls_fit_dki, 2495 'OLLS': ols_fit_dki, 2496 'NLLS': nlls_fit_tensor, 2497 'RT': restore_fit_tensor, 2498 'restore': restore_fit_tensor, 2499 'RESTORE': restore_fit_tensor 2500 } 2501