New algorithms for reconstructing wavefront from slopes data are developed, which exhibit high accuracy over broad spatial-frequency bandwidth. Analyzing wavefront reconstructors in the frequency domain lends new insight into ways to improve frequency response and to understand noise propagation. The mathematical tools required to analyze the frequency domain are first developed for discrete band-limited signals. These tools are shown to improve frequency response in either spatial-or frequency-domain reconstruction algorithms. A new spatial-domain iterative reconstruction algorithm based on the Simpson rule is presented. The local phase estimate is averaged over 8 neighboring points whereas the traditional reconstructors use 4 points. Analytic results and numerical simulations show that the Simpson-rule–based reconstructor provides high accuracy up to 85% of the bandwidth. The previously developed rectangular-geometry band-limited algorithm in frequency domain is adapted to hexagonal geometry, which adds flexibility when applying frequency-domain algorithms. Finally, a generalized analytic expression for error propagation coefficient is found for different reconstructors and compared with numerical simulations.
© 2011 OSA
Frequency-domain wavefront reconstruction methods are as old as the very early wavefront reconstructors [1, 2]. Freischlad placed this subject on solid ground . The rectangular map constraint of the conventional Fourier method has been removed in an iterative Gerchberg-type algorithm dealing with an arbitrary boundary shape . A series of recent papers by Poyneer discusses improvements on handling boundary conditions and applications in extreme adaptive optics [5, 6]. Similar principles have been applied in shearing interferometers . More serious attention has been paid to the accuracy of the reconstruction methods in Refs. [8–10]. The works of Campos and Yaroslavsky presented a one-dimensional solution based on band-limited integration technique in frequency domain. The two-dimensional (2-D) extension of the same method was not discussed. Complementary to their works, Bahk has introduced a full two-dimensional wavefront reconstructor based on the band-limited derivative calculation . Both approaches emphasize the frequency response of the reconstructed signals. The frequency response of wavefront reconstruction has been discussed earlier in the analysis of lateral shearing interferometry . Frequency response characteristic of a reconstructor is important in focal spot diagnostic for high power lasers where the focal spot is indirectly characterized using wavefront information reconstructed from Shack-Hartmann slopes data . The error in the high-frequency components of the reconstructed wavefront results in an error in the estimate of encircled energy.
This paper develops a set of encompassing mathematical tools for wavefront reconstruction problems, where many additional benefits naturally arise interconnecting the results of previous works. The benefits are exemplified by the development of two new wavefront reconstructors which exhibit high accuracy over broad spatial-frequency bandwidth. It also provides a unified analytic expression for the noise propagation coefficients of several well-known traditional wavefront reconstructors as well as the ones developed here.
This paper consists of four parts. Section 2 introduces the mathematical tools and symbols regarding band-limited derivative operations, which are needed for the analyses in the subsequent sections. In Sec. 3, a way to improve the accuracy of the finite difference method is discussed in connection with wavefront reconstruction. The Simpson rule is adopted for developing a new spatial-domain iterative reconstruction algorithm. The exact details of the algorithm and its frequency domain property are described. In Sec. 4, the band-limited reconstruction algorithm developed for rectangular geometry  is extended to hexagonal geometry, which greatly enhances the flexibility of this scheme. Finally, in Sec. 5 the noise propagation curve is analytically derived and compared with numerical simulations.
2. Band-limited derivative
The main results of band-limited derivative techniques in the context of wavefront reconstruction were summarized in . The full derivation of the results will be presented here for the sake of completeness. Additional new notations are introduced which simplify the expressions in Sec. 4.
The motivation for band-limited derivatives, especially for discrete samples, lies in the fact that it provides analytical tool for converting back and forth between slopes measurements and wavefront signal. We start by asking what the exact interpolation formula is for derivatives in discrete samples. According to sampling theorem , a band-limited signal can be exactly reconstructed at any point by convolving a sinc-function with discrete samples.
Using the result
Using the fact that j 1 is an odd function, the derivative interpolation expression at discrete points is
The summation of the left hand side of Eq. (4) for all sample points can be shown to be equal to zero by taking advantage of the right hand side expression and using the periodicity condition of discrete samples:
The coefficient of (k) on the right hand side is a Fourier series which represents a sawtooth wave. Thus Eq. (7) is further simplified as
Equation (8) provides a convenient way of calculating exact derivatives from band-limited signals.
For the sampling points of a derivative signal offset by a half sampling space from the sampling points of the original signal, a slightly different form is derived going through similar steps as in Eqs. (1)–(4),
N can be either an odd or even integer as the floor notation in the above equation will apply to either type.
We also need an interpolation formula for creating a signal shifted by half sample spacing for the Fried geometry,
The DFT of Eq. (14) is
Thus the partial derivative in x-direction for Fried geometry in frequency domain is
Equation (17) has an additional degree of freedom (index p for y-direction) because the reconstructed sample point in the Fried geometry has to be first shifted in y-direction by half-sample size before applying the half-sample shifted derivative operation in x-direction. For consistency, we can verify that the sequential operations of half-pixel shift and derivative operation using 𝖱(k) and 𝖲(k) produce the same result as single operation of 𝖳(k), that is, 𝖳(k) = 𝖱(k) 𝖲(k). This relation, however, does not hold for the value at N/2 for even N, where the left hand side is 0.5 whereas the right hand side is 0. To remove this paradox for even number of samples, we choose to use 𝖲(N/2) = 0.5 and 𝖱(N/2) = 1 or 𝖲(N/2) = −0.5 and 𝖱(N/2) = −1. A similar choice was made in  in defining band-limited integration operators to prevent information loss at the boundary. The redefinition of 𝖲 and 𝖱 at the mid-point is implied hereafter. Using the new definition, the half-pixel operator used in the right hand side of Eq. (15) can be alternatively expressed as
𝖲 can be considered as a discrete angular frequency vector circularly shifted by ⌊N/2⌋ or ⌊N/2⌋ – 1 using the new definition at the mid-point. For an odd number of samples,
For an even number of samples with 𝖲(N/2) = −0.5,
We can collectively define the right hand side of Eqs. (19)-(21) as where the bar over kx implies a circular shift of the discrete angular frequency vector centered around zero by ⌊N/2⌋ or ⌊N/2⌋ – 1 as prescribed above. Using the new notation and Eq. (18) as well as the relation 𝖳(k) = 𝖱(k)𝖲(k), Eqs (8), (12), and (15) can be alternatively expressed as
The notation allows one to see the formal similarity with continuous variable derivative expressions. The continuous variable counterpart for Eq. (22), for example, is (where FT denotes Fourier transform and kx here is a continuous spatial frequency variable. This connection allows one to derive other derivative formulas for discrete band-limited samples by taking advantage of continuous space results.
In many practical situations the band-limited calculations may not produce exact results depending on the nature of signals. The magnitudes of the Fourier coefficients of a linear function, for example, decrease as 1/(spatial frequency) whereas Eq. (22) indicates that the coefficients of the derivative is multiplied by spatial frequency vector. Therefore the highest spatial frequency coefficient does not vanish even if N approaches ∞. The linear terms thus are not band-limited and need to be treated separately using a polynomial fit.
3. Simpson reconstructor
The analysis in Sec. 2 suggests that the accurate derivative calculation at discrete samples requires the superposition sum of the whole available samples. This can be done more conveniently in the frequency domain, which results in a band-limited reconstructor with unity frequency response . On the other hand, it is still worth while to investigate an improved finite difference scheme for purely spatial domain operation. In finite difference methods involving only a few points, a high degree of accuracy can be obtained by distributing the finite difference over both the measured derivative samples (i.e. slopes) and the integrated samples. Denoting the wavefront estimate as and the measured slope as S, one can start from a general finite difference expression such as15] showed a reconstructor based on
The frequency response of Southwell reconstructor is low at high spatial frequency . We enhance the frequency response using the Simpson rule which is
An iterative wavefront reconstructor based on this scheme will be developed in the next section. Other finite schemes which involve more than three samples (slopes) or have asymmetric forms were not considered at this point as the iteration scheme might become unnecessarily complicated.
3.1. Simpson iterator
Casting the local 1-D equation (27) into a least squares form in 2-D, we obtain an error metric (ɛ) as
Δx and Δy are moved around to the side so that squared terms are in units of slopes. This provides equal weight to the differences in x and y directions on the assumption that the magnitude of slopes is comparable in either direction.
The condition ∂ɛ/∂(i, j) = 0 leads to an equation that can be used for an iterative algorithm. At the stationary point, the error is not decreasing anymore which indicates that the phase estimate has reached a solution. It is assumed that the boundary of phase and slope points can be an arbitrary shape. The differentiation of the error metric results in four groups, which are indicated by different colors in Fig. 1. Each group can be used in the equation only if all its elements exist. This strategy is realized by using g-parameters as shown in the following:
gL, gR, gU, gD are flags with values 0 or 1 where L, R, U, D indicate left, right, up, and down directions, respectively. They are zeros if the slopes quantity in the parenthesis next to them are incalculable or ones otherwise. gL at the point (i, j), for example, does not vanish only if all the slopes measurements exist at the additional points (i, j – 1) and (i, j – 2). The scope of each flag is graphically indicated in Fig. 1. For a comparison, the iterative equation for the Southwell reconstructor is written here using the same format.
The successive-over-relaxation technique in  can be applied to the Simpson iterative reconstructor:
3.2. Frequency response and regularization
The frequency response of the Simpson reconstructor will be calculated using a frequency-domain least-squares method. The sum of the squared error in spatial domain in Eq. (28) is equivalent to the sum of the squared error of the Fourier transformed component by the Parseval theorem:
The solution for in Eq. (34) is
As and , the frequency response H defined as the ratio of the reconstructed wavefront amplitude to the true wavefront amplitude associated with the measured slopes at a given spatial frequency point isFig. 1 of ) over all spectrum.
In frequency domain, this transforms into
The denominator of the Simpson frequency response will have an additional term of λ (|Dx, reg|2 + |Dy, reg|2) that removes the singularity. The regularized frequency response is
The 1-D frequency response (ωy = 0) with the regularization term has a second maximum near high spatial frequency for sufficiently small λ (<0.08) [Fig. 2(b)]. The free parameter λ can be fixed to a value such that the second maximum is equal to one. Numerically determined value of λ for such condition is 0.07489. This choice of λ gives only 3% error in wavefront amplitude over 80% of the frequency range. Another choice can be λ = 0.07026 which balances the local maximum and minimum around 1. The second option reduces the maximum deviation below 2.2% within 85% of the spectral range. Figure 2(a) shows three-dimensional view of the frequency response of the Simpson-rule reconstructor with λ = 0.07489. The 1-D response is shown in Fig. 2(b). The solid line is calculated from analytic expression [Eq. (43)] whereas the circles are from numerical simulations. The numerical simulation consists of steps of generating slopes from sinusoid wavefronts at a given spatial frequency using analytic derivatives and of reconstructing the wavefront and comparing the ratio between the original and the reconstructed wavefront amplitude at that frequency. The reconstruction algorithm used in the simulation will be explained in detail in the following section. The result shows good agreement with the analytic curve.
3.3. Iterative algorithm with regularization terms
The frequency-domain analysis does not provide a detailed picture how the successive over-relaxation method can be applied in spatial domain iteration especially around the boundary of slopes map. Resolving the stationary condition with the regularization term introduces additional terms on the left hand side of Eq. (29). These are fully written out using g flags.
gLR or gUD is one only if there exist both points to the right and left or up and down of the point (i, j), respectively and zero otherwise. The iteration formula (32) will be modified to
4. Hexagonal band-limited reconstructor
Iterative Fourier-domain band-limited reconstruction method as proposed in  provides unity frequency response over all spatial bandwidth as long as the signals are band-limited. The band-limited reconstructor was designed for the Southwell geometry; the reconstruction point and the corresponding measured slopes are all at the same point on a rectangular grid. It was shown that band-limited derivative operators are also available for Hudgin and Fried geometries . Table 1 summarizes the three operators belonging to each geometry.
Here we present band-limited reconstructors for hexagonal arrays. Hexagonal geometry may be well suited for adaptive optic systems for large telescopes with hexagonal mirror arrays (e.g., James Webb) . Large deformable mirrors used in laser fusion facilities (National Ignition Facility)  also have hexagonal actuator patterns. The number density of lenslets is slightly higher in hexagonal geometry than square. Figure 3 shows two possible hexagonal arrays. In Fig. 3(a) the unit hexagon is lying on its facet (prostrate hexagon), whereas in Fig. 3(b) the unit hexagon is standing on the apex (standing hexagon). The circles indicate the measured slopes points and the ×’s are reconstruction points. The reconstruction points are set on a rectangular grid because rectangular arrays are easier to visualize in most computer programs.
As will become clear in the discussion of the algorithm, we need procedures for performing DFTs on a hexagonal array. This can be achieved by grouping the slopes into two interleaving groups, which are indicated by red (group 1) and black (group 2) circles in Fig. 3. Starting from Fig. 3(a) geometry, it can be shown that the DFT of the slopes on the rectangular array (on × points) can be expressed as the linear sum of smaller DFTs of the sub-group 1 and 2 (see Appendix).
T hex is a matrix whose size is M by N (i.e. the size of either array 1 or array 2) and “○” denotes entrywise matrix multiplication. The pth row and qth column element of T hex is
The combined total DFT array is therefore a vertical concatenation of the two matrices. On the other hand, the resulting total DFT matrix for Fig. 3(b) geometry is a horizontal concatenation,
The same combination rule applies to y-slope measurements.
The above decomposition technique can be inverted such that each subgroup of the hexagonal array can also be expressed as the linear sum of the block I and II of the rectangular array. This inversion is only used for wavefront points in the algorithm, which is
Using the basic results obtained in Sec. 2 and the DFT procedures for the hexagonal arrays in this section, the band-limited reconstruction algorithm for hexagonal slopes arrays can be implemented as shown in the flowchart in Fig. 4.
Step 1 consists of fitting the slopes over low order polynomials, e.g. third order, which will significantly reduce non-band-limited components of the wavefront. If the regions of interest are disconnected, the fitting has to be performed per each region. Owing to the sum requirement [Eq. (6)], a column and a row are appended to the edge of the measured slope matrices (group 1 and 2 separately), which will satisfy the zero sum conditions in x and y directions.
Step 2 initializes the slopes with measured values. Steps 3-8 form a closed loop required for extrapolating slopes outside the non-rectangular region. The iteration is not required if the region is rectangular.
Slopes in group 1 and 2 are separately Fourier-transformed using Eqs. (47)–(50) in Step 3. In Step 4, wavefront matrices corresponding to each block (I or II) are reconstructed in Fourier domain using the band-limited filter function which is
Step 5 creates wavefront group 1 and 2 by using Eqs. (51)–(52). In Step 3 and 5, the correct T hex must be used according to its geometry. In Step 6, derivative operators are applied to these temporary wavefront matrices to obtain slopes in group 1 and 2, respectively. These new slopes are different from the measured slopes. We leave the values external to the boundary untouched while restoring the internal values to the original measured slopes. The difference between the measured slope and the calculated slopes decreases over the course of iterations. Step 8 determines whether this difference is within tolerance. Once the convergence criterion is met, the wavefront matrices generated in Step 4 ( , ) are combined to form a single matrix either by vertical or horizonal concatenation depending on the hexagon geometry and inverse Fourier-transformed to the spatial domain to produce the final result in Step 9. Small terms in the imagingary part of the solution can be neglected.
The band-limited algorithms as shown in Fig. 2 of  and Fig. 4 can be used together with a non-band-limited filter function, which enables a convenient way to switch between different algorithms. The reconstruction algorithms proposed here are not limited to a specific boundary shape.
5. Error propagation
The wavefront reconstructors have traditionally been characterized with so called error propagation curve. This indicates the sensitivity of the noise in the reconstructed phase to the noise in the slopes measurements. Early numerical and theoretical works show that this sensitivity is a logarithmic function of the number of measurement points [1, 2, 15]. The noise propagation coefficient will be calculated using discrete samples and frequency domain filter functions. We limit the scope to the rectangular area.
Let σφ be the root-mean-square of the reconstructed phase φ. According to the Wiener-Khintchin theorem,
According to linear stochastic system theory, the power spectrum of input and output signal is related by the absolute square of the linear system function. In the case of wavefront reconstruction dictated by the following linear response
Assuming that and are uncorrelated white noise with variance of for each, and since , the noise propagation coefficient η is19] in the case of band-limited operators.
The right hand side of Eq. (61) is inversely proportional to |D 0 ,x|2 + |D 0 ,y|2 for band-limited reconstruction and difficult to visualize in linear scale. We define a ‘noise response function’ SN with the inverse power dependence removed as follows:
It can be shown by the Cauchy-Schwartz inequality that the noise response is always larger than or equal to absolute squared of the frequency response,
The inequality (63) shows that the error propagation is intimately related with the frequency response of a reconstructor. The lower bound of η is
From this one can expect that the Southwell reconstructor will have the lowest lower bound and the Fried reconstructor the highest. It agrees with the result of Zou .
The analytic expression for η can be calculated and fit to a logarithmic curve, although the logarithm dependence is only approximate except for the band-limited reconstructors. The result is summarized in the second column of Table 3. Singularity points were excluded in the summation over spatial-frequency space. The third column shows the simulated η obtained by running actual reconstructors with zero slopes input with Gaussian noise. The Fried reconstructor was not tested. Two hundred realizations were performed at each N, where N 2 is the number of points. N was varied from 10 to 100 by 10. The logarithm fit over the averaged η is shown in the column. The multiplicative coefficients roughly agree with the analytic ones up to the second decimal point, but the additive constants from simulation are always estimated higher than the calculated ones. The offset is about 0.2771 on average. The discrepancy appears to come from the apparent inconsistency in assuming white noise in slopes power spectrum and the use of band-limited derivative formalism. For example, the reconstructed wavefront from white noise spectrum always has some amount of low-order polynomial terms, which cannot be represented by Eq. (59). The constant offset 0.2771 therefore can be considered as the ratio of energy conversion from white noise to non-band-limited signals.
The legacy formulas of noise propagation for each reconstructor are also shown in the fourth column of Table 3, quoted from the three authors’ original publications [1, 2, 15]. The quoted Southwell η is estimated only from the graph in the original paper since no explicit formula was given. Noll’s calculation essentially corresponds to band-limited case. Considering the fact that there is some ambiguity in the determination of the constant offset, at least the multiplicative coefficient of the Fried formula comes close to our analytic result; whereas there is about a factor-of-2 difference in the Southwell and Noll’s expression compared with ours. On the other hand, Hudgin’s formula does not agree with our results. Fried’s formulas is based on comparatively large number of N (≤ 39) compared with Southwell and Hudgin’s calculations (N ≤ 20).
We have presented detailed derivations of band-limited derivative operators in frequency domain. These are important tools for characterizing and improving frequency response of wavefront reconstructors over broad bandwidth. Two new wavefront reconstructors were proposed utilizing these tools. The reconstructors were designed to be accurate up to high spatial frequency. The first one is based on the Simpson integration rule. The bandwidth of the frequency response of this reconstructor, after being regularized, is excellent up to 85% of the spatial frequency range. A successive-over-relaxation iterative solver was presented in detail where the outermost samples are elegantly handled using g flags. The frequency response behavior of the iterative solver agrees well with the predicted frequency response curve. The second reconstructor is an extension of the band-limited reconstruction algorithm previously developed; the measurement points are on a hexagonal array instead of a rectangular array. A Fourier-domain iterative algorithm was proposed for two types of hexagonal arrays. As was pointed before in , the reconstruction process must be preconditioned with low-order polynomial fit. The Simpson-rule based algorithm works purely in spatial domain; therefore it is computationally less complex than band-limited algorithms, whereas the latter provides flexibility against any geometry change. Fourier-domain algorithms have a potential of boosting reduction speed with the help of digital signal processors.
The new wavefront reconstructors are compared with the traditional reconstructors in terms of noise propagation properties through a generalized noise propagation expression. The analytically calculated noise propagation coefficients are consistent with the numerical fit deduced from our own simulations. However we did not find universal agreement with the published results.
The broad-bandwidth wavefront reconstructors developed here are used in wavefront reduction software for characterizing focal spot of OMEGA EP laser beams . The importance of band-limited reconstructor was well illustrated in  for a closed-loop wavefront shaping application. One may also find applications in the study of metrology and atmospheric turbulence .
Equation (47) will be derived for the prostrate hexagon geometry. Supposing the x-slopes are measured on the M by N rectangular grid instead of the hexagonal grid and denoting them with primes, the DFT of S′xEqs. (47), (48) is same as M′ here. The Fourier domain indices are p and q which range from 0 to M – 1 and from 0 to N – 1, respectively.
Denoting the first half range of index p as p′ (p′ = 0, ...,M′ − 1), the above equation becomes
Slopes on the hexagonal grid are and using half-sample shift operator 𝖲. The shift operator and the extra phase factor constitute T hex as defined in the main text [Eq. (48)]. Therefore Eq. (67), (68) can be expressed as a block matrix form as in Eq. (47). The first and second block will be denoted as and , respectively. The DFT of the group 1 or 2 slopes is to be interpreted as operating on matrices of size M′ by N.
The relation between the slopes group 1 and 2 on the hexagonal grid and the first and second blocks on the rectangular grid can be inverted. In other words, the slope group 1 or 2 can be expressed as the linear sum of the first and the second blocks in the Fourier domain.
The block-decomposition technique and the inversion formula can be applied to y-slopes and wavefront points. The results for the standing hexagon geometry can be similarly derived.
This work was supported by the U.S. Department of Energy Office of Inertial Confinement Fusion under Cooperative Agreement No. DE-FC52-08NA28302, the University of Rochester, and the New York State Energy Research and Development Authority. The support of DOE does not constitute an endorsement by DOE of the views expressed in this article.
References and links
1. D. L. Fried, “Least-squares fitting a wavefront distortion estimate to an array of phasedifference measurements,” J. Opt. Soc. Am. 67, 370–375 (1977). [CrossRef]
2. R. H. Hudgin, “Wave-front reconstruction for compensated imaging,” J. Opt. Soc. Am. 67, 375–378 (1977). [CrossRef]
3. K. Freischlad and C. L. Koliopoulos, “Modal estimation of a wavefront from difference measurements using the discrete Fourier transform,” J. Opt. Soc. Am. A 3, 1852–1861 (1986). [CrossRef]
5. L. A. Poyneer, D. T. Gavel, and J. M. Brase, “Fast wave-front reconstruction in large adaptive optics systems with use of the Fourier transform,” J. Opt. Soc. Am. A 19, 2100–2111 (2002). [CrossRef]
6. L. A. Poyneer and B. Macintosh, “Spatially filtered wave-front sensor for high-order adaptive optics,” J. Opt. Soc. Am. A 21, 810–819 (2004). [CrossRef]
8. C. Elster and I. Weingärtner, “High-accuracy reconstruction of a function f(x) when only or is known at discrete measurement points,” in X-Ray Mirrors, Crystals, and Multilayers II, A. K. Freund, A. T. Macrander, T. Ishikawa, and J. L. Wood, eds., Proc. SPIE 4782, 12–20 (2002).
9. J. Campos, L. P. Yaroslavsky, A. Moreno, and M. J. Yzuel, “Integration in the Fourier domain for restoration of a function from its slope: comparison of four methods,” Opt. Lett. 27, 1986–1988 (2002). [CrossRef]
13. J. Bromage, S.-W. Bahk, D. Irwin, J. Kwiatkowski, A. Pruyne, M. Millecchia, M. Moore, and J. D. Zuegel, “A focal-spot diagnostic for on-shot characterization of high-energy petawatt lasers,” Opt. Express 16, 16561–16572 (2008). [PubMed]
14. C. E. Shannon, “Communication in the presence of noise,” Proc. Inst. Radio Eng. 37, 10–21 (1949).
15. W. H. Southwell, “Wavefront estimation from wavefront slope measurements,” J. Opt. Soc. Am. 70, 998–1009 (1980). [CrossRef]
16. B. L. Phillips, “A technique for the numerical solution of certain integral equations of the first kind,” J. ACM 9, 84–97 (1962). [CrossRef]
17. See http://www.jwst.nasa.gov/.
18. S. E. Winters, J. H. Chung, and S. A. Velinsky, “Modeling and control of a deformable mirror,” J. Dyn. Sys. Meas. Control 124, 297–302 (2002). [CrossRef]
19. R. J. Noll, “Phase estimates from slope-type wave-front sensors,” J. Opt. Soc. Am. 68, 139–140 (1978). [CrossRef]
20. W. Zou and J. P. Rolland, “Quantifications of error propagation in slope-based wavefront estimations,” J. Opt. Soc. Am. A 23, 2629–2638 (2006). [CrossRef]
22. T. W. Nicholls, G. D. Boreman, and J. C. Dainty, “Use of a Shack-Hartmann wave-front sensor to measure deviations from a Kolmogorov phase spectrum,” Opt. Lett. 20, 2460–1462 (1995). [CrossRef] [PubMed]