Scholarly article on topic 'Segmentation of Choroidal Boundary in Enhanced Depth Imaging OCTs Using a Multiresolution Texture Based Modeling in Graph Cuts'

Segmentation of Choroidal Boundary in Enhanced Depth Imaging OCTs Using a Multiresolution Texture Based Modeling in Graph Cuts Academic research paper on "Medical engineering"

CC BY
0
0
Share paper
OECD Field of science
Keywords
{""}

Academic research paper on topic "Segmentation of Choroidal Boundary in Enhanced Depth Imaging OCTs Using a Multiresolution Texture Based Modeling in Graph Cuts"

Hindawi Publishing Corporation

Computational and Mathematical Methods in Medicine Volume 2014, Article ID 479268, 9 pages http://dx.doi.org/10.1155/2014/479268

Research Article

Segmentation of Choroidal Boundary in Enhanced Depth Imaging OCTs Using a Multiresolution Texture Based Modeling in Graph Cuts

Hajar Danesh,1 Raheleh Kafieh,1 Hossein Rabbani,1 and Fedra Hajizadeh2

1 Department of Biomedical Engineering, Medical Image and Signal Processing Research Center, Isfahan University of Medical Sciences, Isfahan 81745, Iran 2Noor Ophthalmology Research Center, Tehran 1968653111, Iran

Correspondence should be addressed to Hossein Rabbani; h_rabbani@med.mui.ac.ir

Received 20 August 2013; Revised 30 November 2013; Accepted 19 December 2013; Published 11 February 2014 Academic Editor: William Crum

Copyright © 2014 Hajar Danesh et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

The introduction of enhanced depth imaging optical coherence tomography (EDI-OCT) has provided the advantage of in vivo cross-sectional imaging of the choroid, similar to the retina, with standard commercially available spectral domain (SD) OCT machines. A texture-based algorithm is introduced in this paper for fully automatic segmentation of choroidal images obtained from an EDI system of Heidelberg 3D OCT Spectralis. Dynamic programming is utilized to determine the location of the retinal pigment epithelium (RPE). Bruch's membrane (BM) (the blood-retina barrier which separates the RPE cells of the retina from the choroid) can be segmented by searching for the pixels with the biggest gradient value below the RPE. Furthermore, a novel metho d is proposed to segment the choroid-sclera interface (CSI), which employs the wavelet based features to construct a Gaussian mixture model (GMM). The model is then used in a graph cut for segmentation of the choroidal boundary. The proposed algorithm is tested on 100 EDI OCTs and is compared with manual segmentation. The results showed an unsigned error of 2.48 ± 0.32 pixels for BM extraction and 9.79 ± 3.29 pixels for choroid detection. It implies significant improvement of the proposed method over other approaches like fc-means and graph cut methods.

1. Introduction

Optical coherence tomography (OCT) imaging technique was introduced by Huang et al. in 1991 [1]. This technique is employed for its ability in taking cross-sectional images from microscopic structure of the living tissues. The device has a resolution on the scale of micrometers. In this technique, a number of A-scans or linear scans create B-scan or cross-sectional images [2]. OCT images consist of a great amount of data; therefore, nonautomated and visual analysis of such a great amount of data would be troublesome for the ophthalmologist. The main goal of automatic segmentation is to assist the ophthalmologist in diagnosis and monitoring of eye diseases.

Choroid is one of the structural layers located between the sclera and the retina. This pigmented layer contains many capillaries that supply feeding of the iris and retinal light

receptor cells [3]. Many diseases such as polypoidal choroidal vasculopathy and choroidal tumors cause changes in the structure of this layer [4-8]; therefore, segmentation of this layer has great importance for ophthalmologists. In retinal imaging using conventional OCT, wavelength of the light source is around 800 nanometers, which is not appropriate for imaging of choroid layer, due to signal transmission difficulty through retinal pigment epithelium (RPE) layer and increased depth of imaging [9]. However, increased pixel density and high signal-to-noise ratio in SD-OCT (spectral-domain OCT) in comparison to TD-OCT (time-domain OCT) makes choroidal imaging possible. In SD-OCT, the structures that are closer to "zero delay line" have higher signals than those that are farther away. In conventional OCT, the "zero delay line" is near the inner surface of retina; however, in EDI-OCT, it is placed near the outer retina and choroid, which is the key for EDI imaging. In order to

