Light scattering in the human eye modelled as random phase perturbations

Augusto Arias; Harilaos Ginis; Pablo Artal

doi:10.1364/BOE.9.002664

1. Introduction

Physical and numerical eye models that incorporate intraocular scattering have been developed with two aims: the physical description of optical properties of the ocular media and the evaluation of the scattering effects on the quality of vision under various conditions, such as the effect of glare sources. The model proposed by Navarro [1] replicated an empirical profile of the ocular point spread function (PSF) with a Gaussian diffuser inserted in a schematic eye with aspheric surfaces. Van den Berg et al. [2] simplified the intraocular scattering sources to micro spheres analyzing forward light distribution from donor lenses. Under the same simplification, Chen et al. [3] introduced a biometry-based model with intraocular scattering. In this model, the spheres were located at the cornea and their size and concentration were optimized for adjusting the PSF to a glare equation suggested by the CIE (Commission International d’Eclairage) [4]. A similar method, involving spherical scatterers, was employed by Kelly-Pérez et al. [5]. However, the assumption that ocular scatterers are represented by a distribution of single size particles is not completely realistic. Moreover, the angular distribution of scattering associated to such particles is generally different from the CIE glare function [5].

In addition, static and dynamic optical devices have been reported with the intention of reproducing the straylight in healthy and cataractous eye. Paeglis et al. [6] described a variable diffuser based on polymer dispersive liquid crystal as part of a physical eye model with controllable scattering. However, there is no quantitative assessment of the amount of straylight in this device. Wit et al. [7] evaluated different diffusing filters and concluded that the Black Pro-Mist filters (BPM; The Tiffen Company, NY) provides straylight values corresponding from a minor visual hindrance to mild cataract. These filters, originally designed for decreasing the contrast and sharpness in professional photography, are commonly considered as a reference in the estimation and production of intraocular scattering [8–11].

In this article, the intraocular scattering effects are reproduced by a random phase map -located at the pupil plane- calculated from the weighted combination of wavefronts with different spatial frequencies. This is similar to that presented by Ginis et al. [12], comparing the spatial properties of the confocal microscope image from a rabbit cornea and a fractal surface. This kind of phase maps generate straylight amounts similar to those in the human eye. The amount of straylight can be controlled through the wavefront amplitude [13,14]. However, a model based on these random surfaces requires a complete definition of its parameters. With this motivation, we developed an alternative approach where two parameters are introduced for controlling the wavefront shape and their optimization allows to reproduce the average angular profile of the scattered light distribution.

For the experimental validation of the model, wavefronts were generated using a spatial light modulator (SLM) and the produced amounts of straylight were evaluated using the optical integration method [12]. The capability and usefulness of the model is discussed from both numerical and experimental results.

2. Methodology for wavefront calculation

Our goal is to design phase maps (W) whose angular distribution of their diffraction patterns mimics the intraocular straylight in a large population of eyes. The angular distribution of the wide-angle ocular PSF, has been estimated using various psychophysical methods and it is summarized using the following empirical formula [4]:

\begin{array}{l} P S F_{C I E} = [1 - 0.08 \cdot {(\frac{A}{70})}^{4}] \cdot [\frac{9.2 \times 10^{6}}{{[1 + {(θ / 0.0046)}^{2}]}^{1.5}} + \frac{1.5 \times 10^{5}}{{[1 + {(θ / 0.045)}^{2}]}^{1.5}}] \\ + [1 - 1.6 \cdot {(\frac{A}{70})}^{4}] \cdot {[\frac{400}{1 + {(θ / 0.1)}^{2}} + 3 \times 10^{- 8} . θ^{2}] + p \cdot [\frac{1300}{{[1 + {(θ / 0.1)}^{2}]}^{1.5}} + \frac{0.8}{{[1 + {(θ / 0.1)}^{2}]}^{0.5}}]} \\ + 2.5 \times 10^{- 3} \cdot p \end{array}

where θ is retinal angle in degrees, A is the age and p is a pigmentation coefficient that equals to 0, 0.5 and 1 for black, brown and light eyes respectively. Without loss of generality we assume p = 1 for the following calculations.

