Blind technique for point spread function equalization with application to the Multi-Order Solar Extreme Ultraviolet Spectrograph

Abstract. When combining remote sensing data from multiple instruments or multiple imaging channels, differences in point spread function (PSF) can lead to systematic error. If the PSFs are not well known, then it is difficult to determine which differences in the image data are meaningful for the object being observed and which are artifacts of PSF. Direct PSF measurements can be problematic. For example, in a sounding rocket payload, launch vibrations and acceleration, subsequent operations in micro gravity, and the impact on return to Earth may all affect PSFs. We have developed a blind method to equalize the PSFs of three distinct instrument channels, as found in the Multi-Order Solar Extreme Ultraviolet Spectrograph (MOSES). To validate our technique, we generate three synthetic images with three different PSFs, with some spectrally interesting features. Thence, we demonstrate the successful removal of PSF-induced artifacts is possible, with the genuine spectral features left intact. We also perform blind PSF equalizations on three copies of the same solar image, but with differing PSFs, after applying independent noise to each. The results accurately reproduce corrections performed in the absence of noise, with full knowledge of the PSFs. Finally, we apply PSF equalization to solar images obtained in the 2006 MOSES flight and demonstrate the removal of artifacts.


Introduction
Quantitative image analysis often requires the combination of data from multiple instruments or instrument channels. Differing point spread functions (PSFs) can lead to systematic error in such analyses. 1 For example, taking the difference between two images of the same scene will highlight any differences in PSF size and shape. When image differences are used for tracking asteroids or finding gravitational microlensing, differing PSFs complicate the search for these phenomena. In particular, image differences on scales smaller than the PSF cannot be trusted if the PSFs of the two images are not the same.
PSF effects and the effects of their differences have been addressed in numerous ways, including numerical deconvolution of measured PSFs and other PSF-equalizing techniques. Direct deconvolution is known to be an ill-posed inverse problem. 2 Direct PSF measurements can also be problematic. A more general equalization technique, without a priori knowledge of the PSFs, is therefore useful. For example, Alard and Lupton 3 handle PSF variations as part of an elaborate image subtraction routine. This routine accounts for subtraction of images with different seeing by assuming each frame in a collection of images can be convolved into the best-seeing "reference image" via a convolution kernel. This kernel is modeled as a series of Gaussians plus x-y polynomials and found via the method of least squares. 3 We present a generic technique for blind PSF equalization, valid for imaging systems viewing objects with equal power spectra or objects, whose power spectra differ in a known way. This PSF equalization scheme was born out of necessity for the interpretation of data from our instrument, the Multi-Order Solar Extreme Ultraviolet Spectrograph (MOSES). 4 MOSES is designed to allow simultaneous imaging and spectroscopy over a large 2-D field-of-view (FoV) in a single ultraviolet spectral line, allowing us to study the dynamics of rapid energy release in the solar atmosphere. Transition region explosive events, for example, are the most clearly and frequently observed examples of magnetic reconnection in astrophysics. 5,6 The signatures of magnetic reconnection are large (∼100 km∕s) doppler shifts and line broadenings, which greatly exceed the transition region sound speed (∼30 km∕s). 7 Fortunately, these shifts are large enough to be clearly seen despite the artifacts we describe in this paper. However, the details of the events, including the spectacular explosive event analyzed by Fox et al., 5 can be modified significantly by PSF artifacts. The goal of the present study is to reliably suppress these artifacts and uncover the genuine structure of the explosive events and other transition region flows. The instrument contains three distinct imaging channels, which we compare numerically to obtain estimates of doppler shift and line width in the solar atmosphere. Unfortunately, each channel has a different PSF, causing different distortions of small-scale features. Left untreated, these differences lead to artifacts that could be misidentified as explosive events or other phenomena. The PSFs define minimum spatial and spectral scales, below which we cannot confidently analyze dynamic events in the solar atmosphere. 5

