Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Optimizing near-infrared polariscopic imaging for the living human eye

Open Access Open Access

Abstract

Hardware architectures and image interpretation can be simplified by partial polarimetry. Mueller matrix (MM) polarimetry allows the investigation of partial polarimeter designs for a given scientific task. In this work, we use MM measurements to solve for a fixed polarization illumination and analyzer state that maximize polariscopic image contrast of the human eye. The eye MM image acquisition takes place over 15 seconds which motivates the development of a partial polarimeter that has snapshot operation. Within the eye, the birefringent cornea produces spatially-varying patterns of retardance exceeding half of a wave with a fast-axis varying from linear, to circular, and elliptical states in between. Our closed-form polariscopic pairs are a general solution that maximizes contrast between two non-depolarizing pure retarder MMs. For these MMs, there is a family of polariscopic pairs that maximize contrast. This range of solutions creates an opportunity to use the distance from optimal as a criteria to adjust polarimetric hardware architecture. We demonstrate our optimization approach by performing both Mueller and polariscopic imaging of an in vivo human eye at 947 nm using a dual-rotating-retarder polarimeter. Polariscopic images are simulated from Mueller measurements of 19 other human subjects to test the robustness of this optimal solution.

© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

The birefringence of human eyes offers an interesting and spatially-varying sample for polarimetric imaging. Corneal birefringence is well known, and is the result of anisotropic collagen fibril orientation [14]. Corneal birefringence has been observed using multiple modalities, including polarization-sensitive optical coherency tomography [510]. Of particular relevance is in vivo spatially-resolved polarization imaging between variously-oriented polarizers [1114]. These works as well as Mueller matrix (MM) imaging have been performed in visible wavebands [15]. Much of the literature indicates that the center region of the cornea exhibits linear birefringence, with any circular birefringence being relatively small if at all present [15]. However, our near-infrared MM imaging of 20 human subjects’ entire cornea demonstrates the presence of significant circular birefringence, which has not been reported before to our knowledge [16]. The spatially varying birefringence over the cornea observed in Mueller eye measurements from 20 human subjects is used to demonstrate our approach for optimizing the contrast between two generally elliptical non-depolarizing pure retarder MMs.

MM imaging provides a complete description of a polarized light-matter interaction. When MM imaging is prohibitively complex, expensive, or even unnecessary, partial polarimetry can be used to extract useful polarimetric information [1719]. If some amount of information is known a priori about the subset of MMs to be measured, it becomes possible to compare different partial polarimeter designs. The simplest sample-measuring partial polarimeter is a polariscope, where an object is illuminated by a fixed polarized state generator and imaged through another, generally different, fixed polarization state analyzer. Polariscopic measurements are commonly used for analyzing birefringence in crystals, stress induced by the photoelastic effect in transparent materials, and the concentration of chiral molecules in liquids [20]. Methods for numerical optimization of the polarized illumination and analyzer states can be found in the literature [21,22]. For example, an analytic solution for optimizing the contrast in polariscopic measurements between an object and its background has been studied [23]. These authors analyzed the illumination and analyzer states for MMs with the same retarder fast-axis but varying retardance magnitude and for MMs with the same diattenuation magnitude but varying diattenuation orientation. They point out that optimal illumination states are those which have maximal separation on the Poincaré sphere after the light-matter interaction. They also found that for the retarder MMs, the set of optimal illumination states define a great circle on the Poincaré sphere. Our approach is applicable to any pair of retarder MMs which vary in retardance magnitude, orientation, or both. With respect to corneal birefringence, our approach may be useful for applications such as eye segmentation or eye tracking. The spatially-varying birefringence in the cornea is turned into spatially-varying brightness features, and maximizing the contrast could make such features easier to track.

In our work, polariscopic pairs of a generator and analyzer are optimized for contrast between birefringent MMs of arbitrary fast axes and retardance magnitudes. Over the course of 15 seconds, 25 polarimetric images with different generator and analyzer states were used to reconstruct the initial MM image. The inability to measure an in vivo human eye for an appreciable time duration without motion artifacts motivates our investigation of polariscopic eye imaging. The polariscopic pair optimized from one human subject is used on MM images of 19 other in vivo human eyes for a qualitative comparison of the contrast pattern. The closed-form solution we provide is based on a geometric construction on the Poincaré sphere. As a demonstration of our method, near-infrared polariscopic measurements of the in vivo human eye are simulated from Mueller measurements and compared to polariscopic measurements at nearly-optimal configurations in an existing dual-rotating-retarder polarimeter at 947 nm. This near-infrared wavelength was used because there is less variation in reflectance from the iris than for visible wavelengths among individuals with different eye colors [24]. Exactly optimal polariscopic pairs are not necessarily available in a dual-rotating-retarder polarimeter. Polarimeters with different polarization state generator and analyzer architectures, such as those with liquid-crystal variable retarders, could have the flexibility to access a wider range of optimal polariscopic pairs [2530].

2. Mueller matrix measurements

In Mueller calculus, the polarization state of light is described by the Stokes vector which contains 4 real-valued elements. The $4\times 4$ real-valued MM elements describe the linear equations relating the Stokes vector before and after a light-matter interaction. Knowledge of a MM allows prediction of subsequent polarimetric measurement given arbitrary illumination and analyzer states. A polarization-sensitive measurement $p$ is expressed as

$$p=\mathbf{a}\mathbf{M}\mathbf{g}$$
where the Stokes column vector, $\mathbf {g}$, is the illumination state which is also called the generator, the analyzer is a Stokes row vector, $\mathbf {a}$, quantifying the polarization sensitivity of the imaging optics, and $\mathbf {M}$ is the MM where the dependence on wavelength and scattering geometry have been dropped for brevity. In this work, the MMs of interest represent the near-infrared polarization properties of different regions of an in vivo image of the human eye.

The MM eye images were taken using an existing dual-rotating-retarder Mueller imaging polarimeter named the RGB950 as shown in Fig. 6 in the Appendix [31]. The RGB950 operates at three visible wavebands and one near-infrared waveband. The near-infrared waveband, centered at 947 nm with a full width half max of 40 nm, was used in this work. The eyes were measured in a reflective double-pass configuration where light transmitted through the cornea, reflected off the iris, and transmitted back through the cornea a second time. Over the course of 15 seconds 25 polarimetric images with different generator and analyzer states were used to reconstruct the MM image. Due to the small unconscious random movements of the eye, image registration was performed on the 25 images before the MM reconstruction. The reconstructed MM image from one of 20 human subjects is shown in Fig. 1(a). Remaining motion artifacts on the order of several pixels are to be expected, but more slowly-varying polarimetric patterns are also observed. The inability to measure an in vivo human eye for an appreciable time duration motivates our investigation of snapshot polariscopic eye imaging.

 figure: Fig. 1.

