Abstract
Heliograph imaging is the process to reconstruct the solar image from sparse frequency domain data, and compressed sensing (CS) algorithm has shown potential power to accurately recover images from highly undersampled data. However, the available compressed sensing based models available for heliograph imaging are not able to reconstruct fine solar structures and often suffer from undesired convolutive artifacts. This paper presents an imaging model with total variation constraint and nonlocal operator based on compressed sensing theory, with particular objective to suppress convolutive artifacts and reconstruct fine structures. In particular, an efficient algorithm is presented to solve the formulated model. Finally, a set of simulations has been conducted by using both synthetic and real images, and the results demonstrate that our proposed algorithm has surprisingly lower reconstruction errors than other methods.
1. Introduction
Heliograph is a large-scale antenna array used to provide radio solar images. Chinese Spectral Radio Heliograph (CSRH) is an imaging spectroscopy over centimeter and decimetric wavelength range with high temporal, spatial, and spectral resolutions. The instrument is important for addressing fundamental problems of energy release, particle acceleration, and particle transport. The project of the CSRH has been supported by National Major Scientific Research Facility Program of China and is under construction now, which will open new observational windows on flares and CMEs at radio wavelengths.
According to the program, CSRH consists of 40 antennas of 4.5 m covering 400 MHz–2 GHz and 60 antennas of 2 m covering 2–15 GHz. The array directly probes the Fourier plane associated with the image tangent to the celestial sphere at the solar position and generates the visibilities. Visibility coverage in the Fourier plane is incomplete due to the finite number of antennas. Consequently, recovering solar image from measured visibilities requires solving an ill-posed deconvolution problem.
So far various image reconstruction methods have been proposed to deal with these missing data. For example, linear methods set missing data as zero and then get the solar image, also called solar brightness distribution, through the inverse Fourier transform. These simple methods are generally used to acquire the initial solution of the nonlinear methods. Nonlinear methods extrapolate these missing data based on some physical facts and statistical properties in order to get the “best” estimate of the solar brightness distribution. These nonlinear methods can be divided into two categories: the CLEAN algorithm and the Maximum Entropy Method [1, 2]. Variations and further developments of these two methods have been subsequently proposed by researchers [3, 4]. Recently, compressed sensing has gained great success in image reconstruction [5, 6]. Compressed sensing takes advantage of the signal’s sparseness or compressibility in some domain, allowing the entire signal to be determined from relatively fewer measurements. Wiaux et al. first introduced compressed sensing to radio interferometry and illustrated its superiority compared with imaging methods available [7]. Since then several papers are published on radio image reconstruction using CS methods. In 2011, McEwen and Wiaux generalized the compressed sensing imaging techniques proposed in [7, 8] to recover the interferometric images defined directly on the sphere [9, 10]. Wenger et al. presented a practical method for interferometric image reconstruction that was evaluated on observational data from the Very Large Array (VLA) telescope [11].
Many signals are also sparse or compressible in the magnitude of their gradient space, in which case the total variation (TV) based methods have demonstrated their power in image reconstruction [12]. Wiaux et al. used a reweighted TV minimization approach to recover the signal corrupted by cosmic microwave background distributions [13].
Despite the high effectiveness in image reconstruction, compressed sensing based algorithms often suffer from undesired convolutive artifacts and tend to oversmooth image details. Consenquently, the algorithms fail to recover fine solar structure when used to heliograph imaging. To overcome the drawback, an imaging model based on total variation constraint and nonlocal operator under compressed sensing theroy (hereinafter referred to as CSNL) is proposed in this paper. Besides the total variation constraints, the new designed model introduces nonlocal operator to regularize the ill-posed problem. Minimizing total variation (TV) ensures the sparsity in gradient domain, while the nonlocal operator takes advantage of similarity to improve resolution of fine sturctures [14]. The TV constraint and nonlocal operator are complementary to each other, making the proposed model highly effective in recovering the full disk of solar image and fine structures when used in heliograph imaging.
Our key contributions are the introduction of compressed sensing algorithm to heliograph imaging and the new development of an imaging model with total variation constraint and nonlocal operator. Moreover, an efficient algorithm is presented to solve the formulated model. The rest of this paper is organized as follows. Section 2 describes the basic theory of compressed sensing based signal reconstruction; heliograph imaging model based on TV constraint and nonlocal operator is researched in Section 3. Some simulations are conducted in Section 4, and Section 5 draws some conclusions.
2. Signal Reconstruction Based on Compressed Sensing Theory
In the compressed sensing framework, the reconstruction of a signal that is sparse in a bases system requires just a small number of measured samples , where . This subsampling process can be represented as a projection by measurement matrix . Therefore, the observable can be expressed as follows: The newly defined matrix represents an overcomplete basis.
Equation (1) expresses an underdetermined system of linear equations and the solution is nonunique. To solve this equation, CS assumes that the signal is sparse, which means that the number of the nonzero coefficients of in -domain, that is, tends to be as small as possible:
Apparently, the minimization of (2) is a combinatorial problem that is computationally intractable. Fortunately, the solution of (2) for L0 normal is identical to the solution to a more tractable L1 normal problem when the signal is sparse. So, the CS recovery of from is formulated as the following constrained optimization problem:
3. Heliograph Imaging Model Based on TV Constraint and Nonlocal Operator
Heliograph does not directly measure the brightness distribution of solar image but samples the complex visibility values on the UV plane through the array, and finally the solar image is got through inverse Fourier transform from the obtained data. Therefore, heliograph imaging can be understood as a process of the recovery of the original image from the corresponding complex visibility value and the sampling distribution determined by the given array configuration.
3.1. Imaging Model Based on Compressed Sensing of Minimizing Total Variation
Total variation is widely used to express sparseness in 2D image [15]. According to the compressed sensing theory, the imaging model based on compressed sensing of minimizing total variation for heliograph is as follows: where is a noise-related constant greater than 0, is the two-dimensional discrete Fourier transform matrix, and is the sampling distribution. is the sum of the discrete gradient at each point of the image : where .
3.2. Imaging Model Based on TV Constraint and Nonlocal Operator
The model on compressed sensing of minimizing total variation has good performance in reconstructing of extend sources such as the disk of the sun while it fails to recover the fine solar structure. The reason is that minimizing the total variation tends to oversmooth the gradient of images. To overcome the drawback, the nonlocal operator is introduced. Nonlocal operator is effective in preserving image details by averaging the nearby pixels with similar structures.
For a given pixel in an image , its nonlocal mean, denoted by , is obtained as a weighted average of surrounding pixels within a neighborhood domain, and the nonlocal mean shift of , , is given by the difference between and , namely, as follows: where is the pixel in the neighborhood domain, is the pixel numbers in the neighborhood window, and is the weight of similarity between the pixels and. This similarity is measured as a decreasing function of the weighted Euclidean distance or other-defined distance. To make the algorithm less time consuming, is set a constant scalar determined only by the Euclidean distance between pixels of and in dirty map according to the following equation no matter what the values of and are:
Hence, , together with , can be considered as the weighted sum of . Let with being the similarity weight matrix. As is set to be a constant scalar previously, is a constant matrix for given and , hence lowering the computational burden.
So heliograph imaging based on TV constraint and nonlocal operator under compressed sensing theory can be modeled as follows:
Note that the difference between the model proposed and the one in [16] is as follows.
The model in [16] addressed the image reconstruction under noise-free condition while we relax the problem to that with a noise distribution.
Moreover, we fix a constant matrix instead of updating it from iteration to iteration for the reason of computing, consuming, and convergence.
We introduce auxiliary variables and to transform (7) into
This problem can be recast as second-order cone programs, and (9) can be rewritten as where The standard log-barrier method transforms (10) into a series of linearly constrained programs: where is a slack variable with . The inequality constraint has been incorporated into the functional via a penalty function which is infinite when the constraint is violated and smooth elsewhere. As gets large, the solution to (12) approaches the solution to (10): it can be shown that ; that is, we are within of the optimal value after iteration . The idea here is that each of the smooth subproblems can be solved to fairly high accuracy with just a few iterations of Newton’s method, especially since we can use the solution as a starting point for subproblem .
The complete implementation for this problem follows the following outline.
Step 1 (Initialization). Consider a feasible solution , a tolerance , and zoom factors and . Set , , and and , , and , is the vectorized dirty image.
Step 2. Solve (10) using Newton’s method, using as an initial point to get the solution .
Step 3. If , terminate and return , else , , and go to Step 2.
4. Experimental Results
In the following, we evaluate the performance of our approach by conducting two simulation experiments. A simple synthetic solar image and a more complex image of Cygnus A are used as the simulated sources and considered as “ground truth” image.
The objective fidelity of the reconstructed image is measured by the signal to noise ratio (SNR): and the relative error (RelErr) between the reconstructed image and the “ground truth” image : , where represents the standard deviation.
In our experiments, the antenna array, composed of 100 parabolic antennas with diameter of 3 m, is configured as a three-arm spiral topology with the longest baseline 3 km and the shortest baseline 6 m. Figures 1, 2, and 3 show the topology, UV coverage and, dirty beam (point spread function, PSF), respectively.
The spatial frequencies are just the distances expressed in wavelength units, and the coordinates are obtained from the baseline from the following coordinate transformation: where and are the hour angle and declination of radio source, respectively, and , the wavelength of the observed radio-wave, is set to be 0.3 m.
The synthetic solar brightness image is shown in Figure 4(a) with a size of 512 * 512 pixels in the first experiment. The working frequency is 1 GHz, hence the pixel resolution is 5 arcsec. We generated three different brightness “needle-like” point sources, several different brightness arc sources, and several adjacent point sources with same brightness to verify the performance of dynamic range, shape maintenance, and resolution, respectively.
(a)
(b)
(c)
(d)
Given the test image and UV sampling distribution , we can calculate the simulated measurements and then add Gaussian zero mean white noise with (the dynamic range of the gray level is 0 to 255) to visibilities .
As a reference, Cotton-Schwab clean algorithm (CSC), which is embodied in the AIPS, and compressed sensing based algorithm are used to reconstruct the simulated sources. Figures 4(b), 4(c), and 4(d) give the reconstructed solar images employing the CSC, CS, and CSNL, respectively. The CSC is terminated at 500 iterations and the loop gain is set to be 0.1. In CS, the basis is set to be the identity matrix, and the model is solved by the orthogonal matching pursuit algorithm.
The results indicate that the proposed algorithm has overwhelming advantages over the CSC algorithm and CS algorithm in terms of dynamic range, shape maintenance, and resolution. CSC can only restore several point sources but fails to recover the shape of arc sources and the full disk of solar image. CS can restore not only the shape of arc sources and the full disk, but also most of the point sources. As a contrast, CSNL can reconstruct all the details of the solar image as a result of introduction of the nonlocal operator. The nonlocal operator is very effective in sharpening edges, which helps CSNL in reconstructing the full disk; meanwhile, nonlocal operator gains good denoising performance, which benefits CSNL by restoring the point sources with low brightness.
To simulate images from astronomical radio sources, we just take image of the giant radio galaxy Cygnus A from NRAO image gallery and vectorize it [17]. If the radio source is thought of as an array of point sources, the brightness value of each pixel in the image represents the intensity of each point source. So the second experiment is conducted by using the test image of Cygnus A, shown in Figure 5(a), with the size of 512 * 512 pixels. Figures 5(b), 5(c), and 5(d) give the reconstructed images employing the CSC, CS, and CSNL, respectively. From the images, it is obvious that the CSC suffers from striping artifacts and CS suffers from both artifacts and an uneven background. Benefiting from the nonlocal operator, which exploits the nonlocal self-similarity and eliminates the artifacts, CSNL is superior to CSC and CS in image quality and dynamic range.
(a)
(b)
(c)
(d)
As we have seen, CSC is good at the reconstruction of point sources and CS gives better quality images for images consisting of extended sources. Generally speaking, CS is superior to CSC for it is regularized by the sparse constraint. But both CSC and CS suffer from undesired convolutive artifacts. By incorporating nonlocal operator, CSNL provides much better visual reconstruction quality and SNR than CSC and CS.
5. Conclusions
This paper presents an improved CS based image reconstruction algorithm by incorporating nonlocal operator. The total variations in compressed sensing model and nonlocal operator are complementary to each other, making the new proposed algorithm highly effective in reconstructing the solar disk and preserving fine structure. The results of our experiments demonstrate the potential applications of the proposed method.
Conflict of Interests
The authors declare that they have no financial and personal relationships with other people or organizations that can inappropriately influence their work; there is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, the paper.
Acknowledgment
This work is supported by the National Natural Science Foundation of China under Grants Nos. 61100156, 61070143, and 61201295.