module #include "diplib/deconvolution.h"
Deconvolution Deconvolution algorithms (inverse filtering).
Contents
 Reference
Functions

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.

void dip::
WienerDeconvolution (dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 1e4, dip::StringSet const& options = {S::PAD})  Wiener deconvolution using an estimate of noisetosignal ratio.

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})  TikhonovMiller deconvolution.

void dip::
IterativeConstrainedTikhonovMiller (dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 0.1, dip::dfloat tolerance = 1e6, dip::uint maxIterations = 30, dip::dfloat stepSize = 0.0, dip::StringSet const& options = {S::PAD})  Iterative Constrained TikhonovMiller (ICTM) deconvolution.

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})  RichardsonLucy (RL) deconvolution, also sometimes called the expectation maximization (EM) method.

void dip::
FastIterativeShrinkageThresholding (dip::Image const& in, dip::Image const& psf, dip::Image& out, dip::dfloat regularization = 0.1, dip::dfloat tolerance = 1e6, dip::uint maxIterations = 30, dip::uint nScales = 3, dip::StringSet const& options = {S::PAD})  Fast Iterative ShrinkageThresholding (FISTA) deconvolution.
Function documentation
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 ofin
. 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 noisetosignal
ratio instead of the signal and noise power spectra. Note that can be rewritten as
where is the noisetosignal 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 zeropadded 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 realvalued and scalar, except if the OFT is given instead of the PSF, in which
case psf
could be complexvalued.
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 frequencydomain convolution. "pad"
cannot be combined with "OTF"
.
void
dip::WienerDeconvolution (dip::Image const& in,
dip::Image const& psf,
dip::Image& out,
dip::dfloat regularization = 1e4,
dip::StringSet const& options = {S::PAD})
Wiener deconvolution using an estimate of noisetosignal 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
noisetosignal 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})
TikhonovMiller deconvolution.
Assuming some original image , a known convolution kernel (given by psf
), a noise realization ,
and an observed image (given by in
), the TikhonovMiller 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 TikhonovMiller 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 zeropadded 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 realvalued and scalar, except if the OFT is given instead of the PSF, in which
case psf
could be complexvalued.
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 frequencydomain 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 = 1e6,
dip::uint maxIterations = 30,
dip::dfloat stepSize = 0.0,
dip::StringSet const& options = {S::PAD})
Iterative Constrained TikhonovMiller (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 zeropadded 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 realvalued and scalar, except if the OFT is given instead of the PSF, in which
case psf
could be complexvalued.
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 frequencydomain 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})
RichardsonLucy (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, nonregularized RichardsonLucy 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 staircase effect, as it penalizes slow transitions but allows sharp jumps.
psf
is given in the spatial domain, and will be zeropadded 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 realvalued and scalar, except if the OFT is given instead of the PSF, in which
case psf
could be complexvalued.
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 frequencydomain 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 = 1e6,
dip::uint maxIterations = 30,
dip::uint nScales = 3,
dip::StringSet const& options = {S::PAD})
Fast Iterative ShrinkageThresholding (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 zeropadded 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 realvalued and scalar, except if the OFT is given instead of the PSF, in which
case psf
could be complexvalued.
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 frequencydomain convolution. "pad"
cannot be combined with "OTF"
.