The method, summarized in Fig. 1, is based on the Inverse Discrete Cosine Transform (IDCT) of standard normally distributed random numbers (matrix R) that are weighted by an appropriate power law function (U). In the context of the IDCT, matrix R represents random coefficients in different spatial frequencies while matrix U represents the relative contribution of these different spatial frequencies in the resulting phase map. Matrix R can be re-randomized to create new phase maps that -although different- have the same spectral content as determined by the law in matrix U. In particular, the wavefront W is obtained by applying the IDCT to the element-wise multiplication between U and R:

W_{i', j'} = circ (i', j') \sum_{i = 0}^{N - 1} \sum_{j = 0}^{N - 1} U_{i j} R_{j i} \cos (\frac{π ϕ f_{i}}{N} i') \cos (\frac{π ϕ f_{j}}{N} j')

where N is the dimension of matrices and circ() is the aperture function, valued 1 for

2 \sqrt{{(i^{'} - 0.5 N)}^{2} + {(j^{'} - 0.5 N)}^{2}} \leq N

and zero otherwise, acting as iris with diameter ϕ.

Fig. 1 Phase map calculation.

Download Full Size | PDF

Each position (i,j) in the spectrum represents a mode, i.e., a harmonic phase map in the spatial domain (i’,j’) with single vertical $f_{i} = (i + 1 / 2) / ϕ$ and horizontal $f_{j} = (j + 1 / 2) / ϕ$ frequencies. In this way, the wavefront is the sum of all cosine modes where the amplitude of each one is related to the fraction of the incoming energy into the eye directed to a specific angular position of the PSF.

The modulating matrix U is calculated using a power-law function applied on the radial frequency $f = \sqrt{f_{i}^{2} + f_{j}^{2}}$ (in cycles per millimeter):

F (f) = B f^{β}

where B and β are the model parameters. While β controls the relative weighting of the modes (i.e., the energy distribution between small and large angles in the PSF), B controls the root mean square (RMS) amplitude of the wavefront. Thus, B controls the straylight amount also, due to the linear relationship between logarithmic values of the root mean squared (RMS) and the straylight parameter [13]. For each set of parameters (B,β), the PSF is calculated as the Fraunhofer diffraction pattern associated to the wavefront by applying the fast Fourier transform [15]. The model optimization pertains to the estimation of the parameters B and β, that produce wavefronts for which the PSF approximates optimally the empirical PSF_CIE.

Each calculated PSF was normalized per unit of solid angle (see Appendix). The straylight amount is quantified as the straylight parameter S defined as S = θ² × PSF(θ) at an angle of 6 degrees.

The particular (B,β) couple that reproduces the angular distribution of the ocular straylight was found by an optimization procedure that minimizes the RMS difference between the logarithmic values of the associated PSF and PSF_CIE. This procedure was implemented using the fminunc() function in MATLAB (MathWorks, Natick, MA, USA), which is based on the quasi-Newton method [16]. Once the parameters are optimized, they can be use with different values in R without modifying the average angular profile of the scattered light distribution because those random values slightly change the position of the speckles in the PSF.

Initially the B and β values are optimized for representing the PSF of a young eye, i.e. assuming A≡30 in the Eq. (1). Afterwards, the straylight is evaluated for several B values in order to estimate the relationship between them. The random values in R are preserved in this process. A demonstration of the model performance was developed with the particular features of the experimental setup described in the following Section.

3. Experimental generation and evaluation of straylight

Figure 2 shows the experimental setup used for evaluating the straylight produced by the designed phase maps.

Fig. 2 Experimental setup for the dynamical generation and estimation of straylight.

Download Full Size | PDF

The generation of scattered PSFs is mainly carried out by a Liquid Crystal on Silicon (LCoS) SLM (PLUTO; Holoeye, AG, Germany). Its phase modulation performance was improved by the inclusion of a linear polarizer and a green filter (λ≡540 ± 10 nm).

The pixel pitch of the SLM (8 μm) limits the angular range of the PSF to 1.93 degrees, which is not enough for an appropriated representation of the intraocular scattering. Therefore, a telescope with 6x angular magnification was introduced in the setup. It simultaneously conjugates the SLM plane with the pupil of an artificial eye model, which is composed of a biconvex lens and a 14-bits CCD camera (Luca; Andor, Belfast, UK), acting as ocular optics and retina respectively. In this experience, the size (N) of the matrix displayed on the SLM is 1000 pixels and it produces a pupil size (ϕ) of 1.33 mm. In this configuration, each pixel of the SLM corresponds geometrically to a size of approximately 1.3 μm at the pupil plane.

