Multi-reference global registration of individual A-lines in adaptive optics optical coherence tomography retinal images

Kazuhiro Kurokawa; James A. Crowell; Nhan Do; John J. Lee; Donald T. Miller

doi:10.1117/1.JBO.26.1.016001

6 January 2021 Multi-reference global registration of individual A-lines in adaptive optics optical coherence tomography retinal images

Kazuhiro Kurokawa, James A. Crowell, Nhan Do, John J. Lee, Donald T. Miller

Author Affiliations +

Journal of Biomedical Optics, Vol. 26, Issue 1, 016001 (January 2021). https://doi.org/10.1117/1.JBO.26.1.016001

Abstract

Significance: Adaptive optics optical coherence tomography (AO-OCT) technology enables non-invasive, high-resolution three-dimensional (3D) imaging of the retina and promises earlier detection of ocular disease. However, AO-OCT data are corrupted by eye-movement artifacts that must be removed in post-processing, a process rendered time-consuming by the immense quantity of data.

Aim: To efficiently remove eye-movement artifacts at the level of individual A-lines, including those present in any individual reference volume.

Approach: We developed a registration method that cascades (1) a 3D B-scan registration algorithm with (2) a global A-line registration algorithm for correcting torsional eye movements and image scaling and generating global motion-free coordinates. The first algorithm corrects 3D translational eye movements to a single reference volume, accelerated using parallel computing. The second algorithm combines outputs of multiple runs of the first algorithm using different reference volumes followed by an affine transformation, permitting registration of all images to a global coordinate system at the level of individual A-lines.

Results: The 3D B-scan algorithm estimates and corrects 3D translational motions with high registration accuracy and robustness, even for volumes containing microsaccades. Averaging registered volumes improves our image quality metrics up to 22 dB. Implementation in CUDA™ on a graphics processing unit registers a 512 × 512 × 512 volume in only 10.6 s, 150 times faster than MATLAB™ on a central processing unit. The global A-line algorithm minimizes image distortion, improves regularity of the cone photoreceptor mosaic, and supports enhanced visualization of low-contrast retinal cellular features. Averaging registered volumes improves our image quality up to 9.4 dB. It also permits extending the imaging field of view (∼2.1 × ) and depth of focus (∼5.6 × ) beyond what is attainable with single-reference registration.

Conclusions: We can efficiently correct eye motion in all 3D at the level of individual A-lines using a global coordinate system.

1. Introduction

Adaptive optics-optical coherence tomography (AO-OCT) is a non-invasive, high-resolution method for obtaining three-dimensional (3D) volumetric images of the retina.¹^–⁴ These qualities may support earlier detection and improved monitoring of prevalent blinding retinal diseases, such as glaucoma, diabetic retinopathy, and age-related macular degeneration.¹^–⁴ However, analysis of AO-OCT data is impeded by two factors: First, the time to acquire a single volume is sufficient for it to be distorted by eye motions such as tremors, drifts, and microsaccades,⁵ creating artifacts that are many times larger than the cells being imaged.⁶ These motion artifacts must be measured and compensated for to permit proper analysis of the retinal tissue images. We address this problem by proposing a cascade of two registration algorithms that we call 3D B-scan registration⁷ and global A-line registration. Second, high-resolution AO-OCT systems generate terabytes of image data from a single imaging session, resulting in immense data sets and making motion artifact removal a lengthy process. We decrease processing time by implementing the 3D B-scan registration algorithm with parallel processing on a general-purpose graphics processing unit (GPGPU).⁷ Together, the two registration algorithms allow us to compensate efficiently and accurately for both lateral and axial eye motions. The algorithms correct for $X Y Z$ translational eye movements, torsional eye movements, and variations in image magnification, and generate global motion-free coordinates.

Numerous algorithms have been developed for registering images from adaptive optics scanning laser ophthalmoscopes (AO-SLO)⁸^–¹¹ and AO-OCT.¹²^–¹⁴ These are based universally on strip-wise registration of two-dimensional (2D) en face images. Strip-wise registration specifically corrects motion artifacts that appear in the en face image, which is the most common view of the retina obtained with AO ophthalmoscopes. This approach typically requires that images are decomposed into strips that are sufficiently narrow and acquired sufficiently quickly that they are individually free of motion artifacts (aside from an occasional microsaccade). Motion artifacts are therefore assumed to consist of rigid displacements of these strips that can vary in magnitude and direction. These displacements can be measured and corrected in a target volume by cross-correlating each strip from the en face image of that volume to multiple strips from the en face image of a reference volume.