Fig. 1. Near-infrared MM measurements of an in vivo human eye, masked to the region of the cornea where light is reflected back from the iris. In (a), the image is normalized such that $M_{00}=1$ to visualize the spatially-varying polarization properties. In (b), the MM measurement is approximated as TD using Eq. (2) with (c) being the dominant MJM $\widehat {\mathbf {M}}_0$. In (d), the dominant MJM is further approximated as a pure retarder MJM.

Download Full Size | PDF

MMs have 16 degrees of freedom which are associated with overall reflectance, diattenuation, retardance, and depolarization [20]. A method to reduce the depolarization degrees of freedom from nine to one is the triple-degenerate (TD) model, where the MM is approximated using a non-depolarizing Mueller-Jones matrix (MJM) and an ideal depolarizer which perfectly depolarizes all incident polarization states [32]. The MM is approximated as a convex sum of two terms

$$\mathbf{M}\approx\frac{4M_{00}}{3}\left[\left(\xi_0-\frac{1}{4}\right)\widehat{\mathbf{M}}_0+(1-\xi_0)\mathbf{M}_{ID}\right],$$
where $\widehat {\mathbf {M}}_0$ is the MJM, $\mathbf {M}_{ID}$ is the ideal depolarizer, $M_{00}$ is the overall reflectance, and $\xi _0$ is a scalar-valued depolarization parameter which controls the relative weights in the convex sum. The TD decomposition is a special case of the integral decomposition of a MM, both of which are derived from the Cloude spectral decomposition [3335]. TD models have proven useful for partial polarimetric experiments [36,37]. The TD approximation of the near-infrared eye MM image is shown in Fig. 1(b) with the dominant MJM image shown in Fig. 1(c). The appropriateness of a TD approximation for these measurements is discussed in the Appendix.

To simplify the optimization of generator and analyzer states, the eye MM is further approximated as having a dominant MJM $\widehat {\mathbf {M}}_0$ which is a pure retarder, see Fig. 1(d). The retarder MM is separated from the diattenuation by taking the unitary part from the polar decomposition of the dominant MJM. A general elliptical retarder MM is characterized by its retardance vector $\vec {\boldsymbol {\delta }}=[\delta _H,\delta _{45},\delta _R]$, where the magnitude of the vector is the retardance magnitude and the direction of this vector is the fast-axis on the Poincaré sphere. Due to phase wrapping, the retarder vector calculated from a MM is non-unique. For integers $q$, the retarder vector $\vec {\boldsymbol {\delta }}^{q}=\left (2\pi q+\|\vec {\boldsymbol {\delta }}\|\right )\vec {\boldsymbol {\delta }}/\|\vec {\boldsymbol {\delta }}\|$ will produce the same MM as $\vec {\boldsymbol {\delta }}$. In this work, $q$ is assigned pixel-wise to produce smoothly-varying retardance vectors over the image. The specific value of $q$ does not affect the optimization of the illumination and analyzer states.

The MM for an elliptical retarder is given by

$$\mathbf{M}_{ER}(\vec{\boldsymbol{\delta}})=\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & \frac{\delta_H^2+\left(\delta_{45}^2+\delta_R^2\right)C}{\delta^2} & \frac{\delta_{45}\delta_{H}T}{\delta^2}+\frac{\delta_RS}{\delta} & \frac{\delta_H\delta_RT}{\delta^2}-\frac{\delta_{45}S}{\delta}\\ 0 & \frac{\delta_{45}\delta_{H}T}{\delta^2}-\frac{\delta_RS}{\delta} & \frac{\delta_{45}^2+\left(\delta_R^2+\delta_H^2\right)C}{\delta^2} & \frac{\delta_R\delta_{45}T}{\delta^2}+\frac{\delta_HS}{\delta}\\ 0 & \frac{\delta_H\delta_RT}{\delta^2}+\frac{\delta_{45}S}{\delta} & \frac{\delta_R\delta_{45}T}{\delta^2}-\frac{\delta_HS}{\delta} & \frac{\delta_R^2+\left(\delta_{45}^2+\delta_H^2\right)C}{\delta^2} \end{bmatrix}=\begin{bmatrix} 1 & \vec{\mathbf{0}}^\mathrm{T}\\ \vec{\mathbf{0}} & \mathbf{U}(\vec{\boldsymbol{\delta}}) \end{bmatrix},$$
where $\delta =|\vec {\boldsymbol {\delta }}|$, $C=\cos (\delta )$, $S=\sin (\delta )$, and $T=1-\cos (\delta )$. In this equation, $\vec {\mathbf {0}}$ is a $3\times 1$ vector of zeros and the lower-right $3\times 3$ elements of the MM are written as a unitary matrix $\mathbf {U}(\vec {\boldsymbol {\delta }})$.

Images of each retardance component are shown in Fig. 2 for a single human subject’s eye. The first two retarder vector components, shown in Fig. 2(a) and (b), correspond to the linear retardance. For the majority of the masked area, the $\delta _H$ term has greater magnitude than the $\delta _{45}$ term. Therefore, the linear retardance is primarily horizontally or vertically oriented (where $\delta _H$ is positive or negative, respectively). The circular component, in Fig. 2(c), shows significant circular retardance which has not been previously reported to the authors’ knowledge. The change from positive to negative values in $\delta _R$ indicates a change in the helicity of the fast axis of retardance. The retardance magnitude in Fig. 2(d), which is the norm of the retardance vector, shows two minima on the left and right of the masked region and a maximum at the bottom of the eye. A second maximum likely appears at the top of the eye but is obscured by the eyelid in this image. The circular retardance result highlights the importance of performing the MM characterization as an initial step in the optimization described in the following section. Without first performing the Mueller image characterization, we would have used the purely linear birefringence assumption found in the literature for the optimization. Trying to perform the optimization of the polariscopic pairs based on an assumed MM without the circular birefringence would have resulted in non-optimal pairs.

 figure: Fig. 2.

