Translator Disclaimer
19 November 2018 Advanced wavefront reconstruction methods for segmented Extremely Large Telescope pupils using pyramid sensors
Author Affiliations +
The generation of Extremely Large Telescopes (ELTs) with mirror diameters up to 40 m has thick secondary mirror support structures (also known as spider legs), which cause difficulties in the wavefront reconstruction process. These spider legs create areas where the information of the phase is disconnected on the wavefront sensor detector, leading to pupil fragmentation and a loss of data on selected subapertures. The effects on wavefront reconstruction are differential pistons between segmented areas, leading to poor wavefront reconstruction. The resulting errors make the majority of existing control algorithms unfeasible for telescope systems having spider legs incorporated. A solution, named the split approach, is presented, which suggests to separate reconstruction of segment piston modes from the rest of the wavefront. Further, two methods are introduced for the direct reconstruction of the segment pistons. Due to the separate handling of the piston offsets on the segments, the split approach makes any of the existing phase reconstruction algorithms developed for nonsegmented pupils suitable for wavefront control in the presence of telescope spiders. We present end-to-end simulation results showing accurate, stable, and extremely fast wavefront reconstruction for the first light instrument mid-infrared ELT imager and spectograph of the ELT that is currently under construction.



The pupils on the Extremely Large Telescope (ELT)-generation telescopes are inevitably segmented and partially shaded by thick support structures due to their large sizes. Due to small subapertures, these secondary mirror support structures, also known as spider legs, lead to fragmentation of the pupil into several disjoint segments and cover partially or even completely up to several subapertures. If a wavefront sensor provides only local information about the wavefront, as the Shack–Hartmann (SH) does, differential pistons between the pupil fragments are not seen by the sensor (are in the null space), and therefore, cannot be controlled.1 The pyramid wavefront sensor (PWFS),2 due to its increased sensitivity, has been chosen as part of many instruments currently under the development for ELT-sized telescopes. In contrast to the SH sensor, the pyramid sensor does have a footprint of the segmented piston in its data (see Sec. 5.2.2 and Fig. 5). This suggests that the reconstruction of such modes is possible, if the reconstructor is able to extract this information from the sensor data and efficiently use it.

This paper pursues twofold aims. First, various existing reconstructors are analyzed and compared with respect to their ability to provide a stable control of the fragmented piston modes. This analysis exposes a serious lack of tools for the mentioned task, which makes the majority of available reconstruction algorithms helpless in the presence of telescope spiders. Therefore, the second aim of this paper is to provide such a tool.

Let us drop a few words on our framework. First of all, we distinguish between the so-called low-wind effect (LWE) and pupil fragmentation problem. The LWE is a real low-order distortion in the wavefront caused by the heating of the air in the vicinity of spider legs in particular low-wind conditions as observed on the Very Large Telescope.3 Pupil fragmentation is, on the contrary, a reconstructor-related phenomenon induced by discontinuity of sensor data. In this paper, we are not considering the LWE, but we focus on the differential pistons related to reconstruction methods only. The latter is caused by the inability of the reconstructor to provide an accurate solution of the wavefront based on fragmented data. Such a failure may have two origins. One reason is the inability of the wavefront sensor to provide the measurements of the jumps between the wavefront parts on pupil fragments, as is in the case with the SH sensor. Or it may as well be a deficiency in the reconstructor itself, which is not able to use in a proper way this kind of information that is though provided by the sensor. Moreover, in this paper, we focus on solutions to the pupil fragmentation problem, which do not require bringing any changes into the mechanical setup of the system, i.e., introducing any additional dedicated focal plane sensors that would measure the fragmented low-order modes during the AO loop, as suggested for the LWE in Refs. 45.6.7.

Let us consider the pupil fragmentation problem with PWFSs. We demonstrate in this paper that even if this type of sensor is able to provide the information on the jumps between the fragmented wavefronts, its usage does not automatically lead to correct reconstructions of segmented pistons with an arbitrary reconstruction algorithm designed for nonfragmented annular apertures. Based on numerous end-to-end simulations performed with various reconstruction algorithms (summarized in Secs. 3.5 and 4.2), we draw a conclusion that many of the known methods simply fail, meaning that the residuals contain randomly appearing uncontrolled piston modes on the segments. This leads to a significantly reduced correction quality in terms of Strehl ratio, PSF, and contrast. In addition to pupil fragmentation, another challenge that wavefront reconstruction algorithms for the ELTs have to tackle with is the high number of correcting elements that need to be controlled in real time.

Having these two points in mind, in Secs. 3 and 4, we give an overview of the available algorithms and their readiness/ability to operate both fast and stable with high-quality on ELTs. The analysis is performed in the context of the mid-infrared ELT imager and spectograph (METIS) instrument8 on the ELT containing a single conjugate adaptive optics (SCAO) system with 74×74 WFS subapertures. The geometry of the planned ELT M4 deformable mirror (DM) is taken into account. The thickness of the six secondary mirror support structures is 50 cm and coincides with the subaperture size.

On the one hand, as shown in Sec. 3.5, the current status of the performed simulations shows that among all possible variants of the interaction-matrix-based matrix–vector multiplication (MVM) methods only one, namely, the minimum variance reconstructor using a zonal interaction matrix and Laplacian regularization, is able to provide a stable control of the segmented low-order modes. However, the well-known drawback of any MVM method is the related computational load. Though for the considered SCAO system, the real-time application of MVM is still doable, it is hardly feasible in case of the planned extreme adaptive optics (XAO) system having tens of thousands of correcting elements to control.

On the other hand, as shown in Sec. 4, there exists a variety of interaction-matrix-free, model-based wavefront reconstructors developed for PWFSs in the recent decade, all of them being computationally much more efficient than interaction-matrix-based approaches. However, since those methods were developed such that they are intrinsically using the approximate forward models of the sensor not including segmented pupils, these methods fail when being applied straightforwardly to segmented sensor data. In Sec. 4.2, we additionally mention several attempts of adapting these methods for wavefront reconstruction on segmented pupils that unfortunately do not yield the expected performance. Therefore, in the presence of spiders, model-based wavefront control algorithms need to be adapted in some sophisticated way in order to handle the differential pistons between the pupil segments.

As a solution allowing for both high-quality and high-speed wavefront control, we suggest in Sec. 5 a hybrid scheme. This approach combines the advantage of the interaction-matrix-based methods of being able to handle pupil segmentation and the advantage of the model-based reconstructors of being fast. The solution, named the split approach, treats the reconstruction of segmented pistons separately from the higher-order modes (or frequencies). Here, the piston-free wavefront reconstruction on segments is provided by some model-based algorithm, e.g., the P-CuReD,911 as described in Sec. 5.1. In parallel, the segment pistons are reconstructed from the same sensor data with an interaction-matrix-based MVM approach. Since for the direct segment piston reconstruction we are only interested in the modes of order zero, the computational load can be significantly reduced.

In Sec. 5.2, we demonstrate two possibilities toward the formulation of the direct piston reconstruction for segmented pupils. The first one employs the usual setting of the full zonal interaction matrix using a set of dedicated basis functions representing the wavefront. This big matrix is then inverted via standard techniques of regularization and the resulting intermediate control matrix is afterward reduced to a small sized, control matrix relating the sensor data with the vector of segment pistons having only as many entries as pupil segments. Though still requiring the computationally expensive setting up and inversion of a dense matrix, which can be performed offline, the online calculations have linear complexity O(na) and are very cheap. In the second approach for direct segment piston reconstruction, the initial interaction matrix is formulated in the basis of segment pistons and is, therefore, very small from the start. Resulting in the same number of computations to be performed online, this approach is, in addition, free from the time-consuming (offline) operations involved in the first approach. Sec. 5.4 contains technical details on the implementation of the split approach.

In combination with the P-CuReD, both direct segment piston reconstruction methods have an overall computational effort, which scales linearly. In Sec. 6, we illustrate the performance of the split approach employing the two proposed methods for direct piston reconstruction for segmented pupils with end-to-end closed loop simulation results.



Let us denote the model of the pyramid sensor by an operator P:L2(R2)R2n, which maps real-valued L2 functions (wavefronts and residual phases) to a vector of discrete measurements of length 2n. The measurement process is given by

Eq. (1)

where ϕ describes the incoming phase, s is the pyramid sensor measurements, and η is the noise in the data. The inverse problem is to reconstruct the wavefront ϕ from given noisy sensor data s. Throughout this paper, we assume the pyramid sensor to operate in closed loop AO to be linear.

For a detailed description of the advantages of the PWFS over other sensor types and an overview on the instruments and telescopes, in which this type of sensor is already installed or is planned to be used in the nearest future, we refer the reader to Ref. 12. Here, we only mention briefly that due to the global response, the PWFS is able to sense the differential pistons, which has been successfully demonstrated in the laboratory,13,14 supported by numerical simulations, and validated on sky under seeing-limited conditions.15 The ability of the PWFS to sense the differential pistons of a segmented mirror and correct for it with an inversion based on the singular value decomposition (SVD) of the measured interaction matrix was first demonstrated in numerical simulations in Ref. 16. Among all WFS types tested in Ref. 17, the PWFS takes the most sensitive measurements of the differential pistons on the segments. Apart from that, compared to the SH WFS, the PWFS provides an increased sensitivity, which leads to higher limiting GS magnitudes and higher sky coverage.18