When applied to volumes, strip-wise approaches are computationally inherently efficient in that the entire volume is registered by processing only a small fraction of the volume data, i.e., a 2D en face projection, typically through the photoreceptor layer. This layer has been historically used for registration as it contains the brightest reflections and hence the strongest signals. However, computational efficiency comes at the expense of discarding most of the information in the volume, which could be used for improving the registration. For our 3D B-scan registration algorithm, we take advantage of this additional information by extending the cross-correlation regions into the third dimension (depth), correlating B-scans instead of en face strips. Making use of the additional information outside the photoreceptor layer confers two important advantages. First, it permits us to trade an increased extension in-depth for a decrease in width. Specifically, instead of single-pixel-deep strips we use fast B-scans as single-pixel wide slices extending the full depth of the volume. This allows us to better detect and correct rapidly varying motion artifacts. Second, the depth of focus of AO ophthalmoscopes is typically a small fraction of the retinal thickness. When we image cells outside the photoreceptor layer, we focus the system on those cells, decreasing the sharpness of en face images of the photoreceptors. Cross correlation of B-scans is inherently less sensitive to changes in focus as all cells contribute to the registration and those cells in focus have increased signal and sharpness.

Effective strip-wise and B-scan registration require the reference volume to contain minimal motion artifacts and maximal overlap with the image set to be registered. Considering that no volume is ever completely free of motion artifacts, it follows that volumes registered by this method are still distorted. While several approaches have been proposed to mitigate such residual distortion, e.g., incorporating an eye-tracking system⁶^,¹⁵^–¹⁷ and modifying the scan pattern,¹⁸^,¹⁹ these may require additional hardware and add complexity and cost. Stevenson and Roorda⁸ and others⁹^,¹⁰^,¹⁴ presented a clever method for constructing a synthetic global reference from the complete set of motion traces for each strip within a set of registered images. This permits registration of all images to a global reference and can be done entirely in post-processing. This method is successful, though is most effective within the plane containing the 2D registration strips (e.g., photoreceptor layer) and does not compensate for both scaling and torsional errors even within that plane. The method is also limited to a single reference.

We develop a more rigorous approach to addressing this global registration problem. We use our 3D B-scan registration algorithm to register target volumes to a reference volume. We test this algorithm by examining its effect on a number of image quality metrics. We then repeat the 3D B-scan registration by selecting multiple reference volumes to increase the number of registerable volumes. We use these registration results to derive a single global coordinate system that all individual B-scans from all volumes are mapped into. Finally, we correct for scaling and torsional errors to estimate 3D global coordinates for individual A-lines. We name this approach our global A-line registration algorithm. We test the global A-line algorithm by applying a slightly different set of image quality metrics to averaged volumes registered to the global coordinate system; this second set of metrics are intended to also measure enhancements in the regularity of the resulting images. Figure 1 shows a detailed flowchart of our registration method with steps referenced by section number in the paper.

Fig. 1

Flowchart of our registration method that cascades (a) a 3D B-scan registration algorithm with (b) a global A-line registration algorithm. The 3D B-scan registration algorithm registers multiple target volumes $T = (1, 2, \dots)$ to a single reference volume $R$ . A coarse-to-fine registration approach improves processing speed and accuracy (Sec. 3.1.2.3), accelerated by parallel computing (Sec. A in Supplementary Material). Global A-line registration registers individual A-lines in all volumes to a global coordinate system. The approach consists of three steps: (1) multiple-reference registration to separately register all target volumes to each of multiple reference volumes (Sec. 4.1.1), (2) global registration of individual B-scans (Sec. 4.1.2) to derive a single mapping of individual B-scans to the multiple-reference global coordinate system, and (3) global registration of individual A-lines (Sec. 4.1.3) to correct en face affine transformations among volumes and estimate 3D global coordinates for individual A-lines. Abbreviations: POC, phase-only correlation; NCC, normalized cross-correlation.

Understanding our method and contrasting it against registration methods already in the literature requires an understanding of AO-OCT image acquisition and the impact eye motion has on it. These topics are introduced in Sec. 2. Section 3 presents our 3D B-scan registration algorithm and CUDA™ GPGPU implementation, which allows it to run $\sim 150$ times faster than it would in MATLAB™ on a central processing unit (CPU). We experimentally evaluate our implementation by running it against a dataset of more than 300 volume images and measuring registration accuracy and run time both quantitatively and qualitatively. Section 4 presents our global A-line registration algorithm, which uses the output of the 3D B-scan registration algorithm. We derive the method from theory and experimentally evaluate it by running it against a dataset of more than 900 volume images and measuring the registration accuracy qualitatively and quantitatively. We use the combined registration method (3D B-scan registration plus global A-line registration) to obtain pristine images of retinal cellular structure, to increase the field of view ( $\sim 2 \times$ ) over that of a single acquired AO-OCT volume, and to create through-focus AO-OCT volume images with preserved image sharpness over the full retinal thickness. The latter is made possible by increasing the effective depth of focus ( $\sim 5 \times$ ) over that of a single acquired AO-OCT volume. In Sec. 5, we summarize the performance of our registration method, compare it to registration methods in the literature, and discuss its limitations.