Instrument Concept
Before describing the PSF equalization procedure, we will briefly describe the application that motivates this work. The MOSES sounding rocket payload reflects three distinct images from a concave diffraction grating onto detectors at the m ¼ þ1;0; −1 spectral orders. Each image contains distinct spectral information across the entire FoV. 4 5 The passband also contains less than 10% Si XI 30.33 nm. The optical characteristics on MOSES are outlined in Table 1.
Key parameters for the MOSES instrument design are tabulated in Table 1. Like most solar instruments, MOSES has a high F-ratio (nearly f∕60). The 20 × 10 arc min FoV is unvignetted, with the aperture stop at the primary mirror, which is the concave diffraction grating of Fig. 1. Consequently, the PSF can be expected to be uniform across the FoV. The detectors are rear illuminated CCDs, 2048 × 1024, with the long axis oriented parallel to the grating dispersion. They are placed at the m ¼ −1;0; þ1 orders of the grating. The grating groove profile is lamellar, with the groove depth set to illuminate these three orders in a 1∶2∶1 ratio. Almost no light is diffracted into higher orders. 8 A flat secondary mirror (not shown in Fig. 1) folds the optical system in half. Identical multilayer coatings on both optics give a narrow passband dominated by He II 30.4 nm. Thin film aluminum filters, placed about 200 mm from the focal plane, prevent visible and UV radiation from reaching the detectors. The instrument is described in greater detail by Fox. 10 The primary mirror is in full sunlight, and the secondary is also brightly illuminated, when the instrument shutter opens for exposures during our 5-min observing period on a ballistic trajectory from 160 to 300 km and back. Consequently, the grating and its mount are subject to thermal transients. Despite the low expansion substrates (Zerodur primary, ULE secondary) and semikinematic mount designs, there is evidence of a timevarying PSF during flight. 9 Consequently, ground testing is of little use in combating the PSF effects we address in this paper. Figure 1 displays a simplified conceptual representation of MOSES imaging monochromatic letters "A" and "B," overlapping in the object. In the m ¼ þ1; −1 orders, the "redder," (longer wavelength) light is dispersed farther out from the optical centerline than the "bluer" (shorter wavelength) light by a distance proportional to their difference in wavelength. The central order is a pure imager, with the letters overlapping.
A set of simultaneous images from each order can be viewed as a tomographic projection of a hyperspectral object with two spatial dimensions x and y and one spectral dimension λ. We define x as the axis along which dispersion occurs. An object υðx; y; λÞ forms images at spectral orders m, given as follows: (1) where υ is the object to be imaged, λ is measured from the He II line center in pixels at 29 mÅ∕pixel, and B is the instrument passband. 4 Since dispersion is along the x-axis, the intensity in a given pixel ðx; yÞ is a projection through the x-λ plane at slope. The y coordinate is completely independent of λ and x. This projection is illustrated in Fig. 2, which is situated in a plane of constant y. At the top of the figure, three bright sources are shown at top in the x-λ plane. Their spectra are, respectively, red shifted, broadened, and blue shifted. Their projections onto the coordinates of the three detectors are arrayed along the bottom edge of the figure.
A unique PSF is associated with each order. The systematic error induced, in particular in the m ¼ þ1; −1 orders, is noticeable on small scales. The difficulty is not that there are aberrations, but that they vary in each order.
Our PSF equalization scheme, described in Sec. 3, is intended to solve this issue.

Motivation
The outboard order detectors in MOSES have PSFs with characteristic, oppositely tilted shapes. These PSF shapes appear prominently in unresolved features, that is, features that are similar in scale to the PSF or smaller. Courrier and Kankelborg 11 describe in detail how Doppler shift estimates from MOSES data are influenced by the differing PSFs in our three spectral orders. In the m ¼ −1 order, red shifts move rightward and blue move leftward, and in the m ¼ þ1 order, red moves left while blue moves right. Therefore, these tilts give false apparent Doppler information as if the entire sun were made of small, weak bipolar jets, with a blue-shifted south jet and red-shifted north jet. In the difference image below, subtracting the m ¼ þ1 image from the m ¼ −1 image produces a distinct quadrupole, which is the clearest sign of a bipolar jet (Figs. 3 and 4). This effect has limited MOSES analysis to features whose spectral variations are very large compared to systematic error. 5

