1.

Introduction

Optical coherence tomography (OCT) uses interferometric analysis of low-coherence light to provide depth-resolved micrometer-resolution images of tissue, such as retina and optic nerve head.¹ Clinically, both cross-sectional (B-frame) and three-dimensional (3-D) structural OCT are useful for detecting and monitoring anatomic alterations. With the increase of OCT imaging speed,^2,3 an OCT angiography (OCTA) technique was developed.^4–14 Multiple B-frames can be taken on the same position, and the changes of OCT reflectance properties can be measured to differentiate vasculature from static tissues. These B-frames can also be averaged to generate structural OCT images.

Multiple B-frames taken in the slow transverse direction can create a volumetric scan. A single volumetric scan usually takes 3 to 5 s, during which involuntary eye movements, such as microsaccades and saccades, can introduce motion artifacts.^15–17 They appear as discontinuities on the en face projection of OCT volume and bright lines (saturated signals) on the en face OCTA. These large and rapid eye motions on the transverse direction can be mitigated by adding retinal tracking modality, such as scanning laser ophthalmoscope and infrared fundus video camera that can detect image shifts, pause the OCT scan, and then resume the scan at the proper position.^18–23 Eye tracking strategy can reduce the discontinuities on en face images, but it cannot correct mild transverse or any axial motion. The remaining motion artifacts show as bright lines on en face OCTA. Therefore, additional processing is required to eliminate these lines.

Two-dimensional (2-D) registration algorithms have been developed to reduce motion in en face OCTA—one using the orthogonal method (registration of $x$-fast and $y$-fast en face angiograms)²⁴ and one using the parallel method (registration of two $y$-fast en face angiograms).²⁵ These methods are simple and fast. However, they cannot correct motion on volumetric data, and the distortion of angiographic and structural B-frames on slow axis cannot be removed. Correction of motion artifacts on the volumetric data not only provides clinically useful cross-sectional images along slow axis but also improves the segmentation efficiency on retinal layers. In addition, a motion-corrected OCT volume can improve volume rendering of the retinal vasculature, the optic disc, and the fovea.

One successful algorithm for motion correction on OCT volume is the 3-D orthogonal method.²⁶ It combines the information in $x$-fast and $y$-fast volumetric raster scans and registered A-scans individually by finding a 3-D motion field that minimizes an energy function. This algorithm has been adopted in the commercial Avanti OCT system (Optovue, Inc., Fremont, California) to register angiographic and structural OCT data by detecting the motions calculated from 3-D structural OCT. Although this technology greatly improves the merged scan quality, residual lines can still persist on OCTA.

In this study, we propose a 3-D registration algorithm for OCT volume based on both angiographic and structural OCT information. The vasculature on en face retinal OCTA is used for transverse registration, and the inner limiting membrane (ILM) boundary on structural OCT is taken as the reference for axial registration. Building on our previous work on parallel registration suitable for the ultrahigh-speed OCT (above 200 kHz), this algorithm will provide a comprehensive registration system for automated volume rebuilding in wide-field angiographic and structural OCT.

2.

Method

2.3.

Volumetric Registration Algorithm

This proposed method is composed of four main steps (Fig. 1). First, the transverse motion artifacts on en face retinal OCTA are corrected by an affine registration algorithm on the parallel microsaccade-free strips. This step corrects the motion artifacts between B-frames caused by translational, rotational, and focusing differences (Fig. 2).²⁵ Each A-line of reflectance and decorrelation signal in the volumetric OCT is translated based on the transformation matrix of the corresponding strip. Then, the ILM boundaries on the new transformed slow axis B-frames at the same position in each volume are averaged and smoothed to recover a motion-free ILM boundary. Based on this new boundary, each pair of A-line was initially registered along the axial direction at each B-frame. The final axial position of each A-line was fine-tuned by correlation of reflectance of each pair of A-line. Next, a registration based on cross-sectional OCTA was applied between each pair of fast axis B-frame to refine the final transverse position of each A-line. Lastly, two volumes of structural OCT and OCTA were averaged and merged, respectively, to produce the final motion-corrected volumes.

Fig. 1

The flow chart of the 3-D registration and volume rebuilding algorithm.

Fig. 2