The reconstructors we will be dealing with in Secs. 3 and 4 are using different forward models of the PWFS or its approximations. Here, for the sake of brevity, we intentionally omit describing the forward models of the PWFS since we focus mainly on the achieved performance of various methods under pupil segmentation. For a detailed description of different model approaches and approximations of the WFS data we refer the reader to the devoted works.9,10,12,19,20

Note that the interaction-matrix-based MVM methods are described in Sec. 3 with a much higher level of detail compared to the model-based algorithms in Sec. 4. This is done intentionally by the authors who contributed to the development of model-based algorithms, which turned out to struggle under the circumstances of pupil fragmentation. The focus of this paper is to study, compare, and understand the behavior of different algorithms under pupil fragmentation. As a result of these efforts, this paper serves partially as a review of the currently available wavefront reconstruction algorithms for PWFSs and, in particular, of the present status in the performance those reach on the ELT-era instruments under design.

The connection between incoming wavefronts ϕ, residual wavefronts ϕres, and the mirror shape φ in an AO system is given by

Eq. (2)


Since for ideal compensation the residual wavefront ϕres should be equal to zero, the optimal choice for the DM shape is:

Eq. (3)


Therefore, for the control of DMs, one needs to know either the mirror actuator commands or the shape of the incoming wavefront provided as the solution of the inverse problem (1) for ϕ.


Interaction-Matrix-Based Reconstructors

In this section, we analyze the applicability and performance of the so-called interaction-matrix-based MVM methods for PWFS data and fragmented pupils. These methods involve a registration (or computation) of a WF-to-WFS interaction matrix, its inversion and a subsequent multiplication of the obtained control matrix with a vector of sensor data measurements. In the literature, there have been many variants of the interaction-matrix-based MVM approach presented: statistical estimation or solution in a least-squares sense; zonal or modal control approaches (i.e., the degrees of freedom are actuators/subapertures or modes). The presented overview aims to summarize and compare the performance of the existing interaction-matrix-based MVM methods in case of PWFS data fragmented by spiders. In Sec. 3.1, the generation of a WF-to-WFS interaction matrix is described, and the option of coupling or decoupling this step with the DM is explained in Sec. 3.2. In Sec. 3.3, the simplest least-squares approach and its regularized inversion are specified, whereas Sec. 3.4 deals with more sophisticated statistical approaches. Finally, the quality and speed performance of the interaction-matrix-based MVM algorithms are summarized in Sec. 3.5.


Generating the Interaction Matrix

We introduce (hi) as a set of arbitrary basis functions to represent the wavefront, (him) as a set of modal/global basis functions, and (IF) denotes the DM influence functions. In order to create the interaction matrix of the system, we need to relate the incoming wavefront with the output (measurements) of the PWFS. We represent the incoming wavefront ϕ using a set of basis functions (hi). Thus ϕ can be approximated by

Eq. (4)

where nc indicates the number of used basis functions. The interaction matrix MR2n×nc is then given by

Eq. (5)

i.e., the measurements

Eq. (6)

si=P(hi)for  i=1,2,,nc,
corresponding to the basis function hi build the i’th column of the interaction matrix.

The sensor equation, as already mentioned, reads as


To reconstruct the incoming (residual) wavefront ϕ, the matrix M has to be “inverted” and applied to the measurements represented by

with c=(ci)i=1,,nc.

After the reconstruction step, one has to derive the actuator commands a=(ai)i=1,,na from the reconstruction ϕ, i.e., solve

Eq. (7)


If the chosen basis hi coincides with the DM influence functions, the vectors c and a are equal.


Working with DM Influence Functions

Often, the interaction matrix inversion is coupled with the DM in the sense that for the generation of an interaction matrix one creates a certain (zonal or modal) shape with the DM, which is then sensed by the wavefront sensor. In this approach, one is restricted to wavefront shapes, which can be represented by the DM, i.e., are a linear combination of the DM influence functions:

Eq. (8)

or using DM modes

Eq. (9)

with actuator commands (mlj). This results in

Eq. (10)


Combining Eq. (1) with Eq. (8) or Eq. (9) produces a DM-to-WFS interaction matrix, which relates the sensor measurements s directly with the command vectors:

s=Pϕ+η=P(i=1naaiIFi)+η=i=1naaiP(IFi)+η=:  MIFa+η,
s=Pϕ+η=P(j=1nccjhjm)+η=j=1nccjP(hjm)+η=:  Mmc+η,
assuming a linear response of the pyramid sensor. In our approach, which is more general, the steps of wavefront reconstruction and projection to the DM are decoupled. One is not limited to using only the DM influence functions or modes, which can be provided by those as basis.


Deterministic Setting


Least-squares pseudoinverse

The least-squares problem

of finding the best wavefront fit cLS to the given WFS data vector s is uniquely solved by the least-squares minimum norm solution given as the Moore–Penrose generalized inverse

Such a pseudoinversion of the interaction matrix is considered to be the simplest reconstruction algorithm possible. Historically, least-squares with zonal representation was the first approach applied to the wavefront reconstruction problem involving the SH wavefront sensor.21 For wavefront reconstruction using PWFS measurements, the least-squares approach allows one to reach high correction accuracy without regularization, at least if the number of degrees of freedom is not very big. With the PWFS, the least-squares reconstructor has provided reasonable results on small scale systems having up to 30×30 subapertures within an 8-m telescope diameter as it was demonstrated in Ref. 22 with modal DM control using Karhunen–Loève (KL) polynomials.


Regularized least-squares pseudoinverse

For large-scale AO systems or more sophisticated configurations like multiconjugate adaptive optics (MCAO), the corresponding system matrices have large condition numbers and are difficult to invert. In this case, the conventional least-squares reconstructor performance is not satisfactory, and special treatment is required in the form of a regularization or filtering of unstable modes.

Sensor noise is modeled as a random process obeying zero-mean Gaussian statistics, ηN(0,Cη=σ2I), CηR2n×2n, where σ2 denotes the sensor noise variance. The noise-covariance-weighted least-squares (also known as minimum-norm maximum likelihood23) reconstructor, which minimizes

allows one to take the stochastic measurement uncertainties into account. The solution is given by

The pseudoinverse can be computed using the eigendecomposition of MTCη1M or the SVD of Cη1/2M. However, in practice, such decompositions are related with expensive computations and difficulties in computing the small singular values with the desired numerical precision. Also, small computation errors in eigen- or singular values lead to instabilities in the reconstruction due to noise amplification. As a regularization method, in the truncated SVD, one filters out the modes corresponding to singular values smaller than a given parameter α>0.

Alternatively, a similar effect is achieved by the so-called Tikhonov regularization, well-known in the field of inverse problems.24 For solving the inverse problem (1), we consider the least-squares problem


Then the regularized pseudoinverse is derived as

Eq. (11)

with Tikhonov regularization term consisting of a regularization parameter α>0 and identity matrix I.

For (relatively) small scale systems having up to 30×30 subapertures on an 8-m telescope, high correction quality has been achieved with both an SVD-regularized zonal25 and modal least-squares reconstructor using Zernike polynomials26 or KL polynomials.2729

However, due to noise propagation, the least-squares wavefront reconstruction algorithm performs poorly for large-scale or laser-guide star-based AO applications.30,31 Hence, there has been a tendency observed in the AO community to prefer the regularized variants of the interaction-matrix-based MVM method taking atmospheric statistics into account.


Bayesian Setting

Within a stochastic context, wavefront shapes and sensor noise are independent random processes obeying zero-mean Gaussian statistics, ϕN(0,Cϕ), ηN(0,Cη=σ2I), CϕRnc×nc, and CηR2n×2n, with σ2 denoting the sensor noise variance. Such a point of view allows one to use in the reconstruction the prior knowledge of the atmosphere and measurement noise statistics, expressed with the corresponding covariance matrices Cϕ and Cη, in order to regularize or stabilize the solution.

Two Bayesian statistical approaches, both using a prior probability density assumed on the phase, have been applied to the problem of wavefront reconstruction from sensor data—minimum variance estimation and maximum a posterior (MAP) estimation. The minimum variance [or minimum mean-square error (MMSE)] estimator minimizes the variance of the phase estimation error. The MAP estimator identifies the most likely value of ϕ given the observed data s and prior knowledge on the distribution of ϕ. Since wavefront reconstruction deals with zero-mean Gaussian signal and perturbation, the minimum variance reconstructor providing the minimal MSE coincides with the MAP reconstructor.23,32

In the Bayesian setting, the measurement vector s is a function of the atmospheric turbulence profile. The sensor Eq. (1) using the representation of the wavefront Eq. (4) formulated in terms of arbitrary coefficients ciR, in the stochastic setting is formulated in terms of wavefront coefficients ϕiR as


The aim in this setting is to compensate the turbulence-induced wavefront error.

The minimum variance/MAP reconstructor minimizes the penalized noise-weighted least-squares functional

which can be viewed as an estimator regularized with a Tikhonov term μ0ϕCϕ12. The weighting parameter μ0 allows one to balance between fitting to the data and the prior statistics. The corresponding regularized normal equation is
and the MAP reconstruction is given by
with control matrix MMAP

Eq. (12)


Note that besides the pseudoinverse, there exist several other methods based on the above normal equation for solving the inverse problem. The inverse of the phase covariance Cϕ1 must be chosen such that it is physically realistic. Typically, because of singularities (or ill-conditioning) in the turbulence spectra, it is inevitable to assume some discrete approximation on Cϕ1 and an additional regularization, which results in some loss of accuracy but yields stability.30

