High-resolution wide-field imaging of retinal and choroidal blood perfusion with optical microangiography

The importance of ocular circulation to the function of the eye and in eye diseases has been known for years. There is convincing evidence that vascular pathologies play an important role in the etiology and progression of a number of ophthalmic diseases, including glaucoma, diabetic retinopathy, and age-related macular degeneration (AMD). Clinically, a better visualization of the ocular vasculature could provide important diagnostic information for these diseases. Currently, the most widely used diagnostic tools for imaging ocular blood vessels are fluorescein angiography (FA) and indocyanine green angiography (ICGA). Although these two methods can provide visualization of retinal (using FA) and choroidal (using ICGA) blood vessels, they require injection of a contrast agent (sodium fluorescein dye and indocyanine green dye) into the blood stream. Further, sodium fluorescein may cause nausea, vomiting, itching, dyspnea, and additional adverse allergic reactions,¹ while indocyanine green dye must be used cautiously in patients who are allergic to iodine and shellfish.² Another main disadvantage is that neither of these two methods can provide depth information on ocular blood vessels.

Optical coherence tomography (OCT)^3- 4 has become an important tool in ophthalmology for the diagnosis, management, and high-resolution study of the anterior and posterior segments of the eye. Based on low-coherence interferometry, this noninvasive cross sectional imaging technique offers resolution in the micron range, both latterly and axially, with an ability to image to a depth of several millimeters depending on the opacity of the specimen. Combining OCT with the Doppler principle, phase-resolved Doppler optical coherence tomography (PRDOCT) has successfully been demonstrated for in-vivo imaging of blood vessels and the extraction of functional flow information from blood vessels buried within tissue,^5- 6 particularly after the advent of Fourier domain (FD) OCT.^7- 11 PRDOCT is based on the evaluation of phase difference between adjacent A-lines. However, the accuracy of PRDOCT is deteriorated by two factors: system sensitivity¹² and the heterogeneous texture pattern of tissue.¹³ Due to the degradation of depth-dependent sensitivity, it is hard for PRDOCT to provide blood vessel networks in the choroid, where the signal-to-noise ratio is relatively low compared to the retina. Moreover, the heterogeneous texture pattern artifact reduces the image quality of PRDOCT even in the retina layers, causing difficulty in providing capillary level information. To overcome these problems, several methods have been proposed. Szkulmowski et al. reported a functional Doppler OCT method called joint spectral and time domain OCT.¹² This technology relies on analysis of the amplitude of the OCT signal to extract the flow information buried within the tissue, and is supposed to have higher sensitivity at low signal-to-noise ratio regions than the conventional phase resolved method. A problem with this method is that it needs high A-line density along one B-scan, limiting its suitability for acquiring ocular blood vessel networks for in-vivo applications. Wang and Ma¹³ proposed that by employing a reverse scan pattern, the heterogeneous texture pattern of the tissue can be successfully suppressed. This method showed a great potential to increase the velocity sensitivity of the PRDOCT method to capillary level. However, an additional reverse scan required at the same location restricts its applications for in-vivo experiments. Besides these functional Doppler methods, Yasuno et al. ¹⁴ demonstrated a segmentation method, the scattering optical coherence angiography method (S-OCA), which is based on the analysis of the low intensity region in the OCT structural image of the choroid of the human eye. Though this method can image blood vessel networks in the human choroid, it essentially loses the functional information of blood flow, such as flow velocity and direction. Moreover, the contrast of the S-OCA image is provided by the scattering OCT signal, which makes it difficult to distinguish the flow signals from the structure signals.

