Particle swarm optimization for wavefront correction in ophthalmic applications

Many optical systems require the correction of the cumulative wavefront error of the system for performance optimization. In ophthalmology the wavefront error of the eye corresponds to the visual defect and can be measured up to high-order aberrations today. Lower orders of the wavefront error are usually corrected with spectacles, contact lenses or refractive surgery. In this paper we apply an optimization method called particle swarm optimization to calculate the optimal correction for visual defects based on measured high-order wavefront results. It is shown that an optimized conventional correction and in particular an extended wavefront correction including higher orders yields a better result for the visual Strehl ratio, as compared to a simple conjugate correction scheme.


Introduction
Visual disorders like myopia (nearsightedness) or hypermetropia (farsightedness) are conventional corrected with spectacles, contact lenses or refractive surgery. The optimal refraction of the individual needed correction is determined by the ophthalmologist. The gold standard for this determination is the subjective refraction, where successive different spherical shaped lenses are applied to achieve the best visual acuity. In such an examination the optical aberrations of the eye are not measured, but only the best subjective correction is determined.
Nowadays, the optical aberrations of the eye can be measured fast and precisely with aberrometers. An overview of the available technologies used in these systems is described by Charman in [1]. Regardless of the used technology, the measuring result is always a wavefront error, containing low-and high-order aberrations. The optimal correction should be calculated by a metric for the resulting subjective visual perception. Such a metric should consider all relevant optical and neural factors like image formation and processing for human vision. Studies have shown that the visual Strehl ratio computed by the OTF (VSOTF) describes the visual performance most precisely [2][3][4].
A first attempt to calculate the best correction of the wavefront error was made by Guirao et al [5]. They solved the three-dimensional objective function with an exhaustive method. Unfortunately, it is not described what kind of method they used in detail. Another approach to calculate the optimal correction was introduced by Leube et al in [6]. They trained a deep learning machine with a database of wavefront errors and the respective subjective correction.
Studies have shown that for healthy eyes the wavefront error mainly contains low-order aberrations (LOA) [7], which can be completely corrected with the conventional method. However, the wavefront error of highly irregular corneas, e.g. keratoconus also contains high-order aberrations (HOA) [8], which impairs the visual performance. For an improved visual performance these high-order aberrations need to be included in an optimized calculation procedure.
In this paper, we use a particle swarm optimization algorithm to solve the complex and non-linear optimization problem. The algorithm is flexible to solve corrections with different degrees of freedom. Therefore, a conventional correction of low-order aberration and a wavefront correction including low-and high-order aberrations can be calculated. The wavefront correction, we introduce in this paper, is an extension of the conventional correction with trefoil and coma aberrations. Due to the low-order aberrations cos θ X coma 9 √ 8ρ 3 cos 3θ X trefoil of the conventional correction only quadratic aberrations can be corrected. By including trefoil and coma cubic aberrations can be corrected additionally.
In the first part of this paper we illustrate an algorithm for solving the complex and non-linear objective function. The effects of the conventional, wavefront and simplified conjugate correction on the visual performance are examined in the second part.