increase the quality of EDI-OCT and to reduce "speckles," a high number of images (usually 100 images) are averaged through a software that provides high-quality images with smoother border [9]. An example of this imaging method is compared with conventional OCT in Figure 1. However, another solution for choroidal imaging is a higher wavelength of approximately 1060 nanometer [10,11].

In several studies, EDI-OCT has been used to measure the thickness of the choroid, finding its relation with diseases, and monitoring the treatment process [12-14]. In most of these studies, measurement of choroidal thickness is usually accomplished by manual labeling which is a time-consuming and tedious process. This problem is much more complicated when the number of images is numerous. Therefore, the need for development of an automatic segmentation algorithm on EDI-OCT arises.

A limited number of studies have already been conducted about automatic segmentation of these images. Kajic et al. [15-17] proposed a two-stage statistical model to detect the choroid boundaries in the 1060 nm OCT images in healthy and pathological eyes. This model needed extensive training and the mean error is 13%. Tian et al. [18, 19] found the choroidal boundary by finding the shortest path of the graph formed by valley pixels using dynamic programming (DP). The average of Dice's coefficient on 45 EDI-OCT images was 90.5%. Lu et al. [20] proposed a technique to segment the inner boundary of the choroid using a two-stage fast active contour model. Then a real-time human-supervised automated segmentation on the outer boundary of the choroid was applied. The reported Dice similarity coefficient value on 30 images captured from patients diagnosed with diabetes was 92.7%.

Many algorithms based on wavelet and graph cut are already used in segmentation of retinal layers [21]. However, due to heterogeneity in the choroid layers, such methods cannot be useful in choroid segmentation. However, because of distinct difference between texture of choroid and other retinal layers, an algorithm based on texture classification can be effective. With this idea, we use a combination of graph-based methods and wavelet-domain features for choroidal segmentation. A description ofour method including detection of Bruch's membrane (BM) and segmentation of choroidal-scleral interface (CSI) is provided in Section 2 and the results are presented in Section 3. The results show

the improved robustness of the method compared to other available algorithms. Finally the conclusion and suggestions for future works are presented in Section 4.

2. Material and Methods

The proposed choroid segmentation method is tested on 100 randomly selected two-dimensional EDI-OCT images obtained from 10 eyes of 6 normal subjects by an EDI system of multimodality diagnostic imaging (wavelength: 870 nm; scan pattern: enhanced depth imaging; Spectralis HRA + OCT; Heidelberg Engineering, Heidelberg, Germany) [22]. The study protocol was reviewed in advance by the review board of Noor Ophthalmology Research Center. Each participant was informed of the purpose of the study and provided a written consent to participate. For evaluation and ruling out any ocular pathology all patients underwent thorough ophthalmic examinations including refraction, visual acuity, slit lamp biomicroscopic examination, IOP measurement by Goldmann applanation tonometer, and examination of the fundus with plus 90-D lens. Each dataset consisted of a SLO image and limited number of OCT scans with size of 496 x 768 (namely, for a data with 31 selected OCT slices, the whole data size would be 496 x 768 x 31). The pixel resolutions were 3.9 ^m/pixel in axial direction and 14 ^m/pixel in transversal direction.

2.1. Overview of Algorithm. In the first step, dynamic programming (DP) was applied for boundary detection of RPE layer. The BM is the blood-retina barrier that separates the RPE cells ofthe retina from the choroid and can be segmented by searching for the pixels with the biggest gradient value below the RPE. In this stage, we eliminated the pixels above this boundary form our calculations by making them equal to zero and focused on segmentation of the lower boundary of the choroid.

In the next step, discrete wavelet transform (DWT)— using filter bank and descriptors—was employed to extract appropriate features from EDI-OCT images. The extracted features from training images were used to construct a Gaussian mixture model (GMM). It is assumed that each extracted descriptor using DWT has a multivariate normal

Finding the biggest gradient below the RPE

Extraction of wavelet domain coefficients

Manual segmentation

Training images

BM segmentation

Texture extraction

Gaussian mixture model

Graph cut segmentation

Choroidal boundary segmentation

Figure 2: Choroid segmentation algorithm overview.

c(x™)

Figure 3: Dynamic programming: (a) one step of cost calculation; (b) graph layers and nodes notation. The highlighted pixels show the neighbors.