Recently, based on the principle of full range complex FDOCT^15- 17 with constant carrier frequency modulation along the B-scan direction of an interferogram, Wang and An introduced a new label-free optical microangiography (OMAG) technique to separate the static and moving signals emerged from an illuminated tissue sample. This method was successfully demonstrated for in-vivo imaging of blood vessels at capillary-level resolution.^18- 19 By introducing a phase compensation method to compensate the movement artifact caused by head or eye movements, the OMAG method was also applied successfully for in-vivo imaging of ocular blood vessel networks both in human retina and choroid.²⁰ Although OMAG can deliver superior performance for depth-resolved imaging in both retina and choroid, the first version of the OMAG algorithm encountered certain limitations caused by the introduction of constant modulation frequency. First, the introduction of modulation frequency by the system hardware is difficult to readjust for different experimental situations, such as the density of the A-line scan in one B-scan and the flow direction. Another issue is that this modulation frequency will be coupled into the Doppler frequency shift, which is normally used to evaluate the velocity of the flow, leading to difficulty in retrieving the information regarding the true velocity.¹⁸ To solve the first issue, Wang presented a digital modulation frequency method for OMAG technology.²¹ But this digitally introduced modulation frequency demands computation power to obtain velocity information of OMAG flow images. To solve these two issues, a modified Hilbert transform method was introduced into OMAG.^22- 23 By eliminating the modulation frequency, the speed and flexibility for imaging blood flow can be greatly increased.^23- 24 We have extended that work to perform functional imaging of the choroidal circulation and to obtain wide-field imaging of the ocular vasculature, which is important for successful use in the clinical setting.

In this work, we used OMAG with the modified Hilbert transform method to obtain ocular blood vessel network images. We carried out wide-field imaging of ocular perfusion of both the retina and choroid by reconstructing multiple sequentially registered OMAG projection maps of small areas $(2×2mm)$. For the accurate reconstruction of wide-field OMAG imaging, there are two main issues that need to be resolved in the process of sequential capturing of multiple OMAG images. First, for each of these sequential OMAG datasets, the bulk motion artifacts in each B-scan image need to be suppressed. We solved this issue by directly implementing the phase compensation method.²⁰ In this process, rather than using the phase averaging approach reported in Ref. ²⁰, we instead used the histogram method to estimate the bulk motion phase. The second issue is the displacement between adjacent B-scans. To solve this problem, we used the cross correlation between adjacent 2-D OMAG flow images to compensate for the displacements. After resolving these artifacts, we presented wide-field images of retinal and choroidal blood flow maps, and showed the superior performance of OMAG technology to obtain depth-resolved volumetric ocular perfusion of both retina and choroid. Finally, we compared the wide-field OMAG perfusion map with the clinical FA and ICGA images from the same individual. The well-correlated results show the great potential of OMAG in future clinical applications.

The schematic of the system used to achieve the wide-field ocular blood perfusion image is illustrated in Fig. 1. The light source is a superluminecent diode with a spectral bandwidth of $46nm$ centered at $842nm$ that yields an axial resolution of $∼8μm$ in air. The light emitted from the light source was coupled into a Michelson interferometer, which divides the incident beam into the reference arm and sample arm. The sample light was coupled into a custom-made imaging probe with a pair of galvometer scanners to achieve $x-y$ scanning and to collect the light backscattered from the posterior part of the eye. The optical power of the beam on the cornea was $∼700μW$, which is lower than the ANSI exposure limit. Together with the human eye lens, the imaging probe includes a collimating lens, objective lens, and an ocular lens that provided a lateral resolution of $∼16μm$. In the reference arm, to compensate the dispersion, light was delivered into the stationary mirror through an $∼2-cm$ water cell and a set of optics that is matched with the sample arm. Unlike our previous work, we do not need to introduce a constant modulation frequency into one OCT B-scan. During imaging, the reference arm mirror was kept stationary. The interference signal between the reference arm and sample arm was then collected and delivered into a home-built high-speed spectrometer. This spectrometer had a spectral resolution of $0.055nm$ that achieved a maximum imaging depth of $∼3mm$ on each side of the zero delay line with a line camera of $2048pixels$.

For maintaining high-resolution imaging of blood vessels, the area of each OCT scan was set relatively small. In the lateral direction ($x$ direction), the $x$-galvoscanner was driven by a $27-Hz$ sawtooth waveform with $2-mm$ amplitude, i.e., the imaging speed at 27 frames per second. For elevational direction ($y$ direction), the $y$-galvoscanner was driven by a $0.27-Hz$ sawtooth waveform with $2-mm$ amplitude. In this study, the camera integrating time was set to be $36μs$, allowing $1μs$ for transferring the spectral data from the CCD camera ($2048pixels$, A-scan) to the host computer via Cameralink^™ and a high-speed frame grabber board. For each B-scan we captured 1000 A-scans, and for C-scan (elevational direction) we captured 100 B-scans. Hence, one subarea image dataset contained 2048 (A-scan) by 1000 (B-scan) by 100 (C-scan) pixels.

