Advanced image processing for optical coherence tomographic angiography of macular diseases

,


Introduction
Optical coherence tomography (OCT) provides cross-sectional and three-dimensional (3D) imaging of biological tissues, and is now a part of the standard of care in ophthalmology [1,2]. Conventional OCT, however, is only sensitive to backscattered light intensity and is unable to directly detect blood flow and vascular abnormalities such as capillary dropout or pathologic vessel growth (neovascularization), which are the major vascular abnormalities associated with two of the leading causes of blindness, age-related macular degeneration (AMD) and proliferative diabetic retinopathy (PDR) [3]. Current techniques that visualize these abnormalities require an intravenous dye-based contrast such as fluorescein angiography (FA) or indocyanine green (ICG) angiography.
OCT angiography uses the motion of red blood cells against static tissue as intrinsic contrast. This approach eliminates the risk and reduces the time associated with dye injections [4,5], making it more accessible for clinical use than FA or ICG. A novel 3D OCT angiography technique called split-spectrum amplitude-decorrelation angiography (SSADA) can detect motion-related amplitude-decorrelation on commercially available OCT machines. Using this algorithm, the contrast between static and non-static tissue enables visualization of blood flow, providing high resolution maps of microvascular networks in addition to the conventional structural OCT images [6,7]. En face projection of the maximum decorrelation within anatomic layers (slabs) can produce angiograms analogous to traditional FA and ICG angiography [5,8].
Applying SSADA-based OCT angiography, we and others have quantified vessel density and flow index [9][10][11][12], choroidal neovascularization (CNV) area [4,13], and detected retinal neovascularization (RNV) [12,14] and macular ischemia. Accurate segmentation is necessary for interpretation and quantification of 3D angiograms. However, in the diseased eye, pathologies such as drusen, cystoid macular edema, subretinal fluid, or pigment epithelial detachment distort the normal tissue boundaries. Such distortion increases the difficulty of automated slab boundary segmentation. Although researchers have been working on improving automated segmentation in pathological retina [15][16][17], there is still no fully automated method which guarantees success in all clinical cases and hence manual segmentation or correction is often required. Previously reported manual correction of segmentation is tedious and inefficient [18,19]. In this manuscript, we provide an overview of our advanced image processing of SSADA-based OCT angiography, introduce an automated layer segmentation algorithm with expert correction, which is able to efficiently handle all clinical cases and show results of our technique applied to the processing of OCT angiograms of diseased eyes.

OCT angiography data acquisition
The OCT angiography data was acquired using a commercial spectral domain OCT instrument (RTVue-XR; Optovue). It has a center wavelength of 840 nm with a full-width half-maximum bandwidth of 45 nm and an axial scan rate of 70 kHz. Volumetric macular scans consisted of a 3 × 3 mm or 6 × 6 mm area with a 1.6 mm depth (304 × 304 × 512 pixels). In the fast transverse scanning direction, 304 A-scans were sampled. Two repeated Bscans were captured at a fixed position before proceeding to the next location. A total of 304 locations along a 3 mm or 6 mm distance in the slow transverse direction were sampled to form a 3D data cube. The SSADA algorithm split the spectrum into 11 sub-spectra and detected blood flow by calculating the signal amplitude-decorrelation between two consecutive B-scans of the same location. All 608 B-scans in each data cube were acquired in 2.9 seconds. Two volumetric raster scans, including one x-fast scan and one y-fast scan, were obtained and registered [20].