Point Spread Function Equalization
In each order m, the observed image I 0 m can be understood as a convolution of the "true" image I m with each order's unique PSF, κ m . In Fourier space, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 3 2 6 ; 5 1 5Ĩ 0 m ¼κ mĨm ; (2) where the tilde represents the 2-D Fourier transform with respect to x and y. In this paper, we use integral representation for compactness. In practice, our Fourier transforms and convolutions are carried out discretely on digital data. 12 Instruments on sounding rockets are subject to vibrations during launch and subsequently operate in microgravity. Ground-based measurements of the PSFs are therefore unverifiable. Rather than attempt to ascertain each PSF, we calculate a PSF-compensated image inferred only from flight data. For each order m, we define E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 3 2 6 ; 3 8 2Ĩ 00 m ¼κ 0Ĩ m : We call κ 0 the intermediate kernel. Each image I 00 m can now be directly compared to the others, as they now share a common PSF. For the intermediate kernel, we choose the Fourier-space geometric mean of all the kernels: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 3 2 6 ; 3 0 6κ  Journal of Astronomical Telescopes, Instruments, and Systems 028002-3 Apr-Jun 2018 • Vol. 4 (2) where N ¼ 3 is the number of orders, each with a unique, however, unknown, PSF. Before we proceed further, we must consider the effect of spectral line profiles, which blur each MOSES spectral order differently. The line profiles are not instrumental artifacts but are properties of the object, and we do not wish to remove their effects.
The He II line profile is not constant over space but has a well-defined average. 13 Let us consider the special case in which the spatial and spectral components of the object υ are separable: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 6 3 ; 6 4 2 υðx; y; λÞ ¼ fðx; yÞPðλÞ: If the above equation was actually correct, it would mean that the spectrum is the same in every pixel. While we do not expect this to be true, there is a good reason to believe that the average spectrum taken over any sufficiently large subregion of the image would be nearly the same as a similar subregion. More specifically, when averaging the He II 30.4 nm line profile over the SERTS spectrograph slit, Andretta et al. 13 measure Gaussian, equal-width line profiles in both quiet and active solar regions.
Without loss of generality, we assume ∫ B PðλÞdλ ¼ 1. Separability allows us to cast Eq. (1) in the form: The m ¼ 1 image is then a convolution of the intensity image with the average line profile, P: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 8 ; 6 3 ; 3 7 0 I 1 ¼ The m ¼ −1 case is a convolution with the mirror image of the average line profile: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 9 ; 6 3 ; 3 0 5 Now, the average line profile itself is an intrinsic quantity of the object, but as we see in Eqs. (7)-(9), the way in which this profile affects each channel of our instrument varies. Therefore, it is helpful to define an effective profile for each order: Here, recall that λ is measured from line center. Now, we see Eqs. (7)-(9) can be rewritten as I m ¼ f Ã P m . Consequently, E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 0 ; 3 2 6 ; 7 5 2 jĨ m j jĨ n j ¼ jP m j jP n j : Equation (11) suggests that the power spectra of the images taken in the different orders of the instrument are the same except for a contribution from the line profile. This was derived by assuming a constant line profile everywhere in the image [Eq. (5)]. While we know Eq. (5) is not strictly true, our purpose has been to provide a plausible motivation for Eq. (11), which is the major assumption behind our PSF equalization scheme. The concern about line profiles is unique to MOSES; in a simpler application, such as the PSF equalization of images to find asteroids or microlensing events, Eq. (11) would reduce to an equality between the images' power spectra. Now, we generate a PSF equalization filter C m for each order from only observed quantities: Under the assumption that the power spectra vary only due to the known mean line profile [Eq. (11)], the result applying the filter [Eqs. (11) and (12)] is the same as the results of Eqs. (2) and (3), as can be verified by direct substitution. The mean profile can be estimated from past observations. As mentioned earlier, for He II 30.4 nm, the SERTS spectrograph found very similar average line profiles in an active region and a quiet region. The intent of our filter is to modify the spectra of the images such that the PSFs are equal while leaving the signature of the expected average line profile in the data. Since our filter is defined by looking over a large image with many individual line profiles whose distinct line profiles average to P, we hypothesize that local deviations of the line profile from average will be unaffected by the filter.
The filter described by Eq. (13) is, in principle, neither a sharpening nor a smoothing filter. It aims, rather, to merge the transfer functions of the three images to their geometric mean. Nevertheless, at select high frequencies, where any of the denominator terms contain very small values, the filter can contain large values or floating point exceptions. We replace nonfinite values with zero and apply the following regularization to every element of the filter in Fourier space: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 3 2 6 ; 2 5 7 The regularization maintains the phase of each element of the filter but modifies large magnitudes according to a smooth transfer curve (Note that phase could, in principle, enter into the filter through an asymmetric mean line profile). Magnitudes below a threshold of T ¼ 2 are minimally affected, but none of the magnitudes is allowed to exceed 2, and large magnitudes are zeroed.
In the inverse transform of Eq. (13), C m in spatial coordinates would appear as a small kernel with significant positive and negative contributions in a neighborhood a little larger than the instrument PSFs. Far from this useful kernel, low-level noise extends across the entire image domain. Upon convolution, this noise creates artifacts. Therefore, we cut off the noisy void in our filter arrays, in coordinate space, with the following window: where r is the radial coordinate from the center of the kernel in pixel units. Note that the filter in Eq. (13) is constructed without knowledge of the instrument PSFs.

