1. Introduction

Optical Coherence Tomography (OCT) is a non-invasive optical imaging technique [1], which has been especially successful in ophthalmic applications [2, 3]. By acquiring a series of 1D profile scans (known as A-scans), cross-sectional images of the retina (known as B-scans) are obtained. Composing successive B-scans yields a 3D image of the retina. Typically, B-scans are generated by raster-scanning the retina in one direction, known as the fast-scan direction (Figs. 1(a) & 1(b)). In the 3D volume of the retina, the direction perpendicular to the fast-scan direction is referred to as slow-scan direction. The A-scan acquisition rate in modern spectral domain (SD-OCT) or swept source (SS-OCT) OCT imaging systems [4, 5] is fast enough to assume B-scans are acquired free of motion artifacts [6]. However, during volume acquisition, patient motion artifacts are often observed along the slow-scan direction and may appear as discontinuities in the en face projection of the volume (Fig. 1(b)). These motion artifacts may be due to eye motion such as drifts, saccades or micro saccades, or bulk motion such as heartbeat or respiration [6–8].

 

Fig. 1 Motion artifacts in volumetric OCT imaging of the human retina. (a) Example of a volumetric OCT scan. The B-scan in the fast scan direction shows smooth retinal layers without visible motion artifacts (X-fast scan). The wavy retinal layers in the slow (Y) scan direction is the result of uncompensated motion of the subject’s eye. (b) Example summed voxel projection (SVP) image (Z-direction) of the original volume. The red arrow indicates the clear location of a saccade. (c) SVP after motion compensation. The dark gap region in this image, an artifact of saccades, corresponds to a region not scanned by OCT.

Download Full Size | PDF

Motion artifacts become increasingly problematic with longer duration OCT acquisitions such as those needed for wide field-of-view (FOV) or densely sampled volumetric imaging. Correction of motion artifacts within OCT volumes is important not just for aesthetics but is also critical to the accuracy of OCT based clinical applications such as angiography [9–12] and ocular biometry [13,14].

Several strategies to mitigate motion artifacts in OCT scans have been proposed to ensure data integrity. Referring to the lateral directions as X and Y and to the axial direction as Z, one commonly used strategy is to align the X–Y plane of an OCT volume to a complimentary image acquired with another modality. The use of simultaneously acquired OCT volumes and scanning laser ophthalmoscopy (SLO) images is one such example, where SLO en face retinal images are acquired at a rate that is often the B-scan rate of the OCT system, thus limiting artifacts within the en face image [15,16]. Other methods include finding shifts in the data by maximizing the cross-correlation between B-scans [17], a target tracking approach [18], and assuming a local-symmetry for the shape of the retina [19]. Hardware-based strategies have also been proposed, such as using eye-tracking devices to compensate for the motion [20]. Another method proposed capturing several fast sparsely sampled scans, detecting motion artifacts and removing corresponding scans, and then fusing the artifact-free scans [6,21]. In practice, these methods are most suitable for specific applications due to the need for extra hardware, long image acquisition time, or limited targeted accuracy.

One successful strategy for correcting motion artifacts in OCT is to capture two orthogonally oriented volumetric scans containing uncorrelated motion and subsequently reconstructing a motion-free volume by combining information from both datasets [4,10,22–24]. In such datasets, volumes which are acquired with fast scans oriented in the X-direction are often referred to as X-fast volumes while volumes acquired with fast scans orthogonally oriented in the Y-direction are Y-fast volumes. In [23, 24], the registration is performed by finding a 3D motion field for the A-scans that minimizes an energy function with two terms: one that favors matching of the A-scans between the orthogonally-scanned volumes and one that penalizes irregularities in the motion field, producing smooth motions along the fast scan directions. Recently other redundant scanning patterns have been proposed. Hong et al. [25] utilized a Lissajous scan pattern, so that all locations in the retina are repeatedly acquired, but no discontinuities in the scanning pattern are produced. Zang et al. [12] utilized two OCT volumes with the same fast scan direction to remove en face motion artifacts by registering saccade-free parallel strips.

For registration in the axial direction, one common approach is to segment the retina [26,27] and to fit the segmentations of two orthogonal volumes [22, 28–30]. The retinal boundary (boundary between vitreous and nerve fiber layer) can be segmented as a curve in each B-scan independently [27] or as a surface in the entire volume [26]. We refer the reader to [31] and [32] for extensive reviews on OCT motion correction techniques, denoising and segmentation.

