The non-Fourier image reconstruction method for the STIX instrument

Abstract In this work we aimed to develop the image reconstruction algorithm without any analytical simplifications and restrictions. In our method we abandon Fourier’s approach to image reconstruction, and instead use the number of counts recorded in each detector pixel, and then reconstruct each image using a classical Richardson-Lucy algorithm. Among similar works performed in the past, our approach is based, for the first time, on the real geometry of STIX. We made a preliminary analysis of expected differences in STIX imaging which may occur due to usage of slightly different geometries. The other difference is that we use single-pixel-response maps. Namely, knowing the instrument geometry we are able to calculate the detector response for point sources covering entire the solar disc. Next, we iteratively combine them with varying weights until the best match between reconstructed and observed detector responses is achieved. Preliminary tests revealed that the developed algorithm reproduces high quality images. The algorithm is moderately fast, but the result comparable to CLEAN algorithm is obtained within 20-50 iteration steps which takes less than 2 seconds on typical portable computer configuration. The location, size and intensity of reconstructed sources are very close to simulated ones. Therefore the algorithm is very well suited for the detailed photometry of the solar HXR sources. Moreover, its simplicity allows to improve photon transmission calculation in case of any grids uncertainties measured after the launch.


Introduction
The idea of HXR telescopes using Fourier techniques was given by Oda (1965). The rst observation with such equipment was a rocket experiment which allowed the X-ray imaging of the Crab nebula (Oda et al. 1967). Subsequently, a number of Fourier imagers were constructed (Takakura et al. 1983;Kosugi et al. 1991;Sylwester et al. 2000;Lin et al. 2002).
The main characteristic of Fourier imagers is that the image is not obtained directly. Rather, the imager consists of pairs of grids with various pitches and orientations. Each pair of grids modulates the incoming radiation depending on source size and location. Having sets of intensities measured with dozens of grid pairs, we can reconstruct the spatial distribution of X-ray emission. There are many algorithms developed like those adopted from radio interferometry: BackProjection (Mertz et al. 1986) and CLEAN (Högbom 1974). Other algorithms include Maximum Entropy Method (MEM) (Sato et al. 1999), PIXON (Puetter 1995), Forward Fit (Aschwanden et al. 1999), and Expectation Maximization (EM) (Benvenuto et al. 2013). Each algorithm has its strengths and limitations. There is no one, universal algorithm that is robust for all source con gurations, count statistics etc.
The Spectrometer Telescope for Imaging X-rays (STIX) (Krucker et al. 2020) is one of the instruments installed on board the Solar Orbiter mission (Müller et al. 2013). It will record solar are spectra and images in the energy range 4 − 150 keV. The spatial distribution of hard X-ray (HXR) emission will be found from 30 pairs of grids. Images will be reconstructed with a number of existing algorithms as well as ones developed recently (Duval-Poo et al. 2018;Massa et al 2019). However, at present no algorithm exists which is based on real geometry of STIX telescope.
We decided to develop the image reconstruction algorithm without any analytical simpli cations and restrictions and using parameters for ight model of STIX. The price to pay is computational time, but we are relying on Moore's law which states there will always be an increasing of computing power. This enables us to resolve old problems in more straightforward manner solving numerically basic equations which cannot be solved analytically.
Our algorithm uses transmission functions calculated for single sources smaller than the STIX spatial resolution, giving detector count rates. This is the most straightforwad approach to image reconstruction. At present it allows for obtaining one image in a reasonable time interval of several seconds. In the following sections we brie y describe STIX telescope (Section 2), the grid transmission (Section 3), the reconstruction algorithm (Section 4), the algorithm performance for some typical source distributions (Section 5), and discussion of results (Section 6). STIX (Krucker et al. 2020) is one of the instruments installed on board Solar Orbiter mission (Müller et al. 2013). It will record solar HXR spectra in the energy range 4 − 150 keV divided into 32 bands with widths changing from 1 to 15 keV. STIX will be capable of providing information about the spatial distribution of HXR emission thanks to the system of grid pairs mounted in front of detectors.  Figure 1 is a sketch of the instrument. STIX consists of three main blocks. The rst is the X-Ray window mounted within the thermal shield of the Solar Orbiter. It is responsible for heat shielding and rejecting X-Ray radiation below 4 keV. The next block is system of front and rear grids which constitute the collimator system. Just behind the rear grids the detectors are located together with its electronic module.

