7.1. Mumford-Shah model

In the following we will discuss a seminal work on a variational method for image segmentation by [MS89].

Let \(\Omega \subset \R^d\) be an open, bounded, and connected image domain and let \(f \colon \Omega \rightarrow \R\) be a a given (possibly noisy) image, for which we may assume that \(f \in L^\infty(\Omega)\). The main idea is to partition the image domain \(\Omega\) via a segmentation boundary \(\Gamma \coloneqq \cup_{i=1}^m \partial \Omega_i \setminus \partial \Omega \subset \Omega\) with

\[\Omega \ = \ \Omega_1 \cup \dots \cup \Omega_m \cup \Gamma\]

such that

1. \(f: \Omega \to \R\) varies smoothly within each subregion \(\Omega_i \subset \Omega\)

2. \(f\) has intensity jumps only across the segmentation boundary \(\Gamma\).

The Mumford-Shah energy functional \(E(\Gamma, u)\) models these requirements for image segmentation and the related variational problem is defined for fixed regularization parameters \(\alpha, \beta > 0\) as:

\[\min_{\Gamma \subset \Omega, u \in X} E(\Gamma, u) \ = \ \frac{1}{2} \int_\Omega (f(x)-u(x))^2 \, \mathrm{d}x + \frac{\alpha}{2} \int_{\Omega \setminus \Gamma} \| \nabla u(x) \|^2 \, \mathrm{d}x + \beta |\Gamma|\]

The above functional penalizes the distance between the approximated image \(u\) and the input image \(f\), the lack of smoothness of the approximation within the subregions, and the length of the segmentation boundary \(\Gamma\). This Mumford-Shah model can be classified as a free boundary problem, which is very challenging to solve.

Remark 7.1 (Properties of the Mumford-Shah model)

1. Dropping any term in the energy functional \(E\) leads to trivial solutions and \(\min_{u,\Gamma} E(\Gamma, u) = 0\), due to the following reasons:

   • if one drops the data fidelity term, a trivial solution is \(u \equiv 0\) and \(\Gamma = \emptyset\),

   • if one drops the regularization term, i.e., \(\alpha = 0\), a trivial solution is \(u = f\) and \(\Gamma = \emptyset\),

   • if one drops the geometric prior on the segmentation boundary \(\Gamma\), i.e., \(\beta = 0\), one can compute in the discrete setting a trivial solution for \(u = f\) and setting the segmentation boundary \(\Gamma\) as border between all pixels.

2. Nontrivial solutions \((\hat{u}, \hat{\Gamma})\) of the Mumford-Shah model typically lead to cartoon-like approximations of \(f\) with \(\Gamma\) at the boundary of image regions at high-contrast edges.

3. As we are interested in solutions \(\hat{u}\) which are smooth within each subregion \(\Omega_i \subset \Omega\) of the image domain, but should have jumps across the segmentation boundary \(\Gamma\), a suitable function space for optimization is \(W^{1,2}(\Omega \setminus \Gamma) = H^1(\Omega \setminus \Gamma)\). Note that choosing \(H^1(\Omega)\) for optimization would lead to globally smooth solutions (as in the case of variational denoising) since the segmentation boundary \(\Gamma\) is a null set with respect to the Lebesgue measure, i.e.,

   \[\int_{\Omega \setminus \Gamma} ||\nabla u(x)||^2 \, \mathrm{d}x \ = \ \int_{\Omega} ||\nabla u(x)||^2 \, \mathrm{d}x.\]

4. For \(\alpha \rightarrow \infty\) the optimal approximation \(\hat{u}\) converges to a piecewise constant solution \(\hat{u}_c\), which is constant in each subregion \(\Omega_i \subset \Omega\). Simultaneously, the segmentation boundary \(\Gamma\) becomes the edge set of \(\hat{u}_c\).

Minimizing the Mumford-Shah functional \(E(\Gamma, u)\) in both variables turns out to be a nonconvex problem. We alleviate this problem by keeping the boundary \(\Gamma \subset\Omega\) fixed and optimize for \(u \in W^{1,2}(\Omega \setminus \Gamma) = H^1(\Omega \setminus \Gamma)\), which leads to the following form of our problem:

\[\min_{u \in H^1(\Omega \setminus \Gamma)} E(u) \ = \ \frac{1}{2} \int_\Omega (f(x) - u(x))^2 \, \mathrm{d}x + \frac{\alpha}{2}\int_{\Omega \setminus \Gamma} \| \nabla u \|^2 \, \mathrm{d}x.\]

The energy functional \(E(u)\) is strictly convex and differentiable in \(H^1(\Omega \setminus \Gamma)\) and hence the solution to this problem is unique. Note that it is very similar to the variational denoising problem with generalized Tykonov regularization discussed in Variational denoising problem. The only difference is that the regularization functional has to respect the segmentation boundary \(\Gamma\) to compute piecewise-smooth approximations. Computing the Fréchet derivative of this energy functional leads to the necessary optimality conditions:

\[\begin{aligned} \Delta u(x) \ &= \ \frac{1}{\alpha}(u(x) - f(x)), \qquad \forall x \in \Omega\setminus\Gamma, \\ \frac{\partial u}{\partial \vec{n}}(x) \ &= \ 0, \qquad \forall x \in \partial ( \Omega\setminus\Gamma) = \partial \Omega\cup\Gamma. \end{aligned}\]

Although the above problem for a fixed segmentation boundary \(\Gamma\) is a linear elliptic PDE and not a free‑boundary problem, it still poses significant numerical challenges. The domain \(\Omega \setminus \Gamma\) is a bounded region cut by an internal segmentation boundary \(\Gamma\) that may be highly irregular, contain corners or junctions, and partition the domain into many separate connected components. Using a numerical grid or conforming finite element meshes must resolve \(\Gamma\) exactly, which becomes difficult for complex geometry and leads to large, ill‑conditioned linear systems. The Neumann conditions on both the outer boundary and along the boundary mean that the solution can develop steep gradients near \(\Gamma\), requiring fine local resolution. Moreover, in iterative minimisation schemes the segmentation boundary \(\Gamma\) changes at every outer iteration, demanding repeated remeshing or the use of sophisticated unfitted discretisations (e.g., cut‑cell or extended finite element methods).