Combines two filtered 3D-images at different resolutions and the original
image. Returns the resulting combined image.
Parameters:
imandarray
Original (unfiltered) 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 [1] 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
DWI bias field correction via multi-resolution Legendre polynomial regression.
Estimates the bias field from the mean b0 volume in log space using
coarse-to-fine Legendre polynomial regression, then applies the estimated
field to all DWI volumes.
Parameters:
datandarray
4D DWI data (X, Y, Z, N).
gtabGradientTable
Gradient table.
maskndarray, optional
3D binary brain mask. Auto-computed via median_otsu if None.
orderint, optional
Maximum Legendre polynomial order (terms where i+j+k <= order).
pyramid_levelstuple of int, optional
Downsampling factors for coarse-to-fine pyramid (descending order).
n_iterint, optional
Reweighting iterations per pyramid level.
lambda_regfloat, optional
Ridge regularization strength.
robustbool, optional
Apply Tukey biweight robust reweighting.
gradient_weightingbool, optional
Apply gradient-based edge suppression.
zero_backgroundbool, optional
If True, set the bias field to 1.0 (no correction) outside the brain
mask. If False, the raw extrapolated field values are preserved in
the returned bias_field array. Has no effect on the corrected DWI
data (background voxels are always zeroed by the brain mask).
Top-level DWI bias field correction via regression.
Estimates a smooth multiplicative bias field from the mean b0 volume
using polynomial or B-spline regression in log space, then applies the
correction uniformly to all DWI volumes.
Parameters:
datandarray
4D DWI data (X, Y, Z, N).
gtabGradientTable
Gradient table.
maskndarray, optional
3D binary brain mask. If None, computed via median_otsu.
"bspline": Cubic B-spline regression — more flexible.
"auto": Run both methods and return the one with lower
Coefficient of Variation within the brain mask. The chosen method
is logged at INFO level.
orderint, optional
Maximum Legendre polynomial degree (used only for method=”poly”).
n_control_pointstuple of int, optional
Control grid dimensions at finest level (used only for
method=”bspline”).
pyramid_levelstuple of int, optional
Downsampling factors for coarse-to-fine pyramid (descending order).
n_iterint, optional
Reweighting iterations per pyramid level.
lambda_regfloat, optional
Ridge regularization strength.
robustbool, optional
Apply Tukey biweight robust reweighting at each level.
gradient_weightingbool, optional
Weight regression by edge-suppression map derived from the
image gradient.
return_bias_fieldbool, optional
If True, return the bias field alongside the corrected data.
zero_backgroundbool, optional
If True, set the bias field to 1.0 (no correction) outside the brain
mask. If False, the raw extrapolated field values are preserved in
the returned bias_field array. Has no effect on the corrected DWI
data (background voxels are always zeroed by the brain mask).
Returns:
correctedndarray
Bias-corrected DWI, same dtype as input.
bias_fieldndarray
3D multiplicative bias field (only returned if
return_bias_field=True).
Suppresses Gibbs ringing artefacts of images volumes.
See [2] and [3] for
further details about the method.
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.
n_pointsint, optional
Number of neighbour points to access local TV (see note).
inplacebool, optional
If True, the input data is replaced with results. Otherwise, returns
a new array.
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.
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. The first 3 dimension must have
size >= 2 * patch_radius + 1 or size = 1.
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
[4], [5].
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). E.g. patch_radius=2 gives 5x5x5 patches.
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\)[6]), or automatically calculated using random
matrix theory (in case that \(\tau_{factor}\) is set to None).
return_sigmabool, optional
If true, the Standard deviation of the noise will be returned.
out_dtypestr or dtype, optional
The dtype for the output array. Default: output has the same dtype as
the input.
suppress_warningbool, optional
If true, suppress warning caused by patch_size < arr.shape[-1].
Returns:
denoised_arr4D array
This is the denoised array of the same size as that of the input data,
clipped to non-negative values.
The method follows the algorithm in Manjón et al.[6].
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 not given,
calculate using method in Manjón et al.[6].
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). E.g. patch_radius=2 gives 5x5x5 patches.
gtab: gradient table object (optional if sigma is provided)
gradient information for the data gives us the bvals and bvecs of
diffusion data, which is needed to calculate noise level if sigma is
not provided.
patch_radius_sigmaint, optional
The radius of the local patch to be taken around each voxel (in
voxels) for estimating sigma. E.g. patch_radius_sigma=2 gives
5x5x5 patches.
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 [5]. Default: 2.3
according to Manjón et al.[6].
return_sigmabool, optional
If true, a noise standard deviation estimate based on the
Marcenko-Pastur distribution is returned [5].
correct_biasbool, optional
Whether to correct for bias due to Rician noise. This is an
implementation of equation 8 in [6].
out_dtypestr or dtype, optional
The dtype for the output array. Default: output has the same dtype as
the input.
suppress_warningbool, optional
If true, suppress warning caused by patch_size < arr.shape[-1].
Returns:
denoised_arr4D array
This is the denoised array of the same size as that of the input data,
clipped to non-negative values
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). E.g. patch_radius=2 gives 5x5x5 patches.
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 [4].
out_dtypestr or dtype, optional
The dtype for the output array. Default: output has the same dtype as
the input.
suppress_warningbool, optional
If true, suppress warning caused by patch_size < arr.shape[-1].
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
Non-local means denoising for 3D and 4D images with selectable algorithms.
This implementation provides two different algorithms for non-local means denoising:
the classic approach and an improved blockwise algorithm. The blockwise method
offers better coordinate handling, memory efficiency, and performance through
advanced parallelization and statistical pre-filtering.
See [7] for further details about the classic method.
See [8] and [1] for further details
about the blockwise method.
Parameters:
arr3D or 4D ndarray
The array to be denoised. For 3D arrays, shape should be (height, width, depth).
For 4D arrays, shape should be (height, width, depth, volumes) where the last
dimension represents different volumes (e.g., DWI directions).
sigmafloat, 1D ndarray, or 3D ndarray
Standard deviation of the noise estimated from the data. For 3D arrays,
this should be a scalar. For 4D arrays, this can be either a scalar (same noise
level for all volumes), a 1D array with length equal to the number
of volumes, or a 3D array with shape arr.shape[:3].
mask3D ndarray, optional
Binary mask indicating which voxels to process. Should have shape
(height, width, depth). Voxels with mask value 0 are set to 0 in output.
If None, all voxels are processed.
patch_radiusint, optional
Radius for similarity search neighborhoods. Patches of size
(2*patch_radius + 1)³ around each voxel are compared to find similar
structures.
block_radiusint, optional
Radius for weighted averaging blocks. Each block has size
(2*block_radius + 1)³. Larger blocks provide more smoothing but
increase computational cost. If None, defaults are:
method=’classic’: 5 (11*11*11 blocks)
method=’blockwise’: 2 (5*5*5 blocks)
ricianbool, optional
If True, assumes Rician noise model (appropriate for magnitude MRI data).
If False, assumes Gaussian noise model.
num_threadsint, optional
Number of OpenMP threads to use for parallel processing. If None,
uses all available CPU threads. Set to 1 to disable parallel processing.
methodstr, optional
Algorithm method to use:
‘blockwise’: Improved algorithm with better coordinate handling,
memory efficiency, and statistical pre-filtering
‘classic’: Original algorithm with traditional implementation
Added in version 1.12.0.
Returns:
denoised_arrndarray
The denoised array with the same shape and dtype as the input arr.
Values are clipped to non-negative range for Rician noise model.
Notes
Due to coordinate bug fixes in the blockwise method, equivalent denoising
quality may require different parameters between methods:
- Classic patch_radius=3 ≈ Blockwise patch_radius=2
- Block_radius can be smaller for blockwise due to efficiency improvements
For 4D inputs with method='blockwise', a 3D sigma map is accepted
and passed through for blockwise processing.
References
Examples
>>> importnumpyasnp>>> fromdipy.denoise.nlmeansimportnlmeans>>> # Create synthetic noisy data>>> clean_data=np.random.rand(50,50,50)*100>>> noisy_data=clean_data+np.random.randn(50,50,50)*10>>> # Denoise a 3D image with default blockwise method>>> denoised_data=nlmeans(noisy_data,sigma=10.0)>>> # Use classic method for compatibility:>>> denoised_classic=nlmeans(noisy_data,sigma=10.0,method='classic')>>> # Denoise 4D DWI data with different noise levels per volume:>>> dwi_data=np.random.rand(64,64,40,30)*1000# 30 DWI directions>>> noise_levels=np.linspace(50,100,30)# Varying noise>>> denoised_dwi=nlmeans(dwi_data,sigma=noise_levels)>>> # Denoise 4D DWI data with a 3D sigma map (e.g., from PIESNO):>>> sigma_map=np.ones((64,64,40))*30# one value per spatial voxel>>> denoised_piesno=nlmeans(dwi_data,sigma=sigma_map)
dipy.denoise.noise_estimate.piesno(data, N, *, alpha=0.01, step=100, itermax=100, eps=1e-05, return_mask=False)[source]¶
Probabilistic Identification and Estimation of Noise (PIESNO).
See [9] and [10] for further details
about the method.
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.
stepint
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.
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, optional
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
[11] 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 from
Koay and Basser[11] (see equation 18 in [11]) with
theta = 0. Since this function was introduced in [8] 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.
See [12] for further details about the method.
See [13] for further details about the new method.
Parameters:
datandarray
The 4D noisy DWI data to be denoised.
bvalsarray of shape (N,)
Array of the bvals from the DWI acquisition
patch_radiusint or array of shape (3,), optional
The radius of the local patch to be taken around each voxel (in
voxels).
modelstring, or sklearn.base.RegressorMixin, optional
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.
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.
verbosebool, optional
Show progress of Patch2Self and time taken.
b0_denoisingbool, optional
Skips denoising b0 volumes if set to False.
clip_negative_valsbool, optional
Sets negative values after denoising to 0 using np.clip.
shift_intensitybool, optional
Shifts the distribution of intensities per volume to give
non-negative values.
tmp_dirstr, optional
The directory to save the temporary files. If None, the temporary
files are saved in the system’s default temporary directory.
versionint, optional
Version 1 or 3 of Patch2Self to use.
Returns:
denoised arrayndarray
This is the denoised array of the same size as that of the input data,
clipped to non-negative values.