High-resolution Solar Image Reconstruction Based on Non-rigid Alignment

Suppressing the interference of atmospheric turbulence and obtaining observation data with a high spatial resolution is an issue to be solved urgently for ground observations. One way to solve this problem is to perform a statistical reconstruction of short-exposure speckle images. Combining the rapidity of Shift-Add and the accuracy of speckle masking, this paper proposes a novel reconstruction algorithm-NASIR (Non-rigid Alignment based Solar Image Reconstruction). NASIR reconstructs the phase of the object image at each frequency by building a computational model between geometric distortion and intensity distribution and reconstructs the modulus of the object image on the aligned speckle images by speckle interferometry. We analyzed the performance of NASIR by using the correlation coefficient, power spectrum, and coefficient of variation of intensity profile (CVoIP) in processing data obtained by the NVST (1m New Vacuum Solar Telescope). The reconstruction experiments and analysis results show that the quality of images reconstructed by NASIR is close to speckle masking when the seeing is good, while NASIR has excellent robustness when the seeing condition becomes worse. Furthermore, NASIR reconstructs the entire field of view in parallel in one go, without phase recursion and block-by-block reconstruction, so its computation time is less than half that of speckle masking. Therefore, we consider NASIR is a robust and high-quality fast reconstruction method that can serve as an effective tool for data filtering and quick look.


INTRODUCTION
At present, solar physics research has entered a period of small scale and fine structure. High-resolution observational images are urgently needed. However, diffraction-limited imaging cannot be achieved with ground-based optical telescopes due to the interference of the turbulent atmosphere.
To mitigate the adverse effects of the atmosphere and achieve near-diffraction-limited imaging, the researchers propose the speckle imaging method for high-resolution reconstruction by statistical analysis. The speckle imaging method started from the pioneering work of Labeyrie (Labeyrie 1970), and its processing object is the speckle images. Speckle images are short-exposure images taken with a telescope with an aperture much larger than the atmospheric coherence length. The speckle imaging method can obtain high-resolution reconstructed images close to the diffraction limit of the telescope and has a relatively complete theoretical support.
Speckle imaging methods can roughly be classified into spatial domain and frequency domain reconstruction methods. Spatial reconstruction methods directly reconstruct the spatial structure and intensity of the object. The characteristics of the spatial domain method are simple algorithms and clear spatial descriptions. The typical representative of the spatial domain reconstruction method is Shift-Add Bates & Cady 1980. Frequency domain reconstruction methods usually reconstruct the phase and modulus separately. Phase is reconstructed by speckle masking (Lohmann et al. 1983), while modulus is by speckle interferometry (Labeyrie 1970). Most of the above-mentioned classic algorithms originated in the 1970s and 1980s, and subsequent researchers also have substantial and valuable research works to improve these methods (Jin et al. 2008;yuan Xiang et al. 2016;Huo & Zhou 2010;Liu et al. 1999). However, there are still issues that need to be further studied and resolved in practical applications.The main issue of spatial domain methods is to align the regions containing object diffraction imaging information accurately. Taking Shift-Add as an example, the algorithm assumes that the maximum intensity points of the speckle images coincide with that of the object image and performs rigid translation (uniform displacement of the entire image) and superimposition of each speckle image to obtain a reconstructed image. However, for an extended source like Sun, due to the influence of the surrounding structure, the distribution of maximum intensity points will change significantly under the influence of the atmosphere. Not only exist multiple points of maximum intensity, but their positions will also deviate. The assumption that points of maximum intensity coincide each other no longer holds. In addition, the overall displacement of the image ignores the distortion differences caused by different fields of view and frequency components, which makes it difficult to obtain satisfactory reconstruction results with the traditional rigid Shift-Add.
Speckle imaging frequency-domain reconstruction algorithm is currently the most widely used high-resolution reconstruction method (Liu et al. 1999). Well-known ground-based solar telescopes such as NVST (Liu et al. 2014) and GST (Goode Solar Telescope) (Cao et al. 2010) use this method to reconstruct observational data. A typical representative of the frequency domain reconstruction method is speckle masking.
Since speckle masking estimates the phase of the object based on the triple correlation of the speckle images, that is, recursively from the low-frequency phase to the high-frequency, it is inevitably affected by noise such as photon noise during the recursion process. Due to the recursion process, the noise influence will accumulate in the high-frequency part, and in severe cases, it will distort the fine structure of the object. Therefore, the performance of speckle masking will significantly reduce with the decrease in seeing (Liu et al. 1999).
Therefore, we propose a new high-resolution solar image reconstruction algorithm named NASIR in this paper. NASIR draws on the basic theoretical models of speckle masking and Shift-Add and combines the ideas of the optical flow method (Horn & Schunck 1981) in computer vision. The algorithm obtains the distortion displacement field of each frame of the speckle image relative to the reference image (RI) by establishing a computational model between the intensity change and the displacement field. And the distortion is corrected by pixel-by-pixel non-rigid alignment to obtain the phase estimation of the object image. The modulus of the object image is by speckle interferometry on the aligned speckle images.
Section 2 of the paper introduces the principle and realization of NASIR. Section 3 applies the algorithm to reconstruct the observation data of NVST and compares it with the classical methods. Section 4 discusses the features of NASIR and issues that need to be further studied. Our conclusions are in Section 5.

