denoise

bench([label, verbose, extra_argv])

Run benchmarks for module using nose.

test([label, verbose, extra_argv, doctests, …])

Run tests for module using nose.

Module: denoise.adaptive_soft_matching

adaptive_soft_matching(ima, fimau, fimao, sigma)

Adaptive Soft Coefficient Matching

Module: denoise.denspeed

add_padding_reflection

correspond_indices

cpu_count

Return number of cpus as determined by omp_get_num_procs.

determine_num_threads

Determine the effective number of threads to be used for OpenMP calls

nlmeans_3d

Non-local means for denoising 3D images

remove_padding

Module: denoise.enhancement_kernel

EnhancementKernel

Methods

HemiSphere([x, y, z, theta, phi, xyz, …])

Points on the unit sphere.

Sphere([x, y, z, theta, phi, xyz, faces, edges])

Points on the unit sphere.

ceil(x, /)

Return the ceiling of x as an Integral.

disperse_charges(hemi, iters[, const])

Models electrostatic repulsion on the unit sphere

get_sphere([name])

provide triangulated spheres

gettempdir()

Accessor for tempfile.tempdir.

Module: denoise.gibbs

Version

alias of numpy.lib._version.NumpyVersion

partial

partial(func, *args, **keywords) - new function with partial application of the given arguments and keywords.

Pool([processes, initializer, initargs, …])

Returns a process pool object

deprecated_params(old_name[, new_name, …])

Deprecate a renamed or removed function argument.

determine_num_processes(num_processes)

Determine the effective number of processes for parallelization.

gibbs_removal(vol[, slice_axis, n_points, …])

Suppresses Gibbs ringing artefacts of images volumes.

Module: denoise.localpca

eigh(a[, b, lower, eigvals_only, …])

Solve a standard or generalized eigenvalue problem for a complex Hermitian or real symmetric matrix.

genpca(arr[, sigma, mask, patch_radius, …])

General function to perform PCA-based denoising of diffusion datasets.

localpca(arr, sigma[, mask, patch_radius, …])

Performs local PCA denoising according to Manjon et al.

mppca(arr[, mask, patch_radius, pca_method, …])

Performs PCA-based denoising using the Marcenko-Pastur distribution [1].

Module: denoise.nlmeans

nlmeans(arr, sigma[, mask, patch_radius, …])

Non-local means for denoising 3D and 4D images

nlmeans_3d

Non-local means for denoising 3D images

Module: denoise.nlmeans_block

firdn

Applies the filter given by the convolution kernel ‘h’ columnwise to ‘image’, then subsamples by 2.

nlmeans_block

Non-Local Means Denoising Using Blockwise Averaging

upfir

Upsamples the columns of the input image by 2, then applies the convolution kernel ‘h’ (again, columnwise).

Module: denoise.noise_estimate

convolve(input, weights[, output, mode, …])

Multidimensional convolution.

estimate_sigma(arr[, …])

Standard deviation estimation from local patches

piesno(data, N[, alpha, l, itermax, eps, …])

Probabilistic Identification and Estimation of Noise (PIESNO).

Module: denoise.non_local_means

nlmeans_block

Non-Local Means Denoising Using Blockwise Averaging

non_local_means(arr, sigma[, mask, …])

Non-local means for denoising 3D and 4D images, using

Module: denoise.patch2self

optional_package(name[, trip_msg])

Return package-like thing and module setup for package name

patch2self(data, bvals[, patch_radius, …])

Patch2Self Denoiser

warn(/, message[, category, stacklevel, source])

Issue a warning, or maybe ignore it or raise an exception.

Module: denoise.pca_noise_estimate

pca_noise_estimate

PCA based local noise estimation.

Module: denoise.shift_twist_convolution

EnhancementKernel

Methods

convolve

Perform the shift-twist convolution with the ODF data and the lookup-table of the kernel.

convolve_sf

Perform the shift-twist convolution with the ODF data and the lookup-table of the kernel.

cpu_count

Return number of cpus as determined by omp_get_num_procs.

determine_num_threads

Determine the effective number of threads to be used for OpenMP calls

get_sphere([name])

provide triangulated spheres

sf_to_sh(sf, sphere[, sh_order, basis_type, …])

Spherical function to spherical harmonics (SH).

sh_to_sf(sh, sphere[, sh_order, basis_type, …])

Spherical harmonics (SH) to spherical function (SF).

bench

dipy.denoise.bench(label='fast', verbose=1, extra_argv=None)

Run benchmarks for module using nose.

Parameters
label{‘fast’, ‘full’, ‘’, attribute identifier}, optional

Identifies the benchmarks to run. This can be a string to pass to the nosetests executable with the ‘-A’ option, or one of several special values. Special values are:

  • ‘fast’ - the default - which corresponds to the nosetests -A option of ‘not slow’.

  • ‘full’ - fast (as above) and slow benchmarks as in the ‘no -A’ option to nosetests - this is the same as ‘’.

  • None or ‘’ - run all tests.

  • attribute_identifier - string passed directly to nosetests as ‘-A’.

verboseint, optional

Verbosity value for benchmark outputs, in the range 1-10. Default is 1.

extra_argvlist, optional

List with any extra arguments to pass to nosetests.

Returns
successbool

Returns True if running the benchmarks works, False if an error occurred.

Notes

Benchmarks are like tests, but have names starting with “bench” instead of “test”, and can be found under the “benchmarks” sub-directory of the module.

Each NumPy module exposes bench in its namespace to run all benchmarks for it.

Examples

>>> success = np.lib.bench() 
Running benchmarks for numpy.lib
...
using 562341 items:
unique:
0.11
unique1d:
0.11
ratio: 1.0
nUnique: 56230 == 56230
...
OK
>>> success 
True

test

dipy.denoise.test(label='fast', verbose=1, extra_argv=None, doctests=False, coverage=False, raise_warnings=None, timer=False)

Run tests for module using nose.

Parameters
label{‘fast’, ‘full’, ‘’, attribute identifier}, optional

Identifies the tests to run. This can be a string to pass to the nosetests executable with the ‘-A’ option, or one of several special values. Special values are:

  • ‘fast’ - the default - which corresponds to the nosetests -A option of ‘not slow’.

  • ‘full’ - fast (as above) and slow tests as in the ‘no -A’ option to nosetests - this is the same as ‘’.

  • None or ‘’ - run all tests.

  • attribute_identifier - string passed directly to nosetests as ‘-A’.

verboseint, optional

Verbosity value for test outputs, in the range 1-10. Default is 1.

extra_argvlist, optional

List with any extra arguments to pass to nosetests.

doctestsbool, optional

If True, run doctests in module. Default is False.

coveragebool, optional

If True, report coverage of NumPy code. Default is False. (This requires the coverage module).

raise_warningsNone, str or sequence of warnings, optional

This specifies which warnings to configure as ‘raise’ instead of being shown once during the test execution. Valid strings are:

  • “develop” : equals (Warning,)

  • “release” : equals (), do not raise on any warnings.

timerbool or int, optional

Timing of individual tests with nose-timer (which needs to be installed). If True, time tests and report on all of them. If an integer (say N), report timing results for N slowest tests.

Returns
resultobject

Returns the result of running the tests as a nose.result.TextTestResult object.

Notes

Each NumPy module exposes test in its namespace to run all tests for it. For example, to run all tests for numpy.lib:

>>> np.lib.test() 

Examples

>>> result = np.lib.test() 
Running unit tests for numpy.lib
...
Ran 976 tests in 3.933s

OK

>>> result.errors 
[]
>>> result.knownfail 
[]

adaptive_soft_matching

dipy.denoise.adaptive_soft_matching.adaptive_soft_matching(ima, fimau, fimao, sigma)

Adaptive Soft Coefficient Matching

Combines two filtered 3D-images at different resolutions and the original image. Returns the resulting combined image.

Parameters
imathe original (not filtered) image
fimau3D double array,

filtered image with optimized non-local means using a small block (suggested:3x3), which corresponds to a “high resolution” filter.

fimao3D double array,