To achieve wide-field imaging of fundus vasculature $(7.4×7.4mm)$, we sequentially scanned different areas of the posterior part of the human eye by a scanning order illustrated in Fig. 2. In Fig. 2, the background is an OCT fundus image obtained through the technique described in Ref. ²⁵. To obtain this background image, the scan area of our OCT system was set at $∼9×9mm$. The small blue squares represent the scan area of each smaller OCT dataset with the number indicated referring to the scanning order in the whole experiment. Between every two subarea scans, we introduced $∼0.2-mm$ overlap area to make sure no blank area occurred when reconstructing the wide-field images. For every subscan, it took around $3.7s$ to capture the dataset. Between two subscans, the volunteer was instructed to close his eye for $15s$ for a rest. Hence, it took about $5min$ to acquire the wide-field view of the fundus blood vessel maps. During the whole experiment, the volunteer fixated on a blue light spot emitted by the aiming source to minimize eye movement. The operating software for the whole system, including beam scanning, data acquisition, data storage, and hand-shaking between them, was developed on a Visual C++ platform.

Based on the modified Hilbert transform, OMAG can separate the flow signals from the static tissue signals without introducing any modulation frequency. The detailed description of this method has been reported previously.^23- 24 Here, we briefly review the principle of this method.

The B-scan interference signal captured by the OCT system can be written as:²³Display Formula

In Eq. 1, the first term is the dc signal produced by the light reflected from the reference mirror. The second term is the spatial frequency component of the static tissue sample, which we call the heterogeneous frequency component, and which can provide the static structural information of the sample. Because of the optically heterogeneous nature of the biological tissue sample, this frequency is randomly distributed around zero frequency with a bandwidth of $BW$.¹⁶ The bandwidth of this heterogeneous frequency is determined by the optical characteristics of the tissue sample and the density of A-line scans in a B-scan image.²⁶ The third term is the Doppler beating signal, which is introduced by the moving particles in the tissue sample. Depending on the magnitude of the moving particle velocity, the Doppler beating signal shifts away from the heterogeneous bandwidth $BW$. If we first apply high-pass filtering with cutting frequency $BW∕2$ and then apply the Hilbert transformation along the B-scan direction, we can extract the negative $(v<0)$, Eq. 2, and positive $(v>0)$, Eq. 3, directional flows as follows:¹⁸Display Formula

The movement artifact is one of the most distracting issues for in-vivo imaging of the posterior part of the human eye. This artifact is caused by the involuntary movement of head or eye and can seriously deteriorate the image quality. In our previous work,²³ we proposed a phase compensation method to effectively resolve this issue in the traditional OMAG system. In our modified Hilbert transform-based OMAG system, we directly implemented this phase compensation method to minimize the movement artifact issues.

When the sample moves during imaging, the velocity of this movement $vs$ will introduce a Doppler frequency shift into the light backscattered from the sample with respect to the reference light. If we analyze the frequency components in one B-scan along the $x$ direction, then this Doppler frequency due to the bulk motion will also represent a Doppler beating frequency shift. When this frequency shift is large enough $(>BW∕2)$, the bulk motion signal can be detected through the OMAG method that inevitably gives a motion artifact in the final OMAG flow image. Figure 3 shows typical imaging results obtained from the OMAG method. Figure 3 presents a typical cross sectional structure image of the posterior part of a human eye from a volunteer in our laboratory. The structure information of different layers, such as retina, choroid, and sclera, are clearly delineated from this image. Figure 3 is the corresponding cross sectional flow image, directly obtained from the OMAG method without applying the phase compensation method. Though two large blood vessels are highlighted by the OMAG method, the other part of the image is corrupted by the motion artifacts (white arrow). In this case, the detected moving signal can be written as:Display Formula

As we know that the OMAG image has its origin in FDOCT technology, the phase term $ϕM$ may be obtained from the FDOCT complex signal $B̃OCT(t,z)$ by applying the fast Fourier transform (FFT) along wave number $(k)$:Display Formula

