Principle of Fredholm image reconstruction in the vignetting zone of an externally occulted solar coronagraph
Application to ASPIICS
^{1}
Université Côte d’Azur, Centre National de la Recherche Scientifique, Observatoire de la Côte d’Azur, UMR7293 Lagrange, Parc Valrose, 06108, Nice, France
email: claude.aime@univcotedazur.fr, celine.theys@univcotedazur.fr, henri.lanteri@univcotedazur.fr
^{2}
European Space Research and Technology Center, European Space Agency, Keplerlaan 1, 2201 Noordwijk, The Netherlands
Received:
13
July
2018
Accepted:
13
January
2019
Aims. This study is carried out in the context of data processing of the solar coronagraph ASPIICS of the future formationflying mission Proba3, which is expected to provide images of the corona very close to the limb. There will be a transition zone of the order of 100 arcsec close to the limb, where the telescope aperture suffers a strong vignetting by the external occulter (a disc of 1.42 m at 144 m). The instrument response in this region will vary rapidly both in shape and in integrated intensity, the latter being particular to the external occultation. The aim of this paper is to propose a technique to recover as much as possible of the image of the corona very close to the limb in the vignetting zone.
Methods. The object image relationship in this zone is not defined by the usual convolution but by the more general Fredholm integral of the first kind. Theoretical aspects of the problem are detailed in the context of a matrix formalism for the inversion of the Fredholm integral, formalism that we maintain up to the end of the numerical simulations, which is specific to the present work. The iterative RichardsonLucy algorithm, specially written for the nonconstant integrated intensity of the responses is used here for reconstruction. A study of the effect of noise on a photodetected image is made.
Results. An important part of the work consisted in calculating the elements of the transfer matrix between the object and the image for a simulation on a small region of size 100 × 100 arcsec sampled over 128 × 128 pixels. This is obtained propagating the light through the system using a previously published approach. On a toy object, the reconstruction is excellent down to about 60 arcsec from the limb, corresponding to a vignetting of 50%. The drawback is that the recovery of a N × N object requires the handling of a N^{2} × N^{2} matrix, i.e. a 16384 × 16384 transfer matrix here. However, taking into account radial symmetries of the experiment, we propose the use of a transformation from Cartesian to polar coordinates which allows to apply the same procedure all around the sun as for a small region.
Key words: techniques: high angular resolution / methods: numerical / techniques: image processing / Sun: corona
© ESO 2019
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The visual photometer imagined by Evans (1948) to measure the brightness of the sky near the sun is the precursor of externally occulted solar coronagraphs that are coupled with a classical Lyot (1939)’s coronagraph. The photometer contains some of the principal elements of these hybrid instruments: the external occulter; the internal occulter, the function of which is to intercept some of the brilliant edge of the external occulter; and a diaphragm that operates as a Lyot stop. Evans (1948) was well aware of vignetting effects inherent to externally occulted instruments and the difficulty of observing the sky brightness very close to the solar limb.
Efficient externally occulted coronagraphs were later developed for spaceborne coronagraphy, as described by Purcell & Koomen (1962) and Newkirk & Bohlin (1965). An excellent review of these techniques is made by Koutchmy (1988). The SOHO/LASCO mission (Brueckner et al. 1995) was strongly affected by effects of vignetting and corresponding loss of resolution. Considering geometrical optics, the efficient aperture is what is left inside the telescope aperture and outside the projected image of the occulter; see for example Bayanna et al. (2011).
A clear presentation of the problem can be found in Llebaria et al. (2006). They described the variation in the field of the point spread function (PSF) of stellar sources with the LASCOC2 coronagraph. The response depends on the direction of the point source on the sky, and, as they noticed, is further modified by the propagation inside the Lyot coronagraph itself. For LASCOC2, the effect of vignetting is extremely strong since it occurs for the whole field, due to the closeness of the occulter to the entrance aperture.
A large occulter is needed to perform an observation close to the solar disc, and that is the objective of the ASPIICS coronagraph, which will fly on the Formation Flying mission Proba3 (Lamy et al. 2010). The external occulter will be a disc of 1.42 m located at about 144.3 m from a second spacecraft bearing a Lyot coronagraph. The observation will be made in a narrow domain of wavelength around 550 nm. The entrance aperture of the Lyot coronagraph is 5 cm, giving a theoretical resolution of about 2.3 arcsec on the solar corona. These conditions (a large occulter and a small aperture) are favourable to reduce the effect of vignetting. However, ASPIICS is not exempt from these effects, the vignetting leading to a loss of observation very close to the solar limb.
The response begins to deteriorate at about 100 arcsec from the solar limb, and completely disappears at about 20 arcsec from the solar limb, as was already shown in Fig. 14 of Rougeot et al. (2017). The objective of the present study is to describe a method to recover as much of the image as possible in this severely disrupted region. To establish the coronagraph responses, we use the same computational technique as Rougeot et al. (2017). It is clear that the difficulty of reconstructing the image close to the limb increases with the vanishing of the PSF towards the limb.
The problem to solve is known as the inversion of the Fredholm integral equation of the first kind. It is a delicate problem that has given rise to a significant amount of literature (see e.g. Landweber 1951; Hanson 1971; Strand & Westwater 1968 and others). As described by these authors, the solution to the Fredholm integral, that is the image reconstruction, strongly depends on the illposedness of the problem and, classically, the problem of image reconstruction is solved by minimizing an adequate cost function composed of a data fidelity and a regularization term (Titterington 1985; Demoment 1989). In this preliminary work, the regularization is not considered and the cost function is only composed of the data fidelity term. Due to the illposedness of the problem, it is essential to use an iterative algorithm and to stop the iterations before the occurence of noisy unstable solutions. This procedure is called early stopping and we propose to use the wellknown Richardson Lucy (RL) algorithm, (Richardson 1972; Lucy 1974, 1994), slightly modified for a variant PSF (Lantéri et al. 2006).
The varying PSF, once the problem is discretized, gives rise to a high computational load. A solution proposed in the context of adaptiveoptics imaging, is to partition the image into regions where the PSF is considered to be invariant, to treat them independently, and to synthesize the image at the end (La Camera et al. 2015; Denis et al. 2015; Thiébaut et al. 2016; Flicker & Rigaut 2004). In the present study, the variation of the PSF, both in shape and integrated intensity, is too fast within a small zone to allow for the use of any approximation as we discuss below. A fundamental advantage over adaptive optics work is that here the variation of the response is deterministic. Pointspread functions can be calculated numerically for any point in the field, and it is not necessary to use interpolation. The drawback is that one has to handle very large matrices and it is not possible to use Fourier transforms. Vignetting is quite specific to external occulter coronagraphs; it does not exist in the case of a classical Lyot coronagraph where the response is either fully suppressed or merely equal to the telescope Airy pattern.
The paper is organised as follows: in Sect. 2, we introduce the physical considerations of image formation for an externally occulted coronagraph, and the problem of aperture vignetting. Examples of spacevariant PSFs are given for ASPIICS. Section 3 is devoted to theoretical aspects. The Fredholm relationship is quickly established, the discretized problem and its formulation in a matrix formalism are given. The RL algorithm is expressed in the case of a variant point source response. It is shown that the recovery of an N × N image requires an N^{2} × N^{2} matrix to be handled, and the illustration is made for N = 128 in Sect. 4 on a toy object made of a series of narrow Gaussian spots. A study of noise effect is given for Poisson statistics. In Sect. 5 it is shown how the symmetries of the occultercoronagraph responses make it possible to simplify the problem using a Cartesian to polar transformation. This will be key to processing large, future solar coronal images. Conclusions are given in Sect. 6.
2. Physical considerations: the spacevariant response of ASPIICS
As indicated in the introduction, in the ASPIICS experiment, the vignetting effect is limited to a small angular region of about 100 arcsec around the solar limb. Nevertheless, it is worth the effort to gain a few tens of arcsecs towards the limb, which was one of the main goals of ASPIICS. Since the external occulter is far away (at z = 144.348 m), geometrical considerations are not enough, and the Fresnel diffraction becomes very strong in the vignetting region and must be taken into account, as described in Rougeot et al. (2017).
An illustration of the vignetting effect produced by the external occulter for the observation of the solar corona is presented in Fig. 1. This is a concise illustration that requires some description. Instead of keeping the telescope fixed and scrolling the wave for several directions, we kept the diffraction pattern fixed and move the telescope, which results in the same figures for waves inside the aperture.
Fig. 1. Illustration of the vignetting effect produced by the external occulter. The diffraction pattern of the circular external occulter for an incident plane wave (region: xaxis from 600 mm to 800 mm, yaxis by ±30 mm, origin at the occulter centre) is shown in the background. The geometric edge of the occulter is marked “Oc”, (radius 710 mm), while the stenope image of the Sun is marked “S” (radius of 671 mm). The three circles delimitate intensity beams going through the 50 mm Lyot coronagraph entrance aperture. These correspond to point sources at 914 arcsec (A, occulted), 1000 arcsec (B, vignetting) and 1085 arcsec (C, full illumination) from the solar centre. 