Finally, this study focuses on the development rather than the use of a new registration method. Because of this, we found it most efficient to illustrate and validate our method using just two large datasets acquired on two subjects who had typical fixational ability. We have numerous other datasets that could have been used, and in fact, we now routinely use this method in all our imaging studies. To date, we have successfully registered images acquired on more than 40 subjects of different age, gender, and eye size, with different retinal diseases, and at numerous retinal eccentricities and time points. Some of these results can be found elsewhere.⁴^,²⁰^–²⁵ Also, while we limit the method to AO-OCT images, the method should be applicable to other types of 3D images, most obviously to OCT images (without AO), which are universally used in eye care practice.

2. AO-OCT Image Acquisition and Impact of Eye and Head Motion

Each AO-OCT image volume contains a unique pattern of lateral $(x, y)$ and axial ( $z$ ) distortions and is non-rigidly displaced in all 3D relative to other volumes because of random eye motion during fixation⁵ and subtle movements of the head. The nature of these distortions in each dimension is determined by an interaction between the temporal properties of the component of eye movement and the component of the OCT scan pattern in that dimension. The pixels within a single A-line of the volume (typically oriented in retinal depth, $z$ ) are captured simultaneously. A fast B-scan $(x, z)$ consists of a set of A-lines spaced contiguously along the $x$ axis; in our system, they are typically collected over a period of 0.3 to 1 ms. A complete volume $(x, y, z)$ consists of a set of fast B-scans $(x, z)$ displaced along the $y$ axis and collected over 98 to 414 ms. A slow B-scan $(y, z)$ is constructed from A-lines at the same $x$ position along the $y$ axis, whereas a C-scan $(x, y)$ or en face image is constructed from rows of fast B-scans at the same depth and displaced along the $y$ axis.

When we view a 2D slice through a volume, the appearance of eye motion artifacts depends on which slice is viewed: fast B-scan $(x, z)$ , slow B-scan $(y, z)$ , or C-scan $(x, y)$ . For example, the most commonly seen artifacts in fast B-scan images $(x, z)$ is a tilt induced by a linear variation in the time-of-flight difference of the retinal reflections across the B-scan, as shown in the fast B-scan of Fig. 2 (top left). The time-of-flight variation is caused by displacement of the eye’s optical axis in its pupil plane relative to the imaging beam axis. However, in a slow B-scan image $(y, z)$ , head movement along the $z$ axis is the dominant effect, causing slow (mean wave period $\sim 1 s$ ) shifts or ripples of the image along the $z$ axis, as shown by the slow B-scan and slow B-scan projection of Fig. 2 (top). Microsaccades are more evident in the en face image $(x, y)$ and cause the largest and most abrupt shifts; they can traverse between several dozens and hundreds of photoreceptors over $\sim 25 ms$ ,⁵ as shown in the en face projection in Fig. 2 (bottom). Tremors in the $x$ and $y$ dimensions cause small oscillatory shifts ( $\sim 2$ to $3 μ m$ in amplitude) with a duration of $\sim 11 ms$ ,⁵ resulting in the distorted appearance of cone photoreceptors shown in Fig. 2 (bottom). While it is not evident in Fig. 2, drift movements cause moderate shifts—traversing roughly a dozen photoreceptors—over longer intervals ( $\sim 1 s$ ).⁵

Fig. 2

Eye motion artifacts in AO-OCT volume. (top left) A volume consists of a set number of fast B-scans spaced along the $y$ axis with each fast B-scan consisting of a set number of A-lines spaced along the $x$ axis. (top right) The slow B-scan projection image is the projection of the volume onto the $y z$ plane by averaging all pixels along the $x$ axis. The wavy looking pattern is due to axial motion. (bottom) Similarly, the en face projection image is the projection of the volume onto the $x y$ plane by averaging all pixels along the $z$ axis. The abrupt distortion is due to a microsaccade. The irregular spacing is due to tremor motion.

The process of image registration is also complicated by changes in head position that can give rise to scaling and rotation (torsion) of images. Measurable torsion across look-ins—caused by combination of changes in head position and compensatory torsional eye movements—is common,²⁶^,²⁷ even in subjects experienced in the use of a bite bar. Also, the torsion error drifts over time (standard deviation of torsion error is $\sim 0.15 \deg$ within a 4-s session, but an order of magnitude larger across sessions).²⁸ Image scaling is primarily due to slight changes in position of the bite bar mount between imaging sessions; a slight change in position between two sessions (a year apart) of a longitudinal study²⁵ gave rise to noticeable image scaling (the change in magnification is 2.2%/mm from a pivot conjugate plane). The effects of these motions are considerably more subtle and will be discussed in Sec. 3.

