1 2import logging 3import numpy as np 4import os.path 5from ast import literal_eval 6from warnings import warn 7 8import nibabel as nib 9 10from dipy.core.gradients import gradient_table 11from dipy.data import default_sphere 12from dipy.io.gradients import read_bvals_bvecs 13from dipy.io.peaks import save_peaks, peaks_to_niftis 14from dipy.io.image import load_nifti, save_nifti, load_nifti_data 15from dipy.io.utils import nifti1_symmat 16from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel, 17 auto_response_ssst) 18from dipy.reconst.dti import (TensorModel, color_fa, fractional_anisotropy, 19 geodesic_anisotropy, mean_diffusivity, 20 axial_diffusivity, radial_diffusivity, 21 lower_triangular, mode as get_mode) 22from dipy.direction.peaks import peaks_from_model 23from dipy.reconst.shm import CsaOdfModel 24from dipy.workflows.workflow import Workflow 25from dipy.reconst.dki import DiffusionKurtosisModel, split_dki_param 26from dipy.reconst.ivim import IvimModel 27 28from dipy.reconst import mapmri 29 30from dipy.utils.deprecator import deprecated_params 31 32 33class ReconstMAPMRIFlow(Workflow): 34 @classmethod 35 def get_short_name(cls): 36 return 'mapmri' 37 38 def run(self, data_files, bvals_files, bvecs_files, small_delta, big_delta, 39 b0_threshold=50.0, laplacian=True, positivity=True, 40 bval_threshold=2000, save_metrics=[], 41 laplacian_weighting=0.05, radial_order=6, out_dir='', 42 out_rtop='rtop.nii.gz', out_lapnorm='lapnorm.nii.gz', 43 out_msd='msd.nii.gz', out_qiv='qiv.nii.gz', 44 out_rtap='rtap.nii.gz', 45 out_rtpp='rtpp.nii.gz', out_ng='ng.nii.gz', 46 out_perng='perng.nii.gz', 47 out_parng='parng.nii.gz'): 48 """Workflow for fitting the MAPMRI model (with optional Laplacian 49 regularization). Generates rtop, lapnorm, msd, qiv, rtap, rtpp, 50 non-gaussian (ng), parallel ng, perpendicular ng saved in a nifti 51 format in input files provided by `data_files` and saves the nifti 52 files to an output directory specified by `out_dir`. 53 54 In order for the MAPMRI workflow to work in the way 55 intended either the Laplacian or positivity or both must 56 be set to True. 57 58 Parameters 59 ---------- 60 data_files : string 61 Path to the input volume. 62 bvals_files : string 63 Path to the bval files. 64 bvecs_files : string 65 Path to the bvec files. 66 small_delta : float 67 Small delta value used in generation of gradient table of provided 68 bval and bvec. 69 big_delta : float 70 Big delta value used in generation of gradient table of provided 71 bval and bvec. 72 b0_threshold : float, optional 73 Threshold used to find b0 volumes. 74 laplacian : bool, optional 75 Regularize using the Laplacian of the MAP-MRI basis. 76 positivity : bool, optional 77 Constrain the propagator to be positive. 78 bval_threshold : float, optional 79 Sets the b-value threshold to be used in the scale factor 80 estimation. In order for the estimated non-Gaussianity to have 81 meaning this value should set to a lower value (b<2000 s/mm^2) 82 such that the scale factors are estimated on signal points that 83 reasonably represent the spins at Gaussian diffusion. 84 save_metrics : variable string, optional 85 List of metrics to save. 86 Possible values: rtop, laplacian_signal, msd, qiv, rtap, rtpp, 87 ng, perng, parng 88 laplacian_weighting : float, optional 89 Weighting value used in fitting the MAPMRI model in the Laplacian 90 and both model types. 91 radial_order : unsigned int, optional 92 Even value used to set the order of the basis. 93 out_dir : string, optional 94 Output directory. (default: current directory) 95 out_rtop : string, optional 96 Name of the rtop to be saved. 97 out_lapnorm : string, optional 98 Name of the norm of Laplacian signal to be saved. 99 out_msd : string, optional 100 Name of the msd to be saved. 101 out_qiv : string, optional 102 Name of the qiv to be saved. 103 out_rtap : string, optional 104 Name of the rtap to be saved. 105 out_rtpp : string, optional 106 Name of the rtpp to be saved. 107 out_ng : string, optional 108 Name of the Non-Gaussianity to be saved. 109 out_perng : string, optional 110 Name of the Non-Gaussianity perpendicular to be saved. 111 out_parng : string, optional 112 Name of the Non-Gaussianity parallel to be saved. 113 """ 114 io_it = self.get_io_iterator() 115 for (dwi, bval, bvec, out_rtop, out_lapnorm, out_msd, out_qiv, 116 out_rtap, out_rtpp, out_ng, out_perng, out_parng) in io_it: 117 118 logging.info('Computing MAPMRI metrics for {0}'.format(dwi)) 119 data, affine = load_nifti(dwi) 120 121 bvals, bvecs = read_bvals_bvecs(bval, bvec) 122 if b0_threshold < bvals.min(): 123 warn("b0_threshold (value: {0}) is too low, increase your " 124 "b0_threshold. It should be higher than the first b0 " 125 "value({1}).".format(b0_threshold, bvals.min())) 126 gtab = gradient_table(bvals=bvals, bvecs=bvecs, 127 small_delta=small_delta, 128 big_delta=big_delta, 129 b0_threshold=b0_threshold) 130 131 if not save_metrics: 132 save_metrics = ['rtop', 'laplacian_signal', 'msd', 133 'qiv', 'rtap', 'rtpp', 134 'ng', 'perng', 'parng'] 135 136 if laplacian and positivity: 137 map_model_aniso = mapmri.MapmriModel( 138 gtab, 139 radial_order=radial_order, 140 laplacian_regularization=True, 141 laplacian_weighting=laplacian_weighting, 142 positivity_constraint=True, 143 bval_threshold=bval_threshold) 144 145 mapfit_aniso = map_model_aniso.fit(data) 146 147 elif positivity: 148 map_model_aniso = mapmri.MapmriModel( 149 gtab, 150 radial_order=radial_order, 151 laplacian_regularization=False, 152 positivity_constraint=True, 153 bval_threshold=bval_threshold) 154 mapfit_aniso = map_model_aniso.fit(data) 155 156 elif laplacian: 157 map_model_aniso = mapmri.MapmriModel( 158 gtab, 159 radial_order=radial_order, 160 laplacian_regularization=True, 161 laplacian_weighting=laplacian_weighting, 162 bval_threshold=bval_threshold) 163 mapfit_aniso = map_model_aniso.fit(data) 164 165 else: 166 map_model_aniso = mapmri.MapmriModel( 167 gtab, 168 radial_order=radial_order, 169 laplacian_regularization=False, 170 positivity_constraint=False, 171 bval_threshold=bval_threshold) 172 mapfit_aniso = map_model_aniso.fit(data) 173 174 # for name, fname, func in [('rtop', out_rtop, mapfit_aniso.rtop), 175 # ]: 176 # if name in save_metrics: 177 # r = func() 178 # save_nifti(fname, r.astype(np.float32), affine) 179 180 if 'rtop' in save_metrics: 181 r = mapfit_aniso.rtop() 182 save_nifti(out_rtop, r.astype(np.float32), affine) 183 184 if 'laplacian_signal' in save_metrics: 185 ll = mapfit_aniso.norm_of_laplacian_signal() 186 save_nifti(out_lapnorm, ll.astype(np.float32), affine) 187 188 if 'msd' in save_metrics: 189 m = mapfit_aniso.msd() 190 save_nifti(out_msd, m.astype(np.float32), affine) 191 192 if 'qiv' in save_metrics: 193 q = mapfit_aniso.qiv() 194 save_nifti(out_qiv, q.astype(np.float32), affine) 195 196 if 'rtap' in save_metrics: 197 r = mapfit_aniso.rtap() 198 save_nifti(out_rtap, r.astype(np.float32), affine) 199 200 if 'rtpp' in save_metrics: 201 r = mapfit_aniso.rtpp() 202 save_nifti(out_rtpp, r.astype(np.float32), affine) 203 204 if 'ng' in save_metrics: 205 n = mapfit_aniso.ng() 206 save_nifti(out_ng, n.astype(np.float32), affine) 207 208 if 'perng' in save_metrics: 209 n = mapfit_aniso.ng_perpendicular() 210 save_nifti(out_perng, n.astype(np.float32), affine) 211 212 if 'parng' in save_metrics: 213 n = mapfit_aniso.ng_parallel() 214 save_nifti(out_parng, n.astype(np.float32), affine) 215 216 logging.info('MAPMRI saved in {0}'. 217 format(os.path.dirname(out_dir))) 218 219 220class ReconstDtiFlow(Workflow): 221 @classmethod 222 def get_short_name(cls): 223 return 'dti' 224 225 def run(self, input_files, bvalues_files, bvectors_files, mask_files, 226 b0_threshold=50, bvecs_tol=0.01, save_metrics=[], 227 out_dir='', out_tensor='tensors.nii.gz', out_fa='fa.nii.gz', 228 out_ga='ga.nii.gz', out_rgb='rgb.nii.gz', out_md='md.nii.gz', 229 out_ad='ad.nii.gz', out_rd='rd.nii.gz', out_mode='mode.nii.gz', 230 out_evec='evecs.nii.gz', out_eval='evals.nii.gz', nifti_tensor=True): 231 """ Workflow for tensor reconstruction and for computing DTI metrics. 232 using Weighted Least-Squares. 233 Performs a tensor reconstruction on the files by 'globing' 234 ``input_files`` and saves the DTI metrics in a directory specified by 235 ``out_dir``. 236 237 Parameters 238 ---------- 239 input_files : string 240 Path to the input volumes. This path may contain wildcards to 241 process multiple inputs at once. 242 bvalues_files : string 243 Path to the bvalues files. This path may contain wildcards to use 244 multiple bvalues files at once. 245 bvectors_files : string 246 Path to the bvectors files. This path may contain wildcards to use 247 multiple bvectors files at once. 248 mask_files : string 249 Path to the input masks. This path may contain wildcards to use 250 multiple masks at once. 251 b0_threshold : float, optional 252 Threshold used to find b0 volumes. 253 bvecs_tol : float, optional 254 Threshold used to check that norm(bvec) = 1 +/- bvecs_tol 255 b-vectors are unit vectors. 256 save_metrics : variable string, optional 257 List of metrics to save. 258 Possible values: fa, ga, rgb, md, ad, rd, mode, tensor, evec, eval 259 out_dir : string, optional 260 Output directory. (default current directory) 261 out_tensor : string, optional 262 Name of the tensors volume to be saved. 263 Per default, this will be saved following the nifti standard: 264 with the tensor elements as Dxx, Dxy, Dyy, Dxz, Dyz, Dzz on the 265 last (5th) dimension of the volume (shape: (i, j, k, 1, 6)). If 266 `nifti_tensor` is False, this will be saved in an alternate format 267 that is used by other software (e.g., FSL): a 268 4-dimensional volume (shape (i, j, k, 6)) with Dxx, Dxy, Dxz, Dyy, 269 Dyz, Dzz on the last dimension. 270 out_fa : string, optional 271 Name of the fractional anisotropy volume to be saved. 272 out_ga : string, optional 273 Name of the geodesic anisotropy volume to be saved. 274 out_rgb : string, optional 275 Name of the color fa volume to be saved. 276 out_md : string, optional 277 Name of the mean diffusivity volume to be saved. 278 out_ad : string, optional 279 Name of the axial diffusivity volume to be saved. 280 out_rd : string, optional 281 Name of the radial diffusivity volume to be saved. 282 out_mode : string, optional 283 Name of the mode volume to be saved. 284 out_evec : string, optional 285 Name of the eigenvectors volume to be saved. 286 out_eval : string, optional 287 Name of the eigenvalues to be saved. 288 nifti_tensor : bool, optional 289 Whether the tensor is saved in the standard Nifti format or in an 290 alternate format 291 that is used by other software (e.g., FSL): a 292 4-dimensional volume (shape (i, j, k, 6)) with 293 Dxx, Dxy, Dxz, Dyy, Dyz, Dzz on the last dimension. 294 295 References 296 ---------- 297 .. [1] Basser, P.J., Mattiello, J., LeBihan, D., 1994. Estimation of 298 the effective self-diffusion tensor from the NMR spin echo. J Magn 299 Reson B 103, 247-254. 300 301 .. [2] Basser, P., Pierpaoli, C., 1996. Microstructural and 302 physiological features of tissues elucidated by quantitative 303 diffusion-tensor MRI. Journal of Magnetic Resonance 111, 209-219. 304 305 .. [3] Lin-Ching C., Jones D.K., Pierpaoli, C. 2005. RESTORE: Robust 306 estimation of tensors by outlier rejection. MRM 53: 1088-1095 307 308 .. [4] hung, SW., Lu, Y., Henry, R.G., 2006. Comparison of bootstrap 309 approaches for estimation of uncertainties of DTI parameters. 310 NeuroImage 33, 531-541. 311 312 """ 313 io_it = self.get_io_iterator() 314 315 for dwi, bval, bvec, mask, otensor, ofa, oga, orgb, omd, oad, orad, \ 316 omode, oevecs, oevals in io_it: 317 318 logging.info('Computing DTI metrics for {0}'.format(dwi)) 319 data, affine = load_nifti(dwi) 320 321 if mask is not None: 322 mask = load_nifti_data(mask).astype(bool) 323 324 tenfit, _ = self.get_fitted_tensor(data, mask, bval, bvec, 325 b0_threshold, bvecs_tol) 326 327 if not save_metrics: 328 save_metrics = ['fa', 'md', 'rd', 'ad', 'ga', 'rgb', 'mode', 329 'evec', 'eval', 'tensor'] 330 331 FA = fractional_anisotropy(tenfit.evals) 332 FA[np.isnan(FA)] = 0 333 FA = np.clip(FA, 0, 1) 334 335 if 'tensor' in save_metrics: 336 tensor_vals = lower_triangular(tenfit.quadratic_form) 337 338 if nifti_tensor: 339 ten_img = nifti1_symmat(tensor_vals, affine=affine) 340 else: 341 alt_order = [0, 1, 3, 2, 4, 5] 342 ten_img = nib.Nifti1Image( 343 tensor_vals[..., alt_order].astype(np.float32), 344 affine) 345 346 nib.save(ten_img, otensor) 347 348 if 'fa' in save_metrics: 349 save_nifti(ofa, FA.astype(np.float32), affine) 350 351 if 'ga' in save_metrics: 352 GA = geodesic_anisotropy(tenfit.evals) 353 save_nifti(oga, GA.astype(np.float32), affine) 354 355 if 'rgb' in save_metrics: 356 RGB = color_fa(FA, tenfit.evecs) 357 save_nifti(orgb, np.array(255 * RGB, 'uint8'), affine) 358 359 if 'md' in save_metrics: 360 MD = mean_diffusivity(tenfit.evals) 361 save_nifti(omd, MD.astype(np.float32), affine) 362 363 if 'ad' in save_metrics: 364 AD = axial_diffusivity(tenfit.evals) 365 save_nifti(oad, AD.astype(np.float32), affine) 366 367 if 'rd' in save_metrics: 368 RD = radial_diffusivity(tenfit.evals) 369 save_nifti(orad, RD.astype(np.float32), affine) 370 371 if 'mode' in save_metrics: 372 MODE = get_mode(tenfit.quadratic_form) 373 save_nifti(omode, MODE.astype(np.float32), affine) 374 375 if 'evec' in save_metrics: 376 save_nifti(oevecs, tenfit.evecs.astype(np.float32), affine) 377 378 if 'eval' in save_metrics: 379 save_nifti(oevals, tenfit.evals.astype(np.float32), affine) 380 381 dname_ = os.path.dirname(oevals) 382 if dname_ == '': 383 logging.info('DTI metrics saved in current directory') 384 else: 385 logging.info( 386 'DTI metrics saved in {0}'.format(dname_)) 387 388 def get_tensor_model(self, gtab): 389 return TensorModel(gtab, fit_method="WLS") 390 391 def get_fitted_tensor(self, data, mask, bval, bvec, 392 b0_threshold=50, bvecs_tol=0.01): 393 394 logging.info('Tensor estimation...') 395 bvals, bvecs = read_bvals_bvecs(bval, bvec) 396 gtab = gradient_table(bvals, bvecs, b0_threshold=b0_threshold, 397 atol=bvecs_tol) 398 399 tenmodel = self.get_tensor_model(gtab) 400 tenfit = tenmodel.fit(data, mask) 401 402 return tenfit, gtab 403 404 405class ReconstCSDFlow(Workflow): 406 @classmethod 407 def get_short_name(cls): 408 return 'csd' 409 410 @deprecated_params('nbr_processes', 'num_processes', since='1.4', 411 until='1.5') 412 def run(self, input_files, bvalues_files, bvectors_files, mask_files, 413 b0_threshold=50.0, bvecs_tol=0.01, roi_center=None, roi_radii=10, 414 fa_thr=0.7, frf=None, extract_pam_values=False, sh_order=8, 415 odf_to_sh_order=8, parallel=False, num_processes=None, 416 out_dir='', 417 out_pam='peaks.pam5', out_shm='shm.nii.gz', 418 out_peaks_dir='peaks_dirs.nii.gz', 419 out_peaks_values='peaks_values.nii.gz', 420 out_peaks_indices='peaks_indices.nii.gz', out_gfa='gfa.nii.gz'): 421 """ Constrained spherical deconvolution 422 423 Parameters 424 ---------- 425 input_files : string 426 Path to the input volumes. This path may contain wildcards to 427 process multiple inputs at once. 428 bvalues_files : string 429 Path to the bvalues files. This path may contain wildcards to use 430 multiple bvalues files at once. 431 bvectors_files : string 432 Path to the bvectors files. This path may contain wildcards to use 433 multiple bvectors files at once. 434 mask_files : string 435 Path to the input masks. This path may contain wildcards to use 436 multiple masks at once. (default: No mask used) 437 b0_threshold : float, optional 438 Threshold used to find b0 volumes. 439 bvecs_tol : float, optional 440 Bvecs should be unit vectors. 441 roi_center : variable int, optional 442 Center of ROI in data. If center is None, it is assumed that it is 443 the center of the volume with shape `data.shape[:3]`. 444 roi_radii : int or array-like, optional 445 radii of cuboid ROI in voxels. 446 fa_thr : float, optional 447 FA threshold for calculating the response function. 448 frf : variable float, optional 449 Fiber response function can be for example inputed as 15 4 4 450 (from the command line) or [15, 4, 4] from a Python script to be 451 converted to float and multiplied by 10**-4 . If None 452 the fiber response function will be computed automatically. 453 extract_pam_values : bool, optional 454 Save or not to save pam volumes as single nifti files. 455 sh_order : int, optional 456 Spherical harmonics order used in the CSA fit. 457 odf_to_sh_order : int, optional 458 Spherical harmonics order used for peak_from_model to compress 459 the ODF to spherical harmonics coefficients. 460 parallel : bool, optional 461 Whether to use parallelization in peak-finding during the 462 calibration procedure. 463 num_processes : int, optional 464 If `parallel` is True, the number of subprocesses to use 465 (default multiprocessing.cpu_count()). If < 0 the maximal number 466 of cores minus |num_processes + 1| is used (enter -1 to use as 467 many cores as possible). 0 raises an error. 468 out_dir : string, optional 469 Output directory. (default current directory) 470 out_pam : string, optional 471 Name of the peaks volume to be saved. 472 out_shm : string, optional 473 Name of the spherical harmonics volume to be saved. 474 out_peaks_dir : string, optional 475 Name of the peaks directions volume to be saved. 476 out_peaks_values : string, optional 477 Name of the peaks values volume to be saved. 478 out_peaks_indices : string, optional 479 Name of the peaks indices volume to be saved. 480 out_gfa : string, optional 481 Name of the generalized FA volume to be saved. 482 483 484 References 485 ---------- 486 .. [1] Tournier, J.D., et al. NeuroImage 2007. Robust determination of 487 the fibre orientation distribution in diffusion MRI: Non-negativity 488 constrained super-resolved spherical deconvolution. 489 """ 490 io_it = self.get_io_iterator() 491 492 for (dwi, bval, bvec, maskfile, opam, oshm, opeaks_dir, opeaks_values, 493 opeaks_indices, ogfa) in io_it: 494 495 logging.info('Loading {0}'.format(dwi)) 496 data, affine = load_nifti(dwi) 497 498 bvals, bvecs = read_bvals_bvecs(bval, bvec) 499 print(b0_threshold, bvals.min()) 500 if b0_threshold < bvals.min(): 501 warn("b0_threshold (value: {0}) is too low, increase your " 502 "b0_threshold. It should be higher than the first b0 value " 503 "({1}).".format(b0_threshold, bvals.min())) 504 gtab = gradient_table(bvals, bvecs, b0_threshold=b0_threshold, 505 atol=bvecs_tol) 506 mask_vol = load_nifti_data(maskfile).astype(bool) 507 508 n_params = ((sh_order + 1) * (sh_order + 2)) / 2 509 if data.shape[-1] < n_params: 510 raise ValueError( 511 'You need at least {0} unique DWI volumes to ' 512 'compute fiber odfs. You currently have: {1}' 513 ' DWI volumes.'.format(n_params, data.shape[-1])) 514 515 if frf is None: 516 logging.info('Computing response function') 517 if roi_center is not None: 518 logging.info('Response ROI center:\n{0}' 519 .format(roi_center)) 520 logging.info('Response ROI radii:\n{0}' 521 .format(roi_radii)) 522 response, ratio = auto_response_ssst( 523 gtab, data, 524 roi_center=roi_center, 525 roi_radii=roi_radii, 526 fa_thr=fa_thr) 527 response = list(response) 528 529 else: 530 logging.info('Using response function') 531 if isinstance(frf, str): 532 l01 = np.array(literal_eval(frf), dtype=np.float64) 533 else: 534 l01 = np.array(frf, dtype=np.float64) 535 536 l01 *= 10 ** -4 537 response = np.array([l01[0], l01[1], l01[1]]) 538 ratio = l01[1] / l01[0] 539 response = (response, ratio) 540 541 logging.info("Eigenvalues for the frf of the input" 542 " data are :{0}".format(response[0])) 543 logging.info('Ratio for smallest to largest eigen value is {0}' 544 .format(ratio)) 545 546 peaks_sphere = default_sphere 547 548 logging.info('CSD computation started.') 549 csd_model = ConstrainedSphericalDeconvModel(gtab, response, 550 sh_order=sh_order) 551 552 peaks_csd = peaks_from_model(model=csd_model, 553 data=data, 554 sphere=peaks_sphere, 555 relative_peak_threshold=.5, 556 min_separation_angle=25, 557 mask=mask_vol, 558 return_sh=True, 559 sh_order=sh_order, 560 normalize_peaks=True, 561 parallel=parallel, 562 num_processes=num_processes) 563 peaks_csd.affine = affine 564 565 save_peaks(opam, peaks_csd) 566 567 logging.info('CSD computation completed.') 568 569 if extract_pam_values: 570 peaks_to_niftis(peaks_csd, oshm, opeaks_dir, opeaks_values, 571 opeaks_indices, ogfa, reshape_dirs=True) 572 573 dname_ = os.path.dirname(opam) 574 if dname_ == '': 575 logging.info('Pam5 file saved in current directory') 576 else: 577 logging.info( 578 'Pam5 file saved in {0}'.format(dname_)) 579 580 return io_it 581 582 583class ReconstCSAFlow(Workflow): 584 @classmethod 585 def get_short_name(cls): 586 return 'csa' 587 588 @deprecated_params('nbr_processes', 'num_processes', since='1.4', 589 until='1.5') 590 def run(self, input_files, bvalues_files, bvectors_files, mask_files, 591 sh_order=6, odf_to_sh_order=8, b0_threshold=50.0, bvecs_tol=0.01, 592 extract_pam_values=False, parallel=False, num_processes=None, 593 out_dir='', 594 out_pam='peaks.pam5', out_shm='shm.nii.gz', 595 out_peaks_dir='peaks_dirs.nii.gz', 596 out_peaks_values='peaks_values.nii.gz', 597 out_peaks_indices='peaks_indices.nii.gz', 598 out_gfa='gfa.nii.gz'): 599 """ Constant Solid Angle. 600 601 Parameters 602 ---------- 603 input_files : string 604 Path to the input volumes. This path may contain wildcards to 605 process multiple inputs at once. 606 bvalues_files : string 607 Path to the bvalues files. This path may contain wildcards to use 608 multiple bvalues files at once. 609 bvectors_files : string 610 Path to the bvectors files. This path may contain wildcards to use 611 multiple bvectors files at once. 612 mask_files : string 613 Path to the input masks. This path may contain wildcards to use 614 multiple masks at once. (default: No mask used) 615 sh_order : int, optional 616 Spherical harmonics order used in the CSA fit. 617 odf_to_sh_order : int, optional 618 Spherical harmonics order used for peak_from_model to compress 619 the ODF to spherical harmonics coefficients. 620 b0_threshold : float, optional 621 Threshold used to find b0 volumes. 622 bvecs_tol : float, optional 623 Threshold used so that norm(bvec)=1. 624 extract_pam_values : bool, optional 625 Whether or not to save pam volumes as single nifti files. 626 parallel : bool, optional 627 Whether to use parallelization in peak-finding during the 628 calibration procedure. 629 num_processes : int, optional 630 If `parallel` is True, the number of subprocesses to use 631 (default multiprocessing.cpu_count()). If < 0 the maximal number 632 of cores minus |num_processes + 1| is used (enter -1 to use as 633 many cores as possible). 0 raises an error. 634 out_dir : string, optional 635 Output directory. (default current directory) 636 out_pam : string, optional 637 Name of the peaks volume to be saved. 638 out_shm : string, optional 639 Name of the spherical harmonics volume to be saved. 640 out_peaks_dir : string, optional 641 Name of the peaks directions volume to be saved. 642 out_peaks_values : string, optional 643 Name of the peaks values volume to be saved. 644 out_peaks_indices : string, optional 645 Name of the peaks indices volume to be saved. 646 out_gfa : string, optional 647 Name of the generalized FA volume to be saved. 648 649 References 650 ---------- 651 .. [1] Aganj, I., et al. 2009. ODF Reconstruction in Q-Ball Imaging 652 with Solid Angle Consideration. 653 """ 654 io_it = self.get_io_iterator() 655 656 for (dwi, bval, bvec, maskfile, opam, oshm, opeaks_dir, 657 opeaks_values, opeaks_indices, ogfa) in io_it: 658 659 logging.info('Loading {0}'.format(dwi)) 660 data, affine = load_nifti(dwi) 661 662 bvals, bvecs = read_bvals_bvecs(bval, bvec) 663 if b0_threshold < bvals.min(): 664 warn("b0_threshold (value: {0}) is too low, increase your " 665 "b0_threshold. It should be higher than the first b0 value " 666 "({1}).".format(b0_threshold, bvals.min())) 667 gtab = gradient_table(bvals, bvecs, 668 b0_threshold=b0_threshold, atol=bvecs_tol) 669 mask_vol = load_nifti_data(maskfile).astype(bool) 670 671 peaks_sphere = default_sphere 672 673 logging.info('Starting CSA computations {0}'.format(dwi)) 674 675 csa_model = CsaOdfModel(gtab, sh_order) 676 677 peaks_csa = peaks_from_model(model=csa_model, 678 data=data, 679 sphere=peaks_sphere, 680 relative_peak_threshold=.5, 681 min_separation_angle=25, 682 mask=mask_vol, 683 return_sh=True, 684 sh_order=odf_to_sh_order, 685 normalize_peaks=True, 686 parallel=parallel, 687 num_processes=num_processes) 688 peaks_csa.affine = affine 689 690 save_peaks(opam, peaks_csa) 691 692 logging.info('Finished CSA {0}'.format(dwi)) 693 694 if extract_pam_values: 695 peaks_to_niftis(peaks_csa, oshm, opeaks_dir, 696 opeaks_values, 697 opeaks_indices, ogfa, reshape_dirs=True) 698 699 dname_ = os.path.dirname(opam) 700 if dname_ == '': 701 logging.info('Pam5 file saved in current directory') 702 else: 703 logging.info( 704 'Pam5 file saved in {0}'.format(dname_)) 705 706 return io_it 707 708 709class ReconstDkiFlow(Workflow): 710 @classmethod 711 def get_short_name(cls): 712 return 'dki' 713 714 def run(self, input_files, bvalues_files, bvectors_files, mask_files, 715 b0_threshold=50.0, save_metrics=[], 716 out_dir='', out_dt_tensor='dti_tensors.nii.gz', out_fa='fa.nii.gz', 717 out_ga='ga.nii.gz', out_rgb='rgb.nii.gz', out_md='md.nii.gz', 718 out_ad='ad.nii.gz', out_rd='rd.nii.gz', out_mode='mode.nii.gz', 719 out_evec='evecs.nii.gz', out_eval='evals.nii.gz', 720 out_dk_tensor="dki_tensors.nii.gz", 721 out_mk="mk.nii.gz", out_ak="ak.nii.gz", out_rk="rk.nii.gz"): 722 """ Workflow for Diffusion Kurtosis reconstruction and for computing 723 DKI metrics. Performs a DKI reconstruction on the files by 'globing' 724 ``input_files`` and saves the DKI metrics in a directory specified by 725 ``out_dir``. 726 727 Parameters 728 ---------- 729 input_files : string 730 Path to the input volumes. This path may contain wildcards to 731 process multiple inputs at once. 732 bvalues_files : string 733 Path to the bvalues files. This path may contain wildcards to use 734 multiple bvalues files at once. 735 bvectors_files : string 736 Path to the bvalues files. This path may contain wildcards to use 737 multiple bvalues files at once. 738 mask_files : string 739 Path to the input masks. This path may contain wildcards to use 740 multiple masks at once. (default: No mask used) 741 b0_threshold : float, optional 742 Threshold used to find b0 volumes. 743 save_metrics : variable string, optional 744 List of metrics to save. 745 Possible values: fa, ga, rgb, md, ad, rd, mode, tensor, evec, eval 746 out_dir : string, optional 747 Output directory. (default current directory) 748 out_dt_tensor : string, optional 749 Name of the tensors volume to be saved. 750 out_dk_tensor : string, optional 751 Name of the tensors volume to be saved. 752 out_fa : string, optional 753 Name of the fractional anisotropy volume to be saved. 754 out_ga : string, optional 755 Name of the geodesic anisotropy volume to be saved. 756 out_rgb : string, optional 757 Name of the color fa volume to be saved. 758 out_md : string, optional 759 Name of the mean diffusivity volume to be saved. 760 out_ad : string, optional 761 Name of the axial diffusivity volume to be saved. 762 out_rd : string, optional 763 Name of the radial diffusivity volume to be saved. 764 out_mode : string, optional 765 Name of the mode volume to be saved. 766 out_evec : string, optional 767 Name of the eigenvectors volume to be saved. 768 out_eval : string, optional 769 Name of the eigenvalues to be saved. 770 out_mk : string, optional 771 Name of the mean kurtosis to be saved. 772 out_ak : string, optional 773 Name of the axial kurtosis to be saved. 774 out_rk : string, optional 775 Name of the radial kurtosis to be saved. 776 777 References 778 ---------- 779 780 .. [1] Tabesh, A., Jensen, J.H., Ardekani, B.A., Helpern, J.A., 2011. 781 Estimation of tensors and tensor-derived measures in diffusional 782 kurtosis imaging. Magn Reson Med. 65(3), 823-836 783 784 .. [2] Jensen, Jens H., Joseph A. Helpern, Anita Ramani, Hanzhang Lu, 785 and Kyle Kaczynski. 2005. Diffusional Kurtosis Imaging: The 786 Quantification of Non-Gaussian Water Diffusion by Means of Magnetic 787 Resonance Imaging. MRM 53 (6):1432-40. 788 """ 789 io_it = self.get_io_iterator() 790 791 for (dwi, bval, bvec, mask, otensor, ofa, oga, orgb, omd, oad, orad, 792 omode, oevecs, oevals, odk_tensor, omk, oak, ork) in io_it: 793 794 logging.info('Computing DKI metrics for {0}'.format(dwi)) 795 data, affine = load_nifti(dwi) 796 797 if mask is not None: 798 mask = load_nifti_data(mask).astype(bool) 799 800 dkfit, _ = self.get_fitted_tensor(data, mask, bval, bvec, 801 b0_threshold) 802 803 if not save_metrics: 804 save_metrics = ['mk', 'rk', 'ak', 'fa', 'md', 'rd', 'ad', 'ga', 805 'rgb', 'mode', 'evec', 'eval', 'dt_tensor', 806 'dk_tensor'] 807 808 evals, evecs, kt = split_dki_param(dkfit.model_params) 809 FA = fractional_anisotropy(evals) 810 FA[np.isnan(FA)] = 0 811 FA = np.clip(FA, 0, 1) 812 813 if 'dt_tensor' in save_metrics: 814 tensor_vals = lower_triangular(dkfit.quadratic_form) 815 correct_order = [0, 1, 3, 2, 4, 5] 816 tensor_vals_reordered = tensor_vals[..., correct_order] 817 save_nifti(otensor, tensor_vals_reordered.astype(np.float32), 818 affine) 819 820 if 'dk_tensor' in save_metrics: 821 save_nifti(odk_tensor, dkfit.kt.astype(np.float32), affine) 822 823 if 'fa' in save_metrics: 824 save_nifti(ofa, FA.astype(np.float32), affine) 825 826 if 'ga' in save_metrics: 827 GA = geodesic_anisotropy(dkfit.evals) 828 save_nifti(oga, GA.astype(np.float32), affine) 829 830 if 'rgb' in save_metrics: 831 RGB = color_fa(FA, dkfit.evecs) 832 save_nifti(orgb, np.array(255 * RGB, 'uint8'), affine) 833 834 if 'md' in save_metrics: 835 MD = mean_diffusivity(dkfit.evals) 836 save_nifti(omd, MD.astype(np.float32), affine) 837 838 if 'ad' in save_metrics: 839 AD = axial_diffusivity(dkfit.evals) 840 save_nifti(oad, AD.astype(np.float32), affine) 841 842 if 'rd' in save_metrics: 843 RD = radial_diffusivity(dkfit.evals) 844 save_nifti(orad, RD.astype(np.float32), affine) 845 846 if 'mode' in save_metrics: 847 MODE = get_mode(dkfit.quadratic_form) 848 save_nifti(omode, MODE.astype(np.float32), affine) 849 850 if 'evec' in save_metrics: 851 save_nifti(oevecs, dkfit.evecs.astype(np.float32), affine) 852 853 if 'eval' in save_metrics: 854 save_nifti(oevals, dkfit.evals.astype(np.float32), affine) 855 856 if 'mk' in save_metrics: 857 save_nifti(omk, dkfit.mk().astype(np.float32), affine) 858 859 if 'ak' in save_metrics: 860 save_nifti(oak, dkfit.ak().astype(np.float32), affine) 861 862 if 'rk' in save_metrics: 863 save_nifti(ork, dkfit.rk().astype(np.float32), affine) 864 865 logging.info('DKI metrics saved in {0}'. 866 format(os.path.dirname(oevals))) 867 868 def get_dki_model(self, gtab): 869 return DiffusionKurtosisModel(gtab) 870 871 def get_fitted_tensor(self, data, mask, bval, bvec, b0_threshold=50): 872 logging.info('Diffusion kurtosis estimation...') 873 bvals, bvecs = read_bvals_bvecs(bval, bvec) 874 if b0_threshold < bvals.min(): 875 warn("b0_threshold (value: {0}) is too low, increase your " 876 "b0_threshold. It should be higher than the first b0 value " 877 "({1}).".format(b0_threshold, bvals.min())) 878 879 gtab = gradient_table(bvals, bvecs, b0_threshold=b0_threshold) 880 dkmodel = self.get_dki_model(gtab) 881 dkfit = dkmodel.fit(data, mask) 882 883 return dkfit, gtab 884 885 886class ReconstIvimFlow(Workflow): 887 @classmethod 888 def get_short_name(cls): 889 return 'ivim' 890 891 def run(self, input_files, bvalues_files, bvectors_files, mask_files, 892 split_b_D=400, split_b_S0=200, b0_threshold=0, save_metrics=[], 893 out_dir='', out_S0_predicted='S0_predicted.nii.gz', 894 out_perfusion_fraction='perfusion_fraction.nii.gz', 895 out_D_star='D_star.nii.gz', out_D='D.nii.gz'): 896 """ Workflow for Intra-voxel Incoherent Motion reconstruction and for 897 computing IVIM metrics. Performs a IVIM reconstruction on the files 898 by 'globing' ``input_files`` and saves the IVIM metrics in a directory 899 specified by ``out_dir``. 900 901 Parameters 902 ---------- 903 input_files : string 904 Path to the input volumes. This path may contain wildcards to 905 process multiple inputs at once. 906 bvalues_files : string 907 Path to the bvalues files. This path may contain wildcards to use 908 multiple bvalues files at once. 909 bvectors_files : string 910 Path to the bvalues files. This path may contain wildcards to use 911 multiple bvalues files at once. 912 mask_files : string 913 Path to the input masks. This path may contain wildcards to use 914 multiple masks at once. (default: No mask used) 915 split_b_D : int, optional 916 Value to split the bvals to estimate D for the two-stage process of 917 fitting. 918 split_b_S0 : int, optional 919 Value to split the bvals to estimate S0 for the two-stage process 920 of fitting. 921 b0_threshold : int, optional 922 Threshold value for the b0 bval. 923 save_metrics : variable string, optional 924 List of metrics to save. 925 Possible values: S0_predicted, perfusion_fraction, D_star, D 926 out_dir : string, optional 927 Output directory. (default current directory) 928 out_S0_predicted : string, optional 929 Name of the S0 signal estimated to be saved. 930 out_perfusion_fraction : string, optional 931 Name of the estimated volume fractions to be saved. 932 out_D_star : string, optional 933 Name of the estimated pseudo-diffusion parameter to be saved. 934 out_D : string, optional 935 Name of the estimated diffusion parameter to be saved. 936 937 References 938 ---------- 939 940 .. [Stejskal65] Stejskal, E. O.; Tanner, J. E. (1 January 1965). 941 "Spin Diffusion Measurements: Spin Echoes in the 942 Presence of a Time-Dependent Field Gradient". The 943 Journal of Chemical Physics 42 (1): 288. 944 Bibcode: 1965JChPh..42..288S. doi:10.1063/1.1695690. 945 946 .. [LeBihan84] Le Bihan, Denis, et al. "Separation of diffusion 947 and perfusion in intravoxel incoherent motion MR 948 imaging." Radiology 168.2 (1988): 497-505. 949 """ 950 951 io_it = self.get_io_iterator() 952 953 for (dwi, bval, bvec, mask, oS0_predicted, operfusion_fraction, 954 oD_star, oD) in io_it: 955 956 logging.info('Computing IVIM metrics for {0}'.format(dwi)) 957 data, affine = load_nifti(dwi) 958 959 if mask is not None: 960 mask = load_nifti_data(mask).astype(bool) 961 962 ivimfit, _ = self.get_fitted_ivim(data, mask, bval, bvec, 963 b0_threshold) 964 965 if not save_metrics: 966 save_metrics = ['S0_predicted', 'perfusion_fraction', 'D_star', 967 'D'] 968 969 if 'S0_predicted' in save_metrics: 970 save_nifti(oS0_predicted, 971 ivimfit.S0_predicted.astype(np.float32), affine) 972 973 if 'perfusion_fraction' in save_metrics: 974 save_nifti(operfusion_fraction, 975 ivimfit.perfusion_fraction.astype(np.float32), 976 affine) 977 978 if 'D_star' in save_metrics: 979 save_nifti(oD_star, ivimfit.D_star.astype(np.float32), affine) 980 981 if 'D' in save_metrics: 982 save_nifti(oD, ivimfit.D.astype(np.float32), affine) 983 984 logging.info('IVIM metrics saved in {0}'. 985 format(os.path.dirname(oD))) 986 987 def get_fitted_ivim(self, data, mask, bval, bvec, b0_threshold=50): 988 logging.info('Intra-Voxel Incoherent Motion Estimation...') 989 bvals, bvecs = read_bvals_bvecs(bval, bvec) 990 if b0_threshold < bvals.min(): 991 warn("b0_threshold (value: {0}) is too low, increase your " 992 "b0_threshold. It should be higher than the first b0 value " 993 "({1}).".format(b0_threshold, bvals.min())) 994 995 gtab = gradient_table(bvals, bvecs, b0_threshold=b0_threshold) 996 ivimmodel = IvimModel(gtab) 997 ivimfit = ivimmodel.fit(data, mask) 998 999 return ivimfit, gtab 1000