Wave-aberrations
In wave optics aberrations are deviations of a diffraction limited system in the specific pupil [9]. These deviations are the optical aberrations of the system and are described as a wavefront error. In ophthalmology there are usually listed in microns. Mathematically, the wavefront is related to the phase of a harmonic wave with the wavenumber k = 2π λ , where λ is the wavelength of the light. In ophthalmology the wavefront error of the human eye is measured with an aberrometer [1]. The result of these systems is the wavefront error in the circular exit pupil of the eye. This means that for an eye without any visual defects the wavefront error is equal to 0. The radius of this pupil is described by r 0 .
Analytically, the wavefront W(r, ϕ) is described with the normalized coordinate ρ = r/r 0 as Zernike polynomials and can be written as where Z i (ρ, ϕ) is the ith Zernike polynomial with the specific coefficient c i and n as the maximum number of Zernike polynomials. In ophthalmology, the wavefront is usually specified with 36 Zernike coefficients (n = 35) [10]. The calculation of the Zernike polynomials Z i (ρ, ϕ) is described elsewhere and can be found in [11]. As there are many notations of Zernike polynomials and coefficients, in this paper we use the notation of the Optical Society of America (OSA) [11]. The first 10 Zernike polynomials with the respective description are listed in table 1. The first 6 polynomials are known as low-order aberrations, whereas all other polynomials are declared as high-order aberrations. The Zernike coefficients c i can be interpreted as weighting factor for the particular Zernike polynomial and such for the optical aberration. The sum of equation (1) can also be rewritten as inner product of a row vector ⃗ Z and column vector⃗ c as To correct such a wavefront error, the phase aberrations needs to be compensated so that the resulting wavefront of the eye and the correction is equal to 0. If the distance between the pupils of the correction and the eye is approximately 0, the wave propagation of the light can be neglected. Therefore the resulting wavefront W Res can be calculated as a sum of the correcting wavefront W Corr and the wavefront error of the eye W Eye . Thus, the resulting wavefront can be written with⃗ c as the Zernike coefficients measured from the eye and ⃗ x as variable coefficients of the correction as The resulting wavefront W Res of a wavefront error depends only on the Zernike coefficients ⃗ x of the correction. Equation (4) implies that the optimal correcting wavefront is the conjugate (negative) wavefront of the eye. Therefore, the resulting wavefront is 0 and all aberrations are compensated. However, for several reasons it is not sufficient to use the wavefront directly. One reason is, that the measured wavefront of the eye is just a snapshot of the current error and especially the high-order aberrations varies slightly over multiple measurements [12].
Moreover, note that equation (4) is only valid if the pupil radii of the wavefront error and correction are identical as the pupil diameter of the human eye varies with illumination. Otherwise, the addition of the coefficient vectors⃗ c and ⃗ x is invalid.

Pupil function
The optical aberrations are encoded in the phase and can be expressed with other optical properties like transmission and reflection with the complex pupil function A(r, ϕ) [9].
In 1933, Stiles and Crawford reported in [13] that light rays from the periphery of the pupil are perceived weaker compared to rays from the center with the same intensity. This first Stiles-Crawford-Effect is integrated in the pupil function as amplitude a(r) [14] with On one hand the pupil function contains the aberrations of the optical system represented as the phase and on the other hand the perception sensitivity of the Stiles-Crawford-Effect is included as the amplitude. The pupil function is the base formula for all further calculations of optical properties.

Point spread function and optical transfer function
The monochromatic point spread function (PSF) is the impulse response of an optical system to a light point source of a specific wavelength. In an ideal optical system, where no aberrations exist, the PSF has a bright central region surrounded with concentric rings, well known as Airy disk. If optical aberrations are existing in a system, this Airy disk will be distorted in the image plane. In Fourier optics the PSF is the image of the aperture of the Fraunhofer diffraction and can be calculated by the pupil function [9]. The PSF h is the intensity of the Fourier transform of the pupil function.
In ophthalmology the units of the PSF ψ x , ψ y are usually given in arc minutes [15]. The PSF depends only on the Zernike coefficients of the correction ⃗ x like the resulting wavefront and the pupil function. The PSFs of a wavefront error with three different corrections are shown in figure 1. The left figure 1(a) shows the pure PSF of the wavefront error without any correction. In the middle 1(b) the wavefront error is corrected with low-order aberrations and the right figure 1(c) shows the diffraction limited case, where all aberrations of the eye are compensated.
As mentioned in the beginning the PSF is calculated for one wavelength. However, the visible spectrum is in the range of 400 and 700 nm, whereas the highest sensitivity of the human retina is at 555 nm [16]. For a complete calculation of an optimal correction the full spectrum of the visible light have to be considered. The calculation of the polychromatic PSF can be found in [14]. For the comparison of different correction types, however, the use of monochromatic PSFs is sufficient.
Since the PSF is the impulse response of a point light source, a transferred image through an optical system can be simulated as convolution of the object and the PSF [14]. The objects imaged through an optical system with different PSF from figure 1 are shown in figure 2. In case of an aberrated eye, these images are transferred on the retina for different correction types.
The transmission quality for spatial frequencies is given in the complex optical transfer function (OTF). The function describes the transfer quality through an optical system of specific spatial frequencies in units of cycles per degree (cpd). Mathematically, the OTF H is defined as the Fourier transform of the PSF h [9].
The normalized radial averaged OTFs of the PSFs 1 are illustrated in figure 3. These OTFs reflect the convoluted images from figure 2. If the wavefront error is uncorrected (a), the transmission quality decreases compared to the conventional correction (b) and diffraction limited case (c). The OTF only quantifies the transmission quality of the optical system for single spatial frequencies. For a single-value metric of the visual performance the transmission quality have to be combined with additional neural factors of the visual perception.

