Automated non-rigid registration and mosaicing for robust imaging of distinct retinal capillary beds using speckle variance optical coherence tomography

Variance processing methods in Fourier domain optical coherence tomography (FD-OCT) have enabled depth-resolved visualization of the capillary beds in the retina due to the development of imaging systems capable of acquiring A-scan data in the 100 kHz regime. However, acquisition of volumetric variance data sets still requires several seconds of acquisition time, even with high speed systems. Movement of the subject during this time span is sufficient to corrupt visualization of the vasculature. We demonstrate a method to eliminate motion artifacts in speckle variance FD-OCT images of the retinal vasculature by creating a composite image from multiple volumes of data acquired sequentially. Slight changes in the orientation of the subject’s eye relative to the optical system between acquired volumes may result in non-rigid warping of the image. Thus, we use a B-spline based free form deformation method to automatically register variance images from multiple volumes to obtain a motion-free composite image of the retinal vessels. We extend this technique to automatically mosaic individual vascular images into a widefield image of the retinal vasculature.


Introduction
The primary clinical application of optical coherence tomography (OCT) to date is in retinal imaging.While being used mainly as a diagnostic tool for examining the morphology of the retina, advances have been made in extending the use of OCT to image functional properties, such as blood flow.Doppler OCT has been used to measure blood flow in the vasculature [1][2][3][4][5][6], providing a means to detect the total amount of blood flowing into and out of the inner retina by using either circumpapillary [7] or angle-independent scanning methods [6].However, the information gained from retinal Doppler measurements has not yet been translated to clinical value.Angiography techniques such as fluorescein [8] and indocyanine green angiography [9] allow for detection of areas of vascular perfusion and have a demonstrated clinical utility.However, these methods require the injection of dyes into the patient bloodstream that could cause potential discomfort and possible allergic reactions.Noninvasive optical methods have been developed to image the retinal vasculature.Laser speckle imaging demonstrated the ability to image retinal vessels due to changes in the speckle pattern of moving blood cells [10].Commercially available systems such as the Heidelberg Retinal Flowmeter [11] or the Retinal Function Imager [12] enable imaging and quantification of retinal flow.Adaptive optics scanning laser ophthalmoscopy has been used to image retinal capillaries but only over a very limited field of view [13].However, these imaging modalities are inherently two-dimensional and lack the ability to visualize the location of the vessels in depth.The inner retina is perfused by several distinct vessel beds that lie predominantly in the ganglion cell layer and at the inner and outer borders of the inner nuclear layer [14].These layers contain distinct vasculature that is responsible for the perfusion of the inner retinal tissues to provide nutrients to meet their metabolic demands.The ability to resolve each of these layers in depth could improve understanding of the onset and development of retinal diseases that may affect each vascular bed [15].
Recent advances in FDOCT technology have produced methods including optical microangiography [16], joint time-spectral OCT [17], speckle and phase variance imaging [18,19], and power or variance Doppler techniques [4,20] that enable non-invasive capillary level detection of the retinal vasculature.Recently, a method quantifying the level of perfusion based on the vessel density has been demonstrated [21].The use of increasingly high-speed systems has allowed these methods to overcome the effects of capillary blurring due to subject motion during acquisition of densely sampled B-scans.However, acquiring volume data sets still requires several seconds of imaging time, even with high-speed systems.
For even normal subjects, stable fixation is difficult to achieve over the time course of data acquisition, and motion artifacts corrupt the visualization of portions of the data in a given volume.Eye motion during data acquisition can manifest as a slow shift in gaze (drift), high frequency involuntary motions (tremor), or rapid eye movements (saccades) [22].Subject eye motion results in discontinuities in OCT images, and saccades appear very prominently as streaks along the fast scan axis in angiography data sets due to the high amount of subject motion between successive B-scans, as shown in Fig. 1.Additionally, the effects of motion artifacts in retinal imaging produce not only simple discontinuities throughout the image, but also distortions and image warping in features of interest [23].Retinal OCT systems are typically designed to scan a collimated beam through a pivot location at the center of the pupil of the eye [24].Changes in the orientation of the subject's eye relative to the scanning optics may cause the sample beam to pass through the pupil at different positions or angles, resulting in slight differences in the appearance of structures, such as blood vessels, from one volume to another.Similar types of distortion are also responsible for differences between OCT B-scans in the apparent curvature of the retina, which in some cases may appear flat while in other cases will have either a concave or even a convex curve to the tissue structure [25].These distortions could render rigid transformations insufficient to register images of the microvasculature in the retina from different volumes.
Axial motion artifacts are commonly corrected by using cross-correlation methods and rigid translations of each B-scan [24,26] or through acquisition of orthogonal B-scans to aid alignment by providing a slow scan axis reference frame with no axial motion [27].However, these methods do not correct for lateral motion artifacts.To address this problem, early methods proposed capturing several fast sparse scans, detecting and removing the scans affected by lateral motion artifacts, and fusing the artifact-free scans [28].However, this approach required long image acquisition times because full data sets corrupted by lateral motion were discarded.Another correction method used orthogonal scan patterns to obtain data sets containing uncorrelated motion and subsequently reconstructed a motion-free image by combining information from both data sets [29,30] but required heavy computational power to register even two images.A scanning laser ophthalmoscope (SLO) may be used for retinal tracking [31,32], and other motion correction methods for OCT use the addition of an SLO system to either acquire a reference SLO image [33] or to rescan regions of the volume that suffered from motion [34].Such systems, however, have greatly added complexity due to the additional instrumentation needed.An additional consideration is that slight changes in a subject's gaze and positioning of the optics relative to the pupil can result in slight warping effects in each data set, making it difficult to register images with rigid transformations.Automated mosaicing has been carried out using rigid transforms of the intensity projection data with OCT [35] but did not correct for motion or image warping and has not been applied to images of the microvasculature.It has been demonstrated in other medical imaging modalities that a free-form deformation (FFD) model may account for non-rigid distortions when registering multiple data sets acquired at different points in time [36].Here, we present an automated post-processing technique to remove motion artifacts in speckle variance data sets of the retinal vasculature using localized deformation of the vessel structure to align and register multiple volumes.By using multiple data sets, motion artifacts may appear at different locations in different volumes, as shown in Fig. 1.A motion-free composite image can be constructed by stitching together uncorrupted image segments from each data set.Orthogonally scanned volumes may also be acquired to ensure decorrelated motion artifacts between the different volumes.Using automated segmentation of the retinal layers, the motion correction can be applied to each retinal layer to create images of each of the major vasculature beds in the inner retina that are free of motion artifacts.We also extend the method to automated widefield mosaicing of smaller volumes acquired over a larger retinal surface area.