3. 3D B-Scan Registration

Section 3.1 describes our 3D B-scan registration algorithm that uses a coarse-to-fine registration approach with a combination of phase correlation and normalized cross-correlation (NCC) to improve processing speed and accuracy.⁷ Section A in Supplementary Material summarizes our GPU-based acceleration techniques. Section 3.2 describes the experimental procedure used to evaluate the performance of our approach, and the experimental results.

3.1.

Algorithm

We assume that fast B-scans are rigid, i.e., that they were acquired in a short enough time (0.3 to 1 ms) to contain only negligible motion artifacts. Thus, for each fast B-scan in a target volume we can determine its offset $(Δ x, Δ y, Δ z)$ with respect to the best matching fast B-scan in a reference volume by 3D cross-correlation of individual fast B-scans between reference and target volumes. Figure 3 shows the concept of this 3D B-scan registration method.

Fig. 3

Example of 3D registering a single fast B-scan of an AO-OCT retinal volume to the best-matching fast B-scan in a reference volume. This figure shows the slow B-scan projection (top row) and en face projection (bottom row) images of the target volume (left column), reference volume (middle column), and registered target volume (right column). The red line in each image indicates the location in the target image of the fast B-scan to be registered. The green line indicates the location of the best-matching fast B-scan in the reference volume found by our 3D B-scan registration algorithm, displaced from the target location by $(Δ x, Δ y, Δ z)$ . Target B-scans are shifted by this amount in the right column, registering the B-scans to the reference volume but creating gaps.

We used a coarse-to-fine B-scan search to accelerate image registration. This strategy uses a faster but less accurate method to compute a coarse displacement estimate. This coarse estimate is then used to define a small search space for a subsequent fine estimate using a slower but more accurate technique. The overall 3D B-scan registration approach consists of five steps [Fig. 1(a)]:

1. Prepare AO-OCT volumes for 3D registration by filtering out noise, cropping volumes, and choosing a reference volume (Sec. 3.1.1).
2. Apply high-speed coarse registration with 3D phase-only correlation (POC)²⁹ to obtain an approximate prediction of target B-scan image displacement relative to the best-matching B-scan in the reference volume (Sec. 3.1.2).
3. Define the search space for the subsequent step of fine registration by upsampling the coarsely-predicted B-scan displacements from the preceding step using a linear interpolator (Sec. 3.1.3).
4. Perform fine registration with NCC³⁰ for accurate prediction of individual fast B-scan locations within the reduced search area (Sec. 3.1.4).
5. Iterate (2)–(4) to register other target volumes to the same reference volume and then filter the registration results and visualize the registered averaged high-quality AO-OCT volume (Sec. 3.1.5).

This correlation-based algorithm is well-suited to parallel processing, enabling all computations to be implemented in CUDA (Sec. A in Supplementary Material).

3.1.1.

Data preparation

AO-OCT volumes require six pre-processing operations prior to registration. These are (1) correction of the uneven spacing of the A-lines caused by the nonlinear scan pattern along the fast B-scan direction; (2) correction of tilt introduced by imperfect eye alignment; (3) reduction of noise to avoid spurious correlation peaks; (4) conversion of data to single precision to fit into GPU memory and minimize processing time; (5) cropping each volume; and (6) choosing a reference volume. We summarize these pre-processing steps below.