Implementation of the optical integration method

In parallel, the induced scattering is evaluated using the optical integration method, by recording the projected disks of uniform irradiance and angular wide 2θ from a monitor placed in front of the SLM. The radial profile of the incoherent PSF is calculated from the recorded disks using its relationship with their central intensities I_c (following the notation shown in Fig. 2):

I_{c} (θ) = \int_{0}^{θ} 2 π P S F (ϕ) d ϕ

It means that each point source localized on φ radius has an intensity contribution at the center of the disk that equals to PSF(φ) and therefore the complete annular ring provides a total energy of 2πφPSF(φ). The integration of φ values from 0 to θ allows to predict the effect of a complete disk. Thus, the PSF at θ angle is retrieved calculating the angular derivative of Eq. (4):

P S F (θ) = \frac{1}{2 π θ} \frac{d I_{c} (θ)}{d θ}

A simplified function is proposed by the approximated PSF estimation from the recorded central intensities:

P S F (θ) = C (P S F_{d l}) + D [\frac{1}{{(θ + 0.001)}^{E}}]

where PSF_dl is the diffraction-limited PSF, calculated based on the apertures of the system and radially normalized to 1. The C, D, and E parameters are optimized to minimize the mean squared error between experimental and numerical (calculated replacing Eq. (6) in 4) central intensities data. Moreover, the parameters are adjusted for fulfilling the normalization condition of the encircled energy by the calculated PSF on the considered τ angular range, i.e.,

\int_{0}^{τ} 2 π θ P S F (θ) d θ = 1

. As explained above, due to SLM features and telescope magnification, τ is equivalent to 11.9 degrees. This procedure can fit experimental data with r² higher than 0.990. Finally, straylight amount is calculated from the retrieved PSF.

In this way, the phase maps for the addition of six straylight amounts ranging from 10 to 60 degree²/sr were designed and their experimentally performance was assessed.

4. Results

Initially, the model was implemented considering the above described experimental features for replicating the angular course of the ocular PSF corresponding to a 30 years old eye. The optimized values of B (9.207 µm) and β (−1.214) produce the phase map shown in Fig. 3(a). The angular average of the optimized PSF is compared with the CIE reference, for the considered age, in Fig. 3(b).

Fig. 3 Numerical optimization results.

Download Full Size | PDF

Figure 4 shows the numerical relationship between the generated straylight at six degrees and the RMS of the calculated wavefront for several B values. Following this relationship, the angular course of the PSF with a high amount of straylight (80 degrees²/sr at six degrees) was calculated and shown in the Fig. 3(b). In addition, this numerical profile is compared with the CIE reference that generates the programmed amount of straylight.

Fig. 4 Relationship among B parameter, RMS and the straylight at 6 degrees.

Download Full Size | PDF

The additional straylight amounts experimentally assessed when the designed phase maps were displayed on the SLM are depicted by red dots in Fig. 4. It corresponds to a subtraction between the absolute measurement of straylight and the measured one when no-additional scattering is induced in the setup (i.e., a flat wavefront displayed on the SLM). This background straylight is 9.12 degrees²/sr and corresponds to the straylight from all the optics in the experimental setup including the SLM.

5. Discussion

A methodology for designing phase maps that replicate the angular distribution of intraocular straylight has been developed. The phase maps are calculated as the inverse discrete cosine transform of the modulation of random values by a power-law function. A numerical optimization allows to determine the function parameters that replicates the typical PSF of a healthy and young eye (i.e., with low intraocular scattering). Further straylight amounts can be generated by changing the overall amplitude of the cosine modes. Moreover, the statistical representation of an eye population with near straylight amounts could be simulated by the simple change of the random values while the optimized model parameters are held.