While many variations of orthogonal volume registration have been proposed, they have primarily been limited to correcting narrow field-of-view retinal OCT volumes and have not directly addressed the challenges of wide FOV retinal volumes. The large size of the volumes used in wide FOV OCT can make methods based on A-scan optimization intractable, unless the volumes are subsampled [23,24]. Additionally, in large OCT volumes, saccadic eye motion is more likely to occur due to the longer acquisition times. This can cause relatively large portions of the retina to be missed during acquisitions (Fig. 1(c)), which results in inaccurate registration if not taken into account. By detecting the temporal location of saccades within a volume acquisition, a rigid motion model may be applied to the target between saccades.

In this work, we propose a comprehensive motion correction algorithm, suitable for wide FOV retinal OCT volumes with strong motion artifacts and multiple saccades. Following the overall approach utilized in previous orthogonal scanning methods, when correcting the slow scan axis motion in one volume, the corresponding fast axis B-scans in the other volume are taken as the reference.

Our method is based on two key ideas. One is to exploit the segmentation of the retinal boundary to guide the registration. This segmentation is used twice: to create a projected image of the inner retina that enables lateral plane registration, and as a retinal depth measure to find the axial registration. The second key idea is to include during registration the detection and correction of saccades. This is especially important in the case of wide FOV OCT volumes, since the occurrence of saccades is higher than on narrow FOV volumes. The proposed approach offers a complete robust system for 3D reconstruction for wide FOV OCT scans.

The rest of this article is organized as follows: In Section 2, we describe our approach for X-fast/Y-fast wide FOV OCT registration and motion correction. In Section 3, we show registration results of our proposed method. In Section 4 we present a discussion and conclusions.

The source code corresponding to the proposed algorithm is available online at http://people.duke.edu/~sf59/lezama_BOE_2016.htm.

2. Methods

2.1. Algorithm overview

The proposed algorithm is divided into three main parts (Fig. 2). First, we segment the vitreous/retina boundary in each individual B-scan. We use the segmented boundary to obtain the approximate location of the inner retinal layers, and produce an enhanced summed voxel projection (SVP) en face image of the retina. This SVP image is used for lateral registration in three steps: 1) Detecting saccadic eye motion and splitting the volumes at each saccade location; 2) Stitching the inter-saccadic intervals using orthogonal scans as reference and a locally affine model; and 3) An orientation-aware optical flow algorithm for dense local registration, that results in motion only in the fast-scan direction.

 

Fig. 2 Block diagram of the proposed method. The method is divided in three main blocks. A: retinal surface segmentation, performed in each B-scan separately. B: lateral registration, including saccade detection and correction, and local dense registration using an orientation-aware version of optical flow. C: axial registration, obtained by solving a linear model on the segmented surfaces, equivalent to applying a tilt and an offset to each B-scan.

Download Full Size | PDF

Once the affine transformations and local displacements have been found for each saccade-free region of the SVP, they are applied to the A-scans, or equivalently, to each X–Y slice of the input volumes. In the last stage of the algorithm, the segmented retinal surfaces of each volume are used to correct for motion in the axial direction, which translates into vertical shifts of the A-scans. As a result, we obtain two registered and (generally) artifact-free volumes.

Before providing details on each one of these steps, let us mention how our proposed method improves upon previous approaches:

2.2. Retinal boundary segmentation

The first building block of our method is the segmentation of the nerve fiber layer (NFL)-vitreous boundary. The segmented boundary is then used in the two main steps of the algorithm. For the lateral registration step, we use it to produce SVP images of the approximated inner retinal layers that are used for image registration. For the axial registration step, we match the retinal surfaces obtained from the two volumes by fitting a linear depth distortion model.

We segment the retinal surface on each B-scan image using a method based on finding the shortest-path to traverse the image from left to right (Fig. 3), where the cost of traversing one pixel is given by the magnitude of its gradient. The segmentation procedure consists of the following steps, which we will describe in detail in the rest of this section:

 

Fig. 3 Retinal boundary segmentation. (a) B-scan image (denoised with non-local means [34]). (b) Background pixels binary mask (green=1, blue=0). (c) Normalized square of the gradient magnitude (red is highest, blue is lowest). (d) Resulting segmentation. Red: retinal boundary. Blue: IS/OS layer.

Download Full Size | PDF

2.2.1. Image filtering and background removal

Given a B-scan image b of H pixels height and W pixels width, we first filter the image using the non-local means algorithm [34] (Fig. 3(a)). We normalize the filtered image to have values between 0 and 1. Next, we identify the pixels that belong to the background as those whose value is close to a reference background pixel value, which we find by taking the average pixel value in a region inside the vitreous. For the images used in this work, this region was chosen to be the span of the 10^th to 50^th rows of the B-scan. We create a binary mask indicating the background pixels, whose value is 1 where the difference between the pixel value and the reference background value is small (less than 0.1 in this work) and 0 otherwise (Fig. 3(b)).

2.2.2. Initial estimation of the retinal boundary

