Level-Set Method for Image Analysis of Schlemm's Canal and Trabecular Meshwork

Abstract

Purpose: To evaluate different segmentation methods in analyzing Schlemm's canal (SC) and the trabecular meshwork (TM) in ultrasound biomicroscopy (UBM) images.

Methods: Twenty-six healthy volunteers were recruited. The intraocular pressure (IOP) was measured while study subjects blew a trumpet. Images were obtained at different IOPs by 50-MHz UBM. ImageJ software and three segmentation methods—K-means, fuzzy C-means, and level set—were applied to segment the UBM images. The quantitative analysis of the TM-SC region was based on the segmentation results. The relative error and the interclass correlation coefficient (ICC) were used to quantify the accuracy and the repeatability of measurements. Pearson correlation analysis was conducted to evaluate the associations between the IOP and the TM and SC geometric measurements.

Results: A total of 104 UBM images were obtained. Among them, 84 were adequately clear to be segmented. The level-set method results had a higher similarity to ImageJ results than the other two methods. The ICC values of the level-set method were 0.97, 0.95, 0.9, and 0.57, respectively. Pearson correlation coefficients for the IOP to the SC area, SC perimeter, SC length, and TM width were −0.91, −0.72, −0.66, and −0.61 (P < 0.0001), respectively.

Conclusions: The level-set method showed better accuracy than the other two methods. Compared with manual methods, it can achieve similar precision, better repeatability, and greater efficiency. Therefore, the level-set method can be used for reliable UBM image segmentation.

Translational Relevance: The level-set method can be used to analyze TM and SC region in UBM images semiautomatically.

Introduction

Glaucoma is the world's second-leading ocular disease that causes blindness and is the primary cause of irreversible blindness.¹ Elevated intraocular pressure (IOP) is harmful to the optic nerve and can aggravate glaucoma. Therefore, IOP is the most widely used parameter for evaluating and monitoring glaucoma.² IOP is balanced by the production and outflow of the aqueous humor. Most studies on glaucoma pathogenesis have focused on outflow resistance. The trabecular meshwork (TM) and Schlemm's canal (SC) pathway account for 75% to 80% of the whole outflow,³ making it an important area for study. 

Kagemann et al.⁴ showed that an acute IOP elevation can reduce the SC area and alter the TM configuration in human and animal eyes.⁵ Yan et al.⁶ demonstrated that aerobic exercise can cause TM and SC expansion, which lowers IOP. These findings suggest that TM-SC tissue configurations may determine aqueous outflow and IOP regulation. This conclusion also applies to patients with glaucoma. Swain et al.⁷ reported that SC is collapsed in most patients with primary open-angle glaucoma (POAG). Moreover, clinical studies have shown that canaloplasty is an effective and safe procedure to lower IOP in patients with POAG.⁸ Cagini et al.^9,10 determined that canaloplasty is not as successful in eyes that exhibit an irreversible collapse of outflow pathways. These findings suggest that morphologic changes of TM-SC in patients with glaucoma can lower IOP and improve the disease. 

Since TM-SC morphologic changes can influence IOP regulation, identifying a means to monitor the changes in this area has become an important goal in glaucoma research. Imaging systems are emerging as influential tools for evaluating the TM-SC region in vivo. Optical coherence tomography (OCT) and ultrasound biomicroscopy (UBM) are two of the most predominantly used imaging methods in ophthalmology. The former is considered powerful in distinguishing TM on account of its higher resolution and noninvasive nature.^11,12 The latter can be used in almost all kinds of patients, even those who have a cloudy cornea or arcus senilis and cannot be examined with OCT. 

