Filtering » Deconvolution module

#include "diplib/deconvolution.h"

void dip::WienerDeconvolution( dip::Image const& in, dip::Image const& psf, dip::Image const& signalPower, dip::Image const& noisePower, dip::Image& out, dip::StringSet const& options = {S::PAD})

Wiener deconvolution using estimates of signal and noise power spectra.

Assuming some original image , a known convolution kernel (given by psf), a noise realization , and an observed image (given by in), the Wiener deconvolution filter is the linear filter that, when convolved with , yields an image such that the mean square error between and is minimized.

Finding (returned in out) requires knowledge of the power spectra of the signal and the noise. The Wiener deconvolution filter is defined in the frequency domain as

where is the Fourier transform of in, is the Fourier transform of psf, is signalPower, and is noisePower. These and are typically not known, but:

• signalPower can be estimated as the Fourier transform of the autocorrelation of in. If a raw image is passed for this argument (dip::Image{}), then it will be computed as such.

• noisePower can be estimated as a flat function, assuming white noise. A 0D image can be given here, it will be expanded to the size of the other images. noisePower should not be zero anywhere, as that might lead to division by zero and consequently meaningless results.

The other syntax for dip::WienerDeconvolution takes an estimate of the noise-to-signal ratio instead of the signal and noise power spectra. Note that can be rewritten as

where is the noise-to-signal ratio. If is a constant, then the Wiener deconvolution filter is equivalent to the Tikhonov regularized inverse filter.

psf is given in the spatial domain, and will be zero-padded to the size of in and Fourier transformed. The PSF (point spread function) should sum to one in order to preserve the mean image intensity. If the OTF (optical transfer function, the Fourier transform of the PSF) is known, it is possible to pass that as psf; add the string "OTF" to options.

All input images must be real-valued and scalar, except if the OFT is given instead of the PSF, in which case psf could be complex-valued.

If "pad" is in options, then in is padded by the size of psf in all directions (padded area is filled by mirroring at the image border). This significantly reduces the effects of the periodicity of the frequency-domain convolution. "pad" cannot be combined with "OTF".

void dip::WienerDeconvolution( dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 1e-4, dip::StringSet const& options = {S::PAD})

Wiener deconvolution using an estimate of noise-to-signal ratio.

See the description of the function with the same name above. The difference here is that a single number, regularization, is given instead of the signal and noise power spectra. We then set (the noise-to-signal ratio) to regularization * dip::Maximum(P), with P equal to .

This formulation of the Wiener deconvolution filter is equivalent to the Tikhonov regularized inverse filter.

void dip::TikhonovMiller( dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 0.1, dip::StringSet const& options = {S::PAD})

Tikhonov-Miller deconvolution.

Assuming some original image , a known convolution kernel (given by psf), a noise realization , and an observed image (given by in), the Tikhonov-Miller deconvolution filter is the linear filter that, when convolved with , yields an image that minimizes the Tikhonov functional,

where is the regularization parameter (given by regularization), and is the regularization kernel, for which we use an ideal Laplace kernel here. is returned in out.

In the frequency domain, the Tikhonov-Miller deconvolution filter is defined as

where is the Fourier transform of in, is the Fourier transform of psf, and is the regularization kernel in the frequency domain.

psf is given in the spatial domain, and will be zero-padded to the size of in and Fourier transformed. The PSF (point spread function) should sum to one in order to preserve the mean image intensity. If the OTF (optical transfer function, the Fourier transform of the PSF) is known, it is possible to pass that as psf; add the string "OTF" to options.

All input images must be real-valued and scalar, except if the OFT is given instead of the PSF, in which case psf could be complex-valued.

If "pad" is in options, then in is padded by the size of psf in all directions (padded area is filled by mirroring at the image border). This significantly reduces the effects of the periodicity of the frequency-domain convolution. "pad" cannot be combined with "OTF".

void dip::IterativeConstrainedTikhonovMiller( dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 0.1, dip::dfloat tolerance = 1e-6, dip::uint maxIterations = 30, dip::dfloat stepSize = 0.0, dip::StringSet const& options = {S::PAD})

Iterative Constrained Tikhonov-Miller (ICTM) deconvolution.

Assuming some original image , a known convolution kernel (given by psf), a noise realization , and an observed image (given by in), ICTM deconvolution finds the (returned in out) that minimizes the Tikhonov functional,

where is the regularization parameter (given by regularization), and is the regularization kernel, for which we use an ideal Laplace kernel here. is returned in out.

If stepSize is 0 (the default), ICTM uses the conjugate gradient method to estimate . In this case, it uses the results of Verveer and Jovin to estimate the optimal step size for each step.

