Existing methods for choroid segmentation in OCT are mainly divisible into two categories: traditional methods, and machine learning methods. Zhang et al. (Zhang et al., 2012) first attempted to extract the choroidal layer in 3D SD-OCT by adapting a graph-based method, which produced a relatively accurate choroidal surface. However, this method was tested on normal subjects only, and it is difficult to achieve the expected segmentation performance for some patients, especially those with large changes of choroidal morphology. In order to overcome this limit, Hu et al. (Hu et al., 2013) improved the graph-based multi-layer segmentation method by applying various smoothness and interaction constraints to different choroidal layer structures. This method has been validated on OCT images collected from both healthy subjects and non-neovascular AMD subjects, revealing great similarities with manual segmentation.

With the emergence of Enhanced Depth Imaging OCT (EDI-OCT) technology, its high-resolution imaging made the choroidal layer structure more clearly displayed in OCT B-scans, which is more conducive to choroid segmentation. Tian et al. (Tian et al., 2013) adopted Dijkstra’s algorithm to seek the shortest path, and the choroidal surface was quickly and accurately detected. Similarly, Danesh et al. (Danesh et al., 2014) proposed a segmentation method based on the Gaussian mixture model, to obtain the choroidal structure in EDI-OCT images. However, this still requires handcrafted features, and is sensitive to noise artifacts existing in EDI-OCT images. In addition, Chen et al. (Chen et al., 2015) introduced a new pipeline composed of a progressive intensity distance image generation algorithm and graph search method for the problem of noise and boundary ambiguity. Wang et al. (Wang et al., 2017) used a Markov Random Field (MRF) method to connect adjacent pixels, and a level set method to regularize the distance of uneven textures.

Since graph search technology is greatly affected by manual parameter settings, deep learning-based methods have been developed to obtain the choroidal structure. Sui et al. (2017) proposed a convolutional neural network (CNN)-based method that learns a graph-edge weight directly from raw OCT pixels. The network structure can be divided into two parts: one detects the CSI boundary, and the other detects the BM boundary. This method has revealed good adaptability to 3D EDI-OCT images collected from both healthy subjects and patients with macular edema. He et al. (2021) combined CNN and a l ₂-l _q (0 < q < 1) fitter to segment the outer choroidal surface, in which the CNN is used to generate predicted values, and the l ₂-l _q fitter is employed to maintain the stability of the fitting function. The OCT image is partitioned into small patches to form the input of the CNN, and post-processing is required to discard irrelevant information, which leads to relative inefficiency compared to end-to-end architectures. Similarly, Masood et al. (2019) used deep learning methods to establish a new segmentation structure to obtain the outer choroidal surface. Before being fed into the CNN, the OCT image needs to be divided into patches for data sampling and conversion. This method reduces the average segmentation error, but it is still not as efficient as end-to-end architectures.

As a result, several end-to-end deep learning approaches have been proposed more recently. Zhang et al. (2020) proposed an end-to-end method, consisting of a global multi-layer segmentation block, a choroidal layer segmentation block, and a regularization block. This first segments all the inner retinal layers, and then utilizes global information to detect the choroidal layer. For 3D OCT images collected from healthy subjects, the thickness difference obtained by this method (4.30 ± 0.02 pixels) is more accurate than those obtained by other state-of-the-art methods. Chai et al. (2020) proposed a method that can effectively segment the choroidal boundary by minimizing the differences between different regions, and takes into account the differences between different OCT acquisition equipment. It feeds OCT images from different domains into a U-Net-based network, and uses both adversarial and perceptual loss for domain adaptation.