1. Resample A-lines along the fast B-scan direction: due to nonlinearity of the scan pattern, the spacing of A-lines along the fast ( $X$ ) axis direction is uneven. To correct, we resampled A-lines using a scan pattern measured on a model eye, a procedure that is sometimes referred to as “de-sinusoiding,” though we do not use a sinusoidal waveform.
2. Correct retinal tilt along the fast B-scan direction: varying displacement of the imaging beam axis relative to the center of the ocular pupil causes the images of retinal layers to exhibit differences in tilt between volumes (see Sec. 2); this variation affects both POC and NCC calculations. This tilt is numerically corrected to sub-pixel accuracy using linear interpolation.
3. Remove noise floor and convert data type: to reduce noise, volumes are normalized, scaled and converted to single precision using the following equation:
Eq. (1)
$V^{'} = {\begin{cases} 1 e 5 & (M \leq V) \\ \frac{V - N}{M - N} \times 1 e 5 & (N \leq V < M) \\ 0 & (V < N) \end{cases},$
where $V$ and $V^{'}$ are pixel values before and after normalization, $N$ is the noise offset estimated from regions containing only non-reflective retinal structures (e.g., the vitreous), $M$ is a constant determined by dynamic range and is typically set 35 to 40 dB above the noise offset to capture the strongest retinal reflections, and the scaling factor $1 e 5$ is determined to avoid floating-point underflow and overflow and is a safeguard for squaring numbers and other multiplications in subsequent processing steps. We compute the noise offset $N$ using the 95th percentile of pixel values in the vitreous for each volume.
4. Crop volume: The POC dictates that B-scans in each volume must have the same size. To maximize processing speed, we crop all volumes to the smallest size possible for a specific imaging protocol, which must be no larger than the maximum B-scan size ( $512 \times 512 pixels$ ) that can be processed by our GPGPU implementation.
5. Select reference volume: The reference volume is chosen to be the volume with the highest quality from a series of the same retinal area. Criteria for quality are minimal eye motion artifacts, high image clarity, and high signal-to-noise ratio (S/N). We developed a simple quantitative metric that encapsulates these criteria based on cross-correlation of the en face projections of adjacent fast B-scans. For each pair of adjacent lines, we store the peak correlation value $r$ and the offset $Δ x$ at which the peak occurs. A large eye motion results in a local shearing of the image, leading to systematic variation in $Δ x$ and possibly a decrease in $r$ ; the presence of noise in the image will decrease $r$ and increase the variability in $Δ x$ . We, therefore, define an image metric $m = \overset{´}{r} / σ_{r} σ_{Δ x}$ , where $\overset{´}{r}$ is the mean of the correlation peak heights, $σ_{r}$ is the standard deviation of peak heights, and $σ_{Δ x}$ is the standard deviation of lateral shifts. The volume with the largest $m$ is selected as the reference.

3.1.2.

Coarse registration using 3D phase correlation

We register multiple target volumes to a single reference volume using a coarse-to-fine procedure as shown in Fig. 4. For simplicity, we describe the algorithm to register a single target volume to a single reference volume. The coarse stage applies a quick 3D POC to a set of target sub-volumes (stack of contiguous fast B-scans) and finds the best-matching sub-volume in the reference volume for each target sub-volume (see Secs. 3.1.2.1 and 3.1.2.2).

Fig. 4

Illustration of coarse-to-fine registration. (top) Coarse registration treats each target sub-volume (stack of 16 contiguous fast B-scans) as a rigid body and finds its estimated displacement $(Δ X_{N}, Δ Y_{N}, Δ Z_{N})$ relative to the best-matching reference sub-volume, where $N$ denotes the sample index $(N = 1,2, \dots)$ . Yellow colored areas indicate 20 sampled sub-volumes. Green colored area indicates the search space. (right) To define the search space for the subsequence step of fine registration, we upsample $Δ Y_{n}$ using a linear interpolator to predict the displacement of each target B-scan $Δ {y^{'}}_{i}$ , where $i$ denotes the B-scan index. The pre-registration position of $i$ ’th B-scan is just $i$ , so the predicted position is described as ${y^{'}}_{i} = i - Δ y_{i}^{'}$ . We then set the search space for a given B-scan as $({y^{'}}_{i} \pm w_{y})$ , where $w_{y}$ is a search space size, as indicated by green colored area in the bottom reference image. (bottom) Fine registration finds the accurate estimate of each fast B-scan’s location $(Δ x, Δ y, Δ z)$ within the narrower search space obtained from the preceding step.

Coarse registration sampling

To coarsely estimate the non-rigid displacements between two volumes (mainly due to eye drift), we select several sub-volumes or samples in each target volume that are evenly spaced by $s$ and have a width of $2 d$ along the $y$ axis [see Fig. 4 (top left)]. Selection of optimal values for these parameters requires managing a number of tradeoffs. $s$ must be small enough to accurately represent temporal variation in motion, but decreasing it also increases processing time. Since it is implicitly assumed that no motion occurs within a sample, increasing $d$ runs a greater risk of violating this assumption; on the other hand, decreasing $d$ lowers the S/N of the estimated displacement.

Our B-scan acquisition time typically ranges from 0.3 to 1 ms. We set the sampling interval $s$ to between 5 and 30 fast B-scans, so our sample sub-volumes are temporally spaced $< 30 ms$ apart. This captures the high temporal frequency components of eye motions that are as short as 10 to 60 ms, which we have found through trial and error to be adequate for normal eye movements. Similarly, we set the sample width $2 d$ between 8 and 16 fast B-scans, so each sample sub-volume takes 2.4 to 15 ms. This sample width is sufficiently short to freeze most eye motion⁵ yet sufficiently long for robust POC. In the presence of large, rapid eye movements, on the other hand, we can increase the robustness of POC by decreasing the sampling interval $s$ and creating a partial overlap between sample sub-volumes (at the cost of increased processing time).

3D phase-only correlation