If a positive step size is given (a value in the range (0, 1]), then ICTM uses gradient descent (with steepest descent) and a fixed step size. This is much simpler code, with quicker steps, but converges much more slowly and can even diverge under certain circumstances. It is provided here because this is a common implementation in other software packages; we do not recommend using it.

The iterative algorithm is stopped when the maximum difference of between two steps (ignoring the regularization term) is less than tolerance times the maximum absolute value in .

maxIterations provides an additional stopping condition, in case the algorithm does not converge quickly enough. In a way, providing a small maximum number of iterations is an additional form of regularization. Setting maxIterations to 0 runs the algorithm until convergence.

psf is given in the spatial domain, and will be zero-padded to the size of in and Fourier transformed. The PSF (point spread function) should sum to one in order to preserve the mean image intensity. If the OTF (optical transfer function, the Fourier transform of the PSF) is known, it is possible to pass that as psf; add the string "OTF" to options.

All input images must be real-valued and scalar, except if the OFT is given instead of the PSF, in which case psf could be complex-valued.

If "pad" is in options, then in is padded by the size of psf in all directions (padded area is filled by mirroring at the image border). This significantly reduces the effects of the periodicity of the frequency-domain convolution. "pad" cannot be combined with "OTF".

void dip::RichardsonLucy( dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 0.0, dip::uint nIterations = 30, dip::StringSet const& options = {S::PAD})

Richardson-Lucy (RL) deconvolution, also sometimes called the expectation maximization (EM) method.

Assuming some original image , a known convolution kernel (given by psf), and an observed image (given by in), where is Poisson noise with mean , RL deconvolution finds the (returned in out) with maximal likelihood, given by

This is the basic, non-regularized Richardson-Lucy deconvolution, which requires regularization to be set to 0.

The nIterations parameter serves as regularization, the iterative process must be stopped before the noise gets amplified too much. Even when using the regularization parameter, there is no ideal way to see if the algorithm has converged.

If regularization is positive, total variation (TV) regularization is added, according to Dey et al. In this case, a term is added to the expression above, with equal to regularization. This should be a small value, 0.01 is a good start point. Note that TV regularization tends to introduce a stair-case effect, as it penalizes slow transitions but allows sharp jumps.

psf is given in the spatial domain, and will be zero-padded to the size of in and Fourier transformed. The PSF (point spread function) should sum to one in order to preserve the mean image intensity. If the OTF (optical transfer function, the Fourier transform of the PSF) is known, it is possible to pass that as psf; add the string "OTF" to options.

All input images must be real-valued and scalar, except if the OFT is given instead of the PSF, in which case psf could be complex-valued.

If "pad" is in options, then in is padded by the size of psf in all directions (padded area is filled by mirroring at the image border). This significantly reduces the effects of the periodicity of the frequency-domain convolution. "pad" cannot be combined with "OTF".

void dip::FastIterativeShrinkageThresholding( dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 0.1, dip::dfloat tolerance = 1e-6, dip::uint maxIterations = 30, dip::uint nScales = 3, dip::StringSet const& options = {S::PAD})

Fast Iterative Shrinkage-Thresholding (FISTA) deconvolution.

Assuming some original image , a known convolution kernel (given by psf), a noise realization , and an observed image (given by in), FISTA deconvolution finds the (returned in out) that minimizes the functional

where is the regularization parameter (given by regularization), and is a wavelet transform of . The regularization is applied in some wavelet domain, assuming that the image has a sparse representation in the wavelet domain. We use the Haar wavelet, due to its computational simplicity (it is also the wavelet used by Beck and Teboulle). is returned in out.

The iterative algorithm is stopped when the maximum difference of between two steps (ignoring the regularization term) is less than tolerance times the maximum absolute value in .

maxIterations provides an additional stopping condition, in case the algorithm does not converge quickly enough. In a way, providing a small maximum number of iterations is an additional form of regularization. Setting maxIterations to 0 runs the algorithm until convergence.

nScales determines how many scales of the Haar wavelet to compute. It defaults to 3, as used by Beck and Teboulle. Increasing this value might be useful for very large images.

psf is given in the spatial domain, and will be zero-padded to the size of in and Fourier transformed. The PSF (point spread function) should sum to one in order to preserve the mean image intensity. If the OTF (optical transfer function, the Fourier transform of the PSF) is known, it is possible to pass that as psf; add the string "OTF" to options.

All input images must be real-valued and scalar, except if the OFT is given instead of the PSF, in which case psf could be complex-valued.

If "pad" is in options, then in is padded by the size of psf in all directions (padded area is filled by mirroring at the image border). This significantly reduces the effects of the periodicity of the frequency-domain convolution. "pad" cannot be combined with "OTF".