Numerical Verification
As a first test of our technique, we generate a set of three synthetic images corresponding to MOSES' three instrument channels. As we will discuss in Sec. 5, small-scale artifacts lead us to believe that the outboard MOSES channels have some form of tilted, oppositely oriented PSFs.
To begin, we create three synthetic images with an identical, random field of spots. Since our software infrastructure demands a line profile input, we insert a delta function to convolve with the outboard orders. In view of the relationship between the scene and the MOSES spectral orders [Eq. (6)], this choice guarantees each image remains identical thus far. In the lower left hand corner, we exclude the random spots and place a test pattern consisting of three spots: an unshifted spot, a spot blue-shifted by one pixel, and a two-component event with both red and blue shifts. To imitate MOSES-like artifacts, we convolve each image with the PSFs shown in Fig. 5. We then feed these images to our algorithm and inspect the reduction of the artifacts. Since we have designed all spots (other than the spectral features in the corner) to be identical in each channel prior to convolving the images with our PSFs, the test will be deemed successful if we see mostly empty difference images between channels after correction. However, the blue-shifted spot and the two-component spot in the lower left-hand corner should show appropriate signatures in their difference images.
These PSFs (in particular, the PSFs associated with the outboard orders) blur small-scale details in a way that imitates observed MOSES artifacts. Intensity is blurred upward and "outward" in the sense of Fig. 1, and downward and inward, in the outboard orders, and broadened horizontally in the m ¼ 0 order. The resultant images are shown in Fig. 6.
When studying MOSES-like images, it is instructive to examine difference images. For example, let us consider subtracting the m ¼ 0 order from the m ¼ þ1 in Fig. 1. In the difference image, light from the top of the red "A" and bottom of the blue "B" will stand out brightly while the overlapping core will dim. Similarly, where red and blue shifts occur in our data, image subtraction will highlight brightened spots that have moved outward or inward, respectively, relative to the position of the corresponding spot in the m ¼ 0 order. Let us examine a set of difference images of our synthetic data, prior to PSF equalization Fig. 7.
The true spectral features we have built in appear in the lower left hand corner Figs. 6 and 7. The rest of the images are filled with artifacts due only to PSF differences. Note that in the leftmost image, the artifacts appear as quadrupoles with the same orientation. Let us now look at our PSF-corrected images in Fig. 8.   The systematic tilting visible in the upper panels of Fig. 8 appears to have been eliminated by our PSF equalization. To ensure that our genuine spectral features remain intact, we must examine our difference images after correction (Fig. 9).
Evidently, our technique is effective in removing artifacts due to PSF variation. Samples taken in areas excluding our spectral features indicate the artifact variance is suppressed by a factor of about 230.