We compute the gradient of the filtered image using MATLAB’s imgradient function. We take the square of the gradient magnitude and normalize it so that the values are between 0 and 1. To eliminate noisy gradients from the background noise, we apply the background binary mask computed in the previous step (Fig. 3(b)). The segmentation is finally performed using a shortest path algorithm that cuts the images where the gradient is most prominent (Fig. 3(c)). Our shortest path algorithm is based on the Graylevel-weighted Geodesic Distance Transform [35]. This transform finds the graylevel-weighted distance from every point in the image to any target region in the image. We use MATLAB’s implementation of this function, graydist. We compute the transform twice: once using as target region a column appended to the left side of the image and once using a column appended to the right side of the image. Summing both image transforms yields a distance map where the value of each point is the sum of the weighted shortest paths to the left and right sides of the image. It can be shown [35] that the regional minimum (the set of pixels in an image that are equal to the global minimum of the image) of this map is the weighted shortest path between the left and right sides of the image. We obtain this regional minimum using MATLAB’s imregionalmin function.

2.2.3. Refinement of the initial segmentation

The obtained shortest path cut corresponds to either the vitreous/retina boundary or the IS/OS layer, or in some degenerate cases, a mix of the two. To identify the locations where the segmented boundary was the IS/OS, we look for a bright zone on top of the segmentation. We consider the average pixel values in a band-shaped region between above the segmentation (in this work 94μm to 188μm above it). If the difference between the average value and the reference background value described before is more than a given threshold (0.12 in this work), we mark that location as corresponding to the IS/OS boundary and not the retinal boundary.

To prevent the segmentation from picking the IS/OS boundary, we set the gradient to zero from the bottom of the image to the estimated IS/OS boundary plus a small offset above it (47μm in this work). We then run the shortest path algorithm again in the masked gradient image to obtain a new estimate of the retinal boundary location.

2.2.4. IS/OS layer boundary segmentation

In this step, we segment the B-scan for a third time to find the IS/OS layer. First we mask the location of the retinal boundary found in the previous step. In case our retinal boundary segmentation was already (by chance) the actual IS/OS layer, we only mask the locations where there is a strong brightness change below the current segmentation. To evaluate this condition we use a threshold for the maximum gradient value below the segmentation (in this work 0.3 after normalization). When applied, the mask starts a small distance below the segmentation (47μm in this work) up to the top of the image. The shortest path algorithm is run for the third time on the new masked gradient image. The layer boundary obtained corresponds to the IS/OS layer and it will be used in the next step to correct the first estimate of the retinal boundary.

2.2.5. Layer boundary refinement

To refine the first segmentation of the retinal surface, we impose a prior distance between that boundary and the IS/OS layer. We identify the locations where the retinal boundary and the IS/OS were closer than a given threshold (70.5μm in this work) and we label them as incorrect. In those locations, we correct the retinal boundary segmentation so that the distance to the IS/OS layer is equal to the median separation in the locations not labeled as incorrect. Figure 3(d) shows the final segmentation of the retinal boundary and IS/OS layers.

2.2.6. SVP image of the inner retina

Using the segmented retinal boundary, and the IS/OS layer boundary, we produce a projected image of the inner retina. Our aim is to obtain an image where the major blood vessels of the retina are highly contrasted with respect to the background. To do that, we compute three different depth values z₁, z₂ and z₃. The first value, z₁, is the depth of the retinal boundary plus 19μm. The second value, z₂, is equal to z₁ + d/1.5 + 47μm, where d is the distance between the retinal boundary and the IS/OS layer and. The third value, z₃ is the depth of the IS/OS layer plus 94μm. These values have been found using a training set of 7 pairs of X-fast/Y-fast volumes. The datasets used in the training phase were strictly separate from those used in the testing phase.

We compute the projected image of the inner retina by subtracting the average of the voxel values between z₂ and z₃ from the average of the voxel values between z₁ and z₂. Formally,

 

Fig. 4 Saccade detection based on SVP analysis. (a) & (d): SVP image of the X-fast and (rotated) Y-fast volumes of one subject (Section 2.2). Red arrows indicate saccades. (b) & (e): horizontal gradient of (a) & (d). (c) & (f): Green curve: negative of the correlation between lines (normalized). Blue curve: green curve minus a local average. Peaks in the blue curve determine the location of saccades.

Download Full Size | PDF

 

Fig. 5 Saccades detection and lateral eye motion correction for three example volumes. (a), (e), (c) & (d): original X-fast and Y-fast SVP of two different subjects; red arrows indicate locations of detected saccades. Each SVP image is divided into tiles at the locations of the saccades and the tiles true location and rotation are found using the orthogonal counterpart. (b), (f), (d) & (h): SVP images after correction; green, orange and red arrows indicate successful, partially successful, and unsuccessful corrections, respectively. Notice also how occlusions are revealed as dark gaps between the tiles.