With the MAP/minimum variance estimators, two kinds of errors are related: the approximation error that tells how well the reconstructor approximates the inverse of the sensing operator P and the noise propagation error related to sensor noise. The sources of the model error are the chosen basis representation of the wavefront and the accuracy of the a priori statistical knowledge of the atmosphere.

In contrary to the least-squares approach, a statistical estimation method necessarily needs regularization parameter tuning for an accurate wavefront reconstruction. The numerical simulations indicate that the MAP/minimum variance reconstructor with an optimized parameter μ0 performs better (is more stable) than the (noise-weighted) least squares solution.33,34

In the following, we briefly focus on the implementation details of two variants of Bayesian reconstructors that proved to be efficient in astronomical AO.


MAP reconstructor with modal control of the DM

where the coordinates (x,y) describe a point in the pupil plane, IFl describes the influence function of the l’th actuator, and the integration is performed over the pupil mask domain Ω. Afterward, a KL polynomial-based interaction matrix Mm is constructed by applying to all DM actuators the theoretically predefined mi commands, i.e., they form the closest shape to the desired KL polynomial. The numerical approximation him to the theoretical mode Kidm is represented as linear combination of DM influence functions:
for coefficients mliR.

In this method, Cϕ is the von Karman wavefront covariance matrix restricted to the modal space of the DM. The parameter μ0, which scales as the inverse of the square of the signal-to-noise ratio, allows one to weight the sensor noise and atmospheric priors in a flexible way.

This reconstructor is implemented in Octopus.3436


Zonal minimum variance estimator with regularized sparse approximation of Cϕ

Another variant of MVM that provides high-quality reconstruction is the zonal minimum variance estimator with regularized sparse discrete approximation of Cϕ1, as suggested in Ref. 30, with

where L denotes a discrete Laplacian matrix approximating the Laplacian operator. The constant c0 is physically interpreted as the strength of the turbulence and additionally normalizes the approximation of Cϕ in order to fit the von Karman turbulence spectrum.30,37

Using this covariance approximation corresponds to regularization by the l2-norm of the Laplacian, i.e., to solve the penalized least-squares functional:

with regularization term μ0c0Lϕ2 that removes waffle mode and other high-frequency errors in the phase estimates.37

This reconstructor using DM influence functions as basis functions is incorporated in YAO, an open-source AO simulation tool written in yorick.38 For this study, we have also implemented it in Octopus. In the numerical implementation, stabilization with respect to the wavefront sensor noise is performed by filtering out the columns in the registered interaction matrix corresponding to basis functions whose WFS response was smaller than a certain predefined value (e.g., expressed in a percentage of the maximum registered response). This parameter can be tuned in order to optimize the performance of the reconstructor for different flux settings or spiders’ thickness. Note that in the numerical simulations, the parameter c0 is heuristically tuned as well for various atmosphere strengths and guide star fluxes.


Performance of Interaction-Matrix-Based MVM in the Presence of Spiders



In this section, we summarize and compare the available results of various versions of the interaction-matrix-based MVM methods in the presence of spiders. First, we mention the earlier result achieved in the context of EPICS, the XAO instrument on the ELT in the case of four thick spiders.39 It was reported therein that both the zonal and modal reconstructors provide the same quality in the presence of four spider legs. Apart from that, two important points were underlined: the light behind spiders needs to be used in the reconstruction; and the amount of modulation should not be too large, otherwise, the sensor loses its sensitivity to low-order modes.

Next, we report in Table 1 on the performance we obtained with two variants of Bayesian reconstructors in the case of six thick spiders in the context of the METIS instrument on the ELT. The results are presented for a median atmosphere, on-axis correction and high photon flux. Note that the modal MAP result from Octopus and the zonal reconstructor results from YAO were obtained using a regular Fried geometry, whereas the zonal method results from Octopus imply the real ELT M4 geometry.

Table 1

LE Strehl ratio obtained with three versions of the interaction-matrix-based MVM for pupils with and without spiders.

MethodZonal MMSE (OCT.)Zonal MMSE (YAO)Modal MAP (OCT.)
No spiders, Fried geometry0.89400.89
No spiders, M4 geometry0.89
With spiders, Fried geometry0.89400.86
With spiders, M4 geometry0.89

The zonal minimum variance reconstructor, implemented both in Octopus and YAO, provides the best quality of 0.89 LE Strehl, which barely decreases for fragmented pupil masks. For more YAO simulation results obtained for the METIS case with the zonal minimum variance reconstructor, we refer to Ref. 40.

The modal MAP reconstructor (as implemented in Octopus) achieves the same LE Strehl ratio of 0.89 in the case without spiders. However, in contrary to the zonal approach, its performance in the presence of spiders is not as good. After running multiple tests over a set of regularization parameters, the best LE Strehl ratio that we were able to achieve is 0.859. In the residual screens, we always observed uncompensated random segmented pistons. Note that in order to obtain reasonable results one has to set the illumination parameter to some value lower than 0.5, i.e., to use data from the subapertures partially covered by the spiders. If one does not use this information, the reconstruction gets unstable because of the poorly corrected segmented pistons.

Similar results with the MAP reconstructor have been recently reported in Ref. 41 for an 8-m telescope. First, the authors therein underline the importance of using the light from the subapertures partially covered by the spiders, in the same way that we observed in our tests. Second, the quality of the spider-free case could not be achieved in the presence of spiders. The authors in Ref. 41 conclude that the modal MAP reconstructor as it is known is not able to overcome the fragmentation problem completely. It is presumed that a dedicated change of basis functions, which would be defined on each segment instead of the entire pupil, could potentially be helpful, but no attempts have been performed so far in this direction.

An intermediate conclusion is that only one version of the interaction-matrix-based methods, namely the zonal minimum variance reconstructor, handles the pupil fragmentation problem without almost no loss in quality.



Interaction matrix inversion is related to a high computational load, which can, however, to a large extent be performed offline. Assuming that an interaction-matrix-based MVM algorithm is implemented using conventional matrix inversions and MVMs, the method requires approximately O(na3) operations to compute the control matrix and O(nna) operations to apply the control matrix.30 The computational complexity becomes a limitation of these approaches if the values of na are of order 10,000 as currently under consideration for XAO systems on ELTs.


Advanced Model-Based Reconstructors

Apart from the interaction-matrix-based methods, there exists plenty of recently developed algorithms based on the mathematical analysis of the forward models of the PWFS. The main feature of these algorithms is a low computational complexity resulting in an adequate handling of wavefront reconstruction on AO systems in real time and still guaranteeing high-quality and robustness of the methods. A detailed overview and comparison of various algorithms can be found in Ref. 42. Until lately, these methods have been extensively tested with nonfragmented pupils. The question is: are they also able to provide reasonable reconstructions with fragmented pupil masks?

Keeping the structure similar to the one in the previous section, let us first start with a brief overview of the methods. Note that, for the sake of brevity, these are presented without any mathematical details, which can be found in the provided references. Then we will have a look at quality and speed these methods provide for nonfragmented pupils. Finally, we discuss the performance of the model-based methods in the presence of spiders.

We consider the following algorithms:

  • preprocessed cumulative reconstructor with domain decomposition (P-CuReD)911

  • convolution with the linearized inverse filter (CLIF)9,12,43

  • pyramid Fourier transform reconstructor (PFTR)12,43

  • finite Hilbert transform reconstructor (FHTR)9

  • singular value type reconstructor (SVTR)44

  • a number of linear iterative methods applied to pyramid sensors like conjugate gradient for the normal equation (CGNE), steepest descent (SD), Landweber iteration, and Kaczmarz versions of the previous mentioned algorithms4547

  • nonlinear Landweber and Landweber–Kaczmarz iteration.48

The first three methods, the P-CuReD, the CLIF, and the PFTR are all using Fourier domain analysis of the underlying (modulated and nonmodulated) pyramid sensor operators and their approximations. The CLIF and PFTR suggest different numerical implementations of the same idea of reconstructing the wavefront spectrum from the data spectrum using the known Fourier domain relation between both. The PFTR applies the corresponding inverse Fourier domain filter to the data spectrum, the CLIF algorithm works in the spatial domain. The P-CuReD correlates the pyramid with the SH sensor and by employing the Fourier domain relation between both, suggests a simple convolution-type conversion of pyramid sensor data to SH-like data. Subsequently, the CuReD is applied to the preprocessed pyramid data, a method that was originally developed for the SH-WFS.4951

The FHTR and the SVTR are based on a simplification of the nonmodulated pyramid sensor model represented as the finite Hilbert transform of the incoming wavefront. The FHTR applies a direct inversion formula of the finite Hilbert transform, the SVTR uses a version of the involved operators’ SVD for wavefront estimation. In contrast to MVM approaches, the algorithms consist of analytical derivations of the mentioned tools for inversion.

In the iterative methods, we use well-known mathematical algorithms for solving the inverse problem of wavefront reconstruction from pyramid sensor data, namely CGNE, SD method, SD-Kaczmarz method, Landweber iteration, and Landweber–Kaczmarz iteration. All algorithms are applicable to pyramid sensors with and without modulation, have comparable numerical complexity, and provide similar reconstruction quality. Outstanding are the Kaczmarz versions of the algorithms, which enable a sophisticated combination of the two data sets sx and sy for reconstruction. Although all mentioned model-based reconstructors first were introduced based on a linearization of the pyramid sensor model, for Landweber and Landweber–Kaczmarz iteration, we additionally considered their nonlinear versions in order to develop nonlinear reconstructors.