As we discussed before, a key step for the phase compensation method is how to estimate the phases due to bulk tissue motion between adjacent A—lines $ΔϕSi$. In our previous work,²³ we estimated $ΔϕSi$ by averaging the phase differences between adjacent A-lines along the $z$ direction within the tissue-signal region of an OCT A-scan signal. This method was based on the assumption that the blood vessels within the cross section of the posterior eye segment are quite small compared to the static tissue region itself. This assumption is true for most parts of the posterior segment, except the optic disk, where the blood vessels are relatively large and in which the blood flow is also very fast. However, the phase changes induced by the flow in large blood vessels within the optic disk are not negligible when using the phase-averaging approach,²³ which instead could impose an additional phase offset to the bulk motion estimation. This additional phase offset will lead to a correction artifact in the final output flow image. Such an artifact is illustrated in Fig. 3, resulting from the phase compensation based on the phase averaging method. As we can see, in the left part where there are no large blood vessels, the bulk motion artifacts are eliminated quite successfully. However, in the right part, the existence of the large blood vessel (blue arrow) causes a serious correction artifact above it (pointed by the red arrow). To solve this problem, a more accurate method is therefore needed for the phase estimation due to the bulk tissue motion.

We note that the phase changes solely due to the bulk motion are almost the same at all positions within the tissue region along one OCT A-scan. On the other hand, the blood vessels often occupy a small proportion of the tissue region along one A-scan, and also the velocity of blood flow is parabolically distributed across a wide range. These facts determine that the number of possible positions carrying the phase caused by the bulk motion will be much larger than those of the positions with other phase values along one A-scan. For this reason, we adopted a histogram approach to estimate the phases that are caused by the bulk tissue motion. Next, we briefly describe the procedures of this histogram method, which are similar to that developed by Makita et al. ¹¹

After estimating the phase differences $Δϕi(z)$ between the adjacent A-scans within the tissue region by PRDOCT method, the maximum value of phase difference $Δϕmaxi(z)$ and the minimum value of phase difference $Δϕmini(z)$, as well as the number of phase samples $m$, were calculated. The bin width $h$ for the histogram was then determined by the Freedman and Diaconis rule:²⁷Display Formula

In the case of wide-field imaging, another serious problem caused by the involuntary movement of the human head and eye is the displacement between adjacent B-scans. During imaging, although the volunteer was asked to fixate the blue aiming spot to minimize eye movement, the motion artifacts are still quite dominant and occur occasionally through the whole progress of capturing the 16 OCT datasets, from which the wide-field ocular perfusion map is reconstructed. As previously reported,^11,20 the axial displacement can be compensated for with relative ease by the cross correlation method between two A-lines at the same position of two adjacent B-scans. This is because the retina has well-layered structures along the depth direction that produce a sharp peak for the cross correlation function. Thus, this cross correlation is suitable for determining the displacement by locating its peak position. Using the position of peak value, we can determine the lateral displacement between these two B-scans. For example, if the position index of the peak value is zero, there is no displacement between them. If the position index is positive, the second B-scan moved to the right compared to the previous B-scan and vice versa. However, for several reasons, the cross correlation method based on an OCT structural image often does not work effectively well to compensate the lateral displacements for most parts of the posterior part of the human eye. First of all, except for the optic nerve head region, the structures of the retina are very smooth. Particularly when a small area of $∼2×2mm$ is scanned, the variance of structures is very small in the lateral direction. This means that the structural information presented by the same layer will be very similar along the lateral direction. Because of these reasons, the calculated cross correlation function will behave as if having a very broad peak, which leads to difficulty in locating the peak position for estimating the lateral tissue motion. To give an example, Fig. 5 shows a typical FDOCT B-scan structure image. The cross correlation between the same depth position marked by the red line of two adjacent B-scans is shown in Fig. 5. As it can be seen, the cross correlation function has a very broad peak, and thus the subsequent determination of the peak position can be inaccurate. To obtain accurate determination of the peak position, we propose to use the OMAG flow images to evaluate the cross correlation function, rather than using the structural image. This is reasonable because the OMAG flow images possess more high spatial frequency components than the structural images. According to the information theory as applied to the pattern recognition; the more high spatial frequency components, the sharper the correlation peak. Figure 5 gives the OMAG flow images corresponding to Fig. 5. The cross correlation function at the same position as for Fig. 5 is shown in Fig. 5. It is clear that the use of OMAG flow images produces a cross-correlation function with a very sharp peak, thus delivering an accurate estimation of the displacement between adjacent B-scans.