Visual strehl ratio
The visual perception is a complex process of optical transmission and neural processing [17]. The transmission quality of an optical system can be easily quantified with the OTF. However, the modeling of the neural processing is a challenging task. In 1965, Campbell and Green noted that the visual contrast sensitivity is the product of the optical transmission and a neural contrast sensitivity function (NCSF) [17]. This function is for each eye individual and characterizes the sensitivity of the neural perception of single spatial frequencies. The examination of this sensitivity is a time-consuming process where sinusoidal gratings are projected on the retina, which are not influenced by the optics [17,18]. For a series of spatial frequencies with variable contrast the minimal noticeable contrast is captured. This threshold is equal to the inverse of the neural contrast sensitivity function.
A general description of the neural contrast sensitivity function of 22 eyes was determined by Leube et al in [19]. The equation is given in (11) and the normalized values are plotted in figure 4. The peak of the function is at 6.36 cpd with a value of 138.356 9.
The spatial frequencies with the highest sensitivity are in the range of about 3 to 20 cpd. However, frequencies greater than about 60 cpd are almost not perceived. With the normalized neural contrast   [19]. sensitivity function each spatial frequency of the optical transmission can be weighted. This weighting is used in the single-value metric called visual Strehl ratio computed by the OTF (VSOTF), firstly introduced by Thibos et al [4]. The VSOTF is the ratio of the aberrated OTF H to the diffraction-limited OTF H 0 weighted with the neural contrast sensitivity respectively. Compared to other metrics, the VSOTF has the greatest correlation of objectively calculated and subjectively determined visual performance [2,4]. The VSOTF Q can be written as where The values of the VSOTF are within -1 and 1, whereby higher values correspond to a better visual quality. In our case of a wavefront correction, the VSOTF is a function of the correction Zernike coefficients ⃗ x. Thus, the resulting visual quality of a wavefront error with different corrections can be calculated. In figure 5 the convoluted images from figure 2 are quantified by the visual quality VSOTF. As shown in figure 5, the better the image quality, the higher the VSOTF Q. The different imaging qualities result from different methods of correction.

Correction types
Since the resulting visual quality of a wavefront error depends on the correction, the error will be corrected with four different correction types. These types are defined by their number and calculation method of the Zernike coefficients. The first correction is no correction at all, but the resulting wavefront corresponds to the pure wavefront error.
The second correction is a conventional correction W Conv that only uses low-order aberrations consisting of x-and y-astigmatism as well as defocus (c 3 , · · · , c 5 ). It corresponds to the correction employed today in spectacles and contact lenses. The relations between the Zernike coefficients and the refractive values of the eyeglasses prescription sphere, cylinder and axis are given in [7,20]. Looking at it mathematically, these low-order aberrations only compensate the quadratic terms of the polynomials. The remaining terms (c 6 , c 7 , · · · ) are uncorrected.
Our approach of a wavefront correction W Wave expands the conventional correction with four additional high-order aberrations. These are x-and y-trefoil as well as x-and y-coma (c 6 , · · · , c 9 ), such that cubic terms of the polynomials are also compensated. The remaining terms (c 10 , c 11 , · · · ) are uncorrected. Thus, the number of used coefficients to correct the wavefront error increase from 3 (conventional) to 7.
A full correction corresponds to a conjugate (negative) correction of all Zernike terms. It would correspond to a perfect wavefront. However, this correction is not realistic, as several factors (pupil adaption, shift of contact lens, accuracy of laser ablation) will prohibit it. In order to incorporate this case a last correction is implemented, which is directly calculated from the wavefront. Like the wavefront correction, this simplified conjugate correction W Conj also uses 7 coefficients for the correction. However, these coefficients correspond to the negative values of the wavefront error. Therefore, the respective coefficients of wavefront error and correction cancel each other and the resulting wavefront consists only of the uncorrected high-order aberrations with an index greater than 9.
= Z 10 · (c 10 + 0) + Z 11 · (c 11 + 0) + · · · The used Zernike coefficients of the different correction types are summarized in table 2. The pure and simplified conjugate correction are directly calculated by the wavefront error. However, the conventional and wavefront correction have variable coefficients that need to be optimized by an algorithm. The resulting visual performance is defined with the multidimensional objective function (12) and this results in different optimal solutions for the conventional as well as the wavefront corrections. Due to the complex and non-linear objective function the optimal solution cannot be calculated analytically but iteratively.