THE PRINCIPLE AND IMPLEMENTATION OF NASIR
Similar to speckle masking, NASIR also reconstructs phase and modulus separately. For phase reconstruction, NASIR obtains the distorted displacement field of each frame of the speckle image relative to the reference image (RI) by establishing a computational model between the intensity variation and the displacement field. The phase estimation of the object image is by correcting the distortion by pixel-bypixel non-rigid alignment. The modulus is by speckle interferometry on the aligned speckle images.
High-resolution Solar Image Reconstruction Based on Non-rigid Alignment 3 The above process involves displacement field estimation, band-pass filtering, iterative alignment, and modulus reconstruction of the object image.

Displacement Field Estimation Based on Intensity Changes
Based on the idea of the optical flow, we established the pixel displacement and intensity variation model of the speckle image. With this model, the phase distortion caused by atmospheric turbulence is corrected by reverse displacement correction. For the convenience of analysis and calculation, we first use a polynomial to approximate the relationship between pixel intensity and pixel coordinates in a local region of the image. Taking the quadratic polynomial as an example, the intensity distribution of the local area in a speckle image can express approximately as: Where x is the coordinates of a local region in a speckle image, f (x) is the pixel intensity at x, A is a symmetric matrix, b is a displacement vector, and c is a scalar. These parameters can be estimated by the weighted least square method.
Assume that the approximate expression formula of the corresponding region of the reference image RI (This article uses the ensemble average of the speckle images as the initial reference image) is: If the speckle image has an ideal translation d with respect to RI in this region, the relationship between the intensity of this area and the spatial coordinate of the speckle image can be expressed as: By comparing the coefficients of f 1 (x) and f 2 (x), we can get: In the non-singular case of A 1 , we can obtain the translation d through equation (5): So we established a computational model between pixel intensity and displacement under ideal translation. To more accurately describe the projection of the distortion on the image, a reasonable model of the displacement field can be obtained through optimization strategies such as the affine transformation of the displacement field and the weighting of neighborhood pixels (Farnebäck 2003). The specific method is as follows: for the projection of the general three-dimensional motion on the image plane, we use the affine transformation model to express the displacement field as: Written in matrix form: among them, S = 1 x y 0 0 0 x 2 xy 0 0 0 1 x y xy y 2 (11) P = a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 (12) then: Where i is the pixels in the neighborhood, w i is the weights of the pixel i. Neighborhood size needs to determine according to seeing (r 0 ) (which will discuss in Section 4.2). Solving P in equation (13) and substituting it into equation (10), we can obtain the displacement vector d of each pixel relative to RI. Each speckle image can align to RI by inversely shifting and interpolating each pixel according to the d calculated at each time.