Open with DEXTER 
For regions such as (A), the light coming from the solar disc is fully obscured by the external occulter. For regions such as (C), the aperture is fully illuminated, although a residual effect of diffraction fringes remains there. The transition regions, such as (B), correspond to a strong vignetting. The main variation occurs over ∼70 arcsec, which corresponds to the angle covered by the pupil of the telescope at the distance of 144.3 m.
The wavefront crossing the telescope aperture therefore propagates through ASPIICS’s Lyot coronagraph. Details of the computation of this propagation can be found in Rougeot et al. (2017). The numerical implementation of the propagation is done for the same conditions (number of points 8192 × 8192). The size of the internal occulter used in the present study has an overoccultation of 1.0028 with respect to the external occulter (this corresponds to a disc of 1424 mm onto the occulter) and a Lyot stop reducing the aperture to 48 mm (instead of 50 mm). These diaphragms change the shape of the PSF, but do not introduce a vignetting effect. New numerical computations will have to be performed once the exact parameters of the true experiment are frozen (currently, Proba3 is expected to fly in 2021). We note that possible defaults or misalignments are not taken into account in the present study (see Shestov & Zhukov 2018). Figure 2 gives responses for point sources at three distances from the solar limb (the centre of the Sun being at 960 arcsec from it), over a region limited to 30 × 30 arcsec. Intensities in the Lyot stop plane (referred to as plane C in Rougeot et al. 2017) are also given. A careful examination of these figures, in comparison with the occulter diffraction pattern of Fig. 1, shows how the wave has already been altered in the propagation through the Lyot coronagraph. The leftright inversion is real and due to the propagation inside the coronagraph.
Fig. 2. Left panels: illumination of the image of the aperture in the Lyot stop plane of the coronagraph, with corresponding PSFs (right panels) for unit point sources originating from the corona at 48, 72, and 96 arcsec from the solar limb, or 1008, 1032, and 1056 arcsec from the solar centre. The image in this plane is inverted with respect to the pupil of Fig. 1 and the wave begins to be altered by the propagation through the coronagraph. The PSF degrades faster in the radial direction than in the transverse direction (see Fig. 3 for numerical values). 

