Enhancement and Restoration of Microscopic Images Corrupted with Poisson’s Noise using a Nonlinear Partial Differential Equation-based Filter

An inherent characteristic of the many imaging modalities such as fluorescence microscopy and other microscopic modalities is the presence of intrinsic Poisson noise that may lead to degradation of the captured image during its formation. A nonlinear complex diffusion-based filter adapted to Poisson noise is proposed in this paper to restore and enhance the degraded microscopic images captured by imaging devices having photon limited light detectors. The proposed filter is based on a maximum a posterior approach to the image reconstruction problem. The formulation of the filtering problem as maximisation of a posterior is useful because it allows one to incorporate the Poisson likelihood term as a data attachment which can be added to an image prior model. Here, the Gibb’s image prior model-based on energy functional defined in terms of gradient norm of the image is used. The performance of the proposed scheme has been compared with other standard techniques available in literature such as Wiener filter, regularised filter, Lucy-Richardson filter and another proposed nonlinear anisotropic diffusion-based filter in terms of mean square error, peak signal-to-noise ratio, correlation parameter and mean structure similarity index map. The results shows that the proposed complex diffusion-based filter adapted to Poisson noise performs better in comparison to other filters and is better choice for reduction of intrinsic Poisson noise from the digital microscopic images and it is also well capable of preserving edges and radiometric information such as luminance and contrast of the restored image.

For Poisson’ noise corrupted digital images such as in fluorescence microscopy1-4,9, the sample to be studied is itself the light source and this technique is used to study specimens, which can be made to fluoresce. The fluorescence microscope is based on the phenomenon that certain material emits energy detectable as visible light when irradiated with the light of a specific wavelength. The sample can either be fluorescing in its natural form like chlorophyll and some minerals, or treated with fluorescing chemicals. The basic task of the fluorescence microscope is to let excitation light radiate the specimen and then sort out the much weaker emitted light to make up the image. First, the microscope has a filter that only lets through radiation with the desired wavelength that matches the given fluorescing material. The radiation collides with the atoms in the given specimen and electrons are excited to a higher energy level. When they relax to a lower level, they emit light. To become visible, the emitted light is separated from the brighter excitation light in a second filter. Here, the fact that the emitted light is of lower energy and has a longer wavelength is used. The fluorescing areas can be observed in the microscope and shine out against a dark background with high contrast. Fluorescence microscopy is a rapid expanding technique, both in the medical and biological sciences. The technique has made it possible to identify cells and cellular components with a high degree of specificity. For example, certain antibodies and disease conditions or impurities in inorganic material can be studied with the fluorescence microscopy. In fluorescence microscopy1-3, the captured image may be noisy due to various reasons that may be extrinsic or intrinsic in nature. The extrinsic noise may be induced due to lens miss-focus, environmental factors, instrumental error, light detectors and sensors. The extrinsic noise can be caused by various sources such as dark current, electronic noise, detector readout which are Gaussian distributed, and quantization noise which is uniformly distributed4. The intrinsic noise may be induced by the detection of photons. A light detector is said to be photon-limited if the extrinsic noise is negligibly small compared to the amount of intrinsic noise induced by the detection of photons. Measurements show that scientific CCD cameras and PMTs can be regarded as photon limited2,4.

In the fluorescence microscopy, fluorescent molecules in the illuminated object are excited by an incident light of wavelength λ1 and the excited molecules emit light of wavelength λ2 which is collected by the microscope forming a fluorescent image. The difference of light wavelengths Δ λ=λ1- λ2 between λ1 and λ2 is called the Stokes shift of the fluorescent molecule9. The light may be considered as a series of particles called photons where each photon carries a certain amount of energy E = hc/λ, where c is the speed of light and h is Planck’s constant. Photon production by any light source is a statistical process governed by the laws of quantum physics where the source emits photons at random time intervals and the number of photons in a fixed observation interval results in a number that follows Poisson statistics. Each observation measures a number p with a probability given by the Poisson distribution2,9.

$P\left(p|\rho T\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{\left(\rho T\right)}^{p}{e}^{-pT}}{p!}$ (1)

where T is the observation interval or exposure time and ρ is the photon flux. The average of a large number of observations approximate the expected photon production ρT. This photon detection induced Poisson noise is sometimes referred to as intrinsic noise, and is unavoidable when acquiring an image. A method is proposed to deal with intrinsic Poisson noise for the photon limited light detectors.