System description
A high speed retinal OCT system shown in Fig. 2 was built using a commercial swept-source laser (Axsun, λ o = 1060 nm, Δλ = 100 nm, 100 kHz repetition rate) with 1.8 mW incident at the pupil of the eye in compliance with the ANSI maximum permissible exposure limit [37].A fiber Bragg grating (λ o = 989 nm, Δλ = 0.042 nm) was used to trigger the acquisition of every A-scan to ensure phase stability from one A-scan to the next and to allow for removal of fixed-frequency artifacts.A retinal scanning system consisting of imaging lenses and an X-Y galvanometer pair was used in the sample arm.A fixation target was provided with an LCD display to direct the subject's gaze during imaging sessions.A balanced detection scheme was used to reduce relative intensity noise.Human volunteers were consented and imaged in accordance with institutional review board approval and the declaration of Helsinki.oriented along Y.Each B-scan within a volume was axially registered using cross-correlation with the previous frame to account for subject motion along the axis of the imaging beam [38].For each set of 3 repeated B-scans, the magnitude data after Fourier transformation was averaged to improve the image SNR.The variance was computed across the 3 repeated Bscans at each location [18], resulting in 300 variance B-scans.

Automated image registration procedure
This section describes the method implemented for performing automated image registration for motion correction using multiple speckle variance OCT data sets.The main retinal feature used to guide the registration algorithm is the vascular network.Speckle variance data have greater vessel contrast compared to the reflectance data, which show greater contrast between the retinal layers.Thus, the reflectance OCT B-scans were used to segment the retina into distinct layers while the registration process operated on the variance images.In this section, we first describe image segmentation methods to create images of the distinct vasculature of each retinal layer.Next we describe the process used to detect motion and divide each projection image into motion-free strips.We then describe additional image filtering followed by our global-local registration method and final composite image creation.Figure 3 shows an overview of the steps used to create a motion-free composite image of the retinal vasculature.

Image segmentation
Image segmentation enables division of the retina into layered regions corresponding to the different vascular beds of interest.Variance data sets were collected from the right eye of each subject according to the imaging protocol discussed above, and equal numbers of data sets with scan priorities on the X-and Y-axis were acquired and processed.The averaged intensity B-scans were then automatically segmented to detect 8 boundaries corresponding to 7 distinct retinal layers [39].Three additional boundaries were created relative to the automatically segmented inner nuclear layer (INL), with the added boundaries positioned at the center and at

