denoise
denoise.adaptive_soft_matching
denoise.denspeed
denoise.enhancement_kernel
denoise.gibbs
denoise.localpca
denoise.nlmeans
denoise.nlmeans_block
denoise.noise_estimate
denoise.non_local_means
denoise.patch2self
denoise.pca_noise_estimate
denoise.shift_twist_convolution
EnhancementKernel
HemiSphere
Sphere
Version
partial
Version
denoise

Run benchmarks for module using nose. 

Run tests for module using nose. 
denoise.adaptive_soft_matching

Adaptive Soft Coefficient Matching 
denoise.denspeed
Determine the effective number of threads to be used for OpenMP calls 

Nonlocal means for denoising 3D images 

denoise.enhancement_kernel
Methods 


Points on the unit sphere. 

Points on the unit sphere. 

Return the ceiling of x as an Integral. 

Models electrostatic repulsion on the unit sphere 

provide triangulated spheres 
Accessor for tempfile.tempdir. 
denoise.gibbs
alias of 

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


Returns a process pool object 

Deprecate a renamed or removed function argument. 

Determine the effective number of processes for parallelization. 

Suppresses Gibbs ringing artefacts of images volumes. 
denoise.localpca

This class abstracts handling of a project's versions. 

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

General function to perform PCAbased denoising of diffusion datasets. 

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

Performs PCAbased denoising using the MarcenkoPastur distribution [1]. 

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

Nonlocal means for denoising 3D and 4D images 
Nonlocal means for denoising 3D images 
denoise.nlmeans_block
Applies the filter given by the convolution kernel 'h' columnwise to 'image', then subsamples by 2. 

NonLocal Means Denoising Using Blockwise Averaging 

Upsamples the columns of the input image by 2, then applies the convolution kernel 'h' (again, columnwise). 
denoise.noise_estimate

Multidimensional convolution. 

Standard deviation estimation from local patches 

Probabilistic Identification and Estimation of Noise (PIESNO). 
denoise.non_local_means
NonLocal Means Denoising Using Blockwise Averaging 


Nonlocal means for denoising 3D and 4D images, using 
denoise.patch2self

Return packagelike thing and module setup for package name 

Patch2Self Denoiser. 

Issue a warning, or maybe ignore it or raise an exception. 
denoise.pca_noise_estimate
PCA based local noise estimation. 
denoise.shift_twist_convolution
Perform the shifttwist convolution with the ODF data and the lookuptable of the kernel. 

Perform the shifttwist convolution with the ODF data and the lookuptable of the kernel. 

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


Spherical function to spherical harmonics (SH). 

Spherical harmonics (SH) to spherical function (SF). 
Run benchmarks for module using nose.
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’.
Verbosity value for benchmark outputs, in the range 110. Default is 1.
List with any extra arguments to pass to nosetests.
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” subdirectory 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
Run tests for module using nose.
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’.
Verbosity value for test outputs, in the range 110. Default is 1.
List with any extra arguments to pass to nosetests.
If True, run doctests in module. Default is False.
If True, report coverage of NumPy code. Default is False. (This requires the coverage module).
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.
Timing of individual tests with nosetimer
(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 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 Coefficient Matching
Combines two filtered 3Dimages at different resolutions and the original image. Returns the resulting combined image.
filtered image with optimized nonlocal means using a small block (suggested:3x3), which corresponds to a “high resolution” filter.
filtered image with optimized nonlocal means using a small block (suggested:5x5), which corresponds to a “low resolution” filter.
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
output denoised array which is of the same shape as that of the input
References
Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins. “Multiresolution NonLocal Means Filter for 3D MR Image Denoising” IET Image Processing, Institution of Engineering and Technology, 2011
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.
Desired number of threads to be used.
Nonlocal means for denoising 3D images
The array to be denoised
standard deviation of the noise estimated from the data
patch size is 2 x patch_radius + 1
. Default is 1.
block size is 2 x block_radius + 1
. Default is 5.
If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.
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.
the denoised arr
which has the same shape as arr
.
EnhancementKernel
Bases: object
Methods
Evaluate the kernel at position x relative to position y, with orientation r relative to orientation v. 

Return the computed lookup table. 

Return the orientations. 

Get the sphere corresponding with the orientations 
Compute a lookup table for the contextual enhancement kernel
Spatial diffusion
Angular diffusion
Diffusion time
Always compute the lookup table even if it is available in cache. Default is False.
Specify the number of orientations to be used with electrostatic repulsion, or provide a Sphere object. The default sphere is ‘repulsion100’.
Enable verbose mode.
References
J. Portegies, R. Duits. (2016) Fast implementations of contextual PDE’s for HARDI data processing in DIPY. ISMRM 2016 conference.
on the space of positions and orientations and their application to crossingpreserving smoothing of HARDI images. International Journal of Computer Vision, 92:231264.
(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
G. Girard, and R. Duits. (2015) Improving Fiber Alignment in HARDI by Combining Contextual PDE flow with Constrained Spherical Deconvolution. PLoS One.
HemiSphere
Bases: 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)
Vertices as xyz coordinates.
Vertices as spherical coordinates. Theta and phi are the inclination and azimuth angles respectively.
Vertices as xyz coordinates.
Indices into vertices that form triangular faces. If unspecified, the faces are computed using a Delaunay triangulation.
Edges between vertices. If unspecified, the edges are derived from the faces.
Angle in degrees. Vertices that are less than tol degrees apart are treated as duplicates.
See also
Methods

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

Create instance from a Sphere 

Create a full Sphere from a HemiSphere 

Create a more subdivided HemiSphere 
edges 

faces 

vertices 
Create a HemiSphere from points
Sphere
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)
Vertices as xyz coordinates.
Vertices as spherical coordinates. Theta and phi are the inclination and azimuth angles respectively.
Vertices as xyz coordinates.
Indices into vertices that form triangular faces. If unspecified, the faces are computed using a Delaunay triangulation.
Edges between vertices. If unspecified, the edges are derived from the faces.
Methods

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