The POC algorithm is commonly used for image registration because it is simple and is known to be well-suited for GPU implementation.²⁹ We use the POC algorithm to determine the displacement of each target sub-volume relative to the best-matching reference sub-volume based on the location of the maximum value of the 3D POC [see Fig. 4 (top)].

Prior to the POC calculation, the reference volume is first zero-padded along the $z$ and $x$ axes to the nearest greater power of two and along the $y$ axis to the nearest greater multiple of the sample width $2 d$ because: (1) the POC considers an image to wrap around at the edges, so zero padding or a window function should be applied to eliminate boundary effects; (2) the FFT is faster if the size is a power of two. Then, each target sub-volume ( $a_{N}$ , where $N$ is the sample sub-volume index) is zero-padded to match the size of the entire zero-padded reference volume ( $b$ ), because the POC requires that both input volumes be the same size.

We compute the 3D Fourier transform of each padded target ( $a_{N}$ ) and reference ( $b$ ) volumes as shown in Eq. (2), denoted by $F (a_{N})$ and $F (b)$ , respectively. We calculate the cross-power spectrum by element-wise multiplication of $F (a_{N})$ with the complex conjugate of $F (b)$ and normalize the result (the normalization and division are also element-wise):

Eq. (2)

R_{N} = \frac{F (a_{N}) \circ F {(b)}^{*}}{| F (a_{N}) \circ F {(b)}^{*} |},

We compute the inverse Fourier transform $F^{- 1} (R_{N})$ to obtain a 3D cross-correlation map for each target sub-volume. The location of the maximum value in this map identifies the best-matching identically sized sub-volume within the reference volume and yields the 3D displacement of the target sub-volume $(Δ X_{N}, Δ Y_{N}, Δ Z_{N})$ relative to that reference sub-volume.

3.1.3.

Coarse prediction of individual target B-scan displacements

We use the estimated displacement $(Δ X_{N}, Δ Y_{N}, Δ Z_{N})$ of each target sub-volume relative to the matching reference sub-volume found in the preceding coarse registration step to narrow the search space for the following NCC-based fine registration step, as shown in Fig. 4 (bottom). Reduction of the NCC algorithm search space confers two benefits: (1) it decreases the total computation time; and (2) it increases registration accuracy by decreasing the probability of finding an incorrect match that has significant structural similarity to the target image.

To reduce the search space along the $y$ axis, we first exclude target sub-volumes for which the displacement $Δ Y_{N}$ could not be estimated or deviated significantly from the overall trend ( $| Δ Y_{N} - Δ Y_{N - 1} | > s$ or $| Δ Y_{N + 1} - Δ Y_{N} | > s$ ). We then up-sample $Δ Y_{N}$ using a linear interpolator to predict the displacement of individual target B-scans $Δ {y^{'}}_{i}$ , where $i$ denotes the index of a B-scan within the entire volume. Displacement of target B-scans at the edges of the volume is linearly extrapolated from the nearest pair of $Δ Y_{N} s$ . The pre-registration position of $i$ ’th B-scan is just $i$ , so its predicted position is now described as ${y^{'}}_{i} = i - Δ {y^{'}}_{i}$ . The $y$ axis search space for a given B-scan is defined as $({y^{'}}_{i} \pm w_{y})$ , where the search space half-width $w_{y}$ is typically set to 8: we empirically find that the coarse shift estimate is usually within $\pm 3$ B-scans of the final estimate, so setting $w_{y}$ to eight ensures that the search space contains the matching reference B-scan.

It would be possible to also shrink the search space along the $x$ and $z$ axes for a given coarsely registered target B-scan based on the region of overlap between it and its best-matching reference B-scan. However, by the nature of the FFT-based NCC computation shrinking the search size requires a concomitant decrease in the size of the input data; put simply, the NCC simultaneously computes the correlation for every possible offset along $x$ and $z$ axes allowed by the size of the input. Shrinking the search space thus has two drawbacks. First, the overlapping region depends on the estimated displacement for each B-scan, so each FFT-based NCC computation would require different amount of memory; frequent memory re-allocation slows processing. By reducing the search space only along the $y$ axis, we keep all B-scans the same size as the original, allowing memory to be allocated only once. Second, shrinking the search space along $x$ and $z$ axes would require discarding information. We choose to allow the NCC algorithm to use the full original amount of B-scan information to correct for any prediction failures from the prior coarse registration step.

3.1.4.

Fine registration with normalized cross-correlation