Particle swarm optimization
The optimization problem is to find an optimal correction for a specific visual defect. The correction is defined as Zernike wavefront so that the optimal Zernike coefficients need to be found. The resulting visual performance for a correction is approximated with equation (12). However, the optimal solution of this multidimensional, complex and non-linear objective function (12) has to be calculated iteratively. For this purpose, we apply the stochastic trial and error algorithm particle swarm optimization. This algorithm has already been introduced for classical optical optimizations where parameters of curvature radius, aspherical surface, lens distance etc of an optical system with several lenses are searched for [21][22][23]. An application in the ophthalmic field is so far unknown.
The algorithm was introduced first in 1995 by James Kennedy and Russell Eberhart [24]. Today there are many variants of the original algorithm, in our method we use the original one by Kennedy and Eberhart. The swarm consists of multiple particles and each particle is a possible solution of the problem. The number of particles depends on the specific optimization problem, but an initial number of 40 is suggested [25]. In our case, however, a swarm size of 20 particles has proven to be sufficient. The swarm P (k) with a size of 20 particles at iteration step k is described as with 1 ≤ i ≤ 20. The search space is bounded by the conjugate coefficients of the specific wavefront error with an addition of ± 5 microns for each coefficient.
On each iteration step, all particles move in the search space so that the new position (solution) is within the boundaries ⃗ x min ≤ ⃗ x (k+1) i ≤ ⃗ x max . The position of a particle corresponds to the Zernike coefficients of the correction, which is used for the calculation of the resulting visual quality. The new position of a particle is based on the previous position and is defined as The direction of the movement ⃗ v (k+1) i is calculated individually for each particle. This direction of movement consists only of the direction and is erroneously defined as velocity. But to keep the terminology, the direction of movement is also called velocity in this paper. The velocity is limited to 20% of the search space boundaries so that −2 ≤ v i ≤ 2 [26]. Like the position, the new velocity is based on the previous one and a common (social) ⃗ g (k) i experience. Each individual term is multiplied by a factor so that the cognitive and social influence is weighted differently. Depending on the weighting, the particles are either more obstinate and follow their individual direction or more cooperative and follow the best solution of the swarm. For this kind of optimization problem the values of Trelea (0.6, 1.7 and 1.7) have proven to be successful [27].
The cognitive and social terms are also multiplied with a random number r * cog and r * soc which are uniformly distributed within the range of 0 and 1. The average velocityv of all particles is the stop criteria of the algorithm. Once the average velocity isv ≤ 1.0 × 10 −4 , the swarm has converged and the algorithm stops. The initial (k = 0) positions and velocities of the particles are uniformly distributed within the boundaries ⃗ x min ,⃗ x max and ⃗ v min ,⃗ v max respectively.
The individual best solution ⃗ p (k) i represents the position the particle reached itself with the highest visual quality so far. The common best solution ⃗ g (k) i , on the contrary, is the position of the neighborhood with the highest visual quality. With the last term of equation (30), information is shared between the particles. The type of connection is called topology and the most recent three are illustrated in figure 6 [28]. In our case (in all other cases too), the goal is to use a topology that finds the optimal solution in a short period of time. Like the swarm size, the right choice of a topology depends on the specific optimization problem. A detailed description of the topologies can be found in [28]. To minimize the probability that the algorithm converges in a local maximum, we examine the three topologies: fully connected (global or gbest), ring (lbest) and star. In general, the global maximum of the objective function (12) is unknown, except the wavefront error is corrected completely. In this case, the correction corresponds to the complete conjugated wavefront error. To examine the topologies, an optimization is repeated 50 times and for each optimization the average visual quality of the solutions, standard deviation of the visual qualities, the average number of iterations as well as the success rate are recorded. After each optimization, the initial positions and velocities are recalculated. An optimization will be a success, if the resulting visual quality is 1 ± 0.000 1. As wavefront error, a reduced Zernike coefficients vector up to the 9th coefficient is used. These coefficients⃗ c and the optimal solution ⃗ x Opt as well as the lower and upper boundaries ⃗ x min ,⃗ x max are listed in table 3. The results of the optimizations are shown in table 4. The success rate of 64% of the fully connected topology from table 4 indicates that the objective function consists of several local maxima in which the particles converged. The ring and star topology on the other hand have a success rate of 100%, so that for each optimization the optimal solution was found. A closer look to the average iterations shows the largest number of iterations by the ring topology,   followed by the star topology and the global topology with the lowest number of iterations. This distribution of the number of iterations is also confirmed elsewhere [28]. Thus, the fastest but most inaccurate topology is not applicable for our optimization problem.
To ensure the convergence in the global optimum, we repeated the 50 optimizations, but with a complete Zernike coefficient vector of the wavefront error. For this setup the optimal solution is unknown and the success rate is not meaningful. Therefore, we focus on the highest and average resulting visual quality that are most significantly. For all topologies the highest calculated resulting visual quality is identical to 0.390 5, see table 5. However, the average values show that not for all topologies the solution with a visual quality of 0.390 5 was found for each optimization. This indicates that for these topologies the optimization occasionally converged in a local maximum. This examination is done for a wavefront correction with 7 dimension. To ensure the topology is also applicable for the conventional correction the last setup is repeated, but instead of 7 optimized coefficients only 3 were used. In the case of a conventional correction, all topologies, including the global and star topology, are converged in the optimal solution as shown in table 6. For our simulation, we therefore use the ring topology, which has the longest time for a convergence, but also the lowest probability for a convergence in a local maximum [29]. Table 6. Highest VSOTF Qmax, average VSOTFQ and number of iterations for 50 optimizations of a conventional correction. The initial positions and velocities were recalculated for each optimization.

