*in vivo*. Six healthy subjects were recruited for noninvasive imaging using a snapshot hyperspectral system. The three-dimensional full spatial-spectral data cube was analyzed using non-negative matrix factorization (NMF), wherein the data was decomposed to give spectral signatures and spatial distribution, in search for the MP absorbance spectrum. The NMF was initialized with the

*in vitro*MP spectrum and rank 4 spectral signature decomposition was used to recover the MP spectrum and optical density in vivo. The recovered MP spectra showed two peaks in the blue spectrum, characteristic of MP, giving a detailed

*in vivo*demonstration of these absorbance peaks. The peak MP optical densities ranged from 0.08 to 0.22 (mean 0.15+/−0.05) and became spatially negligible at diameters 1100 to 1760 μm (4 to 6 deg) in the normal subjects. This objective method was able to exploit prior knowledge (the

*in vitro*MP spectrum) in order to extract an accurate

*in vivo*spectral analysis and full MP spatial profile, while separating the MP spectra from other ocular absorbers. Snapshot hyperspectral imaging in combination with advanced mathematical analysis provides a simple cost-effective approach for MP mapping

*in vivo*.

## 1.

## Introduction

The role of macular pigment (MP) in ocular health and disease remains controversial. It has been suggested as a scavenger for a reactive oxygen species and a blue light filter, and, hence, to have a protective role against light-induced damage to the photoreceptor/RPE complex.^{1, 2} Lutein (L) and zeaxanthin (Z), the carotenoids that constitute MP,^{3} have antioxidant properties that appear to be protective against blue light.^{4} The optical filtering properties of MP, to reduce chromatic aberration and glare, have also been described.^{1} Its role in age-related macular degeneration (AMD) remains controversial. While some have proposed a preventative role,^{5, 6} others have disputed this role in their populations of AMD.^{7, 8} Furthermore, the ability to modify the density of macular pigment through dietary supplements remains an important, yet unanswered question. Most studies have shown a significant relationship and this has been comprehensively reviewed.^{9, 10}

MP measures are subject to large amounts of variability in normal subjects.^{11, 12, 13} An obstacle to clearly establishing the role of MP in retinal disease is the lack of a validated, clinically applicable, and widely accessible tool to quantify the MP *in vivo.* Noninvasive MP assessment methods include psychophysical techniques and objective measures. Psychophysical methods such as the matching tasks of heterochromic flicker photometry (HFP) (Refs.
12, 14, 15, 16) or minimum motion photometry^{17, 18} assume normal retinal function and uniform lens density within the area of measurement and rely on the accuracy of patient responses. Patients with advanced ocular disease tend to experience the greatest difficulty with such tests.^{15} Commercially available psychophysical measures of MP generate a single value, rather than a complete profile of measures for each eccentricity. Objective methods include reflectance imaging,^{19, 20, 21} autofluorescence (AF) imaging,^{22, 23, 24} and Raman spectroscopy.^{25, 26} Reflectance and AF imaging give the spatial distribution of MP, but in some cases require accurate fixation and require prebleaching to avoid confounding absorption by photopigments, which involves unpleasant light levels. Reflectance measures are affected by stray light, unless confocal imaging is used. The AF technique assumes that the relative spectral energy of lipofuscin fluorescence is constant across the central retina.^{3} Despite the high chemical specificity of Raman spectroscopy, its signal can be attenuated by lenticular absorbance or scattering and maximal pupil dilation is required for measurement.^{3, 27}

Hence, there is a need for a rapid, objective, and simple method to noninvasively evaluate the MP. In our previous research,^{28, 29} we evaluated a hyperspectral reflectometry camera, which captures a 20 deg field with 76 bands of spectral information, in conjunction with a partially constrained unsupervised data mining approach (blind source separation) to demonstrate the spectra that colocalized with drusen. Here, we applied this technique to quantify the MP in a group of healthy eyes *in vivo*. Our results are comparable to measurements made with other reflectance approaches. There has been recent interest in hyperspectral imaging using a modified fundus camera and the various systems have been briefly reviewed.^{30} Prior multispectral retinal image analysis has been restricted, either in spatial extent to single line scans^{31} or to pointwise acquisition,^{32} or spectrally to a limited number of spectral bands.^{33, 34, 35} While excellent focal analysis of MP has been possible with two wavelength reflectance^{36, 37} or AF,^{36} we wanted to investigate the full 20 deg field with 76 bands of spectral information and explore the use of blind source separation to obtain detailed spectral characterization of MP *in vivo*.

## 2.

## Methods

## 2.1.