Motion detection and image sub-division
We first detected the location of the motion artifacts in each SVP by summing each 2D image along the fast scan axis to create a 1D projection that was normalized according to the mean of the projection.B-scans that were affected by motion resulted in a higher variance across nearly all the pixels corresponding to tissue structure in a variance B-scan.An intensity threshold for values greater than 2.5 times the mean was used to detect points in the 1D projection with a high variance, thus indicating a B-scan that suffered from motion artifacts.This threshold value was determined empirically from data acquired in 2 different subjects and was found to be sufficient for detecting motion artifacts.Each variance SVP was then divided at the location of the detected motion, yielding a sequence of strips that were free of motion artifacts.Image strips that were below a threshold size (< 25 B-scans in width) were discarded, as the amount of information contained within would be too small for accurate registration.Combining information from multiple volumes allowed for recovery of information over the entire retinal area scanned without the need for additional interpolation.Histogram equalization was applied to each strip to normalize intensity differences between variance B-scans.An example of this procedure is shown in Fig. 5. Thus, for all volumes acquired at the same retinal location, a total of N image strips of the GCL + IPL were created for registration.The SVPs from the IPL-INL and INL-OPL junctions were also divided at the same locations as the images of the GCL + IPL.

Gabor filtering
Gabor filtering enabled enhanced visualization and contrast of the retinal vascular network.A Gabor filter is a Gaussian filter modulated by a cosinusoid that acts as a directional low-pass filter and can be used to enhance connectivity in images of the vasculature [40].The standard form of the Gabor filter is given by ( ) In the above equations, m x and m y are the spatial coordinates of each point in the Gabor kernel of size m × m, θ is the angle of orientation of the filter, ν defines the spatial frequency of the filter, and σ x and σ y are the standard deviations of the Gaussian factor in the x and y dimensions.The filter is convolved with the image using varying values of θ and yields a maximum response at each pixel when the filter is aligned locally with the vessels in the image.Keeping the maximum response at each pixel allows for enhanced contrast of the vessels and reduction of background noise.Additionally, to attain the multiresolution property, the size of the Gabor kernel was increased k times, as smaller kernels give a better response to capillaries while larger kernels enhance the response to the larger vessels.Summing the responses of each Gabor kernel yielded images where both capillaries and larger arterioles and venules were preserved.
where ⊗ is the convolution operator.The parameters for the Gabor filter were empirically chosen to give the optimum response in the images shown.For our data, 10 equally spaced angle orientations were used, and the filter kernel was varied in size from 5 × 5 to 50 × 50 pixels with σ x = σ y due to the equal sampling in the x and y dimensions of the image.Gabor filtering was applied to each motion-free strip after motion detection and image sub-division.The Gabor images were cropped by 5 pixels on all sides to avoid edge effects of the filter.An example of Gabor filtering using these parameters of a representative vascular image strip is shown in Fig. 6.Gabor filters were applied to the images of the vasculature at each depth.

Global registration
In order to perform accurate registration of the vasculature, correct global placement of each strip relative to a starting reference image was required.All image strips were zero padded to ensure sufficient ability to register and place each strip in a composite image.In our data sets, each SVP was 300 × 300 pixels and were zero padded to 400 × 400 pixels.A binary mask, b n (x,y), that was equal in size to each image strip, g n (x,y) was created and similarly zero padded.Each b n was transformed using the same parameters as its corresponding g n .The largest image strip was chosen as the starting reference template, g r (x,y), with its binary mask b r (x,y).The largest g n (x,y) orthogonal to the reference was then shifted according to its initial position in its respective SVP.This assumed that the subject was able to re-fixate on the same location after any involuntary motion.To relax that assumption, an additional global x-y translation of g n (x,y) was determined by maximizing the correlation between g r and g n [38,41].A binary mask, described in Eq. ( 4), provided a simple method to constrain the correlation calculation to use only the portion of the image strip that overlapped the reference, which helped to ensure that only sections of the image that contained the same vasculature were used to determine the registration.The effects of noise on the registration were reduced through Gabor filtering.
[ ] [ ] G r (u,v) and G n (u,v) are the 2D fast Fourier transforms (fft) of g r (x,y) and the nth image strip, g n (x,y), respectively, with spatial frequencies u and v. Multiplying G r with the conjugate of G n and taking an inverse Fourier transform yielded a discrete pixel shift (Δx max , Δy max ) in the image that maximized the correlation.This resulted in a translated image.
n n g x y g x x y y → +Δ +Δ (7) An example of this global strip registration is shown in Fig. 7 in the middle column.