Quality and Speed Performance Without Spiders

For a detailed comparison of the methods with respect to the reconstruction quality, we refer to the corresponding publications. Here we only briefly mention that in closed loop simulations all the model-based methods achieve a quality close to that of interaction-matrix-based approaches with some of them even slightly outperforming the latter. In Table 2, we provide a comparison of wavefront reconstruction algorithms for pyramid sensors in terms of their computational complexities. From this table, one can see that all recently developed model-based reconstruction algorithms require much fewer computations to be performed than an MVM implementation. The P-CuReD algorithm reaches the highest quality results in a variety of tested conditions, and at the same time has the lowest, namely linear, computational complexity. In Refs. 52 and 53, AO simulation tools’ users compared the performances of the modal MAP method with the P-CuReD for XAO settings. In Ref. 53, the MVM and P-CuReD algorithm give the same reconstruction quality with only very slight discrepancies in some of the test cases. Moreover, in Refs. 911, it was shown that the P-CuReD provides a significantly improved quality in the low-flux cases compared to an interaction-matrix-based reconstructor, and also that it converges to high Strehl ratios faster than the tested MVM approach.

Table 2

Comparison of the currently existing algorithms for wavefront reconstruction from pyramid sensor data in terms of their flexibility and computational complexity.

Interaction-matrix-based MVMO(na2)
FHTR and SVTRO(na3/2)
Iterative methodsO(na3/2)

Another point we would like to stress here is that all our model-based reconstructors are free from the time-consuming precomputation of matrices and fine tuning of the regularization parameters associated with MVM approaches. Since there are no intrinsic regularization parameters necessary to consider during the reconstruction process, no optimization is needed if atmospheric conditions change. This facilitates the usage of the methods by external users, as confirmed by Refs. 52 and 53 with P-CuReD applications and Ref. 54, where the CuReD method, which was originally developed for the SH sensor49,50 and nowadays also constitutes a part of the P-CuReD method for pyramid sensor, has been successfully tested on-sky.


Adapting Advanced Reconstructors to Segmented Pupils

In this section, we describe our first approaches to overcome the difficulties caused by pupil fragmentation or spider effects when the width of the spiders is no longer negligible compared to the size of the subapertures. Numerical investigations were performed for the ELT having a primary mirror divided into six segments as shown in Fig. 1. As reconstruction method, we use the P-CuReD algorithm. Almost as important as the width of the spider legs is their placement. If the spiders are in parallel with the x and y axes, some basic attempts described in the following succeed. In contrast, for arbitrarily located spiders, we examine poor wavefront reconstruction suffering from differential piston influences. Note that the following described attempts for wavefront reconstruction in the presence of six telescope spiders were only investigated with application of the P-CuReD algorithm. Different reconstruction algorithms will behave differently.

Fig. 1

ELT pupil mask with spiders. The ELT will consist of a 39.3-m-diameter primary mirror and a 11.1-m-diameter secondary mirror, which will be supported by 6 spiders each being 50-cm thick.


One idea is to make the illumination factor necessary for the usage of subapertures for wavefront reconstruction very low, i.e., utilizing (almost) all available subapertures for reconstruction, also those being less illuminated. This method combined with the reconstruction algorithm P-CuReD does not lead to success since the light suffers from obstruction effects especially at the boundaries of the pupil segments. The wavefront reconstruction is adversely influenced on the whole pupil by differential piston effects as seen in Fig. 2 for different illumination factors.

Fig. 2

Quality results for different illumination factors. The SE Strehl ratios for wavefront reconstruction using the P-CuReD without any specific telescope spider handling. The three lines indicate the usage of more or less illuminated pixels for WF reconstruction. The reason for the generally poor performance is segmented piston modes which go out of control. Despite the poor results (with LE Strehl ratios of less than 0.5 instead of   0.9 for spider-free simulations), one can see that an illumination factor slightly under 50% is preferable as already discussed in Sec. 3.5.1 for the modal MAP reconstructor.


Another approach to overcome the effects of pupil segmentation is data interpolation or interpolation of the reconstructed phase under the spider legs. Instead of using defective measurements under the legs provided by the pyramid sensor, we generate data or reconstructions artificially in these areas. For that purpose, we considered bilinear and spline interpolation. Unfortunately, these approaches do not eliminate differential piston effects.

A more sophisticated attempt incorporates the pyramid sensor model in the measurement continuation process by applying an iterative measurement extension method already presented in Refs. 44 and 55. The idea in Ref. 44 was originally developed for wavefront reconstruction on annular, nonsegmented telescope pupils. The basic concept is to generate artificial but pyramid related data under the spider legs by application of the finite Hilbert transform, which is a simplification of the Fourier optics-based pyramid sensor model. Thus the provided data correspond better to pyramid measurements as it is the case for a general interpolation procedure described above. However, the simulation results are quite similar to the approach using bilinear or spline interpolation, and differential piston effects are developing within time and make the reconstruction poor.

An additional experiment is to replace the obstructed data by zeros, i.e., it is assumed that the wavefront is planar in the areas obstructed under the spider legs. Note that zero padding is successfully used for the central obstruction induced by a secondary mirror for some of the reconstructors mentioned in this section, e.g., the SVTR. In contrast to obstruction induced by spiders, a central obstruction does not cause segmented mirrors and hence does not induce differential pistons. Several numerical simulations show that the results of this approach differ significantly for different spider leg locations. In particular, if one considers four spider legs parallel to x and y axes (as those first investigated for SH sensors in Ref. 1), the zero padding approach gives satisfying reconstruction quality with the pyramid sensor. Considering the six ELT legs, the approach again suffers from differential piston development as shown in Fig. 3.

Fig. 3

Residual screen in radians (evaluated in K-band) for the “padding with zeros under the spider legs” approach. Attempts like measurement continuation, interpolation of data, or the reconstructed phase as well as reconstruction using the light under the support structures deliver similar poor results.


Approaches like jump minimization between the segments by boundary integral coupling do bring quality improvements but not as high as we have hoped for. For SH sensors, this method will precisely be described in an upcoming paper of our group. Actuator coupling/slaving of those actuators that are situated at the boundary of spider legs is a hot topic for pyramid sensors on ELTs, especially for wavefront sensing at shorter wavelengths than K-band.

Altogether, the above-described methods do not satisfactorily handle the impact of 1 or 2 subaperture thick spiders on the wavefront correction performance. The negative influences of randomly appearing piston modes on the wavefront reconstruction gave us the motivation to develop ideas to separately examine the segmented pistons. After trying all the described methods, a way for successfully eliminating spider obstruction effects was found. The idea is to reconstruct the segmented piston modes separately from other frequencies in the wavefront. Within these attempts, which we recap as split approach, one still can use the fast interaction-matrix-free wavefront reconstruction algorithms presented in this section and at the same time obtain stable wavefront correction in the presence of spiders.

In the next section, we introduce the basic concept of the split approach and two different algorithms tested for direct segment piston reconstruction. Due to the high-end performance of the P-CuReD, both in terms of quality and speed, we choose this method as one of the components in our split approach aimed at segment-piston-free reconstruction of the wavefront for segmented pupils.


Split Approach

In this section, we describe robust algorithms that allow to compute optimal mirror configurations from signals containing telescope spider obstruction by dividing the wavefront reconstruction into two parts:56

  • 1. Piston-free wavefront reconstruction on each segment: This wavefront reconstruction method handles all modes seen by the PWFS except modes of order zero (piston modes) on each segment.

  • 2. Direct segment piston reconstruction: Here, we focus on modes of order zero on each segment solely.

The reason why we suggest to split piston reconstruction from the full phase reconstruction is twofold. On the one hand, we want to make stable wavefront reconstruction in the presence of spiders feasible with the computationally efficient model-based algorithms presented in Sec. 4. Those algorithms were developed using the forward model of the sensor for annular apertures that do not include segmented pupils and spider effects. Straightforward attempts as described in Sec. 4.2 of applying these algorithms to sensor data “spoiled” by spider obstruction failed so far. On the other hand, it was recognized that the interaction-matrix-based approaches, as described in Sec. 3, are able to correct for differential pistons but are very time-consuming. The related computational effort may be affordable for the METIS system but hardly feasible for XAO systems. Therefore, our goal was to combine the P-CuReD (or any other model-based wavefront reconstruction method) and the advantages of interaction-matrix-based reconstruction to obtain a fast and robust reconstruction approach for segmented pupils having less computational complexity than an MVM.

We recall that P describes the pyramid sensor operator, ϕ is the incoming phase, and s is the pyramid sensor measurements. Usually, all wavefront modes (frequencies) that are seen by the PWFS and afterward corrected are treated within the same wavefront reconstruction process. In the split approach, we separate the incoming wavefront ϕ into the parts:

Eq. (13)

where ϕi indicates the piston-free wavefront reconstruction on segment Ωi. The corresponding reconstruction procedure denoted by P˜ is provided by any of the existing model-based algorithms (mentioned in Sec. 4), which provide high-quality reconstruction on each segment. The term pi describes the corresponding piston information on every segment for i=1,2,,k calculated independently using direct segment piston reconstruction methods that will be described in Sec. 5.2 and are in the following denoted by Π.

Hence, we separate the whole wavefront reconstruction into two parts:

Eq. (14)

which is feasible in closed loop AO for an almost linear pyramid sensor response.19,57


Piston-Free Wavefront Reconstruction on Segments