### Subjects

We recruited six healthy volunteers (mean age 32.8, ±12.1, range 22 to 55 years, three males and three females; Fig. 1) from the University of Southern California. Each patient underwent a complete eye exam followed by pupil dilation and imaging with the hyperspectral camera. Photoreceptor pigment prebleaching was not systematically attempted.

This study was approved by the institutional review boards of the University of Southern California and Columbia University and adhered to the tenets of the Declaration of Helsinki for research involving human subjects. Each participant gave informed consent after explanation of the nature and possible consequences of the study.

## 2.2.

### Hyperspectral Camera

To avoid image coregistration and motion artifacts in spectral space, this camera employed a novel technique, computed tomographic imaging spectrometer (CTIS), for hyperspectral retinal imaging. The hyperspectral imaging device consists of an imaging spectrometer attached to a commercial fundus camera. This device captures a full spatial-spectral image cube (20 deg field, 186 × 186 pixel spatial resolution) in a single camera flash. CTIS uses diffractive optics to multiplex spatial-spectral information onto a focal plane array. With less than 20 msec integration time, the retinal CTIS system acquired 76 contiguous spectral bands between 420 and 720 nm at 4 nm resolution^{38, 39} (Fig. 2), thus obviating the need for subject immobilization for precise image registration.^{40}

## 2.3.

### Image Reconstruction

The images were reconstructed, using the previously described approach.^{39} Briefly, expectation maximization algorithms are employed on the diffraction pattern to reconstruct the image cube. The reconstructed cubes were then anonymously transferred to Columbia University for further data analysis.

## 2.4.

### Data Analysis

We employed a data mining approach to explore and search for constituent pigment spectra of the retina. In particular, we focused on extracting the spatio-spectral MP signatures in healthy subjects. An unsupervised learning algorithm called non-negative matrix factorization (NMF) (Refs. 41 and 42) was employed with a partial constraint in the form of a physiologic meaningful initialization scheme. NMF as usually implemented decomposes the data as a linear sum of sources times abundances; specifically, NMF decomposes a multivariate dataset [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document} $\mathbf{\text{X}}$ into two matrices: a matrix of spectral signatures [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries S}}$\end{document} $\mathbf{\text{S}}$ and their corresponding spatial distribution [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries A}}$\end{document} $\mathbf{\text{A}}$ . NMF assumes the linear model

^{43, 44}Unlike principal component analysis, NMF imposes positivity constraints on the spectra and spatial abundance images, consistent with physical reality. Such constraints yield a parts-based representation of the hyperspectral cube that allows an individual signature to be obtained, i.e., its spatio-spectral characterization, which has a physically meaningful interpretation.

We used NMF in a modified format on the hyperspectral reflectance data cube to recover and investigate absorbance signatures. To do this, the spectral reflectance *R* can be written with the following simplified model (adapted from Refs. 20 and 45):

*R*

_{SC}is the reflectance of the sclera and

*D*is the total optical density (OD) of all absorbers. To convert from reflectance to absorbance the data is log-transformed, assuming

*D*to be

*D*= OD

_{MP}+ OD

_{SubT}, where MP denotes macular pigment and SubT is the subtotal of all optical densities of all absorbers except MP. For log

_{10}of the data for each wavelength, we get

_{MP}, as one of the dominant spectral sources, then at every pixel and for every wavelength λ we have the decomposition

_{MP}is the recovered spectrum and [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries A}}$\end{document} $\mathbf{\text{A}}$

_{MP}is its abundance. It is customary to report the optical density at the absorbance peak, λ = 458 nm, which is the target density mapped. In this study, we reported the OD

_{MP}at λ = 456 nm due to spectral resolution of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document} $\mathbf{\text{X}}$ (i.e., 420 to 720 nm in 4 nm steps). In our case, we alternatively performed NMF on the data model

_{MP}may be constrained and normalized, say with value = 1 at λ = 456 nm, it does not affect the recovery of OD

_{MP}because the abundance [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries A}}$\end{document} $\mathbf{\text{A}}$

_{MP}is automatically scaled so that the product gives the actual data. We note, that besides reporting OD

_{MP}at 456 nm, the MP signature and abundance that NMF returns has a much richer spatio-spectral representation than is possible with other methods.

We performed the following normalization and initialization scheme. First, we divided each spectral band of the cube by its standard deviation to account for reflectance drifts along the spectrum. After normalization, we log-transformed the cube as just described and enforced positivity constraints on [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document} $\mathbf{\text{X}}$ to run NMF on the post-processed data. An example of [TeX:] \documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document} $\mathbf{\text{X}}$ can be seen in Fig. 3.