Open with DEXTER 
Pointspread functions evolve both in shape and integrated intensity, as shown in Fig. 3. An example is given for five point sources plotted together (top of Fig. 3). The intensity transmission curve is obtained computing 128 PSFs equally spaced from 19.5 to 114 arcsec. The integrated PSF, integration computed over a region of 96 × 96 arcsec around the photocentre of each PSF, is normalized to 1 for the shiftinvariant PSF. Below 20 arcsec a point source is almost completely obscured, while the transmission joins the nominal level beyond approximately 95 arcsec. The Strehl ratio, which quantifies the quality of the image, is also shown in this figure.
Fig. 3. Top panel: responses to unit point sources situated at 33.6, 52.3, 70.9, 89.5, and 108 arcsec from the solar limb for ASPIICS. A nonlinear greyscale (intensity to the power 0.2) is used to make the responses more visible. Bottom panel: blue curves, left y axis, resolutions in terms of 1D (radial and transverse) widths, and 2D integrated widths, and red curves, right y axis, integrated PSFs normalized to 1 for the invariant response. The black curve gives the Strehl ratio. Points (black squares) are at positions of the PSFs shown in the top panel. 

Open with DEXTER 
Equivalent widths of the PSFs are also drawn (in blue) in this figure. To account for the anisotropy of the response h(x, y) in the vignetting zone, we compute the 1D equivalent widths Φ_{1Dr} in the radial direction (here along x) and in the transverse direction Φ_{1Dt} (here along y) as a function of the distance from the limb. These quantities are determined as follows, the response h(x, y) is integrated in the y or x directions (here also for a 96 × 96 arcsec region), respectively, to get the 1D functions:
For the 1D functions, equivalent widths are the widths of the rectangles, which have the same surface (width × height) as the integrals of the functions F_{r}, F_{t}, and for heights the maximal values of these functions. For the 2D function, we compute the volume Σ under the function and take the width Φ_{2D} to be the diameter of the cylinder with the same volume Σ. The height of the cylinder is here also the maximal value of the function. The following relations connect these quantities:
Figure 3 clearly shows that the degradation is worse in the radial Φ_{1Dr} direction than in the transverse Φ_{1Dt} direction, as expected from the vignetting of the aperture by the external occulter. The circular averaged equivalent width stands in between. Beyond 90 arcsec from the limb, the response becomes isotropic, and Φ_{1Dr} and Φ_{1Dt} merge to the same value, which is not Φ_{2D}, simply because of the differences between the definitions of the 1D and 2D widths.
A 3D illustration of the Fredholm problem is given in Fig. 4, in opposition to the convolution case. This is another representation of the top panel of Fig. 3, illustrating the rapid variation of the PSF. It would be impossible to partition the image into very small regions to model smooth changes, as evoked by Denis et al. (2015) in their consideration of vignetting in adaptive optics. In its general form, given in Eq. (3) of Sect. 3, the response h(α, ξ, β, η) is a function of four variables. However, in the particular case of the externally occulted Lyot coronagraph, the variation of the response only depends on the radial distance of the point source from the solar centre. This property, that is further used in Sect. 4, makes it possible to represent h(α, ξ, β, η) as a function of three variables, ρ, ξ and η.
Fig. 4. Threedimensional illustration showing PSFs in the two cases of invariant (top panel) and space variant (bottom panel) cases. Vertical axes ( ρ) correspond to positions of points sources in the object at four radial distances from the limb: a = 26.6, b = 48.8, c = 71, and d = 93.3 (in arcsec). The horizontal rectangles show PSF responses in the image plane (a region of η = 43.7 × ξ = 94 arcsec). A nonlinear greyscale is used to make the curves more visible. 

Open with DEXTER 
Figure 5 shows in the Fourier space the rapid variation of the responses. It represents the modulus of the Fourier transform of the space variant PSFs (called pseudoMTF) of Fig. 4 for the positions b, c, and d corresponding to point sources in the solar corona at 48.8, 71, and 93.3 arcsec. All pseudoMTFs are normalized to 1 at the zero frequencies of the u × v plane to better highlight the variation of the spectral coverage of these responses. The transfer function for the PSF at the position d is almost that of the unvignetted aperture, while it practically disappeared for the position a, which was not figured here. In addition to the curves in Fig. 3, this representation is another illustration in the Fourier plane of how fast the variation of the image PSF occurs.
Fig. 5. Modulus of the Fourier transform of the space variant PSFs (called pseudoMTF) of Fig. 4 for the planes b, c, and d corresponding to point sources in the solar corona at 48.8, 71, and 93.3 arcsec. The spatial frequency u is in the limb direction. 