Fig. 2. The horizontal, diagonal, and circular retardance computed from a measured MM are shown in (a), (b), and (c), respectively. The total retardance magnitude is shown in (d). ROIs are designated $i$-$l$. In (e), the retardance vectors within each ROI are plotted in 3D, where the concentric spherical surfaces denote $\pi /2$, $\pi$, and $3\pi /2$ radians of retardance. The point clouds represent the per-pixel retardance vectors and the arrows indicate the average retardance vector for each ROI. Average values are listed in Table 1.

Download Full Size | PDF

Tables Icon

Table 1. Components and magnitudes of the average retarder vectors from the four ROIs shown in Fig. 2 in units of radians. Retardance magnitudes are known up to the order of retardance, where $2\pi q$ for integers $q$ can be added to the magnitude without changing the resulting MM.

To reduce the impacts of measurement noise on the generator and analyzer optimization, the retarder vectors were averaged over each of the regions of interest (ROIs) shown as white boxes in Fig. 2. The average retarder vectors appear in Table 1. The retarder vectors were averaged rather than the MJMs because the summation of different MMs introduces depolarization which is undesirable in this application. ROIs $i$ and $j$ were chosen for use in the following optimization because they have different retardance orientations. The optimization could still be performed for ROIs with similar retardance vectors, such as $j$ and $k$, but the maximum achievable contrast would be lower.

3. Optimizing the polarimeter configuration

3.1 Optimal Illumination and Analyzer States

In this work, the contrast between polarimetric measurements of two MMs is defined as the magnitude of the measurement difference. This contrast $\Delta p$ is written

$$\Delta p = ||\mathbf{a}(\mathbf{M}_i-\mathbf{M}_j)\mathbf{g}||,$$
where $\mathbf {M}_i$ and $\mathbf {M}_j$ are unique, e.g. different pixels or regions of a MM image. The maximum achievable contrast between two polarization states increases as their distance on the Poincaré sphere increases [23]. Therefore, the optimal illumination states $\mathbf {g}$ are those for which, $\mathbf {s}_i=\mathbf {M}_i\mathbf {g}$ and $\mathbf {s}_j=\mathbf {M}_j\mathbf {g}$, are exitant Stokes vectors with the greatest angle possible between them on the Poincaré sphere. Assuming that $\mathbf {M}_i$ and $\mathbf {M}_j$ are pure retarders (see Sect. 2.), maximizing this angle is equivalent to minimizing the inner product
$$\vec{\mathbf{s}}_j^\dagger\vec{\mathbf{s}}_i=\left(\mathbf{U}_j\vec{\mathbf{g}}\right)^\dagger\left(\mathbf{U}_i\vec{\mathbf{g}}\right)=\vec{\mathbf{g}}^\dagger\left(\mathbf{U}_j^\dagger\mathbf{U}_i\vec{\mathbf{g}}\right)$$
where the vector arrow notation indicates a $3\times 1$ vector quantity on the Poincaré sphere, as opposed to a $4\times 1$ Stokes vector (i.e., $\mathbf {g}=S_0[1,s_1,s_2,s_3]=S_0[1,\vec {\mathbf {g}}]$) and $\mathbf {U}$ are $3\times 3$ unitary matrices from Eq. (3). As shown in Eq. (5), the inner product is invariant to multiplying both $\vec {\mathbf {s}}_j$ and $\vec {\mathbf {s}}_i$ by the unitary transform $\mathbf {U}^\dagger _j$, resulting in an inner product between the original generator vector $\vec {\mathbf {g}}$ and that generator vector after transformation by the new matrix $\mathbf {U}_j^\dagger \mathbf {U}_i$. On the Poincaré sphere, a retarder behaves as a rotation. The composition of two rotations is another rotation, so $\mathbf {U}_j^\dagger \mathbf {U}_i$ can be analyzed as any other elliptical retarder matrix.

Points which are further from the rotation axis are rotated a greater distance. States on the surface of the Poincaré sphere which are 90$^\circ$ from the fast axis are rotated the most and therefore are the optimal illumination states. If the composite retarder $\mathbf {U}_j^\dagger \mathbf {U}_i$ has a retarder vector $\vec {\boldsymbol {\delta }}_{j i}=[\delta _H,\delta _{45},\delta _R]$, then the great circle on the Poincaré sphere containing the optimal illumination states $\vec {\mathbf {g}}$ is defined as points satisfying the equation

$$\vec{\boldsymbol{\delta}}_{j i}^\dagger\vec{\mathbf{g}}= \delta_{H}\cos(2\theta)\cos(\eta)+\delta_{45}\sin(2\theta)\cos(\eta)+\delta_{R}\sin(\eta)=0,$$
where $\theta$ and $\eta$ are the coordinates on the Poincaré sphere corresponding to the major axis orientation and ellipticity of polarization states.

For each $\vec {\mathbf {g}}$ on this great circle, there are two antipodal analyzers which maximize $\Delta p$ as calculated in Eq. (4). These solutions can be found with the following geometric construction on the Poincaré sphere. For the two states found by transforming the chosen generator with the different MMs $\vec {\mathbf {s}}_i$ and $\vec {\mathbf {s}}_j$, calculate their unit bisector $\vec {\mathbf {b}}$

$$\vec{\mathbf{b}}=\frac{\vec{\mathbf{s}}_i+\vec{\mathbf{s}}_j}{\|\vec{\mathbf{s}}_i+\vec{\mathbf{s}}_j\|}$$
and a unit vector mutually perpendicular
$$\vec{\mathbf{c}} =\frac{\vec{\mathbf{s}}_i\times\vec{\mathbf{s}}_j}{\|\vec{\mathbf{s}}_i\times\vec{\mathbf{s}}_j\|}.$$

The analyzer states which maximize the measurement difference between $\mathbf {s}_i$ and $\mathbf {s}_j$ are then found using the cross product of the bisector and mutually perpendicular vector

$$\vec{\mathbf{a}}={\pm}\vec{\mathbf{b}}\times\vec{\mathbf{c}}.$$