STIX instrument
The imaging concept is based on pairs of grids which were previously used in several instruments. The example of grid used in STIX is shown in Figure 2. Their characteristic is that slits and slats have di erent widths and they are supported by perpendicular bars (bridges). In previous experiments the front and rear grids were shifted in phase while slats were parallel. In STIX, the front and rear grids are slightly tilted, producing a characteristic Moiré pattern in the detector plane. The pattern obtained via tilted grids is very sensitive to incoming photon directions. To utilize this sensitivity a new type of pixelized detectors, Caliste-SO (Meuris et al. 2015), is used. The pixelated pattern is shown in Figure 3 (left panel). There are 12 pixels: 8 large and four small. They are grouped in four stripes. Stripe B is marked with gray in Figure 3 (left panel). The basic observables are four values measured in stripes A, B, C, and D. The orientation of detector with respect to the Moiré pattern is chosen in such a way that one full period of the pattern is covered with detectors front window (Figure 3,right panel). Such a con guration allows for precise measure of the Moiré pattern phase shift. Moreover, it allows for calibration of the detector orientation and Moiré pattern through comparison of the signal measured in two large pixels from the same stripe. In the case of very strong are, some of the pixels may be switched o to decrease the detector dead time. The whole system should provide us with data for ares of various GOES classes (from A1 to X20) with the angular resolution reaching 7 ′′ . Taking into account that the distance to the Sun will be as small as 0.3 AU (perihelion distance) we expect to achieve unprecedented spatial resolution of about 1500 km for HXR images, which is almost 3 times better than spatial resolution achieved by previous experiments.