Download Full Size | PDF

 

Fig. 6 Vessel likelihood maps. (a) & (c): original SVP images of the inner retina of two different subjects. (b) & (d): vessel likelihood maps obtained using the Gabor filtering process of [10]. These are used during lateral registration to weigh the regions of the image according to the presence of vasculature.

Download Full Size | PDF

2.3. Lateral registration

In this section, we describe how we use the SVP image obtained in the previous section to perform lateral registration of the volumes. To this end, we first perform saccade detection and correction, followed by dense local registration.

2.3.1. Saccades detection

OCT acquired images often suffer from saccadic motion artifacts. These can be seen as abrupt discontinuities in the en face projection of the volume (Fig. 1(b)). As a consequence, the saccadic affected regions of an SVP have a relatively lower correlation with the surrounding regions. We note that even without motion artifacts, the existence of large vessels can also produce low correlation with surrounding regions. Thus, in the X-fast volume, we detect the saccades using the correlation between rows of the horizontal gradient of the SVP (Fig. 4). The horizontal gradient is used to remove horizontal vascular structures that could produce a low correlation between rows that is not due to a saccadic motion. This process is detailed next.

Suppose I is the SVP image of B pixels height and W pixels width, where the B-scans are projected into rows of I. We filter each row of the SVP image I with a 1D Gaussian filter with σ = 6. Then, we compute the negative of the correlation between every two contiguous rows as

Let r be the vector of length B − 1 of correlations between every contiguous rows as computed by (2). We subtract from r a locally averaged version of itself, in order to remove low frequency changes in correlation,

We optimized the parameters σ, λ and r_min on a training set of 14 annotated SVP images containing multiple or no saccades. We use the exact same procedure for finding the saccades in the Y-fast SVP, only using columns instead of rows. Figure 4 shows example SVPs with the saccades detected by our algorithm.

2.3.2. Saccadic eye motion correction

After we detect the saccades, we need to correct for the patient eye motion that caused the discontinuities in the original image (Fig. 5). First, we find a global affine transform between the two orthogonal SVPs, as a pilot (rough) estimate of matching structures. This is done by fixing the Y-fast SVP image as reference and registering the X-fast image to it, using MATLAB’s imregister function. The affine transformed X-fast SVP and the Y-fast SVP images are then divided into tiles cut at the detected saccade locations. Note that the X-fast and Y-fast SVPs will yield horizontal and vertical tiles, respectively. Next, the tiles are matched one at a time, using the vertical tiles as the reference for the horizontal tiles and vice-versa. This process is described next.

We begin by creating two empty reference images, one for the horizontal tiles and one for the vertical tiles. The largest of all the tiles is used to initialize one of the reference images. Next, the tile with the largest overlap with the initialized reference image is chosen as the target. The selected target tile is registered against the reference using only their overlapping part. Next, the reference image of the target’s class is updated by pasting the registered target into it. This process is repeated until all tiles have been registered.

To further improve the registration accuracy, we consider not only the amount of overlap between tiles but also the presence of vascular structures in that area. To measure such presence of vascular structures, we compute a vessel likelihood map from the SVP images. This is done by filtering the SVP image with a bank of multiscale Gabor filters [10]. This filtering process produces a highly contrasted image where thin, elongated structures are highlighted, and the noisy background is filtered out (Fig. 6). By normalizing this image to have values between 0 and 1 we obtain a vessel likelihood map, in which the value of each pixel can be interpreted as the estimated probability of it belonging to a vessel. We use this likelihood map to weigh the area of a region according to the presence of vessels inside it. Consider a rectangular region R of the SVP image, the weighted area A_R of the region is

We use MATLAB’s function imregtform to find affine transforms between images. Since this operation involves a non-convex optimization [36], we try multiple parameters for the initialization (Transformation: affine and similarity. Metric: MSE and mutual information. Pyramid level: 2 and 3. Initial optimization radius: 10^−s, with s = 2, 4, 6) and keep the transformation that yields the highest correlation between the reference and registered parts. We also discard transformations that are too strong by measuring the scale, rotation and skew of the affine transform. Skew is measured as the angle of the lower corner of a transformed rectangle. We set the acceptable limits of these measures to |scale − 1| < 0.1, |rotation| < 3° and |skew − 90°| < 5°. A transformation that is off those limits is considered invalid. We did not set a limit for the translation range. If no valid transformation is found for a tile, that tile is rejected. In practice, this rarely occurs unless the tile is very narrow and there is poor structural information inside it. The gap left by the rejected tile can be filled with information from the orthogonal volume, as will be demonstrated in the experimental section.

