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

Multiscale joint segmentation method for retinal optical coherence tomography images using a bidirectional wave algorithm and improved graph theory

Open Access Open Access

Abstract

Morphology and functional metrics of retinal layers are important biomarkers for many human ophthalmic diseases. Automatic and accurate segmentation of retinal layers is crucial for disease diagnosis and research. To improve the performance of retinal layer segmentation, a multiscale joint segmentation framework for retinal optical coherence tomography (OCT) images based on bidirectional wave algorithm and improved graph theory is proposed. In this framework, the bidirectional wave algorithm was used to segment edge information in multiscale images, and the improved graph theory was used to modify edge information globally, to realize automatic and accurate segmentation of eight retinal layer boundaries. This framework was tested on two public datasets and two OCT imaging systems. The test results show that, compared with other state-of-the-art methods, this framework does not need data pre-training and parameter pre-adjustment on different datasets, and can achieve sub-pixel retinal layer segmentation on a low-configuration computer.

© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Optical coherence tomography is a non-invasive imaging technology with high resolution, which has been widely used in ophthalmology [1]. The retinal layer structure is an important biomarker for diseases such as Cystoid Macular Edema (CME), Alzheimer’s disease (AD), Multiple Sclerosis (MS), Diabetic Macular Edema (DME), and Age-related Macular Degeneration (AMD) [24]. Quantification of retinal layer structure from OCT images is necessary for reliable disease detection, clinical diagnosis, and efficacy monitoring [57]. The retina is usually divided into seven adjacent layers [8] (eight-layer boundaries), as shown in Fig. 1. The seven layers from top to bottom are: the Nerve Fiber layer (i, NFL), the Ganglion Cell and Inner Plexiform complex (ii, GCL-IPL), the Inner Nuclear layer (iii, INL), the Outer Plexiform layer (iv, OPL), the Outer Nuclear layer and Inner Segment complex (v, ONL-IS), the Outer Segment (vi, OS), and the Retinal Pigment Epithelium complex (vii, RPE) [8]. The identification of the seven layers is accomplished by segmenting eight boundaries: Internal Limiting Membrane (I, ILM), NFL-GCI (II), GCI-INL (III), INL-OPL (IV), OPL-ONL (V), IS-OS (VI), Inner Retinal Pigment Epithelium (VII, IRPE) and outer aspect of the Bruch Membrane (VIII, OBM). The thickness of each layer of the retina is different (20µm ∼ 80µm), and the total retinal thickness is generally between 150µm and 300µm.To distinguish each layer clearly, the resolution of OCT equipment is usually less than 10µm. Segmentation of the diseased retina is a challenge because retinopathy can change the normal layer structure of the retina. For example, MS can lead to thinning of the thickness of the retinal layer, AMD can change the structure of the RPE layer, and DME can destroy the layer structure from GCL-IPL to ONL-IS.

 figure: Fig. 1.

Fig. 1. Retinal layers of OCT image centered at the macula. The boundaries from top to bottom are ILM (I), NFL-GCI (II), GCI-INL (III), INL-OPL (IV), OPL-ONL (V), IS-OS (VI), IRPE (VII) and OBM (VIII).

Download Full Size | PDF

Automatic and accurate segmentation of retinal layers is important for ophthalmic disease diagnosis, and many methods have been proposed: Graph [9,10], Active contour [11,12], and Machine learning [1317]. The OCT segmentation method based on active contour was first proposed by Fernández et al. [18]. Yazdanpanah et al. [12] improved Fernández's method and realized the segmentation of intra-retinal layers. In [19], an approximate parallel constraint is integrated into the active contour, which improves the applicability of the method. However, a good initial estimation of the layer boundary is necessary for these methods, otherwise active contour methods may be entrapped in a local energy minimum. Compared with active contour methods, graph methods can find the global optimal solution. In [20], two extensions were added to the graph segmentation method: incorporating varying feasibility constraints and incorporating true regional information, which increased the utility of the graph segmentation method. In [10], a 2D graph theory method based on the shortest path algorithm is proposed, which uses Dijkstra's algorithm to detect the shortest path in an undirected graph. Graph methods can obtain the minimum global error, but graph methods are sensitive to noise and image degradation. Active contour methods and graph theory methods do not require training data, but both require complex constraints and parameter settings, which limit the applicability of the above methods.

Many studies have used machine learning methods to detect retinal layer boundaries. In [21], a random forest classifier is used to estimate boundary probabilities, and these probabilities are used by graph methods to find the optimal segmentation result. In [22], a layer segmentation method combining convolutional neural network (CNN) and generalized U-net network architecture is presented, which is robust to the presence of retinal pathology. In [15], a fully convolutional deep architecture (ReLayNet) is proposed, which is an end-to-end learning method to segment retina into seven layers. In [16], a framework (FCRN) for structured layer surface segmentation was proposed, which can obtain continuous retinal layer surface in an end-to-end deep learning model. Machine learning methods do not have complex constraints, as long as there is enough training data, a satisfactory segmentation result can be obtained. However, the outcome of this approach is dependent on the quantity and quality of labeled data given the difficulty in collecting high-quality data from large and diverse patient populations and corresponding data labeling [14,17].

There are three main challenges in automatic retinal layer segmentation. First, the quality of OCT image is normally degraded by speckle noise [23], which makes layer boundaries unclear and layer segmentation prone to error. Secondly, Pathological alterations [8] in retinal tissue morphology and vessel shadows increase the probability of segmentation failure. Thirdly, there are some differences between OCT images acquired by different instruments [24]. Without manual parameter setting or adding training samples, it is difficult to apply the same segmentation method to different OCT images acquired by different instruments.