A configuration with $-\mathbf {a}$ will have reversed contrast from a configuration with $\mathbf {a}$, i.e., bright regions in the image will become dark and vice versa. The set of optimal illumination states and associated analyzer states for comparing $\mathbf {M}_i$ and $\mathbf {M}_j$ found using Eqs. (6)–(9) are shown in Fig. 3(a). The set of analyzer states associated with the different illumination states trace out a different great circle. A polarimetric image performed using any optimal $\mathbf {g}$ with an optimal $\mathbf {a}$ will have the same $\Delta p$ between the two MM pixels used in the optimization. The family of optimal polariscopic pairs are indistinguishable for measurements comparing two pixels. However, the rest of the polariscopic image will generally be different for this family of polariscopic pairs.

 figure: Fig. 3.

Fig. 3. On the Poincaré sphere, the set of optimal generator illumination states are shown in (a) as a blue great circle perpendicular to the $\vec {\boldsymbol {\delta }}_{ji}$ axis shown with a blue arrow. The red great circle in (a) shows the set of optimal analyzers. The states constrained by the RGB950 hardware are shown in (b), with the set of generators shown on the blue figure-eight and the set of analyzers shown on the red figure-eight. In (c) there are three pairs: I, II, and III, where the optimal and hardware-constrained states are closest, in terms of a root sum of squared distance. In (a) and (b), the I, II, and III polariscopic pairs are shown with a circle, square, and triangle, respectively.

Download Full Size | PDF

3.2 Available illumination and analyzer states

The illumination optics of the RGB950 consist of LEDs, a diffuser, and a fixed linear polarizer which defines the orientation of horizontal polarization, and a rotating linear retarder [31]. The camera optics consist of the lens, bandpass filter, rotating retarder, fixed linear polarizer, and detector. The subject is illuminated with 0.0315 W/m$^2$ at 947 nm. The instrument is shown in Fig. 6 in the Appendix. The rotating retarders have variable orientation, but fixed retardance magnitude for a given wavelength, so the polarization states available for illumination and analysis are constrained.

For regions $i$ and $j$, the compound retarder from Eq. (5) has $\vec {\boldsymbol {\delta }}_{ji}=[-0.080, -3.108, 0.292]$. Optimal polarimetric pairs from Eq. (6) given $\vec {\boldsymbol {\delta }}_{ji}$ are compared to the states constrained by the RGB950 hardware is shown in Fig. 3. The set of polarization states available from rotating the retarders trace out figure-eights on the surface of the Poincaré sphere, see Fig. 3(b). The great circle of optimal illumination states intersects the figure-eight of available illumination states for all of the pairs of birefringent regions considered. However, the figure-eight of analyzer states constrained by the rotating retarder polarimeter will not necessarily intersect the optimal analyzer associated with the selected illumination state.

The spherical distance from the illumination great circle to the illumination figure-eight and the smaller of the spherical distances from the two associated analyzers to the analyzer figure-eight are used to select the polarimeter configurations. The root sum of squares of these distances were used as a figure of merit for how close an available illumination/analyzer pair was to optimal. This was performed numerically by discretizing the great circles into 360 equal segments and discretizing the figure-eights into $0.25^\circ$ step rotations of the retarders. The root sum of squares distance from each (arbitrarily parameterized) point on the great circles to the figure-eights is shown in Fig. 3(c). The three minima at positions 62, 184, and 307 (denoted I, II, and III) were selected for further analysis.

4. Results

The modulation pattern was calculated by applying the chosen polariscopic pairs (e.g. states $\mathbf {g}$ and $\mathbf {a}$) to the pure retarder approximation of $\widehat {\mathbf {M}}_0$ at each pixel according to Eq. (1). The expected modulation patterns for the three optimal pairs discussed in Sect. 3.2 are shown in Fig. 4(a). Region $i$ and region $j$ (see Fig. 2) appear as dark and bright, respectively, in the modulation patterns computed for all three polariscopic pairs. However, the pattern over the rest of the image varies slightly between the three pairs, with I and III having a more rounded central bright region and II having a more distinct "plus sign" bright region. The polariscopic images measured in the RGB950 are taken from Visualization 1, Visualization 2 and Visualization 3 and shown in Fig. 4(b). The rounded bright region expected in the center of images performed with polariscopic pairs I and III, and the bright "plus sign" in the image with pair II are visually apparent. The pattern in measurement III is less obvious, but this is consistent with pair III being further from optimal as compared to I and II. Pixel-wise quantitative comparison is not performed here due to the challenges of low repeatability in positioning. It should be noted that the three different polariscopic pairs shown in Fig. 4 are optimal for the same pair of regions ($i$ and $j$) and therefore are expected to only have minor differences in the modulation pattern over the rest of cornea. Polariscopic pairs optimized for different ROIs would produce different modulation patterns. For example, when $i$ and $j$ are used in the optimization, $i$ and $l$ have similar responses. If the contrast between $i$ and $l$ were optimized, these ROIs would necessarily have different responses.

 figure: Fig. 4.

Fig. 4. Simulated spatial modulation patterns (a) using three optimal illumination/analyzer pairs and measurements (b) with the nearest available pairs in the RGB950. The simulated results are calculated using Eq. (1) and the pure retarder component from the TD approximation of the MM images. The spatially-varying patterns shown in simulation (a) are apparent in the measured results (b), though the presence of depolarization reduced the magnitude of the modulation. For this reason, the dynamic range of the measured polariscopic images was clipped to emphasize the modulation. The images in (b) are still frames taken from Visualization 1, Visualization 2 and Visualization 3.

Download Full Size | PDF

The presence of depolarization reduced the magnitude of modulation for the polariscopic images in Fig. 4(b). The depolarization index (which ranges from 0 for the ideal depolarizer to 1 for a MJM) had an average value of 0.157 over the masked region of the MM measurement. The reduced modulation due to depolarization means that the dark regions of the pattern are not completely black in the measurement. Therefore, the dynamic range of the images in Fig. 4(b) was clipped in post-processing to make the polarimetric modulation pattern more visually apparent. This gives the appearance of overexposure in the sclera to the right of the cornea and on the skin of the lower eyelid. The bright spot to the right of the pupil is glare reflected from the first surface of the eye. During MM acquisition, this spot is overexposed so the reconstructed MM value is inaccurate. Additionally, this light did not transmit into the eye and through the cornea, so the polarimetric response would not be well-represented by the birefringent model used for the cornea.