To construct a model for the image formation of a digital image such as fluorescence microscope2,5,10 both the wave description and quantum nature of light can be used where the incoherent nature of the fluorescence light allows one to model the fluorescence microscope as a linear translation invariant system. The image is in that case a convolution of the object with the point spread function of the fluorescence microscope. The detection of photons in a finite time interval distorts the observed image with Poisson noise. Extrinsic noise sources can further hamper the image. Finally it is common in fluorescence microscopy to measure a non-zero background level arising from auto-fluorescence, inadequate removal of fluorescent staining material, glare and reflections, and/or offset levels associated with the gain of the detector or other electronic sources. Therefore, the image formation of a fluorescence microscope can be modelled as follows

${u}_{0}\left(x,y\right)\text{\hspace{0.17em}}=\eta \text{\hspace{0.17em}}\left[\left(h\left(x,y\right)\otimes u\left(x,y\right)+\beta \left(x,y\right)\right]$ (2)

where u0(x,y) is the observed or recorded fluorescence image, u(x,y) is the original true image without distortion, β(x,y) a background signal, h(x,y) is point spread function (PSF) of the microscope, and η is a general noise distortion function which in the case of scientific grade light detectors is dominated by Poisson noise. For images with a relatively high signal-to-noise ratio and small dynamic range the Poisson process can be approximated by additive Gaussian noise. This approximation of an inhomogeneous Poisson process will in general lead to a different sigma (standard deviation) of the Gaussian distribution for each recorded data point i.e., pixels of image because the variance of a Poisson process is equal to its intensity. However, for images with a relatively high signal-to-noise ratio and small dynamic range Poisson process can be approximated with constant sigma, making the Gaussian noise independent from the pixel intensity. Under these conditions, the image observation model1 given by Eqn. (2) reads

${u}_{0}\left(x,y\right)\text{\hspace{0.17em}}=\left(h\left(x,y\right)\otimes u\left(x,y\right)\right)+\beta \left(x,y\right)$ (3)

$\text{which can be written as}\text{\hspace{0.17em}}{u}_{0}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}Hu+\beta$ (4)

The various restorations methods1,5-8,10-13 available in literature for the model given by Eqn. (4) can be categorised in three broad categories. The first being the linear methods; the second being the iterative nonlinear filters; and the third being the maximum likelihood (ML) estimation based on expectation-maximization (EM) algorithm5,12. The approaches for designing linear filters are based on the concepts of Minimum mean square error (MMSE) and least square restoration. Examples of MMSE and least square filters include the popular Wiener filter5,6. The Wiener filter de-convolves an image using Wiener filter returning de-blurred image. The assumption is that the input image was created by convolving a true image with a point-spread function PSF and possibly by adding noise. The Wiener filter is optimal in a sense of least mean square error between the estimated and the true images, and utilizes the correlation matrices of image and noise. In the absence of noise, the Weiner filter reduces to the ideal inverse filter. Examples of iterative nonlinear filters and ML estimation based on EM includes regularised filter5 and iterative constrained Lucy-Richardson filter5,7,8. The regularised filter5 is a constrained optimum in a sense of least square error between the estimated and the true images under the requirement of preserving image smoothness. The accelerated, damped, Lucy-Richardson filter5,6 maximizes the likelihood that the resulting image, when convolved with the PSF, is an instance of the blurred image, assuming Poisson noise statistics. The Lucy-Richardson filter can be effective when one knows the PSF but know little about the additive noise in the image. The limitations of these filters are: the Wiener1,2 filter is linear and convolution filter and its linear nature makes it incapable of restoring frequencies for which the PSF has a zero response and in general the PSF of a 3-D fluorescence microscope has large regions with zero response known as the missing cones2. Further, linear filters can’t restrict the domain in which the solution should be found and this property is a major drawback since the intensity of an object represents light energy, which is non-negative. Also, linear filters are very sensitive to errors. Further, all the other methods explained as above are based on approximation of the Poisson noise with Gaussian noise. Further all above filters perform better for additive Gaussian noise but do not perform for Poisson noise in fluorescence microscopic images. They are also not capable of preserving edges and radiometric information such as luminance and contrast in the restored image. Therefore, to recover the fluorescence microscopic images an effective method adapted to Poisson noise is required.