The outcome of the saccadic eye motion correction are two SVP images that are the result of pasting affine transformed saccade-free image tiles. While these two SVP images match globally, there are some local minor displacements, which we will correct for as detailed in the next section. The global eye motion displacements have now been corrected, revealing occlusions or repeated imaging of the same location during acquisition (Fig. 5).

2.3.3. Orientation-aware dense local registration

To obtain an accurate and dense lateral registration we use a custom, “orientation-aware” optical flow algorithm that allows independent local displacements to occur inside each B-scan but not between different B-scans. Suppose I(x, y) is the ideal motion artifact free SVP, and I_X (x, y) and I_Y (x, y) are the SVPs of the X-fast and Y-fast volumes, respectively, and assume they are exactly orthogonal. Then,

This is a particular case of the optical flow problem, which is typically solved by a first order decomposition [37,38]. In our case, we first find u(y) in I_X, using I_Y as the reference and then find v(x) to correct I_Y. Our MATLAB based implementation of this particular case of optical flow runs the optical flow algorithm in [38] on the graylevel values of the SVP, and then keeps the components of the flow vectors that are in the direction of the B-scans (the exact orientation is obtained from the affine transforms of Section 2.3.2). Thus, any exchange of information between different B-scans is prevented.

Figure 7 shows the result of the local dense registration for the SVP images of one subject. The overlay between the X-fast and Y-fast SVPs is shown for a cropped region with high detail. In the first step of the “orientation-aware” optical flow algorithm, we correct for the local “horizontal” displacements of the horizontal B-scans. In the second step, we correct for the local “vertical” displacements of the vertical B-scans for a near perfect registration. We note that the terms horizontal and vertical are approximations as B-scans are slightly rotated. Figures 7(e)&(f) show the displacement vectors.

 

Fig. 7 Orientation-aware dense registration. (a) Original overlay between orthogonal SVPs in a cropped region of (d). Green colored vessels correspond to the X-fast SVP and pink to Y-fast. Blue arrows point to regions where displacements need to be corrected. (b) Overlay after horizontal dense motion correction. Blue arrows indicate regions where vertical corrections are required. (c) Overlay after vertical dense motion correction. The two images now match. (d) Resulting overlay of entire orthogonal SVPs. The red box indicates the location of the cropped region in (a)–(c). (e) Flow vectors from (a) to (b) (horizontal motion on X-fast SVP). Note that the orientation of the arrows are those of the original B-scans, so no information is exchanged between B-scans. (f) Flow vectors from (b) to (c) (vertical motion on Y-fast SVP). The resolution of the flow grids has been reduced for improved visualization.

Download Full Size | PDF

2.4. Axial registration

After obtaining the lateral registration, the remaining step to achieve full 3D registration is to correct residual tilts and offsets in the axial direction. To achieve this, we use the segmented NFL-vitreous boundary obtained as described in Section 2.2. Considering this segmentation as a depth measure, we correct the depth of each B-scan in the Z-direction, using the fast-scan directions as reference for the slow-scan directions.

We model the tilt and offset of each B-scan as a linear distortion on the Z-direction. More formally, let ĥ_X (x, y) be the observed depth of the retinal boundary in the X-fast volume at the (x, y) lateral position. Then,

Using equation (9) at each pixel location we obtain a linear system of equations. The solution to the system yields the values of m_X (y), n_X (y), m_Y (x) and n_Y (x), the tilts and offset coefficients for each row and column.

This model is however insufficient in practice because each image tile between saccades has a small but non-negligible independent rotation (see Section 2.3.2). Considering the rotations of the tiles, the equations are now given for each tile:

To make the axial registration process robust to outliers (due for example to errors in the segmentation), after a first registration, we identify outliers as the locations in which the two registered surfaces are separated by more than 23.5μm. To obtain a final estimation of the distortion coefficients, we remove the outliers and re-run the linear solver.

Figure 8 shows an example result of the depth registration. The original segmented retinal surfaces have strong motion artifacts in the slow scan directions. Figures 8(b) and (e) show the tilts and offsets found for each B-scan. By applying them on the original retinal surfaces we obtain two matching surfaces that are smooth on both directions.

 

Fig. 8 Axial registration. (a) Original segmented surface in the X-fast volumes. (b) Linear distortion found by solving (13). The distortions consist of linear displacements of each B-scan (rescaled with respect to (a) for improved visualization). (c) Resulting retinal surface after subtracting (b) from (a). (d) Original segmented surface in the Y-fast volumes. (e) Linear distortions found for the Y-fast volume (rescaled for improved visualization). (f) Resulting retinal surface after subtracting (e) from (d). As a result, the motion in the slow scan direction has been corrected on both volumes, producing matching smooth surfaces.