Open with DEXTER 
3. Theoretical aspects: the matricial formalism and the inversion of the Fredholm integral
With such a variant point response, the relation that connects the object X(α, β), the observation Y(α, β), and the point responses analysed above is no longer the usual convolution relationship, but a more complex integral that is known as the Fredholm integral of the first kind, (Landweber 1951). In the following, we use the notations of Theys & Aime (2016) to derive that relation. For that, we have to denote with different symbols the angular coordinates in the object and image planes. We use (ξ, η) to denote the angular direction in the object plane, and (α, β) in the image plane. Coordinates (ξ, η) play the role of dummy integration variables. An elementary angular region dξdη around the object X(ξ, η) produces an elementary intensity X(ξ, η)dξdη. This intensity is distributed in the image plane according to the instrument response h(α, ξ, β, η), a function of four variables (see Figs. 2 and 3 for illustrations). The object being incoherent, intensities coming from all directions are added together, and we obtain
which is a 2D Fredholm integral of the first kind. Recovering X(α, β) from Y(α, β) is the problem we seek to solve. As already indicated, this problem is more complex to solve than the usual deconvolution problem where the point spread response is shift invariant:
Equation (3) does not allow the use of Fourier transforms as Eq. (4) does.
In a general manner, the discretized model of the Fredholm relation, Eq. (3), leading to the observed process can be written as
where y(m, n), x(i, j) and h(m, i, n, j) are the values of the corresponding arrays.
The values of the image given by Eq. (5) can be put in the matrix form:
where y and x are vectors containing the discretized values of the image and the object. Here, H is the matrix containing elements h(m, i, n, j). A detailed description of the elements of matrix H for both cases of shiftinvariant and shiftvariant PSFs are given for a numerical example in Sect. 4 relative to numerical illustrations.
The estimation of x, given y and H from the matrix form Eq. (6), is a socalled illposed problem and its direct inversion gives poor results since H is very illconditioned. The classical solution is then to use iterative methods to estimate x, which amounts to performing a partial inversion of H. The algorithm used here is a slightly modified version for a variant PSF, Lantéri et al. (2006), of the socalled RL algorithm :
This writing clearly highlights the term related to the flux variation of the PSF, [H^{T}1_{N}]_{r}, that is the sum of the values of the PSF, which varies with r; see the red curve in Fig. 3. In the case of a constant and normalized PSF, the term and Eq. (7) simplifies in the RL algorithm; in the context of illposed problems, instability in the solution appears as the number of iterations increases. The problem is then to determine the optimal iteration number and to stop the iterative process there, for obtaining a physically satisfactory solution, (Lucy 1994). A classical alternative is to introduce an explicit regularization which is left for future work.
4. Numerical illustrations
4.1. Matrix implementation
According to Eqs. (5) and (6), a given value of the image y(m, n) depends on the N × N weighting array :
where y and x are vectors containing the discretized values of the image and the object is lexicographically stored. The term “lexicographic” infers that the values of each table of dimension N × N read from left to right and from top to bottom are arranged in this order in a vector. An image of size N × N then yields a vector y of length N^{2}.
Consequently the matrix H is of dimensions N^{2} × N^{2}.
The array values of each weighting function h(m, i, n, j) are ordered lexicographically on a line of H to form with x the data y(m, n):
We note that each column of H corresponds to one of the PSFs shown in Fig. 2, ordered lexicographically. In the spaceinvariant case, the matrix H has a ToeplitzblockToeplitz structure and some lines can be added to form a circulant matrix. Subsequently, matricial products with H can be advantageously performed in the Fourier space using 2D discrete Fourier Transforms (DFT) on the N × N arrays. This however is not the case here and we have to use very large matrices.
4.2. Illustration on a test object in the vignetting zone
We considered the processing of objects and images of size 96 × 96 arcsec, sampled over N × N = 128 × 128 points. Each PSF is taken as being the same size as the object, that is 96 × 96 arcsec, around the photocentre response. This is not mandatory, and a smaller dimension could have been chosen. These dimensions were a tradeoff between adequate representation of the vignetting zone and an N value of reasonable size, that is, not too large to be handleable with our computers. The size of the matrix H is therefore 16384 × 16384.
The first and most difficult task is to compute 16384 PSFs. This was done using the same computer program and machine (with two 14core Intel Xeon processors and 512 GB of RAM) as in Rougeot et al. (2017). As already indicated in Sect. 2, the computation involves the Fresnel diffraction from the external occulter to the entrance telescope aperture and the propagation through the Lyot coronagraph. Obtaining all 16384 PSFs would have taken days of computing time, which is not impossible but finally not necessarily essential, as we see below. A reasonable approximation was obtained considering that for the small region of 96 × 96 arcsec, the variation of the response is essentially in the radial direction. Indeed, in the transverse direction, responses are affected by a small rotation that is less than 3 degrees at the most.
The procedure we used was the following. We exactly computed 128 radial responses. Each of them was given in an array of 8192 × 8192 points of which we extracted regions of 256 × 256 points. The responses in the transverses directions are obtained by translating subarrays of 128 × 128 points inside the 256 × 256 array. This is also what was previously done to obtain the PSF parameters drawn in Fig. 3. Using a larger array was necessary not to be affected by edge effects. Each of the 128 × 128 point responses was then lexicographically scanned to become a column of the matrices H shown in Fig. 6, for the shiftinvariant response, denoted there H(convolution) and the shiftvariant response, denoted there H(Fredholm).
Fig. 6. Transform matrices H. A central part (600 × 600 points out of 16384 × 16384 points) is represented for (top panel) the shift invariant PSF of the convolution case, denoted H(Convolution), and for (bottom panel) the Fredholm case, denoted H(Fredholm), where the PSF varies from one point to another. 

Open with DEXTER 
As it is not possible to represent the 268 million points of each matrix; in Fig. 6, only a central part of 600 × 600 points is shown in a representation in levels of false colours. We note that the symmetric behaviour of the convolution matrix (except for sparse edge effects) becomes disrupted in the case of the Fredholm matrix. The variable level of intensity already shown in Fig. 3 is also clearly visible here. Figure 7 represents examples for two given rows and columns (here 6655 and 12929) that clearly show the symmetry of the convolution case (to a given shift) and the absence of symmetry of the Fredholm case. In the present computations, the condition numbers for the matrices are 10^{20} for the convolution case and 6.7 × 10^{22} for the Fredholm case.
Fig. 7. An example of rows and columns of the transform matrix H shown in Fig. 6 showing the different behaviour between Fredholm and convolution cases. The symmetry of the rows and columns of the convolution case is lost in the Fredholm case, as expected. Due to the lexicographic order for the image of 128 × 128 points, there is a pseudo periodicity of 128 elements of the matrices H (except for edge effects). Each row and column is made of 16384 points, which makes it difficult to draw. 

