(mumford-shah-model)=
# Mumford-Shah model

In the following we will discuss a seminal work on a variational method
for image segmentation by {cite:p}`mumford_1989`.

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
```{math}
\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:
```{math}
\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.

````{prf:remark} 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.,
    ```{math}
    \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:
```{math}
\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
{ref}`s:variational_denoising`. 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:
```{math}
\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).