Image segmentation is an important task in medical analysis and a critical step in many clinical applications. In addition, different clustering methods are used for medical image segmentation, such as K-means clustering¹³ and fuzzy C-means clustering (FCM).¹⁴ In previous studies, UBM image segmentation was performed freehand or partly assisted by image analysis software.^6,15–18 Among them, ImageJ software (National Institutes of Health, Bethesda, MD) was the most commonly used. However, manual segmentation is time-consuming and depends on the experience of the technician. Moreover, UBM images are usually corrupted with intensity inhomogeneities, which make TM and SC segmentation an inherently difficult task.¹⁹ Recently, level-set methods have been proposed to deal with images with intensity inhomogeneity, such as the local intensity clustering method^20,21 and the edge-based method.²² These approaches have been successfully applied to segment magnetic resonance images of the breast, X-ray images of bones, ultrasound images of the prostate,²⁰ and infrared breast thermography.²³ 

Inspired by the above research, we investigated the application of different segmentation methods for UBM images. Trumpet playing has the advantages of convenience, feasibility, and simplicity. Therefore, we adopted trumpet playing as an intervention to obtain IOP elevations and fluctuations.²⁴ These segmentation methods were further verified in UBM images obtained under different IOPs. With these segmentation results, quantitative analysis of the TM-SC region was done and compared. 

Methods

Image Acquisition

In this study, healthy volunteers were recruited from the staff of Liaocheng People's Hospital, Shandong Province, China. All subjects underwent an ophthalmologic examination and were verified to have met the following requirements: (1) 20 to 40 years old, (2) IOP between 10 and 21 mm Hg, (3) normal anterior chamber depth and open angle, (4) normal structure of TM examined by a gonioscope, (5) no history of inflammatory eye disease or eye surgery, and (6) no family history of glaucoma. The whole process was approved by the ethics committee of Liaocheng People's Hospital and adhered to the tenets of the Declaration of Helsinki. Written informed consent was obtained from each subject before the study began. Participants were examined in a supine position while blowing the trumpet. They were asked to blow the trumpet slowly for as long as possible.²⁴ One eye of each participant was randomly selected for the UBM examination using a 50-MHz UBM (Suoer SW-3200L; Suowei Co., Tianjin, China). The IOP of the other eye before and during the trumpet blowing was measured using the Icare Pro tonometer (Tiolat Oy, Helsinki, Finland). According to the results of our preexperiment and previous studies,⁶ the SC detection rate at the inferior quadrant is the highest. Thus, all the images were obtained from the inferior quadrant of the eye at four time points: before trumpet blowing, 10 seconds after the start of blowing (IOP increasing period), immediately after blowing cessation (IOP peak time), and 10 seconds after blowing cessation. All examinations were conducted by the same ophthalmologists and under the same illumination conditions using identical equipment. 

Algorithms

Image segmentation is the process of partitioning a digital image into multiple regions. In this study, three algorithms were assessed: K-means,¹³ FCM,¹⁴ and the level-set method.²⁰ 

K-Means Clustering

K-means¹³ is a fast and simple clustering algorithm that is used to classify an image into a specific number of disjointed clusters. The general idea is to identify K centroids, one for each cluster, and then associate each data point to the nearest centroid.²⁵ Let Ω be the image domain, I(x): Ω → R be the observed image, and x, y represent the pixel coordinates. In a previous study,²⁶ segmentation of image I(x) into K clusters was achieved by minimizing the following equation:  

\begin{equation}E = \sum\limits_{j = 1}^K {\sum\limits_{i = 1}^M {{{\left\| {I(x_i^{(j)}) - {c_j}} \right\|}^2}} } ,\end{equation}

(1)

where \(\| {I(x_i^{(j)}) - {c_j}} \|\) is the distance between one of the pixels, \(I(x_i^{(j)})\), in cluster j and its cluster centroid, c_j, and M is the pixel number of the image. 

FCM Algorithm

The FCM algorithm was first suggested by Dunn²⁷ and later improved by Bezdek et al.¹⁴ This algorithm is widely used in data clustering and image segmentation.^28,29 By introducing the possibility of partial memberships to clusters, this algorithm attempts to partition every pixel into a collection of fuzzy cluster centroids by minimizing the following objective function:  