Local optimization
Due to geometrical distortions in images between different volumes, local deformation of the image of the vasculature was required to properly align vessels of different sizes.Following the initial global image placement and orthogonal registration, g n (x + Δx,y + Δy) was locally registered to g r (x,y).Free-form deformation (FFD) of the images was carried out using a spline-based registration algorithm based on Rueckert, et al. [36], which uses a cubic B-spline basis to transform an image using a set of control points distributed over the image.B-splines are useful for producing smooth and continuous non-rigid deformations over an image surface [42].Given a set of control points, φ i,j , with p x × p y equally spaced points, the B-spline transformation for each pixel in an image can be described as [36] 3 3 , 0 0 ( , ) ( ) ( ) , ' ( ', ') ( ( , )).
The optimal T(x,y) was found by iteratively adjusting φ using a gradient descent method [36] that minimized the sum of the squared differences between g r (x,y) and g ' n (x + Δx,y + Δy).To locally register g n (x + Δx,y + Δy) to g r (x,y), the area of their overlap after global registration was detected using their binary masks in a similar manner as described for global correlation, and registration was performed using only this region.Using only portions of the images that overlapped helped to ensure a reasonable FFD transformation by removing areas that should not be expected to merge with the reference image.A multi-scale approach was used to refine an initial spline registration control grid of 4 × 4 points to a final grid of 20 × 20 points.Regions of g n (x + Δx,y + Δy) that did not overlap with the reference were constrained against deformation.An example of the result of FFD registration is seen in Fig. 7 in the right most column.

Composite image generation
A new composite reference image was created after the global and local registration of each individual strip.The images g r and g' n were summed together, and pixels with non-zero contributions from both images were averaged and equally weighted.This resulted in a new reference image that was used to register g n + 1 , and the steps from Sections 3.4 and 3.5 were repeated.Once all g n (x,y) had been registered, the final reference image yielded the motion corrected SVP of the retinal vasculature.The transform parameters at each step obtained from using the GCL + IPL images were extracted and applied to the image strips from the retinal layers of the IPL-INL and INL-OPL junctions.Images of the 3 distinct vascular beds were also combined into a single color encoded depth image by writing the intensity values for each layer to a different RGB channel.Pixels in the resulting image with values from multiple color channels revealed regions with overlapping vasculature from different retinal layers.

Automated widefield mosaicing
The previous processing steps allow for registration of multiple volumes acquired over the same area of the retina in order to correct for motion artifacts.An extension of these methods is applicable to automatically register images acquired over different areas of the retina to form a widefield mosaic.The steps below outline the process for automated widefield mosaicing.

Data acquisition and initial motion correction
Data was collected at N different locations of the retina around the foveal avascular zone.Locations for patient fixation were chosen to ensure approximately 50% lateral overlap between data sets.At each location, the same acquisition protocol as in Section 3.1 was used, with two X and two Y oriented volumes acquired at each region.After acquisition, all data sets were processed as described in Section 3 with motion-free composites of the GCL + IPL, IPL-INL, and INL-OPL layers created at each distinct retinal location.Images from the ganglion cell plexus were used to compute the image transforms for widefield mosaic registration.Here, I n (x,y) will denote the nth motion-corrected GCL + IPL image.The motion corrected image corresponding to the first region of the retina that was imaged served as the initial reference frame, I r (x,y), for building up a mosaic.The remaining motion corrected images were registered in the order in which the retinal areas were imaged.This ensured that each image would overlap with preceding data sets in the mosaic, although the amount and location of the overlap were unknown.The transform parameters could then be applied to SVPs of the other retinal layers to provide layer specific widefield mosaics of each vascular bed.

Global registration
I n (x,y) was registered to I r (x,y) using image correlation in a manner similar to that described in Section 3.4.The images were registered in the order of acquisition and so were guaranteed to overlap with previous images.However, because the exact location of I n (x,y) relative to I r (x,y) was unknown, the following process was used to determine a global registration.Binary masks 150 × 150 pixels in size were used to perform a windowed correlation of segments from I n (x,y) with segments from I r (x,y).Each mask was iteratively and independently shifted over each image, and the correlation was calculated over the windowed regions of each image in order to find the optimum image position for I n (x,y).The correlation yielding the maximum peak was used to compute the appropriate x-y translation for I n (x,y).

