scilpy.reconst package

scilpy.reconst.aodf module

scilpy.reconst.aodf.compute_asymmetry_index(sh_coeffs, order, mask)

Compute asymmetry index (ASI) (Cetin-Karayumak et al, 2018) from asymmetric ODF volume expressed in full SH basis.

Parameters:

• sh_coeffs (ndarray (x, y, z, ncoeffs)) – Input spherical harmonics coefficients.

• order (int > 0) – Max SH order.

• mask (ndarray (x, y, z), bool) – Mask inside which ASI should be computed.

Returns:

asi_map – Asymmetry index map.

Return type:

ndarray (x, y, z)

scilpy.reconst.aodf.compute_odd_power_map(sh_coeffs, order, mask)

Compute odd-power map (Poirier et al, 2020) from asymmetric ODF volume expressed in full SH basis.

Parameters:

• sh_coeffs (ndarray (x, y, z, ncoeffs)) – Input spherical harmonics coefficients.

• order (int > 0) – Max SH order.

• mask (ndarray (x, y, z), bool) – Mask inside which odd-power map should be computed.

Returns:

odd_power_map – Odd-power map.

Return type:

ndarray (x, y, z)

scilpy.reconst.bingham module

class scilpy.reconst.bingham.BinghamDistribution(f0, mu_prime1, mu_prime2)

Bases: object

Scaled Bingham distribution, given by:

B(u) = f0 * e^(- k1 * (mu1 * u)**2 - k2 * (mu2 * u)**2), mu1 and mu2 are unit vectors.

Parameters:

• f0 (float) – Scaling parameter of the distribution.

• mu_prime1 (ndarray (3,)) – Axis with highest concentration scaled by the concentration parameter k1.

• mu_prime2 (ndarray (3,)) – Axis with lowest concentration scaled by the concentration parameter k2.

evaluate(vertices)

get_flatten()

peak_direction()

scilpy.reconst.bingham.bingham_fit_sh(sh, max_lobes=5, abs_th=0.0, rel_th=0.0, min_sep_angle=25.0, max_fit_angle=15, mask=None, nbr_processes=None)

Approximate SH field by fitting Bingham distributions to up to max_lobes lobes per voxel, sorted in descending order by the amplitude of their peak direction.

Parameters:

• sh (ndarray (X, Y, Z, ncoeffs)) – SH coefficients array.

• max_lobes (unsigned int, optional) – Maximum number of lobes to fit per voxel.

• abs_th (float, optional) – Absolute threshold for peak extraction.

• rel_th (float, optional) – Relative threshold for peak extraction in the range [0, 1].

• min_sep_angle (float, optional) – Minimum separation angle between two adjacent peaks in degrees.

• max_fit_angle (float, optional) – The maximum distance in degrees around a peak direction for fitting the Bingham function.

• mask (ndarray (X, Y, Z), optional) – Mask to apply to the data.

• nbr_processes (unsigned int, optional) – The number of processes to use. If None, than multiprocessing.cpu_count() processes are executed.

Returns:

out – Bingham functions array.

Return type:

ndarray (X, Y, Z, max_lobes*9)

scilpy.reconst.bingham.bingham_to_peak_direction(bingham_volume)

Compute peak direction for each lobe for a given Bingham volume.

Parameters:

bingham_volume (ndarray (..., max_lobes, 9)) – Bingham volume.

Returns:

peak_dir – Peak direction image.

Return type:

ndarray (…, max_lobes, 3)

scilpy.reconst.bingham.bingham_to_sf(bingham_volume, vertices)

Convert a Bingham distributions volume to a spherical function.

Parameters:

• bingham_volume (Bingham parameters volume.) – A Bingham distributions volume.

• vertices (ndarray (n_vertices, 3)) – Sampling directions.

Returns:

sf – Bingham distribution evaluated at vertices.

Return type:

ndarray (…, n_lobes, N_PARAMS)

scilpy.reconst.bingham.compute_fiber_density(bingham, m=50, mask=None, nbr_processes=None)

Compute fiber density for each lobe for a given Bingham volume.

Fiber density (FD) is given by integrating the Bingham function over the sphere. Its unit is in 1/mm**3.

Parameters:

• bingham (ndarray (X, Y, Z, max_lobes*9)) – Input Bingham volume.

• m (unsigned int, optional) – Number of steps along theta axis for the integration. The number of steps along the phi axis is 2*m.

