Restoration and enhancement of astronomical images using hybrid adaptive nonlinear complex diffusion-based filter

In this paper, a hybrid and adaptive nonlinear complex diffusion based technique is proposed for restoration and enhancement of astronomical images corrupted with additive noise due to diffractions limit, aberrations in telescope camera lens, atmospheric irregularities which are modelled by Gaussian functions. In addition to additive noise removal, the proposed filter is also capable of removing noisy stars and spike noises due to other numerous stars that may have surrounded an object in the astronomical image such as in nebula images. The proposed scheme is completely adaptive in nature in the sense that it estimates all the required filtering parameters from the observed image itself. The performance of the proposed hybrid scheme is compared with other image restoration techniques such as averaging filter, Gaussian filter, median filter, anisotropic diffusion based filter, local variance based adaptive anisotropic diffusion filter and also the adaptive and hybrid version of anisotropic diffusion filter similar to the proposed one in terms of average SNR and BSNR. The results obtained show the efficacy of the proposed scheme.

This paper deals with the restoration of astronomical images where blind restoration techniques are applied for the restoration of the same. In blind image restoration technique, the original image is estimated without the explicit knowledge of the underlying degradation process. This process is difficult since information about the original image or the underlying blurring process is not available in many practical applications such as in space exploration and astronomy. The astronomical images may be obtained with the help of instruments that may use any one of the astronomical imaging techniques such as radio astronomy, infrared astronomy, ultraviolet astronomy, optical astronomy, X-ray astronomy and Gamma-ray astronomy. Optical interferometry is also one of the techniques used in optical astronomy for acquiring images where the data correspond to Fourier Transform of the image. Irrespective of the methods and instruments used to acquire the astronomical images these images may be degraded due to various reasons for example when light from a star or another astronomical object enters the earth's atmosphere, atmospheric turbulence introduced by different temperature layers and different wind speeds interacting can distort and move the image in various ways. Images produced by any telescope larger than a few metres are blurred by these distortions. In addition to atmospheric irregularities, instrument aberrations, diffraction limit and detector noise can also cause the observed image to deviate from the original one. Other type of restoration problem may be related to filtering of noisy stars where for example in a nebula image a mass of stars may be extremely bright but also may be spread randomly in dark space, and the shape of the nebula may therefore appear obscure or deviated from the original one. In addition, there may be some astronomical images such as that of Milky Way and galaxy of stars in which some important information may be required about a major star or some other astronomical object surrounded by mass of stars. To restore the original appearance of a nebula or other astronomical objects noisy stars must be filtered out and the detailed structure of the nebula and the other object is to be enhanced.

The various image restoration techniques available in literature include statistical approaches wavelet based methods and partial differential equation (PDE) based regularization techniques based on total variation (TV) and variation principles1,3-8,13-16. The various statistical estimation approaches8 for restoration of astronomical images include maximum likelihood method (MLE), the maximum entropy (ME), and the Bayesian methods. The various types of commonly used filters based on statistical concepts available in literature for restoration of images are averaging filter Median Filter Gaussian Filter,Wiener filter and relaxed median filter 2,5,6.

The anisotropic diffusion based filter 3 for removal of additive noise from digital images reads: .

$\frac{\partial I}{\partial t}=div\left(C\left(|\nabla I|\right)\nabla I\right)$            (1)

Where the diffusion coefficient is defined as 3

$c\left(|\nabla I|\right)=\frac{1}{1+{\left(\frac{|\nabla I|}{k}\right)}^{2}}$              (2)

In above diffusion coefficient,is the edge threshold parameter and ranges from 20 to 80 for digital gray images as defined and proposed by Persona3. A more general diffusion based filter i.e. complex diffusion based filter1 originally proposed for image enhancement and additive noise removal from general digital gray images reads:

$\frac{\partial I}{\partial t}=div\left(c\left(\mathrm{Im}\left(I\right)\right)\nabla I\right)$                           (3)

For nonlinear complex diffusion, the diffusion coefficient is defined as follows1:

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

Here, k is edge threshold parameter and the value of k ranges from 1 to 1.5 for digital images1. The value of k is fine tuned according to the application in hand and for experimentation purposes1 value of$\theta$ is chosen to be $\frac{\pi }{30}$