Band-pass Filtering
Due to the properties of atmospheric turbulence, the original speckle image is dominated by lowfrequency energy. The displacement field calculated directly from the original speckle image mainly reflects the distortion of the low-frequency structure, which is not conducive to the recovery of the mid and high-frequency phase. Therefore, based on the theoretical speckle interference transfer function (sitf ) (Korff 1973) and the optical transfer function (Otf ) of NVST (yuan Xiang et al. 2016), we constructed a band-pass filter H bp (u) as shown in Eq.14 to filter the original speckle image. This filter enhances mid and high-frequency features and suppresses noise outside the diffraction limit, thereby improving the quality of image alignment.
Where S n is the noise power spectrum estimated from the speckle image sequence, and u is the spatial frequency.

Reference Image Update and Iterative Alignment
Estimating the distortion displacement requires an image with a higher signal-to-noise ratio as the reference image RI (described in 2.1). The initial RI is the ensemble-average image of the speckle images filtered by Eq.14. After the speckle images align to the current RI, the signal-to-noise ratio of their ensemble-averaged image is improved, at this time, NASIR updates the RI with the newer ensembleaveraged image and performs displacement field estimation and alignment again.
After multiple iterations of alignment and updating, the quality of RI and the alignment accuracy of speckle images continuously improved. Then the phase of the ensemble-averaged image of the speckle images is used as the phase estimation of the object image.

Modulus Reconstruction Based on Speckle Interferometry
The basic idea of speckle interferometry is that when the exposure time is less than a certain threshold, the speckle images contain more information than the long-exposure image. The power spectrum of multiple consecutive short exposure speckle images can be time-averaged to restore the power spectrum of the object image (Labeyrie 1970).
In practical observation, speckle images often contain various random noises. Therefore, NASIR uses the Wiener filter method of Eq.15 to reconstruct the modulus of the object image: Where |O(u)| is the modulus of the object image, < |I i (u)| 2 > is the ensemble-average of the power spectrum of the speckle images, S n is the noise power spectrum estimated from the speckle image sequence, sift(u) is calculated based on the seeing (Korff 1973).
Since the distortions in speckle images are corrected, their ensemble-average images contain more mid-and high-frequency energy, which can further improve the contrast of the reconstructed image. So, unlike speckle masking, NASIR uses the aligned speckle image instead of the original speckle image to calculate < |I i (u)| 2 > and then obtains |O(u)| from Eq.(15).

RECONSTRUCTION RESULTS AND ANALYSIS
We compare the reconstruction performance of NASIR with that of correlation Shift-Add and speckle masking by reconstruction experiments on NVST observations. NVST is the largest ground-based vacuum solar telescope in China, and it was put into use in 2010. The NVST has a clear aperture of 985 mm and is mainly used for high-resolution imaging observations of the solar photosphere, chromosphere, and spectrum (Liu et al. 2014). The main instruments of NVST include high-resolution multi-channel imaging systems ) and high-dispersion multi-wavelength spectrometers (Wang et al. 2013).
NVST data products are divided into three levels: Level 0 is the original observation data; Level 1 is based on lucky imaging, corrected by dark current and flat field, and reconstructed by the Shift-Add; Level 1+ is reconstructed by speckle masking , which can produce reconstruction results close to the diffraction limit (Yan et al. 2020).
The experimental data set in this paper is derived from the Level 0 data of the NVST photosphere (TiO, center wavelength 705.8nm, bandwidth 1nm) and chromosphere (Hα, center wavelength 656.28±0.4nm, bandwidth 0.025nm) under different seeing degrees from 2013 to 2019. Each set of experimental data contains 100 short-exposure speckle images (Table 1). The seeing parameter r 0 is calculated from the spectral ratio method(von der Lühe 1984).  Table 1) using correlation Shift-Add, speckle masking, and NASIR, respectively.