The estimated coefficient of the power law in the Eq. (3) (β = −1.214) classifies the estimated wavefront pattern (see Fig. 3(a)) as pink noise which is a common feature in physical and biological systems [17] and, particularly, it was employed for the scattering description of biological tissues [18]. This kind of wavefronts are composed by periodic signals with a wide bandwidth that simultaneously represent the phase perturbations due to high-order Zernike aberrations and typical intraocular scatterers such as cell nuclei, mitochondrias, organelles and protein aggregates. This feature allows the continuous replication of the scattering effects on the ocular PSF along a wide angular range, as it was depicted in Fig. 3(b), unlike the intraocular scattering models based on single size particles [5]. Moreover, as shown in Fig. 4, the amount of straylight can be manipulated in a control way through the amplitude of one calculated phase map, preserving the angular course of the PSF determined by the CIE formula (see Fig. 3(b)). This relationship is linear considering the logarithmic values of RMS and straylight.

The programmed straylight amounts were experimentally verified using a setup with a liquid crystal display as spatial phase modulator. Optical integration method was implemented for the measurement of straylight amount. The spatial properties of the display and a telescope introduced in the system allows to represent phase perturbations with periods ranged from 2.67 to 1333.33 μm, which correspond to an angular domain of the PSF from 0.023 to 11.6 degrees (radius). According to Fig. 4, there is agreement between the programmed amount of straylight and the experimentally measured additions of straylight. However, it is necessary to point out the limitation of the experimental setup due to the addition of a constant offset of straylight on all programmed amounts. This offset, equivalent to the amount of straylight expected for a 55 years-old eye (according to the Eq. (1)), is originated by the diffusion properties of the nematic liquid crystal that composes the SLM [19], scattered light in the others optical components and parasitic light reflections.

Beyond the shown implementation adopting the parameters of the experimental setup, the cosine base may be a powerful mathematical representation for a complete statistic representation of ocular aberrations, including phase perturbations with lowest spatial frequencies. Furthermore, the accuracy generation of the intraocular straylight reveals that the proposed methodology could be useful for further psychophysical experiments where controllable straylight amounts are required, e.g., the cataract simulation.

Appendix normalization of PSF

The units of the PSF_CIE, calculated using the Equation 1, are inverse steradians (sr⁻¹), i.e., the empirical PSF is normalized to the Ω solid angle covered through the propagation [20]. In consequence, the numerical PSFs to be compared with that adopted reference must fulfill the following condition by means of its multiplication with a normalization factor α:

\int_{Ω} α P S F d Ω = 1

The solid angle is defined as the ratio between the area on a sphere and the square of its radius. In the paraxial approximation, the assumption that supports the propagation algorithm, the pixel pitch in the PSF matrix is equivalent to $λ z / ϕ$ being z the propagation distance [15]. Thus, the subtended solid angle by each (x,y) pixel of the PSF matrix is simplified to $λ^{2} / ϕ^{2}$ . Therefore, α is calculated as:

α = {[{(\frac{λ}{ϕ})}^{2} \sum_{x = 1}^{N} \sum_{y = 1}^{N} P S F_{x, y}]}^{- 1}

Funding

European Research Council (ERC) (AdG−2013−339228 SEECAT); Secretaría de Estado e Investigación, Desarrollo e Innovación (SEIDI) (FIS2013-41237-R); Fundación Séneca-Agencia de Ciencia y Tecnología de la Región de Murcia (19897/GERM/15).

Acknowledgments

A. Arias acknowledges the support from the Spanish Government (FPI-BES-2014-070427).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

References and links

1. R. Navarro, “Incorporation of intraocular scattering in schematic eye models,” J. Opt. Soc. Am. A 2(11), 1891–1894 (1985). [CrossRef] [PubMed]

2. T. J. van den Berg, “Light scattering by donor lenses as a function of depth and wavelength,” Invest. Ophthalmol. Vis. Sci. 38(7), 1321–1332 (1997). [PubMed]

3. Y.-C. Chen, C.-J. Jiang, T.-H. Yang, and C.-C. Sun, “Development of a human eye model incorporated with intraocular scattering for visual performance assessment,” J. Biomed. Opt. 17(7), 0750091 (2012). [CrossRef] [PubMed]

4. J. J. Vos and T. J. T. P. van den Berg, “Report on disability glare,” CIE Collect. 135, 1–9 (1999).