Sources of disagreement between the polariscopic image and the modulation pattern computed from MM image are: the assumption of a TD eigenspectrum, the assumption of a pure retarder $\widehat {\mathbf {M}}_0$, the uniformity assumptions on $\xi _0$ and $M_{00}$ (discussed in the Appendix), motion artifacts in the MM measurement because it is reconstructed from 25 images acquired over 15 seconds, blurring in the MM intermediate images and the polariscopic image due to the exposure time, challenges with repeatability in the subject’s eye position between measurements, and measurement noise. Over the masked region of the MJM image, the diattenuation magnitude had an average value of 25.5% which is fairly significant though the polarimetric behavior is dominated by retardance. The differences between the dominant MJM as calculated from the original MM and the pure retarder can be seen Fig. 1(c) and Fig. 1(d).

The modulation patterns for polariscopic pair II were calculated for 19 other human subjects based on their eye MM images from Dataset 1 [16]. These patterns are shown in Fig. 5. The bright central "plus sign" appears for all individuals suggesting that the solution optimized for one human subject is robust to many more, though parts of the pattern may be obscured by the eyelid or pupil.

 figure: Fig. 5.

Fig. 5. Modulation patterns for polariscopic pair II simulated from the MM images of 20 different individuals’ eyes found in Dataset 1. Configuration II was optimized for the first image (top-left), but qualitatively similar modulation patterns appear for all 20 eyes.

Download Full Size | PDF

5. Conclusions

In this work we presented a closed-form solution for maximizing the contrast between two retarder MMs in a polariscope imaging design. The method is general for any pair of elliptical retarders, and can be easily extended to depolarizing MMs that are well-approximated using a TD model. The wide range of observed retardance orientations, from purely linear to purely circular, observed across an individual’s eye suggests that polariscopic imaging can be optimized to effectively identify different regions. To demonstrate the method, we chose two regions from the MM image of an in vivo human cornea and computed the family of optimal illumination and analyzer polarization states. From this family of solutions, we selected three polariscopic pairs that were as close as possible to states available in an existing polarimeter.

We compared the expected modulation pattern of these three optimal pairs to polariscopic using existing polarimeter hardware. The generator and analyzer states for configuration II are both approximately horizontal linear polarization, which is a measurement performed by Sobzcak et al [13]. The modulation pattern observed in our simulation and measurement of configuration II is similar to the patterns reported in their work. The consistency between their measurements and our simulation based on MMs with elliptical retardance suggests that the significant ellipticity of the birefringence human cornea has been overlooked in other works. Notably, elliptical retardance can be created from stacking linear retarders with unaligned fast axes. The stroma is organized in about 200–500 fibrous lamellae and the preferential orientation of collagen fibrils is known to be spatially-varying [38]. In the center of the cornea the lamellae orientations are expected to be one of two perpendicular directions which would not produce circular retardance. This is consistent with the MM measurements in the center of the cornea that appear in the literature [15]. Further from the center as the cornea curves to meet the sclera at the limbus, a greater degree of misalignment in the lamellae is expected which could produce elliptical and circular net birefringence [6]. Further analysis is required to support or refute a hypothesis that the elliptical and circular retardance we observed in MM human eye measurements are the result of fibril orientation.

The optimization is affected by measurement noise to the extent that the initial MM characterization is affected. For this reason, the optimization was performed based on retardance vectors averaged over ROIs. These optimized measurements are invariant to calibration errors in the retardance magnitude or rotational offset of the waveplates because the same instrument was used to perform the initial MM characterization as the subsequent polariscopic images. Even if the reported optimal states do not accurately describe the real illumination and analyzer polarization states, the contrast would still be maximized. However, if we perform the optimization based on MM data with calibration issues and then performed the polariscopic images from a different polarimeter, the calibration issues would manifest as reduced contrast.

It is worth noting that while our method maximizes the contrast between two retarder MMs, the method does not guarantee that the two chosen MMs will be the brightest or darkest in the polarimetric image of an object with spatially varying birefringence. For example, if two MMs are chosen which have similar retardance orientations and magnitudes, the maximum achievable contrast between those two MMs may be lower than the maximum contrast imaged for another pair of regions. Furthermore, the choice of ROIs in the image also affects the specific spatial modulation patterns created by polariscopic imaging. This is an additional design variable as some patterns may be preferable depending on the application.

The use of an existing dual-rotating-retarder polarimeter constrained the polariscopic pairs that could be measured. If a polarimeter is designed for polariscopic eye imaging the illumination and analyzer states can be optimized for this task. Consequently, this polariscopic demonstration does not achieve an upper contrast limit. Increased modulation and different spatial distributions are potentially possible if a near-infrared polariscopic imaging system is designed for in vivo human eyes. Division-of-focal-plane polarization cameras, which acquire images through four different linear polarizers simultaneously at the cost of spatial resolution, have been used in the literature to improve polarimetric measurement speed [39,40]. Simultaneous image capture for the different analyzers is of particular interest for measurements of eyes. However, this would be a fundamentally different optimization problem than is explored in this work.

6. Appendix

MM images and polariscopic images were taken using a dual-rotating-retarder Mueller polarimeter called the RGB950, shown in Fig. 6 [31]. The light source consists of LEDs behind a diffuser, linear polarizer, and rotating retarder. The human subjects are illuminated with a diffuse near-infrared LED light source (945-955 nm) where 0.0315 W/m$^2$ is incident on the subject’s face and/or eye.

 figure: Fig. 6.

Fig. 6. Mueller image measurement of in vivo human eyes being performed in (a) using the RGB950 polarimeter at 947 nm with instrument schematic in (b). The eyes were measured in a reflective double-pass configuration where light transmitted through the cornea, reflected off the iris, and transmitted back through the cornea a second time. The duration of a MM image acquisition is 15 seconds during which 25 unique polarimetric images are captured. When operating at 947 nm, the illumination retarder has a retardance magnitude of 118$^\circ$ and the analyzer retarder has a magnitude of 133$^\circ$.

Download Full Size | PDF

In a TD MM, the three smallest eigenvalues of the coherency matrix (which is linearly related to the MM) are degenerate [32]. If all four eigenvalues are normalized to sum to unity, then taking the smaller three eigenvalues $[\xi _1,\xi _2,\xi _3]$ as 3D rectilinear coordinates forms the natural depolarization space [41]. In this space, the Euclidean distance from $[\xi _1,\xi _2,\xi _3]$ to the TD point $\frac {1-\xi _0}{3}[1,1,1]$ is