Local registration
After global image placement, the free form deformation method used in Section 3.5 was applied to register each image to the current reference mosaic.A similar method of constraining the FFD to operate on regions of overlap between the two images was used.

Mosaic generation
After each iteration, I r (x,y) was updated with the current image and used as the new template for I n + 1 (x,y).I r (x,y) was updated by averaging the contributions of all I n (x,y) to each pixel in the mosaic in a similar manner to the image composite generated in Section 3.6.The steps in Sections 4.2 and 4.3 were repeated for all N images.This resulted in the final motioncorrected widefield mosaic of the ganglion cell plexus.The same image transform parameters were then applied to the SVPs generated from other retinal layers to produce widefield images of the IPL-INL and INL-OPL junctions.A color encoded depth widefield image was also created in a manner similar to that described in Section 3.6 using the widefield mosaics of the 3 retinal layers.

Automated motion artifact correction
Multiple data sets acquired from a healthy 27-year-old human volunteer are shown in Fig. 8.This data set was used to determine the parameters for motion detection, Gabor filtering, and registration in all subsequent data sets.Two X-fast and two Y-fast volumes were acquired, and the SVPs from each layer show evidence of motion artifacts in all volumes.Running the registration algorithm from Section 3 resulted in motion-free composite images for each layer, shown on the right of Fig. 8.An enlarged view of the registered images and SVPs from a single volume is shown in Fig. 9 along with a combined color encoded depth image of the vasculature.After processing, the streak artifacts apparent in the raw image volumes were no longer present, and vessel visualization throughout the image was enhanced in the registered, motion corrected images.Differences in the vasculature of the different layers are apparent, with capillary beds prominent in the layers surrounding the INL and larger vessels present in the GCL + IPL.The distinct vasculature of both the superficial and the deep capillary plexuses was also observed.To demonstrate the robustness of the algorithm in operating on data from different subjects, data from a second volunteer imaged at an area nasal to the fovea is illustrated in Fig. 10.A total of three X-fast and three Y-fast data sets were acquired, and the

Automated widefield mosaic generation
To demonstrate the ability to perform automated widefield mosaicing, data sets from a volunteer were acquired over 12 different regions of the retina.At each location, two X-fast and two Y-fast data sets were recorded.After segmentation, SVPs from each retinal layer were created, and the image registration algorithm of Section 3 was performed.The widefield mosaic algorithm from Section 4 was executed on the motion corrected SVPs from the GCL + IPL.Each of the 12 motion corrected images is shown in Fig. 12 along with the result of the automated mosaic algorithm.The method correctly placed each of the individual images and showed good registration among the overlapping vessels.layer were then applied to the motion corrected images from the INL junctions to generate mosaics of the superficial and deep capillary plexuses in Fig. 13 and Fig. 14, respectively.The mosaics cover approximately a 6 × 4 mm field of view, showing the major vessels and capillaries over a large area as well as the differences between each of the 3 main vasculature beds of the inner retina.A color encoded depth image combining information from each vascular bed mosaic is shown in Fig. 15.

