Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Invariant features-based automated registration and montage for wide-field OCT angiography

Open Access Open Access

Abstract

The field of view of optical coherence tomography angiography (OCTA) images of the retina can be increased by montaging consecutive scans acquired at different retinal regions. Given the dramatic variation in aberrations throughout the entire posterior pole region, it is challenging to achieve seamless merging with apparent capillary continuity across the boundaries between adjacent angiograms. For this purpose, we propose herein a method that performs automated registration of contiguous OCTA images based on invariant features and uses a novel montage algorithm. The invariant features were used to register the overlapping areas between adjacently located scans by estimating the affine transformation matrix needed to accurately stitch them. Then, the flow signal was compensated to homogenize the angiograms with different brightness and the joints were blended to generate a seamless, montaged wide-field angiogram. We evaluated the algorithm on normal and diabetic retinopathy eyes. The proposed method could montage the angiograms seamlessly and provided a wide-field of view of retinal vasculature.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Since the first time-domain optical coherence tomography (OCT) [1] article was published in 1991, OCT technology has developed rapidly. With faster scanning using Fourier-domain modalities, current OCT devices allow acquisition of three-dimensional data representing the retinal layered structure rapidly. Optical coherence tomography angiography (OCTA) detects the variation of OCT signal reflectance over scans repeated at the same location in order to evaluate the microcirculation flow in vivo; allowing sensitive, fast (< 3 seconds), and high resolution volumetric imaging of the retinal flow [2]. This technique has demonstrated utility in key ophthalmic diseases such as diabetic retinopathy, age-related macular degeneration, and glaucoma, allowing visualization of vascular details. In addition, OCTA allows more reliable quantification of the vasculature changes [3]. With the improvement of the hardware and the progress of processing algorithms, a current research endeavor in OCTA is to increase the achievable field of view (FOV) of images. Previous studies have shown that wide-field OCTA can reveal more clinically useful information [4,5].

The availability of faster light sources has allowed acquisition of wide-field OCTA in a single scan [6,7]. This approach to increase the FOV faces two fundamental challenges. First, the retinal curvature covered by the large FOV and limited focusing depth of OCT result in B-scans that are only in focus centrally, and defocused at the periphery [8]. Secondly, propagation of the scanning beam across a larger FOV induces larger wavefront aberrations that are difficult to compensate for without the assistance of adaptive optics [9]. As a consequence, the beam size at the peripheral retina is large and aberrant, deteriorating the lateral resolution of the instrument and degrading the images. An alternative to wide-field scanning is to montage smaller, but high-quality scans with overlapping areas to create a wide-field image. This approach is more cost-effective and amenable to improvement using the existing equipment.

The main task of mosaicking techniques that reproduce wide-FOV from stitching small-FOV scans is to conceal the seam between individual angiograms by ensuring capillary overlap and background level homogeneity. Previous research on montaging retinal images used fundus images with large vessel patterns and intersections to register overlapping areas as well as post-processing of the edges in order to get a seamless, wider fundus image [10–16]. Accomplishing the same task in OCTA requires the use of more challenging image processing techniques, given the finer vessel details visualized by OCTA. Moreover, in order to attain a large FOV by mosaicking, the individual scans need to be large enough to not require an unreasonable number of acquisitions from the subject, making 3 × 3-mm scans impractical. However, the standard 6 × 6-mm angiograms are not isoplanatic. For that FOV, the retinal curvature in the periphery of posterior pole or in near-sighted eyes already produces a significant deformation, and further image processing is required to correct for an adequate blending with adjacent patches.

In our previous research montaging OCTA images, we developed a parallel-strip registration algorithm in 2D [17] and 3D [18]. This method focused on removing motion artifacts and required a large overlapping area. Here, we propose a new mosaicking solution in which patches of the retina with limited overlapping areas are registered by a method based on invariant features [19,20]. The registered patches are seamlessly stitched by a panoramic image seamless stitching routine [21–24] using flow signal compensation. This method allows ultra-wide visualization of the retinal flow by OCTA, with a great potential to better assist clinicians in the assessment of retinal diseases.

2. Data acquisition and preprocessing

We developed an algorithm to montage separate OCTA images in order to achieve a wide field of view. Although we demonstrated the algorithm mainly on scans centered on the optic disc and macula, it can be applied on the scans through to the back of the eyes. The mean distance between disc and fovea was 4.76 ± 0.34 mm [25], and the mean disc-fovea angle was 7.76 ± 3.63° [26]. Assuming the distance from the center of the FAZ to the optic disc is approximately 5 mm horizontally and about 1 mm vertically, if the two scans are perfectly centered the overlap between them would be approximately 1 mm horizontally and 5 mm vertically, with variations between subjects.

Participants were scanned using a 70-kHz commercial AngioVue OCTA system (RTVue-XR; Optovue, Fremont, CA) with the wavelength centered at 840nm. Scans covered a 6 × 6-mm area, using a high-definition scan pattern with 400 × 400 sampling rate (15um/pixel). Two repeated B-scans at the same location were acquired and processed by the commercial version of the split-spectrum amplitude-decorrelation angiography (SSADA), which is sensitive to capillary flow. A horizontal-priority and a vertical-priority scan of the same area were collected and registered to suppress microsaccadic motions artifacts [27]. Ten eyes diagnosed with diabetic retinopathy (DR) and 10 normal eyes were scanned at different retinal landmarks to generate the wide-field OCTA. The DR eyes were scanned centering on the disc, fovea and temporal regions, ensuring a certain overlap between adjacent scans. Normal eyes were scanned centering on the disc and fovea to evaluate the proposed algorithm.