$$\Delta \xi_{TD} =\sqrt{6\left[\left(\xi_1-\frac{1-\xi_0}{3}\right)^2+\left(\xi_2-\frac{1-\xi_0}{3}\right)^2+\left(\xi_3-\frac{1-\xi_0}{3}\right)^2\right]},$$
where the factor of $\sqrt {6}$ is included to normalize the TD distance to range from 0 to 1, with 0 being exactly TD and 1 being as far from TD as possible. Figure 7(a) shows $\Delta \xi _{TD}$ calculated over the MM image of the cornea and Fig. 7(d) shows the same data summarized in a histogram. $\Delta \xi _{TD}$, which has a mean value of 0.151 and a standard deviation of 0.029, is assumed to be low enough for the MM image to be accurately approximated as TD.

 figure: Fig. 7.

Fig. 7. Images (a-c) and histograms (d) of quantities derived from the MM image of an in vivo human cornea. $\Delta \xi _{TD}$, calculated according to Eq. (10), (a) describes how close the coherency matrix eigenspectrum is to being TD. $\Delta \xi _{TD}$ is assumed to be low enough for the MM image to be approximated as TD. The mean and standard deviation of the TD distance are 0.151 and 0.029, respectively. $\xi _0$ (b) is the parameter which controls the amount of depolarization in a TD MM. The mean and standard deviation of $\xi _0$ are 0.355 and 0.023, respectively. $M_{00}$ (c) is the overall reflectance, here normalized by the brightest pixel in the image. The mean and standard deviation of $M_{00}$ are 0.485 and 0.042, respectively. $\xi _0$ and $M_{00}$ are assumed to be uniform over the cornea such that they do not contribute to differences in polarimetric measurement across the image.

Download Full Size | PDF

For any combination of generators and analyzers, the polarimetric measurement of the ideal depolarizer is constant $\mathbf {a}\mathbf {M}_{ID}\mathbf {g}=\frac {1}{2}$. Modulation in a polarimetric measurement of a TD MM is therefore only the result of $\mathbf {a}\widehat {\mathbf {M}}_0\mathbf {g}$, with $M_{00}$ and $\xi _0$ determining the overall brightness and the constant offset, respectively. The average reflectance $M_{00}$ and depolarization parameter $\xi _0$ are shown as images in Fig. 7 b and c, respectively, and summarized as histograms in Fig. 7(d). The standard deviation of $\xi _0$ is 0.023 and the standard deviation of $M_{00}$ is 0.042, so these parameters are assumed to be uniform such that their influence on a polarimetric measurement $p$ is the same over the cornea.

Acknowledgments

The authors would like to thank Adeline Tai for performing the Mueller imaging experiments. This work was made possible thanks to the 20 human subject participants.

Disclosures

The authors declare no conflicts of interest.

Data availability

The 20 MM image measurements discussed in this paper are available in Dataset 1, Ref. [16].

References

1. D. Brewster, “III. Experiments on the depolarisation of light as exhibited by various mineral, animal, and vegetable bodies, with a reference of the phenomena to the general principles of polarisation,” Philos. Trans. R. Soc. London 105, 29–53 (1815). [CrossRef]  

2. S. J. Mccafferty, J. T. Schwiegerling, T. L. Koch, et al., Analysis and Application of Opto-Mechanics to the Etiology of Sub-Optimal Outcomes in Laser Corrective Eye Surgery and Design Methodology of Deformable Surface Accommodating Intraocular Lenses (The University of Arizona., 2015).

3. J. A. Germann, E. Martinez-Enriquez, and S. Marcos, “Quantization of collagen organization in the stroma with a new order coefficient,” Biomed. Opt. Express 9(1), 173–189 (2018). [CrossRef]  

4. J. A. Germann, E. Martínez-Enríquez, M. C. Martínez-García, et al., “Corneal Collagen Ordering After In Vivo Rose Bengal and Riboflavin Cross-Linking,” Invest. Ophthalmol. Vis. Sci. 61(3), 28 (2020). [CrossRef]  

5. A. Stanworth and E. J. Naylor, “The polarization optics of the isolated cornea,” Br. J. Ophthalmol. 34(4), 201–211 (1950). [CrossRef]  

6. R. H. Newton and K. M. Meek, “The integration of the corneal and limbal fibrils in the human eye,” Biophys. J. 75(5), 2508–2512 (1998). [CrossRef]  

7. M. Pircher, E. Goetzinger, R. Leitgeb, et al., “Transversal phase resolved polarization sensitive optical coherence tomography,” Phys. Med. Biol. 49(7), 1257–1263 (2004). [CrossRef]  

8. F. Fanjul-Vélez, M. Pircher, B. Baumann, et al., “Polarimetric analysis of the human cornea measured by polarization-sensitive optical coherence tomography,” J. Biomed. Opt. 15(5), 056004 (2010). [CrossRef]  

9. M. Pircher, C. K. Hitzenberger, and U. Schmidt-Erfurth, “Polarization sensitive optical coherence tomography in the human eye,” Prog. Retinal Eye Res. 30(6), 431–451 (2011). [CrossRef]  

10. J. C. Ramella-Roman, I. Saytashev, and M. Piccini, “A review of polarization-based imaging technologies for clinical and preclinical applications,” J. Opt. 22(12), 123001 (2020). [CrossRef]  

11. M. Sobczak, M. Asejczyk-Widlicka, A. Szafraniec, et al., “Analysis of torsional eye movements using the corneal birefringence pattern,” J. Opt. Soc. Am. A 36(4), B23 (2019). [CrossRef]  

12. M. Sobczak, M. Owczarek, W. A. Wozniak, et al., “In vivo measurements of corneal birefringence properties using the one-way reflective Mueller polarimetry,” Opt. Express 29(10), 15356 (2021). [CrossRef]  

13. M. Sobczak and M. Asejczyk, “Birefringent properties of the cornea measured by a Mueller type polarimeter in healthy adults and children,” Biomed. Opt. Express 12(12), 7872 (2021). [CrossRef]  

14. M. Sobczak, M. Asejczyk, and M. Wilczynski, “The effect of pupil size on the measurement of corneal birefringence properties: preliminary study,” Sci. Rep. 13(1), 17439 (2023). [CrossRef]  

