9  Image Segmentation as Optimization

Segmentation is the task of dividing an image into meaningful regions. In this chapter, segmentation is developed as a mathematical problem in estimation and optimization. We begin with thresholding and clustering, then move to graph-based energies, spectral relaxations, variational models, and neural segmentation. The central theme is simple but powerful: a segmentation is a discrete or continuous label field chosen to balance fidelity to image measurements with spatial regularity.

9.1 Learning Goals

After reading this chapter, students should be able to:

1. define binary, multiclass, semantic, and instance segmentation mathematically;
2. represent a segmentation as a label function, indicator function, partition, or probability field;
3. explain thresholding as a one-dimensional decision rule;
4. derive Otsu-style thresholding from within-class variance minimization;
5. formulate color or feature segmentation using \(k\)-means and mixture models;
6. write segmentation as an energy minimization problem with data and smoothness terms;
7. explain how binary graph cuts solve certain discrete segmentation energies exactly;
8. describe spectral segmentation and the normalized-cut relaxation;
9. connect variational segmentation with perimeter, total variation, and level-set ideas;
10. compute and interpret IoU, Dice, precision, recall, and pixel accuracy;
11. understand how U-Net-style neural segmentation models estimate label fields;
12. design Python experiments comparing classical, statistical, and neural segmentation methods.

9.2 8.1 What Is Image Segmentation?

A digital image is an array

\[ I \in \mathbb{R}^{m\times n\times C}, \]

where \(m\) is the number of rows, \(n\) is the number of columns, and \(C\) is the number of channels. A segmentation assigns labels to pixels. For \(K\) classes, a segmentation is a map

\[ y:\{1,\dots,m\}\times \{1,\dots,n\}\to \{1,\dots,K\}. \]

For binary segmentation, \(K=2\), and it is often convenient to use an indicator function

\[ u_{ij}\in \{0,1\}, \]

where \(u_{ij}=1\) means foreground and \(u_{ij}=0\) means background.

Segmentation can also be represented by regions. If

\[ \Omega = \{1,\dots,m\}\times \{1,\dots,n\}, \]

then a \(K\)-region segmentation is a partition

\[ \Omega = \Omega_1\cup \cdots \cup \Omega_K, \qquad \Omega_a\cap \Omega_b = \emptyset \quad \text{for } a\ne b. \]

The label field \(y\) and the partition \(\{\Omega_k\}\) contain the same information:

\[ \Omega_k = \{(i,j): y_{ij}=k\}. \]

In modern neural segmentation, models often output probabilities rather than hard labels. A soft segmentation is a tensor

\[ p_{ij} \in \Delta^{K-1}, \]

where

\[ \Delta^{K-1} = \left\{p\in \mathbb{R}^K: p_k\ge 0,\; \sum_{k=1}^K p_k=1\right\}. \]

The hard label is then usually obtained by

\[ \widehat y_{ij}=\arg\max_{1\le k\le K} p_{ij,k}. \]

A segmentation is a spatially organized decision rule. Unlike ordinary classification, the decisions are not independent: neighboring pixels are usually expected to have related labels.

9.3 8.2 Four Kinds of Segmentation

Computer vision uses the word segmentation in several related ways.

Binary segmentation separates foreground from background:

\[ y_{ij}\in \{0,1\}. \]

Examples include separating a handwritten digit from the paper, a tumor from healthy tissue, or a road from the rest of a scene.

Multiclass semantic segmentation assigns each pixel a semantic class:

\[ y_{ij}\in \{\text{road},\text{sky},\text{building},\text{car},\dots\}. \]

Pixels belonging to different cars may share the same class label ``car.’’

Instance segmentation separates object instances. If two cars appear in the image, they should receive different instance identities. The output is no longer just a semantic label field; it also contains object identity.

Panoptic segmentation combines semantic segmentation and instance segmentation. Every pixel gets both a semantic category and, when appropriate, an instance identity.

This chapter focuses on the mathematical foundations that support all of these tasks: label functions, energies, probability models, and spatial regularization. Later chapters return to object detection, semantic segmentation, instance segmentation, promptable segmentation, and foundation models.

9.4 8.3 Segmentation as an Inverse Problem

Segmentation is an inverse problem. We observe an image \(I\), and we want to recover an unobserved label field \(y\) or region indicator \(u\). A general statistical formulation is

\[ \widehat y = \arg\max_y p(y\mid I). \]

By Bayes’ rule,

\[ p(y\mid I) \propto p(I\mid y)p(y). \]

Taking negative logarithms gives the optimization problem

\[ \widehat y = \arg\min_y \left\{-\log p(I\mid y)-\log p(y)\right\}. \]

The term \(-\log p(I\mid y)\) measures how well the label field explains the observed image. The term \(-\log p(y)\) encodes prior beliefs, such as spatial smoothness or preference for short boundaries.

A common energy form is

\[ E(y) = \sum_{p\in \Omega} D_p(y_p) + \lambda \sum_{(p,q)\in \mathcal{E}} V_{pq}(y_p,y_q), \]

where \(p\) and \(q\) index pixels, \(D_p\) is a data cost, \(V_{pq}\) is a smoothness cost, \(\mathcal{E}\) is a set of neighboring pixel pairs, and \(\lambda\ge 0\) controls regularization.