During data preprocessing, retinal layers were segmented by a directional graph search algorithm incorporated in the COOL-ART software for OCTA image processing [28] (Fig. 1(A)). Applying the projection-resolved algorithm [29,30], we could remove projection artifacts and visualize the retinal capillary plexuses in three different slabs located at different depths [31] (Fig. 1(B-G)), which are defined as the superficial vascular complex (SVC), the intermediate capillary plexus (ICP), and the deep capillary plexus (DCP). The combination of three plexuses formed the inner retinal angiogram. Since the inner retinal angiogram is rich in microvascular details, it was used to detect the invariant features and to further estimate the optimal affine matrix to montage the en face angiograms of SVC, ICP, DCP, and the inner retina itself.

 figure: Fig. 1

Fig. 1 Illustration of retinal layer segmentation and angiogram generation on SVC (inner 80% of the ganglion cell complex (GCC)), ICP (outer 20% of the GCC and inner 50% of the inner nuclear layer (INL)), DCP (outer 50% of the INL and outer plexiform layer (OPL)), and inner retina (between the inner limiting membrane (ILM) and the OPL).

Download Full Size | PDF

3. Methods

There were three key steps to generate the montaged, wide-field OCTA: accurate estimation of the affine transformation matrix, flow signal compensation, and seamless blending of the edge of overlapping areas. The estimated affine transformation matrix contains the global information used to rotate and shift one image to montage it with the other. The fine registration in the edge area registered the vessels regionally in order to smooth the edges. The flow signal compensation equalized the flow signal to generate a homogenized wide-field OCTA.

3.1 Algorithm overview and data pre-processing

The disc angiogram was selected as the target image (It) and the macular angiogram was the moving image (Im). Then, the invariant features were detected and matched by Speed-Up Robust-Features (SURF) algorithm [32], and further fine matching in local area followed by Random Sampling Consensus (RANSAC) algorithm [33] were applied to filter the outliers. The remaining paired points were used to infer the affine matrix. Then, in order to generate an even and seamless wide-field angiogram, a post-processing step including flow signal compensation and seamless blending was applied. The algorithm framework is described in Fig. 2.

 figure: Fig. 2

Fig. 2 Overview of the invariant features-based automated registration and montage algorithm. SURF – Speed-Up Robust Features.

Download Full Size | PDF

3.2 Invariant features-based registration and montage

3.2.1 Affine transformation and features detection

Stitching the moving angiogram (Im) to the target angiogram (It) consists of re-arranging the pixels of the moving angiogram to match their equivalents in the overlapping area of the target angiogram. The transformation needed for 2D image stitching can be modeled by the FREAaffine transformation matrix:

p1=Ap0,A=[a11a12a13a21a22a23001]
where p0 and p1 are matching points in Im and It. The affine transformation matrix A includes the parameters for rotation, scaling, shifting, flipping and shear mapping. There are 6 degrees of freedom, hence at least 3 matching points are needed to infer the affine transformation matrix.

Scale-invariant feature transform (SIFT) [34] and SURF are popular local feature detection and descriptor methods in computer vision. Both are invariant with rotation, scaling, and intensity variation. Since SURF is more efficient and robust than SIFT, we chose SURF as the invariant features method to select the matching points used to estimate the affine transformation matrix.

SURF detects the features on points of interest in multi-scales, where the points of interest are detected by the determinant of the Hessian matrix following

det(H)=DxxDyy(0.9Dxy)2
where Dxy, Dxy and Dyy are the approximation second derivative of a Gaussian filter. To ensure the invariability of scaling and rotation, the SURF features were detected in multi-scales, and the dominant orientation of the Haar wavelet features were recorded as the SURF orientation descriptor. Furthermore, the summaries of the Haar wavelet features of the neighbors were used to generate the unique and robust SURF descriptors.

The number of SURF features that can be detected is variable, but usually about 800 of SURF features were detected on the OCTA of healthy eyes, which were presented using green crosses in the optic disc angiogram (Fig. 3(A)). Each feature contained location, orientation, scale, and feature descriptors (Fig. 3(B)). In the following step, the SURF descriptors were used to match the detected points, whereas the scale and orientation were used to filter the mismatched pairs.

 figure: Fig. 3

Fig. 3 SURF features detection on OCTAs. (A) The green crosses are the locations of detected features on the angiogram. (B) Some of the SURF features in A were selected randomly for clear representation of orientation and scale. The scale of features is represented by the radius of the green circles. The direction of features and scale is represented by the green line highlighted with white arrows in C and D. (C-D) Magnified insets correspond to the regions indicated by black and blue boxes in B.

Download Full Size | PDF

3.2.2 Angiograms registration and montage

Each detected invariant feature had their unique feature descriptor. The key points with same sign (-/+) of the det(H) would conduct the matching process. We compared the invariant features from It with their corresponding ones in Im using the sum of squared differences (SSD):

paired[Pt(i),Pm(j)]=minj{SSD[Pt(i),Pm(j)]}
where Pm(j) is the position of j-th SURF feature in Im, Pt(i) is the position of the i-th SURF feature in It. The nearest two points were paired.