We next register individual fast B-scans of a target volume to the best-matching B-scans in a reference volume with pixel-level (not subpixel) precision using NCC³⁰ [see Fig. 4 (bottom)]. Each target B-scan is cross-correlated with all reference B-scans within the reduced search space $({y^{'}}_{i} \pm w_{y})$ obtained from the preceding step. That target B-scan’s final displacement estimate $(Δ x_{i}, Δ y_{i}, Δ z_{i})$ is the location of the maximum NCC coefficient. It is worth noting that the $y$ component of displacement is simply the difference between B-scan indices in the target ( $i$ ) and reference ( $j$ ) volumes: $Δ y_{i} = i - j$ .

3.1.5.

Registered volume construction

After iterating the coarse-to-fine procedure described above to register other target volumes to the same reference volume, we estimate the individual fast B-scan displacement for each target volume:

Eq. (3)

(Δ x_{i, T}, Δ y_{i, T}, Δ z_{i, T}),

where

T

is the target volume index. We then construct registered volumes by displacing each fast B-scan by this calculated displacement. During this process, we detect and exclude B-scans that have statistically outlying displacements or low correlation coefficients. Spurious peaks in the NCC for a single B-scan give rise to abrupt change in displacement compared to its neighbors. We detect and exclude these outliers by computing the local standard deviation of displacements

σ_{Δ}

across all B-scans within a

\sim 20 - ms

window centered on the B-scan in question and that B-scan is excluded if

σ_{Δ}

exceeds 3 to 5 pixels or its NCC correlation coefficient is lower than 0.3. These values are determined empirically so as to exclude only outliers under normal conditions and can be easily adjusted in the case of particularly severe eye movements or other abnormal conditions in which there is a breaking point (threshold) where many B-scans are thrown out.

It is worth noting that non-linear distortions caused by variation in eye velocity can result in multiple target B-scans being displaced to the same location, whereas other locations might receive no target B-scans. If multiple B-scans are displaced to the same location, the last acquired B-scan is used. If a location receives no target B-scans, a gap may appear in the registered volume [see Fig. 3 (right column)]. These gaps vary randomly between registered volumes and can be eliminated by averaging volumes as previously demonstrated;¹³ in this study, we typically average more than 100 registered volumes for each location.

3.2.

Experimental evaluation

3.2.1.

Experimental design

We evaluated the performance of our 3D registration using the Indiana AO-OCT system (described by Liu et al.³¹ and Kocaoglu et al.³²). The system ran in two-camera mode to achieve an acquisition speed of 500 K A-lines/s. The light power incident on the subject’s cornea was at or below $430 μ W$ ( $λ_{center} = 790 nm$ ; $Δ λ_{FWHM} = 42 nm$ ) and within American National Standards Institute-defined safe limits for our experimental protocol. The same light was used for both AO-OCT imaging and wavefront sensing. We stabilized the subject’s eye and head position to the AO-OCT system using a dental impression mounted to a motorized $X Y Z$ translation stage. All procedures adhered to the tenets of the Declaration of Helsinki and were approved by the Indiana University Institutional Review Board. Written informed consent was obtained from all subjects prior to the experiment.

As stated in the Introduction, we have successfully used our method to register AO-OCT images acquired from more than 40 subjects. Here, we illustrate the method’s performance on one of our most challenging AO-OCT datasets with images acquired at many time points over a 24-h period on a normal subject with typical fixational ability. This extreme session length introduced a number of challenges to image registration: it required many look-ins by the subject, resulting in increased eye misalignment; it tired the subject, reducing his ability to fixate over the course of the session; it necessitated repeated use of Tropicamide that varied the subject’s refractive state and tear-film state between look-ins; and it included diurnal variation in the state of the eye. These complications collectively generated a wider range of registration errors than that encountered in a typical 1 to 2-h imaging session.

We imaged a 1.5° by 1.5° patch of retina 3° temporal to the fovea in a 29-year-old subject (left eye) who was free of ocular disease and had typical fixational ability. The experiment was designed to observe cone photoreceptor disc shedding,²³ which is a rare event (1/cone/day); we acquired three to four videos every 30 min throughout the evening of one day from 7:30 to 10:00 PM and throughout the following day from 6:00 AM to 6:00 PM. Each 5-s video contained 12 volumes acquired at $2.4 vol / s$ and spaced 0.41 s apart. B-scans were acquired at 1.1 KHz.

We initially registered 300 AO-OCT volumes to a reference volume containing minimal motion artifacts to evaluate registration performance. This set of volumes (34 videos acquired between 6 AM and 10:30 AM) represents a small subset of the total dataset; these are all the videos that could be registered with the 3D B-scan registration algorithm. As we will see in Sec. 4.2, registration to a global coordinate system will permit us to register the entire dataset (a total of 930 AO-OCT volumes in 110 videos.).

3.2.2.

Imaging quality metrics

We quantified the effectiveness of image registration using three metrics: relative power spectral contrast,³³ image sharpness ratio (ISR),³⁴ and mean squared error (MSE)³⁵ computed from en face projection images of cone photoreceptors.

Relative power spectral contrast is defined as the ratio of the power spectra of two images (in this case, an average of multiple targets en face projection versus a single reference en face projection) sampled at the cone spacing period. Larger relative power spectral contrast indicates better registration because correction of eye motion must, from a practical standpoint, increase the cone mosaic regularity, thereby creating more power at the fundamental frequency of the cone mosaic.

ISR is the reduction in image sharpness in the averaged target en face projection image $I_{avg}$ relative to the reference en face projection image $I_{ref}$ ; it is defined as

Eq. (4)

I S R = \sum_{p = 1}^{N} ‖ \nabla I_{avg} ‖ / \sum_{p = 1}^{N} ‖ \nabla I_{ref} ‖,

where

\nabla

denotes the spatial gradient operator (partial derivatives with respect to

x

and

y

),

‖ \cdot ‖

is magnitude, and

\sum_{p = 1}^{N} \dots

denotes summation over pixels common to both images; values range from 0 to 1, with larger values indicating better registration.

MSE is the mean squared difference in pixel value (intensity) between averaged and reference en face projection images, defined as

Eq. (5)

MSE = \sqrt{\frac{1}{N} \sum_{p = 1}^{N} {(I_{avg} - I_{ref})}^{2}} .

Smaller MSE values indicate better registration.

We varied the number of averaged target images and examined the effect of that number on each of these three metrics. We compared: (1) the average registered-target en face projection images versus a single reference en face projection image that we judged to be of high image quality; (2) the average unregistered-target en face projection images versus the same reference image; and (3) simulated en face images composed of random pixel-correlated white noise (uniformly distributed both spatially and temporally) versus the reference. A similar comparison has been reported by Ferguson et al.¹⁵

3.2.3.

Results

Figure 5 shows en face projection images of: Fig. 5(a) the selected reference volume; Fig. 5(b) average of 300 unregistered target volumes; and Fig. 5(c) average of 300 registered target volumes. The substantial difference in visibility of the cone mosaic between Figs. 5(b) and 5(c) show the power of our registration algorithm. This point is demonstrated quantitatively by plotting the circumferential average of the two-dimensional power spectrum of each image in Figure 5(d); power at the cone mosaic fundamental frequency is $\sim 100 \times$ greater after registration and approaches that of the single reference image. The residual degradation relative to the reference is likely caused by residual (sub-pixel) displacement errors; these cannot be removed by our registration algorithm as it registers images to the nearest whole pixel.

Fig. 5

Comparison of en face projection images: (a) reference volume; (b) average of 300 unregistered target volumes; and (c) average of 300 registered target volumes. (d) Circumferential average of two-dimensional power spectrum for the reference image (black curve), average of 300 unregistered target images (blue curve), and average of 300 registered target images (red curve). Fast B-scans are oriented vertically in the images.

Figure 6 shows our image quality metrics (relative power spectral contrast, ISR, and MSE) applied to the same dataset, each plotted as a function of the number of volumes averaged. As a worst-case baseline, we also plot the same metrics for simulated en face images composed of random pixel-correlated white noise (gray circles). Registered en face images have much higher relative contrast and ISR values, and much lower MSE values (better image quality) than unregistered images. Collectively, our results using relative contrast, ISR, and MSE are consistent with our visual observations and show that our 3D B-scan registration algorithm substantially enhances en face image quality by correcting motion artifacts and permitting averaging of multiple volumes.

Fig. 6

Effectiveness of the 3D B-scan registration algorithm for the same dataset used in Fig. 4. Image quality metrics, (a) relative power spectral contrast; (b) ISR; and (c) MSE, plotted as a function of number of volumes averaged.

We also evaluated registration speed for one target volume from the same dataset. Our algorithm runs much more quickly when implemented in CUDA on a GPU than in a high-level language on a CPU as we have done previously: an equivalent MATLAB implementation on a CPU (Intel Xeon CPU E5-1620 v3) registered a $450 \times 450 \times 300 - pixel$ volume in 900 s, whereas our CUDA implementation on an NVIDIA GeForce GTX1080 took only 6.1 s, a $150 \times$ speed increase. Our previous 2D strip-wise registration algorithm (a combination of retinal layer segmentation and strip-wise registration) implemented in MATLAB on a CPU also took several minutes to register the same volume. The lack of a low-level language CPU-based implementation of course precludes a legitimate comparison between CPU and GPU processing. Nevertheless, we expect these speed differences to be primarily attributed to differences in the processing units.

4. Global A-Line Registration

Section 4.1 presents our global A-line registration algorithm for correcting torsional eye movements and image scaling and generating global motion-free coordinates. Section 4.2 describes the experimental procedure used to evaluate the performance of our approach followed by the experimental results. Sections 4.3 and 4.4 provide two applications of the algorithm to extend field of view and depth of focus, respectively.