Based on the prior analysis, we used the cross correlation between adjacent OMAG flow images to obtain an estimation of axial displacement between adjacent B-scans. We also used the same approach to estimate both the axial and lateral displacements directly by 2-D cross correlation method between two neighboring OMAG frames. Figure 6 shows the comparison between the images without and with motion compensation. Figures 6 are OMAG fundus vasculature images for one 3-D OCT dataset before and after displacement compensation, respectively. As we can see in Fig. 6, the two blood vessels in this data volume have broken into several parts because of the eye motion, marked by the blue circles. After the motion compensation, these two vessels were reconnected in Fig. 6. This method works quite well for compensating lateral displacement, even in a small area without any structural landmarks. Figure 6 is the elevational cross sectional image before motion compensation, which was extracted from the position marked by the red line in Fig. 6. Figure 6 is the corresponding image after the motion compensation. The position of this image is also marked by the red line in Fig. 6, which is the same position as in Fig. 6. By comparing these two images, we can see that Fig. 6 is smoother than Fig. 6, where a big strike (the white arrows) is successfully compensated to a level similar to other positions. We further discuss the motion compensation results in Sec. 4, based on the 2-D OMAG $x-y$ projection images.

An important characteristic of OMAG, which is inherited from FDOCT technology, is the depth-resolved imaging capability. Through the OMAG method, blood vessels located at different depths, for example in the retina and choroid, can be separated relatively easily. To achieve the isolation of blood vessel located in retina and choroid, two masks were created through a segmentation process based on the FDOCT B-scan structure image. The detailed procedure of the creation of the masks for both retina and choroid has already been described in our previous work.²⁰ These two masks were then applied onto the OMAG flow image to separate the blood vessels in retina and choroid. After combining 100 B-scans along the elevational direction (C-scan) and doing the maximum projection, the vascular perfusion maps for both the retina and choroid are finally obtained.

After compensating the axial and lateral displacement, we could obtain 3-D vascular images by stacking 100 B-scans together. Through the two masks mentioned before, we could also get the 3-D volumetric image of retinal and choroidal blood vessel networks. The results are presented in the following section.

Because of the high sensitivity to the flow, OMAG delivers superior performance for imaging the blood vessel networks both in retina and choroid. Especially after displacement compensation, even the microvascular networks could be effectively presented through this method. In this section, we give a detailed demonstration of the results, obtained from the previous described procedure for in-vivo experimental 3-D datasets captured by the system described in Fig. 1. Each 3-D FDOCT interferogram is composed of $2048pixels$ in one A-line ($z$ direction), 1000 A-lines in one B-scan ($x$ direction) and 100 B-scans along one C-scan ($y$ direction). This dataset covered $∼2mm$ ($x$ direction) $×2mm$ ($y$ direction) area information of the posterior part of the volunteer eye. The data acquisition time for such a dataset was $<4sec$, with the camera line scan rate at $27kHz$.

Figure 7 is a typical FDOCT B-scan structure image. According to the structural image, two masks were created from the procedure described in Sec. 3d The first one is the anterior boundary presented by the red line and posterior boundary presented by the yellow line, shown in Fig. 7, which is used for vasculature reconstruction in retina. The second one is the anterior boundary presented by the green line and posterior one presented by the blue line in Fig. 7, which is used for choroidal blood vessel reconstruction. Figure 7 is the flow image obtained by applying the OMAG algorithm on the same cross sectional position of Fig. 7. After mirroring the four boundaries in Fig. 7 onto Fig. 7, we successfully segmented out the blood vessels of retina and choroid, respectively. In the retinal region, between the red and yellow line, not only the large vessels but also the microvessels are demonstrated. The choroidal region, between the green and blue lines, presents an abundance of blood vessels. Combining 100 frames along the elevational direction and performing the maximum projection, the ocular vascular projection image for 3-D scanned tissue volume can be obtained. Such results are shown in Fig. 8.

Figures 8 are the ocular perfusion maps of retina and choroid before applying the motion compensation method. Though both of them present abundant blood vessels, lots of them (marked by the ellipse) are not connected well with each other because of eye movements. Figures 8 are the vascular perfusion maps of retina and choroid after applying the displacement compensation. The broken blood vessels presented in Figs. 8 due to motion are reconnected again, with similar areas marked by the ellipses. Especially in the region marked by the green ellipse in the retina and choroid, the blood vessel networks are reasonably well presented compared to those in Figs. 8.