Analysis of Reconstruction Quality When Seeing is Good
It can be seen from Fig.1 (granule outlines, bright point and penumbral filaments of sunspot, etc.) that the reconstruction quality of NASIR is significantly better than that of Shift-Add, and it is closer to the reconstruction quality of speckle masking.
To intuitively compare the reconstruction quality in the frequency domain, we convert the 2-D power spectrum image to polar coordinates. Then, the power spectrum image in polar coordinates is averaged along the polar angular to obtain a 1-D power spectrum curve with the length of the polar diameter as the frequency and the average intensity as the energy.
The power spectrum curve in Fig. 2 (left) shows that the resolution of the NASIR reconstructed image is close to that of speckle masking reconstructed image, about 0.15 arc-second, which is significantly higher than that of Shift-Add reconstructed image of about 0.25 arc-second. The intensity profile in Fig.2 (right) further shows the similarity in intensity distribution between NASIR and speckle masking reconstruction. Fig. 2: Power spectrum curves (left) and intensity profile(right, the red line region in Fig. 1 (a))) of images in Fig1. To further reveal the performance of NASIR, we use the following quantitative metrics to compare and analyze the reconstruction results of speckle masking and NASIR: (1) Image correlation coefficient (r): The correlation coefficient is used to compare the degree of linear correlation between reconstructed images, which is defined as: Here, Cov(X, Y ) is the covariance of the pixel intensity of the two images X and Y , V ar[X] and V ar[Y ] are the variances of X and Y respectively.
(2) Structural SIMilarity (SSIM) (Wang et al. 2004): SSIM is an index that measures the similarity of two images. Given two images X and Y , their structural similarity can be defined in the following way: Where µ x is the mean of X, µ y is the mean of Y , σ x is the variance of X, σ y is the variance of Y , and σ xy is the covariance between X and Y . c 1 = (k 1 L), c 2 = (k 2 L) are constants used to maintain stability. L is the range level of pixel values. k 1 = 0.03, k 2 = 0.03. The range of structural similarity is -1 to 1. When the two images are the same, the value of SSIM is equal to 1. The structural similarity index defines structural information from the perspective of image composition as being independent of brightness and contrast, reflecting the properties of the object structure in the scene, and modeling distortion as a combination of three different factors, brightness, contrast, and structure. The mean is used as an estimate of brightness, the standard deviation is used as an estimate of contrast, and covariance is used as a measure of structural similarity.
(3) Coefficient of Variation of Intensity Profile (CVoIP): We use the coefficient of variation of the intensity profile of the reconstructed image, which is the ratio of the standard deviation of the intensity profile to its mean, to measure the difference in image contrast.
A quantitative comparison of speckle masking with NASIR is shown in Table 2, indicating that the reconstructed images of the two methods have high correlation and structural similarity. Therefore, the visual quality (Fig.1), the consistency of the power spectrum curve (Fig.2), and the quantitative comparison index of Table 2 (first row) show that when seeing is good, NASIR can achieve reconstruction quality close to the speckle masking method.
The reconstruction results of the Hα-band data further verify the effectiveness of the NASIR for high-resolution reconstruction of the speckle images when the seeing is better. Fig.3 shows the reconstructed images of the three algorithms for data set 2 in Table 1. The quantitative comparison result is in the second row in Table 2. The power spectrum curves and the intensity profile curves are shown in Fig.4. Seeing   Fig.5 shows the reconstruction results of dataset 3 (r 0 =5.88cm). It can be seen that as the seeing decreases, the recursive cumulative error of the phase of speckle masking has an aggravated influence on the reconstruction, resulting in stray structures in the reconstructed image, while NASIR can better    6 shows the reconstruction results of data set 4 (r 0 =4.09cm) in Table 1. As we can see, with the further decrease in seeing, the correlation coefficient r and SSIM of the images reconstructed by the two algorithms start to decrease significantly (the fourth row of Table 2), indicating a clear difference in the reconstructed images. The phase recurrence error of speckle masking seriously affects the reconstruction quality, the reconstructed image has been destroyed by the pseudo-structure, while NASIR can still obtain low-frequency structures to a certain degree.