Open with DEXTER 
Computations are made on a toy object with the parameters of ASPIICS. Figure 8 shows, from top to bottom, the true object x_{t} and the result of the product Hx_{t}, giving the image y and a noisy image resulting from a Poisson transform of Hx_{t} with a mean of 163840 photons in the whole image, that is, a mean number of ten photons by pixel. Figure 9 shows the obtained reconstructions at the iteration k̂ corresponding to the minimum of the reconstruction error ϵ_{k} with the RL algorithm, Eq. (7), for the noiseless and the noisy cases. The reconstruction error ϵ_{k} is defined as
where x_{t} is the true object, which is not known in practice. In the noiseless case, it is interesting to observe that there is a substantial improvement to the recovery of information very close to the solar limb, here some 30 or 40 arcsec towards the limb for this toy object. The RL algorithm produces a minimum reconstruction error of 0.40 for an early stopping of k̂ = 2000 iterations. In the noisy case considered, which corresponds to a very high level of noise for a solar observation, the RL algorithm produces a minimum reconstruction error of 0.4888 for an early stopping of k̂ = 27 iterations. It is a wellknown result, the noisier the image, the higher the reconstruction error for a smaller iteration number. It can however be noted that the reconstruction remains very good until 45–50 arcsec towards the limb. For a more accurate view of the reconstructions, a cut of the arrays is given in Fig. 10. The cut is made in the middle, where the dots are circular and equally spaced. The reconstructed object is for an iteration number k̂ = 2000 for the noiseless case and k̂ = 27 for the noisy case.
Fig. 8. From top to bottom panels: toy object made of an ensemble of Gaussian point sources, observed noiseless image, and noisy image with 163840 photons in the whole image, i.e a mean number of 10 photons by pixel. 

Open with DEXTER 
Fig. 9. Reconstructed object using RL at optimal iteration number k̂ = 2000 for the noiseless image (top panel) and at k̂ = 27 for the noisy image (bottom panel). 

Open with DEXTER 
Fig. 10. Cuts of the noiseless observed image y, the object x_{t} and the reconstructed object using the RL algorithm. Top panel: noiseless image. Bottom panel: noisy image. See comments in the body of the paper. 

Open with DEXTER 
In the noiseless case, the reconstruction is excellent down to about 60 arcsec from the limb. There, the integrated intensity has already fallen to 54% of its initial value, the Strehl ratio to 0.3, and the longitudinal PSF width has increased by 83%, as can be seen in Fig. 3. Between 60 and 45 arcsec from the limb, the algorithms start to experience difficulties. In that zone, reconstructed images remain similar to the original object, but the photometry is not preserved, the reconstructed level fluctuates depending on the iteration number. In the example of Fig. 8, one of the peaks is higher, and the other lower than expected. At 45 arcsec from the limb, the transmission has fallen to 27%, the Strehl ratio is less than 0.08, and the longitudinal width has been multiplied by three. Beyond that limit, the reconstruction becomes ineffective and the signal almost disappears below 25 arcsec.
In the noisy case, the behaviour of the algorithm is the same with lower values, the reconstruction is excellent down to about 60 arcsec from the limb with more fluctuating values of the reconstructed intensity. From 60 arcsec the reconstructed intensity falls gradually and the longitudinal width of the peaks increases. The signal disappears below 35 arcsec.
A brief study of the influence of noise on the algorithm is given in Fig. 11. This figure gives the normalized reconstruction error (black curve) and computation time (red curve) for the reconstruction as a function of the mean number of photons by pixel, in logarithmic scale. Results on noisy data, shown in the bottom panel of Fig. 9, correspond to the minimum number of photons of Fig. 11, that is, ten photons per pixel. In this case the reconstruction error is 0.488 and the computation time is less than 10 s. Obviously, the reconstruction error decreases as the number of photons increases and in parallel the computation time increases since the corresponding number of iterations for the early stopping increases. We note that these curves are obtained by stopping the algorithm as soon as the reconstruction error increases. The first three values correspond to low numbers of photons and the results are irregular. It can be noted that the calculation time is at most equal to 10 min, which remains very small compared to that of PSFs, which is spread over several hours but calculated once.
Fig. 11. Black curve, left y axis, reconstruction error as a function of the mean number of photons. Red curve, right y axis, computation time in seconds as a function of the mean number of photons, for a MacBook Pro 2018, i9 processor. 

Open with DEXTER 
5. Reconstructing the entire vignetting zone using a Cartesian to polar transform
The reconstruction of the vignetted zone described above has been done for a small area of N × N points, with N = 128. Since the size of the matrix H is N^{2} × N^{2} points, here with N^{2} = 16384, it is almost impossible to envisage the inverse Fredholm procedure for an image with N much greater than that. Moreover, considering a large region, no approximation can be made for the PSF that must then be computed for all point sources, making the calculation even more difficult.
As we have already discussed it, the method proposed for space varying PSF in the context of adaptive optics imaging, (La Camera et al. 2015; Thiébaut et al. 2016) cannot be used here.
In the particular case of coronagraphy with an external occulter, a very interesting alternative can be imagined considering that a PSF depends solely on the distance to the centre of the occulter and rotates identical to itself around the solar limb. To take advantage of this property we apply to the vignetted region a transformation that we denote as T, giving the value of the observed image from Cartesian coordinates (α, β) into polar coordinates (ρ, θ). This transformation is illustrated in Fig. 12. It changes the ring of the vignetted region into a rectangle. The figure is not to scale and only half of the solar image and corona is represented here so as not to take up space. The inner dashed halfcircle of radius OM corresponds to the solar surface masked by the occulters (external and/or internal). The outer dashed halfcircle of radius OM’ delimits the vignetted region. A simple drawing (in red) is made to schematize the lower solar corona there. This drawing can be followed in the rectangle after the transformation T.
Fig. 12. Schematic drawing of the transformation T from the (α, β) cartesian space to the (ρ, θ) polar space. In this transform, the ring corresponding to the vignetting region becomes a rectangle (here only half the corona region is shown). The drawing is not to scale, since OM corresponds to the Sun radius (about 960 arcsec) and the width (MM’) is the vignetted region (about 100 arcsec). 