• mask (ndarray (X, Y, Z), optional) – Mask to apply to the computation.

• nbr_processes (unsigned int, optional) – The number of processes to use. If None, then multithreading.cpu_count() processes are launched.

Returns:

res – FD per lobe for each voxel.

Return type:

ndarray (X, Y, Z, max_lobes)

scilpy.reconst.bingham.compute_fiber_fraction(fd)

Compute the fiber fraction for each lobe at each voxel.

The fiber fraction (FF) represents the fraction of the current lobe’s FD on the total FD for all lobes. For each voxel, the FF sums to 1.

Parameters:

fd (ndarray (X, Y, Z, max_lobes)) – Fiber density image.

Returns:

ff – Fiber fraction image.

Return type:

ndarray (X, Y, Z, max_lobes)

scilpy.reconst.bingham.compute_fiber_spread(binghams, fd)

Compute fiber spread for each lobe for a given Bingham volume.

Fiber spread (FS) characterizes the spread of the lobe. The higher FS is, the wider the lobe. The unit of the FS is radians.

Parameters:

• binghams (ndarray (X, Y, Z, max_lobes*9)) – Bingham volume.

• fd (ndarray (X, Y, Z, max_lobes)) – Fiber density image.

Returns:

fs – Fiber spread image.

Return type:

ndarray (X, Y, Z, max_lobes)

scilpy.reconst.divide module

scilpy.reconst.divide.fit_gamma(data, gtab_infos, mask=None, fit_iters=1, random_iters=50, do_weight_bvals=False, do_weight_pa=False, do_multiple_s0=False, nbr_processes=None)

Fit the gamma model to data

Parameters:

• data (np.ndarray (4d)) – Diffusion data, powder averaged. Obtained as output of the function reconst.b_tensor_utils.generate_powder_averaged_data.

• gtab_infos (np.ndarray) – Contains information about the gtab, such as the unique bvals, the encoding types, the number of directions and the acquisition index. Obtained as output of the function reconst.b_tensor_utils.generate_powder_averaged_data.

• mask (np.ndarray, optional) – If mask is provided, only the data inside the mask will be used for computations.

• fit_iters (int, optional) – Number of iterations in the gamma fit. Defaults to 1.

• random_iters (int, optional) – Number of random sets of parameters tested to find the initial parameters. Defaults to 50.

• do_weight_bvals (bool , optional) – If set, does a weighting on the bvalues in the gamma fit.

• do_weight_pa (bool, optional) – If set, does a powder averaging weighting in the gamma fit.

• do_multiple_s0 (bool, optional) – If set, takes into account multiple baseline signals.

• nbr_processes (int, optional) – The number of subprocesses to use. Default: multiprocessing.cpu_count()

Returns:

fit_array – Array containing the fit

Return type:

np.ndarray

scilpy.reconst.divide.gamma_fit2metrics(params)

Compute metrics from fit parameters. This is the only function that takes the full brain.

Parameters:

params (np.ndarray) – Array containing the parameters of the fit for the whole brain.

Returns:

• microFA (np.ndarray) – MicroFA values for the whole brain.

• MK_I (np.ndarray) – Isotropic mean kurtosis values for the whole brain.

• MK_A (np.ndarray) – Anisotropic mean kurtosis values for the whole brain.

• MK_T (np.ndarray) – Total mean kurtosis values for the whole brain.

scilpy.reconst.fiber_coherence module

scilpy.reconst.fiber_coherence.compute_coherence_table_for_transforms(directions, values)

Compute fiber coherence indexes for all possible axes permutations/flips (ex, originating from a flip in the gradient table).

The mathematics are presented in : [1] Schilling et al. A fiber coherence index for quality control of B-table orientation in diffusion MRI scans. Magn Reson Imaging. 2019 May;58:82-89. doi: 10.1016/j.mri.2019.01.018.

Parameters:

• directions (ndarray (x, y, z, 3)) – Principal fiber orientation for each voxel.

• values (ndarray (x, y, z)) – Anisotropy measure for each voxel (e.g. FA map).

Returns:

• coherence (list) – Fiber coherence value for each permutation/flip.

• transforms (list) – Transform representing each permutation/flip, in the same order as coherence list.

scilpy.reconst.fiber_coherence.compute_fiber_coherence(peaks, values)

Compute the fiber coherence for peaks and values.

Parameters:

• peaks (ndarray (x, y, z, 3)) – Principal fiber orientation for each voxel.

• values (ndarray (x, y, z)) – Anisotropy measure for each voxel (e.g. FA map).

Returns:

coherence – Fiber coherence value.

Return type:

float

scilpy.reconst.fodf module

scilpy.reconst.fodf.fit_from_model(model, data, mask=None, nbr_processes=None)

Fit the model to data. Can use parallel processing.

Parameters:

• model (a model instance) – It will be used to fit the data. e.g: An instance of dipy.reconst.shm.SphHarmFit.

• data (np.ndarray (4d)) – Diffusion data.

• mask (np.ndarray, optional) – If mask is provided, only the data inside the mask will be used for computations.

• nbr_processes (int, optional) – The number of subprocesses to use. Default: multiprocessing.cpu_count()

Returns:

fit_array – Dipy’s MultiVoxelFit, containing the fit. It contains an array of fits. Any attributes of its individuals fits (of class given by ‘model.fit’) can be accessed through the MultiVoxelFit to get all fits at once.

Return type:

MultiVoxelFit

scilpy.reconst.fodf.get_ventricles_max_fodf(data, fa, md, zoom, sh_basis, fa_threshold, md_threshold, small_dims=False, is_legacy=True)

Compute mean maximal fodf value in ventricules. Given heuristics thresholds on FA and MD values, finds the voxels of the ventricules or CSF and computes a mean fODF value. This is described in Dell’Acqua et al. HBM 2013.

Ventricles are searched in a window in the middle of the data to increase speed. No need to scan the whole image.

Parameters:

• data (ndarray (x, y, z, ncoeffs)) – Input fODF file in spherical harmonics coefficients. Uses sphere ‘repulsion100’ to convert to SF values.

• fa (ndarray (x, y, z)) – FA (Fractional Anisotropy) volume from DTI

• md (ndarray (x, y, z)) – MD (Mean Diffusivity) volume from DTI

• zoom (List of length 3) – The resolution. A total number of voxels of 1000 works well at 2x2x2 = 8 mm3.

• sh_basis (str) – Either ‘tournier07’ or ‘descoteaux07’

• small_dims (bool) – If set, takes the full range of data to search the max fodf amplitude in ventricles, rather than a window center in the data. Useful when the data has small dimensions.

• fa_threshold (float) – Maximal threshold of FA (voxels under that threshold are considered for evaluation). Suggested value: 0.1.

• md_threshold (float) – Minimal threshold of MD in mm2/s (voxels above that threshold are considered for evaluation). Suggested value: 0.003.

• is_legacy (bool, optional) – Whether the SH basis is in its legacy form.

Returns:

mean, mask – Mean maximum fODF value and mask of voxels used.

Return type:

int, ndarray (x, y, z)

scilpy.reconst.fodf.verify_failed_voxels_shm_coeff(shm_coeff)

Verifies if there are any NaN in the final coefficients, and if so raises warnings.

Parameters:

shm_coeff (np.ndarray) – The shm_coeff given by dipy’s fit classes. Of shape (x, y, z, n) with the coefficients on the last dimension.

Returns:

shm_coeff – The coefficients with 0 instead of NaNs.

Return type:

np.ndarray

scilpy.reconst.fodf.verify_frf_files(wm_frf, gm_frf, csf_frf)

Verifies that all three frf files contain four columns, else raises ValueErrors.

Parameters:

• wm_frf (np.ndarray) – The frf directly as loaded from its text file.

• gm_frf (np.ndarray) – Idem

• csf_frf (np.ndarray) – Idem

Returns:

• wm_frf (np.ndarray) – The file. In the case where there was only one line in the file, and it has been loaded as a vector, we return the array formatted as 2D, with the 4 frf values as columns.

• gm_frf (np.ndarray) – Idem

• csf_frf (np.ndarray) – Idem

scilpy.reconst.frf module

scilpy.reconst.frf.compute_msmt_frf(data, bvals, bvecs, btens=None, data_dti=None, bvals_dti=None, bvecs_dti=None, btens_dti=None, mask=None, mask_wm=None, mask_gm=None, mask_csf=None, fa_thr_wm=0.7, fa_thr_gm=0.2, fa_thr_csf=0.1, md_thr_gm=0.0007, md_thr_csf=0.003, min_nvox=300, roi_radii=10, roi_center=None, tol=20)