In recent years, some authors14-19 have proposed partial differential equation (PDE) based filters for the restoration of digital images corrupted with additive noise. These filters are capable of restoring the images and are also well capable of preserving the edges. The PDE based filters can be modified to remove the intrinsic Poisson noise from fluorescence microscopic images in an effective manner. The PDE based filters can be derived in variational framework where the energy of a noisy image may be defined in terms of gradient norm of the image15,17,18 and further, this energy functional is minimized to minimize the variations or noise in an image using Euler-Lagrange minimisation technique combined with gradient descent approach. After minimization, PDE is obtained which acts as the nonlinear filter whose initial condition is the noisy image. The processed or restored image is obtained after certain iterations of the obtained PDE till it converges to the solution.

A nonlinear PDE based filter adapted to Poisson noise is proposed in this paper to remove photon-limited intrinsic noise which has Poisson distribution from the digital microscopic images. Further, the performance of the proposed scheme is evaluated both quantitatively and qualitatively and results are compared against some well known filters in literature such as Wiener filter5,6, Regularised filter5,6, Lucy-Richardson filter5,7,8 and another proposed nonlinear complex diffusion-based filter.

2.1 Anisotropic Diffusion-based Model

To overcome the problems associated with the restoration methods a nonlinear anisotropic diffusion-based filter adapted to Poisson noise is proposed to restore the digital images corrupted with Poisson noise. The proposed filter is based on a maximum a posterior (MAP) approach to the image reconstruction problem. Given an initial noisy image u0, then we reconstruct the filtered image u that maximizes the log-posterior probability

$\mathrm{log}\left(p\left(u/{u}_{0}\right)\right)\mathrm{\alpha }\mathrm{log}\left(p\left({u}_{0}/u\right)\right)+\mathrm{log}\text{\hspace{0.17em}}p\left(u\right)$ (5)

where p(u/u0) is the likelihood term of noise model and p(u) is the prior. The formulation of the filtering problem as maximization of a posterior is useful because it allows one to incorporate the Poisson likelihood term as a data attachment which can be added to an image prior model. In this work, the Gibb’s image prior model based on energy functional defined in terms of gradient norm18 of the image is used.

In fluorescence microscopy and other imaging applications such as medical imaging and astronomical imaging, the elements of noisy image u0 belongs to noisy photon counts which follows Poisson distribution. Therefore, the statistical model that models the error in these systems reads u0 = Poisson(Au+b)Poisson(X). Where Poisson(Au+b) is an independent and identically distributed (iid) Poisson random vector with Poisson parameter vector X.

The probability distribution function of the Poisson noise20 is given as

$p\left(u/{u}_{0}\right)=\frac{{e}^{-u}.{u}^{{u}_{0}}}{\angle {u}_{0}}$ (6)