The number of latent sources in the hyperspectral cube, i.e., the optimum rank of the matrix
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries S}}$\end{document}
$\mathbf{\text{S}}$
of spectral signatures for recovery of the MP spectrum, was experimentally determined to be 4. For other ranks (e.g., 3 or 5) the decomposition of the macular pigment spectra showed redundancies, i.e., similar spectral profiles. In each case, we initialized one source as
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries S}}$\end{document}
$\mathbf{\text{S}}$
_{MP}, the *in vitro* absorbance prior (Fig. 4), where all other sources and abundance images were randomly initialized. However,
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries S}}$\end{document}
$\mathbf{\text{S}}$
_{MP} was allowed to vary during the NMF algorithm to obtain the best fit and most realistic spectrum. The obtained spatio-spectral signatures of the MP were assessed and analyzed by two retinal specialists (RTS, AF).

## 3.

## Results

## 3.1.

### Dataset Description

Figure 1 shows color fundus images of six normal eyes (N1 to N6). The resolution of the color fundus RGB images was 1052 × 914. Figure 3 shows the post-processed hyperspectral cubes of a healthy subject (N6).

## 3.2.

### Non-Negative Matrix Factorization Approach

Table 1 summarizes the employed configuration of NMF and reports the number of iterations used, the reconstruction error in terms of the root-mean-squared (RMS) error, and other convergence specific results. The RMS error is defined as the square root of the mean difference between the data
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document}
$\mathbf{\text{X}}$
and the factor model
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries AS}}$\end{document}
$\mathbf{\text{AS}}$
. We used a rank *r* = 4 decomposition to search for the known MP absorbance spectrum. Thus, if
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document}
$\mathbf{\text{X}}$
is the hyperspectral data set, NMF found a decomposition
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries AS}}$\end{document}
$\mathbf{\text{AS}}$
as an approximation to
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document}
$\mathbf{\text{X}}$
, where
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries S}}$\end{document}
$\mathbf{\text{S}}$
is the matrix of rank 4 of spectral signatures and
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries A}}$\end{document}
$\mathbf{\text{A}}$
is the matrix of abundances. The factors
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries A}}$\end{document}
$\mathbf{\text{A}}$
and
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries S}}$\end{document}
$\mathbf{\text{S}}$
are chosen to minimize the RMS between
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries X}}$\end{document}
$\mathbf{\text{X}}$
and
[TeX:]
\documentclass[12pt]{minimal}\begin{document}$\hbox{\sf{\bfseries AS}}$\end{document}
$\mathbf{\text{AS}}$
, as previously proposed.^{44} For all datasets in these experiments, 1000 iterations of the algorithm resulted in a reconstruction error range of approximately 1 to 4% at each pixel.

## Table 1

Reconstruction error for healthy subjects (N1 to N6). ID denotes an identifier key for a hyperspectral cube or its corresponding color fundus image, N is the number of iterations of the NMF algorithm, RMS is the root-mean-square residual, and D is the maximal change of either factor A or S.

Normals | |||
---|---|---|---|

ID | N | RMS | D |

N1 | 872 | 0.021399 | 0.00009994 |

N2 | 1000 | 0.024398 | 0.00037733 |

N3 | 1000 | 0.025532 | 0.00202895 |

N4 | 1000 | 0.018980 | 0.00013366 |

N5 | 1000 | 0.019192 | 0.00095528 |

N6 | 1000 | 0.021270 | 0.00068574 |

## 3.3.

### Recovered Macular Pigment Spectra and Optical Densities

Figure 5 shows the recovered MP spectra. In all six eyes, the recovered spectra show the two nearby peaks in the blue characteristic of both L and Z, the two carotenoids that constitute MP.^{3} We also found that fundus pigmentation did not have an impact on the performance of this approach, with OD_{MP} recovered in blond (N3) as well as highly pigmented fundi (N1, N2, and N6). The corresponding OD_{MP} maps are shown in Fig. 6. OD_{MP} values ranged from 0.08 to 0.22 (mean 0.15 ± 0.05, Table 2). Figure 7 shows a larger view of the OD_{MP} map for subject N1 and corresponding MP density profile plot along the horizontal meridian.

## Table 2

Macular pigment optical density (ODMP) for healthy subjects (normals). ID denotes an identifier key for a hyperspectral cube or its corresponding color fundus image, Max is the peak optical density value within a two degree annuli centered at the fovea, Mean is the average optical density within the same region as Max, and SD is its standard deviation.