\begin{equation}E = \sum\limits_{j = 1}^K {\sum\limits_{i = 1}^M {u_{_{ij}}^m{{\left\| {I(x_i^{(j)}) - {c_j}} \right\|}^2}} } ,\end{equation}

(2)

where u_ij is the fuzzy membership degree of pixel \(I(x_i^{(j)})\) and cluster centroid c_j, which satisfies u_ij ∈ [0, 1], and \(\sum\limits_{i = 1}^M {{u_{ij}}} = 1,{\rm{ }}j = 1,2, \cdots ,K\). In addition, parameter m(m > 1) is a constant that determines the fuzziness of the resulting partitions. 

Level-Set Method

From the physics of imaging, the observed UBM image I can be modeled as  

\begin{equation}I = BJ + n,\end{equation}

(3)

where J(x) is the real image, B(x) is the bias field that accounts for the intensity inhomogeneity, and n(x) is the noise term.³⁰ The bias field B(x) is assumed to change slowly, and the value B(x) can be considered approximately constant in a neighborhood of O_y = {x||x − y| ≤ ρ}, B(x) ≈ B(y) for x ∈ O_y. Real image J reflects an intrinsic property of the imaging objects, which can be assumed to be a piecewise constant. Moreover, J takes approximately N distinct constant values c₁,c₂,⋅⋅⋅c_N in disjointed regions Ω₁,Ω₂,⋅⋅⋅Ω_N, where \(\Omega = \cup _{i = 1}^N{\Omega _i}\) and Ω_i∩Ω_j = Ø for i ≠ j. Thus, the intensities of points in each subregion Ω_i∩O_y can be approximated as follows:  

\begin{equation}I(x) \approx B(y){c_i} + n(x)\textit{ for }x \in {\Omega _i} \cap {O_y}\end{equation}

(4)

 

Based on the assumption of zero-mean additive Gaussian noise, the intensities in neighborhood O_y can be classified into N distinct clusters with centers m_i ≈ B(y)c_i:  

\begin{equation}I_y^i = \left\{ {I(x):{\rm{ x}} \in {\Omega _i} \cap {O_y}} \right\},{\rm{ }}i{\rm{ = 1,2,}} \cdots {\rm{,}}N{\rm{.}}\end{equation}

(5)

 

The K-means method is used to classify the local intensities in O_y. Then, clustering criterion function ε_y of y in Ω can be written as  

\begin{equation}{\varepsilon _y} = \sum\limits_{i = 1}^N {\int_{{{\Omega _i}}}{{k(y - x){{\left| {I(x) - B(y){c_i}} \right|}^2}dx}}} \end{equation}

(6)

where k(y − x) is the Gaussian kernel function, which is selected as a truncated Gaussian function defined by²⁵ 

