EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 622 (February 2019) A&A, 622 (2019) A212 Full HTML
Open Access
Issue
A&A
Volume 622, February 2019
Article Number A212
Number of page(s) 10
Section Astronomical instrumentation
DOI https://doi.org/10.1051/0004-6361/201833843
Published online 21 February 2019

© ESO 2019

Licence Creative CommonsOpen 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 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 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 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 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 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 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 N2 × N2 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.

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

thumbnail 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: x-axis from 600 mm to 800 mm, y-axis 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.

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.

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

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 95 arcsec. The Strehl ratio, which quantifies the quality of the image, is also shown in this figure.

thumbnail 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 non-linear 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.

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:

(1)

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 Fr, Ft, 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 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 η.

thumbnail Fig. 4.

Three-dimensional 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 non-linear greyscale is used to make the curves more visible.

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.

thumbnail Fig. 5.

Modulus of the Fourier transform of the space variant PSFs (called pseudo-MTF) 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.

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

(3)

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:

(4)

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

(5)

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:

(6)

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 :

(7)

This writing clearly highlights the term related to the flux variation of the PSF, [HT1N]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 ill-posed 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 :

(8)

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

(9)

Consequently the matrix H is of dimensions N2 × N2.

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

(10)

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.

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 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 shift-invariant response, denoted there H(convolution) and the shift-variant response, denoted there H(Fredholm).

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

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 1020 for the convolution case and 6.7 × 1022 for the Fredholm case.

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

Computations are made on a toy object with the parameters of ASPIICS. Figure 8 shows, from top to bottom, the true object xt and the result of the product Hxt, giving the image y and a noisy image resulting from a Poisson transform of Hxt 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

(11)

where xt 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 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 number k̂ = 2000 for the noiseless case and k̂ = 27 for the noisy case.

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

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

thumbnail Fig. 10.

Cuts of the noiseless observed image y, the object xt and the reconstructed object using the RL algorithm. Top panel: noiseless image. Bottom panel: noisy image. See comments in the body of the paper.

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.

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

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 N2 × N2 points, here with N2 = 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 half-circle 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.

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

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

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

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.

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

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 Pk, 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 Pk, 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):

(12)