Wavefront errors
The optimal solution of the objective function for the conventional and wavefront correction of a wavefront error can be calculated with the PSO described in the previous section. Several wavefront errors are needed to examine the different impact of the correction types on the resulting visual performance. In general, the wavefront error for healthy eyes mainly consists of low-order aberrations [7]. Thus, a conventional correction is sufficient in most cases. However, for eye disorders like keratoconus the wavefront error also consists of high-order aberrations that cannot be compensated with the conventional correction [8]. For our simulation we use 50 randomly generated wavefront errors with a coefficient distribution shown in figure 7. Compared to a distribution of healthy wavefront errors, the distribution from figure 7 consist of additional negative affecting high-order aberrations.

Results
The four correction types pure, conjugate, conventional and wavefront are calculated for each of the 50 randomly generated wavefront errors. The conjugate correction is directly calculated by the wavefront error and the conventional as well as the wavefront correction are optimized by the PSO. The pure wavefront error indicates the initial visual performance without any correction. For each wavefront error and correction type, the visual quality (VSOTF) of the optimal solution is recorded. An average of 451 iterations were needed to optimize the conventional correction by the PSO. However, the wavefront correction required more iterations with an average of 682. The distribution of the average visual quality is illustrated in figure 8. The average visual quality of the pure wavefront errors of 0.019 0 indicates the initial bad vision. As expected, the visual quality improves with the conjugate correction to an average of 0.077 5 compared to the pure wavefront error. Unexpectedly, the conventional correction increases the visual performance to an average of 0.151 4. As is well known, the conventional correction uses four coefficients less for the correction. But these three coefficients are optimized to the wavefront error by the PSO. Nevertheless, the best visual quality with an average of 0.352 4 is achieved by the wavefront correction, which is an improvement of more than 130 %.
The effects of the different correction types become clearer when a wavefront error is examined more closely. The Zernike coefficients of this wavefront error are listed in table 7.
The resulting retinal images of the four correcting types with the respective visual quality of the wavefront error from table 7 are shown in figure 9. The size of the letter E is about 5 arc minutes that is equal   The respective Zernike coefficients of the conjugate, conventional and wavefront correction are listed in table 8. The optimized coefficients of the conventional and wavefront correction are totally different from the conjugate values. Not only the absolute values change, also the sign. The conjugate correction simply uses the respective negative coefficients of the wavefront error. All others coefficients as well as neural factors of the visual perception are not involved. If the optimized coefficients of the conventional and wavefront correction are compared, it is noticeable that even the low-order aberrations differ from each other. The high-order aberrations therefore influence the compensating low-order aberrations in a complex way, so that the coefficients cannot be optimized individually, but have to be optimized simultaneously.