filtered image with optimized non-local means using a small block (suggested:5x5), which corresponds to a “low resolution” filter.

sigmathe estimated standard deviation of the Gaussian random variables

that explain the rician noise. Note: In P. Coupe et al. the rician noise was simulated as sqrt((f+x)^2 + (y)^2) where f is the pixel value and x and y are independent realizations of a random variable with Normal distribution, with mean=0 and standard deviation=h

Returns
fima3D double array

output denoised array which is of the same shape as that of the input

References

Coupe11

Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins. “Multiresolution Non-Local Means Filter for 3D MR Image Denoising” IET Image Processing, Institution of Engineering and Technology, 2011

add_padding_reflection

dipy.denoise.denspeed.add_padding_reflection()

correspond_indices

dipy.denoise.denspeed.correspond_indices()

cpu_count

dipy.denoise.denspeed.cpu_count()

Return number of cpus as determined by omp_get_num_procs.

determine_num_threads

dipy.denoise.denspeed.determine_num_threads()

Determine the effective number of threads to be used for OpenMP calls

  • For num_threads = None, - if the OMP_NUM_THREADS environment variable is set, return that value - otherwise, return the maximum number of cores retrieved by openmp.opm_get_num_procs().

  • For num_threads > 0, return this value.

  • For num_threads < 0, return the maximal number of threads minus |num_threads + 1|. In particular num_threads = -1 will use as many threads as there are available cores on the machine.

  • For num_threads = 0 a ValueError is raised.

Parameters
num_threadsint or None

Desired number of threads to be used.

nlmeans_3d

dipy.denoise.denspeed.nlmeans_3d()

Non-local means for denoising 3D images

Parameters
arr3D ndarray

The array to be denoised

mask3D ndarray
sigmafloat or 3D array

standard deviation of the noise estimated from the data

patch_radiusint

patch size is 2 x patch_radius + 1. Default is 1.

block_radiusint

block size is 2 x block_radius + 1. Default is 5.

ricianboolean

If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.

num_threadsint, optional

Number of threads to be used for OpenMP parallelization. If None (default) the value of OMP_NUM_THREADS environment variable is used if it is set, otherwise all available threads are used. If < 0 the maximal number of threads minus |num_threads + 1| is used (enter -1 to use as many threads as possible). 0 raises an error.

Returns
denoised_arrndarray

the denoised arr which has the same shape as arr.

remove_padding

dipy.denoise.denspeed.remove_padding()

EnhancementKernel

class dipy.denoise.enhancement_kernel.EnhancementKernel

Bases: object

Methods

evaluate_kernel

Evaluate the kernel at position x relative to position y, with orientation r relative to orientation v.

get_lookup_table

Return the computed look-up table.

get_orientations

Return the orientations.

get_sphere

Get the sphere corresponding with the orientations

__init__()

Compute a look-up table for the contextual enhancement kernel

Parameters
D33float

Spatial diffusion

D44float

Angular diffusion

tfloat

Diffusion time

force_recomputeboolean

Always compute the look-up table even if it is available in cache. Default is False.

orientationsinteger or Sphere object

Specify the number of orientations to be used with electrostatic repulsion, or provide a Sphere object. The default sphere is ‘repulsion100’.

verboseboolean

Enable verbose mode.

References

[Meesters2016_ISMRM] S. Meesters, G. Sanguinetti, E. Garyfallidis,

J. Portegies, R. Duits. (2016) Fast implementations of contextual PDE’s for HARDI data processing in DIPY. ISMRM 2016 conference.

[DuitsAndFranken_IJCV] R. Duits and E. Franken (2011) Left-invariant diffusions

on the space of positions and orientations and their application to crossing-preserving smoothing of HARDI images. International Journal of Computer Vision, 92:231-264.

[Portegies2015] J. Portegies, G. Sanguinetti, S. Meesters, and R. Duits.

(2015) New Approximation of a Scale Space Kernel on SE(3) and Applications in Neuroimaging. Fifth International Conference on Scale Space and Variational Methods in Computer Vision

[Portegies2015b] J. Portegies, R. Fick, G. Sanguinetti, S. Meesters,

G. Girard, and R. Duits. (2015) Improving Fiber Alignment in HARDI by Combining Contextual PDE flow with Constrained Spherical Deconvolution. PLoS One.

evaluate_kernel()

Evaluate the kernel at position x relative to position y, with orientation r relative to orientation v.

Parameters
x1D ndarray

Position x

y1D ndarray

Position y

r1D ndarray

Orientation r

v1D ndarray

Orientation v

Returns
kernel_valuedouble
get_lookup_table()

Return the computed look-up table.

get_orientations()

Return the orientations.

get_sphere()

Get the sphere corresponding with the orientations

HemiSphere

class dipy.denoise.enhancement_kernel.HemiSphere(x=None, y=None, z=None, theta=None, phi=None, xyz=None, faces=None, edges=None, tol=1e-05)

Bases: dipy.core.sphere.Sphere

Points on the unit sphere.

A HemiSphere is similar to a Sphere but it takes antipodal symmetry into account. Antipodal symmetry means that point v on a HemiSphere is the same as the point -v. Duplicate points are discarded when constructing a HemiSphere (including antipodal duplicates). edges and faces are remapped to the remaining points as closely as possible.

The HemiSphere can be constructed using one of three conventions:

HemiSphere(x, y, z)
HemiSphere(xyz=xyz)
HemiSphere(theta=theta, phi=phi)
Parameters
x, y, z1-D array_like

Vertices as x-y-z coordinates.

theta, phi1-D array_like

Vertices as spherical coordinates. Theta and phi are the inclination and azimuth angles respectively.

xyz(N, 3) ndarray

Vertices as x-y-z coordinates.

faces(N, 3) ndarray

Indices into vertices that form triangular faces. If unspecified, the faces are computed using a Delaunay triangulation.

edges(N, 2) ndarray

Edges between vertices. If unspecified, the edges are derived from the faces.

tolfloat

Angle in degrees. Vertices that are less than tol degrees apart are treated as duplicates.

See also

Sphere
Attributes
x
y
z

Methods

find_closest(xyz)

Find the index of the vertex in the Sphere closest to the input vector, taking into account antipodal symmetry

from_sphere(sphere[, tol])

Create instance from a Sphere

mirror()

Create a full Sphere from a HemiSphere

subdivide([n])

Create a more subdivided HemiSphere

edges

faces

vertices

__init__(x=None, y=None, z=None, theta=None, phi=None, xyz=None, faces=None, edges=None, tol=1e-05)

Create a HemiSphere from points

faces()
find_closest(xyz)

Find the index of the vertex in the Sphere closest to the input vector, taking into account antipodal symmetry

Parameters
xyzarray-like, 3 elements

A unit vector

Returns
idxint

The index into the Sphere.vertices array that gives the closest vertex (in angle).

classmethod from_sphere(sphere, tol=1e-05)

Create instance from a Sphere

mirror()

Create a full Sphere from a HemiSphere

subdivide(n=1)

Create a more subdivided HemiSphere

See Sphere.subdivide for full documentation.

Sphere

class dipy.denoise.enhancement_kernel.Sphere(x=None, y=None, z=None, theta=None, phi=None, xyz=None, faces=None, edges=None)

Bases: object

Points on the unit sphere.

The sphere can be constructed using one of three conventions:

Sphere(x, y, z)
Sphere(xyz=xyz)
Sphere(theta=theta, phi=phi)
Parameters
x, y, z1-D array_like

Vertices as x-y-z coordinates.

theta, phi1-D array_like

Vertices as spherical coordinates. Theta and phi are the inclination and azimuth angles respectively.

xyz(N, 3) ndarray

Vertices as x-y-z coordinates.

faces(N, 3) ndarray

Indices into vertices that form triangular faces. If unspecified, the faces are computed using a Delaunay triangulation.

edges(N, 2) ndarray

Edges between vertices. If unspecified, the edges are derived from the faces.

Attributes
x
y
z

Methods

find_closest(xyz)

Find the index of the vertex in the Sphere closest to the input vector