STIX grid transmission and image formation
A single grid consists of many parallel tungsten bars forming a system of slits with a width d ( Figure 2). The spatial period of the grid p is called 'pitch' and is equal to the sum of the width of the slat and the slit ( Figure 4, right panel). As slats are opaque to incoming radiation (at least in the rst approximation, for photons with energies below 30 − 40 keV) we can describe the single grid transmission with only two values: 0 and 1. Figure 4 (left panel) presents Monte Carlo simulations of grid transmission, where each recorded photon is presented with the black dot. In our convention, the black color in Figure 4 (left panel) corresponds to 1 (full transparency) and white corresponds to 0 (no transmission), thus the black stripes in Figure 4 represent slits. One dimensional cross-section of the grid along the x 1 axis can be described by step function with values of 0 and 1, in the form presented in Figure 4 (right panel). This is a periodic function with period p, known also as a pulse wave, pulse train or rectangular wave. Periodic function can be easily expanded into Fourier series, i.e. a sum of the sine and cosine functions in the form: The origin of the coordinate system can be selected freely, but we put it in the middle of the slit (see Figure 4, right panel). Therefore our function is even. It is f (x) = f (−x), and thus all coe cients bn in the Equation 1 are equal to 0. The a 0 is an average function value over period equal to Ad/p where A is the amplitude. In our case A = 1, so a 0 = d/p. For rectangular wave with pulse width d, the coe cients an are equal (︀ 2/πn )︀ sin (︀ πnd/p )︀ and Equation 1 reduces to: Transformation of 1D pulse train to 2D grid is made by transformation: f (x 1 , x 2 ) = f (x 1 , 0) = f (x 1 ) as explained in Figure 5, therefore: To describe tilted 2D grid we rotate coordinate system from the one with slits parallel to the x 2 axis to the one parallel to x ′ 2 axis as is shown in Figure 6. Corresponding transformations are x ′ 1 = x 1 cos α − x 2 sin α and x ′ 2 = x 1 sin α + x 2 cos α respectively and angle α is calculated clockwise from x 2 axis.
Using Equation 3 for the grid in coordinate system The result of the transformation (Equation 4) applied to the grid with vertical slits from Figure 4 (left panel) and α = 25 ∘ is presented in Figure 6.
Using wave vector we can write Equation 4 in more simple way: There is a di erence between above equations and the STIX grid transmission used in previous works (Giordano et al. 2015;Massa et al 2019). The reason is because we use actual rectangular grids instead of quadratic ones. In the case of a quadratic train the width of slits and slats are equal and d = p/2 therefore d/p = 1/2. This approach is more analytically convenient while may produce signi cant artifacts in reconstructed images which will be discussed further in this Paper. Our grids parametrisation is based on their actual geometry which is not square but rectangular (Krucker et al. 2020).
Having in mind above di erences we can follow the transmission calculations. Each grid pair is described by their pitches p f and pr for front and rear grids respectively. For a given pair, grids have identical slit widths d but slightly di erent slit orientation α f and αr, which results in di erent wave vectors ⃗ k f and ⃗ kr. We introduce four parallel 2D coordinate systems dened by the Sun surface (⃗ x), the two grids (front ⃗ z and rear ⃗ q), and the detector surface ⃗ y, respectively (see Figure 2.2 in Giordano et al. (2015)). The distance of Solar Orbiter front grid from the Sun is S. The separation between front and rear grids is L 1 = 55cm << S. The distance between the rear grids and the detectors is L 2 = 4.7cm << S. Using similar triangles, we can nd a relation between four coordinate systems (Giordano et al. 2015, equations 3.1 and 3.2) and then express the grid transmission as a function of ⃗ x and ⃗ y: where f , r and tr are front, rear and total transmission functions for given STIX subcollimator. Giordano et al. (2015) introduced in their equations the constant vector ⃗ t = (0, t) which accounts for a possible y 2 -translation of the grid pair with respect to the origin of the coordinate system. In Equation 6 and 7, we introduce similar vector ⃗ w, but we use it mainly to the grid translation along y 1 . For example |w| = p 2 means shift of the origin of the coordinate system from center of the slit to the center of the slat. In addition, we calculate transmission functions for bridges, which support nest grids, in a similar way i.e. via Fourier expansion. Therefore, the total transmission of the subcollimator is: where , and rb are front and rear bridges respectively. The bridges should not degrade images (Krucker et al. 2020), but we included them for completeness of the instrument geometry calculation. We calculated tr . We used nite number of harmonics (Equation 6, 7), usually 200. That is enough since we assumed simple grid transmission with only two values: 0 and 1. In order to achieve this, the following approximation was used: if f ≥ 0.5 then f = 1, and if f < 0.5 then f = 0. However, for the entire grid transmission calculations we did not use regular array with equidistant (x i , y i ). Instead 10 5 randomly distributed (x i , y i ) pairs were used. Figure 7 presents the grid pair No. 32 transmission pattern (modulation pattern) on the detector front face, for the point source located on the center of the Sun (left front) and for the point source located at heliographic coordinates 60" W, 60" N (right panel). Grids parameters are: p f = 682.1973 µm, pr = 649.830 µm, α = 51 ∘ .702, α = 48 ∘ .379, and d = 336 µm. The detector has a size 8.8 mm × 9.2 mm and is divided into four large stripes denoted A,B,C and D. The 'observables', i.e. values measured by STIX, are the countrates recorded in these four stripes.