5. I. Kelly-Pérez, N. C. Bruce, L. R. Berriel-Valdos, A. Werner, and J. A. Delgado Atencio, “Computational model of the effect of light scattering from cataracts in the human eye,” J. Opt. Soc. Am. A 30(12), 2585–2594 (2013). [CrossRef] [PubMed]

6. R. Paeglis, M. Ozolinsh, P. Cikmacs, and S. Andersson-Engels, “Eye model with controllable lens scattering,” in (2001), Vol. 4434, pp. 233–238.

7. G. C. de Wit, L. Franssen, J. E. Coppens, and T. J. T. P. van den Berg, “Simulating the straylight effects of cataracts,” J. Cataract Refract. Surg. 32(2), 294–300 (2006). [CrossRef] [PubMed]

8. H. Ginis, O. Sahin, A. Pennos, and P. Artal, “Compact optical integration instrument to measure intraocular straylight,” Biomed. Opt. Express 5(9), 3036–3041 (2014). [CrossRef] [PubMed]

9. H. Ginis, G. M. Pérez, J. M. Bueno, and P. Artal, “The wide-angle point spread function of the human eye reconstructed by a new optical method,” J. Vis. 12(3), 20 (2012). [CrossRef] [PubMed]

10. P. H. B. Kok, H. W. van Dijk, T. J. T. P. van den Berg, and F. D. Verbraak, “A model for the effect of disturbances in the optical media on the OCT image quality,” Invest. Ophthalmol. Vis. Sci. 50(2), 787–792 (2009). [CrossRef] [PubMed]

11. D. Christaras, H. Ginis, A. Pennos, and P. Artal, “Intraocular scattering compensation in retinal imaging,” Biomed. Opt. Express 7(10), 3996–4006 (2016). [CrossRef] [PubMed]

12. H. Ginis, I. Pentari, D. de Brouwere, D. Bouzoukis, I. Naoumidi, and I. Pallikaris, “Narrow angle light scatter in rabbit corneas after excimer laser surface ablation,” Ophthalmic Physiol. Opt. 29(3), 357–362 (2009). [CrossRef] [PubMed]

13. H. Ginis, J. J. Rozema, M.-J. B. R. Tassignon, and P. Artal, “A phase perturbation model of light scattering in the human eye,” Invest. Ophthalmol. Vis. Sci. 55, 2120 (2014).

14. A. Arias, H. Ginis, and P. Artal, “Physical Model of Intraocular Scattering using a Spatial Light Modulator,” Invest. Ophthalmol. Vis. Sci. 56, 1066 (2015).

15. D. Mas, J. Garcia, C. Ferreira, L. M. Bernardo, and F. Marinho, “Fast algorithms for free-space diffraction patterns calculation,” Opt. Commun. 164(4-6), 233–245 (1999). [CrossRef]

16. T. Coleman, M. A. Branch, and A. Grace, “Matlab Optimization Toolbox,” Change, 42 (2008).

17. P. Szendro, G. Vincze, and A. Szasz, “Pink-noise behaviour of biosystems,” Eur. Biophys. J. 30(3), 227–231 (2001). [CrossRef] [PubMed]

18. M. Xu and R. R. Alfano, “Fractal mechanisms of light scattering in biological tissue and cells,” Opt. Lett. 30(22), 3051–3053 (2005). [CrossRef] [PubMed]

19. M. H. Kao, K. A. Jester, A. G. Yodh, and P. J. Collings, “Observation of light diffusion and correlation transport in nematic liquid crystals,” Phys. Rev. Lett. 77(11), 2233–2236 (1996). [CrossRef] [PubMed]

20. T. J. van den Berg, “Analysis of intraocular straylight, especially in relation to age,” Optom. Vis. Sci. 72(2), 52–59 (1995). [CrossRef] [PubMed]

Light scattering in the human eye modelled as random phase perturbations

Abstract

1. Introduction

2. Methodology for wavefront calculation

3. Experimental generation and evaluation of straylight

Implementation of the optical integration method

4. Results

5. Discussion

Appendix normalization of PSF

Funding

Acknowledgments

Disclosures

References and links

Cited By

Figures (4)

Equations (8)

Biomedical Optics Express