Download Full Size | PDF

2.5. System description

2.5.1. Hardware

A custom wide FOV SSOCT system was developed based on a commercially available swept source laser (λ₀ = 1050nm, Δλ = 100nm, Axsun Technologies, Inc.) with an A-scan rate of 100 kHz [39]. A custom retinal sample arm was created using a pair of orthogonal scanning mirrors and two clinically available wide-field, indirect ophthalmic lenses (Volk V20LC and V40LC; Mentor, OH). The retinal FOV of this OCT system relative to the posterior nodal point in a schematic model eye [40] was 80°. Scan time for both orthogonal volumes was just over 8.8 seconds per subject.

2.5.2. Imaging protocol

All human subject research was approved by the Duke University Medical Center Institutional Review Board. Prior to imaging, informed consent was received for all subjects after explanation of the possible consequences and nature of the study. All portions of this research followed the tenets of the Declaration of Helsinki. For each subject, we captured two orthogonally oriented OCT volumes: the first with the fast scan axis in the X-direction (X-fast) and the second with the fast scan axis in the Y- direction (Y-fast). Each volume consisted of 720 B-scans, each with 704 A-scans and 688 pixels per A-scan. The first and last 53 pixels from each A-scan were ignored resulting in an imaging depth of 2.7mm. Three seconds were allowed between orthogonal volumes to allow for the subject to blink.

2.6. Testing protocol

We found the best parameters for our algorithm, written in MATLAB, using a training dataset consisting of seven pairs of X-fast Y-fast volumes. The main parameters of the whole method are: filter window size for B-scan denoising, pixel value thresholds and distance thresholds for retinal segmentation (Section 2.2); filters size and cross-correlation threshold for saccades detection (Section 2.3.1); vessel map weight, initialization parameters and limits for the affine transform for lateral registration (Section 2.3); outlier threshold for the axial registration. The registration algorithm was then tested on a dataset of 24 pairs of X-fast/Y-fast volumes from 12 normal subjects. The entire process took on average 14 ± 2 minutes for a pair of X-fast/Y-fast volumes of 582×704×720 pixels on a Intel Core i7 CPU @3.50GHz with 64GB of RAM. No down-sampling was performed on the original volumes.

3. Evaluation

3.1. Qualitative results

Figure 9 shows the result of the lateral registration stage on the SVP images of the inner retina of three subjects. The column on the right shows the overlay of the images from the X-fast and Y-fast volumes before and after registration. Note the presence of multiple, strong saccades in the original images. After registration, the saccades were detected and the rigid motion of the eye between them has been corrected for. This process reveals regions of the retina that were not acquired in each volume (shown in black). The overlay of the registered images demonstrates the accuracy of the obtained registration.

 

Fig. 9 Results of the lateral registration for three test subjects, A, B and C. The SVP images were computed as described in Section 2.2.6. For each subject, the top and bottom rows show the X-fast and Y-fast SVP images, and their overlay before and after registration, respectively. In the overlay images, the X-fast image is shown in green and the Y-fast in pink. Note the presence of strong saccades in the original images, and the originally poor overlap of the vessels between them. After registration, the saccades have been detected and the rigid motion between them has been corrected for, producing an accurate overlap of the vasculature. Locations where no data was acquired are shown in black. The yellow and cyan lines indicate the locations of the cross-sectional images shown in Fig. 10.

Download Full Size | PDF

For the same subjects, cross-sectional images showing the result of the axial registration stage are presented in Fig. 10. The locations of the cross-sectional cuts are indicated in Fig. 9 with yellow and cyan lines for X and Y-direction, respectively. The black vertical bars represent the locations of missing data revealed after lateral registration. The X-direction in the X-fast volumes is the direction of the B-scans. Cross-sectional images of the Y-fast volume in X-direction correspond to the slow-scan direction and present strong axial motion artifacts, which make the true shape of the retina unrecognizable. After registration, the smooth shape of the retina is recovered along the slow-scan direction, whilst only requiring an independent shift and tilt of each B-scan in the fast-scan direction.

 

Fig. 10 Result of the axial registration step for the three subjects of Figure 9. For each subject, we show cross-sectional images along the X-direction (two columns on the left) and Y directions (two columns on the right). The locations of the X and Y slices are shown with yellow and cyan lines in Fig. 9, respectively. Before axial registration, the X-fast volume presents strong motion artifacts along the Y direction, and the Y-fast volume is corrupted along the X-direction. After registration, the smooth shape of the retina is recovered for the slow-scan directions, whereas for the fast-scan directions, only a tilt and an offset of the image is performed. The black regions correspond to the missing data gaps that have been revealed after lateral registration.