Open with DEXTER 
A block diagram of the operations to apply is given in Fig. 13. We note that after the T transformation, the responses themselves are transformed. It would not be interesting to use this transformation outside the vignetting domain, because then the response, naturally invariant there, would become variable in ρ.
Fig. 13. Block diagram of the proposed procedure to reconstruct the vignetted zone using one and a single H matrix transform. T is the (α, β) to ( ρ, θ) transform, and T^{−1} the ( ρ, θ) to (α, β) transform. 

Open with DEXTER 
An illustration of the rotation of the responses is given in Fig. 14. The direct image Y(α, β) (Fig. 14, top) corresponds to points sources positioned along radial directions from the solar centre. The solar disc is on the left, represented there in greyscale. The curvature is that of a circle of about 960 arcsec. The representation covers a region of 180 × 210 arcsec, sufficient to bring out the rotation of the features of the PSF; it makes use of the parameters of ASPIICS, assuming that all optical parts and settings are perfect.
Fig. 14. Illustration of the transformation T (Cartesian to polar) from the image Y(α, β) (top panel) to its transform Y′ ( ρ, θ). See details in the main text. 

Open with DEXTER 
The transformed image Y′ (ρ, θ) (Fig. 14, bottom) is now a function of ρ, in arcsec, and θ in degrees (or in radians). In the figure the rectangle is represented for θ varying from 0° to 11°. In this representation, the circle corresponding to the solar limb is now a vertical line. The way the transformation has been implemented is by interpolating from a Cartesian to a polar grid.
For that we use a new array of K × L equispaced points in ρ and θ. A point P_{k, l} is defined there by its radial coordinates ρ_{k} = ρ_{0} + kΔρ and θ_{l} = θ_{0} + lΔθ, where the values of ρ_{0}, Δρ, θ_{0}, Δθ, K and L are chosen to fit the region of interest for the transformation. Cartesian coordinates of P_{k, l} are ( ρ_{k} cos(θ_{l}),ρ_{k} sin(θ_{l})) which appear as an interior point of the grid square inside (m, n) and (m + 1, n + 1) of the known array y(m, n). Defining u = ρ_{k} cos(θ_{l})−m and v = ρ_{k} sin(θ_{l})−n the distances to m and n, with 0 ≤ u ≤ 1 and 0 ≤ v ≤ 1, a bilinear interpolation gives the value of y′ (ρ_{k}, θ_{l}):
which has the advantage of preserving positivity. As described by Press et al. (19881992), the bilinear interpolation is the simplest one in two dimensions, is sufficient for an array with smooth variations of values and is currently used in image transformations such as the ones described in Russ (1995). It can be seen that this linear combination will not modify the nature of a Gaussian additive noise. If we now consider a Poisson transformation of the image, as in the numerical examples proposed, since the coefficients are not integers, the noise in the polar coordinate system will change in nature. Another point to note is that the transformation T will introduce a correlation between the pixels and the RL algorithm will no longer be theoretically applicable in this case. However, we can anticipate that solar observations will be high flux and the effect of the noise will be negligible.
We note that if the ρ axis is in arcsec, the θ axis is in degrees. Due to the choice of sampling points used for the transformation T (here we oversampled the image with 1000 × 1000 points), the figure still resembles a direct space.
The main result is that now the responses have become independent of θ. It therefore becomes sufficient to calculate the responses along a radial line. This is exactly what we have done to establish the transform matrix H used in Sect 4. What we have previously considered as an approximation allowed by the treatment of a small area is now the exact computation to perform, thanks to the Cartesian to polar T transform.
Let us estimate in round figures the gain obtained by this operation. Suppose that for the resolution of 2.3 arcsec of the telescope, we take a PSF every 1 arcsec squared. In direct processing, without using the T transformation, it would be necessary to calculate every PSF included in a ring of about 100 arcsec around the limb, which is roughly 2π × R_{Sun} × 100, or about 6 million points sources ( since R_{Sun} = 960 arcsec), instead of just 100 using the T transform. The gain in computing time, simply to get the PSFs is more than a factor of 6000, and also requires less computer memory space. We note that in Sect. 4 we adopted a slightly tighter sampling (128 points for 96 arcsec, or 0.75 arcsec/point), and numbers obtained for the gain in computer time were a lower limit.
Let us denote X′ ( ρ, θ) and Y′ ( ρ, θ) the object and image after application of the transformation T, respectively. The point response h(ξ, η, α, β) then becomes a new function h′ (ξ′,η′,ρ, θ), which is now invariant in θ. The Fredholm equation analogous to Eq. (3) can now be written as
The matrix H, here H′, can be limited to a reasonable size to reconstruct X′ ( ρ, θ) in several pieces that can be stacked on top of each other. An inverse transformation T^{−1} will have finally to be applied to recover the image, as described in the block diagram of Fig. 13.
6. Conclusion
We have proposed a procedure to recover strongly blurred images in the vignetting zone of the externally occulted solar coronagraph ASPIICS. Close to the limb, the PSF varies rapidly in shape and in integrated intensity, which is specific to the external occultation. The object image relationship in this zone is then defined by a Fredholm integral of the first kind. We proposed to invert this problem using a complete matrix formalism and the used method is the nonsimplified form of the RL algorithm for a nonconstant PSF. This PSF was modelled using the approach of Rougeot et al. (2017).
Simulations were made for a simple toy object made of Gaussian point sources. The object is of size 128 × 128 points. Results are very positive since inversion of Fredholm integral gives satisfactory results up to a vignetting of 50% of the aperture. The drawback of the matrix approach is that the recovery of an N × N object requires the handling of an N^{2} × N^{2} matrix, which in the present study is a 16384 × 16384 matrix. Going beyond this dimension becomes difficult.
To overcome this difficulty, we proposed the use of a transformation from Cartesian to polar coordinates that makes it possible to handle the vignetting region of the whole solar image. This was obtained taking into account radial symmetries of the experiment. This transformation changes the ring of the vignetting zone with a variable PSF into a rectangle with a PSF that varies with ρ, the distance to the limb, but which is invariant in θ around the limb. This greatly facilitates both PSF calculations and Fredholm’s inversion.
Future work remains to be done on more realistic images of the solar corona. Such images should show real structures of the solar corona, for example obtained during a true solar eclipse, (Druckmüller et al. 2014), and should be contaminated by the diffraction halo of the solar chromosphere, as computed in Rougeot et al. (2017). This image should also take into account noise sources specific to real observations. Moreover, the transformation T remains to be implemented for an entire image of the solar corona. This will also require a special treatment of the continuity of the image between Fredholm and convolution zones.
A rapid analysis of noise effects is given in Sect. 5. To make these effects visible, a very small number of ten photons per pixel is used in the illustrations. This is certainly a very low number of photons because the area of vignetting corresponds to the brightest zone of the solar corona. The dominant effect in the vignetting zone remains the degradation of the shape of the PSF.
In addition, a reconstruction algorithm adapted to the real data model will have to be developed. A regularization term must be added to overcome early stopping that requires knowledge of the object to be reconstructed.
The technique presented here can be extended to a coronagraph using a serrated external occulter. In this case, the variation of the PSF is no longer only radial, but it also varies angularly with a periodicity of 2π/N_{t}, where N_{t} is the number of teeth. Such coronagraphs have been proposed since precursors projects of Newkirk & Bohlin (1965) and Purcell & Koomen (1962). A recent study by Rougeot & Aime (2018) has shown that these serrated occulters can indeed outperform that the sharpedged disc, provided that the number of teeth is high (e.g. N_{t} = 1000). In this case, the number of PSFs to be computed to establish the matrix H becomes larger, of the order of ten times larger for N_{t} = 1000, which remains quite feasible. Work on this point is being developed and will be the subject of a future contribution.
Finally, a study should be done to investigate whether or not the T technique could find applications in other areas of astronomy where the PSF also suffers from radial variation.
References
 Bayanna, A. R., Mathew, S. K., Sankarasubramanian, K., et al. 2011, Exp. Astron., 29, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, Sol. Phys., 162, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Cash, W. 2006, Nature, 442, 51 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 DaubeWitherspoon, M. E., & Muehllehner, G. 1986, IEEE Trans. Med. Imaging, 5, 61 [CrossRef] [PubMed] [Google Scholar]
 Demoment, G. 1989, IEEE Trans. Acoust. Speech Sig. Proc., 37, 2024 [CrossRef] [Google Scholar]
 Denis, L., Thiébaut, E., Soulez, F., Becker, J. M., & Mourya, R. 2015, Int. J. Comput. Vis., 115, 253 [CrossRef] [Google Scholar]
 DePierro, A. R. 1987, IEEE Trans. Med. Imaging, 2, 328 [Google Scholar]
 Druckmüller, M., Habbal, S. R., & Morgan, H. 2014, ApJ, 785, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Evans, J. W. 1948, J. Opt. Soc. Am., 38, 1083 [NASA ADS] [CrossRef] [Google Scholar]
 Flicker, R. C., & Rigaut, F. J. 2004, J. Opt. Soc. Am. A, 22, 3 [Google Scholar]
 Hanson, R. J. 1971, SIAM. J. Numer. Anal., 8, 616 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Karush, W. 1939, Ph. D. Thesis, Univ. of Chicago, USA [Google Scholar]
 Koutchmy, S. 1988, Space Sci. Rev., 47, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Kuhn, H. W., & Tucker, A. 1951, Proc. 2nd Berkeley Symp., (Berkeley: Univ. of California Press), 481 [Google Scholar]
 La Camera, A., Schreiber, L., Diolaiti, E., et al. 2015, A&A, 579, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lamy, P., Damé, L., Vivès, S., & Zhukov, A. 2010, SPIE, 7731, 18 [Google Scholar]
 Landweber, L. 1951, Am. J. Math., 73, 615 [CrossRef] [MathSciNet] [Google Scholar]
 Lantéri, H., Roche, M., & Aime, C. 2002, Inverse Prob., 18, 1397 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lantéri, H., Roche, M., Cuevas, O., & Aime, C. 2006, Sig. Proc., 81, 945 [CrossRef] [Google Scholar]
 Llebaria, A., Lamy, P., & Danjard, J.F. 2006, Icarus, 182, 281 [NASA ADS] [CrossRef] [Google Scholar]
 Lyot, B. 1939, MNRAS, 99, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1974, AJ, 79, 745 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1994, A&A, 289, 983 [NASA ADS] [Google Scholar]
 Newkirk, Jr., G., & Bohlin, J. D. 1965, IAU Symp., 23, 287 [NASA ADS] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 19881992, Numerical Recipes in C (Cambridge: Cambridge University Press) [Google Scholar]
 Purcell, J. D., & Koomen, M. J. 1962, Coronagraph with Improved ScatteredLight Properties, Report of NRL Progress (Washington, D.C: U.S. GPO) [Google Scholar]
 Richardson, W. H. 1972, J. Opt. Soc. Am., 62, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Rougeot, R., & Aime, C. 2018, A&A, 612, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rougeot, R., Flamary, R., Galano, D., & Aime, C. 2017, A&A, 599, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Russ, J. C. 1995, The Image Processing Handbook, 2nd edn. (Boca Raton: CRCPress) [Google Scholar]
 Shestov, S. V., & Zhukov, A. N. 2018, A&A, 612, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Strand, O. N., & Westwater, E. D. 1968, SIAM J. Numer. Anal., 287 [NASA ADS] [CrossRef] [Google Scholar]
 Theys, C., & Aime, C. 2016, Lect. Notes Phys., 914, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Thiébaut, E., Denis, L., Soulez, F., & Mourya, R. 2016, Proc. SPIE, 9909, 99097N [CrossRef] [Google Scholar]
 Titterington, D. M. 1985, A&A, 144, 381 [NASA ADS] [Google Scholar]
 Vanderbei, R. J., Cady, E., & Kasdin, N. J. 2007, ApJ, 665, 794 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1. Illustration of the vignetting effect produced by the external occulter. The diffraction pattern of the circular external occulter for an incident plane wave (region: xaxis from 600 mm to 800 mm, yaxis by ±30 mm, origin at the occulter centre) is shown in the background. The geometric edge of the occulter is marked “Oc”, (radius 710 mm), while the stenope image of the Sun is marked “S” (radius of 671 mm). The three circles delimitate intensity beams going through the 50 mm Lyot coronagraph entrance aperture. These correspond to point sources at 914 arcsec (A, occulted), 1000 arcsec (B, vignetting) and 1085 arcsec (C, full illumination) from the solar centre. 