which has the advantage of preserving positivity. As described by Press et al. (1988-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π × RSun × 100, or about 6 million points sources ( since RSun = 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

(13)

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 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 N2 × N2 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π/Nt, where Nt 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. Nt = 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 Nt = 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

  1. Bayanna, A. R., Mathew, S. K., Sankarasubramanian, K., et al. 2011, Exp. Astron., 29, 145 [Google Scholar]
  2. Brueckner, G. E., Howard, R. A., Koomen, M. J., et al. 1995, Sol. Phys., 162, 357 [NASA ADS] [CrossRef] [Google Scholar]
  3. Cash, W. 2006, Nature, 442, 51 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  4. Daube-Witherspoon, M. E., & Muehllehner, G. 1986, IEEE Trans. Med. Imaging, 5, 61 [CrossRef] [PubMed] [Google Scholar]
  5. Demoment, G. 1989, IEEE Trans. Acoust. Speech Sig. Proc., 37, 2024 [CrossRef] [Google Scholar]
  6. Denis, L., Thiébaut, E., Soulez, F., Becker, J. M., & Mourya, R. 2015, Int. J. Comput. Vis., 115, 253 [Google Scholar]
  7. DePierro, A. R. 1987, IEEE Trans. Med. Imaging, 2, 328 [Google Scholar]
  8. Druckmüller, M., Habbal, S. R., & Morgan, H. 2014, ApJ, 785, 14 [NASA ADS] [CrossRef] [Google Scholar]
  9. Evans, J. W. 1948, J. Opt. Soc. Am., 38, 1083 [NASA ADS] [CrossRef] [Google Scholar]
  10. Flicker, R. C., & Rigaut, F. J. 2004, J. Opt. Soc. Am. A, 22, 3 [Google Scholar]
  11. Hanson, R. J. 1971, SIAM. J. Numer. Anal., 8, 616 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  12. Karush, W. 1939, Ph. D. Thesis, Univ. of Chicago, USA [Google Scholar]
  13. Koutchmy, S. 1988, Space Sci. Rev., 47, 95 [CrossRef] [Google Scholar]
  14. Kuhn, H. W., & Tucker, A. 1951, Proc. 2nd Berkeley Symp., (Berkeley: Univ. of California Press), 481 [Google Scholar]
  15. La Camera, A., Schreiber, L., Diolaiti, E., et al. 2015, A&A, 579, A1 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
  16. Lamy, P., Damé, L., Vivès, S., & Zhukov, A. 2010, SPIE, 7731, 18 [Google Scholar]
  17. Landweber, L. 1951, Am. J. Math., 73, 615 [NASA ADS] [CrossRef] [Google Scholar]
  18. Lantéri, H., Roche, M., & Aime, C. 2002, Inverse Prob., 18, 1397 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  19. Lantéri, H., Roche, M., Cuevas, O., & Aime, C. 2006, Sig. Proc., 81, 945 [CrossRef] [Google Scholar]
  20. Llebaria, A., Lamy, P., & Danjard, J.-F. 2006, Icarus, 182, 281 [CrossRef] [Google Scholar]
  21. Lyot, B. 1939, MNRAS, 99, 580 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lucy, L. B. 1974, AJ, 79, 745 [NASA ADS] [CrossRef] [Google Scholar]
  23. Lucy, L. B. 1994, A&A, 289, 983 [NASA ADS] [Google Scholar]
  24. Newkirk, Jr., G., & Bohlin, J. D. 1965, IAU Symp., 23, 287 [NASA ADS] [Google Scholar]
  25. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1988-1992, Numerical Recipes in C (Cambridge: Cambridge University Press) [Google Scholar]
  26. Purcell, J. D., & Koomen, M. J. 1962, Coronagraph with Improved Scattered-Light Properties, Report of NRL Progress (Washington, D.C: U.S. GPO) [Google Scholar]
  27. Richardson, W. H. 1972, J. Opt. Soc. Am., 62, 55 [Google Scholar]
  28. Rougeot, R., & Aime, C. 2018, A&A, 612, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Rougeot, R., Flamary, R., Galano, D., & Aime, C. 2017, A&A, 599, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Russ, J. C. 1995, The Image Processing Handbook, 2nd edn. (Boca Raton: CRC-Press) [Google Scholar]
  31. Shestov, S. V., & Zhukov, A. N. 2018, A&A, 612, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Strand, O. N., & Westwater, E. D. 1968, SIAM J. Numer. Anal., 287 [NASA ADS] [CrossRef] [Google Scholar]
  33. Theys, C., & Aime, C. 2016, Lect. Notes Phys., 914, 1 [NASA ADS] [CrossRef] [Google Scholar]
  34. Thiébaut, E., Denis, L., Soulez, F., & Mourya, R. 2016, Proc. SPIE, 9909, 99097N [CrossRef] [Google Scholar]
  35. Titterington, D. M. 1985, A&A, 144, 381 [NASA ADS] [Google Scholar]
  36. Vanderbei, R. J., Cady, E., & Kasdin, N. J. 2007, ApJ, 665, 794 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail 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: x-axis from 600 mm to 800 mm, y-axis 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.
In the text
thumbnail 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).
In the text
thumbnail 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 non-linear 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.
In the text
thumbnail Fig. 4.

Three-dimensional 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 non-linear greyscale is used to make the curves more visible.
In the text
thumbnail Fig. 5.

Modulus of the Fourier transform of the space variant PSFs (called pseudo-MTF) 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.
In the text
thumbnail 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.
In the text
thumbnail 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.
In the text
thumbnail 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.
In the text
thumbnail 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).
In the text
thumbnail Fig. 10.

Cuts of the noiseless observed image y, the object xt and the reconstructed object using the RL algorithm. Top panel: noiseless image. Bottom panel: noisy image. See comments in the body of the paper.
In the text
thumbnail 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.
In the text
thumbnail 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).
In the text
thumbnail 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.
In the text
thumbnail 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.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.