Examining the choroidal layer as extracted from OCT images, ophthalmologists can analyze choroidal variations from different perspectives. In particular, choroidal thickness is of great interest, as it often indicates the presence or even severity of some ophthalmic diseases. Yiu et al. (2015) extracted choroidal thickness from EDI-OCT images collected from subjects with Age-related Macular Degeneration (AMD). Employing on a semi-automatic segmentation method, they analyzed the similarities and differences in choroidal thickness between normal individuals and patients with AMD. Wang et al. (2015) compared choroidal thickness between patients with high myopia and healthy people. By analyzing the experimental results, they found that choroidal thickness in healthy individuals is significantly thicker than that of individuals with high myopia. Regatieri et al. (2012) examined choroidal thickness in diabetic patients and found that the change of choroidal thickness was related to the severity of diabetes. More recently, several studies have shown that choroidal thickness as revealed by retinal OCT images is associated with certain neurodegenerative diseases. Moschos and Chatziralli (2018) extracted the retinal thickness and choroidal thickness of patients with Parkinson’s Disease (PD) from spectral domain OCT, and compared the results with those from healthy individuals. The differences between people with, and without PD were statistically significant. Similarly, Satue et al. (2018) used swept-source OCT to measure retinal and choroidal thickness of patients with PD. They found that the retina of patients with PD became thinner, while choroidal thickness might increase. Similar to our work, Chen et al. (2022) segmented the choroidal layer of highly myopic patients and non highly myopic people and compared the thickness, while they lacked the analysis of three-dimensional features, and the segmentation performance needs to be improved.

To this end, the automatic and accurate quantification of choroidal thickness is potentially crucial to diagnosis of these diseases. However, most quantification approaches of choroidal thickness can only provide two-dimensional measurements at a fixed location, which limits the practicability. Therefore, we proposed to use 3D edge point cloud features to produce a three-dimensional reconstruction of the choroidal layer.

Boundary segmentation in images remains a research hotspot, not only in the fields of medical image analysis but also in many fields of other computer vision such as remote sensing. The mainstream boundary segmentation approaches may be divided into two categories: filtering-based methods, and learning-based methods. Wang et al. (2018) introduced an interactive geodesic method based on CNN into medical image segmentation: a mannual correction of boundary information is required to improve the accuracy of boundary segmentation. Lee et al. (2020) proposed a novel network with boundary preserving blocks to retain the boundary information via learning proper weights of boundary features. Wei et al. (2021) proposed a concentric loop CNN with a boundary detector and a refinement block to improve the effect of boundary segmentation in remote sensing images. Dang and Lee (2021) improved the effect of boundary segmentation in document images by sharing the weights of boundary features and global features, and using adversarial loss to strengthen the learning of boundary information. Wang et al. (2021) combined a transformer with CNN to enhance the segmentation of skin lesions, and used an attention mechanism to boost the performance of boundary segmentation. Recently, Yu et al. (2022) proposed FBCU-Net, which uses boundary semantic features to segment medical images, but it is mainly used for region segmentation and the performance of layer structure segmentation still needs to be improved.

Inspired by Lee et al. (2020), we employed the boundary enhancement point selection algorithm to select several key points, then generated a point map for each B-scan as another ground truth for training. By contrast with the binary disks generated on the selected points in (Lee et al., 2020), we adopted a two-dimensional Gauss function to generate a point map with soft boundaries for more effective guidance with boundary information to segmentation. This modification is mainly based on the following considerations: a binary disk allocates undifferentiated attention to all pixels in the neighborhood of the corresponding selected point, which creates vulnerability to deviation of boundary localization. The soft point map can be expressed as follows: S i , j = max k ∈ 1 , … , K exp − i − x k 2 + j − y k 2 2 σ 2 (1) where (i, j) represent the coordinates of one pixel of the generated soft point map matrix S; (x _k , y _k ) represent the coordinates of kth selected key point (total K selected key points); σ represents the standard deviation of the Gauss function. The differences between the original point map and the proposed soft point map are illustrated in Figure 3.

The proposed choroid segmentation architecture incorporates our novel BEMs into its encoder/decoder layers, as shown in Figure 2. Figure 4 presents the architecture of the BEM, which consists of three branches including FEB, CEB and BAB. The BEM can be embedded in various layers in the segmentation network. The BEM embedded in the ith layer takes the feature maps f i ∈ R w i × h i × c i as input, where w ⁱ , h ⁱ , and c ⁱ respectively represent the width, height and channels of the feature maps at the ith layer. The FEB produces the boundary-enhanced point map M i ∈ R w i × h i × 1 . The CEB generates a channel-wise weighting vector N i ∈ R c i . The BAB outputs a single-channel activation map Q i ∈ R w i × h i × 1 . Then, the final output feature maps v i ∈ R w i × h i × c i is calculated as follows: v i = f i ⊕ f i ⊗ M i ⊗ Q i ⊗ N i , (2) where ⊕ represents element-wise addition; ⊗ represents multiplication (pixel-wise multiplication for single-channel maps M ⁱ and Q ⁱ , and channel-wise multiplication for the weighting vector N ⁱ ).