Open with DEXTER  
In the text 
Fig. 2. Left panels: illumination of the image of the aperture in the Lyot stop plane of the coronagraph, with corresponding PSFs (right panels) for unit point sources originating from the corona at 48, 72, and 96 arcsec from the solar limb, or 1008, 1032, and 1056 arcsec from the solar centre. The image in this plane is inverted with respect to the pupil of Fig. 1 and the wave begins to be altered by the propagation through the coronagraph. The PSF degrades faster in the radial direction than in the transverse direction (see Fig. 3 for numerical values). 

Open with DEXTER  
In the text 
Fig. 3. Top panel: responses to unit point sources situated at 33.6, 52.3, 70.9, 89.5, and 108 arcsec from the solar limb for ASPIICS. A nonlinear greyscale (intensity to the power 0.2) is used to make the responses more visible. Bottom panel: blue curves, left y axis, resolutions in terms of 1D (radial and transverse) widths, and 2D integrated widths, and red curves, right y axis, integrated PSFs normalized to 1 for the invariant response. The black curve gives the Strehl ratio. Points (black squares) are at positions of the PSFs shown in the top panel. 

Open with DEXTER  
In the text 
Fig. 4. Threedimensional illustration showing PSFs in the two cases of invariant (top panel) and space variant (bottom panel) cases. Vertical axes ( ρ) correspond to positions of points sources in the object at four radial distances from the limb: a = 26.6, b = 48.8, c = 71, and d = 93.3 (in arcsec). The horizontal rectangles show PSF responses in the image plane (a region of η = 43.7 × ξ = 94 arcsec). A nonlinear greyscale is used to make the curves more visible. 