Download Full Size | PDF

Figure 11 shows results for two additional test subjects, D, and E, which we consider to be the most challenging volumes in our test dataset. The result is shown in the third column from the left, as the overlay of the two corrected datasets. In the rightmost column, we show how the gaps of missing data in one volume can be completed with the pixel values of the other volume. Both volumes contain multiple saccades, making the inter-saccades regions narrow and difficult to register due to lack of structure. Our algorithm is able to register and correct for lateral motion most of the SVP image of Subject D, except for the rightmost part of the image. For Subject E however, the extreme number of saccades makes the stripes to narrow and the algorithm fails to register most of them.

 

Fig. 11 Examples of shortcomings of our algorithm. Because of the extreme number of saccades, the regions between saccades are too narrow and their overlap contains too little information to achieve a successful registration. From the left, the first two columns show the original SVP images for the X-fast and Y-fast volumes of two test subjects. The third column shows the overlay of the registered images. Green corresponds to the X-fast image and pink to the Y-fast. The rightmost column shows the result of filling the gaps in the Y-fast image with the corresponding pixel values from the X-fast image. Small regions where no data has been acquired in any of the volumes appear in black.

Download Full Size | PDF

Figure 12 shows cross-sectional images for two additional subjects, F and G. The reference image on the leftmost column corresponds to a B-scan. The target image is the corresponding cross-sectional image along the slow scan direction in the orthogonally scanned volume. Note the presence of strong motion artifacts. The third column shows the overlay of the reference, and the target corrected w.r.t. the reference B-scan. The slow-scan direction has been corrected by applying a shift and a tilt to each of the B-scans. Finally, the rightmost column shows the result of filling the gaps of missing data with pixel values from the reference.

 

Fig. 12 Additional results of the axial registration step. The column on the left shows B-scan images for subjects G and F. The second column shows a cross-sectional image of the orthogonal volume in the same location, showing strong axial motion artifacts. The third column shows the overlay image of the registered images. Green corresponds to X-fast and pink to Y-fast. The rightmost column shows the result of registering the target with respect to the reference, plus filling the missing data by copying the pixel values from the reference. Note that the reference image (B-scan) was only rigidly transformed, but the shape was accurately recovered in the motion-corrupted image.

Download Full Size | PDF

The global result of the motion correction algorithm can be better appreciated in a 3D volumetric rendering of the retina, as is shown in Fig. 13. The image on the left shows the original X-fast volume of one of the subjects in the test set. Strong axial motion artifacts can be observed. The image on the right shows the result of registering the X-fast and Y-fast volumes and merging them. After motion correction, the smooth shape of the retina is recovered.

 

Fig. 13 3D rendering of the OCT volume of Subject B (Figs. 9 and 10) before and after motion correction. Left: original X-fast volume, presenting strong motion artifacts. Right merging of the registered and motion-corrected X-fast and Y-fast volumes.

Download Full Size | PDF

3.2. Quantitative results

To assess the performance of the registration algorithm we computed several quantitative performance metrics typically used in the OCT registration literature [22–24], at three stages of the algorithm. The three stages of the registration process considered are: original input volumes, after lateral registration, and after axial registration. We computed the following metrics, all of which show an important improvement after registration: mutual information (MI) between the two volumes; root-mean-squared error (RMSE) between the two volumes; median distance in μm between the retinal boundaries obtained from segmentation; mutual information between the SVPs of the inner retinal layers; and vessel map difference.

The vessel map difference metric [24] consists in measuring the mean absolute difference between the vessel likelihood maps of Section 2.3.2, Fig. 6. (Note that [24] computes the vessel likelihood maps differently.) A lower value of this difference means a better overlap of the vasculature between the two images, thus providing a rough evaluation of the registration performance at a fine scale, even if the vessel segmentation is not exact. Tables 1 and 2 show the average values obtained for these metrics for the 24 pairs of test volumes.

Table 1. Quantitative performance evaluation of the volumetric registration on 24 test volumes. For the RMSE calculation, pixel values are assumed to be between 0 and 255.

View Table | View all tables in this article

Table 2. Quantitative performance evaluation of the lateral registration on 24 test volumes.

View Table | View all tables in this article

3.2.1. Experiments on simulated motion

For further quantitative evaluation, we tried our method on volumes with known motion corruption. We took two X-fast/Y-fast pairs from our test set that do not contain any saccades (but do contain axial motion, Fig. 14). We applied the registration algorithm to each pair to obtain two registered and motion-free X-fast/Y-fast pairs, which we used as ground truth.

 