This decomposition appears repeatedly in computer vision. The notation changes across methods, but the concept remains stable:

\[ \text{segmentation energy} = \text{data fidelity} + \text{spatial regularity}. \]

9.5 8.4 Thresholding: The Simplest Segmentation Rule

Suppose \(I\) is grayscale and we want a binary segmentation. The simplest rule is thresholding:

\[ u_{ij} = \begin{cases} 1, & I_{ij}>t,\\ 0, & I_{ij}\le t. \end{cases} \]

Equivalently,

\[ u_{ij}=\mathbf{1}\{I_{ij}>t\}. \]

This is a one-dimensional classifier applied independently to each pixel. It works well when foreground and background intensities are clearly separated.

Thresholding can be understood as a decision problem. Suppose the two classes have conditional densities

\[ I_{ij}\mid y_{ij}=0 \sim f_0, \qquad I_{ij}\mid y_{ij}=1 \sim f_1. \]

If the class priors are \(\pi_0\) and \(\pi_1\), then the Bayes rule chooses foreground when

\[ \pi_1 f_1(I_{ij}) > \pi_0 f_0(I_{ij}). \]

A threshold rule is optimal when this likelihood-ratio decision boundary occurs at a single intensity value.

9.5.1 Interactive: thresholding a grayscale image

This demo constructs a synthetic grayscale image and compares several thresholds. The threshold is a scalar parameter, but it produces a spatial label field.

9.6 8.5 Otsu Thresholding and Within-Class Variance

Otsu’s method chooses a threshold using the image histogram (Otsu 1979). The method assumes that the histogram is approximately bimodal and chooses a threshold that separates the two intensity populations.

Let intensities take discrete values \(0,1,\dots,L-1\). Let