Subdivides each face of the sphere into four new faces. 
edges 

faces 

vertices 
Find the index of the vertex in the Sphere closest to the input vector
A unit vector
The index into the Sphere.vertices array that gives the closest vertex (in angle).
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
The number of subdivisions to preform.
The subdivided sphere.
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.
Points on a unit sphere.
Number of iterations to run.
Using a smaller const could provide a more accurate result, but will need more iterations to converge.
Distributed points on a unit sphere.
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. Therefore, 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.
provide triangulated spheres
which sphere  one of: * ‘symmetric362’ * ‘symmetric642’ * ‘symmetric724’ * ‘repulsion724’ * ‘repulsion100’ * ‘repulsion200’
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"
Version
partial
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.
The old name of the argument.
None
, optionalThe new name of the argument. Set this to None to remove the
argument old_name
instead of renaming it.
The release at which the old argument became deprecated.
Last released version at which this function will still raise a deprecation warning. Versions higher than this will raise an error.
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.
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
.
Warning to be issued.
Error to be issued
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.
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 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.
Desired number of processes to be used.
Suppresses Gibbs ringing artefacts of images volumes.
Matrix containing one volume (3D) or multiple (4D) volumes of images.
Data axis corresponding to the number of acquired slices. Default is set to the third axis.
Number of neighbour points to access local TV (see note). Default is set to 3.
If True, the input data is replaced with results. Otherwise, returns a new array. Default is set to True.
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.
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 .. [Rfce522a872fd1] 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
Kellner E, Dhital B, Kiselev VG, Reisert M. Gibbsringing artifact removal based on local subvoxelshifts. Magn Reson Med. 2016 doi: 10.1002/mrm.26054.
Version
Bases: _BaseVersion
This class abstracts handling of a project’s versions.
A Version
instance is comparison aware and can be compared and
sorted using the standard Python interfaces.
>>> v1 = Version("1.0a5")
>>> v2 = Version("1.0")
>>> v1
<Version('1.0a5')>
>>> v2
<Version('1.0')>
>>> v1 < v2
True
>>> v1 == v2
False
>>> v1 > v2
False
>>> v1 >= v2
False
>>> v1 <= v2
True
base_version
The “base version” of the version.
dev
The development number of the version.
epoch
The epoch of the version.
is_devrelease
Whether this version is a development release.
is_postrelease
Whether this version is a postrelease.
is_prerelease
Whether this version is a prerelease.
local
The local version segment of the version.
major
The first item of release
or 0
if unavailable.
micro
The third item of release
or 0
if unavailable.
minor
The second item of release
or 0
if unavailable.
post
The postrelease number of the version.
pre
The prerelease segment of the version.
public
The public portion of the version.
release
The components of the “release” segment of the version.
Initialize a Version object.
version – The string representation of a version which will be parsed and normalized before use.
InvalidVersion – If the version
does not conform to PEP 440 in any way then this
exception will be raised.
The “base version” of the version.
>>> Version("1.2.3").base_version
'1.2.3'
>>> Version("1.2.3+abc").base_version
'1.2.3'
>>> Version("1!1.2.3+abc.dev1").base_version
'1!1.2.3'
The “base version” is the public version of the project without any pre or post release markers.
The development number of the version.
>>> print(Version("1.2.3").dev)
None
>>> Version("1.2.3.dev1").dev
1
The epoch of the version.
>>> Version("2.0.0").epoch
0
>>> Version("1!2.0.0").epoch
1
Whether this version is a development release.
>>> Version("1.2.3").is_devrelease
False
>>> Version("1.2.3.dev1").is_devrelease
True
Whether this version is a postrelease.
>>> Version("1.2.3").is_postrelease
False
>>> Version("1.2.3.post1").is_postrelease
True
Whether this version is a prerelease.
>>> Version("1.2.3").is_prerelease
False
>>> Version("1.2.3a1").is_prerelease
True
>>> Version("1.2.3b1").is_prerelease
True
>>> Version("1.2.3rc1").is_prerelease
True
>>> Version("1.2.3dev1").is_prerelease
True
The local version segment of the version.
>>> print(Version("1.2.3").local)
None
>>> Version("1.2.3+abc").local
'abc'
The third item of release
or 0
if unavailable.
>>> Version("1.2.3").micro
3
>>> Version("1").micro
0
The second item of release
or 0
if unavailable.
>>> Version("1.2.3").minor
2
>>> Version("1").minor
0
The postrelease number of the version.
>>> print(Version("1.2.3").post)
None
>>> Version("1.2.3.post1").post
1
The prerelease segment of the version.
>>> print(Version("1.2.3").pre)
None
>>> Version("1.2.3a1").pre
('a', 1)
>>> Version("1.2.3b1").pre
('b', 1)
>>> Version("1.2.3rc1").pre
('rc', 1)
The public portion of the version.
>>> Version("1.2.3").public
'1.2.3'
>>> Version("1.2.3+abc").public
'1.2.3'
>>> Version("1.2.3+abc.dev1").public
'1.2.3'
The components of the “release” segment of the version.
>>> Version("1.2.3").release
(1, 2, 3)
>>> Version("2.0.0").release
(2, 0, 0)
>>> Version("1!2.0.0.post0").release
(2, 0, 0)
Includes trailing zeroes but not the epoch or any prerelease / development / postrelease suffixes.
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 λ (ith entry of w) and its eigenvector vi
(ith 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.
A complex Hermitian or real symmetric matrix whose eigenvalues and eigenvectors will be computed.
A complex Hermitian or real symmetric definite positive matrix in. If omitted, identity matrix is assumed.
Whether the pertinent array data is taken from the lower or upper
triangle of a
and, if applicable, b
. (Default: lower)
Whether to calculate only eigenvalues and no eigenvectors. (Default: both are calculated)
If provided, this twoelement iterable defines the start and the end
indices of the desired eigenvalues (ascending order and 0indexed).
To return only the second smallest to fifth smallest eigenvalues,
[1, 4]
is used. [n3, n1]
returns the largest three. Only
available with “evr”, “evx”, and “gvx” drivers. The entries are
directly converted to integers via int()
.
If provided, this twoelement iterable defines the halfopen 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.
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. The default for standard problems is “evr”. For generalized problems, “gvd” is used for full set, and “gvx” for subset requested cases.
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.
Whether to overwrite data in a
(may improve performance). Default
is False.
Whether to overwrite data in b
(may improve performance). Default
is False.
Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, nontermination) if the inputs do contain infinities or NaNs.
Deprecated since version 1.5.0: eigh keyword argument turbo is deprecated in favour of
driver=gvd
keyword instead and will be removed in SciPy
1.12.0.
Deprecated since version 1.5.0: eigh keyword argument eigvals is deprecated in favour of subset_by_index keyword instead and will be removed in SciPy 1.12.0.
The N (1<=N<=M) selected eigenvalues, in ascending order, each repeated according to its multiplicity.
(if eigvals_only == False
)
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 nonsymmetric 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
occasions 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 respect 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
>>> import numpy as np
>>> 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.69199443e07, 9.11938152e+00])
Request the second smallest 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)
General function to perform PCAbased denoising of diffusion datasets.
Array of data to be denoised. The dimensions are (X, Y, Z, N), where N are the diffusion gradient directions. The first 3 dimension must have size >= 2 * patch_radius + 1 or size = 1.
Standard deviation of the noise estimated from the data. If no sigma is given, this will be estimated based on random matrix theory [1],[R4bed7265c9342]_
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.
The radius of the local patch to be taken around each voxel (in voxels). Default: 2 (denoise in blocks of 5x5x5 voxels).
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.
Thresholding of PCA eigenvalues is done by nulling out eigenvalues that are smaller than:
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.
If true, the Standard deviation of the noise will be returned. Default: False.
The dtype for the output array. Default: output has the same dtype as the input.
If true, suppress warning caused by patch_size < arr.shape[1]. Default: False.
This is the denoised array of the same size as that of the input data, clipped to nonnegative values.
References
Veraart J, Novikov DS, Christiaens D, Adesaron B, Sijbers, Fieremans E, 2016. Denoising of Diffusion MRI using random matrix theory. Neuroimage 142:394406. doi: 10.1016/j.neuroimage.2016.08.016
Veraart J, Fieremans E, Novikov DS. 2016. Diffusion MRI noise mapping using random matrix theory. Magnetic Resonance in Medicine. doi: 10.1002/mrm.26059.
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
Performs local PCA denoising according to Manjon et al. [1].
Array of data to be denoised. The dimensions are (X, Y, Z, N), where N are the diffusion gradient directions.
Standard deviation of the noise estimated from the data.
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.
The radius of the local patch to be taken around each voxel (in voxels). Default: 2 (denoise in blocks of 5x5x5 voxels).
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.
Thresholding of PCA eigenvalues is done by nulling out eigenvalues that are smaller than:
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 MarcenkoPastur distribution [2]. Default: 2.3 (according to [1])
The dtype for the output array. Default: output has the same dtype as the input.
If true, suppress warning caused by patch_size < arr.shape[1]. Default: False.
This is the denoised array of the same size as that of the input data, clipped to nonnegative values
References
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
Veraart J, Novikov DS, Christiaens D, Adesaron B, Sijbers, Fieremans E, 2016. Denoising of Diffusion MRI using random matrix theory. Neuroimage 142:394406. doi: 10.1016/j.neuroimage.2016.08.016
Performs PCAbased denoising using the MarcenkoPastur distribution [1].
Array of data to be denoised. The dimensions are (X, Y, Z, N), where N are the diffusion gradient directions.
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.
The radius of the local patch to be taken around each voxel (in voxels). Default: 2 (denoise in blocks of 5x5x5 voxels).
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.
If true, a noise standard deviation estimate based on the MarcenkoPastur distribution is returned [2]. Default: False.
The dtype for the output array. Default: output has the same dtype as the input.
If true, suppress warning caused by patch_size < arr.shape[1]. Default: False.
This is the denoised array of the same size as that of the input data, clipped to nonnegative values
Estimate of the spatial varying standard deviation of the noise
References
Veraart J, Novikov DS, Christiaens D, Adesaron B, Sijbers, Fieremans E, 2016. Denoising of Diffusion MRI using random matrix theory. Neuroimage 142:394406. doi: 10.1016/j.neuroimage.2016.08.016
Veraart J, Fieremans E, Novikov DS. 2016. Diffusion MRI noise mapping using random matrix theory. Magnetic Resonance in Medicine. doi: 10.1002/mrm.26059.
Nonlocal means for denoising 3D and 4D images
The array to be denoised
standard deviation of the noise estimated from the data
patch size is 2 x patch_radius + 1
. Default is 1.
block size is 2 x block_radius + 1
. Default is 5.
If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.
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.
the denoised arr
which has the same shape as arr
.
References
Descoteaux, Maxime and WiestDaesslé, Nicolas and Prima, Sylvain and Barillot, Christian and Deriche, Rachid Impact of Rician Adapted NonLocal Means Filtering on HARDI, MICCAI 2008
Nonlocal means for denoising 3D images
The array to be denoised
standard deviation of the noise estimated from the data
patch size is 2 x patch_radius + 1
. Default is 1.
block size is 2 x block_radius + 1
. Default is 5.
If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.
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.
the denoised arr
which has the same shape as arr
.
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.
the input image to be filtered
the convolution kernel
NonLocal Means Denoising Using Blockwise Averaging
the input image, corrupted with rician noise
the input mask
similar patches in the nonlocal means are searched for locally, inside a cube of side 2*v+1 centered at each voxel of interest.
the size of the block to be used (2*f+1)x(2*f+1)x(2*f+1) in the blockwise nonlocal means implementation (the Coupe’s proposal).
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
If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.
the denoised output which has the same shape as input image.
References
“An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic Resonance Images” IEEE Transactions on Medical Imaging, 27(4):425441, 2008
“Multiresolution NonLocal Means Filter for 3D MR Image Denoising” IET Image Processing, Institution of Engineering and Technology, 2011
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.
the input image to be filtered
the convolution kernel
Multidimensional convolution.
The array is convolved with the given kernel.
The input array.
Array of weights, same number of dimensions as input
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.
The mode parameter determines how the input array is extended beyond its boundaries. Default is ‘reflect’. Behavior for each valid value is as follows:
The input is extended by reflecting about the edge of the last pixel. This mode is also sometimes referred to as halfsample symmetric.
The input is extended by filling all values beyond the edge with the same constant value, defined by the cval parameter.
The input is extended by replicating the last pixel.
The input is extended by reflecting about the center of the last pixel. This mode is also sometimes referred to as wholesample symmetric.
The input is extended by wrapping around to the opposite edge.
For consistency with the interpolation functions, the following mode names can also be used:
This is a synonym for ‘reflect’.
This is a synonym for ‘constant’.
This is a synonym for ‘wrap’.
Value to fill past edges of input if mode is ‘constant’. Default is 0.0
Controls the origin of the input signal, which is where the filter is centered to produce the first element of the output. Positive values shift the filter to the right, and negative values shift the filter to the left. Default is 0.
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+kj} W_j}\), where W is the weights kernel, j is the ND 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.
>>> import numpy as np
>>> 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]])
Standard deviation estimation from local patches
The array to be estimated
If True, uses all voxels for the estimation, otherwise, only nonzeros voxels are used. Useful if the background is masked by the scanner.
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.
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 coildependent 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 undenoised for spatially varying noise
profiles. Consider using piesno()
to estimate sigma instead if visual
inaccuracies are apparent in the denoised result.
References
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), 31722.
Coupe, P., Yger, P., Prima, S., Hellier, P., Kervrann, C., Barillot,
C., 2008. An optimized blockwise nonlocal means denoising filter for 3D magnetic resonance images, IEEE Trans. Med. Imaging 27, 42541.
Probabilistic Identification and Estimation of Noise (PIESNO).
The magnitude signals to analyse. The last dimension must contain the same realisation of the volume, such as dMRI or fMRI data.
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.
Probabilistic estimation threshold for the gamma function.
number of initial estimates for sigma to try.
Maximum number of iterations to execute if convergence is not reached.
Tolerance for the convergence criterion. Convergence is reached if two subsequent estimates are smaller than eps.
If True, return a mask identifying all the pure noise voxel that were found.
The estimated standard deviation of the gaussian noise.
A boolean mask indicating the voxels identified as pure noise.
Notes
This function assumes two things : 1. The data has a noisy, nonmasked 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
Koay CG, Ozarslan E and Pierpaoli C.
“Probabilistic Identification and Estimation of Noise (PIESNO): A selfconsistent approach and its applications in MRI.” Journal of Magnetic Resonance 2009; 199: 94103.
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: 108119.
NonLocal Means Denoising Using Blockwise Averaging
the input image, corrupted with rician noise
the input mask
similar patches in the nonlocal means are searched for locally, inside a cube of side 2*v+1 centered at each voxel of interest.
the size of the block to be used (2*f+1)x(2*f+1)x(2*f+1) in the blockwise nonlocal means implementation (the Coupe’s proposal).
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
If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.
the denoised output which has the same shape as input image.
References
“An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic Resonance Images” IEEE Transactions on Medical Imaging, 27(4):425441, 2008
“Multiresolution NonLocal Means Filter for 3D MR Image Denoising” IET Image Processing, Institution of Engineering and Technology, 2011
blockwise averaging approach
The array to be denoised
standard deviation of the noise estimated from the data
patch size is 2 x patch_radius + 1
. Default is 1.
block size is 2 x block_radius + 1
. Default is 5.
If True the noise is estimated as Rician, otherwise Gaussian noise is assumed.
the denoised arr
which has the same shape as arr
.
References
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):425441, 2008
Pierrick Coupe, Jose Manjon, Montserrat Robles, Louis Collins. Adaptive Multiresolution NonLocal Means Filter for 3D MR Image Denoising IET Image Processing, Institution of Engineering and Technology, 2011
Return packagelike thing and module setup for package name
package name
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.
TripWire
instanceIf we can import the package, return it. Otherwise return an object raising an error when accessed
True if import for package was successful, false otherwise
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 Denoiser.
The 4D noisy DWI data to be denoised.
Array of the bvals from the DWI acquisition
The radius of the local patch to be taken around each voxel (in voxels). Default: 0 (denoise in blocks of 1x1x1 voxels).
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 ScikitLearn: sklearn.linear_model.LinearRegression, sklearn.linear_model.Lasso or sklearn.linear_model.Ridge and other objects that inherit from sklearn.base.RegressorMixin. Default: ‘ols’.
Threshold for considering volumes as b0.
The dtype for the output array. Default: output has the same dtype as the input.
Regularization parameter only for ridge regression model.
Show progress of Patch2Self and time taken.
Skips denoising b0 volumes if set to False.
Sets negative values after denoising to 0 using np.clip.
Shifts the distribution of intensities per volume to give nonnegative values
This is the denoised array of the same size as that of the input data, clipped to nonnegative values.
References
Denoising Diffusion MRI with Selfsupervised Learning, Advances in Neural Information Processing Systems 33 (2020)
PCA based local noise estimation.
the input dMRI data. The first 3 dimension must have size >= 2 * patch_radius + 1 or size = 1.
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.
The radius of the local patch to be taken around each voxel (in voxels). Default: 1 (estimate noise in blocks of 3x3x3 voxels).
Whether to correct for bias due to Rician noise. This is an implementation of equation 8 in [1].
Radius of a Gaussian smoothing filter to apply to the noise estimate before returning. Default: 2.
The local noise standard deviation estimate.
References
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.
Perform the shifttwist convolution with the ODF data and the lookuptable of the kernel.
The ODF data in spherical harmonics format
The 5D lookup table
Maximal spherical harmonics order
Reduced convolution in one direction only for testing
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.
Apply maxnormalization to the output such that its value range matches the input ODF data.
The ODF data after convolution enhancement in spherical harmonics format
References
J. Portegies, R. Duits. (2016) Fast implementations of contextual PDE’s for HARDI data processing in DIPY. ISMRM 2016 conference.
on the space of positions and orientations and their application to crossingpreserving smoothing of HARDI images. International Journal of Computer Vision, 92:231264.
(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
and R. Duits. (2015) Improving Fiber Alignment in HARDI by Combining Contextual PDE flow with Constrained Spherical Deconvolution. PLoS One.
Perform the shifttwist convolution with the ODF data and the lookuptable of the kernel.
The ODF data sampled on a sphere
The 5D lookup table
Reduced convolution in one direction only for testing
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.
Apply maxnormalization to the output such that its value range matches the input ODF data.
The ODF data after convolution enhancement, sampled on a sphere
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.
Desired number of threads to be used.
Spherical function to spherical harmonics (SH).
Values of a function on the given sphere
.
The points on which the sf is defined.
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.
None
for the default DIPY basis,
tournier07
for the Tournier 2007 [R62f851f19ba82]_[R62f851f19ba83]_ basis,
descoteaux07
for the Descoteaux 2007 [1] basis,
(None
defaults to descoteaux07
).
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.
True to use a legacy basis definition for backward compatibility
with previous tournier07
and descoteaux07
implementations.
Lambdaregularization in the SH fit.
SH coefficients representing the input function.
References
Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Qball Imaging. Magn. Reson. Med. 2007;58:497510.
Tournier J.D., Calamante F. and Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: Nonnegativity constrained superresolved spherical deconvolution. NeuroImage. 2007;35(4):14591472.
Tournier JD, 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:116137.
Spherical harmonics (SH) to spherical function (SF).
SH coefficients representing a spherical function.
The points on which to sample the spherical function.
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.
None
for the default DIPY basis,
tournier07
for the Tournier 2007 [R72146e2b2f522]_[R72146e2b2f523]_ basis,
descoteaux07
for the Descoteaux 2007 [1] basis,
(None
defaults to descoteaux07
).
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.
True to use a legacy basis definition for backward compatibility
with previous tournier07
and descoteaux07
implementations.
Spherical function values on the sphere
.
References
Descoteaux, M., Angelino, E., Fitzgibbons, S. and Deriche, R. Regularized, Fast, and Robust Analytical Qball Imaging. Magn. Reson. Med. 2007;58:497510.
Tournier J.D., Calamante F. and Connelly A. Robust determination of the fibre orientation distribution in diffusion MRI: Nonnegativity constrained superresolved spherical deconvolution. NeuroImage. 2007;35(4):14591472.
Tournier JD, 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:116137.