Effects of Noise
The theory of Sec. 3 was developed without considering the presence of noise in the data. The PSF-corrected images I 00 m are calculated in two steps [Eqs. (11) and (12)]. First, the correction filterC m is calculated from products of amplitude spectra of the noisy dataĨ 0 m ,Ĩ 0 n in the Fourier domain. Then,C m is convolved with the noisyĨ 0 m . The effects of all this noise, especially in the creation ofC m , are difficult to anticipate.
In this section, we test the robustness of the calculation against noise by comparing calculations on synthetic data, performed with and without noise. We generate realistic, PSF-corrected images, with noise, by the following procedure.
1. Arbitrarily choose a single image from a single MOSES exposure; in this case the m ¼ 0 image from the 24th exposure (Fig. 10).   2. Fill in bad/saturated pixels with their surrounding finite pixels via a smoothing convolution kernel.
3. Scale the image to a given mean photon count.
4. Convolve three copies of this image with three different PSFs.
5. Apply a Hanning window to the images to minimize edge effects. 12 6. Apply Poisson noise independently to each of the three images.
7. Generate filters from the noisy images [Eq. (13)], using a delta function for the line profile.
For this numerical experiment, we have chosen the same PSFs shown in Fig. 5 and described in Sec. 4. The symmetry of the PSFs used in this example, and the fact that we generate the three spectral orders by direct convolution with the same source image, guarantees that filters constructed according to the recipe in Eq. (13) would exactly correct the images in the absence of noise. Note that read noise has been omitted from our consideration because it is below the shot noise for our cameras.
The PSF corrections performed by the above procedure are based solely on noisy data, with no prior knowledge of the PSFs. To assess the effect of noise, we compare all three of the PSFequalized images against a single ideal image, generated by a "cheated" technique in which we explicitly generate the intermediate kernel [Eq. (4)] from the known PSFs and convolve this with the original image from Fig. 10. We take this "cheated" image, produced by direct convolution in the absence of noise, as the ideal standard for the PSF-corrected image.
We next generate a reduced chi-squared statistic for each order as a measure of fidelity, comparing the results of the PSF equalization routine to the "cheated" I 00 . A plot is shown in Fig. 11. The technique produces generally satisfactory results for photon counts below 100,000. At high photon counts, the chi squared values transition to a power law.
Beyond this point, our signal-to-noise ratio (SNR; Fig. 12) no longer improves with increased photon count. Evidently, we are limited to an SNR of 400. We suspect that this results from roundoff error propagating through the filter generation and convolution operations. In practice, the systematic error due to PSF reconstruction is well below shot noise. At full well exposure, MOSES gets ∼8000 photons∕pixel (Table 1), leading to a maximum possible instrument SNR of ∼90. Thus, we see two regimes. In the low signal regime, the reconstruction noise is identical to what we would expect if the photon counting noise had simply been added, after the fact, to an ideal PSF-compensated image. In the high signal regime, we find a limiting SNR due to the synthetic error in PSF compensation, but this SNR limit is greater than what could be realized by the instrument. The PSF compensation is therefore robust against all the considered sources of noise.  5 These filters will be applied to the entire dataset.

