Region Based Contour Detection by Dynamic ... - Semantic Scholar
Region Based Contour Detection by Dynamic Programming Xiaoyi Jiang1,2,3 and Daniel Tenbrinck1,2 1
Department of Mathematics and Computer Science, University of M¨ unster, Germany 2 European Institute for Molecular Imaging, University of M¨ unster, Germany 3 Cluster of Excellence EXC 1003, Cells in Motion, CiM, M¨ unster, Germany
Abstract. Dynamic programming (DP) is a popular technique for contour detection, particularly in biomedical image analysis. Although gradient information is typically used in such methods, it is not always a reliable measure to work with and there is a strong need of non-gradient based methods. In this paper we present a general framework for region based contour detection by dynamic programming. It is based on a global energy function which is approximated by a radial ray-wise summation to enable dynamic programming. Its simple algorithmic structure allows to use arbitrarily complex region models and model testing functions, in particular by means of techniques from robust statistics. The proposed framework was tested on synthetic data and real microscopic images. A performance comparison with the standard gradient-based DP and a recent non-gradient DP-based contour detection algorithm clearly demonstrates the superiority of our approach.
1
Introduction
Dynamic programming is a popular technique for contour detection due to its simpleness, efficiency, and guarantee of optimality. One class of detectable contours starts from the left, passes each image column exactly once, and ends in the last column, as exemplified by detecting arterial walls in sonographic artery images [2]. Another, even more important, application class deals with closed contours. Based on a point p in the interior of the contour, a polar transformation with p being the central point brings the original image into a matrix, in which a closed contour becomes one from left to right afterwards. Finally, the detected contour has to be transformed back to the original image space. This technique works well for star-shaped contours, particularly including convex contours. Note that special care must be taken in order to guarantee the closedness of the detected contour [8]. In the following we use the case of closed contours for our discussion, but it applies to the first class of contours as well. Typically, DP-based contour detection assumes strong edges along the contour and is thus based on gradient computation. In the simplest case the contour measure is defined by the sum of gradient magnitudes. More sophisticated formulations can be found for instance in [5], where the object contour is expected
(a)
(b)
(c)
(d)
Fig. 1: (a) Tumor cell ROI; (b) gradient; (c) gradient-based optimal contour; (d) region-based optimal contour from our approach (with L1 norm, see Section 3).
to have approximately uniform edge strength, leading to an alternative contour measure by ratio of the mean and the standard deviation of image gradients. In practice, however, gradient is not always a reliable measure to work with, in particular in biomedical image analysis, due to limiting factors like low contrast, high noise level, etc. Generally, objects may have very smooth or even discontinuous boundaries. One such example is the region-of-interest (ROI) of a tumor cell from microscopic imaging shown in Figure 1. Based on the gradient alone, the contour in (c) is rather suboptimal in the application situation, but maximizes the sum of gradient magnitude (6428 in this case). On the other hand, the contour in (d) is much more reasonable although having a lower sum of gradient magnitude (1917). Looking at the gradient in (b), it indeed can be expected that a small change of intensity values may lead to substantially different contours when maximizing the sum of gradient magnitude. This example clearly demonstrates the need of alternate, non-gradient approaches. There are only very few works on DP-based contour detection using nongradient information. The approach from [9] considers two circular areas, one on each side of a point along the radial ray. The squared difference of their mean intensity values is interpreted as a measure of separability and used as cost for dynamic programming. The circular form for evaluation, however, leads to severe contamination by samples of other distributions, even for many convex shapes, thus substantially limiting the applicability. The authors of [4] present a homogeneity splitting cost for dynamic programming. This method will be further detailed and experimentally compared to our method in Section 4. In the next section we start with a discussion of the fundamental idea of region-based segmentation. Then, Section 3 presents our framework of DP-based contour detection. Our experimental results and comparison with other methods are given in Section 4. Finally, a brief discussion concludes this paper.
2
Fundamentals of region-based segmentation
Given an image f (x, y), a contour C partitions the image into two regions, inside(C) and outside(C). We assume each region is well represented by some model, which can be validated by a model testing function Fi and Fo , respectively. Then, the segmentation task can be generally formulated as one of opti-
mizing the following energy function: Z
Z
E(C) =
Fi (x, y)dxdy +
Fo (x, y)dxdy
inside(C)
(1)
outside(C)
If needed, the two terms can be further weighted. The seminal work by Chan and Vese [1] is a prominent example of this regionbased segmentation scheme with the assumption of two regions with approximately piecewise-constant intensities. Let c1 and c2 denote the average intensity of inside(C) and outside(C), respectively. The energy function in [1] is defined by (additional regularizing terms are ignored): Z
Z
|f (x, y) − c1 |2 dxdy +
E(C) = inside(C)
|f (x, y) − c2 |2 dxdy
(2)
outside(C)
In [1] the optimization is done iteratively by means of a level set solution. Such global optimization schemes suffer from several potential problems: Dependence of the final contour on the initialization, danger of falling into local optima, and high computational burden. For star-shaped contours, which are of interest in many real applications, they are not necessarily the best choice. In this case dynamic programming is applicable and its simpleness, efficiency, and guarantee of optimality make it an attractive option. In this work we explore the potential of region-based dynamic programming for contour detection.
3
Region based segmentation by dynamic programming
A star-shaped contour C can be represented in polar form r(θ), θ ∈ [0, 2π). Given the image boundary B(θ), θ ∈ [0, 2π), the energy function in (1) becomes: Z
Z
E(C) =
Fi (x, y)dxdy +
Fo (x, y)dxdy
inside(C)
Z
2π
outside(C)
r(θ)
Z
=
2π
Z
Z
B(θ)
Fi (θ, r)drdθ + 0
Z
2π
=
0
hZ
0
r(θ)
Z
=
B(θ)
Fi (θ, r)dr +
0
| Z
Fo (θ, r)drdθ
0
r(θ)
i Fo (θ, r)dr dθ
r(θ)
{z
Eθ (C)
}
2π
Eθ (C)dθ
(3)
0
The term Eθ (C) is simply a model testing function applied to one particular radial ray corresponding to angle θ. Note that the model parameters have to be estimated by the entirety of inside(C) and outside(C), respectively.
Fig. 2: Image in polar space. A star-shaped contour in the original image corresponds to a continuous contour from left to right in polar space.
In this formulation the Chan-Vese energy function (2) becomes: Z r(θ) Z B(θ) Eθ (C) = |f (θ, r) − c1 |2 dr + |f (θ, r) − c2 |2 dr 0
R 2π R r(θ) c1 =
0
0 R 2π 0
(4)
r(θ)
f (θ, r)dr
R r(θ) 0
dr
R 2π R B(θ) ,
c2 =
0
f (θ, r)dr r(θ) R 2π R B(θ) dr 0 r(θ)
Note that the two average intensities c1 and c2 depend on the entire inside(C) and outside(C), respectively, even when evaluating Eθ (C) for a single radial ray. If the two average intensity c1 and c2 were known and thus not needed to be estimated, then Eθ (C) in (4) can be easily computed for each candidate r of a given θ. That is, we obtain a cost matrix c(θ, r) := Eθ (C) for all image points in polar space, see Figure 2. A star-shaped contour in original image corresponds to a continuous contour from left to right in this polar space, which can be detected by means of dynamic programming [4] in a time linear to the size of polar space, thus also linear to the original image size. Unfortunately, the two average intensities c1 and c2 are unknown before the contour is completely detected. Therefore, we cannot compute the cost matrix c(θ, r) beforehand and consequently cannot apply dynamic programming to solve the optimization problem (4). The same is true for general Eθ (C) since both model testing functions Fi (θ, r) and Fo (θ, r) depend on global model parameters of inside(C) and outside(C), respectively. The optimization problem (4), however, does have a dynamic programming solution if we model each radial ray separately. Then, a radial ray is partitioned into an inside part r ∈ [0, r(θ)] with the average intensity cθ,1 and an outside part r ∈ [r(θ), B(θ)] with the average intensity cθ,2 . Here both average intensities are locally determined within a specific radial ray. Consequently, Eq.(4) becomes: Z r(θ) Z B(θ) Eθ (C) = |f (θ, r) − cθ,1 |2 dr + |f (θ, r) − cθ,2 |2 dr (5) 0
r(θ)
It can be computed for each r(θ) of a given θ to obtain a cost matrix c(θ, r) := Eθ (C) as follows. For a given θ, each potential value for r(θ) partitions the radial ray into an inside part r ∈ [0, r(θ)] and an outside part r ∈ [r(θ), B(θ)]. This partitioning immediately gives us the respective average intensity cθ,1 and cθ,2 , thus
R enabling the evaluation of the model testing function |f (θ, r) − cθ,k |2 , k = 1, 2. Given the computed cost matrix c(θ, r), dynamic programming can be applied to detect the optimal contour. In general we can formulate: Z Eθ (C) =
r(θ)
Fi∗ (θ, r)dr +
0
Z
B(θ)
Fo∗ (θ, r)dr
r(θ)
In contrast to (3), the representation models and the resulting model testing functions Fi∗ (θ, r) and Fo∗ (θ, r) are bounded to a particular radial ray θ instead of the whole image. Also in this general case we have no trouble to compute a cost matrix c(θ, r) and then to apply dynamic programming for detecting an optimal contour in polar space. This general approach can be flexibly concretized in many different ways: – Any variant of average estimation can be used for the L22 model testing function in (5). In particular, robust statistical estimation methods [6] provide more reliable results for noisy data. – Any model testing function T (f (θ, r), cθ,k ), k = 1, 2, can be used, e.g. variants of Lp norm. The special case in (5) is the squared L2 norm. Of particular interest is the L1 norm, which is provably more robust than L2 . 2 2 – Also variance-based model testing is possible, e.g. E(C) = σθ,1 + σθ,2 , where 2 2 σθ,1 and σθ,2 are variance of inside(C) and outside(C), respectively. So far, the inside and outside part are tested separately. We can go a step further and use combined testing. One example is: Eθ (C) =
|cθ,1 − cθ,2 |2 2 + σ2 σθ,1 θ,2
in accordance with the Fisher linear discriminant. In this most general form the term Eθ (C) cannot be reformulated as a sum of inside and outside term. In summary, the original optimization task (3) cannot be solved by dynamic programming in general. But by modeling each radial ray separately a dynamic programming solution is possible for any representation model and model testing function independent of their form and complexity. This universality gives the rather simple scheme of dynamic programming considerable power for real-world applications. For instance, robust estimation methods such as median-based approaches and L1 norm are highly desired for improved robustness. Also, sophisticated testing criteria like Fisher linear discriminant and others from machine learning theory provides extra useful options for measuring the separability of two distributions. While all such variants are straightforward in our DP-based framework, they are rather complex to implement in calculus, e.g., in the context of level set methods. The proposed approach of modeling each radial ray separately can be interpreted in two ways. First, it can be seen as an approximate solution of the optimization problem. We assume all radial rays are approximately subject to the
(a)
(b)
(c)
(d)
(e)
Fig. 3: (a) Tumor cell ROI; (b) contour by our approach (L1 norm); (c) our approach (L22 ); (d) gradient-based approach; (e) method from [4]. The first two rows are early stage and the last two are late stage cells.
representation model of the whole image. Taking the example of (5), it means cθ,1 ≈ c1 , cθ,2 ≈ c2 , θ ∈ [0, 2π). This assumption may be violated by some of the radial rays. Even in such a case the dynamic programming algorithm still has a good chance of finding a globally optimal solution due to the forced continuity. On the other hand, the dynamic programming solution can also been seen as an extension of the original optimization problem with extra power. The separate modeling of the radial rays allows to deal with inhomogeneous inside and outside regions. As soon as a separation in the individual radial rays can be assumed, an overall segmentation is still possible despite of the inhomogeneity.
4
Experimental results
We have tested our approach on a series of tumor cell images from microscopic imaging. Given an image, cells of interest are marked by a biologist. Then, the approximate cell center and a ROI around the center point are manually selected. Note that we currently refrain to use an automatic detection of cell positions because the microscopic cell images have been acquired with a series of focus levels and the biologists are only interested in well focused cells. Thus, an automatic cell detection must contain a focus assessment which turns out to be not easy, even for an untrained person. The following four algorithm variants have been tested and compared. – Our approach with T (f (θ, r), cθ,k ) = L1 , k = 1, 2. – Our approach with T (f (θ, r), cθ,k ) = L22 , k = 1, 2, see Eq. (5). In both variants of our approach cθ,k is determined by a simple mean computation.
1 0.9 0.8
Dice index
0.7 0.6 0.5 0.4
Gradient−based Malon et al. Our approach (L1) Our approach (L22)
0.3 0.2 0.1 0 5
10
15
20
25
30 35 40 Noise variance σ2
45
50
55
60
Fig. 4: Four examples of synthe- Fig. 5: Segmentation accuracy (Dice index) sized cell images. with increasing noise variance levels.
– Gradient-based contour detection, i.e., Eθ (C) := GradientMagnitude(θ, r(θ)). In this case the dynamic programming is applied to maximize the sum of gradient magnitudes (computed by the Sobel operator) along the contour. – The homogeneity splitting method from [4]: Eθ (C) = σθ (r(θ)) − σθ (r(θ) + δ) + 1. Here σθ (r) represents the standard deviation of inner segment [0, r] on the radial ray θ and δ is a smoothing parameter. In our experiments this parameter was optimized manually. The performance of the four algorithm variants is exemplarily demonstrated on the four images in Figure 3. Our segmentation results appear very convincing to the biologists, in particular using the L1 norm. In contrast, both the standard gradient-based approach and the method from [4] fail in case of early stage cells with complex inner-cell textures which produce substantial ”misleading” gradient hints. Similar reason explains the failure of [4]. The textures induce severe variations of standard deviations along a radial ray, thus leading to a deficient homogeneity splitting measure. On the other hand, late stage cells cause less problems. These images are typically of high contrast to the background and do not have complex inner-cell textures. Even the standard gradient-based approach works well. However, the method from [4] fails in some cases. We have also conducted a study on synthetic data to quantitatively evaluate and compare the performance of the four algorithm variants. Random starshaped objects (similar to cells from microscopic images) were generated with constant intensity value of 60 in a homogeneous background with intensity 0. We added additive Gaussian noise with mean 0 and different variance levels σ 2 ∈ {5, 10, . . . , 55, 60} (see Figure 4 for four examples with increasing level of noise 5, 20, 45, and 60 from top left to bottom right) and evaluated the four discussed methods on 10 arbitrary shapes per noise level. For quantification we compare the computed segmentation result with the (unperturbed) ground truth images using the Dice index D(A, B) = 2|A ∩ B|/(|A| + |B|), see Figure 5. This study indicates that the proposed region-based models are significantly more robust under the impact of noise and give good segmentation results even for very noisy data. In contrast, the gradient-based method cannot cope with the noisy
edge information induced by the perturbations. The approach from [4] is very sensitive to inhomogeneities, since its homogeneity criterion is solely based on the inside region and ignores the outside region. These quantitative results are consistent with our observations on real microscopic images as discussed above.
5
Conclusion
In this paper we have presented a general framework of region based dynamic programming technique for contour detection. This non-gradient approach is efficient (linear to the image size) and allows to use arbitrarily complex region model and model testing functions. In particular, techniques from robust statistics and machine learning theory are helpful when dealing with noisy data. A performance comparison has clearly demonstrated the superiority of our approach. An extension to simultaneously detecting multiple contours is straightforward [7]. The principle behind this work can also be applied to 3D surface detection [3]. The current formulation is based on the fundamental assumption that all the radial rays are approximately subject to the representation model of the whole image. This assumption is always violated to some extent. Although this fact may not effect the contour detection due to the forced continuity, multiple radial rays together as higher-level processing unit instead of single rays will reduce this modeling uncertainty. All these extensions will be studied in future. Acknowledgments: The authors thank Kristen Mills at Max Planck Institute for Intelligent Systems, Stuttgart, for providing the microscopic images.
References 1. Chan, T.F., Vese, L.A.: Active contours without edges. IEEE Trans. on Image Processing 10(2), 266–277 (2001) 2. Cheng, D.C., Jiang, X.: Detections of arterial wall in sonographic artery images using dual dynamic programming. IEEE Trans. on Information Technology in Biomedicine 12(6), 792–799 (2008) 3. Li, K., Wu, X., Chen, D., Sonka, M.: Optimal surface segmentation in volumetric images - a graph-theoretic approach. IEEE Trans. on PAMI 28(1), 119–134 (2006) 4. Malon, C., Cosatto, E.: Dynamic radial contour extraction by splitting homogeneous areas. In: Proc. of CAIP (1). pp. 269–277 (2011) 5. Ray, N., Acton, S.T., Zhang, H.: Seeing through clutter: Snake computation with dynamic programming for particle segmentation. In: Proc. of ICPR. pp. 801–804 (2012) 6. Stewart, C.: Robust parameter estimation in computer vision. SIAM Reviews 41(3), 513–537 (1999) 7. Sun, C., Appleton, B.: Multiple paths extraction in images using a constrained expanded trellis. IEEE Trans. on PAMI 27(12), 1923–1933 (2005) 8. Sun, C., Pallottino, S.: Circular shortest path in images. Pattern Recognition 36(3), 709–719 (2003) 9. Yu, M., Huang, Q., Jin, R., Song, E., Liu, H., Hung, C.C.: A novel segmentation method for convex lesions based on dynamic programming with local intra-class variance. In: Proc. of ACM Symposium on Applied Computing. pp. 39–44 (2012)
Department of Mathematics and Computer Science, University of M¨ unster, Germany 2 European Institute for Molecular Imaging, University of M¨ unster, Germany 3 Cluster of Excellence EXC 1003, Cells in Motion, CiM, M¨ unster, Germany
Abstract. Dynamic programming (DP) is a popular technique for contour detection, particularly in biomedical image analysis. Although gradient information is typically used in such methods, it is not always a reliable measure to work with and there is a strong need of non-gradient based methods. In this paper we present a general framework for region based contour detection by dynamic programming. It is based on a global energy function which is approximated by a radial ray-wise summation to enable dynamic programming. Its simple algorithmic structure allows to use arbitrarily complex region models and model testing functions, in particular by means of techniques from robust statistics. The proposed framework was tested on synthetic data and real microscopic images. A performance comparison with the standard gradient-based DP and a recent non-gradient DP-based contour detection algorithm clearly demonstrates the superiority of our approach.
1
Introduction
Dynamic programming is a popular technique for contour detection due to its simpleness, efficiency, and guarantee of optimality. One class of detectable contours starts from the left, passes each image column exactly once, and ends in the last column, as exemplified by detecting arterial walls in sonographic artery images [2]. Another, even more important, application class deals with closed contours. Based on a point p in the interior of the contour, a polar transformation with p being the central point brings the original image into a matrix, in which a closed contour becomes one from left to right afterwards. Finally, the detected contour has to be transformed back to the original image space. This technique works well for star-shaped contours, particularly including convex contours. Note that special care must be taken in order to guarantee the closedness of the detected contour [8]. In the following we use the case of closed contours for our discussion, but it applies to the first class of contours as well. Typically, DP-based contour detection assumes strong edges along the contour and is thus based on gradient computation. In the simplest case the contour measure is defined by the sum of gradient magnitudes. More sophisticated formulations can be found for instance in [5], where the object contour is expected
(a)
(b)
(c)
(d)
Fig. 1: (a) Tumor cell ROI; (b) gradient; (c) gradient-based optimal contour; (d) region-based optimal contour from our approach (with L1 norm, see Section 3).
to have approximately uniform edge strength, leading to an alternative contour measure by ratio of the mean and the standard deviation of image gradients. In practice, however, gradient is not always a reliable measure to work with, in particular in biomedical image analysis, due to limiting factors like low contrast, high noise level, etc. Generally, objects may have very smooth or even discontinuous boundaries. One such example is the region-of-interest (ROI) of a tumor cell from microscopic imaging shown in Figure 1. Based on the gradient alone, the contour in (c) is rather suboptimal in the application situation, but maximizes the sum of gradient magnitude (6428 in this case). On the other hand, the contour in (d) is much more reasonable although having a lower sum of gradient magnitude (1917). Looking at the gradient in (b), it indeed can be expected that a small change of intensity values may lead to substantially different contours when maximizing the sum of gradient magnitude. This example clearly demonstrates the need of alternate, non-gradient approaches. There are only very few works on DP-based contour detection using nongradient information. The approach from [9] considers two circular areas, one on each side of a point along the radial ray. The squared difference of their mean intensity values is interpreted as a measure of separability and used as cost for dynamic programming. The circular form for evaluation, however, leads to severe contamination by samples of other distributions, even for many convex shapes, thus substantially limiting the applicability. The authors of [4] present a homogeneity splitting cost for dynamic programming. This method will be further detailed and experimentally compared to our method in Section 4. In the next section we start with a discussion of the fundamental idea of region-based segmentation. Then, Section 3 presents our framework of DP-based contour detection. Our experimental results and comparison with other methods are given in Section 4. Finally, a brief discussion concludes this paper.
2
Fundamentals of region-based segmentation
Given an image f (x, y), a contour C partitions the image into two regions, inside(C) and outside(C). We assume each region is well represented by some model, which can be validated by a model testing function Fi and Fo , respectively. Then, the segmentation task can be generally formulated as one of opti-
mizing the following energy function: Z
Z
E(C) =
Fi (x, y)dxdy +
Fo (x, y)dxdy
inside(C)
(1)
outside(C)
If needed, the two terms can be further weighted. The seminal work by Chan and Vese [1] is a prominent example of this regionbased segmentation scheme with the assumption of two regions with approximately piecewise-constant intensities. Let c1 and c2 denote the average intensity of inside(C) and outside(C), respectively. The energy function in [1] is defined by (additional regularizing terms are ignored): Z
Z
|f (x, y) − c1 |2 dxdy +
E(C) = inside(C)
|f (x, y) − c2 |2 dxdy
(2)
outside(C)
In [1] the optimization is done iteratively by means of a level set solution. Such global optimization schemes suffer from several potential problems: Dependence of the final contour on the initialization, danger of falling into local optima, and high computational burden. For star-shaped contours, which are of interest in many real applications, they are not necessarily the best choice. In this case dynamic programming is applicable and its simpleness, efficiency, and guarantee of optimality make it an attractive option. In this work we explore the potential of region-based dynamic programming for contour detection.
3
Region based segmentation by dynamic programming
A star-shaped contour C can be represented in polar form r(θ), θ ∈ [0, 2π). Given the image boundary B(θ), θ ∈ [0, 2π), the energy function in (1) becomes: Z
Z
E(C) =
Fi (x, y)dxdy +
Fo (x, y)dxdy
inside(C)
Z
2π
outside(C)
r(θ)
Z
=
2π
Z
Z
B(θ)
Fi (θ, r)drdθ + 0
Z
2π
=
0
hZ
0
r(θ)
Z
=
B(θ)
Fi (θ, r)dr +
0
| Z
Fo (θ, r)drdθ
0
r(θ)
i Fo (θ, r)dr dθ
r(θ)
{z
Eθ (C)
}
2π
Eθ (C)dθ
(3)
0
The term Eθ (C) is simply a model testing function applied to one particular radial ray corresponding to angle θ. Note that the model parameters have to be estimated by the entirety of inside(C) and outside(C), respectively.
Fig. 2: Image in polar space. A star-shaped contour in the original image corresponds to a continuous contour from left to right in polar space.
In this formulation the Chan-Vese energy function (2) becomes: Z r(θ) Z B(θ) Eθ (C) = |f (θ, r) − c1 |2 dr + |f (θ, r) − c2 |2 dr 0
R 2π R r(θ) c1 =
0
0 R 2π 0
(4)
r(θ)
f (θ, r)dr
R r(θ) 0
dr
R 2π R B(θ) ,
c2 =
0
f (θ, r)dr r(θ) R 2π R B(θ) dr 0 r(θ)
Note that the two average intensities c1 and c2 depend on the entire inside(C) and outside(C), respectively, even when evaluating Eθ (C) for a single radial ray. If the two average intensity c1 and c2 were known and thus not needed to be estimated, then Eθ (C) in (4) can be easily computed for each candidate r of a given θ. That is, we obtain a cost matrix c(θ, r) := Eθ (C) for all image points in polar space, see Figure 2. A star-shaped contour in original image corresponds to a continuous contour from left to right in this polar space, which can be detected by means of dynamic programming [4] in a time linear to the size of polar space, thus also linear to the original image size. Unfortunately, the two average intensities c1 and c2 are unknown before the contour is completely detected. Therefore, we cannot compute the cost matrix c(θ, r) beforehand and consequently cannot apply dynamic programming to solve the optimization problem (4). The same is true for general Eθ (C) since both model testing functions Fi (θ, r) and Fo (θ, r) depend on global model parameters of inside(C) and outside(C), respectively. The optimization problem (4), however, does have a dynamic programming solution if we model each radial ray separately. Then, a radial ray is partitioned into an inside part r ∈ [0, r(θ)] with the average intensity cθ,1 and an outside part r ∈ [r(θ), B(θ)] with the average intensity cθ,2 . Here both average intensities are locally determined within a specific radial ray. Consequently, Eq.(4) becomes: Z r(θ) Z B(θ) Eθ (C) = |f (θ, r) − cθ,1 |2 dr + |f (θ, r) − cθ,2 |2 dr (5) 0
r(θ)
It can be computed for each r(θ) of a given θ to obtain a cost matrix c(θ, r) := Eθ (C) as follows. For a given θ, each potential value for r(θ) partitions the radial ray into an inside part r ∈ [0, r(θ)] and an outside part r ∈ [r(θ), B(θ)]. This partitioning immediately gives us the respective average intensity cθ,1 and cθ,2 , thus
R enabling the evaluation of the model testing function |f (θ, r) − cθ,k |2 , k = 1, 2. Given the computed cost matrix c(θ, r), dynamic programming can be applied to detect the optimal contour. In general we can formulate: Z Eθ (C) =
r(θ)
Fi∗ (θ, r)dr +
0
Z
B(θ)
Fo∗ (θ, r)dr
r(θ)
In contrast to (3), the representation models and the resulting model testing functions Fi∗ (θ, r) and Fo∗ (θ, r) are bounded to a particular radial ray θ instead of the whole image. Also in this general case we have no trouble to compute a cost matrix c(θ, r) and then to apply dynamic programming for detecting an optimal contour in polar space. This general approach can be flexibly concretized in many different ways: – Any variant of average estimation can be used for the L22 model testing function in (5). In particular, robust statistical estimation methods [6] provide more reliable results for noisy data. – Any model testing function T (f (θ, r), cθ,k ), k = 1, 2, can be used, e.g. variants of Lp norm. The special case in (5) is the squared L2 norm. Of particular interest is the L1 norm, which is provably more robust than L2 . 2 2 – Also variance-based model testing is possible, e.g. E(C) = σθ,1 + σθ,2 , where 2 2 σθ,1 and σθ,2 are variance of inside(C) and outside(C), respectively. So far, the inside and outside part are tested separately. We can go a step further and use combined testing. One example is: Eθ (C) =
|cθ,1 − cθ,2 |2 2 + σ2 σθ,1 θ,2
in accordance with the Fisher linear discriminant. In this most general form the term Eθ (C) cannot be reformulated as a sum of inside and outside term. In summary, the original optimization task (3) cannot be solved by dynamic programming in general. But by modeling each radial ray separately a dynamic programming solution is possible for any representation model and model testing function independent of their form and complexity. This universality gives the rather simple scheme of dynamic programming considerable power for real-world applications. For instance, robust estimation methods such as median-based approaches and L1 norm are highly desired for improved robustness. Also, sophisticated testing criteria like Fisher linear discriminant and others from machine learning theory provides extra useful options for measuring the separability of two distributions. While all such variants are straightforward in our DP-based framework, they are rather complex to implement in calculus, e.g., in the context of level set methods. The proposed approach of modeling each radial ray separately can be interpreted in two ways. First, it can be seen as an approximate solution of the optimization problem. We assume all radial rays are approximately subject to the
(a)
(b)
(c)
(d)
(e)
Fig. 3: (a) Tumor cell ROI; (b) contour by our approach (L1 norm); (c) our approach (L22 ); (d) gradient-based approach; (e) method from [4]. The first two rows are early stage and the last two are late stage cells.
representation model of the whole image. Taking the example of (5), it means cθ,1 ≈ c1 , cθ,2 ≈ c2 , θ ∈ [0, 2π). This assumption may be violated by some of the radial rays. Even in such a case the dynamic programming algorithm still has a good chance of finding a globally optimal solution due to the forced continuity. On the other hand, the dynamic programming solution can also been seen as an extension of the original optimization problem with extra power. The separate modeling of the radial rays allows to deal with inhomogeneous inside and outside regions. As soon as a separation in the individual radial rays can be assumed, an overall segmentation is still possible despite of the inhomogeneity.
4
Experimental results
We have tested our approach on a series of tumor cell images from microscopic imaging. Given an image, cells of interest are marked by a biologist. Then, the approximate cell center and a ROI around the center point are manually selected. Note that we currently refrain to use an automatic detection of cell positions because the microscopic cell images have been acquired with a series of focus levels and the biologists are only interested in well focused cells. Thus, an automatic cell detection must contain a focus assessment which turns out to be not easy, even for an untrained person. The following four algorithm variants have been tested and compared. – Our approach with T (f (θ, r), cθ,k ) = L1 , k = 1, 2. – Our approach with T (f (θ, r), cθ,k ) = L22 , k = 1, 2, see Eq. (5). In both variants of our approach cθ,k is determined by a simple mean computation.
1 0.9 0.8
Dice index
0.7 0.6 0.5 0.4
Gradient−based Malon et al. Our approach (L1) Our approach (L22)
0.3 0.2 0.1 0 5
10
15
20
25
30 35 40 Noise variance σ2
45
50
55
60
Fig. 4: Four examples of synthe- Fig. 5: Segmentation accuracy (Dice index) sized cell images. with increasing noise variance levels.
– Gradient-based contour detection, i.e., Eθ (C) := GradientMagnitude(θ, r(θ)). In this case the dynamic programming is applied to maximize the sum of gradient magnitudes (computed by the Sobel operator) along the contour. – The homogeneity splitting method from [4]: Eθ (C) = σθ (r(θ)) − σθ (r(θ) + δ) + 1. Here σθ (r) represents the standard deviation of inner segment [0, r] on the radial ray θ and δ is a smoothing parameter. In our experiments this parameter was optimized manually. The performance of the four algorithm variants is exemplarily demonstrated on the four images in Figure 3. Our segmentation results appear very convincing to the biologists, in particular using the L1 norm. In contrast, both the standard gradient-based approach and the method from [4] fail in case of early stage cells with complex inner-cell textures which produce substantial ”misleading” gradient hints. Similar reason explains the failure of [4]. The textures induce severe variations of standard deviations along a radial ray, thus leading to a deficient homogeneity splitting measure. On the other hand, late stage cells cause less problems. These images are typically of high contrast to the background and do not have complex inner-cell textures. Even the standard gradient-based approach works well. However, the method from [4] fails in some cases. We have also conducted a study on synthetic data to quantitatively evaluate and compare the performance of the four algorithm variants. Random starshaped objects (similar to cells from microscopic images) were generated with constant intensity value of 60 in a homogeneous background with intensity 0. We added additive Gaussian noise with mean 0 and different variance levels σ 2 ∈ {5, 10, . . . , 55, 60} (see Figure 4 for four examples with increasing level of noise 5, 20, 45, and 60 from top left to bottom right) and evaluated the four discussed methods on 10 arbitrary shapes per noise level. For quantification we compare the computed segmentation result with the (unperturbed) ground truth images using the Dice index D(A, B) = 2|A ∩ B|/(|A| + |B|), see Figure 5. This study indicates that the proposed region-based models are significantly more robust under the impact of noise and give good segmentation results even for very noisy data. In contrast, the gradient-based method cannot cope with the noisy
edge information induced by the perturbations. The approach from [4] is very sensitive to inhomogeneities, since its homogeneity criterion is solely based on the inside region and ignores the outside region. These quantitative results are consistent with our observations on real microscopic images as discussed above.
5
Conclusion
In this paper we have presented a general framework of region based dynamic programming technique for contour detection. This non-gradient approach is efficient (linear to the image size) and allows to use arbitrarily complex region model and model testing functions. In particular, techniques from robust statistics and machine learning theory are helpful when dealing with noisy data. A performance comparison has clearly demonstrated the superiority of our approach. An extension to simultaneously detecting multiple contours is straightforward [7]. The principle behind this work can also be applied to 3D surface detection [3]. The current formulation is based on the fundamental assumption that all the radial rays are approximately subject to the representation model of the whole image. This assumption is always violated to some extent. Although this fact may not effect the contour detection due to the forced continuity, multiple radial rays together as higher-level processing unit instead of single rays will reduce this modeling uncertainty. All these extensions will be studied in future. Acknowledgments: The authors thank Kristen Mills at Max Planck Institute for Intelligent Systems, Stuttgart, for providing the microscopic images.
References 1. Chan, T.F., Vese, L.A.: Active contours without edges. IEEE Trans. on Image Processing 10(2), 266–277 (2001) 2. Cheng, D.C., Jiang, X.: Detections of arterial wall in sonographic artery images using dual dynamic programming. IEEE Trans. on Information Technology in Biomedicine 12(6), 792–799 (2008) 3. Li, K., Wu, X., Chen, D., Sonka, M.: Optimal surface segmentation in volumetric images - a graph-theoretic approach. IEEE Trans. on PAMI 28(1), 119–134 (2006) 4. Malon, C., Cosatto, E.: Dynamic radial contour extraction by splitting homogeneous areas. In: Proc. of CAIP (1). pp. 269–277 (2011) 5. Ray, N., Acton, S.T., Zhang, H.: Seeing through clutter: Snake computation with dynamic programming for particle segmentation. In: Proc. of ICPR. pp. 801–804 (2012) 6. Stewart, C.: Robust parameter estimation in computer vision. SIAM Reviews 41(3), 513–537 (1999) 7. Sun, C., Appleton, B.: Multiple paths extraction in images using a constrained expanded trellis. IEEE Trans. on PAMI 27(12), 1923–1933 (2005) 8. Sun, C., Pallottino, S.: Circular shortest path in images. Pattern Recognition 36(3), 709–719 (2003) 9. Yu, M., Huang, Q., Jin, R., Song, E., Liu, H., Hung, C.C.: A novel segmentation method for convex lesions based on dynamic programming with local intra-class variance. In: Proc. of ACM Symposium on Applied Computing. pp. 39–44 (2012)