Open with DEXTER  
In the text 
Fig. 5. Modulus of the Fourier transform of the space variant PSFs (called pseudoMTF) of Fig. 4 for the planes b, c, and d corresponding to point sources in the solar corona at 48.8, 71, and 93.3 arcsec. The spatial frequency u is in the limb direction. 

Open with DEXTER  
In the text 
Fig. 6. Transform matrices H. A central part (600 × 600 points out of 16384 × 16384 points) is represented for (top panel) the shift invariant PSF of the convolution case, denoted H(Convolution), and for (bottom panel) the Fredholm case, denoted H(Fredholm), where the PSF varies from one point to another. 

Open with DEXTER  
In the text 
Fig. 7. An example of rows and columns of the transform matrix H shown in Fig. 6 showing the different behaviour between Fredholm and convolution cases. The symmetry of the rows and columns of the convolution case is lost in the Fredholm case, as expected. Due to the lexicographic order for the image of 128 × 128 points, there is a pseudo periodicity of 128 elements of the matrices H (except for edge effects). Each row and column is made of 16384 points, which makes it difficult to draw. 

Open with DEXTER  
In the text 
Fig. 8. From top to bottom panels: toy object made of an ensemble of Gaussian point sources, observed noiseless image, and noisy image with 163840 photons in the whole image, i.e a mean number of 10 photons by pixel. 

Open with DEXTER  
In the text 
Fig. 9. Reconstructed object using RL at optimal iteration number k̂ = 2000 for the noiseless image (top panel) and at k̂ = 27 for the noisy image (bottom panel). 

Open with DEXTER  
In the text 
Fig. 10. Cuts of the noiseless observed image y, the object x_{t} and the reconstructed object using the RL algorithm. Top panel: noiseless image. Bottom panel: noisy image. See comments in the body of the paper. 

Open with DEXTER  
In the text 
Fig. 11. Black curve, left y axis, reconstruction error as a function of the mean number of photons. Red curve, right y axis, computation time in seconds as a function of the mean number of photons, for a MacBook Pro 2018, i9 processor. 

Open with DEXTER  
In the text 
Fig. 12. Schematic drawing of the transformation T from the (α, β) cartesian space to the (ρ, θ) polar space. In this transform, the ring corresponding to the vignetting region becomes a rectangle (here only half the corona region is shown). The drawing is not to scale, since OM corresponds to the Sun radius (about 960 arcsec) and the width (MM’) is the vignetted region (about 100 arcsec). 

Open with DEXTER  
In the text 
Fig. 13. Block diagram of the proposed procedure to reconstruct the vignetted zone using one and a single H matrix transform. T is the (α, β) to ( ρ, θ) transform, and T^{−1} the ( ρ, θ) to (α, β) transform. 

Open with DEXTER  
In the text 
Fig. 14. Illustration of the transformation T (Cartesian to polar) from the image Y(α, β) (top panel) to its transform Y′ ( ρ, θ). See details in the main text. 

Open with DEXTER  
In the text 