The anisotropic diffusion-based model 3 described by Eqn. (1) and nonlinear complex diffusion based filter1 described by Eqn. (3) have been successfully applied for the restoration of digital gray images and both methods performs better for additive noise removal. The modified and extended forms of these two methods are examined in this paper for restoration of astronomical images. The only modification that is required in these filters is to make it adaptive in nature so that during filtering procedure all required information is calculated automatically for an image. Further, for filtering of noisy stars and spike noises a relaxed median filter is also used. If the proposed scheme is non-adaptive then these methods have to be fine-tuned for various categories of astronomical images and it is a difficult process because we don’t have any exact idea of degradation process of the image or amount of noise added to it based on which the edge threshold parameter k can be calculated which is to be used by the diffusion coefficients. These diffusion based schemes are proposed to be made adaptive by introducing the calculation of the edge threshold parameter k used in calculation of diffusion coefficients described by Eqn. (2) for anisotropic diffusion 3 , and by Eqn. (4) for nonlinear complex diffusion-based filters 1. Normally, in most of the cases the edge threshold parameter k which is a constant is set by hand and fine-tuned according to an image as originally proposed by the authors 1,3. The authors of the paper 4 proposed to make the anisotropic diffusion based scheme an adaptive one by setting the values of edge threshold parameter k to the local variance of the image which reads: .

$k={k}_{0}{\sigma }^{2}$                (5)

where ${\sigma }^{2}$ is the local variance of the image calculated in a window of size 3x3 and k0 is a constant to be fixed by hand. The value of a very large local variance in comparison to magnitude of the gradient of the image is controlled by the value of k0. Depending upon the requirement of filtration of spike noises due to mass of stars from nebula images the value of k0 was set manually which are normally 0.1, 0.5, 1 and 5. After careful examination of this scheme it can be said that it is not totally adaptive in nature and requires manual intervention for setting the values of k0.For best results the value of k should be determined adaptively according to the application in hand rather than fixing it to some constant value throughout the various iterations of the PDE till its convergence.

Therefore, to make the diffusion based scheme overall adaptive in nature, in this paper, the adaptive value of k is proposed to be determined17, 18 using the tools from robust statistics to automatically estimate the robust scale ${\sigma }_{e}$ of an image I. The value of k is set to ${\sigma }_{e}$ which is minimum absolute deviation (MAD) of the gradient of an image. The value of k is estimated as

$k={\sigma }_{e}=1.4826×media{n}_{I}\left[‖\nabla I-media{n}_{I}\left(‖\nabla I‖\right)‖\right]$                           (6)

This method of calculation of k makes the anisotropic diffusion and nonlinear complex diffusion-based PDEs totally adaptive in nature and is well capable of removing the additive noise from any astronomical image in an effective manner.

Further, in astronomical images, the other type of restoration problem may be related to filtering of noisy stars and spike noises. For example, in a nebula image a mass of stars may be extremely bright but also may be spread randomly in dark space, and the shape of the nebula may therefore appear obscure or deviated from the original one. To restore the original appearance of a nebula or other astronomical objects noisy stars must be filtered out and the detailed structure of the nebula and the other object is to be enhanced. Therefore, a complete astronomical image restoration procedure should also account for these types of noises in addition to removal of additive noise as addressed by the proposed adaptive diffusion-based PDEs.To remove the noisy stars or spike or speckle like patterns from the smoothed astronomical image obtained from the application of the proposed diffusion based PDEs, we propose to apply relaxed median filter2 to the smoothed version of original noisy image obtained in every iteration of diffusion-based PDEs till its convergence.

$\frac{\partial I}{\partial t}=RM\left(div\left(C\left(|\nabla I|\right)\nabla I\right)\right)$                                 (7)