In this section, we focus on piston-free wavefront reconstruction on segments. For the estimation of the wavefront ϕpistonfree, we can use any of the existing fast model-based wavefront reconstruction algorithms mentioned in Sec. 4. These methods have shown exceptional wavefront correction quality on nonsegmented pupils, i.e., on the full annular telescope aperture in the following denoted by Ω. The spider legs divide the aperture into segments meaning that for k spider legs we obtain k disjoint segments indicated by Ωi,i=1,2,,k (for instance the six spider legs of the ELT shown in Fig. 1). For segmented pupils, the wavefront reconstruction method P˜—in addition to standard wavefront reconstruction requirements—fulfills two conditions:

  • 1. The method is implemented on segments Ωi, i=1,2,,k instead of the full mask Ω.

  • 2. The reconstruction ϕi on every segment Ωi, i=1,2,,k needs to be piston-free.

Note that condition 2 does not constitute a restriction since one can always compute the local piston information of a full segment from the reconstructed wavefront and subtract it afterward. Condition 1 may possibly be attenuated to an algorithm being implemented on the full aperture but dividing the reconstruction into segments thereafter. In this case also, the elimination of zero-order modes on the segments can be performed separately from the reconstruction process. However, we clearly want to point out that we did not investigate this idea in detail. Until now, we only used a segment-piston-free reconstruction method implemented on the segments. An analysis of these considerations will be part of a subsequent study.

For the segment-piston-free wavefront reconstruction, we use the P-CuReD911 applied to each segment. With linear computational complexity, this algorithm is the fastest method available for wavefront reconstruction from pyramid sensor data and, at the same time, it provides a reconstruction quality close to the theoretical limits.


Direct Segment Piston Reconstructors

Now we introduce two methods for direct segment piston reconstruction denoted by Π providing the piston information on the disjoint segments (pi)i=1,2,k divided by k spider legs. All of them follow the idea of an interaction-matrix-based reconstruction using MVM. In contrast to the conventional approach, we now do not focus on the reconstruction of the complete wavefront, but only on the reconstruction of piston modes on segments. The first method (see Sec. 5.2.1 for details) uses a zonal interaction matrix containing the sensor response to every single zonal basis function. The second one (see Sec. 5.2.2 for details) is rather modal with a very small number of modes used as basis. Namely, we use a segment basis consisting of k modes, where each mode depicts a piston on a given segment.


Direct segment piston reconstructor I: single-poke-approach

We start with wavefront reconstruction using a full zonal interaction matrix as already described in Sec. 3 and transform this algorithm to the reconstruction of segment pistons solely.

For the so-called single-poke-approach, we measure the response of the pyramid sensor for every basis function. More precisely, we use the full interaction matrix M of the system as described in Eqs. (5) and (6) computed for a set of zonal basis functions. The amplitude corresponding to the mentioned basis functions has to be small to ensure a linear response of the pyramid sensor. The control matrix MαRnc×2n is created as in Sec. 3.4.2 using the squared Laplacian as a regularization term.

Since we are only interested in the reconstruction of segment pistons and omit the reconstruction of other modes, the dimension of the problem is drastically reduced. Assume we have the complete wavefront ϕ=Mαs reconstructed. Extraction of the segment piston information p=(pi)1ik from the known wavefront ϕ is obtained by averaging the phase values within each segment. This step is modeled as a multiplication of ϕRnc with a matrix QRk×nc, where the i’th row of Q contains a vector representation of the segment Ωi divided by the number of active subapertures on Ωi (for the averaging). The application of the matrix Q leads to

Eq. (15)

with a dense but small matrix Π1Rk×2n. This means the piston information on the segments is reconstructed from the given sensor data with a linear computational complexity. Hence, we reduce the computationally expensive full interaction matrix approach with complexity O(nna) to a cheap direct segment piston reconstruction method. For usage of the P-CuReD within the split approach, the partition of wavefront reconstruction into separate piston and higher-order frequencies reconstruction only slightly decreases the speed of the reconstruction method, which is still faster as the full interaction matrix approach. Of course, the interaction matrix of the system still needs to be set up for the application of this direct segment piston reconstruction method but these calculations are done offline.

In Fig. 4, we provide an illustration to an Octopus simulation using the split approach with the P-CuReD and the first direct segment piston reconstructor (DSPR). Here we focus on the reconstructions in the initial loop steps, when the sensor is yet far from its linear regime and the loop is only getting closed. One can see that the reconstruction is reasonable and provides a stable convergence. The corresponding closed loop results for a longer simulation are presented in Sec. 6.

Fig. 4

Illustration of wavefront reconstruction with the split approach at the third time step in an AO loop performed in Octopus. The first line shows an incoming screen given in radians (K-band) and the corresponding pistons on the six segments. The second line shows the segment-piston-free wavefront reconstruction using the P-CuReD on segments and the direct segment piston reconstruction described in Sec. 5.2.1. In the last line, we see the combination of both, i.e., the whole wavefront reconstruction, and the difference between the incoming screen and the final reconstruction. Even though the pyramid sensor is not in its linear regime yet, the algorithm already gives reasonable reconstructions of the piston-free wavefront as well as the segment pistons. The simulation parameter corresponds to the ones specified in Sec. 6.1.



Direct segment piston reconstructor II: segment-poke-approach

Inspired by the simple relation Eq. (15) connecting segment piston values p with the sensor data s through a small-size matrix Π1, we want to formulate another direct segment piston reconstruction approach, which will allow us to skip the computationally expensive and time-consuming step of setting up the full interaction matrix of the system. Rewriting Eq. (15) as

we see that, formally, it is possible to define a small interaction matrix of the system in the basis consisting of only a few segment pistons pi, i=1,,k.

We again consider arbitrarily located telescope spiders having k spider legs, i.e., dividing the aperture into k disjoint segments Ωi, i=1,2,,k. The effects of a piston offset with amplitude c on a single segment are described by