Discussion
Nowadays, there are different types of particle swarm optimization, which all are based on the used original PSO like Bare-bones PSO (BBPSO) [30], proactive particles in swarm optimization (PPSO) [31] or an effective fuzzy particle swarm optimization (EFPSO) [32]. The most of them have dynamically controlled boundaries for the search space or velocity. In general, the swarm converges faster and more probably in the global maximum. However, such an implementation is much more complex and will need more time for implementation. Because of the trial and error method of the particle swarm optimization, it cannot be guaranteed that the global maximum is reached for every optimization. But with the ring topology and parameters of Trelea, the probability is minimized that the swarm converges in a local maximum [28]. The VSOTF as objective function is an approximation of the complex visual perception that shows the greatest correlation of objectively calculated and subjectively measured visual performance so far. Fortunately, the objective function can be easily adapted for future possible better models.
The results show that the visual performance can be improved by extending the conventional correction with coma and trefoil aberrations. This effect is not astonishing, but the optimization showed that the visual performance of the optimized wavefront result is much more effective compared to a simple conjugate compensation. Theoretically, the order of correcting aberrations can be increased, but in reality especially the repeatability of high-order aberrations is low. These fluctuations are not due to the measuring systems, but are caused by minimal changes in the tear film, for example. Moreover, the correcting wavefront needs to be converted to a geometrical surface. Depending on the treatment of the eye defect, either a contact lens or an ablation profile for a refractive surgery must be produced. A contact lens or an ablation profile with such a high dynamic surface may be further a challenging task. Nevertheless, if it is necessary to use additional or other aberrations for the correction, the same objective function and optimization method can be used again.
The simulation was done for wavefront errors with comparable large high-order aberrations. Thus, the improvement from the conventional to the wavefront correction is immense. Because regular wavefront errors mainly consists of low-order aberrations, the simulation should lead to the same results, but with a much smaller improvement.
The conventional correction corresponds to a subjective refraction. The optimized coefficients of defocus as well as x-and y-astigmatism can therefore be converted to sphero-cylindrical refractive values of sphere, cylinder and axis. Thus, a subjective evaluation is not required. By optimizing only the 3 coefficients of the conventional scheme, the visual performance can be significantly improved compared to the non-optimized conjugate correction with 7 coefficients. As with subjective refraction, these low-order coefficients are optimized for the entire wavefront error. For example, defocus also partially corrects higher-order aberrations such as spherical aberration [4]. The individual coefficients therefore need to be optimized simultaneously for the entire wavefront error.

Conclusion
The complex and non-linear optimization problem for an objectively calculated correction for high irregular wavefront errors can be solved by the simple particle swarm optimization (PSO) method. The objective function for an optimal resulting visual performance is defined by the visual Strehl ratio computed by the OTF (VSOTF). So, optical and neural factors are considered for the visual perception. Different topologies of the PSO are examined and show that the objective function consists of several local maxima. To avoid a premature convergence, the ring (lbest) topology is used. With the same parameters, the algorithm can optimize corrections with low-as well as high-order aberrations.
The effects on the resulting visual performance of four different correction types are examined for 50 randomly generated high irregular wavefront errors. In this paper, the wavefront correction is an extension of the conventional correction with coma and trefoil aberrations, in x-and y-direction respectively. The coefficients of these corrections are optimized by the PSO.
The main superiority of the optimization process is illustrated by a comparison to the straightforward conjugate compensation, which is just using the negative wavefront coefficients. The comparison shows that the optimized conventional correction is more than 95% better the conjugate correction. We also show that the optimized wavefront correction increases the visual performance by more than 130% compared to the conventional correction. This corresponds to an increase of more than 350% relative to the conjugate correction.
Therefore, this optimization method can be used for an optimized conventional correction (without the need of a subjective evaluation) and it can also be used for an improved high-order correction scheme like free-form contact lenses or laser ablation profiles.