probability density function (pdf) and the number of Gaussian distributions (K) is chosen to be 3, which represent the area above the BM, choroid, and the area beneath the CSI. Finally, the relevance of each data to the model is used in creating the graph for final CSI segmentation. A block diagram of the algorithm could be found in Figure 2.

2.2. BM Detection. We applied boundary detection algorithm using DP to detect boundary of RPE layer which is the brightest part in EDI OCT images. DP is a method based on optimality [23], which seeks for the best functions in which the variables are not simultaneously in conjunction with each other. In simple problem of boundary tracking, the goal is to find the best path (least expensive) between one of the possible starting points and one of the possible ending points.

In this step, upper bound of choroid (lower bound of RPE) can be segmented.

For boundary detection using DP, we are looking for a dark boundary transmitted from one edge to another edge in an image (e.g., from column 1 to M - 1 in Figure 3). DP algorithm associates a cost to each node (brightness of each pixel) and the partial paths are set according to the cost (brightness) of the next neighboring node. The algorithm starts through the columns (column 1 to M - 1) and for each node (from 1 to n), the best path with optimized (lowest) cost is selected among three neighbor pixels of the neighboring column (i = -1,0,1 in (1)). Therefore, a cumulative cost (C(x™+1)) can be made for a pixel in m + 1th column and in the fcth node, according to the cost of previously selected nodes (1).

(a) (b)

Figure 4: (a) EDI-OCT image of the human eye. (b) Bruch's membrane.

The following algorithm describes the DP method [24].

(1) Specify initial costs C(xj) of all nodes in the first graph layer, i = 1,...,n and partial path costs gm(i, k), m = 1,...,M- 1.

(2) Repeat step 3 for all m= 1,...,M-1.

(3) Repeat step 4 for all nodes k = 1,... ,n in the graph layer m.

(4) Let

cixk+1 ) = j=im1in>1 {c(xk+) + gm (i,k)). (i)

Set pointer from node back to node x™*,where * denotes the optimal predecessor. Optimality is defined as minimization of the cost function.

(5) Find an optimal node x1^* in the last graph layer M and obtain an optimal path bybacktracking through the points from xk * to xk*.

The positions of a sample node C(x™+1) and its neighbors in the previous (mth) graph layer are shown in Figures 3(a) and 3(b). The highlighted pixels around C(x™+1) indicate i = -1,0,1 in (1).

It should also be noted that the DP boundary detection is conventionally looking for the darkest boundary in an image and for our application (detection of the brightest boundary) we should inverse the brightness values of the image. After segmentation of RPE layer, BM can be located by searching for the pixels with the biggest gradient below the RPE. Figure 4 shows the result of BM segmentation for a sample EDI-OCT image. From calculative point of view, DP is efficient and flexible and can be considered as a powerful tool even in presence of noise.

2.3. CSI Detection. A novel combination of algorithms is proposed for discrimination of lower boundary of choroid (CSI). Among many possible texture descriptors used for segmentation of images with several textures [24], we show the ability of wavelet texture descriptors in this work.

The algorithm can then be stated in three steps:

(1) calculate wavelet descriptors for each pixel (Section 2.3.1);

(2) create a model of the classes using training images. (Section 2.3.2);

(3) take a new image and segment it using the learned model. Use graph cut segmentation to obtain spatially coherent segmentation (Section 2.3.3).

It is assumed that extracted descriptors using DWT have a multivariate normal pdf. This assumption can be justified by the central limit theorem (CLT). In its classical form, the CLT states that the mean of n independent identically distributed random variables, with finite mean and variance, has a limiting distribution for n ^ >x> that is Gaussian. As an approximation, for a finite number of observations (n ^ >x), it provides a reasonable approximation only close to the peak ofthe normal distribution rather than stretching into the tails (convergence to the limit).

In our datasets of Heidelberg 3D OCT Spectralis, each image is calculated from the mean value of 100 images taken at the same time. The mean process dramatically improves the quality of the information in the images, without introducing an alteration to the original image or adding any distracting electronic noise [22]. Therefore, according to the "convergence to the limit" theorem near the peak areas (needed for GMM assumption), and in accord with linear characteristics of the wavelet transform, we can assume that each extracted descriptor using DWT has a multivariate normal pdf.

