Image Segmentation via Mean Curvature Regularized Mumford-Shah Model and Thresholding

Abstract

Due to the limitations in imaging devices and subject-induced susceptibility effect, general image segmentation is still an open problem. Typical challenges include image noise, intensity inhomogeneity and various image modalities. In this paper, we propose to use a two-step strategy. Specifically, we first utilize a mean curvature regularized Mumford-Shah model to recover an intermediate image with enhanced saliency, and then the segmentation is obtained by a thresholding procedure. For images with intensity inhomogeneity, a bias-corrected fuzzy K-means method is used to correct the bias field before K-means thresholding. The proposed model can be minimized efficiently using the augmented Lagrangian algorithm. Experimental results and comparison analysis demonstrate that the proposed framework is not only able to preserve the geometry of object shapes, especially object corners, but it is also more accurate than state-of-the-art methods.

Appendix: The Derivations in Tables 1 and 2.

In this “Appendix”, we discuss the corresponding functionals for each variable function u, q, p, n, m by fixing the other variable functions and show how to get the associated minimizers or the Euler–Lagrange equations. Similarly as in [28], these functionals can be written as follows.

$$\begin{aligned} \varepsilon _1(u) =&\frac{1}{2}\int _{\Omega }(f-u)^2dx + \frac{\mu }{2}\int _{\Omega \Gamma }|\nabla u|^2dx + \frac{r_2}{2}\int _{\Omega }|p-<\nabla u,1>|^2dx\nonumber \\&+\lambda _2\int _{\Omega }(p-<\nabla u,1>)dx, \end{aligned}$$

(4.1)

$$\begin{aligned} \varepsilon _2(q) =\,&\lambda \int _{\Omega }|q|dx + \frac{r_3}{2}\int _{\Omega }(q-\partial _xn_1-\partial _yn_2)^2dx + \lambda _3\int _{\Omega }(q-\partial _xn_1-\partial _yn_2)dx, \end{aligned}$$

(4.2)

$$\begin{aligned} \varepsilon _3(p) =\,&r_1\int _{\Omega }(|p|-p\cdot m)dx + \lambda _1\int _{\Omega }(|p|-p\cdot m)dx + \frac{r_2}{2}\int _{\Omega }|p-<\nabla u,1>|^2dx\nonumber \\&+\lambda _2\int _{\Omega }(p-<\nabla u,1>)dx, \end{aligned}$$

(4.3)

$$\begin{aligned} \varepsilon _4(n)=\,&\frac{r_3}{2}\int _{\Omega } (q-\partial _xn_1-\partial _yn_2)^2dx + \lambda _3\int _{\Omega }(q-\partial _xn_1-\partial _yn_2)dx + \frac{r_4}{2}\int _{\Omega }|n-m|^2dx\nonumber \\&+\lambda _4\int _{\Omega }(n-m)dx, \end{aligned}$$

(4.4)

$$\begin{aligned} \varepsilon _5(m) =\,&r_1\int _{\Omega }(|p|-p\cdot m)dx + \lambda _1\int _{\Omega }(|p|-p\cdot m)dx\nonumber \\&+\frac{r_4}{2}\int _{\Omega }|n-m|^2dx + \lambda _4\int _{\Omega }(n-m)dx\nonumber \\&+ \delta _{R}(m). \end{aligned}$$

(4.5)

The minimizers of the functionals $\varepsilon _2(q)$, $\varepsilon _3(p)$, $\varepsilon _5(m)$ can be expressed explicitly, while the minimizers of functionals $\varepsilon _1(u)$, $\varepsilon _4(n)$ are determined by the associated Euler–Lagrange equations. For the sake of completeness of presentation, we present the details of how to minimize the subproblems here.

For the functional $\varepsilon _1(u)$, since

$$\begin{aligned} \varepsilon _1(u) =&\frac{1}{2}\int _{\Omega }(f-u)^2dx + \frac{\mu }{2}\int _{\Omega \setminus \Gamma }|\nabla u|^2dx + \frac{r_2}{2}\int _{\Omega }(\partial _xu-d_{1})^2dx\\&+\frac{r_2}{2}\int _{\Omega }(\partial _yu-d_{2})^2dx, \end{aligned}$$

where $(d_{1},d_{2},d_{3})=d=p + \frac{\lambda _2}{r_2}$, then one gets

$$\begin{aligned}{}[I + (\mu +r_2)\Delta ]u = f + r_2 \partial _x^{T}d_{1} +r_2\partial _{y}^{T}d_{2}, \end{aligned}$$

where I is the identity operator, $\partial _{x}^{T}$ is the tranapose of $\partial _{x}$ and $\partial _{y}^{T}$ is the tranapose of $\partial _{y}$