15. J. M. Bueno, “Measurement of parameters of polarization in the living human eye using imaging polarimetry,” Vision Res. 40(28), 3791–3799 (2000). [CrossRef]  

16. A. Tai, Q. Jarecki, and M. Kupinski, “Near-infrared human eye Mueller matrix images,” University of Arizona, Accessed: Dec. 4, 2023https://doi.org/10.25422/azu.data.24722358.

17. T. Novikova and J. C. Ramella-Roman, “Is a complete Mueller matrix necessary in biomedical imaging?” Opt. Lett. 47(21), 5549–5552 (2022). [CrossRef]  

18. M. Kupinski, M. Boffety, F. Goudail, et al., “Polarimetric measurement utility for pre-cancer detection from uterine cervix specimens,” Biomed. Opt. Express 9(11), 5691–5702 (2018). [CrossRef]  

19. M. Kupinski, J. Rehbinder, H. Haddad, et al., “Tasked-based quantification of measurement utility for ex vivo multi-spectral Mueller polarimetry of the uterine cervix,” in Clinical and Preclinical Optical Diagnostics, vol. 10411J. Q. BrownT. G. van Leeuwen, eds., International Society for Optics and Photonics (SPIE, 2017), p. 104110N.

20. R. Chipman, W. Lam, and G. Young, Polarized Light and Optical Systems, Optical Sciences and Applications of Light (CRC Press, 2018).

21. F. Goudail and A. Bénière, “Optimization of the contrast in polarimetric scalar images,” Opt. Lett. 34(9), 1471 (2009). [CrossRef]  

22. M. K. Kupinski, J. Bankhead, A. Stohn, et al., “Binary classification of Mueller matrix images from an optimization of Poincaré coordinates,” J. Opt. Soc. Am. A 34(6), 983 (2017). [CrossRef]  

23. D. Upadhyay, S. Mondal, E. Lacot, et al., “Full analytical solution of adapted polarisation state contrast imaging,” Opt. Express 19(25), 25188 (2011). [CrossRef]  

24. L. Di Cecilia, F. Marazzi, and L. Rovati, “Spectral repeatability of a hyperspectral system for human iris imaging,” in 2018 IEEE 4th International Forum on Research and Technology for Society and Industry (RTSI), (2018), pp. 1–5.

25. A. D. Martino, Y.-K. Kim, E. Garcia-Caurel, et al., “Optimized mueller polarimeter with liquid crystals,” Opt. Lett. 28(8), 616–618 (2003). [CrossRef]  

26. J. Vargas, N. Uribe-Patarroyo, J. A. Quiroga, et al., “Optical inspection of liquid crystal variable retarder inhomogeneities,” Appl. Opt. 49(4), 568–574 (2010). [CrossRef]  

27. N. C. Bruce, O. G. Rodríguez-Herrera, J. M. López-Téllez, et al., “Experimental limits for eigenvalue calibration in liquid-crystal mueller-matrix polarimeters,” Opt. Lett. 43(11), 2712–2715 (2018). ObjectType-Article-1. [CrossRef]  

28. A. Álvarez Herrero, P. G. Parejo, and M. Silva-López, “Fine tuning method for optimization of liquid crystal based polarimeters,” Opt. Express 26(9), 12038–12048 (2018). [CrossRef]  

29. I. Montes-González, O. G. Rodríguez-Herrera, M. A. no alejo, et al., “Effects of typical liquid-crystal retarder errors on optimized stokes polarimeters,” Appl. Opt. 61(35), 10458–10464 (2022). [CrossRef]  

30. I. Montes-Gonzalez, M. A. no Alejo, N. C. Bruce, et al., “Optimized configuration for liquid crystal stokes polarimeters in the presence of fast-axis orientation errors,” Opt. Lett. 49(2), 355–358 (2024). [CrossRef]  

31. J. M. López-Téllez, R. A. Chipman, L. W. Li, et al., “Broadband extended source imaging Mueller-matrix polarimeter,” Opt. Lett. 44(7), 1544–1547 (2019). [CrossRef]  

32. L. Li and M. Kupinski, “Merit functions and measurement schemes for single parameter depolarization models,” Opt. Express 29(12), 18382–18407 (2021). [CrossRef]  

33. R. Ossikovski and O. Arteaga, “Integral decomposition and polarization properties of depolarizing Mueller matrices,” Opt. Lett. 40(6), 954–957 (2015). [CrossRef]  

34. S. R. Cloude and E. Pottier, “Concept of polarization entropy in optical scattering,” Opt. Eng. 34(6), 1599–1610 (1995). [CrossRef]  

35. S. R. Cloude, “Group theory and polarisation algebra,” Optik 75, 26–36 (1986).

36. Q. Jarecki and M. K. Kupinski, “Underdetermined polarimetric measurements for Mueller extrapolations,” Opt. Eng. 61(12), 123104 (2022). [CrossRef]  

37. Q. Jarecki and M. Kupinski, “New polarized representations for depolarization-dominant materials,” Opt. Express (in review).

38. A. Pandolfi and F. Manganiello, “A model for the human cornea: Constitutive formulation and numerical analysis,” Biomech. Model. Mechanobiol. 5(4), 237–246 (2006). [CrossRef]  

39. T. Huang, R. Meng, J. Qi, et al., “Fast mueller matrix microscope based on dual dofp polarimeters,” Opt. Lett. 46(7), 1676–1679 (2021). [CrossRef]  

40. X. Li, F. Goudail, and S.-C. Chen, “Self-calibration for mueller polarimeters based on dofp polarization imagers,” Opt. Lett. 47(6), 1415–1418 (2022). [CrossRef]  

41. R. Ossikovski and J. Vizet, “Eigenvalue-based depolarization metric spaces for Mueller matrices,” J. Opt. Soc. Am. A 36(7), 1173–1186 (2019). [CrossRef]  

Supplementary Material (4)

NameDescription
Dataset 1       Dataset 1, Mueller matrix images of 20 human subjects' eyes.
Visualization 1       Video capture of a human eye under elliptical polarized illumination with another elliptical polarizer in front of the camera. The illumination and analyzer polarizers correspond to pair I in the manuscript.
Visualization 2       Video capture of a human eye under elliptical polarized illumination with another elliptical polarizer in front of the camera. The illumination and analyzer polarizers correspond to pair II in the manuscript.
Visualization 3       Video capture of a human eye under elliptical polarized illumination with another elliptical polarizer in front of the camera. The illumination and analyzer polarizers correspond to pair III in the manuscript.