The schematic diagram of the between-frame registration process between two $y$-fast volumetric scans. (a) First en face OCTA. The white vertical line artifacts are caused by microsaccadic bulk motion. (b) Second en face OCTA. (c, d) The microsaccadic line artifacts are detected and removed, leaving intermicrosaccadic strips for registration. (e) The merged en face OCTA after affine registration of the strips.

2.3.1.

Between-frame registration

We first detected and removed the large motion artifacts shown as bright lines on en face OCTA. The algorithm detected these bright line artifacts as projected B-frames with summed flow signal above the threshold, set as 2.3 standard deviations above the mean. Then en face angiogram was divided into microsaccade-free strips at the detected motion lines [Figs. 2(c) and 2(d)]. We initiated the registration process on two strips, one from each scan that contains the largest overlap. Using the larger of the two strips as the reference strip and the other as the moving strip, we applied the affine registration algorithm. It utilizes a gradient descent to search for a similarity transformation, which allows translation, scaling, rotation, and skewing within the transverse ($x\text{-}y$) plane. The goal is to find a transformation of the moving strip that minimizes the squared difference of the reference vasculature and transformed moving vasculature

Eq. (1)

$$C_{\text{between}}(t)=\sum _{x,y}{\{V_r(x,y)-V_m[T(x,y;t)]\}}^2,$$

where $(x,y)$ is the pixel coordinate, $t$ is a vector of transformation parameters, $V_r$ $(x,y)$ is the reference vasculature, $V_m$ $(x,y)$ is the moving vasculature, and $T$ is the transformation matrix. This minimization can be solved iteratively. The new registered and merged strip was then defined as the new reference strip and the strip with the largest overlap with the new reference strip as the new moving strip. The process was repeated until all strips were registered and merged into one comprehensive en face OCTA image [Fig. 2(e)].

This step provides the transverse transform matrix. Each A-line of the first and second volumetric structural OCT and OCTA was translated, and the between-frame transformation was completed.

2.3.2.

Rebuilding of the inner limiting membrane boundary

Next, the ILM boundary was rebuilt by merging and smoothing ILM boundaries on two transformed volumetric structural OCT.²⁸ This rebuilding process was performed in the slow axis and included three steps: first, a pair of slow axis ILM boundaries from each volume was adjusted to the closest position [Figs. 3(a) and 3(b)]. In this step, the ILM boundary from the reference volume is the reference boundary $B_r$, and the ILM boundary from the moving volume is the moving boundary $B_{m_0}$. The newly moved boundary $B_m$ was obtained when the moving boundary adjusted to the closest position to the reference boundary based on the below equation:

Eq. (2)

$$\begin{gathered} B_m(x)=B_{m_0}(x)+tx=\mathrm{1,2},\dots ,L \\ E=\sum _{i=1}^L|B_r(i)-B_m(i)|B_r(i)>0,B_m(i)>0 \\ {t}^{*}=\arg \min (E|t)t\in [-\mathrm{100,100}], \end{gathered}$$

where $t$ is the displacement of the $B_{m_0}(x)$, $E$ is the cost function, $L$ is the number of the A-line in corresponding slow axis B-frame, and the optimized displacement $t^{*}$ is selected when the cost function $E$ achieves the smallest value. The height profiles of reference boundary and the moved boundary were averaged to generate the initial merged ILM boundary $B_{\mathrm{ave}}$ [Fig. 3(c)]

Eq. (3)

$$B_{\mathrm{ave}}=(B_r+B_m)/2.$$

Fig. 3

The rebuilding process of the ILM boundary height profile along the $x$ (slow)-axis. (a) The original reference (red) and moving (yellow) boundaries. (b) The transformed moving boundary (yellow) is translated to match the reference boundary (red) in average height. (c) The average height profile (green) of the reference (red) and transformed moving (yellow) boundary. (d) The average height profile is smoothed using a moving average filter to obtain the final ILM boundary (green).

Then, the final rebuilt ILM boundary of this B-frame was obtained after a smooth process on $B_{\mathrm{ave}}$

Eq. (4)