OCTA contains a rich vasculature, the image content is monotonic, and the local features are similar. For this reason, if the points of interest were grossly matched based on the feature descriptors only (Eq. (2), many of them could be erroneously paired (Fig. 4). These mismatches would cause failure of the random sample consensus (RANSAC) algorithm in trying to infer the affine transformation matrix A. The RANSAC algorithm was designed to find a homograph matrix (H) that could get the greatest number of pairs that can be properly transformed. Therefore, before applying RANSAC, we imposed two additional restrictions to improve the matching accuracy.

 figure: Fig. 4

Fig. 4 Invariant features matching. The paired points were connected using the cyan dashed lines.

Download Full Size | PDF

For each pair, one point was taken from the moving image (Im) and represented by Pm whereas the other was taken from the target image (It) and represented by Pt. The SURF features of Pm and Pt were represented by (Om, Sm, Dm) and (Ot, St, Dt), where O is the orientation, S is the scale, and D is the feature descriptor. D was already used to grossly pair the SURF points. Then we used the orientation Om and scale Sm to locally rotate and resize the moving image and calculate the correlation with the target image in the paired local area (Fig. 5). The differences of the orientation and scale between the two points of interest were calculated by:

{ΔO=OmOtΔS=SmSt
If -π/8ΔOπ/8 and 1/2ΔS2, the points would be kept for the further descriptor comparison using Eq. (3). Otherwise, the two points would be marked with an “unmatchable” flag, assuming that scale and orientation should not have changed much between two adjacent scans. Some correctly matched and mismatched points were extracted from Fig. 4 to illustrate the differences of scale or orientation between mismatched points (Fig. 5).

 figure: Fig. 5

Fig. 5 Illustration of the scale and orientation of paired points. A1, A2 were correctly paired points of interest; their scale and direction were very similar. B1, B2 were mismatched points, for which both direction and scale were very different. C1, C2 were also mismatched points, where the scale is similar but the difference of direction was great. D is the illustration for the scale and direction restriction for feature matching. The scales of A, B C and D are 1.65 × 1.65-mm, 1.3 × 1.3-mm, 1.65 × 1.65-mm and 0.9 × 0.9-mm, separately.

Download Full Size | PDF

After the first restriction in Eq. (4) is imposed, and for each remaining pair of points, two rectangular areas of size 21 × 21-pixel were extracted from the target and moving images, surrounding the pixels in the pair. The local areas Il surrounding the kth pair were represented by Ilm(k) and Ilt(k). Then, Ilm(k) was resized and rotated to register with the Ilt(k), following

Icl=T[R(Iml,ΔO),ΔS1]
where R is the rotation operation, T is the resize operation, ΔO is the rotation parameter, and ΔS is the rescaling parameter. Then, the second restriction was imposed, consisting of evaluating the correlation of the local areas by:
Cl=Cor(Icl(k),Itl(k))
The areas surrounding truly matching points exhibited high correlation when aligned. The pairs with correlation Clower than a threshold value of 0.5 were also filtered out. By doing this, the SURF descriptors and image content of highly matching points were reserved (Fig. 6).

 figure: Fig. 6

Fig. 6 Illustration of the pairs reserved for the RANSAC input after local area verification by imposition of the two restrictions described in section 3.2.2.

Download Full Size | PDF

After most outliers had been filtered out (Fig. 6), the RANSAC algorithm could be employed on reserved pairs to further remove the remaining mismatching pairs (Fig. 7) and then recognize the registration matrix. The number of available pairs always exceeded the minimum of 3 points required by RANSAC. The least squares approach was used to infer the affine transformation matrix A At this point, the moving angiogram (Im) was well montaged onto the target angiogram (It), but the whole angiogram exhibited an uneven appearance and the seam was visible (Fig. 8).

 figure: Fig. 7

Fig. 7 Illustration of the pairs reserved after the last step of mismatch removal using random sample consensus (RANSAC) algorithm.

Download Full Size | PDF

 figure: Fig. 8

Fig. 8 Angiography montage with visible seams

Download Full Size | PDF

3.3 Flow signal compensation and seamless blending

Aiming at homogenizing the wide-field angiogram, post-processing consisting of flow signal compensation and seamless blending was performed.

Two methods were used to compensate the flow signal: structural reflectance-based compensation and flow signal-based compensation.

The structural reflectance was obtained from the OCT en face projection and used a bias field map to compensate the angiogram:

Ic(A)=IAMean(G(IR))G(IR)
where, IA is the original angiogram image, IR is the en face mean projection of the structural reflectance, and G(·) is the Gaussian filter. The intensity corrected angiogram Ic(A) reduced the effect of structural reflectance variation on flow signal.

Then, flow signal based compensation was applied, inspired by the gain compensation algorithm described in [24]. This task was accomplished by minimizing an error function:

e=1N[αIc(m)(um)βIc(t)(ut)]2
where N is the number of overlapping pixels, um and ut are the overlapping pixels in Im and If with um = Aut, and α and β are gain parameters that need to be optimized. In our application, the α and β were restricted, which were both no less than 0.5 and no more than 2.

After compensation, the flow signal exhibited global uniformity. However, the transition of the flow signal from one angiogram to the other was discontinuous and the vasculatures were not perfectly registered at the seams, and hence the edge was visible. Flow signal blending was necessary to obtain a seamless montage angiogram. The multi-band blending algorithm described in [24] was simplified and applied in this study. The flow signal was fused in multi-scale with a fusion weights map decreasing gradually from 1 to 0 from one angiogram to the other across the seam (Fig. 9). Weights maps for the target It and the moving Im images were represented by Wt and Wm. The center of the weights map was located on the seam, where the weight is 0.5 on both weights maps. The angiograms were filtered by 2D Gaussian kernels with multiple scales; the high-frequency information was extracted by

{Diff(k)=Gσ(k)(I)Gσ(k+1)(I)σ(k+1)=2*σ(k)
When k = 0, we assign to Gσ(k)(I) the original image I

 figure: Fig. 9

Fig. 9 Illustration of the seamless blending routine.

Download Full Size | PDF

The information on the scale was blended by

{Ib(k)=Wt×Diff(tk)+Wm×Diff(mk)Wt+Wm=1
where Ib(k) is the blended image on scale k, Wt is the weights map for the target angiogram, and Wm is the weights map for moving angiogram. The weights maps vary within [0, 1] from the one side of the seam to other. The center of the weights map is 0.5, aligning with the seam. The weights maps were illustrated on Fig. 9. Diff(tk) and Diff(mk) are the detected Diff values in Eq. (9) for the target and moving angiograms separately.

On all scales, the blended images were summed up to get the seamless wide-field angiogram:

Ioutput=k=0NIb(k)
Where N is the number of scales, and Ioutput is final montaged seamless wide-field angiogram. In this study, N = 4.

Comparing the montaged image after seamless blending (Fig. 10(A)) with the montaged angiogram without flow signal compensation in Fig. 8, the transition of the first seems smoother. By inspecting local areas about the seam, the vasculature exhibited a continuous appearance across the seam (Fig. 10(B-C)).

 figure: Fig. 10

Fig. 10 The montaged seamless angiogram applying flow signal compensation and seamless blending. (B) and (C) represent the detail capillary network within the green and blue rectangles in (A).

Download Full Size | PDF

3.4 Bundle processing

As mentioned previously, this invariant features-based automated registration and montage algorithm is not limited to two angiograms only. For the merging of multiple angiograms we designed the bundle processing strategy described below.

First, a target angiogram needs to be pre-assigned and used as the reference for the montage; all other moving angiograms would be stitched to it one by one. Utilizing the proposed registration method in section 3.2, the invariant features of all input angiograms were paired each other. The image connection would be sketched out based on the cluster of the paired points. The pairs only appeared in the overlapping area and were located in a restricted area. The cluster of the paired points in the moving image closer to the target image had priority to be montaged. Once the montage was completed, the resulting angiogram would be updated as the target image, to be further montaged with the following moving angiogram.

4. Results

This proposed method was implemented using Matlab 2018a (Mathworks, Natick, MA) installed on a computer with Intel (R) Core(TM) i7-6800K CPU @ 3.4GHz and DDR4 64GB. It could generate a seamless wide-field OCTA by stitching two OCTAs within 15 seconds.

To evaluate the montage accuracy and the continuity of vasculature across the seams, the difference and correlation of flow signal between opposite sides of the seams were calculated on the angiography montage of 10 normal eyes. For evaluation, disc and macular inner retinal angiograms were montaged. As illustrated in Fig. 11, the green and cyan thin-slabs contained information from the two sets of flow signals (Dt and Dm) on each side of the seam (dotted white line). The difference (E(Dt, Dm)) and correlation (Cor(Dt, Dm)) between Dt and Dm were calculated by

{E(Dt,Dm)=1Li=1L|Mt(i)Mm(i)|Cor(Dt,Dm)=cor(Mt,Mm)Mt(i)=1Wtj=1Wt(Dt(i,j)),Mm(i)=1Wmj=1Wm(Dm(i,j))
where, L is the length of the seams, cor(·) is the correlation calculation, Wt = Wm = 2 pixels are the width of the slab. Mt(i) and Mm(i) are the flow values of the i-th pixel on the target portion and moving portion, separately. While Mt(i) and Mm(i) are paired on the opposite sides of the i-th position on the seam. The difference evaluates the smoothness of angiogram on the mosaic edge, while correlation evaluates the vasculature continuity.

 figure: Fig. 11

Fig. 11 Illustration of the algorithm evaluation. Green and cyan lines illustrated both sides of the seam between the target and moving images.

Download Full Size | PDF

For comparison, we compared the proposed algorithm with the traditional montage method. Previously, the matched points were obtained based on the invariant feature descriptors and removing the outliers using RANSAC, which means the angiograms were montaged without fine registration in local area, flow signal compensation and seamless blending. Using the traditional approach, the angiograms appeared well montaged (Fig. 12 (A)) but vessels were not continuous across the seams (Fig. 12 (B)). In contrast, after applying the proposed method, the montaged angiogram (Fig. 10 (A)) was more homogeneous, the difference and correlation values were improved (Table 1), the vessels were continuous, and the seam was not visible (Fig. 10 vs Fig. 12). The vessels were not absolutely perpendicular to the seams, and the vessels patterns were not completely the same at opposite sides of the seam. Taking these two main factors, which are consequence of the scale over which correlation is computed and the amount of detail available, the correlation of the proposed method cannot be particularly high anyways. We would expect it to be considerably higher if the overlapping area was larger. However, the fact that the correlation is improved with respect to the traditional method indicates that vessel connectivity is improved at the stitched area.

 figure: Fig. 12

Fig. 12 Grossly montaged wide-field angiogram before matching verification in local area, flow signal compensation and seamless blending. Five pixel mismatch was demonstrated.

Download Full Size | PDF

Tables Icon

Table 1. Evaluation of the performance of the automated montage algorithm on 10 normal eyes by the correlation and mean square errors

The affine transformation matrix estimated for the inner retinal angiogram could be used to generate the wide-field angiograms on SVC, ICP, DCP and choriocapillaris (Fig. 13) after projection artifacts were removed by PR-OCTA.

 figure: Fig. 13

Fig. 13 Automated montaged angiogram on four slabs of PR-OCT angiography (2 scans). (A) Superficial vascular complex. (B) Intermediate capillary plexus. (C) Deep capillary plexus. (D) Choriocapillaris.

Download Full Size | PDF

Next, sixteen 6 × 6-mm HD scans were obtained on a normal eye covering an ultra wide-field of view [Fig. 14] and reaching the limit imposed by dilated pupil vignetting. Using the algorithm proposed here, the montaged ultra wide-field angiogram was 15 × 20-mm, and the sampling density was 15 um/pixel. Overviewing the angiogram, the vasculatures were continuous and well connected, which verified the robustness of the bundle of angiograms seamless stitching.

 figure: Fig. 14

Fig. 14 Ultra wide-field OCTA (15 × 20-mm) of a healthy eye registered and montaged by sixteen 6 × 6-mm OCTA scans.

Download Full Size | PDF

For DR cases (Fig. 15), 3 scans centered on the disc, macular, and temporal side of the HD scan pattern were montaged using the proposed algorithm. The montaged wide-field angiogram provided details of the capillary nonperfusion, demonstrating greater extent of abnormalities in the temporal field. The individual scans could not provide equivalent detailed contextualized information for evaluating the progression of vasculopathy.

 figure: Fig. 15

Fig. 15 Montaged wide-field OCTA (6 × 15-mm,3 scans) of a diabetic retinopathy case. Three scans centered on disc, fovea, and temporal were montaged using the algorithm described in the text.

Download Full Size | PDF

This method was also applied on the images acquired from DR patients on a 200-kHz prototype swept source OCT angiography system, with wavelength centered at 1050nm. Four 10 × 8-mm scans, each with 800 × 400 A-scans, centered on the nasal, disc, macular, and temporal areas were merged. Figure 16 demonstrates the four scans stitched seamlessly. The montaged wide-field angiogram presented the detailed pathology, including both capillary loss in nasal and temporal areas as well as the retinal neovascularization.

 figure: Fig. 16

Fig. 16 Montaged OCTA (10 × 25-mm, 4 scans) of the nasal, disc, macula and temporal regions from an image of a diabetic retinopathy case acquired by a prototype OCT angiography system. The white arrow highlights the vitreous neovascularization.

Download Full Size | PDF

5. Discussion and conclusion

In this study, we developed an invariant features-based automated registration and montage algorithm for wide-field OCT angiography, composed of four main steps: (1) detect the invariant features including points of interest and their feature descriptors using the SURF algorithm; (2) match the detected points of interest of the target image with the moving image using their feature descriptors; (3) filter the inaccurate matches and apply the RANSAC algorithm; (4) estimate the affine matrix to stitch the moving image on the target image. The results show that the images were stitched seamlessly and ultra-wide-field OCTAs could be generated.

Clinical studies have validated the benefit of wide-field OCTA over regular field of view [35]. Wide-field OCTA allows clinicians to observe the overall vascular abnormalities, including capillary dropout or neovascularization in the peripheral retina, preserving visualization of the optic disc and macula. To achieve large fields of view by hardware improvements, dynamic compensation of the defocusing on the peripheral retina by adaptive optics would be necessary. The blending could generate a wide FOV angiogram with a very even appearance. For individual scans, each had different brightness. By presenting them with an even appearance, we don’t suggest that some areas have artefactual lower flow speed just because their capillaries show up with lower decorrelation values. While this has been demonstrated in single-scan, wide-field OCT, it has not been accomplished for OCTA thus far. Currently, commercial systems like the RT-Vue XR OCTA have incorporated software solutions to increase the field of view, such as montage of 6 × 6-mm angiograms of the macula and optic disc only, without the processing necessary to eliminate the apparent seams. To the best of our knowledge, no commercial software has incorporated a seamless montage algorithm for generation of clean wide-field angiograms.

Among the feature descriptor algorithms published in the past decade, e.g. BRISK, FREAK, SURF, KAZE, ORB – SURF descriptors generated using the Haar wavelet were more appropriate to describe vasculature patterns. For invariant features detection, both SIFT and SURF are very popular algorithms because they can locate the points of interest and generate the feature descriptors based on the neighborhood. In practice, SIFT could had also been used to perform the task assigned to SURF in our method. As a prior comparison has indicated [36], SIFT is more robust to the scaling and rotation changes, while SURF is more efficient and performs better on illumination changes. In our application, all scans were well focused and acquired over the same area, minimizing changes owing to scaling and rotation between adjacent scans. On the other hand, the strength of the flow signal is intensity dependent and would change considerably owing to reflectance variations between scans. As Fig. 8 demonstrated, the strength of the flow signals on the overlapping area between two scans can produce a very visible seam. Taking into account the factors above, SURF was selected to detect the invariant features.

In the invariant features-based image registration task, RANSAC algorithm is usually used to remove the outliers, but it fails when the number of outliers exceeds the inliers. In this application, the images are composed of very rich vascular patterns and the image features are highly similar, causing plenty of mismatches (as in Fig. 4). In this case, when we set an iteration limitation, it is easy to reach to a wrong solution. To address this problem, we did a further verification of the matches utilizing their scale, rotation and the correlation between the surrounding areas. The number of inaccurate matches was greatly reduced, and only the highly correlated matches were reserved (Fig. 6). Then RANSAC could remove the remaining outliers to estimate the accurate affine matrix used to stitch the moving image onto the target image.

There are a few limitations of the proposed automated seamless montage algorithm. An overlapping area is required for registration to estimate affine matrix. The overlapping area of the acquired data was 5 × 1-mm, but this is not the limitation of the proposed method. For the success montage, the overlapping area could be as small as 2 × 1-mm on normal eye, which was estimated by cropping the images to change the overlapping area. The other limitation is that at least 3 matches were required for the affine transformation estimation. Consequently, an inefficient number of angiograms need to be acquired from the subject in order to generate the wide-field view. For images with low quality on which pair key points could not be detected to match the features automatically, we can manually remove the mismatches or add matches to improve the accuracy of the estimated affine transformation. In this instance, this method would be a semi-automated, the further blending process could also be applied to generate the seamless WF montaged angiogram.

In summary, this invariant features-based automated montage algorithm can stitch multiple angiograms to generate a seamless wide-field angiogram, which provides an efficient way to view the vascular changes in a wide FOV without requiring challenging hardware updates or invasive imaging methods using contrast dye.

Funding

National Institutes of Health (R01 EY027833, R01 EY024544, R01 EY023285, DP3 DK104397, P30 EY010572); unrestricted departmental funding grant and William & Mary Greve Special Scholar Award from Research to Prevent Blindness (New York, NY).

Acknowledgments

We would like to thank Dr. Tristan Hormel for his constructive edit on this paper.

Disclosures

Oregon Health & Science University (OHSU), David Huang and Yali Jia, have a significant financial interest in Optovue, Inc. These potential conflicts of interest have been reviewed and managed by OHSU.

References

1. D. Huang, E. A. Swanson, C. P. Lin, J. S. Schuman, W. G. Stinson, W. Chang, M. R. Hee, T. Flotte, K. Gregory, C. A. Puliafito, and et, “Optical coherence tomography,” Science 254(5035), 1178–1181 (1991). [CrossRef]   [PubMed]  

2. Y. Jia, O. Tan, J. Tokayer, B. Potsaid, Y. Wang, J. J. Liu, M. F. Kraus, H. Subhash, J. G. Fujimoto, J. Hornegger, and D. Huang, “Split-spectrum amplitude-decorrelation angiography with optical coherence tomography,” Opt. Express 20(4), 4710–4725 (2012). [CrossRef]   [PubMed]  

3. Y. Jia, S. T. Bailey, T. S. Hwang, S. M. McClintic, S. S. Gao, M. E. Pennesi, C. J. Flaxel, A. K. Lauer, D. J. Wilson, J. Hornegger, J. G. Fujimoto, and D. Huang, “Quantitative optical coherence tomography angiography of vascular abnormalities in the living human eye,” Proc. Natl. Acad. Sci. U.S.A. 112(18), E2395–E2402 (2015). [CrossRef]   [PubMed]  

4. Q. Zhang, C. S. Lee, J. Chao, C.-L. Chen, T. Zhang, U. Sharma, A. Zhang, J. Liu, K. Rezaei, K. L. Pepple, R. Munsen, J. Kinyoun, M. Johnstone, R. N. Van Gelder, and R. K. Wang, “Wide-field optical coherence tomography based microangiography for retinal imaging,” Sci. Rep. 6(1), 22017 (2016). [CrossRef]   [PubMed]  

5. Q. Zhang, Y. Huang, T. Zhang, S. Kubach, L. An, M. Laron, U. Sharma, and R. K. Wang, “Wide-field imaging of retinal vasculature using optical coherence tomography-based microangiography provided by motion tracking,” J. Biomed. Opt. 20(6), 066008 (2015). [CrossRef]   [PubMed]  

6. T. Hirano, S. Kakihara, Y. Toriyama, M. G. Nittala, T. Murata, and S. Sadda, “Wide-field en face swept-source optical coherence tomography angiography using extended field imaging in diabetic retinopathy,” Br. J. Ophthalmol. 2017, 311358 (2017). [PubMed]  

7. M. Kimura, M. Nozaki, M. Yoshida, and Y. Ogura, “Wide-field optical coherence tomography angiography using extended field imaging technique to evaluate the nonperfusion area in retinal vein occlusion,” Clin. Ophthalmol. 10, 1291–1295 (2016). [CrossRef]   [PubMed]  

8. J. Polans, B. Keller, O. M. Carrasco-Zevallos, F. LaRocca, E. Cole, H. E. Whitson, E. M. Lad, S. Farsiu, and J. A. Izatt, “Wide-field retinal optical coherence tomography with wavefront sensorless adaptive optics for enhanced imaging of targeted regions,” Biomed. Opt. Express 8(1), 16–37 (2016). [CrossRef]   [PubMed]  

9. J. Polans, B. Jaeken, R. P. McNabb, P. Artal, and J. A. Izatt, “Wide-field optical model of the human eye with asymmetrically tilted and decentered lens that reproduces measured ocular aberrations,” Optica 2(2), 124–134 (2015). [CrossRef]  

10. A. Can, C. V. Stewart, B. Roysam, and H. L. Tanenbaum, “A feature-based, robust, hierarchical algorithm for registering pairs of images of the curved human retina,” IEEE Trans. Pattern Anal. Mach. Intell. 24(3), 347–364 (2002). [CrossRef]  

11. A. Can, C. V. Stewart, B. Roysam, and H. L. Tanenbaum, “A feature-based technique for joint, linear estimation of high-order image-to-mosaic transformations: mosaicing the curved human retina,” IEEE Trans. Pattern Anal. Mach. Intell. 24(3), 412–419 (2002). [CrossRef]  

12. A. Can, C. V. Stewart, and B. Roysam, “Robust hierarchical algorithm for constructing a mosaic from images of the curved human retina,” in Computer Vision and Pattern Recognition, 1999. IEEE Computer Society Conference on., (IEEE, 1999)

13. S. Lee, M. D. Abràmoff, and J. M. Reinhardt, “Validation of retinal image registration algorithms by a projective imaging distortion model,” in Engineering in Medicine and Biology Society,2007. EMBS 2007. 29th Annual International Conference of the IEEE, (IEEE, 2007), 6471–6474. [CrossRef]  

14. S. Lee, M. D. Abràmoff, and J. M. Reinhardt, “Feature-based pairwise retinal image registration by radial distortion correction,” in Medical Imaging 2007: Image Processing, (International Society for Optics and Photonics, 2007), 651220.

15. C. V. Stewart, C.-L. Tsai, and B. Roysam, “The dual-bootstrap iterative closest point algorithm with application to retinal image registration,” IEEE Trans. Med. Imaging 22(11), 1379–1394 (2003). [CrossRef]   [PubMed]  

16. J. Chen, J. Tian, N. Lee, J. Zheng, R. T. Smith, and A. F. Laine, “A partial intensity invariant feature descriptor for multimodal retinal image registration,” IEEE Trans. Biomed. Eng. 57(7), 1707–1718 (2010). [CrossRef]   [PubMed]  

17. P. Zang, G. Liu, M. Zhang, C. Dongye, J. Wang, A. D. Pechauer, T. S. Hwang, D. J. Wilson, D. Huang, D. Li, and Y. Jia, “Automated motion correction using parallel-strip registration for wide-field en face OCT angiogram,” Biomed. Opt. Express 7(7), 2823–2836 (2016). [CrossRef]   [PubMed]  

18. P. Zang, G. Liu, M. Zhang, J. Wang, T. S. Hwang, D. J. Wilson, D. Huang, D. Li, and Y. Jia, “Automated three-dimensional registration and volume rebuilding for wide-field angiographic and structural optical coherence tomography,” J. Biomed. Opt. 22(2), 26001 (2017). [CrossRef]   [PubMed]  

19. P. Lukashevich, B. Zalesky, and S. Ablameyko, “Medical image registration based on SURF detector,” Pattern Recognit. Image Anal. 21(3), 519–521 (2011). [CrossRef]  

20. M. Teke and A. Temizel, “Multi-spectral satellite image registration using scale-restricted SURF,” in Pattern Recognition (ICPR),201020th International Conference on, (IEEE, 2010), 2310–2313. [CrossRef]  

21. A. Eden, M. Uyttendaele, and R. Szeliski, “Seamless image stitching of scenes with large motions and exposure differences,” in Computer Vision and Pattern Recognition, 2006 IEEE Computer Society Conference on, (IEEE, 2006), 2498–2505. [CrossRef]  

22. A. Levin, A. Zomet, S. Peleg, and Y. Weiss, “Seamless image stitching in the gradient domain,” in European Conference on Computer Vision, (Springer, 2004), 377–389. [CrossRef]  

23. L. Juan and G. Oubong, “SURF applied in panorama image stitching,” in Image Processing Theory Tools and Applications (IPTA),20102nd International Conference on, (IEEE, 2010), 495–499. [CrossRef]  

24. M. Brown and D. G. Lowe, “Automatic panoramic image stitching using invariant features,” Int. J. Comput. Vis. 74(1), 59–73 (2007). [CrossRef]  

25. R. A. Jonas, Y. X. Wang, H. Yang, J. J. Li, L. Xu, S. Panda-Jonas, and J. B. Jonas, “Optic disc-fovea distance, axial length and parapapillary zones. The Beijing Eye Study 2011,” PLoS One 10(9), e0138701 (2015). [CrossRef]   [PubMed]  

26. R. A. Jonas, Y. X. Wang, H. Yang, J. J. Li, L. Xu, S. Panda-Jonas, and J. B. Jonas, “Optic disc-fovea angle: the Beijing Eye Study 2011,” PLoS One 10(11), e0141771 (2015). [CrossRef]   [PubMed]  

27. M. F. Kraus, B. Potsaid, M. A. Mayer, R. Bock, B. Baumann, J. J. Liu, J. Hornegger, and J. G. Fujimoto, “Motion correction in optical coherence tomography volumes on a per A-scan basis using orthogonal scan patterns,” Biomed. Opt. Express 3(6), 1182–1199 (2012). [CrossRef]   [PubMed]  

28. M. Zhang, J. Wang, A. D. Pechauer, T. S. Hwang, S. S. Gao, L. Liu, L. Liu, S. T. Bailey, D. J. Wilson, D. Huang, and Y. Jia, “Advanced image processing for optical coherence tomographic angiography of macular diseases,” Biomed. Opt. Express 6(12), 4661–4675 (2015). [CrossRef]   [PubMed]  

29. J. Wang, M. Zhang, T. S. Hwang, S. T. Bailey, D. Huang, D. J. Wilson, and Y. Jia, “Reflectance-based projection-resolved optical coherence tomography angiography [Invited],” Biomed. Opt. Express 8(3), 1536–1548 (2017). [CrossRef]   [PubMed]  

30. M. Zhang, T. S. Hwang, J. P. Campbell, S. T. Bailey, D. J. Wilson, D. Huang, and Y. Jia, “Projection-resolved optical coherence tomographic angiography,” Biomed. Opt. Express 7(3), 816–828 (2016). [CrossRef]   [PubMed]  

31. J. P. Campbell, M. Zhang, T. S. Hwang, S. T. Bailey, D. J. Wilson, Y. Jia, and D. Huang, “Detailed vascular anatomy of the human retina by projection-resolved optical coherence tomography angiography,” Sci. Rep. 7(1), 42201 (2017). [CrossRef]   [PubMed]  

32. H. Bay, A. Ess, T. Tuytelaars, and L. Van Gool, “Speeded-up robust features (SURF),” Comput. Vis. Image Underst. 110(3), 346–359 (2008). [CrossRef]  

33. M. A. Fischler and R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography,” Commun. ACM 24(6), 381–395 (1981). [CrossRef]  

34. D. G. Lowe, “Distinctive image features from scale-invariant keypoints,” Int. J. Comput. Vis. 60(2), 91–110 (2004). [CrossRef]  

35. Q. Zhang, C.-L. Chen, Z. Chu, and R. K. Wang, “Wide field OCT angiography by using swept source OCT in living human eye,” in Ophthalmic Technologies XXVII, (International Society for Optics and Photonics, 2017), 100450V.

36. L. Juan and O. Gwun, “A comparison of sift, pca-sift and surf,” International Journal of Image Processing 3, 143–152 (2009).

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (16)

Fig. 1
Fig. 1 Illustration of retinal layer segmentation and angiogram generation on SVC (inner 80% of the ganglion cell complex (GCC)), ICP (outer 20% of the GCC and inner 50% of the inner nuclear layer (INL)), DCP (outer 50% of the INL and outer plexiform layer (OPL)), and inner retina (between the inner limiting membrane (ILM) and the OPL).
Fig. 2
Fig. 2 Overview of the invariant features-based automated registration and montage algorithm. SURF – Speed-Up Robust Features.
Fig. 3
Fig. 3 SURF features detection on OCTAs. (A) The green crosses are the locations of detected features on the angiogram. (B) Some of the SURF features in A were selected randomly for clear representation of orientation and scale. The scale of features is represented by the radius of the green circles. The direction of features and scale is represented by the green line highlighted with white arrows in C and D. (C-D) Magnified insets correspond to the regions indicated by black and blue boxes in B.
Fig. 4
Fig. 4 Invariant features matching. The paired points were connected using the cyan dashed lines.
Fig. 5
Fig. 5 Illustration of the scale and orientation of paired points. A1, A2 were correctly paired points of interest; their scale and direction were very similar. B1, B2 were mismatched points, for which both direction and scale were very different. C1, C2 were also mismatched points, where the scale is similar but the difference of direction was great. D is the illustration for the scale and direction restriction for feature matching. The scales of A, B C and D are 1.65 × 1.65-mm, 1.3 × 1.3-mm, 1.65 × 1.65-mm and 0.9 × 0.9-mm, separately.
Fig. 6
Fig. 6 Illustration of the pairs reserved for the RANSAC input after local area verification by imposition of the two restrictions described in section 3.2.2.
Fig. 7
Fig. 7 Illustration of the pairs reserved after the last step of mismatch removal using random sample consensus (RANSAC) algorithm.
Fig. 8
Fig. 8 Angiography montage with visible seams
Fig. 9
Fig. 9 Illustration of the seamless blending routine.
Fig. 10
Fig. 10 The montaged seamless angiogram applying flow signal compensation and seamless blending. (B) and (C) represent the detail capillary network within the green and blue rectangles in (A).
Fig. 11
Fig. 11 Illustration of the algorithm evaluation. Green and cyan lines illustrated both sides of the seam between the target and moving images.
Fig. 12
Fig. 12 Grossly montaged wide-field angiogram before matching verification in local area, flow signal compensation and seamless blending. Five pixel mismatch was demonstrated.
Fig. 13
Fig. 13 Automated montaged angiogram on four slabs of PR-OCT angiography (2 scans). (A) Superficial vascular complex. (B) Intermediate capillary plexus. (C) Deep capillary plexus. (D) Choriocapillaris.
Fig. 14
Fig. 14 Ultra wide-field OCTA (15 × 20-mm) of a healthy eye registered and montaged by sixteen 6 × 6-mm OCTA scans.
Fig. 15
Fig. 15 Montaged wide-field OCTA (6 × 15-mm,3 scans) of a diabetic retinopathy case. Three scans centered on disc, fovea, and temporal were montaged using the algorithm described in the text.
Fig. 16
Fig. 16 Montaged OCTA (10 × 25-mm, 4 scans) of the nasal, disc, macula and temporal regions from an image of a diabetic retinopathy case acquired by a prototype OCT angiography system. The white arrow highlights the vitreous neovascularization.

Tables (1)

Tables Icon

Table 1 Evaluation of the performance of the automated montage algorithm on 10 normal eyes by the correlation and mean square errors

Equations (12)

Equations on this page are rendered with MathJax. Learn more.

p 1 =A p 0 ,A=[ a 11 a 12 a 13 a 21 a 22 a 23 0 0 1 ]
det(H)= D xx D yy ( 0.9 D xy ) 2
paired[ P t ( i ), P m ( j ) ]= min j { SSD[ P t ( i ), P m ( j ) ] }
{ ΔO= O m O t ΔS= S m S t
I c l =T[ R( I m l ,ΔO ),Δ S 1 ]
C l =Cor( I c l (k), I t l (k) )
I c(A) = I A Mean( G( I R ) ) G( I R )
e= 1 N [ α I c(m) ( u m )β I c(t) ( u t ) ] 2
{ Diff( k )= G σ(k) (I) G σ(k+1) (I) σ(k+1)=2*σ(k)
{ I b ( k )= W t ×Diff( t k )+ W m ×Diff( m k ) W t + W m =1
I output = k=0 N I b ( k )
{ E( D t , D m )= 1 L i=1 L | M t ( i ) M m ( i ) | Cor( D t , D m )=cor( M t , M m ) M t ( i )= 1 W t j=1 W t ( D t (i,j) ) , M m ( i )= 1 W m j=1 W m ( D m (i,j) )
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.