Principle of Fredholm image reconstruction in the vignetting zone of an externally occulted solar coronagraph Application to ASPIICS

Aims. This study is carried out in the context of data processing of the solar coronagraph ASPIICS of the future formation-ﬂying mission Proba-3, which is expected to provide images of the corona very close to the limb. There will be a transition zone of the order of 100arcsec close to the limb, where the telescope aperture su ﬀ ers a strong vignetting by the external occulter (a disc of 1.42m at 144m). 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 deﬁned by the usual convolution but by the more general Fredholm integral of the ﬁrst 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 speciﬁc to the present work.The iterative Richardson-Lucy algorithm, specially written for the non-constant integrated intensity of the responses is used here for reconstruction. A study of the e ﬀ ect 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 × 100arcsec sampled over 128 × 128pixels. 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 60arcsec 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.


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 space-borne 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 out-side 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 LASCO-C2 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 LASCO-C2, 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 Proba-3 (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 ill-posedness 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 ill-posedness 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 well-known Richardson Lucy (RL) algorithm, (Richardson 1972;Lucy 1974Lucy , 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 adaptive-optics 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. Point-spread 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 space-variant PSFs are given for ASPI-ICS. 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 occulter-coronagraph 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.

Physical considerations: the space-variant 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.
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 over-occultation 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, Proba-3 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 left-right inversion is real and due to the propagation inside the coronagraph.
Point-spread 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 shift-invariant PSF. Below 20 arcsec a point source is almost completely obscured, while the transmission joins the nominal level beyond approximately A212, page 2 of 10 95 arcsec. The Strehl ratio, which quantifies the quality of the image, is also shown in this figure.
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: (2) 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 ρ = α 2 + β 2 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 η. 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 pseudo-MTF) 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 pseudo-MTFs 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 un-vignetted 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.

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  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).
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

PSF parameters
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 shift-invariant and shift-variant 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 so-called ill-posed problem and its direct inversion gives poor results since H is very ill-conditioned. 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 so-called RL algorithm : 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.

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 . y = y(1, 1)y(1, 2) . . . y(2, 1) . . . y(N, 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): 1, 1, 1) h (1, 1, 1, 2 (N, 1, N, 1) h (N, 1, N, 2) . . . h (N, 2, N, 1) We note that each column of H corresponds to one of the PSFs shown in Fig. 2, ordered lexicographically. In the space-invariant case, the matrix H has a Toeplitz-block-Toeplitz 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.

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 trade-off 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 14-core 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 sub-arrays 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).
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.
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 iterationk 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 ofk = 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 ofk = 27 iterations. It is a well-known 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 numberk = 2000 for the noiseless case andk = 27 for the noisy case.
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.

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 half-circle 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 A212, page 6 of 10  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.
there. This drawing can be followed in the rectangle after the transformation T.
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 ρ.
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.
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 ): y (ρ k , θ l ) = y(m + u, n + v) = (1 − v)(1 − u) y(m, n)+ uv y(m + 1, n + 1) + (1 − u)v y(m, n + 1) + (1 − v)u y(m + 1, n) which has the advantage of preserving positivity. As described by Press et al. (1988Press et al. ( -1992, 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 S un = 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 T T -1 X(a,b) H' 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.
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.

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 non-simplified form of the RL algorithm for a non-constant 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 sharp-edged 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.