Fig. 14 Result of the algorithm on volumes with simulated motion corruption. (a) & (c): SVPs of two motion-free ground truth volumes. (a) corresponds to the X-fast volume of one subject and (c) to the Y-fast of a different subject. The red lines indicate the location of the cross-sectional images. (b) & (d): Motion-free cross-sectional images of the ground truth volumes. (e) & (g): SVPs of simulated motion corrupted volumes. The simulated saccades are indicated with red arrows. (f) & (h): Cross-sectional images of simulated volumes, showing motion corruption along the slow-scan direction. (i) to (l): Corrected versions of (e)–(h) using the proposed algorithm. Note how the simulated motion has been corrected, producing images that match the ground truth.

Download Full Size | PDF

Using the motion-free volumes, we created 60 simulated pairs of volumes corrupted with random motion. This was performed in a three steps process: 1) We created random saccades in the motion-free volumes (we randomly selected between 1 and 4 saccades per volume, with random locations); 2) We simulated lateral motion by applying random affine transforms to the obtained intra-saccade regions; 3) We simulated axial motion by applying a random offset and rotation to each B-scan.

We ran the registration algorithm on each of the 60 simulated volumes and compared the registration result to the ground-truth motion-free volumes. Figure 14 shows the ground truth motion-free volumes, one simulated motion corrupted volume for each subject, and the result of correcting the simulated X-fast/Y-fast pair.

We computed the following metrics: similarity between the SVP of the corrected volumes and the ground truth as Mutual Information (MI); median retinal surface difference between a corrected volume and the ground truth; precision and recall of saccades detection. The resulting values, averaged over the 30 pairs of simulated motion-corrupted volumes for each subject are presented in Table 3. The results confirm that our algorithm is effectively correcting the motion and producing motion-free images of the retina.

Table 3. Average registration performance of motion simulation/correction.

View Table | View all tables in this article

Note that, as in the real motion scenario, the correction is performed on two distorted volumes. Thus, the motion with respect to the ground truth can only be recovered up to a relative global motion. In the X–Y plane, the SVP of the corrected volume and the ground truth are related by a global affine transformation (Section 2.3.2). In the Z-direction, the corrected volumes and the ground truth are a linear depth correction away, given by the relative position of the two B-scans fixed as references (Section 2.4). Before measuring the registration performance, we found the global affine transformation between the ground truth and the corrected SVPs using MATLAB’s imregister function and the global depth correction by fitting a plane to the difference of the ground truth and the corrected depth.

The precision in saccades detection is lower because a few false positives occur, mostly due to artifacts of the segmentation around the optical nerve. However, the algorithm will typically find the correct location of the resulting regions. Note also that some misdetections can be explained because the random motion could be, by chance, not strong enough to produce a discontinuity.

4. Discussion

The results presented in this article demonstrate the applicability of our approach for 3D registration of wide field-of-view X-fast Y-fast OCT volumes with strong motion artifacts. The key ideas of our method are to separate the X–Y plane registration from the Z-axis registration, to include saccades detection, and to use as a guide the segmented retinal surface for depth registration. One important advantage of our method is that it guarantees minimal data transformation, involving only rigid transformations and small local displacements of each of the original B-scans.

Critical to obtaining good results is a successful segmentation of the retinal boundary. In this work, we used a simple but efficient approach. In terms of computation time, this is the most time consuming stage of the algorithm, as each B-scan has to be individually filtered and segmented. Indeed, in presence of significant retinal layer boundary deformations due to pathology there might be inaccuracies in performance of this algorithm. Note that the full registration process is robust to errors in the segmentation, as long as most parts of each B-scan are segmented correctly. In case of more severe deformations, thanks to the modularity of our approach, more robust to deformation (albeit significantly more time consuming) segmentation methods can be plugged into the algorithm’s pipeline (e.g, [41,42]). Accurate quantification of registration performance in the presence of various types of pathologies is part of our ongoing research.

Another critical part is the detection of saccades and correction of the SVP tiles as described in Section 2.3.2. Our saccade detection method is based on computing the cross-correlation between contiguous lines. The region of the optical nerve can sometimes affect this correlation causing false negatives. Also, the affine registration of the tiles requires enough overlap between them, so that there is sufficient information to compute the affine image registration. If the intervals between saccades are too small, or present no vasculature or other structural information, the registration could fail but still be within the limits of an qualifying transformation. Furthermore, the error of incorrectly registering one tile can be propagated to the other tiles that are yet to be registered. Such a case is exemplified in Subject E of Fig. 11, which presents an extremely high number of saccades and little structural information in large portions of the image.

Finally, the depth registration is obtained by solving a linear system using least squares. The L-2 norm is known to be non-robust to outliers. In this case outliers can occur in regions where the retinal segmentation has failed. The solution we proposed is to detect outliers, remove them and re-iterate the least squares computation. A possible line of future research is to study the performance of more robust linear methods.