OCT segmentation: Integrating open parametric contour model of the retinal layers and shape constraint to the Mumford-Shah functional
In this paper, we propose a novel retinal layer boundary model for segmentation of optical coherence tomography (OCT) images. The retinal layer boundary model consists of 9 open parametric contours representing the 9 retinal layers in OCT images. An intensity-based Mumford-Shah (MS) variational functional is first defined to evolve the retinal layer boundary model to segment the 9 layers simultaneously. By making use of the normals of open parametric contours, we construct equal sized adjacent narrowbands that are divided by each contour. Regional information in each narrowband can thus be integrated into the MS energy functional such that its optimisation is robust against different initialisations. A statistical prior is also imposed on the shape of the segmented parametric contours for the functional. As such, by minimising the MS energy functional the parametric contours can be driven towards the true boundaries of retinal layers, while the similarity of the contours with respect to training OCT shapes is preserved. Experimental results on real OCT images demonstrate that the method is accurate and robust to low quality OCT images with low contrast and high-level speckle noise, and it outperforms the recent geodesic distance based method for segmenting 9 layers of the retina in OCT images.
Optical coherence tomography (OCT) image segmentation to detect retinal layer boundaries is a fundamental procedure for diagnosing and monitoring the progression of retinal and optical nerve diseases. There exist rich literature on approaches for automatic and semi-automatic OCT image segmentation. Common methods include deformable models [1, 2], graph-based and geodesic distance methods [3, 4], statistical shape and appearance models [5, 6], etc. Very recently, deep neural networks are becoming increasingly popular for OCT segmentation, demonstrating excellent performance [7, 8]. However, these deep learning methods usually require the networks to be sufficiently deep to learn all appearance and shape variations of the retinal layers from annotated training sets. Therefore the training set has to be very large and rich to prevent the over-fitting. Large annotated training sets are however difficult to obtain. Existing OCT segmentation algorithms also tend to segment individual retinal layers separately. This form of analysis often fails when there is uncertainty in the image, especially some retinal layers are often difficult to see or missing in OCT images.
We believe that OCT segmentation is more effective by incorporating anatomical shape of the retinal layers and their spatial relations. As such, in this paper we propose a shape-based variational Mumford-Shah (MS) functional for segmentation of up to 9 retinal layer boundaries in OCT images, using only a small training set. We make three distinct contributions to OCT segmentation:
We introduce a new piecewise constant variational MS functional to evolve a pre-defined retinal layer boundary model for OCT image segmentation. It has a region-based data fidelity term and a hybridised first and second regularisation term. The pre-defined retinal layer boundary model consists of 9 retinal layer boundaries, each is an open explicit parametric contour represented by a set of control points. We then construct two narrowbands around each open contour, within which region-based information is derived to aid contour evolution. We show that by incorporating a retinal layer boundary model our method can segment 9 retinal layers simultaneously, and by utilising regional information, the proposed method has a large convergence range and is robust to initialisation.
We introduce to the MS functional a shape constraint learnt from a set of training OCT shapes. We then apply the principal component analysis (PCA) to derive the statistical distribution from the training shapes as well as the resulting irregular contours evolved directly from the MS functional. In this way, the irregular contours are restricted to a manifold of familiar shapes and thereby can be pulled back to appropriate positions to allow a faster convergence.
We apply the proposed method to real OCT dataset acquired from healthy subjects, and demonstrate that the proposed method outperforms the state-of-the-art methods.
2.1 Intensity-based variational MS functional
We start with a new intensity-based MS segmentation functional, and then apply a learnt shape constraint to the functional. We define a retinal layer boundary model as having 9 retinal layers, each of which is an explicit parametric contour . Based on , we propose a new piecewise constant MS functional for each of the 9 retinal layers
where is the input image and are the mean grey value of region . The first regional energy term follows the Chan-Vese model [2, 10]. In the second term, is the first order derivative of the curve with respect to the arc length normalized into the region , and is the second order derivative. and are two regularisation coefficients. To segment a retinal layer in OCT, has two values, namely, . The contour regularisation (last two terms) combines the first and second order derivatives, preventing the contour from bending by introducing elasticity and stiffness to the contour. Each control point of the contour thus can be more equidistant or centred between its neighbourhood points. This makes the functional stable for numerical calculation. is defined as a parametric contour of a set of control points
Here , are the coordinates of the control points to represent one retinal layer boundary. The start and end points are and , respectively.
Variational methods taking into account regional information are generally more robust against noisy and initialisation. In OCT segmentation, we however require the parametric contour to be an open curve . In this case, is no longer equal to . This brings difficulties in estimating parameter in (1) due to the vanished image regions and . To circumvent this limitation, we propose to construct a narrowband around the open contour using the method illustrated in Fig 1:
2.2 Minimisation of MS by gradient descent
We minimise the proposed functional (1) with respect to both the parameter and the contour , where denote the mean grey values on both sides of contour curve and have a simple closed-form
where are the two regions partitioned by the curve and they change as the curve evolves. Since we are able to construct regions from an open contour , it is easy to calculate in (3). After are estimated, we fix and use the gradient descent method to minimise the functional (1) with respect to the open contour . This results in the following contour evolution equation
where , and denotes the outer normal vector of the contour . The first two terms on the right-hand side of (4) minimise the contour length and thereby enforce an equidistant spacing between the control points. The third term maximises the homogeneity in the adjoining regions in the narrowband, which is measured by the energy density (3). This forces move towards a retinal layer boundary in the OCT image.
We now need to discretise (4) for numerical implementation. Specifically, for each control point of the contour , we have the following two semi-implicit iterative schemes
where is an artificial time parameter. The superscript denotes the th iteration, and the subscript the th control point. As compared to explicit iterative methods, semi-implicit relaxation methods allow the use of a relative larger time step , speeding up the convergence rate of contour evolution. We apply the finite difference method with Neumann boundary condition  to discretise the 2nd- and 4th-order derivatives , , and .
To iterate (5) and (6), the open contour normals are calculated in Fig 2. Apart from the two endpoints, for each control point on the contour we first compute its tangent using its neighbours via the centre finite difference scheme. For the start and end points, their tangents are respectively computed with the forward and backward finite differences. The contour normals thereby can be easily derived from the computed tangents because they are perpendicular to each other.
Note that given control points we have to calculate a set of 2 equations as in (5) and (6). Since the 2nd- and 4th-order derivatives are linear differential operators, we can rewrite (5) and (6) to its matrix form as
where is the identity matrix, denotes pointwise multiplication and A is an matrix of the form
2.3 Integrating shape constraint to intensity-based MS
The contour evolution driven by (7) and (8) can become irregular as the MS functional uses only local pixel intensities. To further improve the MS functional we propose to impose some constraint for the functional to control contour evolution. With some manually segmented retinal layer boundaries as training shapes we can derive the statistical distribution of the retinal layer boundaries as a shape constraint on the MS functional. The main idea of this is to project the irregular shapes into a latent space spanned by a few eigenvectors of training shapes with largest eigenvalues. In this way, the resulting irregular shapes are restricted to a manifold of familiar shapes and thereby can be corrected. Note that other dimensionality reduction methods can be also used to derive shape constraint. In this work we found PCA performs reasonably well. Next we detail the use of shape constraint to (1) using PCA. Note that an OCT shape model represents the locations of 9 boundaries.
If one OCT shape model consists of control points, it can be modelled as a -dimensional vector , where ( appears in (2)). Assuming that we have training shapes manually annotated from OCT images. We first align all shapes at the centre of coordinate origin by applying the Procrustes transformations. The mean OCT shape is then computed by . For each shape in the training set, its deviation from the mean shape is . Then the covariance matrix Cov can be calculated by .
We calculate the eigenvectors and corresponding eigenvalues of Cov (sorted so that ). Let be the matrix of the first eigenvectors of . Then we can approximate an OCT shape in the training set by
where defines the parameters for different deformation patterns. Since P is an orthogonal matrix, is given by
By varying the parameters , we can generate new examples of the shape. We can also limit each to constrain the deformation patterns of the shape. Typical limits  are
The use of shape constraint to (1) for segmenting 9 retinal layer boundaries is summarised as follow: P, and are first computed from a set of training OCT shapes. After the original OCT shape model (locations of 9 boundaries) is updated with (7) and (8) (each boundary in the shape model is updated separately), it is transformed to approximate the mean shape with the Procrustes transformation, which forms . This is followed by updating using (10) and then projecting it using (11), i.e. the coefficients obtained with (10) are constrained with equation (11). Next, we correct the original OCT shape using (9) with the projected and warping the corrected version back to the original location. In this way, the original irregular OCT shape is pulled back to a regular one used for the next iteration of (7) and (8). The whole process is repeated until convergence. The overall numerical optimisation of the proposed method is fast as it only computes a few hundred control points.
3 Experimental results
Data: 9 retinal layers boundaries were segmented for each image from in vivo OCT B-scans. 30 Spectralis SDOCT (ENVISU C class 2300, Bioptigen, axial resolution = 3.3, scan depth = 3.4, 32,000 A-scans per second) B-scans from 15 healthy adults were used for the work. The B-scan was imaged from the left and right eye of the healthy adults using a spectral domain OCT device with a chin rest to stabilise the head. The B-scan located at the foveal centre was identified from the lowest point in the foveal pit where the cone outer segments were elongated (indicating cone specialisation). The ground truth boundaries were manually generated by an experienced ophthalmologist. The dataset was randomly split into 20 training and 10 (5 are corrupted with high-level speckle noise using Matlab function with 0.8 variance) validation datasets. For image pre-processing, all images were cropped to extract only region of interest and were flattened using ground truth labels before training. The 9 retinal layer boundaries which can be segmented by the proposed method are shown in Fig 3.
|Notation||Name of retinal boundary/surface||Abbreviation|
|internal limiting membrane||ILM|
|outer boundary of the retinal nerve fibre layer||RNFL|
|inner plexiform layer-inner nuclear layer||IPL-INL|
|inner nuclear layer-outer plexiform layer||INL-OPL|
|outer plexiform layer-outer nuclear layer||OPL-ONL|
|outer nuclear layer-inner segments of photoreceptors||ONL-IS|
|inner segments of photoreceptors-outer segments of photoreceptors||IS-OS|
|outer segments of of photoreceptors-retinal pigment epithelium||OS-RPE|
|retinal pigment epithelium-choroid||RPE-CH|
Parameters: These parameters in our method are the first and second order regularisation parameters and , the artificial time step , the narrowband radius , the number of eigenvectors/eigenvalues , the number of control points for 9 boundaries, and the iteration number (1000).
These parameters are selected as follows: 1) Large and leads to increasingly shorter and finally vanishing segmentation boundaries. Small values may cause control points intersect with each other, thus leading to numerical instabilities. In this work, we fixed ; 2) is bounded by the CFL stability condition. Numerical stabilities can be attained by using ; 3) was selected according to the initial OCT boundaries. If initialisation is close to the true retinal layer boundaries, is small ( pixels for Fig 4 1st row). Otherwise, a large should be used ( pixels for Fig 4 last row). Note that we set the same value for each of 9 OCT boundaries; 4) was confirmed by choosing the first largest eigenvalues such that , where is the total variance of all the eigenvalues. We used the number of training samples as the value of . Eigenvectors corresponding to small eigenvalues do not contribute much to shape variation; 5) 360 control points in total are used for a whole OCT shape (each retinal layer boundary thus has 40). Overall, we only adjusted for different initialisations (see Fig 4).
Comparsion: The performance of the proposed segmentation method was evaluated by computing the Hausdorff distance (HD) metric between the automated and ground truth segmentations for different retinal layer boundaries. We compared our method with the geodesic distance method proposed in . In Fig 5, visual comparison suggests that the proposed method provides significant improvements over intensity-based geodesic distance in OCT segmentation, especially when OCT images are of low contrast and contain high-level speckle noise. In Table 2, we report the HD metric of each boundary as well as the total 9 boundaries over the validation dataset and show that our method outperforms the geodesic distance in terms of the HD metric for all the retinal layer boundaries. The improvements are more evident at IML, RNFL, IPL-INL, INL-OPL and OPL-ONL boundary locations. Note that in order to compare our results with ground truth labels we fitted a spline curve to each segmentation contour such that the resulting contours run across the entire width of the OCT image.
We presented a new segmentation method for optical coherence tomography (OCT) images, which allows the integration of statistical shape models learned from a small OCT dataset. To this end, we developed the Mumford-Shah functional in a way which facilitates a parametric representation of open contours. We then constructed narrowbands around the open contours such that regional information can be derived to assist segmentation. We have shown that integrating such information allows the proposed method to have a large convergence range and thus robust against different initialisations. We have also validated that the proposed method is very accurate even OCT images are of low contrast and contain high-level speckle noise, and that the method outperforms the state-of-the-art geodesic distance segmentation method. The proposed method can be readily extended to other segmentation problems involving open contours.
-  Rossant, F., Bloch, I., Ghorbel, I., Paques, M.: Parallel double snakes. application to the segmentation of retinal layers in 2d-oct for pathological subjects. Pattern Recognition 48(12) (2015) 3857–3870
-  Chan, T.F., Vese, L.A.: Active contours without edges. IEEE Transactions on image processing 10(2) (2001) 266–277
-  Chiu, S.J., Li, X.T., Nicholas, P., Toth, C.A., Izatt, J.A., Farsiu, S.: Automatic segmentation of seven retinal layers in sdoct images congruent with expert manual segmentation. Optics express 18(18) (2010) 19413–19428
-  Duan, J., Tench, C., Gottlob, I., Proudlock, F., Bai, L.: Automated segmentation of retinal layers from optical coherence tomography images using geodesic distance. Pattern Recognition 72 (2017) 158–175
-  Xie, W., Duan, J., Shen, L., Li, Y., Yang, M., Lin, G.: Open snake model based on global guidance field for embryo vessel location. IET Computer Vision 12(2) (2018) 129–137
-  Cootes, T., Baldock, E., Graham, J.: An introduction to active shape models. Image processing and analysis (2000) 223–248
-  Roy, A.G., Conjeti, S., Karri, S.P.K., Sheet, D., Katouzian, A., Wachinger, C., Navab, N.: Relaynet: retinal layer and fluid segmentation of macular optical coherence tomography using fully convolutional networks. Biomedical optics express 8(8) (2017) 3627–3642
-  Fang, L., Cunefare, D., Wang, C., Guymer, R.H., Li, S., Farsiu, S.: Automatic segmentation of nine retinal layer boundaries in oct images of non-exudative amd patients using deep learning and graph search. Biomedical optics express 8(5) (2017) 2732–2744
-  Cremers, D., Tischhäuser, F., Weickert, J., Schnörr, C.: Diffusion snakes: Introducing statistical shape knowledge into the mumford-shah functional. International journal of computer vision 50(3) (2002) 295–313
-  Duan, J., Pan, Z., Yin, X., Wei, W., Wang, G.: Some fast projection methods based on chan-vese model for image segmentation. EURASIP Journal on Image and Video Processing 2014(1) (2014) 7
-  Noye, B.J., Arnold, R.J.: Accurate finite difference approximations for the neumann condition on a curved boundary. Applied Mathematical Modelling 14(1) (1990) 2–13