The gray section in Figure 4 shows the architecture of the FEB, which is designed for extracting a boundary-enhanced point map with different receptive fields. Multiple receptive fields integrate global context information and local detailed information, which is beneficial to accurate localization of boundary key points. We adopted convolution with different dilation rates to obtain the features with different receptive fields. Finally, the features with different receptive fields are concatenated and then fed into a single 1 × 1 convolutional layer with a Sigmoid function. Let d r s ( f i ) be the encoded feature maps of input feature maps f ⁱ using s × s convolution with dilation rate r. Then the generated boundary-enhanced point map can be expressed as: M i = S d 1 1 d 1 1 f i ⊙ d 1 3 f i ⊙ d 2 3 f i ⊙ d 4 3 f i ⊙ d 6 3 f i , (3) where ⊙ and S denote a concatenation operation, and a Sigmoid function, respectively.

Each channel of a feature map may be regarded as a specific-class response. However, there are also some differences in the importance of different feature classes to a specific task (e.g., choroid segmentation). In order to selectively enhance features useful for choroid segmentation, a CEB is designed to calculate channel-wise weighting vectors for extracted feature maps. The detailed structure of CEB is shown in the green section of Figure 4. The branch first adopts global average pooling (gap) to generate channel-wise statistics of the input feature maps f ⁱ , then applies one-dimensional convolution of kernel size 3 to obtain the final channel-wise weighting vectors N ⁱ , which can be denoted as follows: N i = S C 1 D 3 f gap f i , (4) where f _gap denotes global average pooling; C1D ₃ denotes the one-dimensional convolution of kernel size 3; S denotes a Sigmoid function.

In order to further improve detection of the choroidal boundary, we introduced an extra branch called the BAB, as detailed in by the red section of Figure 4. The channel number of the feature maps is reduced to 1 via a 1 × 1 convolutional layer followed by a Sigmoid function: a specific activation function is then applied to the obtained single-channel feature map, which can be formulated as follows: Q i = e x i − 0.5 2 + 1 − e − 0.25 , x i = S C 2 D 1 f i , (5) where f ⁱ represents the input feature maps; C2D ₁ represents the 1 × 1 convolutional layer; S denotes a Sigmoid function.

It is worth noting that the specific activation function was designed based on the observation that choroidal boundary pixels generally have a value around 0.5 in the obtained feature maps, while interior of choroid and background pixels have values near 1 and 0, respectively. To this end, the activation map is utilized to adjust each feature map spatially for the choroid segmentation task, by assigning higher weights for choroidal boundary pixels (close to 2—e ^−0.25), and lower weights for interior of choroid and background pixels (close to 1) (Simonyan and Zisserman, 2014). In this way, the boundary information is highlighted after activation.

In order to effectively train the proposed choroid segmentation network, we introduced a novel loss called boundary perceptual loss (BP-Loss), which embeds the soft point map into the segmentation network. The complete joint loss function is formed after incorporating segmentation loss.

First, we adopted Binary Cross Entropy Loss as the segmentation loss in order to reduce the difference between the ground-truth segmentation map and the predicted segmentation map, which is defined as: L o s s Seg = − ∑ i 1 − S G T i ⋅ log 1 − S ̂ Pred i + S G T i ⋅ log S ̂ Pred i , (6) where S G T i and S ̂ Pred i represent the ith pixel of ground truth segmentation map, and the corresponding predicted segmentation map, respectively.

Unlike general semantic segmentation tasks, medical images require strong expert knowledge to achieve better segmentation performance. Therefore, we concatenated the generated soft point map with the predicted result and ground-truth to form the input of the VGG network, in order to constrain their geometrical relationship. To maintain consistency between the output of the FEB and ground truth. We adopted mean square error (MSE) loss and defined boundary point loss as: L o s s Map i = 1 h i × w i ∑ j = 1 h i ∑ k = 1 w i M j , k i − M G T j , k i 2 , (7) where M ⁱ and M G T i respectively represent the output point map and the ground truth soft point map for the FEB in the ith layer, and h ⁱ and w ⁱ represent the corresponding height and width, respectively. L o s s G F = ∑ i 1 N i ‖ ϕ i S G T ⊙ M G T − ϕ i S ̂ Pred ⊙ M G T ‖ 1 , (8) where ϕ _i (⋅) denotes the feature maps from the ith layer of the VGG-19 network pre-trained on the ImageNet; S _GT denotes ground truth segmentation map; M _GT denotes ground truth soft point map; S ̂ Pred denotes the predicted segmentation map; ⊙ denotes concatenation operation; and N _i denotes the element number of feature maps from the ith layer of VGG-19.