where $C\left(|\nabla I|$ is the diffusion coefficient defined by Eqn. (2); the value of edge threshold parameter k used in diffusion coefficient is calculated adaptively by Eqn. (6) in every iteration of the PDE and RM is the relaxed median filter2 with lower bound and upper bound . Similarly, the proposed adaptive nonlinear complex diffusion based model reads

$\frac{\partial I}{\partial t}=RM\left(div\left(c\left(\mathrm{Im}\left(I\right)\right)\nabla I\right)\right)$                                 (8)

where, c(Im(I)) is the diffusion coefficient defined by eq(4); the value of edge threshold parameter k used in diffusion coefficient is calculated adaptively by Eqn. (6) in every iteration and RM is the relaxed median filter2 with lower bound and upper bound b . The two bounds of the relaxed median filter2 define a sub-list inside the window [W(i, j)] that contains the gray levels that need no filtering. If X(i, j) is the output of a relaxed median filter then X(i, j) is given by:

$X\left(i,j\right)=R{M}_{\alpha ,\beta }\left\{W\left(i,j\right)\right\}$                                 (9)

$=X\left(i,j\right)$ if $X\left(i,j\right)\in \left[{\left[W\left(i,j\right)\right]}_{\left(\alpha \right)}{\left[W\left(i,j\right)\right]}_{\left(\beta \right)}\right]$

$={\left[W\left(i,j\right)\right]}_{\left(m\right)}$otherwise;

where [W(i, j)](m) is the median value of the samples inside the window W(i, j).The sliding window is

$W\left(i,j\right)=\left\{{X}_{\left(i+r,j+r\right)}:r\in W\right\}$                                 (10)

located at position i .

For digital implementations, the proposed PDE based schemes can be discretized using finite difference scheme. The discretized form of the proposed hybrid and adaptive anisotropic diffusion-based PDE given by Eqn. (7) reads

${I}^{n+1}\left(x,y\right)={I}^{n}\left(x,y\right)+\Delta t.\left[RM\left(\stackrel{⇀}{\nabla }.c\stackrel{\to }{\nabla }I\right)\right]$                                 (11)

The second term in R.H.S. of Eq. (11) can be further discretized using centred difference scheme as proposed by Persona3.

Similarly, the discretized form of the proposed hybrid and adaptive nonlinear complex diffusion-based PDE given by Eq. (8) reads

${I}^{n+1}\left(x,y\right)={I}^{n}\left(x,y\right)+\Delta t.\left[RM\left(\stackrel{⇀}{\nabla }.c\left(\mathrm{Im}\left(I\right)\right)\stackrel{\to }{\nabla }I\right)\right]$                                 (12)

For the numerical scheme, given by Eqs. (11-12) to be stable, the von Neumann analysis 19, shows that we require $\frac{\Delta t}{{\left(\Delta x\right)}^{2}}$ < ¼. If the grid size is set to $\Delta x=1$ , then $\Delta t$ < ¼ i.e. $\Delta t$ < 0.25. Therefore, the value of $\Delta t$ is set to 0.25 for stability of Eqs. (11-12).

Out of the two proposed hybrid and adaptive schemes, the hybrid and adaptive version of nonlinear complex diffusion-based technique is outperforming all the schemes in consideration which is validated by the results.

The performance measurement metrics used are blurred signal-to-noise ratio (BSNR) and average signal-to-noise ratio (SNR) which are defined as follows:

$Avg.SNR=\frac{mean\left(I\right)}{std.deviation\left(I\right)}=\frac{\mu }{\sqrt{{\sigma }_{n}^{2}}}$                                 (13)

A higher value of avg. SNR indicates less noise i.e. better restoration.

$BSNR=10\mathrm{log}10\left[\frac{\frac{1}{\left(M×N\right)}\sum \sum \left[{I}_{orignoisy}-{I}_{restored}{\right]}^{2}}{{\sigma }_{n}^{2}}\right]$                                 (14)

Where M×N is the size of image,${I}_{orignoisy}$ is the original noisy image; ${I}_{restored}$ is the restored filtered image and ${\sigma }_{n}^{2}$ is the variance of noise in original noisy image. A low value of BSNR indicates less blurring or less noise or better restoration. The main issue with the use of BSNR is the estimation of noise variance ${\sigma }_{n}^{2}$ in the image. There are various schemes that exist in literature for the estimation of noise variance from an image20. In this paper, an adaptive noise estimation technique20 is used which is defined as follows:

${\sigma }_{n}^{2}=\frac{1}{M×N}\sum \sum$ mode $\left({\sigma }_{v\left(i,j\right)}^{2}\right)$                                 (15)

The mode of a sample data can be calculated by following relation in terms of mean and median of the sample image data:

(16)

The results of the proposed schemes and its comparison with other standard techniques are presented for three astronomical images21. The two proposed schemes described by Eqs. (11-12) have been implemented using MATLAB 7.0. The values of edge threshold parameter for both schemes were calculated adaptively and is described by Eq. (6). The values of k were calculated during each iteration of the evolution of PDEs described by Eqs. (11-12) and in subsequent iterations its value decreases as the restored image evolves in every iteration till its convergence. It has been tested experimentally that the proposed schemes converges after 50 iterations for most of the astronomical images and gives acceptable quality of the restored image. After 50 iterations the image starts blurring and fine structure may be lost. The value of $\Delta t$ is set to 0.25 for stability of Eqs. (11-12). For implementation of Eq. (12) the value of $\theta$ is chosen to be $\frac{\pi }{30}$ . The relaxed median filter 2 is applied to the output produced in each iteration of the proposed PDEs for the removal of spike noise and noisy stars from astronomical images. The two proposed schemes have also been compared against some well known image restoration techniques such as averaging filter5,6, Median Filter5,6, Gaussian Filter5,6, Wiener filter5,6, anisotropic diffusion based filter described by Eq. (8) 3 and adaptive anisotropic diffusion filter based on local variance described by Eqs. (1) and,(5)4 and their efficacy is examined for astronomical images. During implementations of the anisotropic diffusion filter 3 and adaptive anisotropic diffusion filter based on local variance 4 the values of k and k0 were set to 10 and 0.1 respectively. The numbers of iterations for all diffusion-based filters were set to 50. The proposed scheme was tested for several astronomical images and the results of three sample astronomical images are shown in this paper. The performance trend is almost same for all the images. Table 1 lists the average SNR comparison and Table 2 lists the BSNR comparison of the various restoration schemes, respectively, for the three sample images in consideration.

Table 1. Average SNR comparison of proposed scheme with other schemes.

Table 2. BSNR comparison of proposed scheme with other schemes.

The SNR for the three noisy sample images shown in this paper, Table 1, are 1.1957, 1.2721, and 0.7799 respectively. After applying the proposed scheme there are significant improvement in the SNR of the sample images which are 1.9079, 1.4011, and 1.6503 respectively over the sample noisy images in comparison to the SNR of another recent research

Sample image 1 is aat111h.jpg of size 101x126; sample image 2 is orionnebula.jpg of size 227x200 and sample image 3 is m52rb.jpg of size 259x20021.

Figure 1. Results of various image restoration schemes for Sample image 1: aat111h.jpg, 101x126 21.

Figure 2. Results of various image restoration schemes for Sample image 2: Orionnebula.jpg ,227x200 21.

paper i.e. local variance based adaptive anisotropic diffusion filter4 which are 1.8364, 1.3474, and 1.2878 respectively. The BSNR for the three noisy sample images shown in this paper, Table 2, are 18.817, 13.705, and 24.947 respectively. After applying the proposed scheme there are significant reduction in the BSNR of the sample images which are 11.5150, 12.082, and 11.2850 respectively in comparison to the BSNR of the recent research paper i.e. local variance based adaptive anisotropic diffusion filter4 which are 12.7965, 12.867, and 14.9159 respectively. The performance trend remained same for most of the test images in consideration.

Sample image 1 is aat111h.jpg of size 101x126; sample image 2 is orionnebula.jpg of size 227x200 and sample image 3 is m52rb.jpg of size 259x20021.

Figs. (1-2) show the visual results for various image restoration schemes for two sample images. From Figs. 1 and 2, it can be observed that among all eight restoration schemes that have been used for comparison the Gaussian, averaging, Wiener, median, anisotropic diffusion 3 and adaptive anisotropic diffusion based filter 4 are not capable of filtering the noisy stars and also not effective for removing additive noises. Out of the two proposed filters, the adaptive and hybrid nonlinear complex diffusion based scheme produces best results and also it associated with larger value of average SNR and smaller value of BSNR in comparison to all schemes (see Table 1 and 2). Therefore, it can be observed that both the proposed schemes, which are hybrid adaptive anisotropic diffusion-based and hybrid adaptive nonlinear complex diffusion-based filters, are performing better than all other standard restoration schemes for all sample images. Further, out of the two proposed schemes hybrid adaptive nonlinear complex diffusion filter is outperforming all schemes in consideration and producing the better and required results for astronomical images. The performance trend remained the same for all other sample images.

In this paper, two hybrid and adaptive diffusion-based filters were proposed for blind restoration of astronomical images. Their efficacy and performance were examined in comparison to other standard image restoration schemes in terms of average SNR and BSNR. Out of the two proposed schemes which are hybrid adaptive anisotropic diffusion-based filter and hybrid adaptive nonlinear complex diffusion-based filter, the later outperforms all other schemes and is a better and suitable choice for restoration of astronomical images. It is capable of removing all possible additive noises from astronomical images that may have been introduced during its acquisition due to various reasons that may include atmospheric irregularities, instrument aberrations, diffraction limit and detector noise. In addition to, these noises the proposed scheme is also capable of removing spike noises due to presence of noisy stars in astronomical images. In addition to the one application in astronomical image processing as discussed in this paper, the proposed scheme may also be utilized for restoration and enhancement of images, degraded in a manner as described in this paper, in other application areas including defence applications where we are interested in removal of arbitrarily scattered granular noises and spike noises.

1. Gilboa, G.; Sochen, N. & Zeevi, Y.Y. Image enhancement and denoising by complex diffusion processes. IEEE Trans. Pattern Anal. Mach. Intel., 2004, 25(8), 1020-36.

2. Hamza, A.B.; Escamillia, P.L.; Aroza, J.M. & Roldan, R. Removing noise and preserving details with relaxed median filter. J. Math. Imaging Vis., 1999, 11(2), 161-177.

3. Perona P & Malik J. Scale space and edge detection using anisotropic diffusion. IEEE Trans Pattern Anal Mach Intell, 1990, 12(7), 629–39.

4. Chao, S.-M Tsai, & D.-M. Astronomical image restoration using an improved anisotropic diffusion, Pattern Recognition Letters, 2006, 27(5), 335–344.

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. Starck, J.L.; Pantin, E. & Murtagh, F. Deconvolution in astronomy: A review. Astron. Soc. Pacific, 2002, 114, 1051–1069.

8. Molina R. et. al. Image restoration in astronomy- a Bayesian perspective. IEEE signal Processing magazine, 2001,18(2), 11-29.

9. Molina, R.; del Olmo, A.; Perea, J. & Ripley, B.D. Bayesian deconvolution in optical astronomy. Astron. J., 1992, 103(689), 666-75.

10. Molina, R.; de la Blanca, Pérez N. & Ripley, B.D. Statistical restoration of astronomical images. Data Analysis in Astronomy III, V. di Gesu, L. Scarsi, P. Crane, J.H. Friedman, S. Levialdi and M.C. Maccarone, Eds. New York. 1988. pp.75-82.

11. Molina, R. & Ripley, B.D. Using spatial models as priors in astronomical images analysis. J. Appl. Statist., 1989,16(2), 193-206.

12. Bouanno, R.; Buscema, G.; Corsi, C.E. et al. Automated photographic photometry of stars in globular clusters. Astron. Astrophys., 1983, 126(2), 278-82.

13. Starck, J.L.; Candes, E. & Donoho, D.L.. Astronomical image representation by the curvelet transform. Astron. Astrophys., 2003, 398(2), 785–800.

14. Mihcak, M.K.; Kozintsev, I.; Ramchandran, K. & Moulin, P. Low complexity image denoising based on statistical modeling of wavelet coefficients. IEEE Signal Process. Lett., 1999, 6(12), 300–303.

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

16. 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.

17. Black, M.J. & Sapiro, G.. Robust anisotropic diffusion. IEEE Trans Image Process 1998, 7(3), 421-432.

18. Voci, F.; Eiho, S.; Sugimoto, N. & Sekiguchi, H. Estimating the gradient threshold in Perona-Malik equation. IEEE signal processing magazine. 2004, 21(3), 39-46.

19. 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.

20. Aja-Fernández, S. et al. Automatic noise estimation in images using local statistics: Additive and multiplicative cases. Image and Vis. Comput., 2009, 27(6), 756–770.

21.National Optical Astronomy Observatory/Association of Universities for Research in Astronomy/National Science Foundation. U.S.A. http://www.noao.edu/image_gallery.

Dr Rajeev Srivastava received his PhD (Computer Engineering) from University of Delhi, Delhi, India. Presently working as an Associate Professor in the Department of Computer Engineering, Indian Institute of Technology, Banaras Hindu University, Varanasi, India. He has around 15 years of teaching and research experience. He has published around 32 research papers in reputed journals and conferences; 04 book chapters; edited one research reference book. His research interests include : Image processing and computer vision, pattern recognition, soft computing, bioinformatics algorithms and e-learning.