After motion compensation, we could obtain the 3-D volumetric flow image of both retina and choroid by combining 100 B-scan frames in each 3-D subdataset. The results are shown in Fig. 9. Figures 9 demonstrate the 3-D volumetric blood vessel networks in the retina and choroid. In Fig. 9, the retina blood vessels (coded with green) and choroidal blood vessels (coded in red) are merged together to present their relative positions. Together with a cross section B-scan structure image, Fig. 9 also demonstrates the positions of blood vessels relative to the structure.

As we discussed in the first section, the wide field of view ocular perfusion maps of retina and choroid may provide a broad spectrum of information for ophthalmic diagnosis and eye disease studies. To obtain the wide-field imaging of ocular perfusion maps while maintaining the high resolution OMAG imaging of blood perfusion, we captured 16 OCT data cubes in a predesigned order (see Fig. 2). Each of these small datasets covered around a $2×2-mm$ area. In both lateral and elevational directions, we introduced a $∼0.2-mm$ overlap between adjacent small datasets to eliminate blank space during the whole image stitching process. In this case, all 16 datasets covered an entire area of about $7.4×7.4mm$ over the posterior segment of the eye. For every OCT raw data, we first obtained all cross sectional structural images, i.e., OCT images, to segment out the tissue signal region from which the retina and choroid layers, through the creation of two masks for retina and choroid, were separated. Subsequently, the motion compensation methods described in Secs. 3b,3c were applied to remove/minimize the motion artifacts due to the inevitable eye movement during imaging. Next, the maximum projection approach along the $z$ direction was used to obtain the projection images $(x-y)$ of the separate vasculature maps of retinal and choroidal blood vessels. Finally, all the projection images of the 16 areas were manually stitched together to form two wide-field vascular perfusion maps of retina and choroid. The results are illustrated in Figs. 10, respectively. Figure 10 is the retina blood perfusion map and Fig. 10 is the choroid blood perfusion map. Both images demonstrate clear, well-connected blood vessel networks, not only the large blood vessels but also the capillary vessels.

At last, to compare the performance of OMAG technology with standard clinical ocular perfusion imaging modalities, we performed FA and ICGA procedures for the same volunteer on whom we performed OMAG scans. The results are presented in Figs. 10 for FA and ICGA, respectively. In Figs. 10, the light-blue squares correspond to the scan area of OMAG imaging shown in Figs. 10. By comparison, the OMAG results demonstrated very good agreement with the clinical FA and ICGA results. Through scrutinizing the appearances of the blood vessels, almost all noticeable vessels presented by FA and ICGA images are demonstrated by OMAG imaging. Furthermore, in the retinal perfusion map [Fig. 10], the OMAG results demonstrated more microcirculation information in the fovea region than that of the FA image [marked by the ellipse in Fig. 10]. The main advantages of OMAG compared to FA and ICGA are that the images are obtained without injected contrast agents, and the vessels can be correlated with the depth of location in the retina or choroid. These two advantages of OMAG wide field of view perfusion maps indicate this imaging modality will potentially be a very powerful tool for in-vivo clinical application in the future.

We use OMAG technology to achieve wide-field imaging of ocular perfusion maps, both in retina and choroid, by manually stitching 16 small area datasets. The datasets provide high-resolution blood vessel maps down to the capillary level. To achieve this, we develop methods to solve the problems associated with subject movements. We use a histogram-based method to estimate the bulk motion phase by which the bulk motion artifacts are eliminated/minimized through a phase compensation method. To remove the displacement between adjacent frames caused by involuntary head and eye movements, we demonstrate a method that uses the correlation function between two adjacent B-scan OMAG flow frames to estimate the displacement. We show that these motion compensation methods work well to remove/minimize the subject movement artifacts presented in OMAG flow images. Finally, we compare the wide-field OMAG perfusion maps with the currently available FA and ICGA images, and the well-correlated results show the great potential of OMAG in future clinical applications.

This study was supported in part by research grants from the National Institutes of Health (R01 HL093140-01, R01 EB009682-01, and R01 DC010201-01), the American Heart Association (0855733G), and an unlimited grant from Research to Prevent Blindness. The content is solely the responsibility of the authors and does not necessarily represent the official views of grant-giving bodies.