Discussion
The results described in this paper have demonstrated the ability to segment retinal data to perform layer specific angiography using speckle variance OCT.The proposed motion correction algorithm has been demonstrated to give good results for creating motion-free composite images from multiple volumes of data and can also be used to generate widefield visualization of the vasculature.
One of the key processing steps in the algorithm is the application of Gabor filters to enhance the appearance of the vessels in the SVPs.Applying the Gabor filter causes additional smoothing of the vessel edges that may cause vessels to appear larger than their actual size.However, the aberrations in the eye already provide a limiting factor on the lateral resolution.While variance processing methods in OCT are able to visualize the location of capillary structures, this does not imply that their size is accurately represented.Retinal capillaries are approximately 5 µm in diameter, yet due to the resolution limitations of the optics of the eye, the lateral resolution of standard OCT retinal systems is 10-20 µm.Vessels with a diameter smaller than the lateral resolution will appear to be blurred to the beam spot size.Thus, the use of Gabor filtering is not likely the limiting factor in image resolution.Measuring the exact capillary dimensions is not possible without the use of adaptive optics [43].However, the presence of the capillaries may still be detecting using variance based processing if the sampling density is sufficiently high for each B-scan.By oversampling with a high number of A-scans, image pixels containing even small blood vessels will show a higher variance due to the presence of moving scatterers within each image voxel.For detecting motion in the variance SVPs, a threshold of 2.5 times the mean variance over all B-scans was used to detect B-scans that were affected by patient motion.This value was chosen empirically based on data acquired from the 2 subjects presented here.Though the threshold value worked well for detecting motion in the data from this report, an optimal threshold value may be determined if a larger set of subjects were to be analyzed.As signal from the vasculature also produces a high variance, and it may be possible for B-scans containing large vessels oriented parallel to the fast scan axis to be falsely identified as being affected by motion if the vessel spanned much of the length of a given B-scan.However, vessels typically possess some level of curvature that would cause a B-scan to partially bisect only a segment of the vessel.Thus, it would be expected that false positives due to large vasculature would be minimally present in most cases.
Imaging subjects with poor fixation ability may pose a problem.The motion correction method described here requires that subjects be able to re-fixate upon a target after any motion event.This allows multiple data sets to be acquired over the same retinal region, however, a subject who is unable to consistently focus or rapidly refocus would prevent the necessary data from being acquired.
The A-scan spacing used in our acquisition protocol was 8.3 µm in both x and y dimensions.While this sampling density is still larger than the diameter of capillary structures, the presence of microvessels within a given image voxel was still sufficient to cause detectable differences in the intensity variance compared to voxels lacking vasculature.Scanning a smaller area would allow for denser sampling to enhance detection of capillary vessels.However, the field of view would be reduced.Visualization of a wide area over the retinal surface is important because pathologies may affect specific retinal regions while leaving other areas appearing normal.Thus, there is a trade-off in imaging a larger field of view for a single volume while still trying to maintain a high sampling density for capillary visualization.Use of the widefield mosaic algorithm could offer a method to increase the field of view by acquiring densely sampled volumes to enhance the detection of microvessels.
In certain areas of the motion corrected images, the connectivity of the smaller vessels may be difficult to observe.This artifact may be due to insufficient sampling of the vessel network.Increasing the A-scan sampling rate as well as the lateral resolution of the system would improve visualization of these vessels; however, this would come at a cost to either the field of view or an increase in imaging time.
The larger vessels in the retina are expected to lie in the GCL + IPL with smaller vessels feeding into the superficial and deep capillary plexuses.The superficial capillary plexus region was defined to be 16 µm above the detected anterior boundary of the INL.However, in some cases, the larger vessels appeared to enter into this region, as can be seen from the yellow coloring of some of the major vessels in Fig. 15.This may indicate that the superficial capillary plexus either varies in width over the retinal surface or is thinner than the 16 µm region chosen to capture the vessels about the INL.
The image transform parameters for registration of the GCL + IPL images were directly applied to the registration of the IPL-INL and INL-OPL images.It was assumed that the acquisition rate of an A-scan was instantaneous relative to the time of a B-scan.However, slight registration artifacts may still manifest in the images of the vasculature IPL-INL and INL-OPL, which may indicate that an additional correction factor is necessary.It is possible that a direct application of the image transform parameters from the GCL + IPL image is not the most optimal solution for all layers.
Use of orthogonal scan patterns in this work was motivated by the observation that multiple volumes scanned along the same axis from the same retinal area do not always guarantee that motion artifacts will be randomly distributed from one volume to the next.It is possible that regularly occurring events may cause the subject to move at the same point in each volume.For example, when fixating during imaging, subjects are still able to see the scanning sample beam before it is incident on the fovea using peripheral vision.As the scanning beam passes near the periphery of the fovea, many subjects tend to temporarily shift their fixation to the location of the moving beam, thus causing a motion artifact that is reproduced in successive volumes.These consistent motion artifacts would prevent correction of motion at that location and motivated acquisition of both X-and Y-fast scans.
The computation time of the registration algorithm in Section 3 for registering data from 4 volumes was 3.5 minutes using a Pentium i5 2.8 GHz single core processor in Matlab®.Generation of the widefield mosaic from 12 motion-corrected images required an additional 14 minutes of processing time.It may be possible to optimize the algorithm and implement parallel processing to increase the speed of the registration.However, this does not include the time to process and segment the data sets after acquisition, which increases the amount of time before useful data can be presented.For clinical use, it would be desirable to reduce the processing time for the entire process.Use of graphics processing units (GPU) have shown much promise for increasing the rate of data processing [44,45] and can assist in making this method useful for immediate analysis of data in the clinic.
Improved image acquisition protocols may be used to optimize data collection for widefield mosaicing.It may be possible to mosaic image with less overlap than was used in this study, if the more accurate (yet more computationally expensive) multi-frame joint image registration algorithms are utilized [46][47][48].By increasing the distance between retinal locations imaged, total imaging time per subject may be reduced or a wider field may be imaged.
Imaging about the peripapillary region would pose additional difficulty with our method.Segmenting the retina near the optic nerve is difficult, as the distinct retinal layers narrow and appear to merge together as they approach the nerve head.This would prevent accurate separate of the different vascular layers.Additionally, the presence of large blood vessels would cause shadowing of deeper retinal layers regions beneath the large vessel, thus preventing detection of the deeper vasculature.Thus, it would be difficult to create segmented vascular maps of the distinct retinal layers.However, the presented method of global and local registration should still be applicable if properly segmented and clearly detected vessel maps could be generated.
Although speckle variance OCT was used in this study, the motion correction algorithm operates on the projection images of the vasculature.Thus, other OCT techniques used to create angiograms, such as phase variance or power Doppler methods, may also be used with this method of motion removal.