si=P[c·XΩi(x,y)]for  i=1,2,,k,
XΩi(x,y):={1,for  (x,y)Ωi0,else,
denoting the characteristic functions of the telescope aperture segments Ωi. Again, we assume that the pyramid sensor fulfills the linearity assumption if c is chosen small enough.

With the measurement vectors siR2n for every segment Ωi, i=1,2,,k, we obtain a piston interaction matrix MpR2n×k consisting of


The pyramid response to a segment poke is shown in Fig. 5. While in the previous method a very time-consuming precomputation of the full interaction matrix is required, the matrix necessary for this approach only has dimension 2n×k, and hence its computation is much cheaper. Therefore, setting up the matrix containing the WFS response to segment piston modes can be recomputed quickly for changing seeing conditions and readily performed online.

Fig. 5

Pyramid sensor measurements for a single-segment piston mode of height 5×108  m. The first column shows the segment piston, the second and third column illustrate the corresponding sensor measurements in x and y directions. The sensor measurements are shown in the range [1.5×108  m, 1.5×108  m].


Now the piston reconstruction on every segment is described as minimization of

for some suitable chosen Tikhonov matrix Γ and regularization parameter α>0. Solving the equation in a least squares sense leads to the normal equation:
where the right-hand side represents a projection of the measurements containing all modes seen by the pyramid sensor s to data MpTs including piston information only. Using Tikhonov regularization, we obtain the piston control matrix Π2Rk×2n for the direct piston reconstruction on segments. The segment piston pRk are then obtained by


Fast and Robust Wavefront Reconstruction Under Pupil Segmentation Using the Split Approach

The general scheme of the split approach for wavefront reconstruction using one of the above-introduced direct segment piston reconstruction methods is described by Algorithm 1.

Algorithm 1

Split approach.

choose segment-piston-free wavefront reconstruction method P˜
choose direct segment piston reconstruction method Π
for all time steps do
 get measurements s
for segments i=1,2,,k
  substract global piston for every ϕi if necessary
end for
end for

The segment piston control matrices Πi, i=1,2 are dense but only of dimension k×2n, which leads to an optimization of the proposed approaches with respect to computational complexity. The expensive steps can be precomputed offline, thus, the algorithms scale linearly. This is a clear advantage to the computationally expensive wavefront reconstruction using full interaction-matrix-based MVM approaches as indicated in Fig. 6. For a moderately large-scale SCAO system such as METIS, the gain in the computational efficiency provided by the split approach is of order 104. For an extremely large-scale SCAO system as the planned XAO system, the corresponding gain is of order 105.

Fig. 6

Logarithmic plot of the approximate complexities of the MVM methods using the full control matrix (dotted line) and of the split approach using any of the described methods for direct segment piston reconstruction applying the small-size control matrices (solid line).



Details on the Realization

At the end, we will focus on a few numerical implementation details of the split approach, namely the incorporated pyramid sensor model, illumination factor, loop gain, and phase ambiguity.

For the proposed direct segment piston reconstruction methods, the nonlinear PWFS model including interference effects can be taken into account for the calculation of the piston control matrices Πi, i=1,2.

Our numerical simulations and those performed in Ref. 41 show that for wavefront estimation under pupil fragmentation it is important to choose an appropriate illumination factor determining that subapertures are active and therefore used for wavefront reconstruction. As already discussed, we perform the segment-piston-free reconstruction on disjoint segments, hence we only use subapertures that are illuminated at least 75%. In contrast, for the DSPR, usage of measurements on less illuminated subapertures is crucial.

The loop gains for the two parts of the split approach have to be considered separately. We identify an optimal loop gain for the segment-piston-free wavefront reconstruction using the P-CuReD and another one for the direct segment piston reconstruction. At the beginning of a closed loop simulation, it is crucial to avoid phase ambiguity caused by the sinusoidal part of the pyramid data, meaning that piston offsets of size 2π radians in the phase cannot be distinguished by the pyramid sensor but heavily influence the image quality.58 Due to the nonlinearity of the PWFS (“optical gain”) and the need to correct low frequencies fast, we use a higher integrator control loop gain for the DSPR in the first iterations. This results in a stronger emphasis on the correction of low-order modes and provides an adequate control of piston offsets for data corresponding to larger phases.

The last step of the algorithm contains the projection of the reconstructed wavefront on the DM. One can, therefore, either solve Eq. (7) or simply evaluate ϕ at the actuator positions. We did choose the second approach. To be able to control also actuators outside the reconstruction area, we smoothly extend the reconstruction ϕ to a larger domain covering all used DM actuator positions.


Numerical Validation


Simulation Environment and Parameters

We use the simulation tool Octopus35,36 provided by the European Southern Observatory and study the METIS instrument having an SCAO system incorporated. As one of the three first light instruments of the ELT, METIS will allow investigations of exoplanets concerning physical and chemical properties like weather, temperature, seasons, or the composition of their atmospheres. It will, among others, focus on proto-planetary disks, the formation of planets, and the Solar System as well as the growth of supermassive black holes. The science wavelengths of the instrument range from 3 to 19  μm, the wavefront sensing is performed at 2.2  μm. For the METIS instrument, an annular mask is currently considered with an outer diameter of 37 m and an inner diameter of 11.1 m (30%), i.e., all edges of the real ELT primary mirror having a diameter of 39 m are cropped such that we obtain a circular area. As a sensing device, we simulate a PWFS having 74×74 subapertures corresponding to a subaperture size of 0.5 m. The simulation grid size is selected as 740×740 pixels on the aperture resulting in a resolution of 0.05 m per pixel for a 37-m telescope.

The phase screens of our simulation are based on the von Karman atmospheric model having 35 layers and an outer scale of L0=25  m. The layer model of the atmosphere ranging from 30 m to 26.5 km above the observatory’s platform was provided by the European Southern Observatory. Two seeing conditions were simulated: median with the Fried parameter r0 at 0.5  μm=0.157  m and very bad with r0 at 0.5  μm=0.097  m. We consider two different photon fluxes of 600 and 100 incident photons per subaperture per frame corresponding to a guide star magnitude of 8 or 10, respectively.

A key element of every AO system is the DM. It provides fast steering capabilities to compensate for wavefront aberrations caused by atmospheric turbulence and telescope perturbations in real time and hence allows one to optimize the telescope performance. The actuator positions of the M4 DM used in our Octopus simulations are shown in Fig. 7. We use the M4 influence functions internally incorporated in Octopus. Altogether, we have 3874 active subapertures and 5190 active actuators in use. All test case parameters are summarized in Table 3.

Table 3

Simulation parameters used in the tests.

Simulation parameters
Telescope diameter37 m
Central obstruction11.1 m, i.e., 30%
Pupil maskELT pupil with 50-cm wide spiders
Pupil segments6
Science targetOn-axis (SCAO)
Sensing bandK (2.2  μm)
Evaluation bandsK (2.2  μm)
L (3.7  μm)
N (10.0  μm)
Atmospheric modelvon Karman
Number of simulated layers35
Outer scale L025 m
Fried radius r0 at λ=500  nm0.157 m, 0.097 m
Coherence time τ0 at λ=500  nm5.35 ms
Number of subapertures74×74
Minimum subaperture illumination45%
Number of active subapertures3874 out of 5476
Linear size of simulation grid740 pixels
DM geometryELT M4 model
DM delay1 frame
Number of active actuators5190
Detector read-out noise1 electron/pixel
Background flux0.431 photons/pixel/frame or 4.19  ADU/m2/ms/arcsec2
Frame rate500 Hz
Photon flux or guide star brightness[100, 600] photons/subaperture/frame or [8, 10] mag
Zenith angle[0 deg,30 deg,60 deg]
Simulation time2 s (1000 iterations)

Fig. 7

Actuator positions of the ELT M4 geometry implemented in Octopus. The spider boundaries are indicated as lines.



Numerical Results

A numerical performance analysis and a comparison of the approaches in terms of reconstruction quality for the above-described simulation environment are provided in Tables 4Table 56 and Fig. 8. The three tables present the long-exposure Strehl ratios evaluated at three different wavelengths λ=2.2  μm, λ=3.7  μm, and λ=10.0  μm correspondingly. In each of these three tables, the test cases cover two photon fluxes of 600 and 100 photons/subaperture/frame and three zenith angles of 0 deg, 30 deg, and 60 deg. Each table demonstrates the results obtained with the P-CuReD for a spider-free simulation as benchmark and the spider effected results obtained within the split approach using each of the DSPR.

Table 4

Long-exposure Strehl ratios in K-band (2.2  μm) after 1000 closed loop simulation steps. As reconstruction method, we use the P-CuReD. We compare the results for a benchmark simulation without spiders using the P-CuReD only and employing the split approach combined with the proposed DSPRs in the presence of telescope spiders.

Photon fluxZenith angle (deg)No spidersDSPR IDSPR IITheoretical limit

Table 5

Long-exposure Strehl ratios in N-band (10  μm) after 1000 closed loop simulation steps using the P-CuReD in all simulations.

FluxZenith angle (deg)No spidersDSPR IDSPR IIESO goalESO requirementTheoretical limit

Table 6

Long-exposure Strehl ratios in L-band (3.7  μm) after 1000 closed loop simulation steps using the P-CuReD in all simulations.

FluxZenith angle (deg)No spidersDSPR IDSPR IIESO goalESO requirementTheoretical limit

Fig. 8

Robustness of the DSPR I and the instabilities of the DSPR II method illustrated by the corresponding short-exposure Strehl ratios in K-band for 1000 time steps of test case 1 (600 photons/subaperture/frame and zenith angle 0 deg).


From the tables, we see that, compared to the case without spiders, some small loss in quality is present when pupil segmentation due to spiders is taken into account and solved with the split approach. Moreover, one can see that the loss increases slightly for larger zenith angles and that the DSPR I provides slightly higher Strehl ratios compared to the DSPR II. Figure 8 shows the corresponding short-exposure Strehl ratios over the simulation time demonstrating that the DSPR I is indeed more stable compared to the DSPR II. More precisely, the performance of the two methods is not identical because phase averaging over the segments and matrix inversion are interchanged in the DSPR II compared to the DSPR I. As confirmed in Fig. 9 showing the residual pistons on each segment for the DSPR I, the reconstructions obtained with the split approach are free from large differential piston effects.

Fig. 9

Residual pistons in radians (K-band) on the disjoint segments for the split approach using the DSPR I. The photon flux corresponds to 100 photons/subaperture/frame using a zenith angle of 30 deg. There is almost no residual piston development. Additionally, the reconstructions do not suffer from the phase ambiguity of the pyramid sensor since the piston offsets between the segments are very much smaller than 2π radians. The global piston is subtracted from each of the segment pistons for better visibility.


Based on the provided reconstruction quality and stability features, the conclusion is that the DSPR I method is preferable. That is why we focused on the DSPR I in further numerical studies with Figs. 10 and 11 presenting acceptable METIS results obtained at a zenith angle 60 deg under very bad seeing conditions (nominal r0 at 0.5  μm=0.064  m), which corresponds to r0 at 2.2  μm=0.379  m.

Fig. 10

Simulations without telescope spiders and with telescope spiders and the split approach under very bad seeing conditions illustrated by the corresponding short-exposure Strehl ratios in K-band for 1000 time steps (600 photons/subaperture/frame and zenith angle 60 deg).


Fig. 11

Residual pistons in radians (K-band) on the disjoint segments for the split approach using the DSPR I. The photon flux corresponds to 600 photons/subaperture/frame using a zenith angle of 60 deg and very bad atmosphere. There appear slight differential pistons, which are continuously corrected during the AO loop. The global piston is subtracted from each of the segment pistons for better visibility.


However, one should keep in mind that in case the atmospheric or observational conditions change, registration and inversion of an interaction matrix will be required. This process is very time-consuming and may lead to a loss of precious operational time on the telescope. On the contrary, in the DSPR II method, the underlying interaction matrix is much smaller, therefore, significantly less operations are needed for its registration and inversion. Since the DSPR II method provides acceptable reconstruction quality as well, it can be used at least as a substitute of the DSPR I method when the latter requires some updates to be performed.

In addition, Tables 5 and 6 contain the ESO goals and ESO requirements. Note that both are defined as performance evaluations over at least 15 min of telescope operations under nominal conditions. These include some additional error sources, e.g., wind induced vibrations or noncommon path aberrations (NCPAs), which were not taken into account in the presented simulations. Therefore, we cannot straightforwardly compare our results that evaluate only 2 s of operation, with the ESO goals and requirements, but we can still see that there is a big safety gap allowed for vibrations and other error terms reducing the quality. We would like to point out that we are aware of those additional error sources and consider them in separate dedicated tests. For instance, the influence of wind shake and NCPAs on the METIS performance has been previously analyzed in Ref. 40. In this paper, we focus on eliminating the impact of spiders in the first turn. Analyzing and reducing different error terms one-by-one allows one to achieve confidence in the final performance of the instrument.

Moreover, in each of the tables, we show the roughly estimated theoretical limit of the achievable long-exposure Strehl ratio. These rough estimates are obtained using the assumption that in the high-flux case the reachable quality is limited from above by the two main error sources, the fitting, and the temporal delay error. Following the approach in Ref. 34, we evaluate the fitting error by

where d denotes the average actuator’s distance, r0 is the Fried radius, and the delay error by
where τ is the delay and τ0 is the coherence time of the atmospheric turbulence.

Furthermore, we want to remark that the results presented in the tables were obtained with the same loop gains for different guide star magnitudes and zenith angles, which clearly underline the stability of the algorithms. The DSPR II was experienced to be more sensitive to the choice of the loop gain than the DSPR I. The quality results of several numerical simulations can even be slightly improved by applying a loop gain being optimized with respect to the special parameter choices of the individual test cases.

In order to quantify the achieved raw contrast, we plot in Fig. 12 the simulated radially averaged PSF obtained with the DSPR I method. The area of interest, where the correction actually happens, is defined by the number of actuators and the telescope diameter as na(λ/D). For λ=3.7  μm, the best corrected central area of the PSF has a square shape with a side length of 1546 mas. The maximum achieved contrast is of order 106 within the control area of the DM. Figure 13 shows the difference in the simulated radially averaged point spread function (PSF) at a wavelength of λ=3.7  μm obtained with the DSPR II method compared to the DSPR I in the case of a bright guide star and zenith angle 0 deg.

Fig. 12

Simulated radially averaged PSF at a wavelength of λ=3.7  μm obtained with the DSPR I method in the case of a bright guide star and zenith angle 0 deg.


Fig. 13

Difference in the simulated radially averaged PSF at a wavelength of λ=3.7  μm obtained with the DSPR II method compared to the DSPR I in the case of a bright guide star and zenith angle 0 deg.




The first part of this paper provided a brief overview on the ability of existing algorithms to reconstruct wavefronts from PWFS data in the presence of support structures for the secondary mirror, the so-called spiders. It is known that a partial or complete shading of subapertures, depending on the thickness of the spider legs, leads to difficulties in controlling piston modes on disjoint pupil segments. We have shown that on the one hand, at least some variants of the interaction-matrix-based MVM approaches are able to handle spider obstruction effects successfully, though they are known to be computationally demanding. On the other hand, there exist several fast interaction-matrix-free model-based algorithms, which provide high-quality reconstruction in case of annular apertures without spiders. Unfortunately, they run into problems when dealing with sensor data (partially) shaded by spider legs. In end-to-end simulations, we observe in this case random uncontrolled segment pistons in the residuals. For the high-contrast large-scale AO systems in design, a combination of the advantages of both reconstructor types, i.e., the ability to reconstruct wavefronts with high precision and speed in the presence of spiders, is highly desirable.

In the second part of this paper, we presented a solution, the so-called split approach, which combines the advantages of the interaction-matrix-based and free methods and provides a high-quality wavefront reconstruction in the presence of spiders with little computational demand. The solution is based on the idea to split the reconstruction of pistons on the pupil segments from the reconstruction of higher-order modes (frequencies). For the latter part, one can use any of the available fast model-based algorithms. Especially, beneficial here is the P-CuReD, which has the smallest computational load and also provides the best reconstruction quality in case no spider legs shade the pupil.

For segment piston reconstruction, we presented two approaches. One of them requires the registration and inversion of a full interaction matrix. However, the control matrix is afterward projected to the space of segment pistons only, hence reduced in size. Another one uses a small size interaction matrix instead, which is registered in the basis of segment pistons. The resulting control matrices are in both methods of small size and the direct segment piston reconstruction step has a linear complexity.

Combined with the P-CuReD used for reconstruction of higher-order modes, the number of computations required for the complete wavefront reconstruction scales linearly with the number of controlled actuators. This represents a big advantage of the split approach compared to the usual interaction-matrix-based MVM, whose complexity scales quadratically with the number of actuators. Although MVM approaches may still be computationally doable for relatively large-scale systems like METIS having a 74×74 pyramid sensor in 2026, they are hardly feasible for extremely large systems like XAO having 200×200 subapertures and a corresponding number of actuators to be controlled in real time. The presented split approach causes, in contrast, no difficulties for the real-time implementation even on the largest AO systems of ELT-sized telescopes.

Moreover, the split approach makes existing model-based phase reconstruction algorithms developed for nonsegmented pupils suitable for wavefront control in the presence of telescope spiders. Alternatively, another idea to overcome the effects of pupil fragmentation with the model-based algorithms consists in an appropriate adjustment of the underlying forward models. A corresponding extension of the algorithms will be necessary and may allow for a stable control of segmented piston offsets without using interaction-matrix-based attempts. The question whether such an extension is possible and for which methods it is applicable needs further investigations.

The analysis of the split approach was done in the context of the instrument METIS for sensing in the K-band. Further investigations about DSPRs using pyramid sensors in different sensing wavelengths are of high interest for several instruments in development for ELTs and therefore planned future work. Although, the theory of the split approach is nonspecific to the instrument and therefore applicable to any AO system, the wavefront reconstruction may possibly suffer from phase ambiguity or other wavelength dependent effects of the pyramid sensor.


The authors are grateful to Miska Le Louarn from the European Southern Observatory (ESO) and the Austrian Adaptive Optics team for support. This work was partially funded by the Austrian Federal Ministry of Science and Research (HRSM) and the Austrian Science Fund (FWF) F68-N36, project 5.



S. Bonnefond et al., “Wavefront reconstruction with pupil fragmentation: study of a simple case,” Proc. SPIE, 9909 990972 (2016). PSISDG 0277-786X Google Scholar


R. Ragazzoni, “Pupil plane wavefront sensing with an oscillating prism,” J. Mod. Opt., 43 (2), 289 –293 (1996). JMOPEW 0950-0340 Google Scholar


J.-F. Sauvage et al., “Tackling down the low wind effect on sphere instrument,” Proc. SPIE, 9909 990916 (2016). PSISDG 0277-786X Google Scholar


M. J. Wilby et al., “A ‘fast and furious’ solution to the low-wind effect for SPHERE at the VLT,” Proc. SPIE, 9909 99096C (2016). PSISDG 0277-786X Google Scholar


K. El Hadi et al., “Pupil phase discontinuity measurement: comparison of different wavefront sensing concepts,” Proc. SPIE, 9909 990967 (2016). PSISDG 0277-786X Google Scholar


M. N’Diaye et al., “Mitigate the impact of ELT architecture on AO performance: learn from today’s telescopes to characterize and prevent the island effect,” in Proc. AO4ELT5, (2017). Google Scholar


M. N’Diaye et al., “Calibration of the island effect: experimental validation of closed-loop focal plane wavefront control on Subaru/SCExAO,” Astron. Astrophys., 610 A18 (2018). AAEJAF 0004-6361 Google Scholar


B. R. Brandl et al., “Status of the mid-infrared E-ELT imager and spectrograph METIS,” Proc. SPIE, 9908 990820 (2016). PSISDG 0277-786X Google Scholar


I. Shatokhina, “Fast wavefront reconstruction algorithms for extreme adaptive optics,” Johannes Kepler University Linz, (2014). Google Scholar


I. Shatokhina et al., “Preprocessed cumulative reconstructor with domain decomposition: a fast wavefront reconstruction method for pyramid wavefront sensor,” Appl. Opt., 52 (12), 2640 –2652 (2013). APOPAI 0003-6935 Google Scholar


I. Shatokhina, A. Obereder and R. Ramlau, “Fast algorithm for wavefront reconstruction in XAO/SCAO with pyramid wavefront sensor,” Proc. SPIE, 9148 91480P (2014). PSISDG 0277-786X Google Scholar


I. Shatokhina and R. Ramlau, “Convolution- and Fourier-transform-based reconstructors for pyramid wavefront sensor,” Appl. Opt., 56 (22), 6381 –6390 (2017). APOPAI 0003-6935 Google Scholar


S. Esposito et al., “Pyramid sensor for segmented mirror alignment,” Opt. Lett., 30 2572 –2574 (2005). OPLEDP 0146-9592 Google Scholar


S. Esposito et al., “Cophasing of segmented mirrors using the pyramid sensor,” Proc. SPIE, 5169 72 –78 (2003). PSISDG 0277-786X Google Scholar


S. Esposito et al., “Wavefront sensor design for the GMT natural guide star AO system,” Proc. SPIE, 8447 84471L (2012). PSISDG 0277-786X Google Scholar


S. Esposito and N. Devaney, “Segmented telescopes co-phasing using pyramid sensor,” in Proc. Beyond Conventional Adaptive Optics, (2001). Google Scholar


I. Surdej, “Co-phasing segmented mirrors: theory, laboratory experiments and measurements on sky,” (2011). Google Scholar


S. Esposito and A. Riccardi, “Pyramid wavefront sensor behavior in partial correction adaptive optic systems,” Astron. Astrophys., 369 L9 –L12 (2001). AAEJAF 0004-6361 Google Scholar


C. Vérinaud, “On the nature of the measurements provided by a pyramid wave-front sensor,” Opt. Commun., 233 27 –38 (2004). OPCOB8 0030-4018 Google Scholar


V. Korkiakoski et al., “Comparison between a model-based and a conventional pyramid sensor reconstructor,” Appl. Opt., 46 (24), 6176 –6184 (2007). APOPAI 0003-6935 Google Scholar


D. L. Fried, “Least-square fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am., 67 370 –375 (1977). JOSAAH 0030-3941 Google Scholar


S. Esposito et al., “Laboratory characterization and performance of the high-order adaptive optics system for the Large Binocular Telescope,” Appl. Opt., 49 G174 –G189 (2010). APOPAI 0003-6935 Google Scholar


C. Béchet, M. Tallon and E. Thiébaut, “Comparison of minimum-norm maximum likelihood and maximum a posteriori wavefront reconstructions for large adaptive optics systems,” J. Opt. Soc. Am. A, 26 497 –508 (2009). JOAOD6 0740-3232 Google Scholar


C. R. Vogel and Q. Yang, “Multigrid algorithm for least-squares wavefront reconstruction,” Appl. Opt., 45 (4), 705 –715 (2006). APOPAI 0003-6935 Google Scholar


S. Esposito, O. Feeney and A. Riccardi, “Laboratory test of a pyramid wavefront sensor,” Proc. SPIE, 4007 416 –422 (2000). PSISDG 0277-786X Google Scholar


S. Esposito, A. Riccardi and O. Feeney, “Closed-loop performance of pyramid wavefront sensor,” Proc. SPIE, 4034 184 –189 (2000). PSISDG 0277-786X Google Scholar


S. Esposito et al., “Large binocular telescope adaptive optics system: new achievements and perspectives in adaptive optics,” Proc. SPIE, 8149 814902 (2011). PSISDG 0277-786X Google Scholar


S. Esposito et al., “LBT AO on-sky results,” in Proc. of the Second AO4ELT Conf., 3 (2011). Google Scholar


E. Pinna et al., “XAO at LBT: current performances in the visible and upcoming upgrade,” in Proc. of the Fourth AO4ELT Conf., (2015). Google Scholar


B. L. Ellerbroek, “Efficient computation of minimum-variance wave-front reconstructors with sparse matrix techniques,” J. Opt. Soc. Am. A, 19 (9), 1803 –1816 (2002). JOAOD6 0740-3232 Google Scholar


R. H. Hudgin, “Wave-front reconstruction for compensated imaging,” J. Opt. Soc. Am., 67 375 –378 (1977). JOSAAH 0030-3941 Google Scholar


E. Thiébaut and M. Tallon, “Fast minimum variance wavefront reconstruction for extremely large telescopes,” J. Opt. Soc. Am. A, 27 1046 –1059 (2010). JOAOD6 0740-3232 Google Scholar


R. M. Clare and R. G. Lane, “Comparison of wavefront sensing with the Shack–Hartmann and pyramid sensors,” Proc. SPIE, 5490 1211 –1222 (2004). PSISDG 0277-786X Google Scholar


I. Montilla et al., “Performance comparison of wavefront reconstruction and control algorithms for extremely large telescopes,” J. Opt. Soc. Am. A, 27 A9 –A18 (2010). JOAOD6 0740-3232 Google Scholar


M. Le Louarn et al., “Parallel simulation tools for AO on ELTs,” Proc. SPIE, 5490 705 –712 (2004). PSISDG 0277-786X Google Scholar


M. Le Louarn et al., “Adaptive optics simulations for the European extremely large telescope,” Proc. SPIE, 6272 627234 (2006). Google Scholar


J. M. Bardsley, “Wavefront reconstruction methods for adaptive optics systems on ground-based telescopes,” SIAM J. Matrix Anal. Appl., 30 67 –83 (2008). SJMAEL 0895-4798 Google Scholar


F. Rigaut and M. van Dam, “ELT AO simulations on a laptop with YAO,” in Proc. of AO4ELT3, (2013). Google Scholar


V. Korkiakoski and C. Vérinaud, “Extreme adaptive optics simulations for EPICS,” in Proc. of the First AO4ELT Conf., 03007 (2010). Google Scholar


S. Hippler et al., “Single conjugate adaptive optics for the ELT instrument METIS,” (2017). Google Scholar


B. Engler et al., “Effects of the telescope spider on extreme adaptive optics systems with pyramid wavefront sensors,” Proc. SPIE, 10703 107035F (2018). PSISDG 0277-786X Google Scholar


V. Hutterer et al., “Wavefront reconstruction for elt-sized telescopes with pyramid wavefront sensors,” Proc. SPIE, 10703 1070344 (2018). PSISDG 0277-786X Google Scholar


I. Shatokhina, V. Hutterer and R. Ramlau, “Two novel algorithms for wavefront reconstruction from pyramid sensor data: convolution with linearized inverse filter and pyramid Fourier transform reconstructor,” in Proc. of AO4ELT5, (2017). Google Scholar


V. Hutterer and R. Ramlau, “Wavefront reconstruction from non-modulated pyramid wavefront sensor data using a singular value type expansion,” Inverse Prob., 34 (3), 035002 (2018). INPEEY 0266-5611 Google Scholar


V. Hutterer, R. Ramlau and I. Shatokhina, “Real-time adaptive optics with pyramid wavefront sensors: accurate wavefront reconstruction using iterative methods,” (2018). Google Scholar


V. Hutterer, R. Ramlau and I. Shatokhina, “Real-time adaptive optics with pyramid wavefront sensors: a theoretical analysis of the pyramid sensor model,” (2018). Google Scholar


V. Hutterer, “Model-based wavefront reconstruction approaches for pyramid wavefront sensors in adaptive optics,” (2018). Google Scholar


V. Hutterer and R. Ramlau, “Non-linear wavefront reconstruction methods for pyramid sensors using Landweber and Landweber–Kaczmarz iteration,” Appl. Opt., 57 (30), 8790 –8804 (2018). APOPAI 0003-6935 Google Scholar


M. Rosensteiner, “Cumulative reconstructor: fast wavefront reconstruction algorithm for extremely large telescopes,” J. Opt. Soc. Am. A, 28 2132 –2138 (2011). JOAOD6 0740-3232 Google Scholar


M. Rosensteiner, “Wavefront reconstruction for extremely large telescopes via CuRe with domain decomposition,” J. Opt. Soc. Am. A, 29 2328 –2336 (2012). JOAOD6 0740-3232 Google Scholar


M. Zhariy et al., “Cumulative wavefront reconstructor for the Shack-Hartman sensor,” Inverse Prob. Imaging, 5 893 –913 (2011). Google Scholar


R. Clare and M. L. Louarn, “Numerical simulations of an Extreme AO system for an ELT,” in Proc. AO4ELT2, 1100 –1107 (2011). Google Scholar


R. M. Clare et al., “Numerical evaluation of pyramid type sensors for extreme adaptive optics for the European extremely large telescope,” in Proc. of the Fifth AO4ELT Conf., (2017). Google Scholar


U. Bitenc et al., “On-sky tests of the CuReD and HWR fast wavefront reconstruction algorithms with CANARY,” Mon. Not. R. Astron. Soc., 448 (2), 1199 –1205 (2015). MNRAA4 0035-8711 Google Scholar


C. Z. Bond et al., “Iterative wave-front reconstruction in the fourier domain,” Opt. Express, 25 11452 –11465 (2017). OPEXFF 1094-4087 Google Scholar


A. Obereder et al., “Dealing with spiders on ELTs: using a pyramid wfs to overcome residual piston effects,” Proc. SPIE, 10703 107031D (2018). PSISDG 0277-786X Google Scholar


A. Burvall et al., “Linearity of the pyramid wavefront sensor,” Opt. Express, 14 (25), 11925 –11934 (2006). OPEXFF 1094-4087 Google Scholar


E. Pinna et al., “Phase ambiguity solution with the pyramid phasing sensor,” Proc. SPIE, 6267 62672Y (2006). PSISDG 0277-786X Google Scholar


Victoria Hutterer received her BS and MS degrees in mathematics from Johannes Kepler University Linz, Austria, in 2014 and 2015, respectively. She is a PhD student at the Johannes Kepler University Linz. Her current research interests include inverse problems, adaptive optics, pyramid wavefront sensors, wavefront reconstruction algorithms, and iterative and nonlinear algorithms.

Iuliia Shatokhina received her BS and MS degrees in physics from the National University of Kyiv-Mohyla Academy in 2006 and 2007, respectively, and her PhD in mathematics from Johannes Kepler University Linz, Austria, in 2014. She is a research assistant at Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences. Her current research interests include inverse problems, adaptive optics, pyramid wavefront sensors, and wavefront reconstruction algorithms.

Andreas Obereder received his MS degree in mathematics in 1998 and his PhD in mathematics from Johannes Kepler University Linz, Austria, in 2003. He is a research assistant at Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences. His current research interests include inverse problems, astronomical adaptive optics, and wavefront reconstruction algorithms.

Ronny Ramlau received his MS degree in mathematics in 1994 and his PhD in mathematics from the University of Potsdam, Germany, in 1997. He is a university professor at Johannes Kepler University Linz, Austria, where he leads the Industrial Mathematics Institute. Also, he is a scientific director at Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austrian Academy of Sciences, leading the Transfer Group. He is the author of more than 50 journal papers and has written 3 book chapters. His current research interests include nonlinear inverse problems, medical imaging, regularization methods, and iterative methods.

© The Authors. Published by SPIE under a Creative Commons Attribution 4.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Victoria Hutterer, Iuliia Shatokhina, Andreas Obereder, and Ronny Ramlau "Advanced wavefront reconstruction methods for segmented Extremely Large Telescope pupils using pyramid sensors," Journal of Astronomical Telescopes, Instruments, and Systems 4(4), 049005 (19 November 2018).
Received: 2 July 2018; Accepted: 25 October 2018; Published: 19 November 2018

Back to Top