Data availability

The 20 MM image measurements discussed in this paper are available in Dataset 1, Ref. [16].

16. A. Tai, Q. Jarecki, and M. Kupinski, “Near-infrared human eye Mueller matrix images,” University of Arizona, Accessed: Dec. 4, 2023https://doi.org/10.25422/azu.data.24722358.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1.
Fig. 1. Near-infrared MM measurements of an in vivo human eye, masked to the region of the cornea where light is reflected back from the iris. In (a), the image is normalized such that $M_{00}=1$ to visualize the spatially-varying polarization properties. In (b), the MM measurement is approximated as TD using Eq. (2) with (c) being the dominant MJM $\widehat {\mathbf {M}}_0$. In (d), the dominant MJM is further approximated as a pure retarder MJM.
Fig. 2.
Fig. 2. The horizontal, diagonal, and circular retardance computed from a measured MM are shown in (a), (b), and (c), respectively. The total retardance magnitude is shown in (d). ROIs are designated $i$-$l$. In (e), the retardance vectors within each ROI are plotted in 3D, where the concentric spherical surfaces denote $\pi /2$, $\pi$, and $3\pi /2$ radians of retardance. The point clouds represent the per-pixel retardance vectors and the arrows indicate the average retardance vector for each ROI. Average values are listed in Table 1.
Fig. 3.
Fig. 3. On the Poincaré sphere, the set of optimal generator illumination states are shown in (a) as a blue great circle perpendicular to the $\vec {\boldsymbol {\delta }}_{ji}$ axis shown with a blue arrow. The red great circle in (a) shows the set of optimal analyzers. The states constrained by the RGB950 hardware are shown in (b), with the set of generators shown on the blue figure-eight and the set of analyzers shown on the red figure-eight. In (c) there are three pairs: I, II, and III, where the optimal and hardware-constrained states are closest, in terms of a root sum of squared distance. In (a) and (b), the I, II, and III polariscopic pairs are shown with a circle, square, and triangle, respectively.
Fig. 4.
Fig. 4. Simulated spatial modulation patterns (a) using three optimal illumination/analyzer pairs and measurements (b) with the nearest available pairs in the RGB950. The simulated results are calculated using Eq. (1) and the pure retarder component from the TD approximation of the MM images. The spatially-varying patterns shown in simulation (a) are apparent in the measured results (b), though the presence of depolarization reduced the magnitude of the modulation. For this reason, the dynamic range of the measured polariscopic images was clipped to emphasize the modulation. The images in (b) are still frames taken from Visualization 1, Visualization 2 and Visualization 3.
Fig. 5.
Fig. 5. Modulation patterns for polariscopic pair II simulated from the MM images of 20 different individuals’ eyes found in Dataset 1. Configuration II was optimized for the first image (top-left), but qualitatively similar modulation patterns appear for all 20 eyes.
Fig. 6.
Fig. 6. Mueller image measurement of in vivo human eyes being performed in (a) using the RGB950 polarimeter at 947 nm with instrument schematic in (b). The eyes were measured in a reflective double-pass configuration where light transmitted through the cornea, reflected off the iris, and transmitted back through the cornea a second time. The duration of a MM image acquisition is 15 seconds during which 25 unique polarimetric images are captured. When operating at 947 nm, the illumination retarder has a retardance magnitude of 118$^\circ$ and the analyzer retarder has a magnitude of 133$^\circ$.
Fig. 7.
Fig. 7. Images (a-c) and histograms (d) of quantities derived from the MM image of an in vivo human cornea. $\Delta \xi _{TD}$, calculated according to Eq. (10), (a) describes how close the coherency matrix eigenspectrum is to being TD. $\Delta \xi _{TD}$ is assumed to be low enough for the MM image to be approximated as TD. The mean and standard deviation of the TD distance are 0.151 and 0.029, respectively. $\xi _0$ (b) is the parameter which controls the amount of depolarization in a TD MM. The mean and standard deviation of $\xi _0$ are 0.355 and 0.023, respectively. $M_{00}$ (c) is the overall reflectance, here normalized by the brightest pixel in the image. The mean and standard deviation of $M_{00}$ are 0.485 and 0.042, respectively. $\xi _0$ and $M_{00}$ are assumed to be uniform over the cornea such that they do not contribute to differences in polarimetric measurement across the image.

Tables (1)

Tables Icon

Table 1. Components and magnitudes of the average retarder vectors from the four ROIs shown in Fig. 2 in units of radians. Retardance magnitudes are known up to the order of retardance, where 2 π q for integers q can be added to the magnitude without changing the resulting MM.

Equations (10)

Equations on this page are rendered with MathJax. Learn more.

p = a M g
M 4 M 00 3 [ ( ξ 0 1 4 ) M ^ 0 + ( 1 ξ 0 ) M I D ] ,
M E R ( δ ) = [ 1 0 0 0 0 δ H 2 + ( δ 45 2 + δ R 2 ) C δ 2 δ 45 δ H T δ 2 + δ R S δ δ H δ R T δ 2 δ 45 S δ 0 δ 45 δ H T δ 2 δ R S δ δ 45 2 + ( δ R 2 + δ H 2 ) C δ 2 δ R δ 45 T δ 2 + δ H S δ 0 δ H δ R T δ 2 + δ 45 S δ δ R δ 45 T δ 2 δ H S δ δ R 2 + ( δ 45 2 + δ H 2 ) C δ 2 ] = [ 1 0 T 0 U ( δ ) ] ,
Δ p = | | a ( M i M j ) g | | ,
s j s i = ( U j g ) ( U i g ) = g ( U j U i g )
δ j i g = δ H cos ( 2 θ ) cos ( η ) + δ 45 sin ( 2 θ ) cos ( η ) + δ R sin ( η ) = 0 ,
b = s i + s j s i + s j
c = s i × s j s i × s j .
a = ± b × c .
Δ ξ T D = 6 [ ( ξ 1 1 ξ 0 3 ) 2 + ( ξ 2 1 ξ 0 3 ) 2 + ( ξ 3 1 ξ 0 3 ) 2 ] ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.