Normals | |||
---|---|---|---|

ID | Max | Mean | SD |

N1 | 0.14504 | 0.14243 | 0.00531 |

N2 | 0.22397 | 0.22117 | 0.00987 |

N3 | 0.10200 | 0.07576 | 0.01001 |

N4 | 0.16015 | 0.13787 | 0.01581 |

N5 | 0.18881 | 0.16846 | 0.01847 |

N6 | 0.17816 | 0.17613 | 0.00531 |

## 4.

## Discussion

The hyperspectral tool we have utilized herein is an integration of highly novel imaging technology and the blind source separation method of NMF. The high spectral resolution (4 nm) of the CTIS multiplexing device coupled with the ability of NMF to precisely detect a spectral signal and its abundance, have together achieved a detailed demonstration of the two absorbance peaks of MP in normal eyes *in vivo*. To our knowledge, this is the first demonstration of the bifid absorbance spectrum of MP *in vivo*. Other researchers have attempted to fit their recovered spectral signatures to known MP spectra,^{34} which is rather different to our approach. The flexibility allowed by having multiple spectral points in the region of interest (21 measurements at each pixel in the region from 420 to 500 nm), and the use of MP signature only to initialize NMF search, give this approach the distinct advantage of reaching MP spectral signature based on individual eye's spatial/spectral data encoded in the hyperspectral cube. The accuracy of this method is demonstrated by the highly resolved recovered spectral signal in normal subjects. Moreover, the method is versatile in its ability to recover the spectra of known or unknown reflectors and absorbers. While in this project we used the known *in vitro* spectrum of MP to guide the recovery of its abundance, we have previously utilized the same blind source separation in a completely unsupervised fashion to extract spectral signatures that consistently co-localized with drusen.^{28, 29} Suggesting this approach can be used to achieve important tasks in the study of the normal and diseased fundus.

Advantages of this approach include its relative freedom from dependence on visual acuity or fixation, which may be highly valuable when evaluating patients with unstable and/or eccentric fixation. This is one of the main limitations of psychophysical approaches to MP assessments, such as HFP and customized HFP. Another limitation of the psychophysical approach relates to its basic assumption that OD_{MP} is zero outside the central 5 deg, which leads to its systematic underestimation of OD_{MP} when compared with other approaches.^{36, 46, 47, 48} It was previously shown that MP can still be present at low levels at 1 mm eccentricity (≈4 deg) in primate eyes^{49} and out to 8 deg in human donor eyes.^{50} A residual signal in our data at more eccentric locations might be attributed to this low level. A further advantage of the snapshot hyperspectral method is the capability to map a complete spatial profile of MP. Other hyperspectral methods, although they may use careful fits to models of fundus reflectance, only acquire MP density at a single point, for example, a measure of the averaged density in the central 1.9 deg.^{51} The snapshot method acquires a full three-dimensional spectral-spatial cube, which is simultaneously analyzed to measure MP density at every point of a 4 mm square field at a resolution of 22 μm.

The limitations of this approach include the inherent restrictions of measuring MP with reflectance. Delori have provided a detailed analysis showing that while MP measurements are highly correlated among three different approaches, reflectance values are much lower than other approaches, and offered several explanations, including the presence of anterior reflectors in the retina and lens scatter.^{36} The low OD_{MP} values obtained in our subjects are consistent with these previous findings obtained by reflectometry.^{36} Methods to correct these errors include baseline spectral measures of lens scatter^{36} and model-based estimates of (anterior reflectance) foveal reflectance using four spectral lines.^{45} These corrections could be implemented in future work. Another limitation may be related to the possible chromatic aberrations and defocus inherent in the spectral range at issue (420 to 500), although our dataset showed essentially uniform image quality in this range (Fig. 3). The optical density of the individual eye in principle should not confound the NMF hyperspectral method for OD_{MP}, because it is contained in OD_{SubT}, the separate term for other absorbers. Absorption by incompletely bleached cone photoreceptor pigment could also influence the hyperspectral data, but again, in principle, these effects should be separated from the MP by the NMF method, which specifically isolates the contribution from the MP spectrum. In fact, there could be imperfect unmixing of these signals. In the case of cones, the peak absorbances of the short and medium wavelength cones are centered about 50 nm on either side of the peaks of the MP profile,^{52} hence, they would have only a relatively small effect, similar to that of photoreceptor absorbance on HFP depending on differential absorption of blue and green light. Finally, the internal limiting membrane reflectance signal is apparent in the clinical photos of the younger normal subjects N2 and N6, and to a lesser extent in other images, overall, the effect of scattering due the internal limiting membrane is small.^{51}