Computes a multi-shell, multi-tissue single Fiber Response Function from a DWI volume. A DTI fit is made, and voxels containing a single fiber population are found using a threshold on the FA and MD, inside a mask of each tissue type.

Parameters:

• data (ndarray) – 4D Input diffusion volume with shape (X, Y, Z, N)

• bvals (ndarray) – 1D bvals array with shape (N,)

• bvecs (ndarray) – 2D bvecs array with shape (N, 3)

• btens (ndarray 1D) – btens array with shape (N,), describing the btensor shape of every pair of bval/bvec.

• data_dti (ndarray 4D) – Input diffusion volume with shape (X, Y, Z, M), where M is the number of DTI directions.

• bvals_dti (ndarray 1D) – bvals array with shape (M,)

• bvecs_dti (ndarray 2D) – bvals array with shape (M, 3)

• btens_dti (ndarray 1D) – bvals array with shape (M,)

• mask (ndarray, optional) – 3D mask with shape (X,Y,Z) Binary mask. Only the data inside the mask will be used for computations and reconstruction.

• mask_wm (ndarray, optional) – 3D mask with shape (X,Y,Z) Binary white matter mask. Only the data inside this mask will be used to estimate the fiber response function of WM.

• mask_gm (ndarray, optional) – 3D mask with shape (X,Y,Z) Binary grey matter mask. Only the data inside this mask will be used to estimate the fiber response function of GM.

• mask_csf (ndarray, optional) – 3D mask with shape (X,Y,Z) Binary csf mask. Only the data inside this mask will be used to estimate the fiber response function of CSF.

• fa_thr_wm (float, optional) – Use this threshold to select single WM fiber voxels from the FA inside the WM mask defined by mask_wm. Each voxel above this threshold will be selected. Defaults to 0.7

• fa_thr_gm (float, optional) – Use this threshold to select GM voxels from the FA inside the GM mask defined by mask_gm. Each voxel below this threshold will be selected. Defaults to 0.2

• fa_thr_csf (float, optional) – Use this threshold to select CSF voxels from the FA inside the CSF mask defined by mask_csf. Each voxel below this threshold will be selected. Defaults to 0.1

• md_thr_gm (float, optional) – Use this threshold to select GM voxels from the MD inside the GM mask defined by mask_gm. Each voxel below this threshold will be selected. Defaults to 0.0007

• md_thr_csf (float, optional) – Use this threshold to select CSF voxels from the MD inside the CSF mask defined by mask_csf. Each voxel below this threshold will be selected. Defaults to 0.003

• min_nvox (int, optional) – Minimal number of voxels needing to be identified as single fiber voxels in the automatic estimation. Defaults to 300.

• roi_radii (int or array-like (3,), optional) – Use those radii to select a cuboid roi to estimate the FRF. The roi will be a cuboid spanning from the middle of the volume in each direction with the different radii. Defaults to 10.

• roi_center (tuple(3), optional) – Use this center to span the roi of size roi_radius (center of the 3D volume).

• tol (int) – tolerance gap for b-values clustering. Defaults to 20

Returns:

• reponses (list of ndarray) – Fiber Response Function of each (3) tissue type, with shape (4, N).

• frf_masks (list of ndarray) – Mask where the frf was calculated, for each (3) tissue type, with shape (X, Y, Z).

Raises:

ValueError – If less than min_nvox voxels were found with sufficient FA to estimate the FRF.

scilpy.reconst.frf.compute_ssst_frf(data, bvals, bvecs, b0_threshold=20, mask=None, mask_wm=None, fa_thresh=0.7, min_fa_thresh=0.5, min_nvox=300, roi_radii=10, roi_center=None)

Computes a single-shell (under b=1500), single-tissue single Fiber Response Function from a DWI volume. A DTI fit is made, and voxels containing a single fiber population are found using either a threshold on the FA, inside a white matter mask.

Parameters:

• data (ndarray) – 4D Input diffusion volume with shape (X, Y, Z, N)

• bvals (ndarray) – 1D bvals array with shape (N,)

• bvecs (ndarray) – 2D bvecs array with shape (N, 3)

• b0_threshold (float) – Value under which bvals are considered b0s.

• mask (ndarray, optional) – 3D mask with shape (X,Y,Z) Binary mask. Only the data inside the mask will be used for computations and reconstruction. Useful if no white matter mask is available.