\begin{equation}k(y - x) = \left\{ {\begin{array}{@{}l@{\quad}l@{}} {\frac{1}{a}{e^{ - \frac{{{{\left| {y - x} \right|}^2}}}{{2{\sigma ^2}}}}}},&{x \in {O_y},}\\ \\ {0,}&{x \notin {O_y}} \end{array}} \right.\end{equation}

(7)

where a is a normalization constant, such that ∫k(u)du = 1, and σ is the standard deviation of the function. The smaller the value of ε_y, the better the classification of y in Ω. Therefore, the optimal partition of the entire domain Ω can be realized by joint-minimizing ε_y, which can be written as the following local clustering criterion function:  

\begin{eqnarray}\begin{array}{@{}l@{}} \displaystyle\varepsilon = \int\!{{{\varepsilon _y}dy}} = \int\!\!{{\left( {\sum\nolimits_{i = 1}^N {\int_{{{\Omega _i}}}{{k(y - x)}}} {{\left| {I(x) - B(y){c_i}} \right|}^2}dx} \!\right)\!dy}}.\end{array}\!\!\!\nonumber\\\end{eqnarray}

(8)

 

However, function ε is difficult to solve. Therefore, ε is converted into a level-set formulation with several level-set functions. Let φ: Ω → R be a level-set function, and function ε can be written as the function of Φ = (φ₁,φ₂,⋅⋅⋅φ_k), c = (c₁,c₂,⋅⋅⋅c_N) and the bias field b:  

\begin{equation}\varepsilon \left( {\Phi ,c,b} \right) = \int{{\sum\nolimits_{i = 1}^N {\int_{{{\Omega _i}}}{{{e_i}(x)}}} {m_i}(\Phi (x))dx}},\end{equation}

(9)

where e_i(x) = ∫k(y − x)|I(x) − B(y)c_i|²dy and the membership functions m_i=1 for y ∈ Ω_i, and m_i= 0 for y∉Ω_i. For the case of two phases, the membership functions are defined by m₁(φ) = H(φ) and m₂(φ) = 1 − H(φ). The energy function in the two-phase level-set formulation is defined by  

\begin{equation}F\left( {\phi ,c,b} \right) = \varepsilon \left( {\phi ,c,b} \right) + \nu L(\phi ) + \mu {R_p}(\phi )\end{equation}

(10)

where L(φ) and R_p(φ) are the regularization terms. Energy minimization is achieved by an iterative process. By minimizing this energy, the level-set method³⁰ can segment the image and estimate the bias field that can be applied for bias correction. When the above energy function F(φ, c, b) obtains the minimum value or the maximum number of iterations is reached, the iteration is terminated. 

Using ImageJ software, each image was segmented three times by the same ophthalmologist. The average value of the measurements was regarded as the ground truth. The paired t-test was used to compare the mean differences when the measurement data were obtained by ImageJ and the three segmentation methods (P < 0.05 was considered statistically significant).^31,32 The relative error and the interclass correlation coefficient (ICC) were employed to quantify the accuracy and repeatability of the three methods.³³ Because systematic differences are part of the measurement error, a two-way random-effects model was used to calculate the ICC.³⁴ 

Image segmentation and analysis were conducted on images with clear TM-SC structures that were identified by the ophthalmologist. Owing to the low resolutions of the 50-MHz UBM images, it is almost impossible to obtain perfect quality UBM images. If the TM-SC region in the UBM image was primarily not deformed or slightly deformed but not significantly different from the physiologic anatomical structure, then the UBM image was considered of good quality and of use for segmentation. If the TM-SC region in the UBM image was severely deformed, was very different from the physiologic anatomy, or could not be recognized by the ophthalmologist, then the UBM image was considered of poor quality and excluded. Four UBM images of one subject were used to illustrate the segmentation effects of the different methods. The TM-SC regions of all UBM images were extracted using the above three segmentation methods. By measuring the TM-SC area, the reliability and repeatability of the segmentation results were quantified. Furthermore, the correlations between the measurements and the IOP were obtained. All experiments were carried out in MATLAB (MathWorks, Natick, MA) on a PC with a 3.6-GHz Intel core processor and 8 GB of memory. 

Results

A total of 26 volunteers completed the whole experimental process, including 10 males and 16 females. The average age was 34.53 years. Four UBM images were collected from each subject at different IOPs, and 104 images were collected in total. The average IOP before trumpet blowing and immediately after was 17.5 mm Hg and 28.8 mm Hg, respectively. The UBM image and its region of interest (ROI) are shown in Figure 1. The size of the UBM image was 1024 × 655 pixels, and the size of the ROI was 150 × 100 pixels. 

Figure 1.

 

UBM image and ROI.

Image Segmentation with Different Methods

To illustrate the segmentation effects of the different methods, four UBM images (Fig. 2) from one subject were segmented. ImageJ, K-means, FCM, and level set were used to obtain the boundary curves of the different gray regions in the UBM images. The segmentation results are shown in Figures 3 to 6. 

Figure 2.

 

Preprocessed UBM images from one subject.

Figure 3.

 

UBM image segmentation result of the ROI at the IOP of 15 mm Hg using (a) ImageJ, (b) K-means, (c) FCM, and (d) level set.

When the IOP is 15 mm Hg, the SC region segmentation result in Figure 3d is similar to that in Figure 3a. The TM and SC in Figures 3b and 3c are blended, resulting in an unrecognizable SC boundary, whereas the TM boundary in Figure 3d is easy to identify. The TM boundary in Figure 3b is discontinuous, and the boundary in Figure 3c is overlapped and unclear. In terms of the TM and SC segmentation effects in the UBM image, the level-set method is thus better than the other two methods. 

Figure 4 shows the UBM image segmentation results of the ROI at the pressure of 22 mm Hg. It is observed that the SC region in Figure 4d is more consistent with Figure 4a. The SC boundaries in Figures 4b and 4c are rougher than those in Figure 4a. There are some disturbances on the TM boundary in Figure 4d. However, there are two TM boundaries in Figure 4b, and it is difficult to choose the correct one. The TM boundary in Figure 4c is oversegmented. Therefore, we can conclude that the level-set method can obtain results similar to those of ImageJ when the IOP is elevated. 

Figure 4.

 

UBM image segmentation result of the ROI at the IOP of 22 mm Hg using (a) ImageJ, (b) K-means, (c) FCM, and (d) level set.

The UBM image segmentation results at the pressure of 27 mm Hg are shown in Figure 5. The SC boundaries in Figures 5b and 5c obtained by the FCM and K-means methods are discontinuous. Some interference occurs between the TM and SC, which makes identifying the TM area difficult. The SC region segmented using the level-set method has a clear outline, and the TM region boundary can be easily identified. It is obvious that the level-set method produces more accurate segmentation results than the other two methods. 

Figure 5.

 

UBM image segmentation result of the ROI at the IOP of 27 mm Hg with (a) ImageJ, (b) K-means, (c) FCM, and (d) level set.

To verify the universality of the segmentation methods, we tested the three segmentation methods using the UBM image at the pressure of 33 mm Hg. From Figure 6, we can observe that the geometric deformation of SC becomes more obvious as the IOP increases. There are some small differences in the geometry of the TM-SC area in Figures 6a and 6d. The TM and SC in Figures 6b and 6c are connected on account of speckle noise, which also leads to the interruption of the TM boundary. The TM-SC region extracted by the level-set method is more similar to that of the ImageJ result than the other two methods. 

Figure 6.

 

UBM image segmentation result of the ROI at the IOP of 33 mm Hg using (a) ImageJ, (b) K-means, (c) FCM, and (d) level set.

Quantitative Analysis of Segmentation Results

Eighty-four images that were adequately clear were selected from a total of 104 images. After the UBM images were segmented, the TM-SC region was extracted and measured. The measurement was performed as follows. The Canny edge detection algorithm was applied to convert the segmentation result to a binary form. The number of pixels inside the SC region was used as the SC area, and the number of pixels on the boundary was regarded as the SC perimeter. The average of the three maximum numbers of pixels per row in the SC region was used as the SC length. The Sobel edge detection algorithm was employed to binarize the boundary curve of the TM-SC region. The average of the three maximum numbers of pixels between the TM and SC per column was used as the TM width. 

Table shows the measurement data of the four methods as the mean ± SD. There were no statistically significant differences between the measurement data obtained by the level-set method and ImageJ (P = 0.663, P = 0.071, P = 0.755, and P = 0.117 for SC area, SC perimeter, SC length, and TM width, respectively). There were no statistically significant differences in the SC area and SC perimeter measurements between the K-means method and by ImageJ (P = 0.103 and P = 0.901 for SC area and SC perimeter, respectively), while the differences for the SC length and TM width between the two methods were statistically significant (P < 0.001 for SC length and TM width). There was no statistically significant difference between the SC area measured by the FCM method and the corresponding result measured by ImageJ (P = 0.662), while the differences in the SC perimeter, SC length, and TM width between the two methods were statistically significant (P = 0.032, P < 0.001, and P < 0.001 for SC perimeter, SC length, and TM width, respectively). 

Table.

 

Measurement Data (Mean ± Standard Deviation) for Segmentation Results of the Four Methods (Pixels, n = 84)

 

As can be seen from the above results, the level-set method showed better similarity to the ImageJ method than the FCM and K-means methods. Among the four parameters measured by the ImageJ and three segmentation methods, the SC area was the most consistent parameter. 

Figure 7 presents schematic diagrams of the relative errors and ICC values of the three segmentation methods. It is evident that the relative errors of the level-set method are less than 0.08, whereas the other two methods have large relative errors. The ICC values of the level-set method are 0.97, 0.95, 0.9, and 0.57, respectively, whereas the corresponding ICC values of the other two methods are less than 0.4. From these results, we can conclude that the measurements of the level-set method have a higher reliability and better repeatability than the other two methods. 

Figure 7.

 

Schematic diagrams of relative errors and ICC values of the three segmentation methods.

Correlation Analysis Between IOP and Measurements

Owing to the poor performance of the K-means and FCM methods for segmenting the TM-SC region, only the level-set method was used to perform the correlation analysis. The correlations among the SC area, SC perimeter, SC length, TM width, and IOP were analyzed using the Pearson correlation coefficient and linear regression analysis. The results are shown in Figure 8. It can be observed that, as the IOP increases, the SC area, perimeter, and length tend to decrease. The TM width likewise decreases. Pearson correlation coefficients for IOP to the SC area, SC perimeter, SC length, and TM width are −0.91, −0.72, −0.66, and −0.61, respectively. Thus, a negative correlation relationship between the IOP and the geometrical measurement of TM and SC can be inferred. 

Figure 8.

 

IOP versus measurements of the (a) SC area, (b) SC perimeter, (c) SC length, and (d) TM width. For reference, regression equations, linear regression lines (red), Pearson correlation coefficients (R), and P values are provided.

Discussion

Image segmentation plays an important role in many medical imaging applications. In previous studies, the TM-SC region was manually obtained from UBM images.^6,15–18 However, manual outlining is a time-consuming and tedious task. In this article, we showed that the level-set method can be used to extract the TM-SC region from the UBM image and useful features can be obtained from the segmentation. Compared with manual segmentation, the level-set method produced similar segmentation results while providing better repeatability and efficiency. In addition, the level-set method had higher accuracy than the classical FCM and K-means methods. 

The negative correlation between IOP and the measurements of TM and SC inferred from our results were similar to those of previous studies.^4,5 The reason for TM-SC region collapse may be the compression force caused by acute IOP elevation to the elastic structure. Assuming that accurate measurements of the TM and SC response to IOP fluctuation in patients in vivo are realized, mathematical models can be used to calculate the TM stiffness, which has recently been shown to be associated with resistance to outflow.^35–37 Among the four measurement indicators in this study, the Pearson correlation coefficient for IOP to the SC area is the highest. Thus, we may speculate that the SC area can be used as a sensitive indicator for measuring the TM-SC region. 

However, our study had the following limitations. First, the resolution of the 50-MHz UBM was not high; therefore, approximately 20% of the images could not be used. Second, since the UBM images were continuously obtained while the measurement of IOP was discontinuous, the mismatch between the two may have affected the correlation analysis results. Third, the subjects recruited for this experiment were young healthy adults. We are therefore unsure whether this method can be used with elder subjects and patients with glaucoma. Therefore, the effectiveness of the level-set method for UBM image segmentation of different subject groups will be our future work. 

In conclusion, changes in the TM-SC region can be detected by UBM and extracted by image segmentation methods. The level-set method can accurately and efficiently segment UBM images of TM and SC. Therefore, the level-set method is an effective technique for UBM image segmentation. 

Acknowledgments

Supported by National Key Research and Development Program for Stem Cell and Translational Research Foundation of China (2018YFA0109500). 

Disclosure: X. Wang, None; Y. Zhai, None; X. Liu, None; W. Zhu, None; J. Gao, None