Computational Time-consuming
NASIR is a spatially parallel algorithm that avoids the multi-path phase recursion, block reconstruction, and stitching, so its average computation time is about 48.7% of speckle masking.

DISCUSSION
The core idea of NASIR is to recover the phase distortion of different frequency components through non-rigid (non-linear) geometric correction. It is a fast and robust method for high-resolution statistical reconstruction.

Features of NASIR
NASIR performs phase correction through pixel-by-pixel non-rigid alignment, which is the affine correspondence between pixel feature descriptors. Feature descriptors are polynomial fit coefficients of the intensities of pixels and their neighborhoods. So non-rigid alignments are not limited by spatial invariance (isoplanatic angle). In addition, the continuous update of the reference image RI by NASIR can make the low-frequency structure tend to be consistent during the iterative alignment process, so mitigating the geometric distortions between isoplanatic angles. Therefore, NASIR can directly reconstruct the field of view at a time without block-wise reconstruction.
NASIR avoids the computation of high-order statistics and complex multi-path phase recurrence, and the computation time is only about half that of speckle masking.
NASIR estimates the displacement field for each pixel of the field of view in parallel, which can be immune to recursive error accumulation.
Band-pass filtering of the raw speckle images suppresses noise interference beyond the diffraction limit, allowing NASIR to reconstruct the modulus and phase of the mid-frequency structure more efficiently.

Further Improvements to NASIR
While achieving successful experiments, NASIR deserves further optimization and improvement in the following areas: (1) Optimization of displacement field-intensity distribution model. Considering the convenience of analysis and calculation, NASIR adopts polynomial expansion to approximate the intensity distribution of the image area. It may differ from physical models of actual intensity changes caused by atmospheric turbulence. Therefore, establishing a speckle image intensity distribution-displacement representation model more in line with physics is one of the issues worthy of further study.
(2) Adaptive selection of neighborhood size. NASIR uses the intensity distribution model of the neighborhood as a feature descriptor and estimates the displacement through the matching between feature descriptors. A larger neighborhood size can improve the matching stability, but it will produce more smoothing effects and reduce the reconstruction quality of high-frequency structures. Conversely, a smaller neighborhood size can preserve high-frequency details. But in the case of poor vision and high noise, it is easy to cause mismatch and reduce the accuracy of displacement estimation. Therefore, a more appropriate size should be chosen according to the atmospheric coherence length for better alignment performance.
(3) Design and implementation of seeing-based band-pass filters. NASIR uses band-pass filtered speckle images for alignment, essentially phase correction based on the alignment of object structures within the diffraction limit. Therefore, seeing-based adaptive band-pass filter design is also one of the fundamental issues to consider.
(4) Improve alignment and update iteration strategies. Non-rigid alignment based on speckle imaging requires a reference image. Iterative improvement of reference image quality is crucial for phase reconstruction. Empirically, the number of iterations for NASIR to update the reference image is between 3∼5 times, therefore, it is necessary to make more reasonable choices and judgments on iteration and convergence based on seeing r 0 .

CONCLUSION
NASIR is a novel high-resolution reconstruction method for speckle imaging that combines classical statistical reconstruction methods and computer vision techniques. Reconstruction results from groundbased telescope observations demonstrate the potential of NASIR in several ways: Flexibility. NASIR can directly reconstruct the full field of view without sub-block reconstruction and stitching, avoiding stitching errors and discontinuities between stitched sub-blocks.
Robustness. When seeing is good, the reconstruction quality of NASIR is close to speckle masking. When seeing decreases, NASIR is more stable.
Fastness. NASIR does not perform phase recursion and avoids block reconstruction/stitching operations. Computation time is reduced by more than 50% compared to speckle masking.
In conclusion, although NASIR does not have theoretical accuracy guarantees like speckle masking. However, experimental results show that NASIR is still a convenient, robust, and fast method for highresolution solar image reconstruction. We believe that it will play an important role in many applications such as quick look and data filtering.