Quantitative Reduction of False Doppler Shifts
To quantitatively examine reduction of the artificial Doppler shift, we examine the small feature from Fig. 3 (Figs. 13 and  14). An additional example of PSF-derived artifacts is shown in the left panel of Fig. 15. The two quadrupolar features are similar to those seen in our synthetic data (Fig. 9) and to the example in Fig. 14. In fact, similar features appear in many places across the FoV in the 2006 MOSES flight data. If these were genuine spectral features, they would not all have the same orientation. After correction, the artifacts are strongly suppressed.
To quantify the difference in Doppler shifts, we employ a program that background subtracts around the feature and takes the centroid of each row in the array from bottom to top, finding the "center of brightness" (CoB). The program then calculates the difference of CoB between the plus and zero orders to infer the false Doppler shift and repeats this procedure with the PSF-equalized data. A CoB-shift of one pixel corresponds to a Doppler shift of 29 km∕s. To aid the eye in interpreting the plot (Fig. 16), we calculate the integrated intensity of each row, going from bottom to top, and place it directly under the plot of the Doppler shift. Vertical dashed lines mark off the FWHM of the intensity distribution, to illustrate the section of the window where we have an event, rather than background.  The pattern of positive-to-negative Doppler shift, which we identify as an artifact, has been reduced to well below the local sound speed (∼28 km∕s) by applying the filter. A smaller residual blue shift remains. This may represent a real, physical blue shift; but we cannot categorically eliminate the possibility that some artifact remains.

Preservation of True Solar Features
It is important that our correction technique corrects only artifacts while leaving true solar features intact. To assess this, we examine a blue-shifted jet with difference images from the m ¼ þ1; −1 orders, before and after processing. Figure 17 shows the difference between the m ¼ −1 and m ¼ þ1 images of a blue shifted jet, before and after PSF equalization. While the obvious artifacts in the MOSES data are suppressed, genuine spectral features remain largely unaffected.
To quantify preservation of Doppler shifts in observing true solar features, we also apply our CoB finding technique (see the previous section) to the bipolar explosive event described by Fox et al. 5 The event in m ¼ 0 is shown in Fig. 18. Difference images, m ¼ þ1 minus m ¼ 0, are shown in Fig. 19. The bright kernel at the center of Fig. 18 shows up as a dark center with bright wings to the left and right in both of the difference images. This is a clear signature of line broadening. The jets to the north and south show bright-dark and dark-bright patterns in the difference images, indicating red and blue shifts, respectively. After PSF correction (right panel of Fig. 19), the   qualitative Doppler signature of the explosive event remains intact. However, the CoB plot (Fig. 20) shows a less extreme redshift in the lower jet, indicating that some of the apparent Doppler shifts were likely caused by PSF variation (compare Fig. 18 in Fox et al. 5 ).

Conclusions and Future Work
We have developed a technique that permits a series of images of the same scene, but with different PSFs, to be filtered so that they have the same PSF. This PSF equalization is a blind technique in the sense that no prior knowledge of the image PSFs is assumed. We have developed this for application to images taken by the MOSES objective grating spectrograph, which forms images of the solar transition region in the He II   Journal of Astronomical Telescopes, Instruments, and Systems 028002-9 Apr-Jun 2018 • Vol. 4 (2) 304 A line at multiple orders. The technique also leaves room for a known, prior, systematic difference in the image power spectra due to the average spectral line profile of the object. We have validated our PSF equalization scheme numerically on synthetic data. Difference images generated in this test clearly illustrated the types of artifacts generated by differing PSFs, which we see repeated in cookie-cutter fashion across the MOSES images. The test with synthetic data demonstrates that those artifacts can be suppressed by a large factor (image variance is reduced by a factor of more than 200). A second test on images derived from solar data shows that the method is practically immune to noise.
Next, we applied the method to actual MOSES data. Artifacts that stand out in interorder image differences are visibly suppressed, whereas true solar features of size scales larger than the artifacts appear to retain their essential Doppler and spectral characteristics. Our PSF equalization filter adjusts only the image power spectra and not the phases. This implicitly assumes symmetric PSFs. While estimates of the MOSES PSFs indicate slight asymmetry, 9 the visual results of PSF equalization of MOSES data look promising.
We still need to verify that spatially varying line profiles (the very thing we wish to observe) do not unduly influence the generation of PSF equalization filters. Other possible concerns include the presence of contaminant lines in the data. Hyperspectral data that could be used to generate plausible synthetic datasets including such effects are available from a number of sources, including the NASA IRIS mission. 14 The method we have developed has other potential applications. Interimage PSF differences analogous to those encountered in the MOSES data also afflict datasets used for multiobservatory astrophysical observations, 1 studies of gravitational microlensing, 3 and asteroid finding. 15