Conclusion
We have demonstrated an automated post-processing method to correct for motion artifacts and generate widefield mosaics in speckle variance OCT images with layer specific visualization of the retinal vasculature.We created motion-free composite images from multiple orthogonally scanned volumes using a non-rigid free form deformation technique to correct for image warping between volumes.We have shown the ability to visualize the 3 distinct vascular layers in the inner retina as well as the ability to automatically generate widefield mosaic images of the retinal vasculature from multiple volumes acquired over different regions of the retina.This method may help improve the clinical viability of OCT by enhancing the ability to visualize and diagnose diseases of the retinal vasculature.

Fig. 1 .
Fig. 1.Example of retinal motion artifacts in speckle variance OCT data sets.Two volumes acquired from the same subject at different times show horizontal streaks along the fast scan axis due to subject motion.

Fig. 2 .
Fig. 2. System schematic of the swept-source optical coherence tomography retinal scanner.Data sets consisted of a 2.5 × 2.5 mm volume scan with 300 A-scans/B-scan, 3 repeated B-scans, and 300 unique B-scan positions, thus yielding a total of 900 B-scans in each volume.This yielded a lateral sample spacing of 8.3 μm in both scan dimensions.For scanning stability, 75 A-scans were added to each B-scan to allow for flyback and repositioning of the fast galvo, yielding a B-scan time of 3.75 ms and a duty cycle of 80% with a total volume acquisition time of 3.4 s.Multiple data sets were acquired from each subject, and half the volumes were oriented with the fast scan axis along X and the other half
16 µm anterior and posterior to the INL.The size of these boundary regions was chosen to ensure that the capillary plexuses were fully captured in the SVPs.All boundaries were #187382 -$15.00USD Received 19 Mar 2013; revised 25 Apr 2013; accepted 27 Apr 2013; published 7 May 2013 (C) 2013 OSA 1 June 2013 | Vol. 4, No. 6 | DOI:10.1364/BOE.4.000803 | BIOMEDICAL OPTICS EXPRESS 808 constrained so as not to cross neighboring boundaries near the fovea where the retinal layers merge together.Summed volume projections (SVPs) were automatically created for each of the retinal layers from the variance B-scans using the segmentation data obtained from the averaged intensity B-scans.An example of the segmentation and layer specific SVP creation from a sample volume data set is shown in Fig. 4. The SVP corresponding to the ganglion cell plexus was created by summing the variance B-scans over the ganglion cell + inner plexiform layer (GCL + IPL in Fig. 4B).At the INL junctions, 2 SVPs were created spanning 16 µm above the INL to the mid-point of the INL (superficial capillary plexus) and from the mid-point of the INL to 16 µm below the INL (deep capillary plexus).Image registration was performed on the SVPs of the ganglion cell plexus to take advantage of the larger vessels located there, which were helpful in image alignment.The transform parameters were then applied directly to the SVPs of the superficial (IPL-INL junction) and deep (INL-OPL junction) capillary plexuses.

Fig. 4 .
Fig. 4. Automatic segmentation of retinal layers and corresponding variance SVPs.(A) Automated segmentation of 7 retinal layers according to [39].NFL: Nerve Fiber Layer, GCL + IPL: Ganglion Cell Layer + Inner Plexiform Layer, INL: Inner Nuclear Layer, OPL: Outer Plexiform Layer, ONL: Outer Nuclear Layer, OS: Outer Segments, RPE: Retinal Pigment Epithelium.(B) Expanded view of the inner retina from (A).Dashed lines indicate added boundaries used to create SVPs of the vascular beds in the GCL + IPL, IPL-INL and INL-OPL junction.Boundaries were created at 16 µm above (green dash), in the middle (red dash), and 16 µm below (yellow dash) the INL.(C) SVPs of each segmented retinal layer.The GCL + IPL (green), IPL-INL (red), and INL-OPL (yellow) contain the distinct vessel layers of interest and correspond to the shaded regions in (B).The choroid image was summed over a depth of 30 µm below the choroid boundary.Each SVP was 2.5 × 2.5 mm in size.