Overview of advanced image processing
Segmentation of OCT angiography 3D flow data allows visualization and analysis of isolated vascular beds. OCT structural images provide reference boundaries for the segmentation of 3D OCT angiograms. Useful reference boundaries ( Fig. 1(A)) include, but are not limited to, the inner limiting membrane (ILM), outer boundary of the inner plexiform layer (IPL), inner nuclear layer (INL), outer boundary of the outer plexiform layer (OPL), outer nuclear layer (ONL), retinal pigment epithelium (RPE), and Bruch's membrane (BM). Vascular layers or "slabs" are identified by two relevant tissue boundaries. For example, retinal circulation is between the boundaries Vitreous/ILM and OPL/ONL. B-scan images can be automatically segmented by a graph search technique [16,21]. We used a conceptually simple directional graph search technique and simplified the complexity of the graph to reduce the computation time. When pathology severely disrupts normal tissue anatomy, manual correction is required. Manual correction of B-scan segmentation was propagated forward and backward across multiple B-scans, expediting image processing and reducing manpower cost. In evaluation, our approach shows good efficiency and accuracy in clinical cases.
We created composite cross-sectional OCT images by combining color-coded angiogram B-scans (flow information) superimposed on gray-scale structural B-scans ( Fig. 1(A)), presenting both blood flow and retinal structure together. This provided detailed information on the depth of the microvasculature network.
OCT angiograms are generated by summarizing the maximum decorrelation within a slab encompassed by relevant anatomic layers [6]. The 3D angiogram slabs are then compressed and presented as 2D en face images so they can be more easily interpreted in a manner similar to traditional angiography techniques. Using the segmentation of the vitreous/ILM, IPL/INL, OPL/ONL, and RPE/BM, five slabs can be visualized as shown in Figs. 1(B)-1(D), 1(G) and 1(H).

Advanced image processing: healthy eye
In a healthy eye, the vitreous is avascular, and there is no flow above the vitreous/ILM boundary. Therefore, the en face image will appear black ( Fig. 1(B)). The superficial inner retinal angiogram (between vitreous/ILM and IPL/INL) shows healthy retinal circulation with a small foveal avascular zone ( Fig. 1(C)). The deep inner retina angiogram (between IPL/INL and OPL/ONL) shows the deep retinal plexus which is a network of fine vessels ( Fig. 1(D)). The inner retina angiogram ( Fig. 1(E)) is a combination of two superficial slabs ( Fig. 1(C) and 1(D)).
Blood flow from larger inner retinal vessels casts a fluctuating shadow, inducing signal variation in deeper layers. This variation is detected as decorrelation and results in a shadowgraphic flow projection artifact. Signal characteristics alone cannot distinguish this shadowgraphic flow projection from true deep-tissue blood flow, but it can be recognized by its vertical shadow in the cross-sectional OCT angiogram (white arrows in Fig. 1(A)). A comparison of Figs. 1(F) and 1(E) reveals projection and replication of the vascular patterns from superficial slabs in the deeper layers. This is particularly evident in the outer retinal slab, where the retinal pigment epithelium (RPE) is the dominant projection surface ( Fig. 1(F)). Subtracting the angiogram of the inner retina from that of the outer retina can remove this artifact, producing an outer retinal angiogram devoid of flow, as would be expected in a healthy retina ( Fig. 1(G)). Flow detected in the outer retinal angiogram after removal of projection artifact is pathologic [13,22]. The choriocapillaris angiogram (RPE/BM to 15 µm below) shows nearly confluent flow ( Fig. 1(H)). Figure 1(I) shows an en face thickness map of the retina, segmented from vitreous/ILM to RPE/BM.
After removal of flow projection, the outer retinal en face angiogram ( Fig. 1(G)) is then used as the reference for removing shadowgraphic flow projection on cross-sectional images. This produces composite B-scan images with color-coded flow corresponding to various slabs, without vertical shadowgraphic artifacts in the outer retina ( Fig. 1(J) compared to Fig.  1(A)). Similarly, we can generate composite C-scan images. Because of the curved nature of the retina, the volume data is flattened using RPE/BM to produce flat C-scan images ( Fig.  1(K)).