2.3.1. Wavelet Descriptors. DWT using filter bank and descriptors is employed to extract features for each pixel of the image. In other applications of wavelets transform for texture analysis, most prevalent features are wavelet energy signatures and their second-order statistics [25, 26]. In this case, Haar wavelet filters (H(z) = (1 + z)/2 and G(z) = (z - 1)/2) are utilized for ease of use. Filter bank stages are repeated (multiresolution approach) for four levels. Since we are using wavelet frame instead of the standard DWT, the number of levels can optionally grow. But, with trial and error, we found that, for levels higher than 4, the results did not improve and we only faced with more complexity. The energy of each high-pass filter at different stages of subbands and the energy of the last stage of the low-pass filter are extracted as image features (5 features). Extracted features and a training image with a known segmentation (mask) are then used to construct GMM for the image.

2.3.2. Gaussian Mixture Model. GMM is a parametric pdf represented as a weighted sum of Gaussian component distributions. The idea for clustering with GMMs is the same as for fc-means. In conventional GMM, we make an initial guess for mean and covariance rk and then try to iteratively refine them. The number of Gaussian distributions (K) is chosen to be 3, which represent the area above the BM, choroid, and the area beneath CSI.

Assume that n-dimensional data constituted K Gaussian distribution, N2,..., Nfc, in which Rk with the mean and covariance rk is considered as follows:

Considering the weights of each relevant K as nk,

Œfc=i nk = 1)> the final pdf is defined as follows:

p(xj) = £nkp(xj I Nfc),

in which

p(xj I*k)

(2 7t)

n/2\r r

exP (-\(xj - ^k)Trk1 (xj - Fk)) •

Since we are working on wavelet features described in Section 2.3.1, n should be set to 5. We now assume that actual values of the observed EDI-OCT images X = {X1, X2,..., Xt,... XM] and value of pixels in each image

are independent. We also

suppose that EDI images have equal size of AxB. With defining N = MxAxB,we will have the following equations:

P(X) = UP(xj)

— n£

j=1k=1 (2^)n/2\r

X eXP (-2(Xj - ^k)Trk1 (xj - Fk)) :

where M represents the number of EDI-OCT images used for training (10 for our case). For each training image, a mask was made manually which classifies each image into three labels: the area above the BM, choroid, and the area beneath CSI.

In conventional GMM, the parameters are estimated using iterative EM algorithm to fit a mixture model to a set of available data [24,27].

EM algorithm is rapidly repeated and at each stage, Gaussian expectation is estimated for each data sample, and then Gaussian element estimation (maximization) is altered. If ^k and rk are in hand as an estimation, we can calculate the fcth Gaussian probability for Xy