$$B_{\text{rebuilt}}=\{\begin{array}{ll}{B}_{\text{rebuilt}}(n)=B_{\mathrm{ave}}(1)+\cdots +B_{\mathrm{ave}}(n)+\cdots +B_{\mathrm{ave}}(n+n-1)& n\in [1,\mathrm{sp}]\\ B_{\text{rebuilt}}(n)=B_{\mathrm{ave}}(n-\mathrm{sp})+\cdots +B_{\mathrm{ave}}(n)+\cdots +B_{\mathrm{ave}}(n+\mathrm{sp})& n\in [\mathrm{sp}+1,L-\mathrm{sp}]\\ B_{\text{rebuilt}}(n)=B_{\mathrm{ave}}[n-(L-n)]+\cdots +B_{\mathrm{ave}}(n)+\cdots +B_{\mathrm{ave}}(L)& n\in [L-\mathrm{sp}+1,L]\end{array}\mathrm{sp}=\lfloor L/20\rfloor ,$$

where sp is the smooth parameter and $\lfloor \rfloor$ means round down [Fig. 3(d)].

2.3.3.

Axial registration

After the rebuilt, the ILM boundary was obtained, and each A-line of two OCT and OCTA volumes was translated based on the displacement between $B_{\text{rebuilt}}$ and corresponding ILM boundary

Eq. (5)

$$\begin{gathered} D_r(x,y)=B_{\text{rebuilt}}(x,y)-B_r(x,y) \\ {D}_m(x,y)=B_{\text{rebuilt}}(x,y)-B_m(x,y), \end{gathered}$$

where $D_r(x,y)$ and $D_m(x,y)$ are the displacement of each A-line in reference and moving volumes, respectively, and $(x,y)$ is the coordinate of each A-line. Each A-line was then translated based on the corresponding displacements. To reduce the registration error caused by ILM segmentation, another registration between these two A-lines was applied to fine-tune the axial final positions. Before the registration, the reflectance of two A-lines [Figs. 4(a) and 4(b)] was first normalized to [0, 1] and filtered by a $1\times N_z$ pixel Gaussian filter with 5 pixel standard deviation [Fig. 4(c)]. The registration was optimized by a cost function $C(p)$

Eq. (6)

$$C(p)=\frac{1}{N_z}\sum _{z=\mathrm{1,2},\dots ,N_z}{[A_r(z)-A_m(z+p)]}^{2}p=[-\mathrm{10,10}],$$

where $C(p)$ is the cost function of the registration, $A_r(z)$ is the normalized reflectance of the reference A-line, $A_m(z+p)$ is the normalized reflectance of the moving A-line after a displacement $p$, and $N_z$ is the number of pixel in each A-line. The optimized displacement $P^{*}$ was selected when the cost function $C(p)$ achieved the smallest value. In the final axial registration step, $A_r(z)$ and $A_m(z)$ were translated to $A_r[z-(p^{*}-\lfloor p^{*}/2\rfloor )]$ and, respectively, [Fig. 4(d)].

Fig. 4

The axial registration process. (a) The reference A-line of the normalized reflectance signal (red). (b) The moving A-line of the normalized reflectance signal (green). (c) The reference and moving A-line low-pass filtered by the Gaussian function. (d) The moving and reference A-line were registered by minimizing their sum squared difference and shifted to their mean position.

2.3.4.

Within-frame registration

To further improve the transverse registration accuracy, a within-frame registration based on the cross-sectional decorrelation signal profile was performed.

The decorrelation signal on each pair of B-frames along $y$ (fast)-axis was registered by an affine registration (Fig. 5). Since there were just some small distortions between two initially registered volumes, the ranges of translation, scale, rotation, and skew were limited to $|\text{translation}|<5$, $|\text{scale}-100\%|<5\%$, $|\text{rotation}|<2\mathrm{deg}$, and $|\text{skew}-90\mathrm{deg}|<3\mathrm{deg}$.

Fig. 5

An example of the within-frame registration between a pair of OCTA B-frames. (a) The reference B-frame (flow signal in red). (b) The moving B-frame (green). (c) The moving B-frame is overlaid on the reference B-frame without registration. (d) The two B-frames after affine registration to minimize sum squared difference. The enlarged regions are outlined in blue. (e) Enlargement of panel c showing mismatch in the position of two large vessels (blue arrows). (f) Enlargement of panel d showing registration of the two mismatched vessels.

After the within-frame registration, the final merged volumetric scans were obtained by averaging two registered OCT and OCTA together. A postprocessing, including Gabor filter and multiscale vessel enhancement filter, was performed to improve the image quality of the en face OCTA.^29,30

3.

Evaluations

The automated volumetric registration of two scans took an average of 7 min. The test was performed on a workstation with Inter(R) Xeon(R) CPU E3-1226 v3 at 3.30 GHz and 16.0 GB RAM using MATLAB^® 2014b (Mathworks, Natick, Massachusetts).

Figure 6 shows the registration results of two $y$-fast scans on central macula. The volumetric rendering of the scan volumes (first row) before and after motion correction shows the surface ripples due to axial motion on the original $y$-fast volume [Figs. 6(a) and 6(b)]. The correction recovers the smoothness of volume surface [Fig. 6(c)]. B-frames on slow axis, seen on the second row, also show the performance of axial registration. The microsaccades along axial direction distort the retinal and choroidal layers in the original B-frames [Figs. 6(d) and 6(e)]. Axial registration recovers the smooth boundaries of retina and choroid, including the outer boundary of the choroid [Fig. 6(f)]. No motion artifact is shown on B-frames along the $y$ (fast)-axis [Figs. 6(g) and 6(h)] due to the high A-line rate on the vertical priority. Our results show no registration error is introduced on fast axis in the merged B-frame [Fig. 6(i)], although the axial registration is performed on the slow axis. In addition, the transverse motion artifacts and the vasculature incoherence are all corrected on OCTA [Figs. 6(j)–6(l)].

Fig. 6

3-D registration of two $y$-fast volumetric OCT and OCTA ($6\times 10\times 7\mathrm{mm}$) of macular region on a normal eye: (a) first volumetric OCT, (b) second volumetric OCT, (c) merged volumetric OCT after the 3-D registration of a and b, (d) the B-frame along $x$ (slow)-axis in first volume (red line in j), (e) the B-frame along $x$-axis in second volume (red line in k), (f) the B-frame along $x$-axis in the merged volume (red line in l), (g) the B-frame along $y$ (fast)-axis in first volume (yellow line in j), (h) the B-frame along $y$-axis in second volume (yellow line in k), (i) the B-frame along $y$-axis in the merged volume (yellow line in l), (j) first en face OCTA, (k) second en face OCTA, and (l) the merged en face OCTA after 3-D registration of two volumes.

Figure 7 shows that the registered and merged structural B-frame on slow axis crossing fovea mimics the B-frame acquired by $x$-fast scan on the same position of the same eye, indicating that the axial registration based on ILM boundary rebuilding on slow axis resembles the anatomic retinal contour.

Fig. 7

(a) A comparison of the slow axis B-frame of the merged OCT and (b) the $x$-fast B-frame crossing the fovea, which represents the true contour of the retina.

To further demonstrate the stability and adaptability of this method, we also show the registration results on a scan that includes the optic nerve head, which is more complex than the fovea. The registration performance can be globally reviewed by volumetric rendering of the scan volumes before and after registration [Figs. 8(a)–8(c)]. The registered volume shows the even and regular appearance of the optic nerve head and surrounding retina structures. The distorted cross-sectional B-frames on slow axis in the original volumes were corrected in the registered volume [Figs. 8(d)–8(f)]. The image contrast of retinal layers and its microstructures are enhanced. Similar to the original fast axis B-frames, the registered and merged B-frame shows no translation artifacts [Figs. 8(g)–8(i)]. En face OCTA after registration shows more continuous and distinctive vasculature than that in the nonregistered volume. Again, the retinal layer profiles and the regular anatomic microstructures of disc shown in the merged B-frame are highly correlated to the corresponding $x$-fast B-frame on the same position (Fig. 9).

Fig. 8

3-D registration of two $y$-fast volumetric OCT and OCTA ($6\times 10\times 7\mathrm{mm}$) of disc region on a normal eye: (a) first volumetric OCT, (b) second volumetric OCT, (c) merged volumetric OCT after the 3-D registration of a and b, (d) the B-frame along $x$ (slow)-axis in first volume (red line in j), (e) the B-frame along $x$-axis in second volume (red line in k), (f) the B-frame along $x$-axis in the merged volume (red line in l), (g) the B-frame along $y$ (fast)-axis in first volume (yellow line in j), (h) the B-frame along $y$-axis in second volume (yellow line in k), (i) the B-frame along $y$-axis in the merged volume (yellow line in l), (j) first en face OCTA, (k) second en face OCTA, and (l) the merged en face OCTA after 3-D registration of two volumes.

Fig. 9

(a) A comparison of the slow axis of the merged OCT and (b) the $x$-fast B-frame across the optic nerve head, which represents the true contour of the retina.

To quantitatively evaluate the registration accuracy, the average pixel distance (APD) between the ILM boundary of fast axis and the slow axis was calculated using the below equation:

Table 1

The comparison of APD to fast axis B-frames on slow axis B-frames before and after registration.

To quantitatively measure the improvement of contrast on en face OCTA, the root-mean-square (RMS) contrast³¹ was used and calculated using the below equation:

These quantitative analyses demonstrate significant improvement of the retinal structural contour (Table 1) on the structural OCT and the RMS contrast (Table 2) on the en face OCTA.

Table 2

The comparison of RMS contrast on en face OCTA before and after registration.

4.

Discussion

We have developed a 3-D registration method that removes motion artifacts and merges two volumetric angiographic and structural OCT volumes in the same transverse priority. This is well suited for OCTA using ultrahigh-speed OCT. The algorithm effectively eliminates axial and transverse motion artifacts, resulting in smooth volume rendering and high-contrast angiography with excellent vessel continuity.

The main interest of this study is the use of OCTA and structural OCT information separately for transverse and axial registration. To the best of our knowledge, no previous work has used a motion indicator on OCTA to guide volumetric registration. Compared to shadow-graphic pattern of large vessels on en face structural OCT, the detailed vasculature on OCT angiograms provides a better reference for transverse registration. Another advantage of using en face OCTA is that the microsaccade motions have been naturally shown as the bright motion artifact, eliminating the need to compute the cross-correlation between contiguous lines (projected B-frames) to detect this type of motion on en face structural OCT.

Unlike the previously reported orthogonal 3-D registration methods²⁶ using full A-scan profiles, our method performs transverse and axial registration separately to achieve 3-D registration. This strategy allows a more efficient processing, especially on the large scan volume with large number of A-lines.

The transverse registration method in this study differs from our previous work²⁵ in that we did not apply free-form deformation, which is a type of local nonrigid registration, because the translation of A-lines based on this deformation requires much computation time and memory space. Although the registration on the en face direction may not be as accurate as the work reported on 2-D OCTA,²⁵ by adding the within-frame affine registration based on cross-sectional flow profile between each pair of $y$-fast B-frames, we were able to further improve the registration accuracy in 3-D.

To correct the motion in axial direction, we needed to rebuild the anatomic structure on the slow axis. Although we do not have a reference scan with different scan priority, we have demonstrated that the convergence of two B-frames on the slow axis from different OCT volumes can effectively recover a smooth retinal surface. However, the notable limitation of this method is that the registered retinal surface is still a few pixels different from the true retinal contour. This deficit may not be clinically significant because the recovered contour allows for accurate segmentation of retinal sublayers and generation of high-quality en face structural and angiographic images. Moreover, in this step, it is critical to accurately trace ILM boundaries. The large gradient of reflectance between the vitreous and the ILM usually allows accurate segmentation of this boundary. In our study, we applied the directional graph search method²⁸ to further improve the segmentation efficacy.

Finally, this registration method is complementary to a real-time eye tracking function. Even with the acquisition of two volumetric scans, there is still a small probability of gaps occurring in the same place at both volumes, causing registration error and stripe defect in the merged volume. We have recently shown that we can reduce this failure probability to near zero by combining tracking assisted scanning with registration.¹⁸ Therefore, the registration algorithms developed by us and other groups would be useful in the correction of residual tracking error.

5.

Conclusion

We propose an automated 3-D registration algorithm to remove motion artifacts and merging two volumetric OCT and OCTA scanned by an ultrahigh-speed swept-source OCT. We removed the motion artifacts in the “transverse” direction using the between-frame registration and corrected the “axial” location based on the ILM boundary of the retina. Two volumetric scans were merged after the axial registration and within-frame registration. We demonstrated our algorithm on OCT and OCTA volume including the macula and the disc in healthy volunteers. By improving the volume rendering of the retina and the quality of OCTA, this method may improve the utility of OCT-based ocular imaging.