The Boundary Perceptual Loss is defined as: L o s s B P = L o s s G F + ∑ i n L o s s Map i (9) where n indicates the number of BEMs in the proposed choroid segmentation network.

Finally, the total loss function is defined as: L o s s Total = λ Seg L o s s Seg + λ B P L o s s B P (10) where λ _Seg and λ _BP are set as 0.5 and 0.5 in our task.

In this work, a new Choroidal OCT image for SegmenTAtion (COSTA) dataset, which consists of three subsets named COSTA-H, COSTA-T and COSTA-B, was constructed for our proposed approach. These subsets were acquired from different devices or adopted different bit depths, as illustrated in Figure 5.

• COSTA-H consists of 10 OCT volumes from 10 healthy subjects. Each volume was captured by the Heidelberg Spectrails system, and contains 384 non-overlapping B-scans, each covering a 3 × 3 × 2 mm ³ region. Two groups of ophthalmologists were invited to independently make manual annotations of the upper and lower boundaries of the choroidal layer (BM and CSI), and their consensus were used as ground truth after discussion. In order to reduce the number of manual annotations, we asked these ophthalmologists to annotate one B-scan every six consecutive B-scans, due to the high similarity between adjacent B-scans in an OCT volume. Finally, we obtained a total of 384/6 × 10 = 640 B-scans with manually annotated choroid boundaries.

• COSTA-T was captured by the Topcon DRI-OCT-1 system, containing a total of 20 OCT volumes from 20 healthy human eyes. Each volume contains 256 B-scans with a resolution of 512 × 992 pixels covering a 6 × 6 × 2 mm ³ region. This dataset was also annotated by the same protocol as COSTA-H, yielding a total of 256/4 × 20 = 1280 B-scans with annotated choroid boundaries.

Both COSTA-H and COSTA-T datasets were used for training and testing, where the ratio of data volume between the training and testing sets is 3:1. To more accurately and credibly evaluate the proposed network, we adopted the 4-fold cross-validation strategy, i.e., randomly dividing all samples into 4 equal pieces and taking each piece as the validation set and others as the training set in turn. After 4 groups of tests, different validation sets are replaced each time. That is, the results of four groups of models are obtained, and the average value is taken as the final result.

• COSTA-B was captured by a homemade 70-Khz SD-OCT system with different bit depths, and all this data was selected from Hao et al. (Hao et al., 2020). It contains 199 annotated B-scans with a resolution of 270 × 450 pixels from one normal subject. By contrast with COSTA-H and COSTA-T, COSTA-B was only used to test the robustness of the proposed method with respect to imaging quality. For the same B-scan, we also made a comparison of segmentation results based on different bit depths. The higher bit depth represents the better image quality, which makes for easier choroid segmentation.

All deep learning approaches in the experiments were implemented with PyTorch (Paszke et al., 2019) and ran on a single NVIDIA GeForce GTX 3090 GPU with 24 GB memory under an Ubuntu 16.04 system. The proposed network was trained with 400 epochs, and some hyper-parameters were set as follows: Adam optimization, with an initial learning rate of 0.0005 and batch size of 8. For other comparison methods, we adopted the same training strategy in the original paper.

In order to compare the performance of the proposed method with other state-of-the-art deep learning networks, the following routine metrics for image segmentation were adopted and calculated:

• Dice Coefficient (Dice) = 2 × T P 2 × T P + F P + F N ;

• Intersection over Union (IoU) = T P T P + F P + F N ;

• Accuracy (Acc) = T P + T N T P + T N + F P + F N ;

• Sensitivity (Sen) = T P T P + F N ;where TP is true positive, FP is false positive, TN is true negative, and FN is false negative.

In addition, we adopted Average Unsigned Surface Detection Error (AUSDE) (Xiang et al., 2018) based on BM and CSI: A U S D E = 1 m ∑ i = 1 m | y i − y ̂ i | (11)