subdivide([n])

Subdivides each face of the sphere into four new faces.

edges

faces

vertices

__init__(x=None, y=None, z=None, theta=None, phi=None, xyz=None, faces=None, edges=None)

Initialize self. See help(type(self)) for accurate signature.

edges()
faces()
find_closest(xyz)

Find the index of the vertex in the Sphere closest to the input vector

Parameters
xyzarray-like, 3 elements

A unit vector

Returns
idxint

The index into the Sphere.vertices array that gives the closest vertex (in angle).

subdivide(n=1)

Subdivides each face of the sphere into four new faces.

New vertices are created at a, b, and c. Then each face [x, y, z] is divided into faces [x, a, c], [y, a, b], [z, b, c], and [a, b, c].

      y
      /\
     /  \
   a/____\b
   /\    /\
  /  \  /  \
 /____\/____\
x      c     z
Parameters
nint, optional

The number of subdivisions to preform.

Returns
new_sphereSphere

The subdivided sphere.

vertices()
property x
property y
property z

ceil

dipy.denoise.enhancement_kernel.ceil(x, /)

Return the ceiling of x as an Integral.

This is the smallest integer >= x.

disperse_charges

dipy.denoise.enhancement_kernel.disperse_charges(hemi, iters, const=0.2)

Models electrostatic repulsion on the unit sphere

Places charges on a sphere and simulates the repulsive forces felt by each one. Allows the charges to move for some number of iterations and returns their final location as well as the total potential of the system at each step.

Parameters
hemiHemiSphere

Points on a unit sphere.

itersint

Number of iterations to run.

constfloat

Using a smaller const could provide a more accurate result, but will need more iterations to converge.

Returns
hemiHemiSphere

Distributed points on a unit sphere.

potentialndarray

The electrostatic potential at each iteration. This can be useful to check if the repulsion converged to a minimum.

Notes

This function is meant to be used with diffusion imaging so antipodal symmetry is assumed. Therefor each charge must not only be unique, but if there is a charge at +x, there cannot be a charge at -x. These are treated as the same location and because the distance between the two charges will be zero, the result will be unstable.

get_sphere

dipy.denoise.enhancement_kernel.get_sphere(name='symmetric362')

provide triangulated spheres

Parameters
namestr

which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’

Returns
spherea dipy.core.sphere.Sphere class instance

Examples

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') 
Traceback (most recent call last):
    ...
DataError: No sphere called "not a sphere name"

gettempdir

dipy.denoise.enhancement_kernel.gettempdir()

Accessor for tempfile.tempdir.

Version

dipy.denoise.gibbs.Version

alias of numpy.lib._version.NumpyVersion

partial

class dipy.denoise.gibbs.partial

Bases: object

partial(func, *args, **keywords) - new function with partial application of the given arguments and keywords.

Attributes
args

tuple of arguments to future partial calls

func

function object to use in future partial calls

keywords

dictionary of keyword arguments to future partial calls

Methods

__call__(*args, **kwargs)

Call self as a function.

__init__(*args, **kwargs)

Initialize self. See help(type(self)) for accurate signature.

args

tuple of arguments to future partial calls

func

function object to use in future partial calls

keywords

dictionary of keyword arguments to future partial calls

Pool

dipy.denoise.gibbs.Pool(processes=None, initializer=None, initargs=(), maxtasksperchild=None)

Returns a process pool object

deprecated_params

dipy.denoise.gibbs.deprecated_params(old_name, new_name=None, since='', until='', version_comparator=<function cmp_pkg_version>, arg_in_kwargs=False, warn_class=<class 'dipy.utils.deprecator.ArgsDeprecationWarning'>, error_class=<class 'dipy.utils.deprecator.ExpiredDeprecationError'>, alternative='')

Deprecate a renamed or removed function argument.

The decorator assumes that the argument with the old_name was removed from the function signature and the new_name replaced it at the same position in the signature. If the old_name argument is given when calling the decorated function the decorator will catch it and issue a deprecation warning and pass it on as new_name argument.

Parameters
old_namestr or list/tuple thereof

The old name of the argument.

new_namestr or list/tuple thereof or None, optional

The new name of the argument. Set this to None to remove the argument old_name instead of renaming it.

sincestr or number or list/tuple thereof, optional

The release at which the old argument became deprecated.

untilstr or number or list/tuple thereof, optional

Last released version at which this function will still raise a deprecation warning. Versions higher than this will raise an error.

version_comparatorcallable

Callable accepting string as argument, and return 1 if string represents a higher version than encoded in the version_comparator, 0 if the version is equal, and -1 if the version is lower. For example, the version_comparator may compare the input version string to the current package version string.

arg_in_kwargsbool or list/tuple thereof, optional

If the argument is not a named argument (for example it was meant to be consumed by **kwargs) set this to True. Otherwise the decorator will throw an Exception if the new_name cannot be found in the signature of the decorated function. Default is False.

warn_classwarning, optional

Warning to be issued.

error_classException, optional

Error to be issued

alternativestr, optional

An alternative function or class name that the user may use in place of the deprecated object if new_name is None. The deprecation warning will tell the user about this alternative if provided.

Raises
TypeError

If the new argument name cannot be found in the function signature and arg_in_kwargs was False or if it is used to deprecate the name of the *args-, **kwargs-like arguments. At runtime such an Error is raised if both the new_name and old_name were specified when calling the function and “relax=False”.

Notes

This function is based on the Astropy (major modification). https://github.com/astropy/astropy. See COPYING file distributed along with the astropy package for the copyright and license terms.

Examples

The deprecation warnings are not shown in the following examples. To deprecate a positional or keyword argument:: >>> from dipy.utils.deprecator import deprecated_params >>> @deprecated_params(‘sig’, ‘sigma’, ‘0.3’) … def test(sigma): … return sigma >>> test(2) 2 >>> test(sigma=2) 2 >>> test(sig=2) # doctest: +SKIP 2

It is also possible to replace multiple arguments. The old_name, new_name and since have to be tuple or list and contain the same number of entries:: >>> @deprecated_params([‘a’, ‘b’], [‘alpha’, ‘beta’], … [‘0.2’, 0.4]) … def test(alpha, beta): … return alpha, beta >>> test(a=2, b=3) # doctest: +SKIP (2, 3)

determine_num_processes

dipy.denoise.gibbs.determine_num_processes(num_processes)