To solve the above problems, a joint segmentation framework combining the local optimal solution and the global optimal solution is proposed. This method is insensitive to speckle noise, and can adaptively and accurately segment 8-layer boundaries from retinal images without data training and manual parameter adjustment. This method is convenient to use on OCT devices and has lower performance requirements on computers. This method has the potential to be widely applied to retinopathy images and retinal images acquired by different imaging systems. The main contributions of this paper can be concluded as follows:

  • 1. The bidirectional wave algorithm with two-way boundary detection is proposed, which improves the layer structure recognition ability and layer detection speed.
  • 2. We improved the shortest path algorithm and designed a specially constructed undirected graph with dynamic weights, which can greatly reduce the amount of data and effectively reduce the impact of speckle noise on layer segmentation.
  • 3. We propose a multiscale joint framework based on bidirectional wave algorithm and improved graph theory, which shows better robustness and adaptability in retinal layer segmentation. Compared with other methods, the method proposed can achieve the same or better retinal segmentation results without data pre-training and parameter pre-adjustment.

2. Method

The multiscale joint segmentation framework in this paper consists of three parts: preprocessing, bidirectional wave algorithm, and improved graph theory, as shown in Fig. 2. The main function of preprocessing are speckle noise suppression, variance calculation and image multiscale transformation, which is the necessary basis for image segmentation. This segmentation framework uses the bidirectional wave algorithm to detect retinal edge information of multiscale images and then uses the improved graph theory to modify the edge information globally, so that the accurate 8-layer boundaries can be segmented.

 figure: Fig. 2.

Fig. 2. A schematic of the proposed method. Gray represents the preprocessing part, blue represents the bidirectional wave algorithm part, and green represents the improved graph theory part.

Download Full Size | PDF

2.1 Preprocessing

Figure 3 shows the preprocessing process. First, the original OCT image is enhanced by background noise removal (If the pixel gray value is less than the overall average gray value of the image, the pixel gray value is set to zero) and contrast stretching (If the pixel gray value is greater than the overall average gray value of the image, the pixel gray value should be linearly stretched to 0-255). The enhanced images (image size is W × H) were then downsampled by 2 and 4 to obtain two new images of sizes W/2×H/2 and W/4×H/4, respectively. Multiscale transform generates images that show structural information of retina at different resolution levels, which is helpful to improve the robustness of layer segmentation. Multiscale images are simultaneously filtered by Gaussian low-pass filter (the size of Gaussian kernel is W × H, and the filtering threshold is W/10) to smooth the image, which helps segment a more complete layer edge. Finally, the variance of the smoothed image is calculated to obtain the target area (the kernel size for calculating the variance is 3 × 3), which largely reduced the calculation workload.

 figure: Fig. 3.

Fig. 3. Diagram of the preprocessing part. Red represents (W/4×H/4) size images, orange represents (W/2×H/2) size images and pink represents the original size image. W is the width of the image and H is the height of the image.

Download Full Size | PDF

2.2 Bidirectional wave algorithm

The wave algorithm [25] is a new edge detection method proposed by our laboratory. Wave algorithm regards the bright and dark changes in OCT image as the change of potential energy in fluid, thus transforming the edge detection problem into a potential energy calculation problem. Although the wave algorithm can detect layer boundaries well, it can only detect the boundary that changes from dark to bright (gray value from low to high). To detect the boundary from bright to dark (gray value from high to low), the image needs to be rotated 90 degrees [26]. The wave algorithm is composed of the potential energy correction equation (PEC) and the wave potential energy equation (WPE). The basic reason why the wave algorithm can only detect the boundary unidirectionally is that the wave potential energy equation lacks the ability to detect the direction of potential energy change. The original wave potential energy equation [25] is:

$$WPE\textrm{ = }v \times {v_q} \times \sigma + {\varphi _g}$$
where φg is the gravitational potential energy, v is the fluid velocity in the Z template, σ is the scale factor, vqis the fluid velocity in the Q template. To enable the wave algorithm to bidirectionally detect the boundary, the wave potential energy equation is improved. In this paper, the Q template and the Z template are merged into the M template, v and vq are redefined as V. V is the normalized difference value in the M template, and the mathematical expression of V is shown in Fig. 4. Even gray change direction is different from the M template movement direction, as long as the gray change trend (gray changes from low to high or high to low) in the M template is more obvious, the value of the kinetic energy V2 will be larger.

V2 can better represent the value of the kinetic energy in the template, but it cannot determine the direction of the kinetic energy. To enable the algorithm to automatically identify the kinetic energy direction, a kinetic energy direction equation is presented in this paper:

$$\left\{ {\begin{array}{{c}} {s = \sum\limits_{i \in M} {\sum\limits_{j \in M} {\frac{{I(i,j) - I(i,j - 1)}}{{|{I(i,j) - I(i,j - 1)} |}}} } }\\ {\tanh (s) = \frac{{{e^s} - {e^{ - s}}}}{{{e^s} + {e^{ - s}}}}} \end{array}} \right.$$
where s represents the number of pixels with the same gray change trend in the template. s is positively related to the gray change from low to high and negatively related to the gray change from high to low. V is related to the value of the kinetic energy, and s is related to the direction of the kinetic energy. The kinetic energy direction factor tanh(s) can be obtained by substituting s as a variable into the hyperbolic tangent function. According to the properties of the hyperbolic tangent function, tanh(s) ranges from -1 to 1 as s varies from -12 to 12. The improved wave potential energy equation is:
$$BWPE\textrm{ = }\tanh (s) \times ({V^2} \times \sigma + {\varphi _g})$$

 figure: Fig. 4.

Fig. 4. Improvement of velocity definition. (i, j) is the pixel to be detected in the image, and I(i, j) is the gray value of the pixel in the image. The 3 × 3 window centered on (i, j) is the Z template. The 3 × 2 window overlap with the left end of the Z template is the Q template. The 3 × 5 window centered on (i, j) is the M template. All the templates in this paper are moved from left to right (direction of movement). The normalized difference value in the Z template is defined as v, and the normalized difference value in the Q template is defined as vq.

Download Full Size | PDF

It is obvious from the above formula that the kinetic energy direction equation adds the sign to the unsigned wave potential energy equation. When the template kinetic energy direction is the same as the template movement direction, tanh(s)→1, that is, BWPE is a positive number. When the template kinetic energy direction is opposite to the template movement direction, tanh(s)→-1, that is, BWPE is a negative number. Finally, the boundary line can be obtained through the improved potential energy correction equation:

$$BPEC\textrm{ = }\left\{ \begin{array}{@{}cc@{}} + 1,& ({(K > 1) \cap (C > 1)} )\cap ({BWPE \ge \alpha } )\\ - 1,&({(K < 1) \cap (C < 1)} )\cap ({BWPE \le ( - \alpha )} )\\ 0,& others \end{array} \right.$$
where K is the geometric discriminant equation (K = Kb/Kc, Kb is the slope of a straight line composed of three points behind the measuring point, and Kc is the slope of a straight line composed of three points before the measuring point), C is the gray discriminant function (C = Sc/Sb, Sc is the average gray value of three points before the measuring point, and Sb is the average gray value of three points behind the measuring point), and α is the potential energy threshold (α=meanB/meanT, meanB is the average gray value of the background, and meanT is the average gray value of the target). When BPEC=+1, it means that the detected point is the positive edge. When BPEC = -1, it means that the detected point is the negative edge. When BPEC = 0, it means that the detected point is not an edge. This design not only enables the algorithm to detect the boundary in both directions but also improves the detection speed (edge detection part) by about 50%.

In this paper, the three sizes of images are processed using the improved wave potential energy equation, and three different maps are obtained as shown in Fig. 5. Among them, the yellow area represents the area where BWPE is positive (positive area), and the cyan area represents the area where BWPE is negative (negative area). The positive edge (yellow line) can be detected from the positive area by the improved potential energy correction equation, and the negative edge (cyan line) can also be detected from the negative area. The potential energy map (Fig. 5(a).) derived from the original image has rich detail, from which all texture information of retina can be segmented, but it may also contain false edges caused by noise. The potential energy map (Fig. 5(b).) derived from the W/2×H/2 size image has a large Receptive Field, from which the structural framework of retina can be segmented. The potential energy map (Fig. 5(c).) derived from the W/4×H/4 size image only has the general outline of retina, from which useful edge information is difficult to segment, but it can be used as a segmentation guide map. The segmentation guide map can divide the retinal layer structure into four areas (Fig. 5(f).): the first area (Green) contains ILM and NFL-GCI, the second area (Blue) contains GCI-INL, INL-OPL, and OPL-ONL, the third area (Pink) contain IS-OS and IRPE, and the fourth area (Red) contains OBM.

 figure: Fig. 5.

Fig. 5. Schematic diagram of improved wave algorithm. (a∼c) Potential energy map, (d∼e) Boundary Map, (f) Guide map. The improved wave potential energy equation can effectively detect two types of areas, namely the area where the gray value changes from dark to light (yellow) and the area where the gray value changes from light to dark (cyan). The improved potential energy correction equation can detect two types of boundary lines from two types of areas.

Download Full Size | PDF

2.3 Graph theory with dynamic weights

The bidirectional wave algorithm is an edge detection algorithm based on local information, and it is difficult to acquire the global optimal solution. To acquire the global best segmentation effect, the Dijkstra algorithm [27] is used to modify the segmentation results of the bidirectional wave algorithm. To deeply fuse the Dijkstra algorithm with the bidirectional wave algorithm, a specially constructed undirected graph with dynamic weights is designed. In this paper, each image is represented as a network (shown in Fig. 6.), where each node corresponds to a pixel. However, not all pixels are nodes, only the pixel corresponding to the layer boundary (yellow and green lines in Fig. 5(d). and Fig. 5(e).) detected by the bidirectional wave algorithm are nodes. A link connecting two nodes is called an edge. The weight of the edge is the weighted sum of all pixels on the edge:

$$W\textrm{ = }\sum\limits_{i \in L} {{w_i}} $$
where L is the edge, W is the weight of the edges, i is the pixel that L passes through, and wi is the weight of the pixel. A set of connected edges form a pathway, and the layer boundary is the pathway with the minimum weight sum. The network defined in this paper has many advantages. First, nodes are no longer adjacent pixels, but edge pixels selected by the bidirectional wave algorithm, which can greatly reduce the amount of data and speed up the segmentation. Secondly, the weighted sum of pixels is the weight of edges, which can reduce the negative effect of speckle noise on layer segmentation.

 figure: Fig. 6.

Fig. 6. Network definitions in improved graph theory. Green represents nodes (BPEC≠0), Yellow represents the pathway (edges), Gray represents the background area, and Blue represents the gray value change areas (target area) detected by the bidirectional wave algorithm. “+” represents positive areas and positive edges, and “-” represents negative areas and negative edges.

Download Full Size | PDF

Edge weights are important to segmentation accuracy, and common metrics for the weight value include pixel-to-pixel distances or differences between intensity values [12]. However, the edge weights in this paper are the weighted sum of all pixels corresponding to the edge, and the pixel weight is set according to the bidirectional wave algorithm. In the weight map, only the pixels in the potential energy area (positive area and negative area) detected by the bidirectional wave algorithm are given weights (0∼1), and the weights of pixels in other areas are set to the maximum value (255). In this paper, the weights of pixels are not constant, and the weights of pixels will change dynamically according to BPEC. If the boundary to be detected is a positive boundary (ILM, INL-OPL, IS-OS, and IRPE), the pixel weight of positive areas and positive edges is set to be lower (0∼1), and that of negative areas and negative edges to be higher (255). If the boundary to be detected is a negative boundary (NFL-GCI, GCI-INL, OPL-ONL, and OBM), the pixel weight of positive areas and positive edges is set to be higher (255), and that of negative areas and negative edges to be lower (0∼1).

2.4 Automatic segmentation of eight-layer boundaries

The nodes map of the retina (Fig. 7(a).) can be obtained by applying the OR operation to the W/2×H/2 size boundary image and the W × H size boundary image. Nodes map consists of positive boundaries and negative boundaries. The frame map of the retina (Fig. 7(b).) can be obtained by applying the OR operation to the W/2×H/2 size potential energy image and the W × H size potential energy image. Frame map consists of the positive area and the negative area, the positive area contains four positive boundaries, and the negative area contains four negative boundaries. The W/4×H/4 size potential energy image (Fig. 7(c).) can be used as a guide map to divide the retina into four areas. Before obtaining the Nodes map, Frame map, and Guide map, the images of three sizes are uniformly transformed to the W/2×H/2 size. Combining the frame map and the guide map can divide the retinal image into 8 sub-region, and each sub-region contains a retinal layer boundary. In the 8 sub-region, the nodes map and the frame map can assign weights to the pixels in each sub-region, and get the weight map (the weight value of the nodes map is 0, the weight value of the frame map is 1, and the weight value of other areas is 255) of eight boundaries respectively. Finally, eight complete retinal layer boundary lines (Fig. 7(e).) can be obtained by using the improved Dijkstra algorithm to calculate the optimal path of the eight weight maps (Fig. 7(d).) respectively.

 figure: Fig. 7.

Fig. 7. Automatic segmentation of the 8-layer boundary of the retina. (a) Nodes map, (b) Frame map, (c) Guide map. Before obtaining the Nodes map, Frame map and Guide map, the images of three sizes are uniformly transformed to the W/2×H/2 size.

Download Full Size | PDF

3. Experiments and results

3.1 Data

The segmentation framework is evaluated on two public datasets that contain healthy controls (HC) as well as patients with AMD and MS (MS causes mild thinning of retinal layers). Reference [28] contains 14 healthy controls and 21 people with MS. The data acquired with a Heidelberg Spectralis scanner come with 9-layer boundaries manually delineated in each B-scan. The lateral and axial resolutions are 5.8 µm and 3.87 µm, respectively. There are 49 B-scans from each of the 35 subjects and the size of each B-scan is 496 × 1024. Reference [29] consists of 384 macular scans (115 normal eyes and 269 eyes with AMD) and 3-layer boundaries of each B-scan are manually delineated. The dataset was acquired using Bioptigen imaging systems at 4 different clinical sites. The lateral and axial resolutions are 6.54µm and 3.23µm (Each B-scan had a size of 512 × 1000), respectively.

3.2 Segmentation results

The higher the image quality, the more layers the algorithm can segment. To ensure the universality of the algorithm, the proposed method only segments the seven layers (eight-layer boundaries) that are most widely used in clinical medicine. The difference between the segmentation algorithm and the manual delineation is evaluated using mean absolute distance (MAD) and standard deviation (Std). In [28], eight retinal layer boundaries of ILM (I), NFL-GCI (II), GCI-INL (III), INL-OPL (IV), OPL-ONL (V), IS-OS (VI), IRPE(VII) and OBM (VIII) were segmented (Table 1 and Fig. 8. are the test results). In [29], since experts only manually labeled ILM, IRPE, and OBM, the proposed method only segmented ILM, IRPE, and OBM (Table 1 and Fig. 8. are the test results). Detailed experimental data has been published on GitHub [30]. The proposed method has low requirements for computers, and the implementation environment of our method in this paper is: a personal computer without GPU, CPU is Intel i5-7500, Random Access Memory is 16 G, and the programming software is Visual Studio Community 2017. Many high-performance OCT systems are equipped with GPU, but this increases hardware costs and system complexity. In many underdeveloped areas or community hospitals, most of them use low-cost OCT equipment without GPU. The proposed method can be easily applied to low-cost OCT equipment without adding new hardware, which is beneficial to improving medical care.

 figure: Fig. 8.

Fig. 8. Automatic segmentation results of the proposed method in [28] and [29]. The proposed method (red) and expert manual tracings (green).

Download Full Size | PDF

Tables Icon

Table 1. The mean absolute distance (MAD(Std) in pixel) in [28] and [29]. Since only 3-layer in [29] are manually delineated (ILM, IRPE, and OBM), the MAD (Std) of the unlabeled layers in the table is “–”

Each B-scan in [28] was averaged over at least 12 images at the same location [28], while the B-scan in [29] was not averaged, so it has lower speckle noise and better SNR than [29]. From Table 1, we can see that the overall segmentation error of the proposed method in [28] and [29] are within 1 pixel and 1.5 pixels, respectively. From this, it is not difficult to find that the segmentation accuracy of this framework will improve with the improvement of image quality. For both datasets, the accuracy of retinal layer segmentation in patients (AMD and MS) is lower than that in healthy controls. This is because the presence of pathologies changes the normal physiological structures of the retina, making the segmentation method more error-prone. Overall, without additional training data and manual parameter adjustment, the proposed method is feasible for both high-quality images and low-quality images from patients and normal subjects. From Fig. 8, we can see that the auto-segmentation result (red line) of the proposed method overlaps well with the expert's manual delineation result (green line), which shows that our results are basically the same as manual segmentation.

3.3 Comparison of results

In [28], the proposed method is compared with FCRN [16] and ReLayNet [15]. FCRN and ReLayNet are one of the state-of-the-art deep learning methods. Since FCRN and ReLayNet require training before they can be used, the last 6 HC scans and the last 9 MS scans were used for training, and the other 20 scans are used for testing [13,16]. The proposed method can be used directly without training and manual parameter settings. For direct comparison, the proposed method also only uses the other 20 scans for testing. Table 2 shows the segmentation errors of different methods in [28]. From Table 2, we can see that FCRN, ReLayNet and the proposed method all have a sub-pixel error on average. Machine learning methods trained on a specific dataset generally achieve better segmentation results in that dataset, but such methods are more prone to be affected by limited training data, inaccurate labels, or outlier images.

Tables Icon

Table 2. The mean absolute distance (MAD) and root mean square error (RMSE) in [28] for three of the tested methods (in pixel)

In [29], our method is compared with the 3D optimal surfaces segmentation method (OSCS) [31] and the convolution-based neural network method (CNN-S) [32]. Since CNN-S requires training before it can be used and OSCS needs to set different parameters to achieve the best performance, 330 volumes (90 HC, 240 AMD) are randomly selected in [29] to train the network and optimize the parameters, and 50 volumes (25 HC, 25 AMD) are used for testing [32]. For direct comparison, the proposed method also only uses 50 volumes (25 HC, 25 AMD) for testing. Table 3 shows the segmentation errors of different methods in [29]. It can be seen from Table 3 that the overall segmentation error of three methods in [29] is less than 2 pixels. Compared with other methods, our method has better segmentation accuracy for both normal subjects and AMD patients. Our method does not have data pre-training and parameter pre-tuning, so the segmentation effect of our method may not be the best on some datasets (such as [28]). However, our method can be easily applied to existing commercial equipment (no data pre-training, no parameter pre-adjustment, no GPU), which is very meaningful for the intelligence of medical instruments.

Tables Icon

Table 3. The mean absolute distance (MAD(Std) in pixel) in [29] for three of the tested methods

4. Discussion

CNN-S, FCRN, and ReLayNet are deep learning methods, which require training before they can be used. Due to the diversity of organisms, if deep learning methods are used to segment biological tissues, there needs to be enough sample data for training. In fact, it is difficult to acquire and manually delineate large volumes of diverse medical data. The limitation of deep learning methods to segment diseased retina is that the test image must have corresponding training samples in the training dataset. If retinopathy is not in the training dataset, the probability of segmentation errors is higher. OSCS does not require training but requires manual setting of different parameters for different data to achieve the best performance. The limitation of OSCS to segment diseased retina is that layer segmentation for retinopathy can only be achieved by manually adjusting parameters based on prior knowledge of this retinopathy. The proposed method is an adaptive retinal layer segmentation method that can be used directly without training and manual parameter settings (except for the overall deviation of the boundary line caused by different manual delineations). In Ref. [29], and other OCT images, the proposed method will adaptively adjust parameters (without manual intervention) to achieve the best layer segmentation effect. However, the proposed method also has a limitation, that is, it cannot effectively segment images with severely damaged retinal layer structures. The bidirectional wave algorithm in this paper is directional, so the segmentation effect of the diseased retina with a steep boundary is poor. The DME dataset [33] used the standard Spectralis to collect retinal images of patients with DME [33]. The lateral and axial resolutions are 11.07–11.59µm and 3.87 µm (Each B-scan had a size of 496 × 768), respectively. In [33], the proposed method can segment the retinal layer structure of mild and intermediate DME cases (Fig. 9(c). and Fig. 9(b)., the layer boundary changes gently), but cannot segment the layer structure of advanced DME cases (Fig. 9(a)., retinopathy makes the layer boundary changes steep). If the layer structure of a severely damaged area needs to be segmented, some parameters in the proposed method have to be adjusted manually. In the future, we will try to improve the proposed method by using the thermal diffusion equation so that our method can also effectively segment the diseased retina with steep boundaries. To encourage research in this field, we have published the code of this algorithm [30] for reference.

 figure: Fig. 9.

Fig. 9. Automatic segmentation results of the proposed method in [33]. (a) Advanced DME cases, (b) Intermediate DME cases, (c) Mild DME cases. The proposed method (red) and expert manual tracings (green).

Download Full Size | PDF

5. Additional experiments

To further verify the universality and robustness of our method, we tested the algorithm on two sets of OCT imaging systems developed by Nanyang Technological University with a center wavelength of 1060 nm. One is a point-scan imaging system with an axial resolution of 5µm, and the imaging results are shown in Fig. 10(a). Another is a line-scan imaging system with an axial resolution of 15µm, the imaging result is shown in Fig. 10(b). We use the proposed method to segment the retinal images obtained by two systems respectively (without training data and manual parameter adjustment), and Fig. 10 shows the segmentation results. Even though the image quality of the line-scan image is low, the proposed method can still segment the 8-layer boundaries of retina. From Fig. 10, we can see that our framework can adaptively obtain satisfactory retinal layer segmentation results for images with different resolutions or different imaging modalities. Moreover, the proposed method can be easily applied to these OCT imaging systems and does not require additional GPU configuration or pre-collected data for training. This method has good robustness to different OCT images, and the main reason is that this method skillfully fuses prior knowledge, local information, and global information.

 figure: Fig. 10.

Fig. 10. Tests of the proposed method on different OCT imaging systems. (a) point scan images (size of B-scan is 400 × 400), (b) line-scan images (size of B-scan is 400 × 510). The red line is the segmentation result of the proposed method.

Download Full Size | PDF

6. Conclusion

To adaptively segment retinal OCT images obtained from different OCT instruments without data pre-training and parameter pre-adjustment, a multiscale joint segmentation framework based on the bidirectional wave algorithm and improved graph theory is proposed, which combines the local optimal solution and the global optimal solution. Our framework is evaluated on two public datasets and two OCT imaging systems. Compared with the other four methods, our framework can adaptively segment retinal images acquired by different OCT instruments without training and manual parameter adjustment, which will advance 3D-OCT, OCTA, and other fields that require retinal layer segmentation.

Funding

National Key Research and Development Program of China (2017YFC0109901); Ministry of Education - Singapore (MOE-T2EP30120-0001); Ministry of Health -Singapore (MOH-OFIRG19may-0009).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Data Availability

Data is available in [28,29,33]. Other data may be obtained from the authors upon reasonable request. Codes underlying the results presented in this paper are available in Code 1 [30].

References

1. J. F. D. Boer, R. Leitgeb, and M. Wojtkowski, “Twenty-five years of optical coherence tomography: the paradigm shift in sensitivity and speed provided by Fourier domain OCT,” Biomed. Opt. Express 8(7), 3248–3280 (2017). [CrossRef]  

2. V. T. T. Chan, Z. Sun, S. Tang, L. J. Chen, A. Wong, C. C. Tham, T. Y. Wong, C. Chen, M. K. Ikram, H. E. Whitson, E. M. Lad, V. C. T. Mok, and C. Y. Cheung, “Spectral-domain OCT measurements in Alzheimer’s disease: a systematic review and meta-analysis,” Ophthalmology 126(4), 497–510 (2019). [CrossRef]  

3. A. Rothman, O. Murphy, K. Fitzgerald, J. Button, E. GordonLipkin, J. Ratchford, S. Newsome, E. Mowry, E. Sotirchos, S. SycMazurek, J. Nguyen, N. GonzalezCaldito, L. Balcer, E. Frohman, T. Frohman, D. Reich, C. Crainiceanu, S. Saidha, and P. Calabresi, “Retinal measurements predict 10 year disability in multiple sclerosis Ann,” Ann. Clin. Transl. Neurol. 6(2), 222–232 (2019). [CrossRef]  

4. L. Yihao, A. Carass, H. Yufan, B. J. Antony, A. Filippatou, S. Saidha, S. D. Solomon, P. A. Calabresi, and J. L. Prince, “Layer boundary evolution method for macular OCT layer segmentation,” Biomed. Opt. Express 10(3), 1064–1080 (2019). [CrossRef]  

5. S. A. Read, D. Alonso-Caneiro, and S. J. Vincent, “Longitudinal changes in macular retinal layer thickness in pediatric populations: Myopic vs non-myopic eyes,” PLoS One 12(6), e0180462 (2017). [CrossRef]  

6. W. J. Lee, Y. K. Kim, Y. W. Kim, J. W. Jeoung, S. H. Kim, J. W. Heo, H. G. Yu, and K. H. Park, “Rate of macular ganglion cell-inner plexiform layer thinning in glaucomatous eyes with vascular endothelial growth factor inhibition,” Journal of glaucoma 26(11), 980–986 (2017). [CrossRef]  

7. J. Hamwood, D. Alonso-Caneiro, S. A. Read, S. J. Vincent, and M. J. Collins, “Effect of patch size and network architecture on a convolutional neural network approach for automatic segmentation of OCT retinal layers,” Biomed. Opt. Express 9(7), 3049–3066 (2018). [CrossRef]  

8. A. Chakravarty and J. Sivaswamy, “A Supervised Joint Multi-layer Segmentation Framework for Retinal Optical Coherence Tomography Images using Conditional Random Field,” Comput. Methods Programs Biomed. 165, 235–250 (2018). [CrossRef]  

9. J. Tian, B. Varga, G. M. Somfai, W. H. Lee, W. E. Smiddy, and D. C. DeBuc, “Real-time automatic segmentation of optical coherence tomography volume data of the macular region,” PLoS One 10(8), e0133908 (2015). [CrossRef]  

10. S. J. Chiu, X. T. Li, P. Nicholas, C. A. Toth, J. A. Izatt, and S. Farsiu, “Automatic segmentation of seven retinal layers in SDOCT images congruent with expert manual segmentation,” Opt. Express 18(18), 19413–19428 (2010). [CrossRef]  

11. L. D. Sisternes, G. Jonna, J. Moss, M. F. Marmor, T. Leng, and D. L. Rubin, “Automated intraretinal segmentation of SD-OCT images in normal and age-related macular degeneration eyes,” Biomed. Opt. Express 8(3), 1926–1949 (2017). [CrossRef]  

12. A. Yazdanpanah, G. Hamarneh, B. R. Smith, and M. V. Sarunic, “Segmentation of Intra-Retinal Layers From Optical Coherence Tomography Images Using an Active Contour Approach,” IEEE Trans. Med. Imaging 30(2), 484–496 (2011). [CrossRef]  

13. Y. He, A. Carass, Y. Liu, B. M. Jedynak, S. D. Solomon, S. Saidha, P. A. Calabresi, and J. L. Prince, “Deep learning based topology guaranteed surface and MME segmentation of multiple sclerosis subjects from retinal OCT,” Biomed. Opt. Express 10(10), 5042–5058 (2019). [CrossRef]  

14. J. Kugelman, D. Alonso-Caneiro, S. A. Read, S. J. Vincent, and M. J. Collins, “Automatic segmentation of OCT retinal boundaries using recurrent neural networks and graph search,” Biomed. Opt. Express 9(11), 5759–5777 (2018). [CrossRef]  

15. A. G. Roy, S. Conjeti, S. P. K. Karri, D. Sheet, A. Katouzian, C. Wachinger, and N. Navab, “Relaynet: retinal layer and fluid segmentation of macular optical coherence tomography using fully convolutional networks Biomed,” Biomed. Opt. Express 8(8), 3627–3642 (2017). [CrossRef]  

16. Y. He, A. Carass, Y. Liu, B. M. Jedynak, S. D. Solomon, S. Saidha, P. A. Calabresi, and J. L. Prince, “Structured layer surface segmentation for retina OCT using fully convolutional regression networks,” Med. Image Anal. 68, 101856 (2021). [CrossRef]  

17. J. Wang, Z. Wang, F. Li, G. Qu, Y. Qiao, H. Lv, and X. Zhang, “Joint retina segmentation and classification for early glaucoma diagnosis,” Biomed. Opt. Express 10(5), 2639–2656 (2019). [CrossRef]  

18. D. C. Fernández, N. Villate, C. A. Puliafito, and P. J. Rosenfeld, “Comparing total macular volume changes measured by Optical Coherence Tomography with retinal lesion volume estimated by active contours,” Invest. Ophth. Vis. Sci. 45(4), 3072 (2004).

19. F. Rossant, I. Bloch, I. Ghorbel, and M. Paques, “Parallel Double Snakes. Application to the segmentation of retinal layers in 2D-OCT for pathological subjects,” Pattern Recogn. 48(12), 3857–3870 (2015). [CrossRef]  

20. M. K. Garvin, M. D. Abramoff, X. Wu, S. R. Russell, T. L. Burns, and M. Sonka, “Automated 3-D Intraretinal Layer Segmentation of Macular Spectral-Domain Optical Coherence Tomography Images,” IEEE Trans. Med. Imaging 28(9), 1436–1447 (2009). [CrossRef]  

21. A. Lang, A. Carass, A. K. Bittner, H. S. Ying, and J. L. Prince, “Improving graph-based OCT segmentation for severe pathology in Retinitis Pigmentosa patients,” Proc. SPIE 10137, 426–433 (2017).

22. F. G. Venhuizen, B. V. Ginneken, B. Liefers, M. J. J. P. van Grinsven, S. Fauser, C. Hoyng, T. Theelen, and C. I. Sánchez, “Robust total retina thickness segmentation in optical coherence tomography images using convolutional neural networks,” Biomed. Opt. Express 8(7), 3292–3316 (2017). [CrossRef]  

23. S. L. Lou, X. Chen, J. Liu, Y. Shi, H. Qu, Y. Wang, and H. Cai, “Fast OCT image enhancement method based on the sigmoid-energy conservation equation,” Biomed. Opt. Express 12(4), 1792–1803 (2021). [CrossRef]  

24. R. D. Kinkelder, D. M. D. Bruin, F. D. Verbraak, T. G. V. Leeuwen, and D. J. Faber, “Comparison of retinal nerve fiber layer thickness measurements by spectral-domain optical coherence tomography systems using a phantom eye model,” J. Biophoton. 6(4), 314–320 (2013). [CrossRef]  

25. S. L. Lou, X. Chen, X. Han, J. Liu, Y. Wang, and H. Cai, “Fast Retinal Segmentation Based on the Wave Algorithm,” IEEE Access 8, 53678–53686 (2020). [CrossRef]  

26. J. Liu, S. L. Lou, X. Chen, H. Cai, and Y. Wang, “Fast Segmentation Algorithm for Cystoid Macular Edema Based on Omnidirectional Wave Operator,” Appl. Sci. 11(14), 6480 (2021). [CrossRef]  

27. E. W. Dijkstra, “A note on two problems in connexion with graphs,” Numer. Math. 1(1), 269–271 (1959). [CrossRef]  

28. Y. He, A. Carass, S. D. Solomon, S. Saidha, P. A. Calabresi, and J. L. Prince, “Retinal layer parcellation of optical coherence tomography images: Data resource for multiple sclerosis and healthy controls,” Data Brief 22, 601–604 (2019). [CrossRef]  

29. S. Farsiu, S. J. Chiu, R. V. O’Connell, F. A. Folgar, E. Yuan, J. A. Izatt, and C. A. Toth, “Quantitative Classification of Eyes with and without Intermediate Age-related Macular Degeneration Using Optical Coherence Tomography,” Ophthalmology 121(1), 162–172 (2014). [CrossRef]  

30. S. L. Lou, “WAVE-GT-OCTimage,” GitHub2022, https://github.com/loushiliang/WAVE-GT-OCTimage.

31. Q. Song, J. Bai, M. K. Garvin, M. Sonka, J. M. Buatti, and X. Wu, “Optimal multiple surface segmentation with shape and context priors,” IEEE Trans. Med. Imaging 32(2), 376–386 (2013). [CrossRef]  

32. A. Shah, L. Zhou, M. D. Abrámoff, and X. Wu, “Multiple surface segmentation using convolution neural nets: application to retinal layer segmentation in OCT images,” Biomed. Opt. Express 9(9), 4509–4526 (2018). [CrossRef]  

33. J. C. Stephanie, J. A. Michael, S. M. Priyatham, W. C. Scott, A. I. Joseph, and F. Sina, “Kernel regression based segmentation of optical coherence tomography images with diabetic macular edema,” Biomed. Opt. Express 6(4), 1172–1194 (2015). [CrossRef]  

Supplementary Material (1)

NameDescription
Code 1       WAVE-GT-OCTimage

Data Availability

Data is available in [28,29,33]. Other data may be obtained from the authors upon reasonable request. Codes underlying the results presented in this paper are available in Code 1 [30].

28. Y. He, A. Carass, S. D. Solomon, S. Saidha, P. A. Calabresi, and J. L. Prince, “Retinal layer parcellation of optical coherence tomography images: Data resource for multiple sclerosis and healthy controls,” Data Brief 22, 601–604 (2019). [CrossRef]  

29. S. Farsiu, S. J. Chiu, R. V. O’Connell, F. A. Folgar, E. Yuan, J. A. Izatt, and C. A. Toth, “Quantitative Classification of Eyes with and without Intermediate Age-related Macular Degeneration Using Optical Coherence Tomography,” Ophthalmology 121(1), 162–172 (2014). [CrossRef]  

33. J. C. Stephanie, J. A. Michael, S. M. Priyatham, W. C. Scott, A. I. Joseph, and F. Sina, “Kernel regression based segmentation of optical coherence tomography images with diabetic macular edema,” Biomed. Opt. Express 6(4), 1172–1194 (2015). [CrossRef]  

30. S. L. Lou, “WAVE-GT-OCTimage,” GitHub2022, https://github.com/loushiliang/WAVE-GT-OCTimage.

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 (10)

Fig. 1.
Fig. 1. Retinal layers of OCT image centered at the macula. The boundaries from top to bottom are ILM (I), NFL-GCI (II), GCI-INL (III), INL-OPL (IV), OPL-ONL (V), IS-OS (VI), IRPE (VII) and OBM (VIII).
Fig. 2.
Fig. 2. A schematic of the proposed method. Gray represents the preprocessing part, blue represents the bidirectional wave algorithm part, and green represents the improved graph theory part.
Fig. 3.
Fig. 3. Diagram of the preprocessing part. Red represents (W/4×H/4) size images, orange represents (W/2×H/2) size images and pink represents the original size image. W is the width of the image and H is the height of the image.
Fig. 4.
Fig. 4. Improvement of velocity definition. (i, j) is the pixel to be detected in the image, and I(i, j) is the gray value of the pixel in the image. The 3 × 3 window centered on (i, j) is the Z template. The 3 × 2 window overlap with the left end of the Z template is the Q template. The 3 × 5 window centered on (i, j) is the M template. All the templates in this paper are moved from left to right (direction of movement). The normalized difference value in the Z template is defined as v, and the normalized difference value in the Q template is defined as vq.
Fig. 5.
Fig. 5. Schematic diagram of improved wave algorithm. (a∼c) Potential energy map, (d∼e) Boundary Map, (f) Guide map. The improved wave potential energy equation can effectively detect two types of areas, namely the area where the gray value changes from dark to light (yellow) and the area where the gray value changes from light to dark (cyan). The improved potential energy correction equation can detect two types of boundary lines from two types of areas.
Fig. 6.
Fig. 6. Network definitions in improved graph theory. Green represents nodes (BPEC≠0), Yellow represents the pathway (edges), Gray represents the background area, and Blue represents the gray value change areas (target area) detected by the bidirectional wave algorithm. “+” represents positive areas and positive edges, and “-” represents negative areas and negative edges.
Fig. 7.
Fig. 7. Automatic segmentation of the 8-layer boundary of the retina. (a) Nodes map, (b) Frame map, (c) Guide map. Before obtaining the Nodes map, Frame map and Guide map, the images of three sizes are uniformly transformed to the W/2×H/2 size.
Fig. 8.
Fig. 8. Automatic segmentation results of the proposed method in [28] and [29]. The proposed method (red) and expert manual tracings (green).
Fig. 9.
Fig. 9. Automatic segmentation results of the proposed method in [33]. (a) Advanced DME cases, (b) Intermediate DME cases, (c) Mild DME cases. The proposed method (red) and expert manual tracings (green).
Fig. 10.
Fig. 10. Tests of the proposed method on different OCT imaging systems. (a) point scan images (size of B-scan is 400 × 400), (b) line-scan images (size of B-scan is 400 × 510). The red line is the segmentation result of the proposed method.

Tables (3)

Tables Icon

Table 1. The mean absolute distance (MAD(Std) in pixel) in [28] and [29]. Since only 3-layer in [29] are manually delineated (ILM, IRPE, and OBM), the MAD (Std) of the unlabeled layers in the table is “–”

Tables Icon

Table 2. The mean absolute distance (MAD) and root mean square error (RMSE) in [28] for three of the tested methods (in pixel)

Tables Icon

Table 3. The mean absolute distance (MAD(Std) in pixel) in [29] for three of the tested methods

Equations (5)

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

W P E  =  v × v q × σ + φ g
{ s = i M j M I ( i , j ) I ( i , j 1 ) | I ( i , j ) I ( i , j 1 ) | tanh ( s ) = e s e s e s + e s
B W P E  =  tanh ( s ) × ( V 2 × σ + φ g )
B P E C  =  { + 1 , ( ( K > 1 ) ( C > 1 ) ) ( B W P E α ) 1 , ( ( K < 1 ) ( C < 1 ) ) ( B W P E ( α ) ) 0 , o t h e r s
W  =  i L w i
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.