where m represents the width of the B-scan, y ⁽ⁱ⁾ and y ̂ ( i ) represent vertical coordinates of the ith point in the horizontal direction of BM or CSI in the predicted segmentation map and ground truth, respectively. Based on y ⁽ⁱ⁾ and y ̂ ( i ) of BM and CSI (respectively denoted as y BM ( i ) , y ̂ BM ( i ) and y CSI ( i ) , y ̂ CSI ( i ) ), average Thickness Difference (TD) can also be calculated as: T D = 1 m ∑ i = 1 m ‖ y CSI i − y BM i | − | y ̂ CSI i − y ̂ BM i ‖ (12)

Figure 6 shows the visualization of the results of training and testing on Heidelberg, and we compare them with other popular methods that use deep learning to segment the choroid layer, including the U-Net (Ronneberger et al., 2015), FCN (Long et al., 2015), DeepLab v3+ (Chen et al., 2018), SegNet (Badrinarayanan et al., 2017), CE-Net (Gu et al., 2019), CS-Net (Mou et al., 2019), and SCA-CENet (Mao et al., 2020). In Figure 7, the BM and CSI are located by blue and red lines, respectively. In both the overall segmentation and the boundary extraction (shown by the zoomed-in part of Figure 7), the proposed BEM performed better than its counterparts, which shows that our reinforcement of boundary characteristics is useful and efficient.

The testing results on the COSTA-B dataset are presented in Figure 8, which shows three OCT B-Scans with different bit depths, including 6, 8, and 12 (the best image quality) bits. The model with the best segmentation results on the COSTA-H dataset was used to segment these images. The segmentation results show the excellent robustness of the proposed model. Common segmentation methods that focus on global information and lack detailed features have difficulty in fully segmenting all the choroid layers. In contrast, the proposed method, which benefits from its ability to enhance boundary characteristics and extract different features from different perceptual fields and channels, shows good robustness across images of differing qualities.

In order to verify the advantages of our method from a quantitative perspective, we selected the Dice, IoU, AUSDE and TD as evaluation metrics. For the COSTA-B test, only two metrics, Dice and IoU, were selected for evaluation, as the segmentation results of many methods could not form clear boundary lines (as shown in Figure 8).

Table 1 shows the quantitative comparison of various deep learning methods applied to the COSTA-H dataset. As is shown in Table 1, our method achieved 97.06% of the Dice coefficient and 94.31% of the IoU value, with a Boundary Error of 0.9496 for AUSDE of the BM and a Detection Error of 3.0029 for the lower boundary CSI, outperforming other methods. After adding BEM and BP-Loss to UNet, the value of AUSDE of CSI is improved by 2.57 pixels and that of TD is improved by 2.68 pixels, demonstrating the value of boundary extraction.

In order to further evaluate our proposed method, we conducted additional tests on the COSTA-T dataset and the results are shown in Table 2. With the help of BEM, our method achieves 92.87% of the Dice coefficient and 86.91% of the IoU value: after adding the BEM and BP-Loss, this improves to 2.02% and 2.39% on the Dice and IoU values, respectively.

We also select 6, 8, 10, and 12 bit depth images in the COSTA-B dataset for evaluation on the DICE and IoU metrics. For each network, both the best trained and validated models on the COSTA-H dataset were tested. It may be seen from Figure 9, by contrast with our proposed method, the segmentation results of the other networks degrade significantly over decreased depths of the same image. Specifically, all of the tested algorithms obtain better segmentation results on the 12-bit depth map than those on the 6 and 8-bit depths. This further validates the robustness of the proposed method.

In order to verify that each branch in the BEM and the BP-Loss are effective, we conducted ablation experiments by removing each branch separately and performing the experiments on the same dataset, and the results are shown in Figure 10. It can be seen that the dice and IoU metrics gradually increase with the addition of different branches, both in COSTA-T and COSTA-H datasets. It is obvious that the segmentation result benefits from every branch of the BEM.

After adding BEM, the IoU of the segmentation result on COSTA-T dataset reaches 86.53%, which is 3.01% higher than the baseline. With the help of BP-Loss, the segmentation result reaches 86.91%, which further improves the segmentation effect. Similarly, each module and branch plays a role on the COSTA-H dataset.