Image restoration with RL algorithm
In a standard approach, the uxes measured in STIX stripes A, B, C and D on each detector are converted into visibilities, i.e. amplitudes and phases in Fourier space (Krucker et al. 2020). These visibilities are then used to reconstruct image with a use of several numerical algorithms. In this paper we omit the Fourier approach and limit ourselves only to the uxes measured at each detector stripe. We look for distribution of the emission source which matches the best measured uxes.
For this purpose we have adopted well known, classical Richardson-Lucy (RL) algorithm based on Bayes theorem of the statistical distribution of the detection noise, and used extensively in astrophysical imaging. Richardson (1972) and Lucy (1974) gave independently formulas for iterative deconvolution of images blurred by instrumental Point Spread Function (PSF). In this case, the image element (pixel) u j is smeared by PSF array p ij to the observed image element o i = ∑︀ j p ij u j . Looking for the maximum likelihood distribution of u j that gives o i under known p ij they found that next, k + 1 approximation to u j can be obtained from previous solution k as: j is the image value calculated from k-th solution. In other words, k + 1 solution is obtained from previous one by multiplying it by the sum of the observed to calculated ratio weighted by p ij . The summation is done over all image elements i which contributed to the observed image element j. Sheep and Vardi (1982) showed that above iteration process converges to the maximum likelihood solution for Poisson statistics in the data. In general the total sum of weighted functions is not equal 1 as in case of PSF (︀ p ij )︀ , and we should use normalization instead: Even if, in various applications, the method does not meet entire Bayesian assumptions, its idea is very simple and intuitive. If the observed values are greater than calculated ones, we should increase all elements that give a contribution to the data, proportionally to these contributions. If the observed values are less than calculated ones, we should decrease all elements correspondingly. A nal correction to given element at the next k +1 step of iteration is the sum of observed to calculated data ratios weighted by corresponding contributions. This is the simplest method and anyone can experience this marvelous algorithm using for example max_likelihood.pro function from IDL Solar Soft Ware (SSW) library.
There are many examples of the application of the iterative RL algorithms in 1, 2 and 3 dimensions. In the 1D case we can use it to calculate the di erential emission measure distribution (DEM) as function of temperature from set of uxes observed in X-ray emission spectral lines (Withbroe 1975;Sylwester et al. 1980;Kępa 2008). A new approximation of the DEM distribution is obtained from the previous one using observed to calculated uxes ratio, weighted by the set of temperature dependent emission function of the lines used. In 2D case we modi ed and accelerated RL method for deconvolution of soft X-ray images from the Yohkoh/SXT telescope (Sylwester and Sylwester 1998). The weighted function used is the PSF of the telescope. In 3D case we deconvolved X-ray light curves of eclipsing binaries into 3D distribution of emitting X-ray stellar coronae (Siarkowski 1992;Siarkowski et al. 1996;Preś et al. 1997). New iteration of 3D corona is obtained from previous one using observed to calculated ux for given binary phase angle, weighted by contribution of those coronal elements which are not occulted by both stars.
Using Equation 8 and Figure 7 we can describe transmission of any STIX grid as transformation of emitting point (x1, x 2) at solar surface to four uxes measured at four detector stripes. The spatial distribution of emission sources is represented in the form of discrete image im(i,j). Therefore the ux F cal at pixel A, B, C and D is:  N, N) and we have to normalize it before using in RL algorithm. We de ne w (m) = ∑︀ i,j tr1(m, i, j)/N 2 and w1 (m, i, j) = tr1 (m, i, j) /w (m). An iterative procedure for "images" observed by STIX (Equation 9) can be now expressed as: where F obs (m) and F cal (

Algorithm performance
We started from a simple test of the RL method assuming 1 × 1 arcsec pixel on the Sun. The image was 601 × 601 array covering 10 × 10 arcmin area. The source was a single Gaussian ellipse (Figure 8, top-left panel). The counts measured in each of 120 stripes is presented in the right panel of Figure 8. The STIX transmission function has the form tr1 (120, 601, 601). We calculated STIX grids transmission using Equations 6 to 8 and 10 for each detector (30), each subpixel (4), and for each single-pixel source of 601 × 601 map. All these transmissions were folded into nal STIX "PSF" -tr1(120, 601, 601) array. Sources were simulated with Poisson noise included. As a rst approximation we assumed a at image with im 0 (i, j) = const ≈ ∑︀ m F obs (m) /N 2 and run iteration as described by Equation 12. Figure 8 shows the results of iterations after 10, 50 and 500 steps. The presented maps are cuts (65 × 65 pixels) from the entire reconstructed image. After 500 steps no background outside this area was detected. Mean computation time for 100 steps was 2.5 s on the average laptop con guration (i5-6300U 2.40 GHz CPU, 8 GB RAM). At this point we do not have any parameter that stops the iteration process. The result is presented in Figure 8. The source ts best the simulated one after 500 iterations steps. We conducted reconstruction until 10000 steps, and observed small degradation of resulting image. To avoid this problem in further tests we assumed that the iteration is stopped when the change of χ 2 is less than 10 −4 .
Next, we performed tests on various emission sources con gurations, and compared results with CLEAN-Vis algorithm existing in the SSW tree. Assumed total counts are related to moderate signals (less than 10000 counts/detector).
The exception is low signal case where 1000 counts in total were simulated. The testing con gurations were: -Two Gaussians, dynamic range 3 : 5, 130000 counts, -Two Gaussians, higher dynamic range 1 : 5, 100000 counts, -Two Gaussians, low signal -the total number of counts recorded by all 30 detectors was less than 1000, -Faint structure, small sources located close to each other, each source is smaller than instrument angular resolution, 230000 counts, -Typical aring geometry: two foot points, loop-top source, above the loop-top source, 300000 counts. Figure 9 presents simulated con gurations ( rst column), images reconstructed with CLEAN-Vis algorithm (middle column), and results of reconstruction made with RL algorithm (right column). The images reconstructed with both algorithms show the same number of sources. The exception is a ne structure, but in that case size of model sources were below the STIX angular resolution. For each model and reconstructed source de ned by a 50% isophote we tted ellipse. From the t we took following parameters for a comparison: -center coordinates, -diameter described by semi-major and semi-minor axes of the ellipse t to the source, -brightness which is a source ux measured in the 50% of the maximum isophote and normalized to the total ux.
The obtained parameters are collected in the Table 1. For a better understanding of the results we present them in three plots ( Figure 10). Each plot presents comparison of the model source parameters obtained with CLEAN-Vis algorithm (red triangles), and RL algorithm (blue asterisks). The Figure 10. The plots present parameters collected in the Table 1. Red triangles represents results from CLEAN-Vis, and blue asterisks are for RL algorithm (see text for further details). equality of parameters is represented by a dashed line. The top panel contains both x and y, coordinates of the source center. It is clear that the locations are reconstructed with very good accuracy. Usually reconstructed sources were shifted not farther than 1 pixel from simulated source. We did not observe any directivity in reconstructed sources position shift. However, we performed simulations for small number of sources, therefore we cannot de nitely conclude that CLEAN-Vis or RL does not have a tendency for systematical shift of sources in x or y directions. The middle panel of Figure 10 presents semi-axes of tted ellipses. The result for CLEAN-Vis is what we expected from RHESSI experience. Namely, sizes are systematically larger than simulated ones. Both semi-axes were larger up to 2 times in comparison with simulated sources. The RL algorithm reconstructs sources sizes comparable to simulated ones. A similar situation holds for source brightness. The CLEAN-Vis underestimated the brightness which can be seen Figure 10 (middle column) as signi cantly higher background. The photometry from images reconstructed with RL algorithm is very good. At this point we are not able to compare with PIXON algorithm as it is not implemented for STIX data yet. However, these results shows the strength of our very simple algorithm: very good intensity and size reconstruction. This is important from physical interpretation point of view. Namely, a are energy budget analysis strongly depends on the density derived for the source (Ne = √︀ EM/V). When emission measure (EM) is underestimated, and volume (V) is overestimated the density is very biased. Thus, a reliable size and emission measure obtained from image are essential.
The parameters obtained with our algorithm were estimated from the images reconstructed after few hundreds of iteration steps. We do not use any regularization. The iteration was stopped when the change in the calculated χ 2 was less than 10 −4 . This value was estimated empirically. We observed that when the best similarity between the simulated and reconstructed images is achieved then the χ 2 changes very slightly. However, after next 10000 iteration steps the image quality might decrease, and χ 2 becomes worse. We estimated empirically that the quality of the image is good enough (less than 1% di erence with regard to the best image) when the χ 2 changes to less than 10 −4 . Usually the criterion was achieved after 300-500 iteration steps. However, it should be mentioned that the images which were closest to the CLEAN-Vis ones were obtained after 20-50 iteration steps only.

Square and rectangular approach differences
Our reconstruction algorithm, Maximum Likelihood, and Expectation Maximization (Lay and Katsaggelos 1990) are based on the same Richardson-Lucy approach. Therefore, we can investigate the degree of simpli cation introduced by assumption of the square wave functions. For this purpose we reconstructed images with our algorithm but separately for square and rectangular transmission functions. The result for a single source is presented in Figure 11. Top row of top panel presents images obtained with real STIX geometry after 10,1000, and 2000 reconstruction steps. The same is shown in bottom row of top panel, but assuming square transmission functions. There is signi cant di erence comparing arti cial sources which are much more signi cant when using simpli ed square transmissions. About 10% of signal is spread outside source when square transmission is used. Moreover, the source size is larger for square transmission. It may be connected to the fact that the convergence stops at some level which is illustrated in the χ 2 behaviour (Figure 11, top row of bottom panel). In addition, the reconstructed source ux stops at a xed level when square transmission is used (Figure 11, bottom row of bottom panel). We observed similar e ects for other sources con gurations. STIX, like any experiment in space, is a "living" instrument. Namely, the geometry of construction elements may change during entire mission. The imaging philosophy of the STIX is based on extremely precise development and alignment of grids. Grids were chosen to be rectangular to achieve the best possible angular resolution, and they were aligned precisely. The image reconstruction algorithm has to have a capability for utilization of data which comes from real instrument, not idealized one. For that reason we checked how the counts measured in stripes A,B,C, and D changes when alignment of grids is not perfect.
Again, we performed number of tests, but here we present only one example for clarity. Figure 12 presents changes of count number in the stripe C in a case of misalignment of front grid of collimator #18 (see Figure 4 in Krucker et al. 2020). We rotated front grid by an angle changing from −0. ∘ 5 to +0. ∘ 5. Change in registered counts is clearly visible. Moreover we saw large di erence in counts measured in square and rectangular approximations, which comes mainly from di erent transmissions. However, there is another source of di erence which is seen when square and rectangular counts are divided (red dashed line in Figure 12). The red line reveals that square approximation may be a source of additional errors which do not scale linearly with grid misalignment degree.
The misalignment analysis is easy and straightforward with our algorithm. It allows to calculate transmissions taking into account possible shifts, rotations, and skewness of individual grids. Moreover, we can include thickness of grids which is important for energy-dependent transmission, and/or for photons that hit grids at extremely large angles.  Changes of signal measured in the stripe C taking into account misalignment (rotation) of the first grid of collimator#18 (finest one). Black line presents rectangular grids, blue line is for square grids, and red line presents signals ratio.

Conclusions
In this paper we presented the simple image reconstruction method based on the Richardson-Lucy algorithm. Similar method exists (Massa et al 2019), but their methodology is di erent. We calculate photon source transmission (images) to counts directly (single-pixel response maps) while transmission functions used by Massa et al (2019) are approximated via visibilities. This additional step, from count rates to visibilities and back, can be a source of additional image degradation because it is fragile to noise in measured counts. Moreover, we have to keep in mind that the work of Massa et al (2019) is based on a square wave train approximation (Giordano et al. 2015), which is not a real STIX grids geometry.
With pixel count data calculated for point sources we were able to reconstruct images that shows reliable sources. Their positions, sizes and intensities are very similar to simulated sources. At this point we were able to compare our method with CLEAN-Vis algorithm which is the only other one implemented for STIX data. The main conclusions about presented algorithm performance are as follow: -The algorithm gave stable results with small tendency for over-resolution which may be easily avoided by using the stopping parameter related to the rate of χ 2 change. -The algorithm is moderately fast, but the result comparable to CLEAN algorithm is obtained within 20-50 iteration steps which takes less than 2 seconds on typical portable computer con guration. -Sources are reconstructed very well. The location, size and intensity are very close to simulated ones. Therefore the algorithm is very well suited for the detailed photometry of the solar HXR sources. -The algorithm used with real STIX geometry gives better results in comparison to simpli ed, square transmission analytical approach used elsewhere. -Simple approach we used, allows to improve photon transmission calculation in case of grids thickness, and uncertainties.
The algorithm presented here is the most straightforward approach to the image reconstruction problem. For the rst time we use real set of parameters describing STIX subcollimators. Despite the algorithm simplicity, it is very well constrained for detailed photometry of solar HXR sources.
Istvan Etesi for developing are simulation software and implementation of the CLEAN-Vis algorithm. This work was supported by the Polish National Science Centre grant number 2015/19/B/ST9/02826.

Conflict of Interest:
The authors have no con icts of interest to declare. All co-authors have seen and agree with the contents of the manuscript and there is no nancial interest to report.