Fig. 5 .
Fig. 5. Example of image sub-division at locations of motion artifacts.Image strips on right have undergone histogram equalization to remove slight intensity differences between B-scans.

Fig. 6 .
Fig. 6.Example of multiresolution Gabor filtering on speckle variance image of retinal vessels.

Fig. 7 .
Fig. 7. Examples of image registration process.Top row shows images after global and local registration.Bottom row shows magnification of the region within the red box on each of the top images.Yellow arrows point out several vessels that show improved registration after free form deformation.
performs a flooring operation.The B-spline functions are defined as for interpolating each pixel in g n (x,y) to a new location, g' n (x',y'), based on the shift of each local control point in φ.A bilinear interpolation method was then used to #187382 -$15.00USD Received 19 Mar 2013; revised 25 Apr 2013; accepted 27 Apr 2013; published 7 May 2013 (C) 2013 OSA 1 June 2013 | Vol. 4, No. 6 | DOI:10.1364/BOE.4.000803 | BIOMEDICAL OPTICS EXPRESS 812 compute the new intensity values at each pixel based on the deformation of g n according to T(x,y).

Fig. 8 .
Fig. 8. Registration result from subject 1.Two X-fast and two Y-fast data sets were acquired.The original SVPs for each of the 3 main vessel layers are shown in the 4 left columns.The right most column shows results of image registration for each of the 3 layers.Each image covers a 2.5 × 2.5 mm scan area.

Fig. 9 .
Fig. 9. Enlarged view of images in Fig. 8 comparing the registered images to SVPs from a single volume.Motion artifacts are removed and visualization of the vasculature is enhanced in the registered images.A color encoded depth image is shown on the right, combining information from the registered images of the 3 vessel layers.Red indicates more superficial vessels while blue indicates deeper vessels.

#
187382 -$15.00USD Received 19 Mar 2013; revised 25 Apr 2013; accepted 27 Apr 2013; published 7 May 2013 (C) 2013 OSA 1 June 2013 | Vol. 4, No. 6 | DOI:10.1364/BOE.4.000803 | BIOMEDICAL OPTICS EXPRESS 815 original variance images as well as the results of the motion correction are shown.The expanded view in Fig. 11 shows the improvements in vasculature visualization gained from Gabor filtering and motion artifact correction.Once again, the distinct vasculature from the 3 main vessel beds can be observed and are also shown in a color encoded depth image.

Fig. 10 .
Fig. 10.Registration results of a second subject from a region nasal to the fovea.Three x-fast and three y-fast data sets were acquired.Each image covers a 2.5 × 2.5 mm scan area.

Fig. 11 .
Fig. 11.Enlarged view of images in Fig. 10 comparing the registered images to SVPs from a single volume.A color encoded depth image is shown on the right combining information from the 3 registered images.

Fig. 12 .
Fig. 12. Automated widefield mosaic showing vasculature of the ganglion cell plexus as generated from segmentation of the GCL + IPL.The individual motion corrected images are shown on the left.Representative B-scans from different volumes are shown corresponding to the dashed lines on the widefield mosaic.The nasal retina is on the left side of the mosaic, temporal on the right.The widefield image spans approximately a 6 × 4mm field of view.

Fig. 13 .
Fig. 13.Widefield mosaic of retinal layers showing vasculature of the superficial capillary plexus based on the segmentation of the IPL-INL junction.The individual motion corrected images are shown on the left.The nasal retina is on the left side of the mosaic, temporal on the right.

Fig. 14 .
Fig. 14.Widefield mosaic of retinal layers showing vasculature of the deep capillary plexus based on the segmentation of the INL-OPL junction.The individual motion corrected images are shown on the left.The nasal retina is on the left side of the mosaic, temporal on the right.

Fig. 15 .
Fig. 15.Widefield mosaic of retinal layers showing a color encoded depth image of the widefield mosaic with information from the registered mosaics of the 3 main vessel layers.Red indicates superficial vessels while blue indicates deeper vessels.The individual, color encoded depth, motion corrected images are shown on the left.The nasal retina is on the left side of the mosaic, temporal on the right.