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

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.

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][11][12][13][14][15][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][22][23][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.

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 × 6mm 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. 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).

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.

Algorithm overview and data pre-processing
The disc angiogram was selected as the target image (I t ) and the macular angiogram was the moving image (I m ). 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.

Affine transformation and features detection
Stitching the moving angiogram (I m ) to the target angiogram (I t ) consists of re-arranging the pixels of the moving angiogram to match their equivalents in the overlapping area of the

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 I t with their corresponding ones in I m using the sum of squared differences (SSD): where P m (j) is the position of j-th SURF feature in I m , P t (i) is the position of the i-th SURF feature in I t . 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.
For each pair, one point was taken from the moving image (I m ) and represented by P m whereas the other was taken from the target image (I t ) and represented by P t . The SURF features of P m and P t were represented by 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 O m and scale S m 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: If -π/8≤ΔO ≤π/8 and 1/2≤ΔS≤2, 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).

Flow sign
Aiming at hom compensation Two meth compensation The struct field map to c where, I A is th reflectance, an effect of struc Then, flow algorithm des where N is th with u m = Au the α and β w After com of the flow si were not per blending was algorithm des fused in mult angiogram to I m images wer F nal compensa mogenizing th n and seamless hods were use n and flow sign tural reflectanc compensate the he original ang nd G(·) is the G ctural reflectan w signal based scribed in [24]. e number of ov t , and α and β ere restricted, w mpensation, the ignal from one rfectly register necessary to scribed in [24] ti-scale with a the other acros re represented n I m and I f plication, transition sculatures ow signal blending ignal was from one he moving ted 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 When k = 0, we assign to G σ(k) (I) the original image I The information on the scale was blended by ( ) where I b (k) is the blended image on scale k, W t is the weights map for the target angiogram, and W m 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(t k ) and Diff(m k ) 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: Where N is the number of scales, and I output 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)).

Bundle p
As mentioned algorithm is n designed the b First, a ta montage; all proposed regi were paired e the paired po restricted area had priority t would be upd angiogram.

This propose installed on a It could gener
To evalua difference and on the angiog angiograms w information fr white line). T calculated by 10 Fig. 12) patterns were n factors, which nt of detail avai nyways. We w However, the indicates that matching verificat el mismatch was de on, W t = W m = i-th pixel on t paired on the moothness of an uity. rated both sides of he traditional n the invarian ns the angiogra pensation and well montaged (B)). In contr . 10 (A)) w (  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.