™k P(xj 1 Nk) a= 1 ni P(xj 1

in which Xj (probability ratio for Nk) is granted based on sum of Xj (regardless of Gaussian production) balanced with nt. Therefore, the following equation is defined:

1 n — Ñ£pjk

In this equation, mean pjk is calculated in data series. Similarly, one can estimate corrected values of and rk

new _ ¿j=1 PjkXj Vk — VN ^ ,

Zj=1 Pjk

rnew — IN=1 ^jk (Xj -Fr)(*j -FfcneW)T Zj=1 Pjk

With a slight modification to conventional GMM, which uses EM algorithm to calculate the parameters of each Gaussian function (including mean and covariance in (8)), we use the training step for finding the parameters of GMM. Namely, we use training images with known segmentation (mask) and calculate the mean and covariance of each section (for the training data). These parameters are used for all images in our database (i.e., there is no need to use EM algorithm for each image separately, which speeds up the algorithm). Then we only calculate the responsibility factor by (6) to construct a learned model to be fed to graph cut segmentation.

When a new test image is considered for segmentation, we calculate wavelet descriptors for each pixel (and accordingly, 5-dimensional vectors of Xj and The probability (p) of a pixel belonging to a particular class of the learned model can be obtained, afterward. The value of probability (p) is then used in the construction of a graph cut segmentation.

2.3.3. Graph Cut Segmentation. Graph cut segmentation is constructed based on a learned model to obtain spatially coherent segmentation.

The direct use of optimization algorithms of minimumcut/maximum flow for graph partitioning was first introduced by Greig et al. [28] in image processing of binary images. Using graph optimization algorithms, a powerful method of optimal bordering and region classification in N-dimensional image data was proposed [29, 30].

This method starts by one or more points representing "object" and one or more points representing "background," determined using interactive or automatic identification. The overall shape of cost function C is represented as follows [31]:

C (/) — cdata (/) + csmooth (/).

To minimize C(f), a special class of arc-weighted graphs Gst = (V U {terminal nodes}, E) is employed. In addition to the set of nodes V corresponding to pixels of the image I (Figure 5(a)), terminal nodes (shown by 1, 2, and 3 in (6) Figure 5(b)) are also added to Gst. These terminals are hard-linked with the segmentation seed points (bold links in

Figure 5: Graph cut classification—example of simple classification. (a) Image with 3 classes. (b) Related graph and terminal nodes. Green lines show «-links, while red, yellow, and blue curves represent the i-links.

Figure 5(b)) and represent the segmentation labels (1, 2, and

The arcs E in Gst can be classified into two categories: n-links and i-links. The n-links connect pairs of neighboring pixels whose costs are derived from the smoothness term csmooth(f). The i-links connect pixels whose costs are derived from the data term cdata(f). Green lines in Figure 5(b) show n-links, while red, yellow, and blue curves represent the t-links. A s - t cut in Gst is a set of arcs whose removal partitions the nodes into three disjoint subsets (1, 2, and 3 in Figure 5). The cost of a cut is the total cost of arcs in the cut, and a minimum s - t cut is a cut whose cost is minimal. The minimum s - t cut problem or its dual (the maximum flow problem) can be solved by various algorithms. In maximum flow algorithm, maximum amount of water to sink is delivered using directed arc graphs and the amount of water flow through each separate arc is determined using arc capacity or cost. The greatest amount of maximum flow from s to t saturates a set of graph arcs. These saturated arcs divide nodes into two separate sections of S and T, related to minimum cuts [24, 32, 33].

Let each image pixel ik take a binary label Lk e {1,2,3}. The labeling vector L = (L 1,L 2,...,L |7|) defines the resulting binary segmentation.

In this paper, Cdata(f) is the distance between each image pixel and initial class made by GMM in Section 2.3.2. For a test image, wavelet descriptors are calculated for each pixel and probability (p) of a pixel belonging to a particular class of the learned model is obtained. Then Cdata(f) is created by

Cdata (f) = - log ip + eps).

Csmooth(f) is a matrix of costs which is related to adjacent pixel values. In this work, we assigned a square matrix with fixed values (r = 4) and with zeros on diagonal as Csmooth (f) (9). The increase in fixed term (r) would enhance adjacency

smooth

constraint constant and therefore would result in a finer separation:

0r r ...... r

r 0 r r '• .

r r 0 r r .

. r r 0 r r

. '• r r 0 r r ...... r r 0

3. Results

One hundred two-dimensional EDI-OCT images were obtained from Heidelberg 3D OCT Spectralis and were used in statistical analysis to produce the results. For our 100 two-dimensional dataset, we chose 10 images to train the learned model. Actually, we can choose any 10 or less images for training and calculate means and covariances and weights of GMM and there would be no considerable difference between these parameters. This step will expedite our algorithm instead of using EM algorithm for each image. The performance of the method is reported based on its ability in correct segmentation of the test images. For evaluation of the proposed method, the manual segmentation of two observers was used as the gold standard.

For validation purpose, the mean signed and unsigned border positioning errors for each border were compared with other algorithms such as graph cut, fc-means, and DP and the results are presented in Tables 1 and 2 for each boundary. We implemented each of these algorithms and tested them on our 90 two-dimensional test set. In fc-means algorithm, we selected k = 3 and applied fc-means algorithm on the image, directly. In graph cut method, we used the result of fc-means algorithm to create initial prototype for each class and Cdata(f) was calculated by the distance between each image pixel to initial class made by fc-means algorithm.

Table 1: Summary of mean signed border positioning errors (mean ± std).

Avg. obs. versus our alg. Avg. obs. versus graph cut alg. Avg. obs. versus fc-means alg. Avg. obs. versus DP alg. Obs. 1 versus obs. 2

BM 0.44 ± 1.18 4.30 ± 0.72 7.73 ±3.35 0.44 ± 1.18 1.12 ± 1.29

Choroid 5.77 ± 2.77 57.69 ± 6.93 -31.25 ± 11.51 9.65 ± 5.41 7.65 ± 2.63

Table 2: Summary of mean unsigned border positioning errors (mean ± std).

Avg. obs. versus our alg. Avg. obs. versus graph cut alg. Avg. obs. versus fc-means alg. Avg. obs. versus DP alg. Obs. 1 versus obs. 2

BM 2.48 ± 0.32 4.51 ± 0.71 7.73 ±3.35 2.48 ± 0.32 2.64 ± 0.98

Choroid 9.79 ± 3.29 65.69 ± 7.53 33.73 ± 12.23 12.48 ± 5.41 8.78 ± 2.38

Table 3: Improvement of the proposed method compared with other algorithms.

P value, our alg. versus graph cut alg. P value, our alg. versus fc-means alg. P value, our alg. versus DP :

BM «0.001 «0.001 —

Choroid «0.001 «0.001 «0.01

Csmooth(/) was calculated in the same method described in Section 2.3.3 by small change of selecting r = 3 which could give the best results. For DP method, the first step is similar to the proposed method in this paper (we applied boundary detection algorithm using dynamic programming to detect RPE layer boundary and eliminated the pixels above this boundary by making them equal to zero). Then we eliminated a region beneath the RPE (7 pixels below) and applied DP to search the image for another time. The results were based on the whitest route available after elimination of RPE layer.

According to Tables 1 and 2, the signed border positioning errors were 0.44 ± 1.18 pixels for BM extraction and 5.77 ± 2.77 pixels for choroid segmentation, and the unsigned border positioning errors were 2.48 ± 0.32 pixels for BM extraction and 9.79 ± 3.29 pixels for choroid segmentation, respectively. The errors between the proposed algorithm and the reference standard were similar to those computed between the observers. For example, the overall observer error was 2.64 ± 0.98 pixels and 8.78 ± 2.38 pixels for BM and CSI, respectively, which is comparable to the results of the algorithm. The border positioning errors of the proposed method showed significant improvement over other algorithms, compared in both of the tables.

To show the statistically significant improvement of the proposed method over the compared algorithms, Table 3 shows the obtained P values. The values confirm that our results have a significant improvement. For instance, the algorithm's overall P value against fc-means and graph cut was less than 0.001 and against dynamic programming, it was less than 0.01.

The pixel resolution of our datasets in axial direction was 3.9 ^m/pixel. Therefore, mean signed positioning errors for localization of BM and choroid are 1.71 and 22.50 ^m, respectively. In accord with repeatability measurements for choroidal thickness of EDI-OCT images [12], it can be

concluded that the positioning error has acceptable accuracy. Figure 6 demonstrates two samples showing the performance of the proposed method. Furthermore, Figure 7 shows the results of segmentation using the proposed algorithm in comparison with other methods.

The computational complexity of the proposed algorithm is around 13 seconds in the case of using the saved parameters for constructed models; otherwise, if the computation consists of the model construction for one image, it takes around 21 seconds on a PC with Microsoft Windows XP x32 edition, Intel core 2 Duo CPU at 3.00 GHz, 4 GB RAM.

4. Conclusion

In this paper, a new method for choroid segmentation in EDI-OCT is introduced. In the first step, RPE boundary was segmented using DP method and the area above RPE was eliminated from next steps. Then, the wavelet descriptors were calculated for each pixel of the training images and assuming that the pdf of these descriptors are normal, their parameters were calculated to construct a model.

When a new test image is considered for segmentation, wavelet descriptors were calculated for each pixel and probability (_p) of a pixel belonging to a particular class of the learned model was obtained. The value of probability (_p) was then used in the construction of a graph cut segmentation.

The main limitation of this method is the need for manual segmentation to construct the model. Despite the fact that only a few training images are sufficient to produce good results, the manual labeling may be troublesome and replacing this step with an automatic method can be studied in future works. For example, in this work we choose several images as known data and extracted their mixture model's parameters for using as parameters of all new data

in the database. For each database, arbitrary sample data can be used for constructing this mixture model (with known parameters) and then used for all data in database. However, the automatic version of this method can also been considered without using a predefined mixture model and by using EM algorithm for each image.

The new method has better accuracy than previous methods and the algorithm is tested on a larger dataset compared to older algorithms. The performance of the method is also fast and the implementation is relatively simple in comparison to more complicated algorithms.

As a future work, the proposed method should be tested on EDI images taken from other modalities to prove its robustness to the imaging technique. Furthermore, a 3D choroidal thickness map can be constructed using 3D OCT dataset which can assist the ophthalmologist in the

diagnosis of the choroidal diseases. We also like to work on segmentation of the blood vessels in choroidal layer and produce a 3D vascular structure to give more information about distribution of the vessels.

Conflict of Interests

The authors disclose that there is no conflict of interests including any financial and personal relationships with other people or organizations that could inappropriately influence (bias) their work.

References

[1] D. Huang, E. A. Swanson, C. P. Lin et al., "Optical coherence tomography," Science, vol. 254, no. 5035, pp. 1178-1181,1991.

[2] D. Huang, E. A. Swanson, C. P. Lin et al., Optical Coherence Tomography, Massachusetts Institute of Technology, Whitaker College of Health Sciences and Technology, 1993.

[3] G. A. Cioffi, E. Granstam, and A. Alm, "Ocular circulation," in Adler's Physiology of the Eye, Mosby, St. Louis, Mo, USA, 2003.

[4] G. Barteselli, J. Chhablani, S. El-Emam et al., "Choroidal volume variations with age, axial length, and sex in healthy subjects: a three-dimensional analysis," Ophthalmology, vol. 119, no. 12, pp. 2572-2578, 2012.

[5] S. E. Chung, S. W. Kang, J. H. Lee, and Y. T. Kim, "Choroidal thickness in polypoidal choroidal vasculopathy and exudative age-related macular degeneration," Ophthalmology, vol. 118, no. 5, pp. 840-845, 2011.

[6] P. Jirarattanasopa, S. Ooto, A. Tsujikawa et al., "Assessment of macular choroidal thickness by optical coherence tomography and angiographic changes in central serous chorioretinopathy," Ophthalmology, vol. 119, no. 8, pp. 1666-1678, 2012.

[7] J. R. Ehrlich, J. Peterson, G. Parlitsis, K. Y. Kay, S. Kiss, and N. M. Radcliffe, "Peripapillary choroidal thickness in glaucoma measured with optical coherence tomography," Experimental Eye Research, vol. 92, no. 3, pp. 189-194, 2011.

[8] J. Yeoh, W. Rahman, F. Chen et al., "Choroidal imaging in inherited retinal disease using the technique of enhanced depth imaging optical coherence tomography," Graefe's Archive for Clinical and Experimental Ophthalmology, vol. 248, no. 12, pp. 1719-1728, 2010.

[9] R. F. Spaide, H. Koizumi, and M. C. Pozonni, "Enhanced depth imaging spectral-domain optical coherence tomography," The American Journal of Ophthalmology, vol. 146, no. 4, pp. 496500, 2008.

[10] R. F. Spaide, "Enhanced depth imaging optical coherence tomography of retinal pigment epithelial detachment in age-related macular degeneration," The American Journal of Ophthalmology, vol. 147, no. 4, pp. 644-652, 2009.

[11] B. PovaZay, B. Hermann, A. Unterhuber et al., "Three-dimensional optical coherence tomography at 1050 nm versus 800 mn in retinal pathologies: enhanced performance and choroidal penetration in cataract patients," Journal ofBiomed-ical Optics, vol. 12, no. 4, Article ID 041211, 2007.

[12] W. Rahman, F. K. Chen, J. Yeoh, P. Patel, A. Tufail, and L. Da Cruz, "Repeatability of manual subfoveal choroidal thickness measurements in healthy subjects usingthe technique of enhanced depth imaging optical coherence tomography," Investigative Ophthalmology and Visual Science, vol. 52, no. 5, pp. 2267-2271, 2011.

[13] T. Fujiwara, Y. Imamura, R. Margolis, J. S. Slakter, and R. F. Spaide, "Enhanced depth imaging optical coherence tomography of the choroid in highly myopic eyes," The American Journal of Ophthalmology, vol. 148, no. 3, pp. 445-450, 2009.

[14] R. L. Wong, P. Zhao, and W. W. Lai, "Choroidal thickness in relation to hypercholesterolemia on enhanced depth imaging optical coherence tomography," Retina, vol. 33, no. 2, pp. 423428, 2013.

[15] V. Kajic, M. Esmaeelpour, B. PovaZay, D. Marshall, P. L. Rosin, and W. Drexler, "Automated choroidal segmentation of1060 nm OCT in healthy and pathologic eyes using a statistical model," Biomedical Optics Express, vol. 3, no. 1, pp. 86-103, 2012.

[16] M. Esmaeelpour, S. Brunner, S. Ansari-Shahrezaei et al., "Cho-roidal thinning in diabetes type 1 detected by 3-dimensional 1060 nm optical coherence tomography," Investigative Ophthalmology and Visual Science, vol. 53, no. 11, pp. 6803-6809, 2012.

[17] V. Kajic, M. Esmaeelpour, C. Glittenberg et al., "Automated three-dimensional choroidal vessel segmentation of 3D, 1060 nm OCT retinal data," Biomedical Optics Express,vol. 4, no. 1, pp. 134-150, 2013.

[18] J. Tian, P. Marziliano, M. Baskaran, T. A. Tun, and T. Aung, "Automatic measurements of choroidal thickness in EDI-OCT images," in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC '12), pp. 5360-5363, San Diego, Calif, USA, September 2012.

[19] J. Tian, P. Marziliano, M. Baskaran, T. A. Tun, and T. Aung, "Automatic segmentation of the choroid in enhanced depth imaging optical coherence tomography images," Biomedical Optics Express, vol. 4, no. 3, pp. 397-411, 2013.

[20] H. Lu, N. Boonarpha, M. T. Kwong, and Y. Zheng, "Automated segmentation of the choroid in retinal optical coherence tomography images," in Proceedings of the 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC '13), pp. 5869-5872, Osaka, Japan, July 2013.

[21] M. K. Garvin, M. D. Abràmoff, R. Kardon, S. R. Russell, X. Wu, and M. Sonka, "Intraretinal layer segmentation of macular optical coherence tomography images using optimal 3-D graph search," IEEE Transactions on Medical Imaging, vol. 27, no. 10, pp. 1495-1505, 2008.

[22] http://www.heidelbergengineering.com/us/wp-content/ uploads/hra2-aquire-the-perfect-image.pdf.

[23] R. Bellman, "The theory of dynamic programming," Bulletin of the American Mathematical Society, vol. 60, pp. 503-515,1954.

[24] M. Sonka, V. Hlavac, and R. Boyle, Image Processing, Analysis, and Machine Vision, Brooks/Cole, Albany, NY, USA, 1999.

[25] S. Arivazhagan and L. Ganesan, "Texture classification using wavelet transform," Pattern Recognition Letters, vol. 24, no. 910, pp. 1513-1521, 2003.

[26] G. van de Wouwer, P. Scheunders, and D. van Dyck, "Statistical texture characterization from discrete wavelet representations," IEEE Transactions on Image Processing, vol. 8, no. 4, pp. 592598, 1999.

[27] C. E. Rasmussen, "The infinite Gaussian mixture model," Advances in Neural Information Processing Systems, vol. 12, pp. 554-560, 2000.

[28] D. Greig, B. Porteous, and A. H. Seheult, "Exact maximum a posteriori estimation for binary images," Journal of the Royal Statistical Society B, vol. 51, no. 2, pp. 271-279,1989.

[29] Y. Boykov and G. Funka-Lea, "Graph cuts and efficient N-D image segmentation," International Journal of Computer Vision, vol. 70, no. 2, pp. 109-131, 2006.

[30] Y. Boykov and V. Kolmogorov, "Computing geodesics and minimal surfaces via graph cuts," in Proceedings of the 9th IEEE International Conference on Computer Vision, pp. 26-33, October 2003.

[31] S. Geman and D. Geman, "Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images," IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 6, no. 6, pp. 721-741, 1984.

[32] A. V. Goldberg and R. E. Tarjan, "New approach to the maximum-flow problem," Journal of the ACM, vol. 35, no. 4, pp. 921-940, 1988.

[33] L. R. Ford and D. R. Fulkerson, "Maximal flow through a network," Canadian Journal of Mathematics, vol. 8, pp. 399-404, 1956.

Copyright of Computational & Mathematical Methods in Medicine is the property of Hindawi Publishing Corporation and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.