Directional graph search
Graph search is a common technique for image segmentation [23][24][25]. We designed a directional graph search technique for retinal layer segmentation. Since retinal layers are primarily horizontal structures on B-scan structural images, we first defined an intensity gradient in depth along the A-line, with each pixel assigned a value G x,z, where and I x,z is the intensity of the pixel, and I x,z-1 is the intensity of the previous pixel located within the A-line. From this, we established a gradient image by normalizing each G x,z value with the function where C (1) is a normalized value between 0 and 1, and min(G) and max(G) are the minimum and maximum G, respectively, for the entire B-scan structural image containing W columns and H rows. An example of a gradient image is displayed in Fig. 2(B) and assigns light-todark intensity transitions as having low C (1) values, shown here with the dark line at the NFL/GCL, IPL/INL, OPL/ONL and RPE/BM tissue boundaries. Because retinal tissue boundaries displayed on structural OCT B-scans ( Fig. 2(A)) have two types of intensity transitions (i.e. light-to-dark and dark-to-light [21]), an inverse gradient image was also generated using the function (2) (1) ,, 1 x z x z CC  thereby defining dark-to-light intensity transitions with a low C (2) value, demonstrated by the horizontal black lines in Fig. 2(C) at vitreous/ILM, IS/OS boundary, and INL/OPL. Graph search segments an image by connecting C values with the lowest overall cost. Typically, a graph search algorithm considers all 8 surrounding neighbors when determining the next optimal connection (Fig. 3(A)). Our directional graph search algorithm considers only 5 directional neighbors, as illustrated by the 5 dashed lines of Fig. 3(B). Because retinal layers are nearly flat, it can be assumed tissue boundaries will extend continuously across structural OCT B-scans, unlikely to reverse in direction. Since we perform the graph search directionally, starting left and extending right, the neighbor to C x,z that is likely to have the lowest connection cost will be on the right side. Therefore, we exclude from the search the left side neighbors and the upward C x,z-1 and downward C x,z+1 neighbors. We can then include C x+1,z-2 and C x+1,z+2 positions to make our directional graph search sensitive to stark boundary changes. We assign a weight of 1 to C x+1,z-1 , C x+1,z , and C x+1,z+1 and a weight of 1.4 to C x+1,z-2 and C x+1,z+2 , thus giving extra cost to curvy paths.  In order to automatically detect the start point of a retinal layer boundary, the directional graph search starts from a virtual start point located outside the graph, such that all adjacent neighbors are in the first column ( Fig. 3(B)). The lowest cost of connecting C values then ends at the rightmost column. This directional graph search method reduces computation complexity since fewer neighbors are considered, and therefore improves segmentation efficiency.
Automated image segmentation of retinal layers using graph search is a common practice with image processing and has been described at length in the literature [16,21,[26][27][28]. Similar to previous demonstrations [18,21], our directional graph search detects seven boundaries of interest one by one on a B-scan image (Fig. 2(A)).

Propagated 2D automated segmentation
In the clinic, OCT images often contain pathological abnormalities such as cysts, exudates, drusen, and/or layer separation. These abnormalities are difficult to account for on conventional 2D and 3D segmentation algorithms [16,28]. Figure 4(A1) shows an example of layer segmentation attracted to strong reflectors, exudates in this case. Our propagated 2D automated segmentation takes into consideration the segmentation result of the previous Bscan, assuming that boundaries do not change much in adjacent B-scans. Specifically, it first segments a B-scan frame with relatively few pathological structures, which is chosen by the user. To segment the remaining B-scans using directional graph search, we further confine the boundary to be within a range that is 15 µm above and below the same boundary in the previous B-scan frame. The segmentation is propagated to the rest of the volume frame by frame. Figure 4(A2) shows the propagated automated segmentation provides accurate segmentation even in tissue disrupted by exudates.
The en face images Fig. 4(B1) and 4(B2) map the distance between the segmented OPL/ONL position and the bottom of the image, with each horizontal line corresponding to a B-scan frame. The conventional 2D algorithm (Fig. 4(B1)) shows segmentation errors, while an accurate segmentation from propagated 2D algorithm generates a continuous map (Fig.  4(B2)). This map facilitates monitoring and identification of possible segmentation errors. Fig. 4. Comparison of performance on pathologic tissue using 2D automated segmentation (A1, B1), and propagated 2D automated segmentation (A2, B2). En face images B1 and B2 maps the position of the OPL/ONL boundary of A1 and A2. Red arrows in A1 and A2 point to the segmentation differences. The colorbar of B1 and B2 is the same as Fig. 1(I).

Propagated 2D automated segmentation with intelligent manual correction
When propagated 2D automated segmentation fails, expert manual correction is required. In the manual correction mode, the user pinpoints several landmark positions with red crosses (Fig. 5(B), 4 red crosses). An optimal path through these landmarks is automatically determined using directional graph search. After manual corrections are made on B-scans within a volume, the corrected boundary curve is propagated to adjacent frames. For example, in Fig. 5, only frame n was manually corrected, and the manual correction successfully propagated to frame n + 30, as shown in Fig. 5 (propagated correction), identified by the red arrow. Fig. 5. Illustration of 2D automated segmentation with and without intelligent manual correction. Manual correction (middle image, red crosses) was performed on frame n, and the correction propagated to frame n + 30. Red arrows identify the segmentation differences.

Semi-automatic segmentation
For cases where the retina is highly deformed and the automated segmentation completely fails, we devised a semi-automatic segmentation method. Similar to intelligent scissors [23], while the users moves the cursor along the boundary path, directional graph search is applied and displayed locally in real time (Fig. 6(A)).

Interpolation mode
Automated/manual segmentation can also be applied at regular intervals ( Fig. 6(B)), followed by interpolation across the entire volume. This greatly reduces the segmentation workload while maintaining reasonable accuracy. The frame interval for manual segmentation is determined according to the variation among B-scan frames, usually 10 to 20 for 3 × 3 mm scans and 5 to 10 for 6 × 6 mm scans. This provides a balance between segmentation accuracy and required manual segmentation workload.

Volume flattening
On 6 × 6 mm images, the automated segmentation using directional graph search may fail due to significant tissue curvature, as we will show in the result section ( Fig. 9(A), inside the yellow box). A flattening procedure was utilized to solve this problem (Fig. 7). We first found the center of mass (pixel intensity) of each A-scan, represented as blue dots in Figs. 7(B1) and 7(C1). We then fitted a polynomial plane to these centers of mass. A shift along the depth (z) was performed to transform the volume, so that the curved plane transformed to a flat plane ( Fig. 7(B2)), flattening the curved volume (compare Figs. 7(A1) and 7(A2)). By using the center of mass instead of an anatomic tissue plane, the volume flattening procedure is not subject to boundary distortion caused by pathology. Note that this volume flattening procedure is only used to aid the segmentation. For visualization, the volume is flattened using the RPE/BM boundary after segmentation.

Standard procedure
As a first step, volume flattening is performed on and only on 6 × 6 mm scans (2.3.6). The user then chooses a frame with few pathological structures, runs 2D automated segmentation (2.3.2), and corrects segmentation errors by providing several key points (2.3.3) or using semi-automatic segmentation (2.3.4). Propagated 2D automated segmentation (2.3.3) then segments the rest of the frames. If the user observes a segmentation error, he or she performs manual correction and then rerun the propagation. In rare cases when propagation fails to correct errors, the user can manually segment selected frames and perform interpolation (2.3.5). It should be noted that often the user only needs to use interpolation to segment one or two of the boundaries, and automatic segmentation is able to work out the rest of the boundaries. Segmentation is always performed under the supervision of the user to minimize errors.

Study population
We systematically tested our segmentation technique in eyes with DR and AMD. In the DR study, 5 normal cases, 10 non-proliferative DR (NPDR), and 10 proliferative DR (PDR) were studied. The layers of interest for segmentation include vitreous/ILM, IPL/INL, OPL/ONL, RPE/BM. In the AMD study, 4 normal, 4 dry AMD and 4 wet AMD eyes were examined. The layers of interest for segmentation included vitreous/ILM, OPL/ONL, IS/OS, RPE/BM. Table 1 summarizes the average number of layers corrected and processing time.

Automated segmentation of pathology
The 2D automated algorithm correctly segmented the tissue boundaries despite disruption caused from RNV (Fig. 8(A)), small exudates ( Fig. 8(B)), small intraretinal cysts (Fig. 8(C)), or drusen with strong boundary (Fig. 8(D)). However, the algorithm failed in some severe pathological cases. During segmentation of epiretinal membrane (Fig. 8(E) 8(H)) caused the segmentation of the IS/OS and RPE to not accurately follow the more elevated drusens, and as a consequence, NFL/GCL, IPL/INL, and INL/OPL were also segmented incorrectly. In these cases, propagated 2D automated segmentation with intelligent manual correction was applied. And in rare cases, interpolation mode was used (e.g. IPL/INL and OPL/ONL in Fig. 11, PDR with edema).

Segmentation processing time and accuracy
During segmentation, we recorded the number of manual corrections made on each type of boundary for both DR and AMD. The average number of manual corrections is given in Table 1. The automated segmentation of vitreous/ILM was highly accurate in DR and AMD cases. In severe DR cases, edema sometimes caused tissue boundaries to be located outside of the searching regions, and therefore required manual corrections of both the IPL/INL and OPL/ONL boundaries. Similarly, AMD with large drusen caused segmentation failure of OPL/ONL and IS/OS. RPE/BM needed to be manually corrected in AMD cases where the boundary became unclear. In general, an increase in severity of either DR or AMD required a longer average processing time. Compared to a purely manual segmentation approach (typically taking 3-4h to complete [19],), our intelligent manual correction method efficiently segmented tissue boundaries in eyes with DR and AMD, only taking 15 minutes to complete, including the most difficult case. To evaluate segmentation accuracy, we compared the results of manual segmentation (using intelligent scissors) with those from our propagated automated segmentation with manual correction. For each case, 2 subjects were randomly chosen and 20 B-scans were randomly selected for evaluation. Three graders independently performed manual segmentation of each tissue boundary, with the help of intelligent scissors. The manually segmented boundaries were averaged among the three graders and taken as gold standard. The absolute errors of our propagated automated segmentation with manual correction was determined (mean ± std in unit of pixels). The result is given in Table 2. In more than 62% of images, the segmentation error is less than 1 pixels (3.1 µm).

Volume flattening
The aforementioned volume flattening procedure was able to solve stark curvature segmentation errors. The yellow box in Fig. 9(A) demonstrates segmentation failure at multiple tissue boundaries in an area of stark curvature. By flattening the volumetric data, our automated segmentation algorithm was able to accurately segment all seven tissue boundaries inside the yellow box as shown in Fig. 9(B). When the image was restored to its original curvature, the corrected segmentation remained ( Fig. 9(C)). This automated volume flattening allows for efficient image processing of large area OCT scans (e.g. 6 × 6 mm).

Age related macular degeneration (3 × 3 mm scans)
CNV, which is the pathologic feature of wet AMD, occurs when abnormal vessels grow from the choriocapillaris and penetrate Bruch's membrane into the outer retinal space [13]. Detection of CNV depends on the segmentation of three reference planes (Vitreous/ILM, OPL/ONL, and RPE/BM) used to generate three slabs: inner retina, outer retina, and choriocapillaris. Figure 10 shows representative images of dry and wet AMD cases. Structural information from the OCT angiography scans was used to create retinal thickness (Figs. 10(B1) and 10(B2)) and RPE-drusen complex (RPEDC) maps (Figs. 10(D1) and 10(D2), distance between IS/OS and RPE/BM). The retinal thickness map is clinically useful in determining atrophy and exudation. The RPEDC map, representing the size and the volume of drusen, has been correlated with risk of clinical progression [29].
Because of shadowgraphic flow projection, true CNV is difficult to identify in both the composite B-scan and en face angiogram. We used an previously published automated CNV detection algorithm to removed projection from the outer retinal [22]. A composite en face angiogram displaying the two retinal slabs in different colors shows the CNV in relation to the retinal angiogram (Figs. 10(C2) and 12(A1)). An overlay of this composite angiogram on the cross-sectional angiogram shows the depth of the CNV in relation to retinal structures ( Fig. 10(A2)). The size of CNV can be quantified by calculating the area of the vessel in the outer retinal slab.

Diabetic retinopathy (3 × 3 mm scans)
RNV, or growth of new vessels above the ILM, is the hallmark of proliferative diabetic retinopathy (PDR). The presence of RNV is associated with high risk of vision loss and is an indication for treatment with panretinal photocoagulation, which reduces the risk of vision loss [30].
Segmenting along the vitreous/ILM border reveals the RNV in the vitreous slab, distinguishing it from intra-retinal microvascular abnormalities (IRMA), which can be difficult to distinguish clinically from early RNV (Fig. 11, case 3, PDR without edema and Fig. 12(A2)). By quantifying the RNV area, one can assess the extent and activity of PDR. Figure 11 shows representative images of DR. The first row shows color-coded B-scans without flow projection artifacts and the boundaries for en face projection. Presenting the structural and flow information simultaneously clarifies the anatomic relationship between vessels and tissue planes. En face composite angiograms of the superficial and deep plexus (second and third rows, respectively) disclose vascular abnormalities including RNV, IRMA, thickening/narrowing of vessels, and capillary dropout as with typical dye-based angiography.
Capillary nonperfusion is a major feature of DR that is associated with vision loss and progression of disease [31,32]. Using an automated algorithm, we identified and quantified capillary nonperfusion [4,13] and created a nonperfusion map (Fig. 11(D)) showing blue areas with flow signal lower than 1.2 standard deviations above the mean decorrelation signal in the foveal avascular zone. Assessing the distance from vitreous/ILM to RPE/BM across the volume scan created the retinal thickness map (Fig. 11(E)). This allows the clinician to assess the central macula for edema, atrophy, and distortion of contour.

Clinical evaluation of 6 × 6 mm scans
The pathology in AMD and DR can extend beyond the central macular area. While OCT angiography cannot match the field of view of the current dye-based widefield techniques [33, 34], 6 × 6 mm OCT angiography scans cover a wider area and can reveal pathology not shown in 3 × 3 mm scans. Figures 12(A1) and 12(A2) show examples of 6 × 6 mm scans of the wet AMD case in Fig. 10 and PDR without edema case seen in Fig. 11, respectively. Although these scans are of lower resolution, 6 × 6 mm scans captured areas of capillary nonperfusion not present in the 3 × 3 mm scan area (black areas outside of the blue square). Fig. 12. A1, B1, results of 6 × 6 mm scan for the wet AMD case in Fig. 10. A2, B2, results of 6 × 6 mm scan for the PDR without edema case in Fig. 11. The blue square marks the 3 × 3 mm range corresponding to Fig. 10 and Fig. 11.

Conclusion
We have described in detail advanced image processing for OCT angiography quantification and visualization. Our proposed segmentation method shows good accuracy and efficiency in clinical applications. In the current phase of development, segmentation still requires manual correction in a minority of cases, but its frequency and associated workload has been highly reduced with techniques such as semi-automatic segmentation, propagated manual corrections, and interpolation as compared to previous reports [18,19]. We also showed innovative ways of visualizing OCT angiography data including composite B-scan images and composite en face angiograms. Integration of these methods into commercial OCT angiography instruments can potentially improve the utility and diagnostic accuracy of OCT angiography.