• mask_wm (ndarray, optional) – 3D mask with shape (X,Y,Z) Binary white matter mask. Only the data inside this mask and above the threshold defined by fa_thresh will be used to estimate the fiber response function. If not given, all voxels inside mask will be used.

• fa_thresh (float, optional) – Use this threshold as the initial threshold to select single fiber voxels. Defaults to 0.7

• min_fa_thresh (float, optional) – Minimal value that will be tried when looking for single fiber voxels. Defaults to 0.5

• min_nvox (int, optional) – Minimal number of voxels needing to be identified as single fiber voxels in the automatic estimation. Defaults to 300.

• roi_radii (int or array-like (3,), optional) – Use those radii to select a cuboid roi to estimate the FRF. The roi will be a cuboid spanning from the middle of the volume in each direction with the different radii. Defaults to 10.

• roi_center (tuple(3), optional) – Use this center to span the roi of size roi_radius (center of the 3D volume).

Returns:

full_response – Fiber Response Function, with shape (4,)

Return type:

ndarray

Raises:

ValueError – If less than min_nvox voxels were found with sufficient FA to estimate the FRF.

scilpy.reconst.frf.replace_frf(old_frf, new_frf, no_factor=False)

Replaces the 3 first values of old_frf with new_frf. Formats the new_frf from a string value and verifies that the number of shells corresponds.

Parameters:

• old_frf (np.ndarray) – A loaded frf file, of shape (N, 4), where N is the number of shells.

• new_frf (str) – The new frf, to be interpreted with a 10**-4 factor. Ex: 15,4,4. With multishell: all values, concatenated into one string. Ex: 15,4,4,13,5,5,12,5,5.

• no_factor (bool) – If true, the fiber response function is evaluated without the 10**-4 factor.

Returns:

response – Formatted new frf, of shape (n, 4)

Return type:

np.ndarray

scilpy.reconst.mti module

scilpy.reconst.mti.adjust_B1_map_header(B1_img, slope)

Fixe issue in B1 map header, by applying the scaling (slope) and setting the slope to 1.

Parameters:

• B1_img (nifti image object) – The B1 map.

• slope (float) – The slope value, obtained from the image header.

Return type:

Adjusted B1 nifti image object (new_b1_img)

scilpy.reconst.mti.adjust_B1_map_intensities(B1_map, nominal=100)

Adjust and verify the values of the B1 map. We want the B1 map to have values around 1.0, so we divide by the nominal B1 value if it is not already the case. Sometimes the scaling of the B1 map is wrong, leading to weird values after this process, which what is verified here.

Parameters:

• B1_map (B1 coregister map.)

• nominal (Nominal values of B1. For Philips, should be 1.)

Return type:

Ajusted B1 map in 3d-array.

scilpy.reconst.mti.apply_B1_corr_empiric(MT_map, B1_map)

Apply an empiric B1 correction on MT map.

see Weiskopf et al., 2013 https://www.frontiersin.org/articles/10.3389/fnins.2013.00095/full

Parameters:

• MT_map (3D-Array MT map.)

• B1_map (B1 coregister map.)

Return type:

Corrected MT matrix in 3D-array.

scilpy.reconst.mti.apply_B1_corr_model_based(MTsat, B1_map, R1app, fitvalues_file, B1_ref=1)

Apply a model-based B1 correction on MT map.

see Rowley et al., 2021 https://onlinelibrary.wiley.com/doi/10.1002/mrm.28831

Parameters:

• MTsat (3D-Array MTsat map.)

• B1_map (B1 coregister map.)

• R1app (Apparent R1 map, obtained from compute_saturation_map.)

• fitvalues_file (Path to the fitValues_*.mat file corresponding to the) – MTsat map. This file is computed with Christopher Rowley’s Matlab library.

• B1_ref (Reference B1 value used for the fit (usually 1).)

Return type:

Corrected MTsat matrix in 3D-array.

scilpy.reconst.mti.compute_ratio_map(mt_on_single, mt_off, mt_on_dual=None)

Compute magnetization transfer ratio (MTR), and inhomogenous magnetization transfer ratio (ihMTR) if mt_on_dual is given. - MT ratio (MTR) is computed as the percentage difference of single frequency mean image from the reference. - ihMT ratio (ihMTR) is computed as the percentage difference of dual from single frequency.

For ihMTR, see Varma et al., 2015 www.sciencedirect.com/science/article/pii/S1090780715001998

Parameters:

• mt_on_single (MT-on image with single frequency pulse: can be one) – single positive/negative frequency MT image or the average of both images (positive and negative).

• mt_off (MT-off image: the reference image without MT) – preparation.

• mt_on_dual (MT-on image with dual frequency pulse: can be one) – dual alternating positive/negative frequency MT image or the average of both images. Optional. If given, will compute the ihMTR also.

Return type:

MT ratio (MTR), ihMT ratio (ihMTR).

scilpy.reconst.mti.compute_saturation_map(MT, PD, T1, a, TR)

Compute saturation of given contrast (MT). Saturation is computed from apparent longitudinal relaxation time (T1app) and apparent signal amplitude (Aapp), as the difference of longitudinal relaxation times of bound to macromolecules pool from free water pool saturation. The estimation of the saturation includes correction for the effects of excitation flip angle and longitudinal relaxation time, and remove the effect of T1-weighted image.

see Helms et al., 2008 https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.21732

Parameters:

• MT (Contrast of choice (can be a single contrast or the) – mean of positive and negative for example)

• PD (Reference proton density weighted image (MToff PD))

• T1 (Reference T1 weighted image (MToff T1))

• a (List of two flip angles corresponding to PD and T1. If B1) – correction is on, this should be two 3D-array of flip angles varying spatially with respect to the B1 map.

• TR (List of two repetition times corresponding to PD and T1.)

Return type:

MT saturation matrice in 3D-array (sat), computed apparent T1 map (T1app).

scilpy.reconst.mti.process_contrast_map(merged_images, single_echo=False, filtering=False)

Average echoes of a contrast map and apply gaussian filter.

Parameters:

• 4d_image (4D image, containing every echoes of the contrast map)

• single_echo (Apply when only one echo is used to compute contrast map)

• filtering (Apply Gaussian filtering to remove Gibbs ringing) – (default is False).

Return type:

Contrast map in 3D-Array.

scilpy.reconst.mti.py_fspecial_gauss(shape, sigma)

Function to mimic the ‘fspecial gaussian’ MATLAB function Returns a rotationally symmetric Gaussian lowpass filter of shape size with standard deviation sigma. see https://www.mathworks.com/help/images/ref/fspecial.html

Parameters:

• shape (Vector to specify the size of square matrix (rows, columns).) – Ex.: (3, 3)

• sigma (Value for standard deviation (Sigma value of 0.5 is) – recommended).

Return type:

Two-dimensional Gaussian filter h of specified size.

scilpy.reconst.mti.smooth_B1_map(B1_map, wdims=5)

Apply a smoothing to the B1 map.

Parameters:

• B1_map (B1 coregister map.)

• wdims (Window dimension (in voxels) for the smoothing.)

Return type:

Smoothed B1 map in 3d-array.

scilpy.reconst.mti.threshold_map(computed_map, in_mask, lower_threshold, upper_threshold, idx_contrast_list=None, contrast_maps=None)

Remove NaN and apply different threshold based on - maximum and minimum threshold value - T1 mask - combination of specific contrast maps

idx_contrast_list and contrast_maps are required for thresholding of ihMT images.

Parameters:

• computed_map (3D-Array data.) – Myelin map (ihMT or non-ihMT maps)

• in_mask (Path to binary T1 mask from T1 segmentation.) – Must be the sum of GM+WM+CSF.

• lower_threshold (Value for low thresold <int>)

• upper_thresold (Value for up thresold <int>)

• idx_contrast_list (List of indexes of contrast maps corresponding to) – that of input contrast_maps ex.: [0, 2, 4] Altnp = 0; Atlpn = 1; Negative = 2; Positive = 3; PD = 4; T1 = 5

• contrast_maps (List of 3D-Array. File must containing the) – 5 or 6 contrast maps.

Return type:

Thresholded matrix in 3D-array.

scilpy.reconst.sh module

scilpy.reconst.sh.compute_rish(sh, mask=None, full_basis=False)

Compute the RISH (Rotationally Invariant Spherical Harmonics) features of the SH signal [1]. Each RISH feature map is the total energy of its associated order. Mathematically, it is the sum of the squared SH coefficients of the SH order.

Parameters:

• sh (np.ndarray object) – Array of the SH coefficients

• mask (np.ndarray object, optional) – Binary mask. Only data inside the mask will be used for computation.

• full_basis (bool, optional) – True when coefficients are for a full SH basis.

Returns:

• rish (np.ndarray with shape (x,y,z,n_orders)) – The RISH features of the input SH, with one channel per SH order.

• orders (list(int)) – The SH order of each RISH feature in the last channel of rish.

References

[1] Mirzaalian, Hengameh, et al. “Harmonizing diffusion MRI data across

multiple sites and scanners.” MICCAI 2015. https://scholar.harvard.edu/files/hengameh/files/miccai2015.pdf

scilpy.reconst.sh.compute_sh_coefficients(dwi, gradient_table, b0_threshold=20, sh_order=4, basis_type='descoteaux07', smooth=0.006, use_attenuation=False, mask=None, sphere=None, is_legacy=True)

Fit a diffusion signal with spherical harmonics coefficients. Data must come from a single shell acquisition.

Parameters:

• dwi (nib.Nifti1Image object) – Diffusion signal as weighted images (4D).

• gradient_table (GradientTable) – Dipy object that contains all bvals and bvecs.

• b0_threshold (float) – Threshold for the b0 values. Used to validate that the data contains single shell signal.

• sh_order (int, optional) – SH order to fit, by default 4.

• basis_type (str) – Either ‘tournier07’ or ‘descoteaux07’

• smooth (float, optional) – Lambda-regularization coefficient in the SH fit, by default 0.006.

• use_attenuation (bool, optional) – If true, we will use DWI attenuation. [False]

• mask (nib.Nifti1Image object, optional) – Binary mask. Only data inside the mask will be used for computations and reconstruction.

• sphere (Sphere) – Dipy object. If not provided, will use Sphere(xyz=bvecs).

• is_legacy (bool, optional) – Whether or not the SH basis is in its legacy form.

Returns:

sh_coeffs – Spherical harmonics coefficients at every voxel. The actual number of coefficients depends on sh_order.

Return type:

np.ndarray with shape (X, Y, Z, #coeffs)

scilpy.reconst.sh.convert_sh_basis(shm_coeff, sphere, mask=None, input_basis='descoteaux07', output_basis='tournier07', is_input_legacy=True, is_output_legacy=False, nbr_processes=None)

Converts spherical harmonic coefficients between two bases

Parameters:

• shm_coeff (np.ndarray) – Spherical harmonic coefficients

• sphere (Sphere) – The Sphere providing discrete directions for evaluation.

• mask (np.ndarray, optional) – If mask is provided, only the data inside the mask will be used for computations.

• input_basis (str, optional) – Type of spherical harmonic basis used for shm_coeff. Either descoteaux07 or tournier07. Default: descoteaux07

• output_basis (str, optional) – Type of spherical harmonic basis wanted as output. Either descoteaux07 or tournier07. Default: tournier07

• is_input_legacy (bool, optional) – If true, this means that the input SH used a legacy basis definition for backward compatibility with previous tournier07 and descoteaux07 implementations. Default: True

• is_output_legacy (bool, optional) – If true, this means that the output SH will use a legacy basis definition for backward compatibility with previous tournier07 and descoteaux07 implementations. Default: False

• nbr_processes (int, optional) – The number of subprocesses to use. Default: multiprocessing.cpu_count()

Returns:

shm_coeff_array – Spherical harmonic coefficients in the desired basis.

Return type:

np.ndarray

scilpy.reconst.sh.convert_sh_to_sf(shm_coeff, sphere, mask=None, dtype='float32', input_basis='descoteaux07', input_full_basis=False, is_input_legacy=True, nbr_processes=2)

Converts spherical harmonic coefficients to an SF sphere

Parameters:

• shm_coeff (np.ndarray) – Spherical harmonic coefficients

• sphere (Sphere) – The Sphere providing discrete directions for evaluation.

• mask (np.ndarray, optional) – If mask is provided, only the data inside the mask will be used for computations.

• dtype (str) – Datatype to use for computation and output array. Either float32 or float64. Default: float32

• input_basis (str, optional) – Type of spherical harmonic basis used for shm_coeff. Either descoteaux07 or tournier07. Default: descoteaux07

• input_full_basis (bool, optional) – If True, use a full SH basis (even and odd orders) for the input SH coefficients.

• is_input_legacy (bool, optional) – Whether the input basis is in its legacy form.

• nbr_processes (int, optional) – The number of subprocesses to use. Default: multiprocessing.cpu_count()

Returns:

shm_coeff_array – Spherical harmonic coefficients in the desired basis.

Return type:

np.ndarray

scilpy.reconst.sh.maps_from_sh(shm_coeff, peak_dirs, peak_values, peak_indices, sphere, mask=None, gfa_thr=0, sh_basis_type='descoteaux07', nbr_processes=None)

Computes maps from given SH coefficients and peaks

Parameters:

• shm_coeff (np.ndarray) – Spherical harmonic coefficients

• peak_dirs (np.ndarray) – Peak directions

• peak_values (np.ndarray) – Peak values

• peak_indices (np.ndarray) – Peak indices

• sphere (Sphere) – The Sphere providing discrete directions for evaluation.

• mask (np.ndarray, optional) – If mask is provided, only the data inside the mask will be used for computations.

• gfa_thr (float, optional) – Voxels with gfa less than gfa_thr are skipped for all metrics, except rgb_map. Default: 0

• sh_basis_type (str, optional) – Type of spherical harmonic basis used for shm_coeff. Either descoteaux07 or tournier07. Default: descoteaux07

• nbr_processes (int, optional) – The number of subprocesses to use. Default: multiprocessing.cpu_count()

Returns:

nufo_map, afd_max, afd_sum, rgb_map, gfa, qa

Return type:

tuple of np.ndarray

scilpy.reconst.sh.peaks_from_sh(shm_coeff, sphere, mask=None, relative_peak_threshold=0.5, absolute_threshold=0, min_separation_angle=25, normalize_peaks=False, npeaks=5, sh_basis_type='descoteaux07', is_legacy=True, nbr_processes=None, full_basis=False, is_symmetric=True)

Computes peaks from given spherical harmonic coefficients

Parameters:

• shm_coeff (np.ndarray) – Spherical harmonic coefficients

• sphere (Sphere) – The Sphere providing discrete directions for evaluation.

• mask (np.ndarray, optional) – If mask is provided, only the data inside the mask will be used for computations.

• relative_peak_threshold (float, optional) – Only return peaks greater than relative_peak_threshold * m where m is the largest peak. Default: 0.5

• absolute_threshold (float, optional) – Absolute threshold on fODF amplitude. This value should be set to approximately 1.5 to 2 times the maximum fODF amplitude in isotropic voxels (ex. ventricles). scil_fodf_max_in_ventricles.py can be used to find the maximal value. Default: 0

• min_separation_angle (float in [0, 90], optional) – The minimum distance between directions. If two peaks are too close only the larger of the two is returned. Default: 25

• normalize_peaks (bool, optional) – If true, all peak values are calculated relative to max(odf).

• npeaks (int, optional) – Maximum number of peaks found (default 5 peaks).

• sh_basis_type (str, optional) – Type of spherical harmonic basis used for shm_coeff. Either descoteaux07 or tournier07. Default: descoteaux07

• is_legacy (bool, optional) – If true, this means that the input SH used a legacy basis definition for backward compatibility with previous tournier07 and descoteaux07 implementations. Default: True

• nbr_processes (int, optional) – The number of subprocesses to use. Default: multiprocessing.cpu_count()

• full_basis (bool, optional) – If True, SH coefficients are expressed using a full basis. Default: False

• is_symmetric (bool, optional) – If False, antipodal sphere directions are considered distinct. Default: True

Returns:

peak_dirs, peak_values, peak_indices

Return type:

tuple of np.ndarray

scilpy.reconst.sh.verify_data_vs_sh_order(data, sh_order)

Raises a warning if the dwi data shape is not enough for the chosen sh_order.

Parameters:

• data (np.ndarray) – Diffusion signal as weighted images (4D).

• sh_order (int) – SH order to fit, by default 4.

scilpy.reconst.utils module

scilpy.reconst.utils.find_order_from_nb_coeff(data)

scilpy.reconst.utils.get_maximas(data, sphere, b_matrix, threshold, absolute_threshold, min_separation_angle=25)

scilpy.reconst.utils.get_sh_order_and_fullness(ncoeffs)

Get the order of the SH basis from the number of SH coefficients as well as a boolean indicating if the basis is full.

scilpy.reconst.utils.get_sphere_neighbours(sphere, max_angle)

Get a matrix of neighbours for each direction on the sphere, within the min_separation_angle.

min_separation_angle: float

Maximum angle in radians defining the neighbourhood of each direction.

Returns:

neighbours – Neighbour directions for each direction on the sphere.

Return type:

ndarray