As the functional $\varepsilon _2(q)$ can be reformulated as

$$\begin{aligned} \varepsilon _2(q) = \lambda \int _{\Omega }|q|dx + \frac{r_3}{2}\int _{\Omega }(q-\tilde{q})^2dx, \end{aligned}$$

where $\tilde{q} = \partial _{x}n_1+\partial _{y}n_2-\frac{\lambda _3}{r_3}$, then by Lemma 1 in [28], one gets

$$\begin{aligned} \arg \min \limits _{q}\varepsilon _2(q) = \max \big (|\tilde{q}|-\frac{\lambda }{r_3}, 0\big )\frac{\tilde{q}}{|\tilde{q}|}. \end{aligned}$$

Similarly, as

$$\begin{aligned} \varepsilon _3(p) =&(r_1+\lambda _1)\int _{\Omega }|p|dx + \frac{r_2}{2}\int _{\Omega }(p-\tilde{p})^2dx,\\ \varepsilon _5(m) =&\int _{\Omega }(m-\tilde{m})^2dx + \delta _{R}(m), \end{aligned}$$

where $\tilde{p} = <\nabla u,1> - \frac{\lambda _2}{r_2} + \frac{m(r_1+\lambda _1)}{r_2}$,  $\tilde{m} = n+\frac{\lambda _4}{r_4}+\frac{(r_1+\lambda _1)p}{r_4}$. By using Lemma 1 in [28], one obtains

$$\begin{aligned} \arg \min \limits _{p}\varepsilon _3(p) = \max \big (|\tilde{p}|-\frac{r_1+\lambda _1}{r_2}, 0\big )\frac{\tilde{p}}{|\tilde{p}|}. \end{aligned}$$

By using Lemma 2 in [28], we get

$$\begin{aligned}&\arg \min \limits _{m}\varepsilon _5(m) = {\left\{ \begin{array}{ll} &{}\tilde{m}, ~~~~~|\tilde{m}|\le 1,\\ &{}\frac{\tilde{m}}{|\tilde{m}|}, ~~~ |\tilde{m}|>1. \end{array}\right. }\ \end{aligned}$$

As for the functional $\varepsilon _4(n)$, standard procedures lead to the following Euler–Lagrange equations:

$$\begin{aligned} \Big (\frac{r_4}{r_3}+\partial _x^{T}\partial _x\Big )n_1 + \partial _x^{T}\partial _yn_2&= \frac{r_4}{r_3}h_1+\partial _x^{T}q+\frac{1}{r_3}\partial _x^{T}\lambda _3,\\ \Big (\frac{r_4}{r_3}+\partial _y^{T}\partial _y\Big )n_2 + \partial _{y}^{T}\partial _{x}n_1&= \frac{r_4}{r_3}h_2+\partial _{y}^{T}q+\frac{1}{r_3}\partial _y^{T}\lambda _3,\\ n_3&=h_3. \end{aligned}$$

where $n=(n_1,n_2,n_3)$, $m=(m_1,m_2,m_3)$, $\lambda _2=(\lambda _{21},\lambda _{22},\lambda _{23})$, $\lambda _4=(\lambda _{41},\lambda _{42},\lambda _{43})$, $h=(h_1,h_2,h_3)=m-\frac{\lambda _4}{r_4}$. Both u and n can be solved by the fast Fourier transform [28]. Moreover, based on the above formulations, we may update all the Lagrange multipliers:

$$\begin{aligned} \lambda _1^{new} =\,&\lambda _1^{old} + r_1(|p^{old}|-p^{old}\cdot m^{old}),\\ \lambda _{21}^{new} =\,&\lambda _{21}^{old} + r_2(p_{1}^{old}-\partial _x u^{old}),\\ \lambda _{22}^{new} =\,&\lambda _{22}^{old} + r_2(p_{2}^{old}-\partial _y u^{old}),\\ \lambda _{23}^{new} =\,&\lambda _{23}^{old} + r_2(p_{3}^{old}-1),\\ \lambda _3^{new} =\,&\lambda _3^{old} + r_3(q^{old}-\partial _xn_{1}^{old}-\partial _yn_{2}^{old}),\\ \lambda _{41}^{new} =\,&\lambda _{41}^{old} + r_4(n_{1}^{old}-m_{1}^{old}),\\ \lambda _{42}^{new} =\,&\lambda _{42}^{old} + r_4(n_{2}^{old}-m_{2}^{old}),\\ \lambda _{43}^{new} =\,&\lambda _{43}^{old} + r_4(n_{3}^{old}-m_{3}^{old}), \end{aligned}$$

where $|p|=\sqrt{(p_{1})^2+(p_{2})^2+(p_{3})^2}$.