MP/BME 574 Lecture 16: Deconvolution: The Inverse Filter, Iterative, and Expectation Maximization

Learning Objectives:

  • 2D deconvolution concepts
  • Introduction to the inverse filter and iterative solutions to the inverse problem

Assignment:

  1. Read Chapter 5of Gonzalez, Woods, and Eddins, titled “Image Restoration.”
  2. Read “An iterative technique for the rectification of observed distributions,” by Lucy et al. available on the website in pdf.
  1. The “Forward” and “Inverse” problem
  1. Forward problem
  2. Creating an object and simulating images of this object so as to test an imaging system’s performance without actually building it.
  3. For example, say we have a new reconstruction (or image compression method, etc...) method we wish to test out.
  4. Create or choose a phantom
  5. Apply the experimental reconstruction to noiseless and noisy data
  6. Compare with “conventional method” using some accepted image metric of performance
  7. Inverse problem
  8. Given noisy real-world data, we desire to reconstruct the object that produced the data
  9. For example,
  10. Establish boundaries on the problem, i.e. is the data projection data? Is it sufficiently sampled? Do we know the PSF of the system used to acquire the data?
  11. Choose the appropriate mathematical reconstruction method
  12. May be entirely deterministic but may have to be modified if data are very noisy or only partially known
  13. Restore the image and compare with “conventional method” using some accepted image metric of performance.
  1. Deconvolution is an example of the Inverse problem
  1. Typically used in imaging applications involving astronomy and microscopy (optical).
  2. Application in modeling of biological processes such as perfusion.
  3. Has all of the elements of the inverse problem used in parametricimage reconstruction.
  1. Inverse Filter
  2. Consider the 2D convolution, , and it Fourier transform

. This is a filter process such that the original data are filtered by the point response function of the system. It is desired that we obtain the deconvolved image, , where we define the inverse filter .

  1. The inverse filter is unstable for near-zero values of the modulation transfer function, MTF = , andwhen contains significant noise, particularly at higher spatial frequencies. The resulting noise amplification is unacceptable in most practical imaging applications.

  1. Stabilizing the inverse filter can be achieved by applying an empirical noise threshold, such that the inverse filter values are de-emphasized for lower values of the MTF:

(1)

This will mitigate noise amplification at the expense of degraded performance for the deconvolution at higher spatial frequencies. Given that the PSF is likely a low pass filter process; this is also less than desirable.

  1. If the noise process is known, then Weiner deconvolution becomes an option:

This filter can be shown to minimize the mean square error (MSE) between the true object and the estimate:

  1. Iterative approaches are a robust way to both stabilize the deconvolution process and improve performance at higher frequencies. The iterative deconvolution algorithm can be generally represented mathematically below:

(2)

where the ‘^’ symbol indicates an estimate of the function and the subscript indicates the iteration number. Note that the initial estimate,, for the image is simply a scaled version of the measured image, . The choice for the initial guess is not as critical in the present implementation where the point response function, , is known.

The monotonic convergence for this iterative approach to the inverse filter is assured assuming is continuous and appropriately normalized. The term  controls the convergence rate and can be used to slow the rate of convergence to prevent large departures from the optimal estimate of the image. The convergence to the inverse filter can be shown by solving Equations 2 above recursively. Solving for after Fourier transformation of Eqs. 2:

.

For k iterations then,

Therefore as , then for .

Note that for large the estimate converges more rapidly to the inverse filter. For noisy data the faster rate of convergence can lead to oscillations and errors. For small , the convergence rate is slowed with increased stability at the expense of increased computation time.

For so-called “Blind” deconvolution methods where both and must be estimated, the optimal solution can be highly dependent on the initial guess. This will be addressed later in this lecture.

  1. Blind Deconvolution
  1. In the blind deconvolution problem, both f and h are need to be estimated simultaneously. If nothing is known about either function this is not possible.
  1. Therefore a strategy is to combine the robust convergence properties of iterative techniques with apriori assumptions about the form of the data including statistical models of uncertainty in the measurements.
  2. General assumptions about the physical boundary conditions and uncertainty in the data
  3. e.g. Non-negativity and compact support.
  4. Statistical models of variation in the measured data:
  5. e.g. Poisson or Gaussian distributed.
  6. This leads to estimates of expected values for the measured data for Maximum Likelihood (ML) optimization.
  7. Physical parameter constrain the solution, while the ML approach provides a criterion for evaluating convergence.
  8. Maximum Likelihood:
  9. Consider a data estimation problem in which the uncertainty in the measured data are assumed to be governed by a Gaussian probability density function (pdf):

It is acceptable to evaluate the log-likelihood since log(Pr) is a monotonically increasing function. Therefore, we maximize the total probability by maximizing:

.

Therefore, the log likelihood of the measured data is maximized for a model in which, is minimized.

  1. Now consider the iterative ML approach to the blind deconvolution problem:

then at each k,

and

is minimized and used to optimize the convergence. The conditions for convergence are similar to the iterative procedure when is known except that the convergence to the inverse filter is no longer guaranteed and is sensitive to noise, choice of , and the initial guess, .

  1. Statistical model of the convolution process:
  2. Derives from the ML concept applied to a statistical model of the convolution process. In this approach, x, is a random variable, and represents an estimate of a probability density function, , that models the missing or unknown data.
  3. Intuitively is the superposition of multiple random processes used as probes (i.e. individual photons, or molecules of dye) use to measure the response of the system. The physical system must adhere to mass balance and, for finite counting statistics, non-negativity.
  4. For example, consider:

, where is our measured image data, is the desired corrected image, and is a conditional probability density function kernel that relates the expected value of the data to the measured data, e.g. assuming the photon counts in our x-ray image are approximately Gaussian distributed about there expected value  For this example, then

becomes our familiar convolution process.

  1. Expectation maximization
  2. Expectation in the sense that we use a statistical model of the data measurement process to estimate the missing data
  3. Maximization in the maximum likelihood sense, where iteration is used within appropriate physical constraints to maximize the likelihood of the probability model for the image.
  4. Consider an “inverse” conditional probability function given by Beyes’ Theorem:

.

Then it is possible to estimate the value of the inverse probability function iteratively from current guesses of at the kth iteration of the measured image and deconvolved image respectively.

Our iterative estimate of the inverse filter is then:

,

where,

, and

Putting this all together starting with the last result and substituting, then the iterative estimate of the image is:

.

This is guaranteed to converge if the x,  are non-negative and the respective areas of and are conserved. This is because will approach and the will then approach in these circumstances. Note that the model has remained general. As long as the model follows the requirements of a probability density function (pdf), its form can depend on the desired application. This is not to say that the algorithm is guaranteed to converge to the global maximum likelihood result although in practice the algorithm is very robust in applications where there is sufficient SNR.