\[ p_\ell = \frac{\#\{(i,j): I_{ij}=\ell\}}{mn} \]

be the normalized histogram. For a threshold \(t\), define two classes

\[ C_0(t)=\{0,1,\dots,t\}, \qquad C_1(t)=\{t+1,\dots,L-1\}. \]

The class weights are

\[ w_0(t)=\sum_{\ell=0}^t p_\ell, \qquad w_1(t)=\sum_{\ell=t+1}^{L-1}p_\ell. \]

The class means are

\[ \mu_0(t)=\frac{1}{w_0(t)}\sum_{\ell=0}^t \ell p_\ell, \qquad \mu_1(t)=\frac{1}{w_1(t)}\sum_{\ell=t+1}^{L-1}\ell p_\ell. \]

The within-class variance is

\[ \sigma_W^2(t) = w_0(t)\sigma_0^2(t)+w_1(t)\sigma_1^2(t). \]

Otsu’s threshold minimizes \(\sigma_W^2(t)\). Equivalently, because the total variance decomposes into within-class and between-class parts, it maximizes the between-class variance

\[ \sigma_B^2(t) = w_0(t)w_1(t)\left(\mu_0(t)-\mu_1(t)\right)^2. \]

Thus Otsu thresholding is not merely a heuristic. It is a variance decomposition method for one-dimensional unsupervised classification.

NoteStatistical interpretation

Otsu’s method chooses the threshold that makes the two induced intensity classes as internally concentrated and externally separated as possible. For students in statistics, this is closely related to one-dimensional clustering and ANOVA-style variance decomposition.

9.7 8.6 Multilevel Thresholding

For \(K>2\) intensity classes, thresholding uses \(K-1\) thresholds

\[ t_1<t_2<\cdots<t_{K-1}. \]

The label rule is

\[ y_{ij}=k \quad \text{if} \quad t_{k-1}< I_{ij}\le t_k, \]

with \(t_0=-\infty\) and \(t_K=\infty\). The same variance idea can be extended by minimizing the sum of within-class variances:

\[ \sum_{k=1}^K w_k \sigma_k^2. \]

However, pure intensity thresholding ignores spatial relationships and texture. Two visually distinct regions may share the same intensity range. Conversely, one object may contain shadows, highlights, and texture that produce a wide intensity distribution.

This limitation motivates feature-based segmentation.

9.8 8.7 Feature-Based Segmentation

Instead of segmenting using intensity alone, we can construct a feature vector for each pixel:

\[ z_p \in \mathbb{R}^d, \qquad p\in \Omega. \]

Features may include:

• color values such as RGB or Lab;
• image coordinates \((i,j)\);
• texture responses from filters;
• gradient magnitude and orientation;
• local statistics such as mean and variance in a window;
• learned embeddings from a neural network.

Then segmentation becomes clustering or classification in feature space.

For example, a simple feature vector may be

\[ z_p = \begin{bmatrix} R_p\\G_p\\B_p\\ \alpha i_p\\ \alpha j_p \end{bmatrix}, \]

where the parameter \(\alpha\) controls how strongly spatial proximity influences the clustering.

The inclusion of coordinates is important. Without coordinates, two far-away regions with similar colors may be grouped together. With coordinates, the algorithm is encouraged to form spatially coherent regions.

9.9 8.8 \(k\)-Means Segmentation

Given feature vectors \(z_1,\dots,z_N\in \mathbb{R}^d\), where \(N=mn\), \(k\)-means solves

\[ \min_{c_1,\dots,c_K,\; y_1,\dots,y_N} \sum_{p=1}^N \|z_p-c_{y_p}\|_2^2, \]

where \(c_k\in \mathbb{R}^d\) is the centroid of cluster \(k\) and \(y_p\in \{1,\dots,K\}\) is the cluster assignment for pixel \(p\).

For fixed centroids, the optimal label is

\[ y_p = \arg\min_{1\le k\le K}\|z_p-c_k\|_2^2. \]

For fixed labels, the optimal centroid is

\[ c_k = \frac{1}{|C_k|}\sum_{p:y_p=k} z_p. \]

The standard \(k\)-means algorithm alternates between these two steps. Each step decreases the objective, so the algorithm converges to a local minimum.

9.9.1 Interactive: feature-space clustering for segmentation

This demo shows a synthetic color/spatial feature cloud and its cluster assignment. Changing the feature representation changes the implied segmentation geometry.

9.10 8.9 Gaussian Mixture Models and Soft Segmentation

\(k\)-means produces hard clusters. A probabilistic alternative is the Gaussian mixture model:

\[ p(z) = \sum_{k=1}^K \pi_k \mathcal{N}(z\mid \mu_k,\Sigma_k), \]

where \(\pi_k\ge 0\) and \(\sum_k \pi_k=1\).

The posterior probability that pixel \(p\) belongs to class \(k\) is

\[ \gamma_{pk} = P(y_p=k\mid z_p) = \frac{\pi_k \mathcal{N}(z_p\mid \mu_k,\Sigma_k)} {\sum_{\ell=1}^K \pi_\ell \mathcal{N}(z_p\mid \mu_\ell,\Sigma_\ell)}. \]

The soft segmentation is the matrix of posterior probabilities \(\gamma_{pk}\). The hard segmentation is

\[ \widehat y_p=\arg\max_k \gamma_{pk}. \]

Parameters can be estimated by the expectation-maximization algorithm. The E-step computes the responsibilities \(\gamma_{pk}\). The M-step updates the parameters:

\[ N_k=\sum_{p=1}^N \gamma_{pk}, \]

\[ \pi_k=\frac{N_k}{N}, \qquad \mu_k=\frac{1}{N_k}\sum_{p=1}^N \gamma_{pk}z_p, \]

and

\[ \Sigma_k = \frac{1}{N_k}\sum_{p=1}^N \gamma_{pk}(z_p-\mu_k)(z_p-\mu_k)^T. \]

This probabilistic viewpoint is useful because segmentation uncertainty matters. Boundary pixels, shadows, low contrast regions, and ambiguous textures often should not be forced into overconfident labels too early.

9.11 8.10 Why Clustering Alone Is Not Enough

Thresholding and clustering treat pixels as independent feature vectors. However, images are spatial objects. Neighboring pixels are often related. A cluster assignment may be noisy if it ignores the image grid.

A segmentation method should balance two goals:

1. data fidelity: assign labels that agree with the observed features;
2. spatial regularity: prefer simple, coherent regions unless the data strongly suggest a boundary.

This leads to discrete energy models.

9.12 8.11 Discrete Segmentation Energies

Let \(G=(V,E)\) be a graph whose vertices are pixels and whose edges connect neighboring pixels. Usually \(E\) connects 4-neighbors or 8-neighbors on the image grid.

A common binary segmentation energy is

\[ E(u) = \sum_{p\in V}D_p(u_p) + \lambda\sum_{(p,q)\in E}w_{pq}|u_p-u_q|, \]

where \(u_p\in \{0,1\}\). The data term \(D_p(u_p)\) measures the cost of assigning label \(u_p\) to pixel \(p\). The pairwise term penalizes label changes between neighboring pixels.

The weight \(w_{pq}\) can be chosen to preserve image edges. A common choice is

\[ w_{pq} = \exp\left(-\frac{\|I_p-I_q\|^2}{2\sigma^2}\right). \]

If neighboring pixels have similar intensities, then \(w_{pq}\) is large, so cutting between them is expensive. If neighboring pixels have very different intensities, then \(w_{pq}\) is small, so a boundary is allowed.

This model encodes the following principle:

\[ \text{good boundaries should pass through weak similarity links.} \]

9.13 8.12 Maximum A Posteriori Interpretation

The same energy arises from a Markov random field prior. Suppose

\[ p(u\mid I) \propto \exp\left(-E(u)\right). \]

Then minimizing \(E(u)\) is equivalent to maximum a posteriori estimation:

\[ \widehat u = \arg\max_u p(u\mid I). \]

The unary terms act like pixelwise likelihoods. The pairwise terms act like a prior favoring spatial smoothness.

For binary labels, one may set

\[ D_p(1)=-\log P(I_p\mid u_p=1)-\log P(u_p=1), \]

and

\[ D_p(0)=-\log P(I_p\mid u_p=0)-\log P(u_p=0). \]

The pairwise Potts penalty

\[ V_{pq}(u_p,u_q)=w_{pq}\mathbf{1}\{u_p\ne u_q\} \]

encourages neighboring pixels to have equal labels.

9.14 8.13 Graph Cuts for Binary Segmentation

Certain binary segmentation energies can be minimized exactly by an \(s\)-\(t\) graph cut (Boykov et al. 2001; Kolmogorov and Zabih 2004). The idea is to construct a graph with two special nodes: a source \(s\) representing foreground and a sink \(t\) representing background.

Each pixel is connected to the source and sink using capacities derived from unary costs. Neighboring pixels are connected to each other using capacities derived from pairwise smoothness weights.

An \(s\)-\(t\) cut separates the graph into two parts, one containing \(s\) and one containing \(t\). A pixel assigned to the source side receives foreground label. A pixel assigned to the sink side receives background label.

For the energy

\[ E(u) = \sum_p D_p(u_p)+\lambda\sum_{(p,q)} w_{pq}|u_p-u_q|, \]

the cost of the cut can be made equal to \(E(u)\) up to an additive constant. Therefore the minimum cut gives the globally optimal binary segmentation.

9.14.1 Interactive: graph-cut intuition

This demo uses a small grid with unary foreground/background evidence and pairwise smoothing. The display shows how a boundary prefers to pass through weak similarity edges.

ImportantExactness is special

Graph cuts are powerful because they provide exact global minimizers for important binary energies. However, not every discrete segmentation energy can be solved exactly by a graph cut. The pairwise terms must satisfy appropriate submodularity conditions.

9.15 8.14 Submodularity in Binary Pairwise Energies

Consider a binary pairwise energy term

\[ V_{pq}(u_p,u_q), \qquad u_p,u_q\in \{0,1\}. \]

For graph-cut minimization, a key condition is submodularity:

\[ V_{pq}(0,0)+V_{pq}(1,1) \le V_{pq}(0,1)+V_{pq}(1,0). \]

This says that equal labels should not be more expensive than unequal labels. The Potts penalty

\[ V_{pq}(u_p,u_q)=w_{pq}\mathbf{1}\{u_p\ne u_q\} \]

satisfies this condition when \(w_{pq}\ge 0\), since

\[ 0+0\le w_{pq}+w_{pq}. \]

For multiclass segmentation, exact minimization is generally harder. Approximate methods such as alpha-expansion can be used for certain metric pairwise costs.

9.16 8.15 Spectral Segmentation and Normalized Cuts

Graph cuts minimize a discrete cut energy with unary terms. Spectral methods focus on graph partitioning using the graph Laplacian. Let \(W\) be a symmetric affinity matrix, where \(W_{pq}\) measures similarity between pixels or superpixels. Let

\[ D_{pp}=\sum_q W_{pq} \]

be the degree matrix. The unnormalized graph Laplacian is

\[ L=D-W. \]

For a vector \(f\in \mathbb{R}^N\),

\[ f^T L f = \frac{1}{2}\sum_{p,q}W_{pq}(f_p-f_q)^2. \]

Thus the Laplacian quadratic form penalizes label variation across strong similarity edges.

In normalized cuts (Shi and Malik 2000), the goal is to partition the graph into sets \(A\) and \(B\) while accounting for the total connection volume of each set. The normalized cut objective is

\[ \operatorname{Ncut}(A,B) = \frac{\operatorname{cut}(A,B)}{\operatorname{assoc}(A,V)} + \frac{\operatorname{cut}(A,B)}{\operatorname{assoc}(B,V)}, \]

where

\[ \operatorname{cut}(A,B)=\sum_{p\in A,q\in B}W_{pq} \]

and

\[ \operatorname{assoc}(A,V)=\sum_{p\in A,q\in V}W_{pq}. \]

The exact discrete problem is difficult, but a continuous relaxation leads to a generalized eigenvalue problem involving \(L\) and \(D\):

\[ Lz = \lambda Dz. \]

The eigenvectors provide low-dimensional coordinates for graph partitioning.

Graph-cut segmentation and spectral segmentation both use graphs, but they optimize different objectives. Graph cuts are often used with unary data terms and can solve certain binary energies exactly. Spectral methods use eigenvectors of graph Laplacians to reveal global grouping structure.

9.17 8.16 Variational Segmentation

The discrete graph model has a continuous analogue. Let \(\Omega\subset \mathbb{R}^2\) be an image domain and let \(u:\Omega\to \{0,1\}\) be a binary segmentation. A variational segmentation energy has the form

\[ E(u)=\int_\Omega D(x,u(x))\,dx + \lambda \operatorname{Length}(\partial \{x:u(x)=1\}). \]

The first term is data fidelity. The second term penalizes boundary length. This model prefers regions that explain the data with short, smooth boundaries.

A classical model in this direction is the Mumford–Shah functional (Mumford and Shah 1989). It seeks a piecewise-smooth approximation \(v\) and an edge set \(K\) by minimizing

\[ E(v,K) = \int_{\Omega\setminus K} |I(x)-v(x)|^2\,dx + \alpha\int_{\Omega\setminus K}|\nabla v(x)|^2\,dx + \beta\operatorname{Length}(K). \]

Here \(v\) approximates the image away from edges, \(K\) represents the edge set, and the last term penalizes edge length.

9.18 8.17 Chan–Vese Segmentation

The Chan–Vese model (Chan and Vese 2001) is a region-based simplification of Mumford–Shah. For a binary region \(R\subset \Omega\), it minimizes

\[ E(R,c_1,c_0) = \int_R |I(x)-c_1|^2\,dx + \int_{\Omega\setminus R}|I(x)-c_0|^2\,dx + \lambda \operatorname{Length}(\partial R), \]

where \(c_1\) and \(c_0\) are the average intensities inside and outside the region.

For fixed \(R\), the optimal constants are

\[ c_1 = \frac{1}{|R|}\int_R I(x)\,dx, \qquad c_0 = \frac{1}{|\Omega\setminus R|}\int_{\Omega\setminus R} I(x)\,dx. \]

The model is powerful because it can segment objects even when edges are weak, provided that the inside and outside regions have different average intensities.

9.19 8.18 Level-Set Representation

Instead of representing a region \(R\) directly, one may represent it as the positive set of a function \(\phi:\Omega\to \mathbb{R}\):

\[ R=\{x\in \Omega: \phi(x)>0\}. \]

The boundary is the zero level set

\[ \partial R=\{x\in \Omega: \phi(x)=0\}. \]

Level-set methods allow the boundary to move, split, and merge naturally. This is useful for topology changes during segmentation.

The binary indicator can be written using a Heaviside function:

\[ u(x)=H(\phi(x)). \]

The boundary length can be expressed formally as

\[ \operatorname{Length}(\partial R) = \int_\Omega \delta(\phi(x))|\nabla \phi(x)|\,dx, \]

where \(\delta\) is the Dirac delta distribution. In numerical algorithms, smooth approximations of \(H\) and \(\delta\) are used.

9.19.1 Interactive: region energy and boundary regularization

This demo compares a noisy data term with a smoothed boundary preference. Increasing regularity suppresses isolated noisy label decisions.

9.20 8.19 Total Variation and Convex Relaxation

For binary segmentation, the perimeter of a region is closely related to the total variation of its indicator function:

\[ \operatorname{TV}(u) = \int_\Omega |\nabla u(x)|\,dx. \]

When \(u\) is an indicator function, \(\operatorname{TV}(u)\) corresponds to boundary length in a generalized sense. This motivates relaxed energies of the form

\[ E(u) = \int_\Omega D_1(x)u(x)+D_0(x)(1-u(x))\,dx + \lambda \operatorname{TV}(u), \]

where now

\[ 0\le u(x)\le 1. \]

The constraint has been relaxed from binary values to the interval \([0,1]\). Under suitable conditions, thresholding the relaxed solution gives a binary segmentation. This connection is one reason total-variation methods are important in mathematical imaging.

9.21 8.20 Segmentation Losses and Evaluation Metrics

For supervised segmentation, we compare a predicted mask with a ground-truth mask. Let \(A\) be the predicted foreground set and \(B\) the ground-truth foreground set.

The intersection-over-union score is

\[ \operatorname{IoU}(A,B) = \frac{|A\cap B|}{|A\cup B|}. \]

The Dice coefficient is

\[ \operatorname{Dice}(A,B) = \frac{2|A\cap B|}{|A|+|B|}. \]

For binary masks, these are related by

\[ \operatorname{Dice} = \frac{2\operatorname{IoU}}{1+\operatorname{IoU}}, \qquad \operatorname{IoU} = \frac{\operatorname{Dice}}{2-\operatorname{Dice}}. \]

Pixel accuracy is

\[ \frac{\operatorname{TP}+\operatorname{TN}} {\operatorname{TP}+\operatorname{TN}+\operatorname{FP}+\operatorname{FN}}. \]

However, pixel accuracy can be misleading when foreground is rare. In medical imaging, for example, a model may achieve high pixel accuracy by predicting background everywhere.

9.21.1 Interactive: IoU, Dice, and threshold choice

This demo shows how changing a probability threshold changes the predicted mask and the resulting IoU and Dice values.

9.22 8.21 Cross-Entropy and Dice Loss

Neural segmentation models usually output a probability \(p_{ij}\in (0,1)\) for binary foreground. Given ground-truth labels \(y_{ij}\in \{0,1\}\), the binary cross-entropy loss is

\[ \operatorname{BCE}(y,p) = -\frac{1}{mn}\sum_{i,j} \left[y_{ij}\log p_{ij}+(1-y_{ij})\log(1-p_{ij})\right]. \]

For multiclass segmentation with probabilities \(p_{ij,k}\), the cross-entropy loss is

\[ \operatorname{CE}(y,p) = -\frac{1}{mn}\sum_{i,j}\log p_{ij,y_{ij}}. \]

Dice loss is based on the soft Dice score. For binary segmentation,

\[ \operatorname{SoftDice}(y,p) = \frac{2\sum_{i,j}y_{ij}p_{ij}+\epsilon} {\sum_{i,j}y_{ij}+\sum_{i,j}p_{ij}+\epsilon}, \]

and

\[ \operatorname{DiceLoss}=1-\operatorname{SoftDice}. \]

Dice loss directly addresses foreground overlap and is often useful when classes are imbalanced.

9.23 8.22 Neural Segmentation Models

A neural segmentation model is a function

\[ f_\theta:\mathbb{R}^{m\times n\times C}\to [0,1]^{m\times n\times K}. \]

Unlike image classification, segmentation requires an output at every pixel. Therefore the model must preserve or recover spatial resolution.

A fully convolutional network replaces dense classification layers with convolutional layers so that the output remains spatial (Long et al. 2015). An encoder-decoder network compresses the image into feature maps and then upsamples those features back to image resolution.

U-Net is a particularly influential encoder-decoder architecture (Ronneberger et al. 2015). It uses skip connections from the encoder to the decoder:

\[ \text{decoder feature at scale }s = \operatorname{concat}(\text{upsampled feature},\text{encoder feature at scale }s). \]

The skip connections help recover fine spatial details that may be lost during downsampling.

9.24 8.23 Segmentation as Dense Classification with Spatial Structure

At first glance, semantic segmentation is pixelwise classification. For each pixel \(p\), the model predicts

\[ p_\theta(y_p=k\mid I). \]

But this description is incomplete because \(p_\theta(y_p=k\mid I)\) depends on the whole image or a large receptive field. The label decision at one pixel is informed by context.

Classical segmentation used explicit regularizers such as graph cuts or total variation. Neural segmentation often learns spatial regularity implicitly from data. However, many modern architectures still contain mathematical structures that encourage spatial consistency:

• convolutional layers enforce local translation-equivariant processing;
• pooling and downsampling build multiscale context;
• skip connections preserve localization;
• attention layers capture long-range interactions;
• loss functions such as Dice emphasize region overlap;
• post-processing can use conditional random fields or connected components.

9.25 8.24 AI Component: Segmentation with Foundation Models

Modern segmentation increasingly uses large pretrained models. Promptable segmentation models take an image plus prompts such as points, boxes, or masks and return segmentation candidates. This changes the mathematical problem from

\[ I \mapsto y \]

to

\[ (I,\text{prompt})\mapsto y. \]

The prompt acts as side information or a constraint. For example, a positive point prompt says that a particular pixel belongs to the target object. A box prompt says that the target object lies inside a rectangular region.

From an optimization viewpoint, prompts modify the feasible set or data term. From a statistical viewpoint, prompts condition the posterior distribution:

\[ p(y\mid I,\text{prompt}). \]

This idea will be developed further in the chapter on promptable segmentation and foundation models.

9.26 8.25 AI-Assisted Mathematical Study

Large language models and multimodal AI tools can be useful for learning segmentation, but they should be used with mathematical verification.

Try the following activities.

1. Ask an AI assistant to explain the difference between thresholding, \(k\)-means segmentation, graph cuts, and U-Net segmentation. Then rewrite the answer as optimization problems.
2. Ask for a derivation of Otsu thresholding. Check whether the assistant correctly distinguishes within-class variance and between-class variance.
3. Ask why Dice loss is useful for class imbalance. Verify the answer using a synthetic example with rare foreground pixels.
4. Ask for pseudocode for binary graph-cut segmentation. Identify which parts require a max-flow/min-cut algorithm.
5. Ask an AI assistant to generate a Python experiment comparing thresholding and U-Net segmentation. Check whether the proposed evaluation metrics are appropriate.

The goal is not to outsource understanding. The goal is to use AI to create drafts, explanations, and tests that you then verify mathematically.

9.27 8.26 Python Example: Thresholding and Otsu’s Method

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(7243)

m, n = 160, 220
y, x = np.mgrid[0:m, 0:n]

background = 0.25 + 0.05 * rng.normal(size=(m, n))
object_mask = ((x - 115)**2 / 50**2 + (y - 80)**2 / 35**2) <= 1
image = background.copy()
image[object_mask] = 0.68 + 0.08 * rng.normal(size=object_mask.sum())
image = np.clip(image, 0, 1)

# Otsu threshold from scratch
hist, edges = np.histogram(image.ravel(), bins=256, range=(0, 1), density=False)
p = hist / hist.sum()
levels = (edges[:-1] + edges[1:]) / 2

omega = np.cumsum(p)
mu = np.cumsum(p * levels)
mu_total = mu[-1]

sigma_b2 = (mu_total * omega - mu)**2 / (omega * (1 - omega) + 1e-12)
idx = np.argmax(sigma_b2)
threshold = levels[idx]

pred = image > threshold

fig, axes = plt.subplots(1, 3, figsize=(11, 3.5))
axes[0].imshow(image, cmap="gray")
axes[0].set_title("synthetic image")
axes[1].plot(levels, p)
axes[1].axvline(threshold, color="black", linestyle="--")
axes[1].set_title(f"Otsu threshold = {threshold:.3f}")
axes[2].imshow(pred, cmap="gray")
axes[2].set_title("binary segmentation")
for ax in axes:
    ax.set_xticks([])
    ax.set_yticks([])
plt.tight_layout()
plt.show()

9.28 8.27 Python Example: \(k\)-Means Segmentation in Color and Space

import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import KMeans

rng = np.random.default_rng(7243)

m, n = 140, 180
y, x = np.mgrid[0:m, 0:n]

image = np.zeros((m, n, 3), dtype=float)
image[..., 0] = 0.20 + 0.10 * x / n
image[..., 1] = 0.35 + 0.08 * y / m
image[..., 2] = 0.55

mask1 = (x - 60)**2 + (y - 65)**2 < 32**2
mask2 = ((x - 125)**2 / 36**2 + (y - 78)**2 / 25**2) < 1
image[mask1] = np.array([0.90, 0.25, 0.20])
image[mask2] = np.array([0.20, 0.75, 0.35])
image += 0.04 * rng.normal(size=image.shape)
image = np.clip(image, 0, 1)

# Feature vector: RGB plus scaled coordinates
coords = np.stack([y / m, x / n], axis=-1)
alpha = 0.45
features = np.concatenate([image, alpha * coords], axis=-1).reshape(-1, 5)

km = KMeans(n_clusters=3, n_init=10, random_state=0)
labels = km.fit_predict(features).reshape(m, n)

fig, axes = plt.subplots(1, 2, figsize=(9, 4))
axes[0].imshow(image)
axes[0].set_title("input image")
axes[1].imshow(labels, cmap="tab10")
axes[1].set_title("k-means labels")
for ax in axes:
    ax.set_axis_off()
plt.tight_layout()
plt.show()

9.29 8.28 Python Example: Discrete Energy with Iterated Conditional Modes

The following example uses a simple coordinate-descent heuristic called iterated conditional modes. It is not a graph-cut solver. Its purpose is pedagogical: it shows how data terms and smoothness terms compete.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(7243)
m, n = 120, 160
y, x = np.mgrid[0:m, 0:n]
true_mask = ((x - 80)**2 / 42**2 + (y - 60)**2 / 28**2) <= 1

image = 0.30 + 0.35 * true_mask.astype(float) + 0.12 * rng.normal(size=(m, n))
image = np.clip(image, 0, 1)

# Unary costs from simple Gaussian class models.
mu0, mu1 = 0.30, 0.65
sigma = 0.13
D0 = 0.5 * ((image - mu0) / sigma)**2
D1 = 0.5 * ((image - mu1) / sigma)**2

u = (D1 < D0).astype(int)
lmbda = 0.9

for sweep in range(8):
    for i in range(1, m - 1):
        for j in range(1, n - 1):
            neighbors = [u[i-1, j], u[i+1, j], u[i, j-1], u[i, j+1]]
            cost0 = D0[i, j] + lmbda * sum(v != 0 for v in neighbors)
            cost1 = D1[i, j] + lmbda * sum(v != 1 for v in neighbors)
            u[i, j] = int(cost1 < cost0)

fig, axes = plt.subplots(1, 3, figsize=(11, 3.6))
axes[0].imshow(image, cmap="gray")
axes[0].set_title("image")
axes[1].imshow(D1 < D0, cmap="gray")
axes[1].set_title("unary-only")
axes[2].imshow(u, cmap="gray")
axes[2].set_title("with smoothness heuristic")
for ax in axes:
    ax.set_axis_off()
plt.tight_layout()
plt.show()

9.30 8.29 TensorFlow Example: A Tiny U-Net-Style Model

import tensorflow as tf

inputs = tf.keras.Input(shape=(128, 128, 3))

# Encoder
c1 = tf.keras.layers.Conv2D(32, 3, padding="same", activation="relu")(inputs)
c1 = tf.keras.layers.Conv2D(32, 3, padding="same", activation="relu")(c1)
p1 = tf.keras.layers.MaxPooling2D()(c1)

c2 = tf.keras.layers.Conv2D(64, 3, padding="same", activation="relu")(p1)
c2 = tf.keras.layers.Conv2D(64, 3, padding="same", activation="relu")(c2)
p2 = tf.keras.layers.MaxPooling2D()(c2)

# Bottleneck
b = tf.keras.layers.Conv2D(128, 3, padding="same", activation="relu")(p2)
b = tf.keras.layers.Conv2D(128, 3, padding="same", activation="relu")(b)

# Decoder with skip connections
u2 = tf.keras.layers.UpSampling2D()(b)
u2 = tf.keras.layers.Concatenate()([u2, c2])
c3 = tf.keras.layers.Conv2D(64, 3, padding="same", activation="relu")(u2)
c3 = tf.keras.layers.Conv2D(64, 3, padding="same", activation="relu")(c3)

u1 = tf.keras.layers.UpSampling2D()(c3)
u1 = tf.keras.layers.Concatenate()([u1, c1])
c4 = tf.keras.layers.Conv2D(32, 3, padding="same", activation="relu")(u1)
c4 = tf.keras.layers.Conv2D(32, 3, padding="same", activation="relu")(c4)

outputs = tf.keras.layers.Conv2D(1, 1, activation="sigmoid")(c4)

model = tf.keras.Model(inputs, outputs)
model.compile(
    optimizer="adam",
    loss="binary_crossentropy",
    metrics=[tf.keras.metrics.BinaryIoU(target_class_ids=[1], threshold=0.5)]
)

model.summary()

9.31 8.30 Conceptual Exercises

1. Explain why thresholding is a pixelwise classifier but not a spatial segmentation model.
2. Give an example where thresholding fails but \(k\)-means with color and coordinate features succeeds.
3. Explain the difference between semantic segmentation and instance segmentation.
4. Why can pixel accuracy be misleading for rare-object segmentation?
5. Explain why Dice loss is often useful for medical image segmentation.
6. Compare graph-cut smoothness and total-variation regularization.
7. Explain why U-Net skip connections are helpful for segmentation.
8. Describe how a point prompt changes the segmentation problem.

9.32 8.31 Mathematical Exercises

1. Let \(A\) and \(B\) be two finite masks. Prove that

   \[ \operatorname{Dice}(A,B) = \frac{2\operatorname{IoU}(A,B)}{1+\operatorname{IoU}(A,B)}. \]

2. Derive the optimal \(k\)-means centroid formula

   \[ c_k=\frac{1}{|C_k|}\sum_{p\in C_k}z_p. \]

3. Show that the graph Laplacian satisfies

   \[ f^TLf=\frac{1}{2}\sum_{p,q}W_{pq}(f_p-f_q)^2. \]

4. Verify that the binary Potts penalty

   \[ V(u_p,u_q)=w\mathbf{1}\{u_p\ne u_q\} \]

   satisfies the submodularity condition when \(w\ge 0\).

5. For the Chan–Vese energy with fixed region \(R\), derive the formulas for \(c_1\) and \(c_0\).

6. Suppose a segmentation model predicts probability \(p\) for every foreground pixel and probability \(q\) for every background pixel. Write the binary cross-entropy loss in terms of \(p\), \(q\), and the foreground fraction.

9.33 8.32 Python Exercises

1. Implement Otsu thresholding from scratch and compare it with a manually chosen threshold.
2. Segment a color image using \(k\)-means with RGB features only. Then add coordinate features and compare the result.
3. Implement IoU and Dice for binary masks. Test your functions on synthetic masks.
4. Create a noisy binary image and compare unary-only segmentation with a smoothed segmentation heuristic.
5. Train a tiny U-Net on synthetic circles, rectangles, and ellipses. Evaluate Dice on a held-out test set.
6. Compare binary cross-entropy and Dice loss when foreground occupies only 1 percent of the image.
7. Use a pretrained segmentation model and evaluate how input resizing affects boundary accuracy.

9.34 8.33 AI-Assisted Exercises

1. Ask an AI assistant to explain graph cuts for segmentation. Convert the explanation into the energy notation used in this chapter.
2. Ask for an implementation of Otsu thresholding. Check whether the implementation uses normalized histogram probabilities correctly.
3. Ask for a U-Net architecture. Identify the encoder, decoder, bottleneck, skip connections, and output layer.
4. Ask an AI assistant to generate a segmentation metric function. Test it on masks where the correct IoU and Dice can be computed by hand.
5. Ask for a comparison between Dice loss and cross-entropy. Create a numerical example that confirms or refutes the explanation.

9.35 8.34 Summary

Segmentation transforms an image into a structured label field. The clean mathematical formulation is not one method but a family of related models:

\[ \widehat y = \arg\min_y \left[\text{data fidelity}(I,y)+\text{regularity}(y)\right]. \]

Thresholding and Otsu’s method solve simple intensity-based versions. \(k\)-means and Gaussian mixtures segment in feature space. Graph cuts add spatial smoothness and solve certain binary energies exactly. Spectral methods use graph Laplacian eigenvectors to reveal grouping structure. Variational models use boundary length, total variation, and level sets. Neural models learn dense image-to-label maps from data, often using encoder-decoder architectures and losses such as cross-entropy and Dice.

For applied mathematics and statistics students, segmentation is a central example of how geometry, probability, optimization, and learning interact in modern computer vision.

9.36 References

This chapter builds on classical and modern ideas in image processing, graph-based vision, variational segmentation, and neural segmentation; see Gonzalez and Woods (2018), Szeliski (2022), Otsu (1979), Boykov et al. (2001), Kolmogorov and Zabih (2004), Shi and Malik (2000), Mumford and Shah (1989), Chan and Vese (2001), Long et al. (2015), and Ronneberger et al. (2015).

Boykov, Yuri, Olga Veksler, and Ramin Zabih. 2001. “Fast Approximate Energy Minimization via Graph Cuts.” IEEE Transactions on Pattern Analysis and Machine Intelligence 23 (11): 1222–39.

Chan, Tony F., and Luminita A. Vese. 2001. “Active Contours Without Edges.” IEEE Transactions on Image Processing 10 (2): 266–77.

Gonzalez, Rafael C., and Richard E. Woods. 2018. Digital Image Processing. 4th ed. Pearson.

Kolmogorov, Vladimir, and Ramin Zabih. 2004. “What Energy Functions Can Be Minimized via Graph Cuts?” IEEE Transactions on Pattern Analysis and Machine Intelligence 26 (2): 147–59.

Long, Jonathan, Evan Shelhamer, and Trevor Darrell. 2015. “Fully Convolutional Networks for Semantic Segmentation.” IEEE Conference on Computer Vision and Pattern Recognition, 3431–40.

Mumford, David, and Jayant Shah. 1989. “Optimal Approximations by Piecewise Smooth Functions and Associated Variational Problems.” Communications on Pure and Applied Mathematics 42 (5): 577–685.

Otsu, Nobuyuki. 1979. “A Threshold Selection Method from Gray-Level Histograms.” IEEE Transactions on Systems, Man, and Cybernetics 9 (1): 62–66.

Ronneberger, Olaf, Philipp Fischer, and Thomas Brox. 2015. “U-Net: Convolutional Networks for Biomedical Image Segmentation.” Medical Image Computing and Computer-Assisted Intervention.

Shi, Jianbo, and Jitendra Malik. 2000. “Normalized Cuts and Image Segmentation.” IEEE Transactions on Pattern Analysis and Machine Intelligence 22 (8): 888–905.

Szeliski, Richard. 2022. Computer Vision: Algorithms and Applications. 2nd ed. Springer. https://szeliski.org/Book/.