Determine the effective number of processes for parallelization.

  • For num_processes = None` return the maximum number of cores retrieved

by cpu_count().

  • For num_processes > 0, return this value.

  • For num_processes < 0, return the maximal number of cores minus

|num_processes + 1|. In particular num_processes = -1 will use as many cores as possible.

  • For num_processes = 0 a ValueError is raised.

Parameters
num_processesint or None

Desired number of processes to be used.

gibbs_removal

dipy.denoise.gibbs.gibbs_removal(vol, slice_axis=2, n_points=3, inplace=True, num_processes=1)

Suppresses Gibbs ringing artefacts of images volumes.

Parameters
volndarray ([X, Y]), ([X, Y, Z]) or ([X, Y, Z, g])

Matrix containing one volume (3D) or multiple (4D) volumes of images.

slice_axisint (0, 1, or 2)

Data axis corresponding to the number of acquired slices. Default is set to the third axis.

n_pointsint, optional

Number of neighbour points to access local TV (see note). Default is set to 3.

inplacebool, optional

If True, the input data is replaced with results. Otherwise, returns a new array. Default is set to True.

num_processesint or None, optional

Split the calculation to a pool of children processes. This only applies to 3D or 4D data arrays. Default is 1. If < 0 the maximal number of cores minus |num_processes + 1| is used (enter -1 to use as many cores as possible). 0 raises an error.

Returns
volndarray ([X, Y]), ([X, Y, Z]) or ([X, Y, Z, g])

Matrix containing one volume (3D) or multiple (4D) volumes of corrected images.

Notes

For 4D matrix last element should always correspond to the number of diffusion gradient directions.

References

Please cite the following articles .. [Rfce522a872fd-1] Neto Henriques, R., 2018. Advanced Methods for Diffusion MRI Data

Analysis and their Application to the Healthy Ageing Brain (Doctoral thesis). https://doi.org/10.17863/CAM.29356

2

Kellner E, Dhital B, Kiselev VG, Reisert M. Gibbs-ringing artifact removal based on local subvoxel-shifts. Magn Reson Med. 2016 doi: 10.1002/mrm.26054.

eigh

dipy.denoise.localpca.eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=True, eigvals=None, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None)

Solve a standard or generalized eigenvalue problem for a complex Hermitian or real symmetric matrix.

Find eigenvalues array w and optionally eigenvectors array v of array a, where b is positive definite such that for every eigenvalue λ (i-th entry of w) and its eigenvector vi (i-th column of v) satisfies:

              a @ vi = λ * b @ vi
vi.conj().T @ a @ vi = λ
vi.conj().T @ b @ vi = 1

In the standard problem, b is assumed to be the identity matrix.

Parameters
a(M, M) array_like

A complex Hermitian or real symmetric matrix whose eigenvalues and eigenvectors will be computed.

b(M, M) array_like, optional

A complex Hermitian or real symmetric definite positive matrix in. If omitted, identity matrix is assumed.

lowerbool, optional

Whether the pertinent array data is taken from the lower or upper triangle of a and, if applicable, b. (Default: lower)

eigvals_onlybool, optional

Whether to calculate only eigenvalues and no eigenvectors. (Default: both are calculated)

subset_by_indexiterable, optional

If provided, this two-element iterable defines the start and the end indices of the desired eigenvalues (ascending order and 0-indexed). To return only the second smallest to fifth smallest eigenvalues, [1, 4] is used. [n-3, n-1] returns the largest three. Only available with “evr”, “evx”, and “gvx” drivers. The entries are directly converted to integers via int().

subset_by_valueiterable, optional

If provided, this two-element iterable defines the half-open interval (a, b] that, if any, only the eigenvalues between these values are returned. Only available with “evr”, “evx”, and “gvx” drivers. Use np.inf for the unconstrained ends.

driver: str, optional

Defines which LAPACK driver should be used. Valid options are “ev”, “evd”, “evr”, “evx” for standard problems and “gv”, “gvd”, “gvx” for generalized (where b is not None) problems. See the Notes section.

typeint, optional

For the generalized problems, this keyword specifies the problem type to be solved for w and v (only takes 1, 2, 3 as possible inputs):

1 =>     a @ v = w @ b @ v
2 => a @ b @ v = w @ v
3 => b @ a @ v = w @ v

This keyword is ignored for standard problems.

overwrite_abool, optional

Whether to overwrite data in a (may improve performance). Default is False.

overwrite_bbool, optional

Whether to overwrite data in b (may improve performance). Default is False.

check_finitebool, optional

Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.

turbobool, optional

Deprecated since v1.5.0, use ``driver=gvd`` keyword instead. Use divide and conquer algorithm (faster but expensive in memory, only for generalized eigenvalue problem and if full set of eigenvalues are requested.). Has no significant effect if eigenvectors are not requested.

eigvalstuple (lo, hi), optional

Deprecated since v1.5.0, use ``subset_by_index`` keyword instead. Indexes of the smallest and largest (in ascending order) eigenvalues and corresponding eigenvectors to be returned: 0 <= lo <= hi <= M-1. If omitted, all eigenvalues and eigenvectors are returned.

Returns
w(N,) ndarray

The N (1<=N<=M) selected eigenvalues, in ascending order, each repeated according to its multiplicity.

v(M, N) ndarray

(if eigvals_only == False)

Raises
LinAlgError

If eigenvalue computation does not converge, an error occurred, or b matrix is not definite positive. Note that if input matrices are not symmetric or Hermitian, no error will be reported but results will be wrong.

See also

eigvalsh

eigenvalues of symmetric or Hermitian arrays

eig

eigenvalues and right eigenvectors for non-symmetric arrays

eigh_tridiagonal

eigenvalues and right eiegenvectors for symmetric/Hermitian tridiagonal matrices

Notes

This function does not check the input array for being hermitian/symmetric in order to allow for representing arrays with only their upper/lower triangular parts. Also, note that even though not taken into account, finiteness check applies to the whole array and unaffected by “lower” keyword.

This function uses LAPACK drivers for computations in all possible keyword combinations, prefixed with sy if arrays are real and he if complex, e.g., a float array with “evr” driver is solved via “syevr”, complex arrays with “gvx” driver problem is solved via “hegvx” etc.

As a brief summary, the slowest and the most robust driver is the classical <sy/he>ev which uses symmetric QR. <sy/he>evr is seen as the optimal choice for the most general cases. However, there are certain occassions that <sy/he>evd computes faster at the expense of more memory usage. <sy/he>evx, while still being faster than <sy/he>ev, often performs worse than the rest except when very few eigenvalues are requested for large arrays though there is still no performance guarantee.

For the generalized problem, normalization with respoect to the given type argument:

type 1 and 3 :      v.conj().T @ a @ v = w
type 2       : inv(v).conj().T @ a @ inv(v) = w

type 1 or 2  :      v.conj().T @ b @ v  = I
type 3       : v.conj().T @ inv(b) @ v  = I

Examples

>>> from scipy.linalg import eigh
>>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
>>> w, v = eigh(A)
>>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
True

Request only the eigenvalues

>>> w = eigh(A, eigvals_only=True)

Request eigenvalues that are less than 10.

>>> A = np.array([[34, -4, -10, -7, 2],
...               [-4, 7, 2, 12, 0],
...               [-10, 2, 44, 2, -19],
...               [-7, 12, 2, 79, -34],
...               [2, 0, -19, -34, 29]])
>>> eigh(A, eigvals_only=True, subset_by_value=[-np.inf, 10])
array([6.69199443e-07, 9.11938152e+00])

Request the largest second eigenvalue and its eigenvector

>>> w, v = eigh(A, subset_by_index=[1, 1])
>>> w
array([9.11938152])
>>> v.shape  # only a single column is returned
(5, 1)

genpca

dipy.denoise.localpca.genpca(arr, sigma=None, mask=None, patch_radius=2, pca_method='eig', tau_factor=None, return_sigma=False, out_dtype=None)

General function to perform PCA-based denoising of diffusion datasets.

Parameters
arr4D array

Array of data to be denoised. The dimensions are (X, Y, Z, N), where N are the diffusion gradient directions.

sigmafloat or 3D array (optional)

Standard deviation of the noise estimated from the data. If no sigma is given, this will be estimated based on random matrix theory [1],[R4bed7265c934-2]_

mask3D boolean array (optional)

A mask with voxels that are true inside the brain and false outside of it. The function denoises within the true part and returns zeros outside of those voxels.

patch_radiusint or 1D array (optional)

The radius of the local patch to be taken around each voxel (in voxels). Default: 2 (denoise in blocks of 5x5x5 voxels).

pca_method‘eig’ or ‘svd’ (optional)

Use either eigenvalue decomposition (eig) or singular value decomposition (svd) for principal component analysis. The default method is ‘eig’ which is faster. However, occasionally ‘svd’ might be more accurate.

tau_factorfloat (optional)

Thresholding of PCA eigenvalues is done by nulling out eigenvalues that are smaller than:

\[\tau = (\tau_{factor} \sigma)^2\]

tau_{factor} can be set to a predefined values (e.g. tau_{factor} = 2.3 [3]), or automatically calculated using random matrix theory (in case that tau_{factor} is set to None). Default: None.

return_sigmabool (optional)

If true, the Standard deviation of the noise will be returned. Default: False.

out_dtypestr or dtype (optional)

The dtype for the output array. Default: output has the same dtype as the input.

Returns
denoised_arr4D array

This is the denoised array of the same size as that of the input data, clipped to non-negative values

References

1

Veraart J, Novikov DS, Christiaens D, Ades-aron B, Sijbers, Fieremans E, 2016. Denoising of Diffusion MRI using random matrix theory. Neuroimage 142:394-406. doi: 10.1016/j.neuroimage.2016.08.016

2

Veraart J, Fieremans E, Novikov DS. 2016. Diffusion MRI noise mapping using random matrix theory. Magnetic Resonance in Medicine. doi: 10.1002/mrm.26059.

3

Manjon JV, Coupe P, Concha L, Buades A, Collins DL (2013) Diffusion Weighted Image Denoising Using Overcomplete Local PCA. PLoS ONE 8(9): e73021. https://doi.org/10.1371/journal.pone.0073021

localpca

dipy.denoise.localpca.localpca(arr, sigma, mask=None, patch_radius=2, pca_method='eig', tau_factor=2.3, out_dtype=None)

Performs local PCA denoising according to Manjon et al. [1].

Parameters
arr4D array

Array of data to be denoised. The dimensions are (X, Y, Z, N), where N are the diffusion gradient directions.

sigmafloat or 3D array

Standard deviation of the noise estimated from the data.

mask3D boolean array (optional)

A mask with voxels that are true inside the brain and false outside of it. The function denoises within the true part and returns zeros outside of those voxels.

patch_radiusint or 1D array (optional)

The radius of the local patch to be taken around each voxel (in voxels). Default: 2 (denoise in blocks of 5x5x5 voxels).

pca_method‘eig’ or ‘svd’ (optional)

Use either eigenvalue decomposition (eig) or singular value decomposition (svd) for principal component analysis. The default method is ‘eig’ which is faster. However, occasionally ‘svd’ might be more accurate.

tau_factorfloat (optional)

Thresholding of PCA eigenvalues is done by nulling out eigenvalues that are smaller than:

\[\tau = (\tau_{factor} \sigma)^2\]

tau_{factor} can be change to adjust the relationship between the noise standard deviation and the threshold tau. If tau_{factor} is set to None, it will be automatically calculated using the Marcenko-Pastur distribution [2]. Default: 2.3 (according to [1])

out_dtypestr or dtype (optional)

The dtype for the output array. Default: output has the same dtype as the input.

Returns
denoised_arr4D array

This is the denoised array of the same size as that of the input data, clipped to non-negative values

References

1(1,2)

Manjon JV, Coupe P, Concha L, Buades A, Collins DL (2013) Diffusion Weighted Image Denoising Using Overcomplete Local PCA. PLoS ONE 8(9): e73021. https://doi.org/10.1371/journal.pone.0073021

2

Veraart J, Novikov DS, Christiaens D, Ades-aron B, Sijbers, Fieremans E, 2016. Denoising of Diffusion MRI using random matrix theory. Neuroimage 142:394-406. doi: 10.1016/j.neuroimage.2016.08.016

mppca

dipy.denoise.localpca.mppca(arr, mask=None, patch_radius=2, pca_method='eig', return_sigma=False, out_dtype=None)

Performs PCA-based denoising using the Marcenko-Pastur distribution [1].

Parameters
arr4D array

Array of data to be denoised. The dimensions are (X, Y, Z, N), where N are the diffusion gradient directions.

mask3D boolean array (optional)

A mask with voxels that are true inside the brain and false outside of it. The function denoises within the true part and returns zeros outside of those voxels.

patch_radiusint or 1D array (optional)

The radius of the local patch to be taken around each voxel (in voxels). Default: 2 (denoise in blocks of 5x5x5 voxels).

pca_method‘eig’ or ‘svd’ (optional)

Use either eigenvalue decomposition (eig) or singular value decomposition (svd) for principal component analysis. The default method is ‘eig’ which is faster. However, occasionally ‘svd’ might be more accurate.

return_sigmabool (optional)

If true, a noise standard deviation estimate based on the Marcenko-Pastur distribution is returned [2]. Default: False.

out_dtypestr or dtype (optional)

The dtype for the output array. Default: output has the same dtype as the input.

Returns
denoised_arr4D array

This is the denoised array of the same size as that of the input data, clipped to non-negative values

sigma3D array (when return_sigma=True)

Estimate of the spatial varying standard deviation of the noise

References

1(1,2)

Veraart J, Novikov DS, Christiaens D, Ades-aron B, Sijbers, Fieremans E, 2016. Denoising of Diffusion MRI using random matrix theory. Neuroimage 142:394-406. doi: 10.1016/j.neuroimage.2016.08.016

2

Veraart J, Fieremans E, Novikov DS. 2016. Diffusion MRI noise mapping using random matrix theory. Magnetic Resonance in Medicine. doi: 10.1002/mrm.26059.

nlmeans

dipy.denoise.nlmeans.nlmeans(arr, sigma, mask=None, patch_radius=1, block_radius=5, rician=True, num_threads=None)

Non-local means for denoising 3D and 4D images

Parameters
arr3D or 4D ndarray

The array to be denoised

mask3D ndarray
sigmafloat or 3D array

standard deviation of the noise estimated from the data

patch_radiusint

patch size is 2 x patch_radius + 1. Default is 1.

block_radiusint

block size is 2 x block_radius + 1. Default is 5.

ricianboolean

If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.

num_threadsint, optional

Number of threads to be used for OpenMP parallelization. If None (default) the value of OMP_NUM_THREADS environment variable is used if it is set, otherwise all available threads are used. If < 0 the maximal number of threads minus |num_threads + 1| is used (enter -1 to use as many threads as possible). 0 raises an error.

Returns
denoised_arrndarray

the denoised arr which has the same shape as arr.

References

Descoteaux08

Descoteaux, Maxime and Wiest-Daesslé, Nicolas and Prima, Sylvain and Barillot, Christian and Deriche, Rachid Impact of Rician Adapted Non-Local Means Filtering on HARDI, MICCAI 2008

nlmeans_3d

dipy.denoise.nlmeans.nlmeans_3d()

Non-local means for denoising 3D images

Parameters
arr3D ndarray

The array to be denoised

mask3D ndarray
sigmafloat or 3D array

standard deviation of the noise estimated from the data

patch_radiusint

patch size is 2 x patch_radius + 1. Default is 1.

block_radiusint

block size is 2 x block_radius + 1. Default is 5.

ricianboolean

If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.

num_threadsint, optional

Number of threads to be used for OpenMP parallelization. If None (default) the value of OMP_NUM_THREADS environment variable is used if it is set, otherwise all available threads are used. If < 0 the maximal number of threads minus |num_threads + 1| is used (enter -1 to use as many threads as possible). 0 raises an error.

Returns
denoised_arrndarray

the denoised arr which has the same shape as arr.

firdn

dipy.denoise.nlmeans_block.firdn()

Applies the filter given by the convolution kernel ‘h’ columnwise to ‘image’, then subsamples by 2. This is a special case of the matlab’s ‘upfirdn’ function, ported to python. Returns the filtered image.

Parameters
image: 2D array of doubles

the input image to be filtered

h: double array

the convolution kernel

nlmeans_block

dipy.denoise.nlmeans_block.nlmeans_block()

Non-Local Means Denoising Using Blockwise Averaging

Parameters
image3D array of doubles

the input image, corrupted with rician noise

mask3D array of doubles

the input mask

patch_radiusint

similar patches in the non-local means are searched for locally, inside a cube of side 2*v+1 centered at each voxel of interest.

block_radiusint

the size of the block to be used (2*f+1)x(2*f+1)x(2*f+1) in the blockwise non-local means implementation (the Coupe’s proposal).

hdouble

the estimated amount of rician noise in the input image: in P. Coupe et al. the rician noise was simulated as sqrt((f+x)^2 + (y)^2) where f is the pixel value and x and y are independent realizations of a random variable with Normal distribution, with mean=0 and standard deviation=h

ricianboolean

If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.

Returns
fima: 3D double array

the denoised output which has the same shape as input image.

References

[1] P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, C. Barillot,

“An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic Resonance Images” IEEE Transactions on Medical Imaging, 27(4):425-441, 2008

[2] Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins.

“Multiresolution Non-Local Means Filter for 3D MR Image Denoising” IET Image Processing, Institution of Engineering and Technology, 2011

upfir

dipy.denoise.nlmeans_block.upfir()

Upsamples the columns of the input image by 2, then applies the convolution kernel ‘h’ (again, columnwise). This is a special case of the matlab’s ‘upfirdn’ function, ported to python. Returns the filtered image.

Parameters
image: 2D array of doubles

the input image to be filtered

h: double array

the convolution kernel

convolve

dipy.denoise.noise_estimate.convolve(input, weights, output=None, mode='reflect', cval=0.0, origin=0)

Multidimensional convolution.

The array is convolved with the given kernel.

Parameters
inputarray_like

The input array.

weightsarray_like

Array of weights, same number of dimensions as input

outputarray or dtype, optional

The array in which to place the output, or the dtype of the returned array. By default an array of the same dtype as input will be created.

mode{‘reflect’, ‘constant’, ‘nearest’, ‘mirror’, ‘wrap’}, optional

The mode parameter determines how the input array is extended beyond its boundaries. Default is ‘reflect’. Behavior for each valid value is as follows:

‘reflect’ (d c b a | a b c d | d c b a)

The input is extended by reflecting about the edge of the last pixel.

‘constant’ (k k k k | a b c d | k k k k)

The input is extended by filling all values beyond the edge with the same constant value, defined by the cval parameter.

‘nearest’ (a a a a | a b c d | d d d d)

The input is extended by replicating the last pixel.

‘mirror’ (d c b | a b c d | c b a)

The input is extended by reflecting about the center of the last pixel.

‘wrap’ (a b c d | a b c d | a b c d)

The input is extended by wrapping around to the opposite edge.

cvalscalar, optional

Value to fill past edges of input if mode is ‘constant’. Default is 0.0

originint or sequence, optional

Controls the placement of the filter on the input array’s pixels. A value of 0 (the default) centers the filter over the pixel, with positive values shifting the filter to the left, and negative ones to the right. By passing a sequence of origins with length equal to the number of dimensions of the input array, different shifts can be specified along each axis.

Returns
resultndarray

The result of convolution of input with weights.

See also

correlate

Correlate an image with a kernel.

Notes

Each value in result is \(C_i = \sum_j{I_{i+k-j} W_j}\), where W is the weights kernel, j is the N-D spatial index over \(W\), I is the input and k is the coordinate of the center of W, specified by origin in the input parameters.

Examples

Perhaps the simplest case to understand is mode='constant', cval=0.0, because in this case borders (i.e., where the weights kernel, centered on any one value, extends beyond an edge of input) are treated as zeros.

>>> a = np.array([[1, 2, 0, 0],
...               [5, 3, 0, 4],
...               [0, 0, 0, 7],
...               [9, 3, 0, 0]])
>>> k = np.array([[1,1,1],[1,1,0],[1,0,0]])
>>> from scipy import ndimage
>>> ndimage.convolve(a, k, mode='constant', cval=0.0)
array([[11, 10,  7,  4],
       [10,  3, 11, 11],
       [15, 12, 14,  7],
       [12,  3,  7,  0]])

Setting cval=1.0 is equivalent to padding the outer edge of input with 1.0’s (and then extracting only the original region of the result).

>>> ndimage.convolve(a, k, mode='constant', cval=1.0)
array([[13, 11,  8,  7],
       [11,  3, 11, 14],
       [16, 12, 14, 10],
       [15,  6, 10,  5]])

With mode='reflect' (the default), outer values are reflected at the edge of input to fill in missing values.

>>> b = np.array([[2, 0, 0],
...               [1, 0, 0],
...               [0, 0, 0]])
>>> k = np.array([[0,1,0], [0,1,0], [0,1,0]])
>>> ndimage.convolve(b, k, mode='reflect')
array([[5, 0, 0],
       [3, 0, 0],
       [1, 0, 0]])

This includes diagonally at the corners.

>>> k = np.array([[1,0,0],[0,1,0],[0,0,1]])
>>> ndimage.convolve(b, k)
array([[4, 2, 0],
       [3, 2, 0],
       [1, 1, 0]])

With mode='nearest', the single nearest value in to an edge in input is repeated as many times as needed to match the overlapping weights.

>>> c = np.array([[2, 0, 1],
...               [1, 0, 0],
...               [0, 0, 0]])
>>> k = np.array([[0, 1, 0],
...               [0, 1, 0],
...               [0, 1, 0],
...               [0, 1, 0],
...               [0, 1, 0]])
>>> ndimage.convolve(c, k, mode='nearest')
array([[7, 0, 3],
       [5, 0, 2],
       [3, 0, 1]])

estimate_sigma

dipy.denoise.noise_estimate.estimate_sigma(arr, disable_background_masking=False, N=0)

Standard deviation estimation from local patches

Parameters
arr3D or 4D ndarray

The array to be estimated

disable_background_maskingbool, default False

If True, uses all voxels for the estimation, otherwise, only non-zeros voxels are used. Useful if the background is masked by the scanner.

Nint, default 0

Number of coils of the receiver array. Use N = 1 in case of a SENSE reconstruction (Philips scanners) or the number of coils for a GRAPPA reconstruction (Siemens and GE). Use 0 to disable the correction factor, as for example if the noise is Gaussian distributed. See [1] for more information.

Returns
sigmandarray

standard deviation of the noise, one estimation per volume.

Notes

This function is the same as manually taking the standard deviation of the background and gives one value for the whole 3D array. It also includes the coil-dependent correction factor of Koay 2006 (see [1], equation 18) with theta = 0. Since this function was introduced in [2] for T1 imaging, it is expected to perform ok on diffusion MRI data, but might oversmooth some regions and leave others un-denoised for spatially varying noise profiles. Consider using piesno() to estimate sigma instead if visual inaccuracies are apparent in the denoised result.

References

1

Koay, C. G., & Basser, P. J. (2006). Analytically exact correction

scheme for signal extraction from noisy magnitude MR signals. Journal of Magnetic Resonance), 179(2), 317-22.

2

Coupe, P., Yger, P., Prima, S., Hellier, P., Kervrann, C., Barillot,

C., 2008. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images, IEEE Trans. Med. Imaging 27, 425-41.

piesno

dipy.denoise.noise_estimate.piesno(data, N, alpha=0.01, l=100, itermax=100, eps=1e-05, return_mask=False)

Probabilistic Identification and Estimation of Noise (PIESNO).

Parameters
datandarray

The magnitude signals to analyse. The last dimension must contain the same realisation of the volume, such as dMRI or fMRI data.

Nint

The number of phase array coils of the MRI scanner. If your scanner does a SENSE reconstruction, ALWAYS use N=1, as the noise profile is always Rician. If your scanner does a GRAPPA reconstruction, set N as the number of phase array coils.

alphafloat

Probabilistic estimation threshold for the gamma function.

lint

number of initial estimates for sigma to try.

itermaxint

Maximum number of iterations to execute if convergence is not reached.

epsfloat

Tolerance for the convergence criterion. Convergence is reached if two subsequent estimates are smaller than eps.

return_maskbool

If True, return a mask identifying all the pure noise voxel that were found.

Returns
sigmafloat

The estimated standard deviation of the gaussian noise.

maskndarray (optional)

A boolean mask indicating the voxels identified as pure noise.

Notes

This function assumes two things : 1. The data has a noisy, non-masked background and 2. The data is a repetition of the same measurements along the last axis, i.e. dMRI or fMRI data, not structural data like T1/T2.

This function processes the data slice by slice, as originally designed in the paper. Use it to get a slice by slice estimation of the noise, as in spinal cord imaging for example.

References

1

Koay CG, Ozarslan E and Pierpaoli C.

“Probabilistic Identification and Estimation of Noise (PIESNO): A self-consistent approach and its applications in MRI.” Journal of Magnetic Resonance 2009; 199: 94-103.

2

Koay CG, Ozarslan E and Basser PJ.

“A signal transformational framework for breaking the noise floor and its applications in MRI.” Journal of Magnetic Resonance 2009; 197: 108-119.

nlmeans_block

dipy.denoise.non_local_means.nlmeans_block()

Non-Local Means Denoising Using Blockwise Averaging

Parameters
image3D array of doubles

the input image, corrupted with rician noise

mask3D array of doubles

the input mask

patch_radiusint

similar patches in the non-local means are searched for locally, inside a cube of side 2*v+1 centered at each voxel of interest.

block_radiusint

the size of the block to be used (2*f+1)x(2*f+1)x(2*f+1) in the blockwise non-local means implementation (the Coupe’s proposal).

hdouble

the estimated amount of rician noise in the input image: in P. Coupe et al. the rician noise was simulated as sqrt((f+x)^2 + (y)^2) where f is the pixel value and x and y are independent realizations of a random variable with Normal distribution, with mean=0 and standard deviation=h

ricianboolean

If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.

Returns
fima: 3D double array

the denoised output which has the same shape as input image.

References

[1] P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, C. Barillot,

“An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic Resonance Images” IEEE Transactions on Medical Imaging, 27(4):425-441, 2008

[2] Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins.

“Multiresolution Non-Local Means Filter for 3D MR Image Denoising” IET Image Processing, Institution of Engineering and Technology, 2011

non_local_means

dipy.denoise.non_local_means.non_local_means(arr, sigma, mask=None, patch_radius=1, block_radius=5, rician=True)
Non-local means for denoising 3D and 4D images, using

blockwise averaging approach

Parameters
arr3D or 4D ndarray

The array to be denoised

mask3D ndarray
sigmafloat

standard deviation of the noise estimated from the data

patch_radiusint

patch size is 2 x patch_radius + 1. Default is 1.

block_radiusint

block size is 2 x block_radius + 1. Default is 5.

ricianboolean

If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.

Returns
denoised_arrndarray

the denoised arr which has the same shape as arr.

References

Coupe08

P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, C. Barillot, An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic Resonance Images, IEEE Transactions on Medical Imaging, 27(4):425-441, 2008

Coupe11

Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins. Adaptive Multiresolution Non-Local Means Filter for 3D MR Image Denoising IET Image Processing, Institution of Engineering and Technology, 2011

optional_package

dipy.denoise.patch2self.optional_package(name, trip_msg=None)

Return package-like thing and module setup for package name

Parameters
namestr

package name

trip_msgNone or str

message to give when someone tries to use the return package, but we could not import it, and have returned a TripWire object instead. Default message if None.

Returns
pkg_likemodule or TripWire instance

If we can import the package, return it. Otherwise return an object raising an error when accessed

have_pkgbool

True if import for package was successful, false otherwise

module_setupfunction

callable usually set as setup_module in calling namespace, to allow skipping tests.

Examples

Typical use would be something like this at the top of a module using an optional package:

>>> from dipy.utils.optpkg import optional_package
>>> pkg, have_pkg, setup_module = optional_package('not_a_package')

Of course in this case the package doesn’t exist, and so, in the module:

>>> have_pkg
False

and

>>> pkg.some_function() 
Traceback (most recent call last):
    ...
TripWireError: We need package not_a_package for these functions, but
``import not_a_package`` raised an ImportError

If the module does exist - we get the module

>>> pkg, _, _ = optional_package('os')
>>> hasattr(pkg, 'path')
True

Or a submodule if that’s what we asked for

>>> subpkg, _, _ = optional_package('os.path')
>>> hasattr(subpkg, 'dirname')
True

patch2self

dipy.denoise.patch2self.patch2self(data, bvals, patch_radius=[0, 0, 0], model='ols', b0_threshold=50, out_dtype=None, alpha=1.0, verbose=False, b0_denoising=True, clip_negative_vals=False, shift_intensity=True)

Patch2Self Denoiser

Parameters
datandarray

The 4D noisy DWI data to be denoised.

bvals1D array

Array of the bvals from the DWI acquisition

patch_radiusint or 1D array, optional

The radius of the local patch to be taken around each voxel (in voxels). Default: 0 (denoise in blocks of 1x1x1 voxels).

modelstring, or initialized linear model object.

This will determine the algorithm used to solve the set of linear equations underlying this model. If it is a string it needs to be one of the following: {‘ols’, ‘ridge’, ‘lasso’}. Otherwise, it can be an object that inherits from dipy.optimize.SKLearnLinearSolver or an object with a similar interface from Scikit-Learn: sklearn.linear_model.LinearRegression, sklearn.linear_model.Lasso or sklearn.linear_model.Ridge and other objects that inherit from sklearn.base.RegressorMixin. Default: ‘ols’.

b0_thresholdint, optional

Threshold for considering volumes as b0.

out_dtypestr or dtype, optional

The dtype for the output array. Default: output has the same dtype as the input.

alphafloat, optional

Regularization parameter only for ridge regression model. Default: 1.0

verbosebool, optional

Show progress of Patch2Self and time taken.

b0_denoisingbool, optional

Skips denoising b0 volumes if set to False. Default: True

clip_negative_valsbool, optional

Sets negative values after denoising to 0 using np.clip. Default: True

shift_intensitybool, optional

Shifts the distribution of intensities per volume to give non-negative values Default: False

Returns
denoised arrayndarray

This is the denoised array of the same size as that of the input data, clipped to non-negative values.

References

[Fadnavis20] S. Fadnavis, J. Batson, E. Garyfallidis, Patch2Self:

Denoising Diffusion MRI with Self-supervised Learning, Advances in Neural Information Processing Systems 33 (2020)

warn

dipy.denoise.patch2self.warn(/, message, category=None, stacklevel=1, source=None)

Issue a warning, or maybe ignore it or raise an exception.

pca_noise_estimate

dipy.denoise.pca_noise_estimate.pca_noise_estimate()

PCA based local noise estimation.

Parameters
data: 4D array

the input dMRI data.

gtab: gradient table object

gradient information for the data gives us the bvals and bvecs of diffusion data, which is needed here to select between the noise estimation methods.

patch_radiusint

The radius of the local patch to be taken around each voxel (in voxels). Default: 1 (estimate noise in blocks of 3x3x3 voxels).

correct_biasbool

Whether to correct for bias due to Rician noise. This is an implementation of equation 8 in [1].

smoothint

Radius of a Gaussian smoothing filter to apply to the noise estimate before returning. Default: 2.

Returns
sigma_corr: 3D array

The local noise standard deviation estimate.

References

1

Manjon JV, Coupe P, Concha L, Buades A, Collins DL “Diffusion Weighted Image Denoising Using Overcomplete Local PCA”. PLoS ONE 8(9): e73021. doi:10.1371/journal.pone.0073021.

EnhancementKernel

class dipy.denoise.shift_twist_convolution.EnhancementKernel

Bases: object

Methods

evaluate_kernel

Evaluate the kernel at position x relative to position y, with orientation r relative to orientation v.

get_lookup_table

Return the computed look-up table.

get_orientations

Return the orientations.

get_sphere

Get the sphere corresponding with the orientations

__init__()

Compute a look-up table for the contextual enhancement kernel

Parameters
D33float

Spatial diffusion

D44float

Angular diffusion

tfloat

Diffusion time

force_recomputeboolean

Always compute the look-up table even if it is available in cache. Default is False.

orientationsinteger or Sphere object

Specify the number of orientations to be used with electrostatic repulsion, or provide a Sphere object. The default sphere is ‘repulsion100’.

verboseboolean

Enable verbose mode.

References

[Meesters2016_ISMRM] S. Meesters, G. Sanguinetti, E. Garyfallidis,

J. Portegies, R. Duits. (2016) Fast implementations of contextual PDE’s for HARDI data processing in DIPY. ISMRM 2016 conference.

[DuitsAndFranken_IJCV] R. Duits and E. Franken (2011) Left-invariant diffusions

on the space of positions and orientations and their application to crossing-preserving smoothing of HARDI images. International Journal of Computer Vision, 92:231-264.

[Portegies2015] J. Portegies, G. Sanguinetti, S. Meesters, and R. Duits.

(2015) New Approximation of a Scale Space Kernel on SE(3) and Applications in Neuroimaging. Fifth International Conference on Scale Space and Variational Methods in Computer Vision

[Portegies2015b] J. Portegies, R. Fick, G. Sanguinetti, S. Meesters,

G. Girard, and R. Duits. (2015) Improving Fiber Alignment in HARDI by Combining Contextual PDE flow with Constrained Spherical Deconvolution. PLoS One.

evaluate_kernel()

Evaluate the kernel at position x relative to position y, with orientation r relative to orientation v.

Parameters
x1D ndarray

Position x

y1D ndarray

Position y

r1D ndarray

Orientation r

v1D ndarray

Orientation v

Returns
kernel_valuedouble
get_lookup_table()

Return the computed look-up table.

get_orientations()

Return the orientations.

get_sphere()

Get the sphere corresponding with the orientations

convolve

dipy.denoise.shift_twist_convolution.convolve()

Perform the shift-twist convolution with the ODF data and the lookup-table of the kernel.

Parameters
odfsarray of double

The ODF data in spherical harmonics format

kernelarray of double

The 5D lookup table

sh_orderinteger

Maximal spherical harmonics order

test_modeboolean

Reduced convolution in one direction only for testing

num_threadsint, optional

Number of threads to be used for OpenMP parallelization. If None (default) the value of OMP_NUM_THREADS environment variable is used if it is set, otherwise all available threads are used. If < 0 the maximal number of threads minus |num_threads + 1| is used (enter -1 to use as many threads as possible). 0 raises an error.

normalizeboolean

Apply max-normalization to the output such that its value range matches the input ODF data.

Returns
outputarray of double

The ODF data after convolution enhancement in spherical harmonics format

References

[Meesters2016_ISMRM] S. Meesters, G. Sanguinetti, E. Garyfallidis,

J. Portegies, R. Duits. (2016) Fast implementations of contextual PDE’s for HARDI data processing in DIPY. ISMRM 2016 conference.

[DuitsAndFranken_IJCV] R. Duits and E. Franken (2011) Left-invariant diffusions

on the space of positions and orientations and their application to crossing-preserving smoothing of HARDI images. International Journal of Computer Vision, 92:231-264.

[Portegies2015] J. Portegies, G. Sanguinetti, S. Meesters, and R. Duits.

(2015) New Approximation of a Scale Space Kernel on SE(3) and Applications in Neuroimaging. Fifth International Conference on Scale Space and Variational Methods in Computer Vision

[Portegies2015b] J. Portegies, R. Fick, G. Sanguinetti, S. Meesters, G.Girard,

and R. Duits. (2015) Improving Fiber Alignment in HARDI by Combining Contextual PDE flow with Constrained Spherical Deconvolution. PLoS One.

convolve_sf

dipy.denoise.shift_twist_convolution.convolve_sf()

Perform the shift-twist convolution with the ODF data and the lookup-table of the kernel.

Parameters
odfsarray of double

The ODF data sampled on a sphere

kernelarray of double

The 5D lookup table

test_modeboolean

Reduced convolution in one direction only for testing

num_threadsint, optional

Number of threads to be used for OpenMP parallelization. If None (default) the value of OMP_NUM_THREADS environment variable is used if it is set, otherwise all available threads are used. If < 0 the maximal number of threads minus |num_threads + 1| is used (enter -1 to use as many threads as possible). 0 raises an error.

normalizeboolean

Apply max-normalization to the output such that its value range matches the input ODF data.

Returns
outputarray of double

The ODF data after convolution enhancement, sampled on a sphere

cpu_count

dipy.denoise.shift_twist_convolution.cpu_count()

Return number of cpus as determined by omp_get_num_procs.

determine_num_threads

dipy.denoise.shift_twist_convolution.determine_num_threads()

Determine the effective number of threads to be used for OpenMP calls

  • For num_threads = None, - if the OMP_NUM_THREADS environment variable is set, return that value - otherwise, return the maximum number of cores retrieved by openmp.opm_get_num_procs().

  • For num_threads > 0, return this value.

  • For num_threads < 0, return the maximal number of threads minus |num_threads + 1|. In particular num_threads = -1 will use as many threads as there are available cores on the machine.

  • For num_threads = 0 a ValueError is raised.

Parameters
num_threadsint or None

Desired number of threads to be used.

get_sphere

dipy.denoise.shift_twist_convolution.get_sphere(name='symmetric362')

provide triangulated spheres

Parameters
namestr

which sphere - one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’

Returns
spherea dipy.core.sphere.Sphere class instance

Examples

>>> import numpy as np
>>> from dipy.data import get_sphere
>>> sphere = get_sphere('symmetric362')
>>> verts, faces = sphere.vertices, sphere.faces
>>> verts.shape == (362, 3)
True
>>> faces.shape == (720, 3)
True
>>> verts, faces = get_sphere('not a sphere name') 
Traceback (most recent call last):
    ...
DataError: No sphere called "not a sphere name"

sf_to_sh

dipy.denoise.shift_twist_convolution.sf_to_sh(sf, sphere, sh_order=4, basis_type=None, full_basis=False, legacy=True, smooth=0.0)

Spherical function to spherical harmonics (SH).

Parameters
sfndarray

Values of a function on the given sphere.

sphereSphere

The points on which the sf is defined.

sh_orderint, optional

Maximum SH order in the SH fit. For sh_order, there will be (sh_order + 1) * (sh_order + 2) / 2 SH coefficients for a symmetric basis and (sh_order + 1) * (sh_order + 1) coefficients for a full SH basis.

basis_type{None, ‘tournier07’, ‘descoteaux07’}, optional

None for the default DIPY basis, tournier07 for the Tournier 2007 [R62f851f19ba8-2]_[R62f851f19ba8-3]_ basis, descoteaux07 for the Descoteaux 2007 [1] basis, (None defaults to descoteaux07).

full_basis: bool, optional

True for using a SH basis containing even and odd order SH functions. False for using a SH basis consisting only of even order SH functions.

legacy: bool, optional

True to use a legacy basis definition for backward compatibility with previous tournier07 and descoteaux07 implementations.

smoothfloat, optional

Lambda-regularization in the SH fit.

Returns
shndarray

SH coefficients representing the input function.

References

1

Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510.

2

Tournier J.D., Calamante F. and Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution. NeuroImage. 2007;35(4):1459-1472.

3

Tournier J-D, Smith R, Raffelt D, Tabbara R, Dhollander T, Pietsch M, et al. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage. 2019 Nov 15;202:116-137.

sh_to_sf

dipy.denoise.shift_twist_convolution.sh_to_sf(sh, sphere, sh_order=4, basis_type=None, full_basis=False, legacy=True)

Spherical harmonics (SH) to spherical function (SF).

Parameters
shndarray

SH coefficients representing a spherical function.

sphereSphere

The points on which to sample the spherical function.

sh_orderint, optional

Maximum SH order in the SH fit. For sh_order, there will be (sh_order + 1) * (sh_order + 2) / 2 SH coefficients for a symmetric basis and (sh_order + 1) * (sh_order + 1) coefficients for a full SH basis.

basis_type{None, ‘tournier07’, ‘descoteaux07’}, optional

None for the default DIPY basis, tournier07 for the Tournier 2007 [R72146e2b2f52-2]_[R72146e2b2f52-3]_ basis, descoteaux07 for the Descoteaux 2007 [1] basis, (None defaults to descoteaux07).

full_basis: bool, optional

True to use a SH basis containing even and odd order SH functions. Else, use a SH basis consisting only of even order SH functions.

legacy: bool, optional

True to use a legacy basis definition for backward compatibility with previous tournier07 and descoteaux07 implementations.

Returns
sfndarray

Spherical function values on the sphere.

References

1

Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Q-ball Imaging. Magn. Reson. Med. 2007;58:497-510.

2

Tournier J.D., Calamante F. and Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution. NeuroImage. 2007;35(4):1459-1472.

3

Tournier J-D, Smith R, Raffelt D, Tabbara R, Dhollander T, Pietsch M, et al. MRtrix3: A fast, flexible and open software framework for medical image processing and visualisation. NeuroImage. 2019 Nov 15;202:116-137.