It is important to note that approaches using differential reflectance or AF images to calculate MP density maps rely on major assumptions. For reflectance, the entire difference in reflectance between a single blue and green wavelength is attributed to MP. The NMF method makes no such assumptions and instead uses the full spectral information in the 27 bands between 420 and 530 nm to fit the reflectance data to a detailed prior model of the MP spectrum. Similarly, the AF method makes several assumptions. First, the difference between AF intensity centrally and peripherally is entirely due to the absorbance by MP of blue light. Second, the relative spectral energy of lipofuscin fluorescence is constant across the central retina and, last, that MP density is zero at a peripheral point, another potential source of error for all comparison methods. Again, the strength of the NMF method is that it directly searches for the absorber of interest based on its own signature and specifically separates the distribution of this signature from the other components of the data. No assumptions are made or required about the distribution of other absorbers and reflectors.

There are potential advantages and disadvantages to the general NMF method and specifically to the use of an initialized choice of spectral basis vector. The power of NMF lies in its ability to extract meaningful spectra from data without prior knowledge. For example, this would apply to the reflectance spectra of drusen, which are unknown. In such cases, the initial basis vectors are specified in number only (the rank of the decomposition) and otherwise, the initial spectra can be virtually random. This approach is the broadest but with no guarantee that meaningful spectra will be recovered. However, NMF is also used in the context where a specific approximate spectral signal is known to be present in the data. Thus, the power of the NMF approach with prior knowledge is that, by specifying one initial basis vector for the *in vivo* spectrum of MP, the algorithm is already “close” to a solution. During the solution, the algorithm modifies the initial basis vector and calculates its spatial distribution as necessary to find the best objective fit to real closely-spaced spatio-spectral data. Recovery of both a physiologically reasonable MP spectrum and density profile in turn, provides substantial evidence that the NMF method has effectively identified the MP and separated it from other ocular absorbers, even though no other absorbers were explicitly modeled.

Realistically, no spectral separation will be perfect. Some of the abundance images for MP in normal eyes show signal in the retinal vessels, suggesting that hemoglobin absorbance has not been completely separated from the MP spectrum. Indeed, because we did not explicitly include hemoglobin absorbance spectra in our other basis vectors, a slight hemoglobin admixture would be expected. Moreover, the MP distributions in other published methods highlight the vessels,^{20, 23, 27} suggesting that a complete separation of hemoglobin and MP has not yet been achieved. A simple future solution would be an initial segmentation of the vessels, either manually or by automated techniques. A more sophisticated solution, directly amenable to NMF, would be to search for hemoglobin basis vectors.

A future goal will be to separately measure L and Z. Their spectra are similar, with Z shifted to the red by 10 nm. We believe that this should be possible to achieve with our method, as suggested by the fact that separate L and Z optical densities have been measured using fundus reflectometry with 5.8 nm resolution *in vivo*
^{32} (although the actual spectra were not recovered in this experiment). Our method has shown the ability to resolve peaks 10 nm apart, hence, it should also be able to recover the L and Z spectra.

In summary, we have shown that the combination of a novel hyperspectral imaging tool with advanced mathematical analysis can, to the best of our knowledge, for the first time, recover detailed spectral absorption curves for MP *in vivo* that correspond to physically realistic retinal distributions. The calculated optical densities are consistent with other reflectometry methods. With regard to MP specifically, the ease, speed, and precision of MP density acquisition suggests that further effort with known methods is warranted to correct for error sources such as lens scatter and pre-retinal reflectors. More generally, success with MP suggests that this hyperspectral method will prove effective in recovering known and unknown spectra and densities of important ocular reflectors and absorbers. Among these, which should be investigated, are drusen in age-related macular degeneration, whose complex biochemical composition is critical to understanding the disease process.

## Acknowledgments

The authors would like to thank Jennifer Dalberth for editorial assistance, the New York Community Trust (RTS), NEI R01 EY015520 (RTS), and R01 EY021470 (RTS, AAF) for Research support, and Reichert, Inc. (AAF) for CR: Grant support.

## References

*Identification of spectral phenotypes in age-related macular degeneration patients*, M. Fabrice, G. S. Per, H. Bruce, E. S. Michael, and B. Arthur, Eds., p. 64261I, SPIE (2007). Google Scholar

*in vivo*using hyperspectral image analysis," Journal of Biomedical Optics 16(10), 106008 (1 October 2011). https://doi.org/10.1117/1.3640813