AN INTERFERENCE FRINGE REMOVAL METHOD BASED ON MULTI-SCALE DECOMPOSITION AND ADAPTIVE PARTITIONING FOR NVST IMAGES
2College of Computer and Information Technology, China Three Gorges University, Yichang 443002, China ctguhuangyao@163.com
Published under Creative Commons license CC BY-SA 4.0
Abstract
The New Vacuum Solar Telescope (NVST) is the largest solar telescope in China. When using CCDs for imaging, equal-thickness fringes caused by thin-film interference can occur. Such fringes reduce the quality of NVST data but cannot be removed using standard flat fielding. In this paper, a correction method based on multi-scale decomposition and adaptive partitioning is proposed. The original image is decomposed into several sub-scales by multi-scale decomposition. The region containing fringes is found and divided by an adaptive partitioning method. The interference fringes are then filtered by a frequency-domain Gaussian filter on every partitioned image. Our analysis shows that this method can effectively remove the interference fringes from a solar image while preserving useful information.
Keywords:
Sun: chromosphere, techniques: image processing1. INTRODUCTION
The one-meter New Vacuum Solar Telescope (NVST) is the primary equipment for optical and near-infrared observations of the Sun in the Chinese solar physics and space science communities (Li et al. 2014). The main goals of NVST are high resolution imaging and spectral observations, including measurements of the solar magnetic field. When the telescope is observing, the solar chromosphere imaging system is sometimes affected by interference fringes caused by thin film interference in the CCD system (Liu 2016; Wei 2011). The spacing of the interference fringes is constant over time, but the intensity of the fringes varies significantly. Instead of decreasing, some of the interference fringes in the NVST image are enhanced after traditional flat fielding as shown in Figure 1a,b (Wang et al. 2016). In Figure 1c,d, the interference fringes appearing in the shown region of the solar image are more apparent. The presence of interference fringes can seriously affect the quality of solar observation data products and interfere with subsequent solar physics research. This paper studies a method to improve the quality of the data by removing interference fringes from the perspective of image processing.
In image processing, a multi-scale wavelet transform can decompose the original image into multiple scale maps, and the fringe features will be clearly displayed on some scale maps. This approach has achieved good results in preserving image detail and fringe removal (Liang et al. 2009; Sun et al. 2016). Previously, He Shuan et al. processed the Huairou solar images using a multi-scale wavelet decomposition and fast fourier transformation (FFT), and achieved positive results. However, the method works best when applied to straight line fringes with uniform directions. The interference fringes in the NVST image data present themselves as curves which cannot be directly removed using existing methods, as shown in Figure 1 (He et al. 2014).
At present, there is no effective method for removing NVST interference fringes. Given that the curvature of the fringes in the NVST image does not vary rapidly, the area containing the fringes can be divided into small blocks. The direction of the fringes in each of the small blocks is approximately uniform and linearly distributed, and hence can be filtered in the frequency domain. Therefore, this paper proposes a new fringe removal method based on multi-scale decomposition and adaptive partitioning for NVST images. First, the multi-scale wavelet transform is utilized to divide the original image into multiple scales. Second, the Radon transform is utilized to detect the line direction, and identify all blocks containing fringes adaptively. Third, a frequency domain Gaussian filter is utilized to remove the fringes in all blocks.
2. OBSERVATIONAL DATA
The NVST is a vacuum solar telescope with a clear aperture of 985 mm which began operations in 2012 (Xiang et al. 2016). Its main structure is a multi-channel high resolution imaging system which consists of one channel for the chromosphere and two channels for the photosphere. All channels can observe and record images simultaneously. The band used for observing the chromosphere is Hα(6563Å). The Hα filter is a tunable Lyot filter with a bandwidth of 0.25 Å (Liu et al. 2014). The Hα light is recorded using a 14-bit CCD camera pco4000 manufactured by the PCO company. The exposure time of each frame is 10 or 20 ms. The camera has 2K×2K square pixels which cover a field of view (FOV) of 168′′×168′′. The image scale of each pixel is 0.082′′. In order to improve the signal-to-noise ratio (SNR) of the data and shorten the sampling interval between images, the image can be binned to 1K×1K by using the binning mode of the camera. After binning, the scale of each pixel is 0.164′′.
In this paper, we use Hα data of the chromosphere for our fringe removal experiments. To illustrate the effectiveness of our algorithm, we select NVST imaging data with obvious fringes. We selected two datasets taken on 2013 June 12 and 2014 December 23. Additional data taken on 2014 August 27, 2015 December 1, 2016 November 20, 2017 August 6, and 2018 August 20 were used to verify the efficiency of the method in simulations.
3. ALGORITHM
3.1. Multi-Scale Wavelet Transform
The multi-scale decomposition makes use of a multi-scale wavelet transform. From the large number of existing multi-scale wavelet transform methods, we select the “à trous” wavelet algorithm to decompose the image into wavelet planes (Núñez et al. 1999; Chibani 2005; Holschneider et al. 1989). We construct a sequence of approximations from a given image p:
$$$\begin{array}{ccc}{F}_{1}\left(\mathbf{p}\right)={\mathbf{p}}_{1},& {F}_{2}\left(\mathbf{p}\right)={\mathbf{p}}_{2},& {F}_{3}\left(\mathbf{p}\right){=\mathbf{p}}_{3},\end{array}\cdot \cdot \cdot \cdot $$$ | (1) |
F_{1}, F_{2}, ... represent the sub-components of the multi-scale wavelet transform to different frequency segments, and the sub-components correspond to the output of different frequency segment filters. To construct the sequence, the algorithm performs successive convolutions with a filter obtained from an auxiliary “scaling function”. We use a scaling function which has a B_{3} cubic spline profile. The use of a B_{3} cubic spline leads to a convolution with a fixed 5 × 5 mask of the form
$$$\frac{1}{256}\left(\begin{array}{ccccc}1& 4& 6& 4& 1\\ 4& 16& 24& 16& 4\\ 6& 24& 36& 24& 6\\ 4& 16& 24& 16& 4\\ 1& 4& 6& 4& 1\end{array}\right).$$$ | (2) |
The wavelet planes are computed as the differences between two consecutive approximations p_{l−1} and pl. Defining w_{l} = p_{l−1} −p_{l} (l = 1, ..., n), in which p_{0} = p, we can write the reconstruction formula like
$$$\mathrm{p}=\sum _{l=1}^{n}{\mathrm{p}}_{l}+{\mathrm{p}}_{r}.$$$ | (3) |
In this representation, the images p_{l} are versions of the original image p at increasing scales (decreasing resolution levels), w_{l} are the multi-resolution wavelet planes, and p_{r} is a residual image. The image p represents each scale component, and the subscript represents the number of the component.
After the multi-scale wavelet transform, the original image is decomposed into multiple scale images, in which the interference fringes are present only at a few specific scales. In other words, the interference fringes are separated out at this step. Figure 2 shows a Hα image and its decomposed images containing the interference fringes. The image was decomposed into 19 scale images by the wavelet toolkit.
3.2. Radon Transform and Adaptive Partitioning
To find the region containing fringes in the scale maps, the Radon transform is a good choice (van Ginkel et al. 2004). This method has the advantage of good anti-noise performance, and is widely utilized for image processing in medical, optical, and total interferometric metrics.
When mapping a line ρ = x cos θ + y sin θ in a plane space (x, y) to a point (ρ, θ) in Radon space, the Radon transform of the continuous image is
$$$R\left(\rho ,\theta \right)={\iint}_{D}f\left(x,y\right)\delta \left(\rho -x\mathrm{c}\mathrm{o}\mathrm{s}\theta -y\mathrm{s}\mathrm{i}\mathrm{n}\theta \right)dxdy$$$ | (4) |
where D is the entire image plane, f(x, y) is the gray value of the image, and the characteristic function is the Dirac function δ(t) = 1 for t = 0, δ(t) = 0 for t ≠ 0. ρ and θ are the distance and angle of a point in ρ−θ space, respectively. The eigenfunction δ(t) limits the integral to be non-zero only along the straight line ρ = x cos θ + y sin θ. The Radon transform can be understood as the projection of an image in space ρ−θ, and each point of the space ρ−θ corresponds to a straight line of the image space. The Radon transform is the integral of the image pixel values on each line and can also be understood as the projection of the image onto the horizontal axis after clockwise rotation.
To define the scale map f, the length of each small block is s, so the small block size is s × s. Then the following matrix can be defined:
$$$f=\left[\begin{array}{cccc}{f}_{11}& {f}_{12}& \cdots & {f}_{1m}\\ {f}_{21}& {f}_{22}& \cdots & {f}_{2m}\\ \vdots & \vdots & \ddots & \vdots \\ {f}_{n1}& {f}_{n2}& \cdots & {f}_{nm}\end{array}\right]$$$ | (5) |
$$$\begin{array}{l}{f}_{ij}=\left\{\begin{array}{l}1:\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{}\mathrm{f}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{e}\mathrm{s},\\ 0:\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{}\mathrm{b}\mathrm{l}\mathrm{o}\mathrm{c}\mathrm{k}\mathrm{}\mathrm{d}\mathrm{o}\mathrm{e}\mathrm{s}\mathrm{}\mathrm{n}\mathrm{o}\mathrm{t}\mathrm{}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{}\mathrm{f}\mathrm{r}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{e}\mathrm{s},\end{array}\right.\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}1\le i\le m,1\le j\le n.\end{array}$$$ | (6) |
In this paper, we utilize the Radon transform to determine the peak points of the first-order differential image to realize adaptive partitioning. If the positions of the highest intensity peak points are distributed along a column line in a certain block image, where the largest and smallest difference of the first four peak points is less than 10 degrees, we assert that the block image contains fringes. After the Radon line detection in the scale map, the component f_{ij} will be set to unity because the block contains fringes.
The adaptive partitioning is done primarily to automatically calculate the position of each block containing fringes within the full-block scale image. An example is presented in Figure 3. In order to approximate the fringes in each block as straight lines with the same slope, the length of each block is 50 pixels. The blue lines represent all blocks, and the red lines represent the blocks with fringes selected using the Radon line detection method. As can be seen from Figure 3, the method correctly selects the blocks in which fringes are located. The number of blocks with fringes per scale is calculated using the adaptive partitioning method. The 4th, 7th, and 10th scale images have the largest number of fringed blocks. On other scales, the number of fringed blocks is small (less than 40 blocks) and they do not conform to the space size of the fringes. Therefore, the 4th, 7th, and 10th scale images are the scales that must be processed. Figure 4 shows the number of fringed blocks at each scale. In order to check the robustness of the selection of fringed blocks, we calculate the relative error from the ratio of the number of fringed blocks not selected to the total number of blocks in the 4th, 7th, and 10th scale images; errors are less than 10%.
In Figure 5 we define x′ = ρ. Figure 5a shows the original block image with fringes, Figure 5b shows the angular projected intensity image generated with the Radon transform, where the position of the column represents the angles. Figure 5c shows the first-order differential of the projected intensity image, and Figure 5d shows the distribution of the first few peak points in the projected intensity image. Similarly, Figure 5e–h shows the processed image obtained from the non-fringed blocks. The following rules can be determined based on Figure 5. When the first few peak points of the first-order differential image are distributed around a certain angle, the block can be counted as a fringed block. Thus, by positioning we can identify all fringed blocks in the 7th scale. In all the fringed scales, the same adaptive partitioning can be undertaken.
3.3. Frequency Domain Filtering through a Block Image
In order to remove the fringes in the block image, the filter in the frequency domain should account for the periodicity of the fringes. The fringes in the block can be approximated as sine waves along the normal direction. Due to the frequency spectrum distribution of the sinusoidal signal, it can be removed using a bimodal Gaussian filter.
The bimodal Gaussian filter function is:
$$$f\left(x\right)={{A}_{1}e}^{-\frac{{\left(x-{x}_{1}\right)}^{2}}{2{\delta}_{1}^{2}}}+{{A}_{2}e}^{-\frac{{\left(x-{x}_{2}\right)}^{2}}{2{\delta}_{2}^{2}}}$$$ | (7) |
where A_{1},A_{2} are the amplitudes of each Gaussian component; x_{1}, x_{2} is the position of each peak; δ_{1}, δ_{2} is the standard deviation of each Gaussian curve. The filter is shown in Figure 6a.
In order to treat the fringe as sinusoidal signals, the fringed block image is rotated such that the fringes become parallel to the horizontal axis. After rotating, the column direction can be approximated as a sine wave signal. For each block containing fringes, the bimodal Gaussian function (Figure 6a) is utilized to remove the interference fringes which affect the observed image along column direction in the frequency domain.
Figure 7 shows that the positions of the fringes in the 4th and 10th scales are not as clear as in the case of the 7th scale. If the position of the fringes is again determined at the 4th and 10th scales using the Radon transformation, the detection precision is lower and the position of the fringes in each scale is approximately the same. Therefore, we choose to perform Gaussian filtering on each selected fringed block in the 7th scale and record the position of the peak point of each fringe in the frequency spectrogram, and then remove fringes using these positions in the 4th and 10th scales.
The fringe direction is rotated into the horizontal direction in the block, and data are fitted along the vertical direction. In Figure 6b, the solid line is the spectrum of a certain column and the dashed line is the result of fitting by the bimodal Gaussian filter function in Figure 6a. The peak positions of the two Gaussian curves correspond to the peak points in the spectrum of every column. The parameters δ_{1}, δ_{2} are taken to be 1.5. Each fringed block must be filtered column by column and then joined with the unfringed blocks. All blocks filtered on the 7th scale are shown in Figure 7a,b. We can see that there are almost no fringes in the blocks and the scale image.
In Figure 7c–f, the 4th and 10th scales are likewise filtered with the Gaussian filter and the fringes of the corresponding scale are mostly removed.
4. INTERFERENCE FRINGE REMOVAL METHOD
The flow of our fringe removal method is shown in Figure 8.
There are four main steps in our new interference fringe removal method for NVST data. First, the multi-scale decomposition is utilized to divide the original image into multiple scales to realize the separation between signals and fringes. Second, the Radon transform is utilized to detect the line direction, and furthermore identify all blocks containing fringes. Third, a frequency domain Gaussian filter is utilized to remove the fringes in all blocks. Finally, the image reconstruction can be completed with an inverse multi-scale wavelet transform using the corrected scales.
The key problem of the method is finding the block area containing fringes in the scale image adaptively. Since the scale images containing fringes are divided into small blocks, the direction of the fringes in each of the small blocks is basically uniform. Therefore, this paper utilizes the idea of applying the Radon transform to find every block containing fringes adaptively by the method outlined in Section 3. After Gaussian filtering is performed on each fringed scale image, all scales are synthesized by multi-scale inverse wavelet transform using MATLAB, and the final reconstructed image is generated.
5. ANALYSIS RESULTS AND CONCLUSIONS
In this work we have shown that the fringes existing in the original image have been largely eliminated, and the original details of the image are well preserved. In order to demonstrate the effectiveness of the method, we select several NVST chromosphere observation data for testing. Our tests indicate that the method is very effective for removing fringes as shown in Figure 9.
In a further simulation experiment, the extracted fringes are superimposed onto an image without fringes. We use two image quality evaluation methods, the Standard Error (SE) and the Structural Similarity Index (SSIM), to assess the quality of the image restoration (Chae 2004). The smaller the value of the SE, or the closer the value of the SSIM is to unity, the better the performance of the algorithm.
We now define f(x, y) as the image before processing and $$ \stackrel{~}{f}$$(x, y) as the image after processing. The formulas for calculating SE and SSIM are:
$$$SE=\sqrt{\frac{1}{M\times N}\sum _{i=1}^{M}\sum _{j=1}^{N}{\left|f\left(x,y\right)-\stackrel{~}{f}\left(x,y\right)\right|}^{2}}$$$ | (8) |
$$$SSIM\left(a,b\right)=\frac{\left(2{\mu}_{a}{\mu}_{b}+{C}_{1}\right)\left(2{\sigma}_{ab}+{C}_{2}\right)}{\left({\mu}_{a}^{2}+{\mu}_{b}^{2}+{C}_{1}\right)\left({\sigma}_{a}^{2}+{\sigma}_{b}^{2}+{C}_{2}\right)}$$$ | (9) |
where a is f(x, y), b is $$ \stackrel{~}{f}$$(x, y), C_{1},C_{2} are constants, μ_{a}, μ_{b} are the averages of all pixel values in an image, and σ_{a}, σ_{b} are variances of image pixel values.
Five sets of data are selected for the simulation experiment, and the results for the first set of data obtained with adaptive partitioning are shown in Figure 10. After applying three different correction methods – adaptive partitioning, frequency domain filtering, and median filtering (Wang et al. 2016) – the results of calculating SE for all data are shown in Table 1, and the results of calculating SSIM for all data are shown in Table 2.
If we define the image without fringes as f_{1}, the image with overlayed fringes as f_{2}, the restored image from adaptive partitioning as f_{3}, the restored image from frequency domain filtering as f_{4}, and the restored image from median filtering as f_{5}, then se1, se2, se3, se4, ssim1, ssim2, ssim3, and ssim4 represent SE(f_{1}, f_{2}), SE(f_{1}, f_{3}), SE(f_{1}, f_{4}), SE(f_{1}, f_{5}), SSIM(f_{1}, f_{2}), SSIM(f_{1}, f_{3}), SSIM(f_{1}, f_{4}), and SSIM(f_{1}, f_{5}) respectively. It can be seen from Table 1 that compared to the reference value se1, the SE value of the adaptive partitioning method se2 is much lower than than that of the other two methods se3 and se4. Similarly, as seen from Table 2, the SSIM value of the adaptive partitioning method ssim2 is closer to unity than that of the the other two methods ssim3, ssim4. The intensity of the fringes is significantly reduced after treatment by the adaptive partitioning method. From Figure 10, we see that the superimposed fringes are substantially reduced through adaptive partitioning processing, and the recovered data is consistent with the original reference image.
Through multiple observations and analysis of NVST chromosphere data, this paper proposes a method to remove interference fringes from the observed images by multi-scale decomposition and adaptive partitioning. Our results show that the method can effectively remove the interference fringes while retaining the original image information. There remains some areas for improvement of the method. Although relatively obvious fringe areas can be selected in the scale map, a small number of fringe-containing areas are greatly affected by the target image, resulting in no selection processing. The texture features are difficult to test by line detection in the block, and some streaks that are not clearly noticeable will remain.
Acknowledgments
The authors thank the NVST team, especially the observing colleagues for their excellent observations. The data used in this paper were obtained with the New Vacuum Solar Telescope in Fuxian Solar Observatory of Yunnan Astronomical Observatory, CAS. The research work is also supported by the National Natural Science Fund Committee of the Chinese Academy of Sciences astronomical union funds no. U1731124.
References
- Chae, J., (2004), Flat-Fielding of Solar Hα Observations Using Relatively Shifted Images, Sol. Phys., 221, p1. [https://doi.org/10.1023/b:sola.0000033357.72303.89]
- Chibani, Y., (2005), Selective Synthetic Aperture Radar and Panchromatic Image Fusion by Using the à Trous Wavelet Decomposition, EURASIP J. Adv. Signal Process., 14, p2207.
- Holschneider, M., & Tchamitchian, P., (1990), Regularite Locale de la fonction “Non-Differentiable” de Riemann, in: Les Ondelettes en 1989: Seminaire d’Analyse Harmonique, Lemariè, P. G., ed.Lemariè, P. G. (Berlin: Springer-Verlag), p102.
- van Ginkel, M., Luengo Hendriks, C. L. L., & van Vliet, L. J., (2004), A Short Introduction to the Radon and Hough Transforms and How They Relate to Each Other, Quantitative Imaging Group Technical Report Series, QI-2004-01.
- He, S., Zheng, S., & Huang, Y., (2014), Removing the Stripe Noises Interference on the Hα Full-Disk Solar Image Based on Multiscale Transform, Applied Mechanics & Materials, 596, p365. [https://doi.org/10.4028/www.scientific.net/amm.596.365]
- Liang, W., Cui, G., Chen, G. Q., (2009), Feature Edge Extraction Based on Wavelet Multiscale Correlation, Electronics Optics & Control, 16, p86.
- Li, X.B., Wang, F., Xiang, Y.Y., et al. , (2014), Parallel Image Reconstruction for New Vacuum Solar Telescope, JKAS, 47, p43. [https://doi.org/10.5303/jkas.2014.47.2.43]
- Liu, Y.-M., (2016), Hyperspectral Image Debanding Method Based on Adaptive Univariate Variation, Laser & Optoelectronics Progress, 9, p80.
- Liu, Z., Xu, J., Gu, B.-Z., et al. , (2014), New Vacuum Solar Telescope and Observations with High Resolution, Res. Astron. Astrophys., 14, p705.
- Núñez, J., Otazu, X., Fors, O., et al. , (1999), Multiresolution-Based Image Fusion with Additive Wavelet Decomposition, IEEE Transactions on Geoscience and Remote Sensing, 37, p1204.
- Sun, Q., & Qian, Z., (2016), Multi-Scale Wavelet Transform Filtering of Non-Uniform Pavement Surface Image Background for Automated Pavement Distress Identification, Measurement, 86, p26. [https://doi.org/10.1016/j.measurement.2016.02.044]
- Wang, S. B., Xu, Z., Xiang, Y. Y., & Jin, Z. Y., (2016), Study and Removal of the Interference Fringes in Images Observed by the Imaging System of NVST, Acta Astron. Sinica, 57, p608.
- Wei, C., (2011), Airspace Processing and Application Research of Interference Fringe Image, Nanjing University of Science and Technology.
- Xiang, Y. Y., Liu, Z., & Jin, Z. Y., (2016), High Resolution Reconstruction of Solar Prominence Images Observed by the New Vacuum Solar Telescope, New Astron., 49, p8. [https://doi.org/10.1016/j.newast.2016.05.002]