Given image data ouarising from model u0=Poisson(Au+b) Poisson(X), whose pdf is described by Eqn (6), the maximum likelihood estimate of u is obtained by maximising p(u/ u0)with respect to subject to u ≥ 0. Alternatively, we can calculate the maximum likelihood estimate of u by minimising the negative log likelihood of Poisson pdf given by ${u}_{\text{ML}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\underset{u\ge 0}{\mathrm{arg}\mathrm{min}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left\{-\text{In}\text{\hspace{0.17em}}\text{\hspace{0.17em}}p\left(u/{u}_{0}\right)\right\}=\text{\hspace{0.17em}}\underset{u\ge 0}{\mathrm{arg}\mathrm{min}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left\{u-{u}_{0}\text{\hspace{0.17em}}\text{In}u\right\}$ to restore the image data u.

The log likelihood of Poisson pdf reads

$\mathrm{log}\left(p\left(u/{u}_{0}\right)\right)=\mathrm{log}\left[\frac{{e}^{-u}.{u}^{{u}_{0}}}{\angle {u}_{0}}\right]=\text{\hspace{0.17em}}\left(-u+{u}_{0}\mathrm{log}u-\mathrm{log}{u}_{0}\right)$ (7)

For obtaining the maximum likelihood of u, the derivative of log likelihood of Poisson pdf wrt u reads

$\frac{\partial \mathrm{log}\left(p\left(u/{u}_{0}\right)\right)}{\partial u}=\frac{\left({u}_{0}-u\right)}{u}$ (8)

$\text{Let},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}B=\frac{\left({u}_{0}\text{-}u\right)}{u}$ (9)

The term B acts as the data attachment term or data likelihood term i.e. B is the maximum likelihood estimate of u assuming Poisson noise.

The Gibbs image prior model is considered as the image prior model for MAP based restoration method. The Gibbs prior model is based on the energy functional, which is defined in terms of the gradient norm of the image, related to anisotropic diffusion-based PDE18 for additive noise removal. The Gibb’s prior model reads

$p\left(u\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{\text{z}}\mathrm{exp}\left(-E\left(u\right)\right)$ (10)

Where energy functional18 is defined as

$E\left(u\right)=\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}\left(\text{λ}\underset{\mathrm{\Omega }}{\int }||\nabla u|{|}^{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}d\mathrm{\Omega }\right)$ (11)

Using the Euler-Lagrange minimization technique combined with gradient descent approach the Eqn. (11) leads to following PDE for anisotropic regularisation of the image data18

$\frac{\partial u}{\partial t}=\text{λ}\nabla .\left(c\left(||\nabla u||\right)\nabla u\right)$ (12)

By adding the data likelihood term B given by Eqn. (9) with the regularisation term given by Eqn. (11), the update equation for the filtered image reads

$\frac{\partial u}{\partial t}=B+\text{λ}\nabla .\left(c\left(||\nabla u||\right)\nabla u\right)$ (13)

Where the diffusion coefficient is defined as18,

$c\left(||\nabla u||\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{1+{\left(\frac{\nabla u}{k}\right)}^{2}}$ (14)

$\frac{\partial u}{\partial t}=\frac{\left({u}_{0}-u\right)}{u}+\lambda \nabla .\left(c\left(||\nabla u||\right)\nabla u\right)$ (15a)

with initial condition being the observed noisy image,

$\text{i.e}\text{.,}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}u\left(t=0\right)={u}_{0}$ (15b)

The first term, which is first derivative of log likelihood of Poisson probability distribution function (pdf) with respect to estimated image, acts as the data attachment term and measures the dissimilarities at a pixel between observed image and its estimated value obtained during filtering process there by making the whole filtering process adapted to noise. The second term is responsible for regularisation and smoothing of the image data by minimising the variance of pixels.

The image restoration models given by Eqns (15a) and (15b) derived using MAP approach can also be alternatively derived in minimisation framework using calculus of variations and Euler-Lagrange minimisation. Therefore, the anisotropic diffusion penalized Poisson maximum likelihood estimation can also be obtained by minimising ${u}_{\lambda }\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\underset{u\ge 0}{\mathrm{arg}\text{\hspace{0.17em}}\mathrm{min}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left\{{T}_{\lambda }\left(u\right)=\left(u-{u}_{o}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{ln}\text{\hspace{0.17em}}\text{\hspace{0.17em}}u\right)\text{\hspace{0.17em}}+\frac{\lambda }{2}\text{\hspace{0.17em}}|\nabla {\text{u|}}^{2}\right\}$ resulting in Eqns (15a) and (15b).

2.2 Proposed Complex Diffusion-based Model

Another filter-based on the concepts of complex diffusion-based processes19 is proposed and its efficacy is also examined. The major advantage of complex diffusion-based filter over anisotropic diffusion filter is that more generalized and efficient19. To derive the filter based on complex diffusion-based processes, if in second term of the anisotropic diffusion-based model given by Eqn. (15a), real time factor t is replaced by the complex time factor it and the diffusion coefficient c(|∆u|) by c(Im(u)), where Im(u) is the imaginary part of the image as computed in paper19, then it leads to following proposed nonlinear complex diffusion-based model19 adapted to Poisson noise:

$\frac{\partial I}{\partial t}=\frac{\left({u}_{0}-u\right)}{u}+\mathrm{\lambda }\cdot \text{\hspace{0.17em}}\text{\hspace{0.17em}}div\left(c\left(\mathrm{Im}\left(u\right)\right)\nabla u\right)$ (16a)

(16b)

The diffusion coefficient c(Im(u)) is defined as follows19

$c\left(\mathrm{Im}\left(I\right)\right)=\frac{{e}^{i\mathrm{\theta }}}{1+{\left(\frac{\mathrm{Im}\left(I\right)}{k\mathrm{\theta }}\right)}^{2}}$ (17)

where k is edge threshold parameter and the value of k ranges from 1 to 1.5 for digital images19. For nonlinear complex diffusion process defined by Eqns (16) and (17) the evolution of real part of the image is controlled by the linear forward diffusion, whereas evolution of imaginary part of the image is controlled by both the real and imaginary equations. A qualitative property of edge detection, i.e., second smoothed derivative is described by the imaginary part of the image for small value of θ, whereas real values depict the properties of ordinary Gaussian scale-space. For large values of θ, the imaginary part feeds back in to the real part creating the wave like ringing effect which is an undesirable property. Here, for experimentation purposes the value of θ is chosen to be $\frac{\ne }{30}$ .

2.3 Discretisation of the Proposed Model

For digital implementations, the Eqn. (15) can be discretised using finite differences schemes21. The discretised form of the anisotropic diffusion-based model, given by Eqn. (15), reads

${u}^{n+1}={u}^{n}+\mathrm{\Delta }t.\left[\frac{\left({u}_{0}-{u}^{n}\right)}{{u}^{n}}+\mathrm{\lambda }\nabla \text{.(}c\left(‖\nabla {u}^{n}‖\right)\nabla {u}^{n}\right)\right]$ (18a)

$u\left(t=0\right)={u}_{0}$ (18b)

And the discretised form of the proposed complex diffusion-based model reads

${u}^{n+1}={u}^{n}+\mathrm{\Delta }t.\left[\frac{\left({u}_{0}-{u}^{n}\right)}{{u}^{n}}+\mathrm{\lambda }\nabla \text{.(}c\left(\mathrm{Im}\left({u}^{n}\right)\nabla {u}^{n}\right)\right]$ (19a)

$u\left(t=0\right)={u}_{0}$ (19b)

For the numerical scheme, given by Eqns (18) and (19) to be stable, the von Neumann analysis21, shows that we require $\frac{\varnothing t}{{\left(\varnothing x\right)}^{2}}<1/4\text{\hspace{0.17em}}$. If the grid size is set to Δx =1, then $\mathrm{\Delta }t<1/4\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{i}\text{.e}\text{.}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{\Delta }t<0.25.$ i.e ∆t < 0.25. Therefore, the value of ∆t is set to 0.25 for stability of Eqn. (18). In the similar fashion, the nonlinear complex diffusion model given by Eqn. (16) can be discretised using finite difference scheme. Out of the two proposed schemes, the complex diffusion-based model adapted to Poisson noise is performing better to all other schemes in consideration which is validated by the results.

2.4 Estimation of Regularisation Parameter λ

The regularisation parameter λ, used in Eqns (15), (16) and (19) are computed dynamically from the information present in the fluorescence microscopic image. Initially it was set to 0 and in successive iterations of the proposed PDE based model the specific value of λ is calculated. The regularisation parameter λ is responsible for making a balance between the two terms which are data attachment term and regularisation term in the proposed model. The regularisation parameter is defined in terms of the inverse of average signal-to-noise ratio (SNR) computed from mean and standard deviation of the image in specific iterations. The average SNR is defined as follows:

$Avg\text{\hspace{0.17em}}.\text{\hspace{0.17em}}SNR\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{mean\left(I\left(t\right)\right)}{std.\text{\hspace{0.17em}}deviation\left(I\left(t\right)\right)}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\mathrm{\mu }}{\sqrt{{\mathrm{\sigma }}_{n}^{2}}}$ (20)

$\mathrm{\lambda }\text{\hspace{0.17em}}\text{=}\text{\hspace{0.17em}}\frac{1}{Avg\text{\hspace{0.17em}}.\text{\hspace{0.17em}}SNR}$ (21)

The lower value of SNR denotes that more noise is present in image whereas the higher value of SNR denotes that the image is of better quality and contains less noise in comparison to its previous more noisy version. Therefore, for lower SNR value, the value of λ increases allowing more weight to the second term in the proposed model which is responsible for regularisation or smoothing of data. For higher SNR value, the value of λ decreases allowing less weight to the second term in the proposed model which is responsible for regularisation or smoothing of data.

The performance of the proposed anisotropic diffusion- based scheme adapted to Poisson noise and other noise reduction schemes available in literature such as Wiener filter5,6, Lucy-Richardson filter5,7,8, Regularised filter5,6 and another proposed complex diffusion-based filter have been evaluated both qualitatively and quantitatively and comparative study of their performances is presented for four sample microscopic images.

3.1 Performance Metrics

The metrics for comparing the performance of various schemes in considerations for digital images corrupted with Poisson noise are defined as follows:

Mean square error5:

$MSE\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{m\text{\hspace{0.17em}}×\text{\hspace{0.17em}}n}\text{\hspace{0.17em}}\sum _{i\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{m}\sum _{j\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{n}{\left[{I}^{\prime }\left(i\text{\hspace{0.17em}},\text{\hspace{0.17em}}j\right)\text{\hspace{0.17em}}-\text{\hspace{0.17em}}I\left(i,\text{\hspace{0.17em}}j\right)\right]}^{2}$ (22)

where is I′ the original image without noise, I′ is the filtered noise reduced image, m x n is the size of the image and i = 1……..m, j = 1……..n.

Peak signal-to-noise ratio5:

$PSNR\text{\hspace{0.17em}}=\text{\hspace{0.17em}}20\text{\hspace{0.17em}}{\mathrm{log}}_{10}\left[\frac{255}{RMSE}\right]$ (23)

Here RMSE is the root mean square error. For optimal performance, measured values of MSE should be small and that of peak signal-to-noise ratio (PSNR) should be large.

Correlation parameter22: Correlation parameter (CP) is a qualitative measure for edge preservation. To evaluate the performance of the edge preservation or sharpness, the correlation parameter is defined as follows:

$CP\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\sum _{i\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{m}\sum _{j\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{n}\left(\mathrm{\Delta }I\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\mathrm{\Delta }\overline{I}\right)\text{\hspace{0.17em}}×\text{\hspace{0.17em}}\left(\mathrm{\Delta }\stackrel{^}{I}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\mathrm{\Delta }\overline{\stackrel{^}{I}}\right)}{\sqrt{\sum _{i\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{m}\sum _{j\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{n}{\left(\mathrm{\Delta }I\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\mathrm{\Delta }\overline{I}\right)}^{2}×}\text{\hspace{0.17em}}\sum _{i\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{m}\sum _{j\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{n}{\left(\mathrm{\Delta }\stackrel{^}{I}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\mathrm{\Delta }\overline{\stackrel{^}{I}}\right)}^{2}}$ (24)

where ΔΙ and ΔÎ are high pass filtered versions of original imageand I filtered image Î obtained via a 3 × 3 pixel standard approximation of the Laplacian operator. The ΔĪ and $\mathrm{\Delta }\overline{\stackrel{^}{I}}$ are the mean values of I and Î respectively. The correlation parameter should be closer to unity for an optimal effect of edge preservation.

Structure similarity index map23: Structure similarity index map (SSIM) is used to compare luminance, contrast and structure of two different images. It can be treated as a similarity measure of two different images. This similarity measure is a function of luminance, contrast and structure. The SSIM of two images X and Y can be calculated as:

$SSIM\left(X,\text{\hspace{0.17em}}Y\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\left(2{\mathrm{\mu }}_{x}{\mathrm{\mu }}_{y}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{C}_{1}\right)\text{\hspace{0.17em}}×\text{\hspace{0.17em}}\left(2{\text{σ}}_{xy}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{C}_{2}\right)}{\left({\mathrm{\mu }}_{x}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\mathrm{\mu }}_{y}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{C}_{1}\right)\text{\hspace{0.17em}}×\text{\hspace{0.17em}}\left({\mathrm{\sigma }}_{x}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\mathrm{\sigma }}_{y}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{C}_{2}\right)}$ (25)

where ${\propto }_{i}\left(i=X\text{\hspace{0.17em}}\mathrm{or}\text{\hspace{0.17em}}Y\right)$ is the mean intensity, σi (i = X or Y) is the standard deviation, σxyxy and Ci (i=1 or 2) is the constant to avoid instability when ${\mu }_{x}^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\mu }_{y}^{2}$ is very close to zero and is defined as Ci = (kiL)2 in which ki << 1 and L is the dynamic range of pixel values e.g. L = 255 for 8-bit gray scale image. To have an overall quality measurement of the entire image, mean SSIM is defined as

$MSSIM\left(X,\text{\hspace{0.17em}}Y\right)\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{1}{mn}\text{\hspace{0.17em}}\sum _{i\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{m}\sum _{j\text{\hspace{0.17em}}=\text{\hspace{0.17em}}1}^{n}SSIM\left({X}_{ij},\text{\hspace{0.17em}}{Y}_{ij}\right)$ (26)

The MSSIM value should be closer to unity for optimal measure of similarity.

3.2 Results and Discussions

The performances of all schemes in consideration have been compared in terms of MSE, PSNR, CP and MSSIM for the all four sample microscopic images. The Poisson noise was introduced artificially in the test images. It has been tested through experimentation that after 50 iterations the proposed diffusion-based PDEs in consideration converge to the desired level of solution i.e., produces acceptable quality of de-noised images. After 50 iterations the PSNR values of the proposed PDE based schemes start decreasing. Therefore total numbers of iterations were fixed to 50 for the proposed schemes. The value of was set to 0.25 for stability reasons. The initial condition of the proposed PDE based filter is the noisy image and the final solution is the de-noised image. The initial value of the dynamic regularisation parameter λ(0) is set to zero and the value of λ was calculated in each iteration of the proposed PDE-based model by Eqns (19) and (20). All schemes have been implemented using MATLAB 7.0 software. For implementing the proposed method, the Eqn. (19) was used for evolving the solution. The diffusion coefficient C(Im(u))is defined by Eqn. (17). The value of threshold parameter k used in Eqn. (17) lies in between 1 and 1.5 for de-nosing the image corrupted by Gaussian noise19. For images corrupted with Poisson noise, the value of edge threshold parameter k for the nonlinear complex diffusion-based scheme was set to 0.0001 and is the optimal choice as tested by experimentations. For other values of k, the performance decreases. For implementing the anisotropic diffusion-based model Eqn. (18) was used. The diffusion coefficient C(║∇I║) defined by Eqn. (14) was used, where the value of threshold parameter k used in Eqn. (14) was set to 20 as proposed18. It has been tested through experimentation that, for other values of k the performance decreases. For implementing nonlinear complex diffusion-based technique, the discretised versions of Eqn. (16) was used.

The results for four test cases are shown in this paper though the performances of all schemes in considerations were evaluated for several other test images and the performance trend almost remained same. Tables 1-4 lists the performance results for the four sample microscopic images cheek3oilmicroscopic.jpg of size 480 × 640, microscopicobject.jpg of size 80 × 80, cos12.jpg, 504 × 700, and bpae3.jpg, 504 × 700 respectively for the proposed scheme and other schemes in consideration in terms of MSE, PSNR, CP and MSSIM. Figures 2-5 show comparison of visual results of various filters for the four sample microscopic images in consideration. From Tables 1-4, it can be observed that the proposed complex diffusion-based model for Poisson noise reduction is associated with minimum MSE, maximum PSNR, CP and MSSIM values in comparison to other schemes for all four sample microscopic images in consideration. The maximum values of CP and MSSIM associated with the proposed method, which are very close to one, indicate that the proposed scheme is well capable of preserving edges and structures of microscopic images in addition to effective reduction of intrinsic Poisson noise.

Therefore, from the results obtained it can be concluded that the proposed complex diffusion-based model is a better choice for reduction of the intrinsic Poisson noise from microscopic images and it also preserves the edges and other radiometric information such as luminance and contrast of the image.

A nonlinear PDE based filter, i.e., nonlinear complex diffusion-based filter adapted to Poisson noise is proposed in this paper to restore and enhance the degraded microscopic images corrupted by intrinsic Poisson noise.

The proposed filter is based on a maximum a posterior (MAP) approach to the image reconstruction problem. The first term in the proposed model is the log likelihood term that makes the overall filtering procedure adapted to Poisson noise and the second term is responsible for the regularisation and smoothing of the image data. In the similar framework, another nonlinear anisotropic diffusion filter is proposed and its efficacy was also examined. For digital implementations, the proposed PDE based filters were discretised using finite difference scheme. The value of regularisation parameter λ was calculated dynamically in each iteration. Performance and efficacy of the proposed scheme were examined for several microscopic images and results and performance analysis for four sample microscopic images are presented here. The performance of the proposed scheme has also been compared with other standard techniques available in literature such as Wiener filter, regularised filter, Lucy-Richardson filter and another proposed nonlinear anisotropic diffusion-based filter in terms of mean square error (MSE), peak signal-to-noise ratio (PSNR), correlation parameter (CP) and mean structure similarity index map (MSSIM). From the results and performance analysis it can be concluded that the proposed complex diffusion- based filter adapted to Poisson noise performs better in comparison to other filters and is a better choice for reduction of intrinsic Poisson noise from the microscopic images and it is also well capable of preserving edges and radiometric information such as luminance and contrast of the restored image.

1. van Kempen, G.M.P.; van Vliet, L.J.; Verveer, P.J. & van der Voort, H.T.M. A quantitative comparison of image restoration methods for confocal microscopy. Journal Microscopy, 1997, 185, 354.

2. van Kempen, G.M.P. Image restoration in fluorescence microscopy. Delft University, The Netherlands, 1999. PhD Thesis.

3. Goodman, J.W. Statistical optics. Wiley-Interscience, New York, USA, 1985.

4. Castleman, K.R. Digital image processing. Prentice Hall Englewood Cliffs, 1996.

5. Jain, A. K. Fundamentals of digital image processing. PHI, India, 2006.

6. Gonzalez, R.C. & Wintz, P. Digital image processing. Academic Press, New York, USA, 1987.

7. Lucy, L.B. Image restoration of high photometric quality. In Proceedings of the Restoration of HST Images and Spectra, edited by R.J. Hanisch & R.L. White, STSci. 1974, pp. 79-85.

8. White, R.L. Image restoration using the damped Lucy- Richardson method. In Astronomical Data Analysis Software and System III, ASP Conference series, 61, 1994.

9. Lakowicz, J.R. Principles fluorescence spectroscopy. Plenum Press, New York,1983.

10. Snyder, D.L.; Hammoud, A.M. & White, R.L. Image recovery from data acquired with a charge-coupled-device camera. J. Opti. Soc. Ame. A – Opt. Ima. Sci., 1993, 10, 1014.

11. Agard, D.A.; Hiraoka, Y. & Sedat, J.W. Three dimensional microscopy: Image processing for high resolution subcellular imaging. In New Methods in Microscopy and Low Light Imaging, edited by J. E. Wampler, SPIE, USA, 1989, 1161, pp. 24-30.

12. Krishnamurthy, V.; Liu, Y.H.; Bhattacharyya, S.; Turner, J.N. & Holmes, T.J. Blind deconvolution of fluorescence micrographs by maximum likelihood estimation. Applied Optics, 1995, 34, 6633.

13. Tao, L. & Nicholson, C. The 3-D point spread functions of a microscope objective in image and object space. Journal Microscopy, 1995, 178, 267.

14. Caselles, V.; Morel, J. & Sapiro, G. Introduction to the special issue on partial differential Equations and geometry driven diffusions in image processing. IEEE Trans. Image Proce., 1998, 7, 269.

15. Witkin, A.P. Scale space filtering. In Proceedings of the International Joint Conference on Artificial Intelligence, Germany, 1983, pp. 1019-023.

16. B. ter Harr Romeny (ed.), Geometry driven diffusion in computer vision. Kluwer, Boston, MA, 1994.

17. Rudin, L.I.; Osher, S. & Fatemi, E. Non linear total variation based noise removal algorithm. In Proceedings of the ModeLIZ. Mat. Traittement d’Images, INRIA, 1992. pp. 149-79.

18. Perona, P. & Malik, J. Scale-space and edge detection using anisotropic diffusion. IEEE Trans. Pattern Anal. Mach. Intel., 1990, 12, 629.

19. Gilboa, G.; Sochen, N. & Zeevi, Y.Y. Image enhancement and denoising by complex diffusion processes. IEEE Trans. Pattern Anal. Mach. Intel., 2004, 25, 1020.

20. Papoulis. Probability, random variables, and stochastic processes. New York, 1991.

21. Press, W.H.; Teukolsky, S.A.; Vetterling, W.T. & Flannery, B.P. Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, USA, 1992.

22. Salinas, H.M. & Fernandez, D.C. Comparison of PDE- based nonlinear diffusion approaches for image enhancement and de-noising in optical coherence tomography. IEEE Trans. Med. Imag., 2007, 26, 761.

23. Wang, Z.; Bovik, A.C.; Sheikh, H.R. & Simon celli, EP. Image quality assessment: from error visibility to structural similarity. IEEE Trans. Image Proc., 2004, 13, 1.

 Dr Rajeev Srivastava received his ME (Computer Technology and Applications) and PhD (Computer Engg) both from University of Delhi, India, in 2005 and 2011 respectively. Currently working as an Associate Professor in the Department of Computer Engineering, Institute of Technology, Banaras Hindu University (IT-BHU), Varanasi, India. He has 26 research publications in national/ international conferences and journals and 04 book chapters to his credit. His research interests include image processing and vision, and algorithms. Dr J.R.P. Gupta received his BSc (Electrical Engineering) and PhD from the University of Bihar, India, in 1972 and 1983 respectively. Curently working as a Head, Department of Instrumentation and Control Engineering, University of Delhi. He is a senior member of IEEE and recipient of IETE K.S. Krishnan memorial award for the best system oriented paper. His research interest includes signal processing, power electronics and control systems. Dr Harish Parthasarthy received his BTech (Electrical Engineering) from Indian Institute of Technology (IIT), Kanpur, India, in 1990 and PhD from (IIT), Delhi, India, in 1994. He has also pursued Post-doctorate for a period of one year at Indian Institute of Astrophysics (IIAP), Bengaluru, India. Currently he is a Professor in the Department of Electronics and Communication Engineering, Netaji Subhas Institute of Technology (NSIT), New Delhi, India. His research interest include signal processing.