Digital Image Processing: Image Segmentation

Theoretical Foundations and Mathematical Definition of Image Segmentation

Image segmentation partitions a digital image into its constituent, non-overlapping regions or objects. Formally, let $R$ represent the entire spatial region occupied by an image. Image segmentation is a process that partitions $R$ into $n$ sub-regions, $R_1, R_2, \dots, R_n$, such that the following five fundamental conditions are satisfied:

1. Completeness:

$\bigcup_{i=1}^n R_i = R$

Every pixel in the image must belong to a region; no pixels are left unassigned.

2. Connectivity:

$R_i \text{ is a connected set for } i = 1, 2, \dots, n$

The pixels within each partitioned sub-region must form a connected component in terms of spatial adjacency.

3. Disjointness (Mutually Exclusive Sub-regions):

$R_i \cap R_j = \emptyset \quad \text{for all } i \neq j$

Regions must be strictly disjoint, meaning no pixel can simultaneously belong to more than one region.

4. Regional Homogeneity (Predicate Satisfaction):

$P(R_i) = \text{TRUE} \quad \text{for } i = 1, 2, \dots, n$

A logical predicate $P(R_i)$ evaluates the uniformity of pixels in region $R_i$ (e.g., uniformity in intensity, texture, or color) and evaluates to $\text{TRUE}$ if all pixels within $R_i$ satisfy the specified homogeneity condition.

5. Regional Distinctness:

$P(R_i \cup R_j) = \text{FALSE} \quad \text{for any adjacent regions } R_i \text{ and } R_j$

Adjacent regions must differ significantly with respect to the predicate $P$. Merging two adjacent distinct regions creates a composite region that fails the homogeneity criterion.

Segmentation algorithms generally follow one of two basic approaches based on intensity properties:

• Discontinuity-based segmentation: Partitions an image based on abrupt, local changes in intensity (such as isolated points, lines, and edges).

• Similarity-based (region-based) segmentation: Partitions an image into regions that are similar according to a set of predefined criteria (such as thresholding, region growing, and region splitting/merging).

When segmenting a smooth region of constant intensity, edge detection based on intensity discontinuities reliably recovers the closed outer boundary. However, when segmenting textured or noisy regions, edge-based operations produce many spurious internal and external edges, preventing the formation of a clean, closed boundary. In such cases, segmentation based on region properties (e.g., texture descriptors) produces a coherent, segmented partition.

Digital Derivatives and Intensity Discontinuity Profiles

Abrupt local transitions in gray-level values are modeled mathematically using discrete first-order and second-order derivatives along 1D spatial profiles.

1D Discrete Derivatives

For a discrete 1D intensity function $f(x)$, derivatives are approximated via finite differences:

• First-order derivative (Forward difference):

$\frac{\partial f}{\partial x} = f'(x) = f(x + 1) - f(x)$

• Second-order derivative (Central second difference):

$\frac{\partial^2 f}{\partial x^2} = f''(x) = f(x + 1) + f(x - 1) - 2f(x)$

Numerical Analysis of Transition Profiles

The behavior of first and second derivatives across fundamental image features—such as ramps, isolated points, flat segments, lines, and steps—is illustrated by the following discrete profile:

• Image Strip Intensity Values:

$[5,\, 5,\, 4,\, 3,\, 2,\, 1,\, 0,\, 0,\, 0,\, 6,\, 0,\, 0,\, 0,\, 0,\, 1,\, 3,\, 1,\, 0,\, 0,\, 0,\, 0,\, 7,\, 7,\, 7,\, 7]$

• First-Order Derivative Values:

$[-1,\, -1,\, -1,\, -1,\, -1,\, 0,\, 0,\, 6,\, -6,\, 0,\, 0,\, 0,\, 1,\, 2,\, -2,\, -1,\, 0,\, 0,\, 0,\, 7,\, 0,\, 0,\, 0]$

• Second-Order Derivative Values:

$[-1,\, 0,\, 0,\, 0,\, 0,\, 1,\, 0,\, 6,\, -12,\, 6,\, 0,\, 0,\, 1,\, 1,\, -4,\, 1,\, 1,\, 0,\, 0,\, 7,\, -7,\, 0,\, 0]$

Derivative Behaviors and Comparative Properties

1. Edge Thickness: First-order derivatives produce thicker, constant-amplitude responses along ramps. Second-order derivatives produce a zero response along flat slopes of a ramp and non-zero responses only at the onset and termination of the ramp.

2. Sensitivity to Fine Detail: Second-order derivatives have a significantly stronger response to fine details (isolated points, thin lines, impulse noise) than first-order derivatives.

3. Double-Edge Response: Second-order derivatives produce a characteristic double-line effect at ramp and step transitions (a positive spike followed by a negative spike, or vice versa).

4. Sign and Edge Transition Direction: The sign of the second derivative determines whether an intensity transition is from light to dark (positive then negative) or from dark to light (negative then positive).

5. Zero-Crossing: The point where the second derivative changes sign passes through zero precisely at the center of a ramp edge, providing a reliable marker for edge localization.

Detection of Isolated Points and Lines

Point Detection

An isolated point differs significantly from its local background. Detecting an isolated point requires a spatial operator that responds strongly to sharp localized discontinuities in all directions. This is achieved using the 2D Laplacian operator.

Continuous and Discrete Formulation of the 2D Laplacian

The 2D continuous Laplacian is defined as:

$\nabla^2 f(x, y) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}$

Using discrete finite differences along the $x$ and $y$ axes:

$\frac{\partial^2 f}{\partial x^2} = f(x + 1, y) + f(x - 1, y) - 2f(x, y)$

$\frac{\partial^2 f}{\partial y^2} = f(x, y + 1) + f(x, y - 1) - 2f(x, y)$

Combining these yields the standard 4-neighbor discrete Laplacian:

$\nabla^2 f(x, y) = f(x + 1, y) + f(x - 1, y) + f(x, y + 1) + f(x, y - 1) - 4f(x, y)$

Including diagonal terms yields the isotropic 8-neighbor Laplacian mask:

$\begin{bmatrix} 1 & 1 & 1 \\ 1 & -8 & 1 \\ 1 & 1 & 1 \end{bmatrix}$

Spatial Mask Filtering and Thresholding

For a $3 \times 3$ neighborhood centered at $(x, y)$ with mask weights $w_k$ and pixel intensities $z_k$:

$R(x, y) = \sum_{k=1}^9 w_k z_k$

A point at $(x, y)$ is classified as an isolated point if the absolute value of the response exceeds a specified non-negative threshold $T$:

$g(x, y) = \begin{cases} 1 & \text{if } |R(x, y)| \ge T \\ 0 & \text{otherwise} \end{cases}$

In industrial non-destructive testing, such as examining X-ray images of turbine blades, manufacturing porosity often appears as a single dark pixel. Convolving the image with the point detection mask and thresholding the output isolates the porosity.

Line Detection

Second derivatives produce a stronger response and thinner lines than first derivatives. However, the double-line effect of the second derivative must be handled by taking either the absolute value or only the positive values of the Laplacian.

Direction-Specific Line Detectors

To detect lines oriented in specific directions, four $3 \times 3$ masks are used:

• Horizontal Mask ($R_1$):

$\begin{bmatrix} -1 & -1 & -1 \\ 2 & 2 & 2 \\ -1 & -1 & -1 \end{bmatrix}$

• $+45^\circ$ Diagonal Mask ($R_2$):

$\begin{bmatrix} 2 & -1 & -1 \\ -1 & 2 & -1 \\ -1 & -1 & 2 \end{bmatrix}$

• Vertical Mask ($R_3$):

$\begin{bmatrix} -1 & 2 & -1 \\ -1 & 2 & -1 \\ -1 & 2 & -1 \end{bmatrix}$

• $-45^\circ$ Diagonal Mask ($R_4$):

$\begin{bmatrix} -1 & -1 & 2 \\ -1 & 2 & -1 \\ 2 & -1 & -1 \end{bmatrix}$

Each mask sums to zero, ensuring a zero response over regions of constant intensity.

Directional Selection Rule

Let $R_1, R_2, R_3, R_4$ be the responses of the four masks at a given point $(x, y)$. The point is associated with the direction of mask $k$ if:

$|R_k| > |R_j| \quad \text{for all } j \neq k$

To find lines of a specific orientation (e.g., $+45^\circ$), the image is filtered with that mask, negative values are set to zero, and the result is thresholded at $g \ge T$.

Edge Detection Principles and Models

An edge is a set of connected pixels that lie on the boundary between two regions characterized by distinct intensity levels.

Edge Models

Edges are classified into three geometric models:

1. Step Edge: An abrupt, instantaneous change in intensity occurring over a distance of 1 pixel. Step edges occur at clean physical boundaries and borders in range imaging.

2. Ramp Edge: An intensity transition that increases or decreases gradually over a finite distance (several pixels). Real-world optical blurring, lighting transitions, and finite sensor sampling render most edges as ramps.

3. Roof Edge: An intensity profile that rises and falls over a narrow region. Roof edges appear in range imaging when objects are closer to the sensor than the background, in thin line digitizations, and as roadways in satellite imagery.

In natural images, edge widths vary considerably. For instance, in a $1508 \times 1970$ pixel scene, practical ramp edge profiles can span 9 pixels, step-like transitions can span 2 pixels, and the base of a roof edge may span 3 pixels.

Mathematical Behavior Across a Ramp Edge

Across an ideal ramp transition:

• Intensity Profile $f(x)$: Constant low value, linear ramp slope, constant high value.

• First Derivative $f'(x)$: Constant zero on both flat sides, constant positive value across the ramp transition (producing a thick pulse proportional to edge width).

• Second Derivative $f''(x)$: Zero in flat regions, a positive impulse at the onset of the ramp, zero along the linear slope, a negative impulse at the end of the ramp, and a zero crossing precisely at the inflection point.

The Three-Stage Edge Detection Pipeline

Edge detection requires three sequential operations:

1. Image Smoothing for Noise Reduction: High-frequency noise corrupts derivative operators. Smoothing filters suppress noise prior to differentiation.

2. Detection of Edge Points: Local spatial operators (gradient or Laplacian) identify all candidate edge points.

3. Edge Localization: Thresholding and nonmaxima suppression extract true edge points while discarding false responses.

Impact of Additive Gaussian Noise

Noise degrades derivative operations. When zero-mean additive Gaussian noise is added to a ramp edge:

• At $\sigma = 0.0$, derivatives are clean and unambiguous.

• At $\sigma = 0.1$, slight ripples appear in derivative profiles.

• At $\sigma = 1.0$, noise peaks match derivative amplitudes.

• At $\sigma = 10.0$, derivative profiles are dominated by noise, obscuring true zero-crossings and first-derivative peaks.

lFirst-Order Gradient Operators and Vector Analysis

First-order edge detection relies on the gradient of a 2D continuous or discrete image function $f(x, y)$.

Vector Definition, Magnitude, and Direction

The gradient is a 2D column vector pointing in the direction of the maximum rate of change of $f$:

$\nabla f \equiv \text{grad}(f) = \begin{bmatrix} g_x \\ g_y \end{bmatrix} = \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{bmatrix}$

• Magnitude of the Gradient:

$M(x, y) = \text{mag}(\nabla f) = \sqrt{g_x^2 + g_y^2}$

• Computational Approximations:

$\text{mag}(\nabla f) \approx |g_x| + |g_y|$

$\text{mag}(\nabla f) \approx \max(|g_x|,\, |g_y|)$

• Direction of the Gradient Vector:

$\alpha(x, y) = \tan^{-1}\left[\frac{g_x}{g_y}\right]$

• Direction of the Edge:

The physical direction of an edge at $(x, y)$ is perpendicular to the gradient vector:

$\phi = \alpha - 90^\circ$

Spatial Gradient Masks

For a $3 \times 3$ neighborhood labeled as follows:

$\begin{bmatrix} z_1 & z_2 & z_3 \\ z_4 & z_5 & z_6 \\ z_7 & z_8 & z_9 \end{bmatrix}$

1D Finite Difference Masks

• Vertical component $g_x$: $\begin{bmatrix} -1 \\ 1 \end{bmatrix}$

• Horizontal component $g_y$: $\begin{bmatrix} -1 & 1 \end{bmatrix}$

Roberts Cross-Gradient Operators ($2 \times 2$)

$g_x = (z_9 - z_5) \implies \begin{bmatrix} -1 & 0 \\ 0 & 1 \end{bmatrix}$

$g_y = (z_8 - z_6) \implies \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}$

Prewitt Operators ($3 \times 3$)

• Horizontal ($g_x$):

$\begin{bmatrix} -1 & -1 & -1 \\ 0 & 0 & 0 \\ 1 & 1 & 1 \end{bmatrix}$

• Vertical ($g_y$):

$\begin{bmatrix} -1 & 0 & 1 \\ -1 & 0 & 1 \\ -1 & 0 & 1 \end{bmatrix}$

Sobel Operators ($3 \times 3$)

Sobel operators apply a weight of 2 to the central coefficient to smooth along the direction orthogonal to differentiation, reducing noise sensitivity:

• Horizontal ($g_x$):

$\begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ 1 & 2 & 1 \end{bmatrix}$

• Vertical ($g_y$):

$\begin{bmatrix} -1 & 0 & 1 \\ -2 & 0 & 2 \\ -1 & 0 & 1 \end{bmatrix}$

Diagonal Edge Operators

• Prewitt Diagonal Masks:

$\begin{bmatrix} 0 & 1 & 1 \\ -1 & 0 & 1 \\ -1 & -1 & 0 \end{bmatrix} \quad \text{and} \quad \begin{bmatrix} -1 & -1 & 0 \\ -1 & 0 & 1 \\ 0 & 1 & 1 \end{bmatrix}$

• Sobel Diagonal Masks:

$\begin{bmatrix} 0 & 1 & 2 \\ -1 & 0 & 1 \\ -2 & -1 & 0 \end{bmatrix} \quad \text{and} \quad \begin{bmatrix} -2 & -1 & 0 \\ -1 & 0 & 1 \\ 0 & 1 & 2 \end{bmatrix}$

Practical Filtering Pipeline on Complex Scenes

In complex scenes (such as an $834 \times 1114$ pixel image scaled to $[0, 1]$):

1. Raw gradient filtering via Sobel masks yields directional components $|g_x|$ and $|g_y|$, and composite gradient $|g_x| + |g_y|$.

2. Smoothing with a $5 \times 5$ averaging filter prior to computing gradients eliminates micro-texture (e.g., mortar lines in brick walls) while preserving structural edges.

3. Setting a global gradient threshold at $33\%$ of the maximum value isolates dominant structural boundaries.

Advanced Edge Detectors: Marr-Hildreth and Canny

The Marr-Hildreth Edge Detector

The Marr-Hildreth algorithm filters an image with a Gaussian to control detail scale, followed by a second-order Laplacian operator to locate edges via zero-crossings.

The Laplacian of Gaussian (LoG)

The 2D Gaussian function is:

$G(x, y) = e^{-\frac{x^2 + y^2}{2\sigma^2}}$

where $\sigma$ is the space constant (standard deviation).

Applying the Laplacian operator yields the Laplacian of Gaussian (LoG):

$\nabla^2 G(x, y) = \frac{\partial^2 G}{\partial x^2} + \frac{\partial^2 G}{\partial y^2} = \left[\frac{x^2 + y^2 - 2\sigma^2}{\sigma^4}\right] e^{-\frac{x^2 + y^2}{2\sigma^2}}$

The cross-section exhibits a central lobe flanked by negative lobes ("Mexican hat"). The distance between its zero-crossings is:

$w = 2\sqrt{2}\sigma$

A discrete $5 \times 5$ integer approximation to the negative of this shape is:

$\begin{bmatrix} 0 & 0 & -1 & 0 & 0 \\ 0 & -1 & -2 & -1 & 0 \\ -1 & -2 & 16 & -2 & -1 \\ 0 & -1 & -2 & -1 & 0 \\ 0 & 0 & -1 & 0 & 0 \end{bmatrix}$

Marr-Hildreth Algorithm Steps

1. Filter the input image $f(x, y)$ with an $n \times n$ Gaussian lowpass filter, where $n$ is chosen as the smallest odd integer greater than or equal to $6\sigma$:

$n \ge 6\sigma$

2. Compute the Laplacian of the smoothed image:

$g(x, y) = \nabla^2 [G(x, y) \star f(x, y)] = [\nabla^2 G(x, y)] \star f(x, y)$

3. Determine the zero-crossings of $g(x, y)$. A pixel is marked as a zero-crossing if the signs of opposing neighbors differ and their absolute difference exceeds a threshold.

Thresholding zero-crossings at $0$ yields closed-loop edges ("spaghetti effect"). Setting the threshold to $4\%$ of the maximum value of the LoG response suppresses weak zero-crossings and preserves thin, meaningful edges.

The Canny Edge Detector

The Canny detector satisfies three key performance criteria for step edges corrupted by white noise:

• Low Error Rate: All true edges should be detected, and no false edges should be created.

• Localization: The distance between detected edge pixels and actual physical edge centers must be minimized.

• Single Edge Response: The detector must return only a single point for each edge pixel, suppressing false multiple responses across broad profiles.

Canny Edge Detection Algorithm

Step 1: Gaussian Smoothing

Smooth the input image $f(x, y)$ using a Gaussian filter $G(x, y)$:

$f_s(x, y) = G(x, y) \star f(x, y)$

Step 2: Gradient Vector Computation

Compute the gradient magnitude $M(x, y)$ and direction $\alpha(x, y)$ using Sobel masks:

$M(x, y) = \sqrt{g_x^2 + g_y^2}$

$\alpha(x, y) = \arctan\left(\frac{g_y}{g_x}\right)$

Step 3: Nonmaxima Suppression (NMS)

Gradient magnitude arrays display wide ridges around local maxima. Nonmaxima suppression thins these ridges to 1-pixel-wide lines.

Define four basic discrete edge directions in a $3 \times 3$ neighborhood:

• $d_1$: Horizontal ($0^\circ$)

• $d_2$: $-45^\circ$

• $d_3$: Vertical ($90^\circ$)

• $d_4$: $+45^\circ$

The normal direction $\alpha$ is mapped into one of four sectors:

1. If $\alpha \in [-22.5^\circ, 22.5^\circ]$ or $\alpha \in [157.5^\circ, -157.5^\circ]$, the edge normal is horizontal, corresponding to a vertical edge ($d_3$).

2. If $\alpha \in [22.5^\circ, 67.5^\circ]$ or $\alpha \in [-157.5^\circ, -112.5^\circ]$, the edge direction is $-45^\circ$ ($d_2$).

3. If $\alpha \in [67.5^\circ, 112.5^\circ]$ or $\alpha \in [-112.5^\circ, -67.5^\circ]$, the edge normal is vertical, corresponding to a horizontal edge ($d_1$).

4. If $\alpha \in [112.5^\circ, 157.5^\circ]$ or $\alpha \in [-67.5^\circ, -22.5^\circ]$, the edge direction is $+45^\circ$ ($d_4$).

Let $p$ be the current pixel. Inspect its two neighbors along the mapped normal direction $d_k$. If $M(x, y)$ is strictly less than at least one of these two neighbors, set $g_N(x, y) = 0$; otherwise, retain $g_N(x, y) = M(x, y)$.

Step 4: Hysteresis Thresholding

Suppress false edge points using two thresholds: a high threshold $T_H$ and a low threshold $T_L$ (typically with a ratio $T_H / T_L$ between $2:1$ and $3:1$):

$g_{NH}(x, y) = g_N(x, y) \ge T_H$

$g_{NL}(x, y) = g_N(x, y) \ge T_L$

Subtract strong edges from weak edges to form an isolated weak-edge image:

$g_{NL}(x, y) = g_{NL}(x, y) - g_{NH}(x, y)$

Step 5: Connectivity Analysis and Edge Linking

1. Locate the next unvisited strong edge pixel $p$ in $g_{NH}(x, y)$.

2. Mark as valid edge pixels all weak pixels in $g_{NL}(x, y)$ that are connected to $p$ via an 8-connected path.

3. Repeat until all non-zero pixels in $g_{NH}(x, y)$ have been visited.

4. Set all unmarked pixels in $g_{NL}(x, y)$ to zero.

Edge Linking and Boundary Detection

Edge detection often yields broken boundaries due to noise, uneven illumination, and non-uniform material properties. Edge linking algorithms assemble edge pixels into continuous, closed boundaries.

Local Processing

Local processing inspects a small neighborhood $S_{xy}$ centered at every declared edge point $(x, y)$. A neighboring edge pixel at $(s, t) \in S_{xy}$ is linked to $(x, y)$ if both magnitude and directional similarity criteria are met:

• Magnitude Similarity:

$|M(s, t) - M(x, y)| \le E$

• Angle Similarity:

$|\alpha(s, t) - \alpha(x, y)| \le A$

Alternatively, a binary linked image $g(x, y)$ is constructed directly using thresholds on gradient magnitude and angle:

$g(x, y) = \begin{cases} 1 & \text{if } M(x, y) > T_M \text{ and } \alpha(x, y) = A \pm T_A \\ 0 & \text{otherwise} \end{cases}$

where $T_M$ is the magnitude threshold, $A$ is the specified angle, and $T_A$ defines an acceptable angular band about $A$.

Rows or columns are subsequently scanned to fill isolated gaps of zeros that do not exceed a length $K$.

Regional Processing: Iterative Polygonal Fit

When edge points delimit a closed or open boundary of a known region, polygonal approximations simplify the boundary representation.

The Polygonal Fit Algorithm

Let $P$ be an ordered sequence of distinct 1-valued contour pixels. Define two stacks: $\text{OPEN}$ and $\text{CLOSED}$, and specify a tolerance threshold $T$.

1. Initialization:

   • If the curve is closed: select starting vertices $A$ and $B$. Push $A$ into $\text{OPEN}$, and push $B$ into both $\text{OPEN}$ and $\text{CLOSED}$.

   • If the curve is open: push $A$ into $\text{OPEN}$ and $B$ into $\text{CLOSED}$.

2. Line Fitting: Compute the parameters of a straight line connecting the top vertex of $\text{CLOSED}$ to the top vertex of $\text{OPEN}$.

3. Distance Evaluation: Compute the perpendicular distances from this line to all points in $P$ located between these two vertices. Find the point $V_{\max}$ with the maximum distance $D_{\max}$.

4. Branching Condition:

   • If $D_{\max} > T$, push $V_{\max}$ onto $\text{OPEN}$ as a new vertex. Return to Step 2.

   • If $D_{\max} \le T$, pop the top vertex from $\text{OPEN}$ and push it onto $\text{CLOSED}$.

5. Termination: If $\text{OPEN}$ is not empty, return to Step 2. If $\text{OPEN}$ is empty, exit. The points in $\text{CLOSED}$ form the vertices of the polygonal approximation.

Execution Trace for a Closed Curve

Global Processing: The Hough Transform

When broken edge points cannot be grouped locally, global arrangements are extracted using the Hough Transform, developed by Paul Hough in 1962.

Parameterization of Lines

Representing lines in Cartesian form ($y = ax + b$) fails for vertical lines because the slope $a \rightarrow \infty$. The normal (polar) parameterization resolves this issue:

$x \cos\theta + y \sin\theta = \rho$

where:

• $\rho$ is the perpendicular distance from the coordinate origin to the line.

• $\theta$ is the angle of the normal vector relative to the horizontal axis.

Properties of Parameter Space

1. A single point $(x_i, y_i)$ in the spatial plane maps to a sinusoidal curve in the $\rho\theta$-parameter plane:

$\rho = x_i \cos\theta + y_i \sin\theta$

2. Collinear spatial points correspond to sinusoidal curves that intersect at a single parameter point $(\rho', \theta')$.

3. The parameter space is discretized into accumulator cells spanning $\theta \in [-90^\circ, 90^\circ]$ and $\rho \in [-\sqrt{M^2 + N^2},\, \sqrt{M^2 + N^2}]$ for an $M \times N$ image.

Edge Linking Algorithm Based on the Hough Transform

1. Compute a binary edge image (e.g., using Canny or Sobel operators).

2. Divide the $\rho\theta$-plane into accumulator cells $A(\rho, \theta)$ initialized to zero.

3. For each non-zero edge pixel $(x_k, y_k)$, compute $\rho$ for every quantized $\theta$ value and increment the corresponding cell: $A(\rho, \theta) = A(\rho, \theta) + 1$.

4. Find the local maxima in $A(\rho, \theta)$ that exceed a count threshold.

5. Map peak cells back to lines in image space and link collinear edge segments by bridging gaps that are smaller than a specified pixel distance.

Intensity Thresholding and Noise Degradation

Formulations of Thresholding

• Single Global Thresholding:

$g(x, y) = \begin{cases} 1 & \text{if } f(x, y) > T \quad \text{(object point)} \\ 0 & \text{if } f(x, y) \le T \quad \text{(background point)} \end{cases}$

• Multiple Thresholding:

$g(x, y) = \begin{cases} a & \text{if } f(x, y) > T_2 \\ b & \text{if } T_1 < f(x, y) \le T_2 \\ c & \text{if } f(x, y) \le T_1 \end{cases}$

Role of Noise and Illumination Non-Uniformity

In ideal noiseless images, histograms of bimodal scenes exhibit distinct spikes separated by broad valleys. However, two phenomena compromise global thresholding:

1. Additive Gaussian Noise: As noise variance increases (from $\sigma = 0$ to $\sigma = 10$ and $\sigma = 50$), the valley separating the foreground and background modes broadens, flattens, and eventually merges into a single unimodal distribution.

2. Spatial Illumination Variations: Multiplying a scene by an intensity gradient (e.g., a spatial ramp from $0.2$ to $0.6$) shifts the local reflectance values, causing the modes to spread across the intensity spectrum and obscuring the valley.

Basic Global Thresholding (Iterative Algorithm)

1. Select an initial estimate for the global threshold $T$ (typically the average intensity of the image).

2. Segment the image using $T$, producing two pixel groups:

   • $G_1$: pixels with $f(x, y) > T$

   • $G_2$: pixels with $f(x, y) \le T$

3. Compute the average intensities $m_1$ and $m_2$ of regions $G_1$ and $G_2$:

$m_1 = \frac{1}{|G_1|} \sum_{(x, y) \in G_1} f(x, y)$

$m_2 = \frac{1}{|G_2|} \sum_{(x, y) \in G_2} f(x, y)$

4. Compute an updated threshold value:

$T_{\text{new}} = \frac{1}{2}(m_1 + m_2)$

5. Repeat Steps 2 through 4 until the difference between successive thresholds satisfies:

$|T_{\text{new}} - T_{\text{old}}| < \Delta T$

where $\Delta T$ is a predefined convergence threshold.

Optimal Global Thresholding: Otsu's Method

Otsu's method (1979) selects a global threshold by maximizing the between-class variance $\sigma_B^2$, which measures the statistical separability of classes.

Statistical Formulation

Let an image of size $M \times N$ have $L$ distinct gray levels $\{0, 1, 2, \dots, L - 1\}$. Let $n_i$ denote the count of pixels with intensity $i$.

The normalized histogram represents the probability distribution of gray levels:

$p_i = \frac{n_i}{MN} \quad \text{such that} \quad \sum_{i=0}^{L-1} p_i = 1$

Partition the pixels into two classes using a threshold $k$ ($0 \le k \le L - 1$):

• $C_1$: gray levels in $[0, k]$

• $C_2$: gray levels in $[k + 1, L - 1]$

Cumulative Probabilities of Occurrence

$P_1(k) = \sum_{i=0}^k p_i$

$P_2(k) = \sum_{i=k+1}^{L-1} p_i = 1 - P_1(k)$

Class Means

$m_1(k) = \sum_{i=0}^k i P(i \mid C_1) = \frac{1}{P_1(k)} \sum_{i=0}^k i p_i$

$m_2(k) = \sum_{i=k+1}^{L-1} i P(i \mid C_2) = \frac{1}{P_2(k)} \sum_{i=k+1}^{L-1} i p_i$

Global Mean and Cumulative Mean

The cumulative mean up to intensity level $k$ is:

$m(k) = \sum_{i=0}^k i p_i = P_1(k) m_1(k)$

The global mean intensity across the entire image is:

$m_G = \sum_{i=0}^{L-1} i p_i = P_1 m_1 + P_2 m_2$

Between-Class Variance Derivation

The between-class variance $\sigma_B^2$ is defined as:

$\sigma_B^2(k) = P_1(m_1 - m_G)^2 + P_2(m_2 - m_G)^2$

Substituting $m_G = P_1 m_1 + P_2 m_2$ and $P_2 = 1 - P_1$ yields:

$\sigma_B^2(k) = P_1 P_2 (m_1 - m_2)^2 = \frac{[m_G P_1(k) - m(k)]^2}{P_1(k)[1 - P_1(k)]}$

Threshold Optimization and Separability Measure

The optimal threshold $k^*$ is the index that maximizes $\sigma_B^2(k)$:

$\sigma_B^2(k^*) = \max_{0 \le k < L-1} \sigma_B^2(k)$

If the maximum is not unique, the optimal threshold is calculated as the average of the maximizing indices.

The global variance of the image is:

$\sigma_G^2 = \sum_{i=0}^{L-1} (i - m_G)^2 p_i$

The separability measure $\eta$ evaluates the quality of the thresholding ($0 \le \eta \le 1$):

$\eta = \frac{\sigma_B^2(k^*)}{\sigma_G^2}$

Algorithmic Summary of Otsu's Method

1. Compute the normalized histogram $p_i$ for $i = 0, 1, \dots, L - 1$.

2. Compute the cumulative sums $P_1(k)$ for $k = 0, 1, \dots, L - 1$.

3. Compute the cumulative means $m(k)$ for $k = 0, 1, \dots, L - 1$.

4. Compute the global intensity mean $m_G$.

5. Compute the between-class variance $\sigma_B^2(k)$ for all $k$.

6. Find $k^*$ that maximizes $\sigma_B^2(k)$.

7. Segment the image into background and foreground using $k^*$:

$g(x, y) = \begin{cases} 1 & \text{if } f(x, y) > k^* \\ 0 & \text{if } f(x, y) \le k^* \end{cases}$

Enhancements for Challenging Conditions

Image Smoothing to Restore Bimodality

When heavy additive noise blurs histogram peaks into a unimodal distribution, Otsu's method on the raw image often fails. Pre-filtering with a $5 \times 5$ spatial averaging mask attenuates noise variance, re-establishing a bimodal profile and enabling successful thresholding.

Edge-Assisted Global Thresholding

When an object is small relative to the image area, background pixels dominate the histogram, suppressing the object mode into an imperceptible tail and causing Otsu's method to fail.

To correct for skewed class proportions:

1. Compute an edge image using either the gradient magnitude or the absolute Laplacian of $f(x, y)$.

2. Threshold the edge image at a high percentile (e.g., the $99.7^{\text{th}}$ percentile) to produce a sparse binary edge mask.

3. Compute an intensity histogram using only the pixels in $f(x, y)$ that coincide with the non-zero mask pixels.

4. Apply Otsu's algorithm to this edge-conditioned histogram to determine the optimal threshold (e.g., $T = 134$), which will lie midway between the object and background peaks.

5. Segment the entire original image globally using this threshold.

Multiple Thresholds for Multimodal Distributions

For $K$ distinct classes $C_1, C_2, \dots, C_K$, the between-class variance generalized to $K - 1$ thresholds $\{k_1, k_2, \dots, k_{K-1}\}$ is:

$\sigma_B^2 = \sum_{k=1}^K P_k (m_k - m_G)^2$

where:

$P_k = \sum_{i \in C_k} p_i \quad \text{and} \quad m_k = \frac{1}{P_k} \sum_{i \in C_k} i p_i$

The optimal thresholds satisfy:

$\sigma_B^2(k_1^*, k_2^*, \dots, k_{K-1}^*) = \max_{0 \le k_1 < k_2 < \dots < L-1} \sigma_B^2(k_1, k_2, \dots, k_{K-1})$

Spatially Adaptive and Variable Thresholding

When background illumination varies across an image, global thresholding fails. Spatially adaptive methods vary the threshold $T(x, y)$ as a function of local coordinates.

Image Partitioning

The image is subdivided into non-overlapping rectangular blocks chosen small enough that illumination across each block is approximately uniform. Otsu's method is then applied to each subimage individually.

Variable Thresholding Based on Local Image Properties

Let $\sigma_{xy}$ and $m_{xy}$ denote the local standard deviation and local mean within a neighborhood $S_{xy}$ centered at $(x, y)$. The local threshold is defined as:

$T_{xy} = a\sigma_{xy} + bm_{xy}$

If the background intensity is nearly constant across the scene, the global mean $m_G$ can replace the local mean:

$T_{xy} = a\sigma_{xy} + bm_G$

The image is thresholded according to:

$g(x, y) = \begin{cases} 1 & \text{if } f(x, y) > T_{xy} \\ 0 & \text{if } f(x, y) \le T_{xy} \end{cases}$

Dual Predicate Formulation

Alternatively, local parameters can be evaluated using a logical predicate $Q$:

$g(x, y) = \begin{cases} 1 & \text{if } Q(\text{local parameters}) \text{ is TRUE} \\ 0 & \text{otherwise} \end{cases}$

For example:

$Q(\sigma_{xy}, m_{xy}) = \begin{cases} \text{TRUE} & \text{if } f(x, y) > a\sigma_{xy} \text{ AND } f(x, y) > bm_{xy} \\ \text{FALSE} & \text{otherwise} \end{cases}$

Variable Thresholding Using Moving Averages

For text document processing and scenes containing small features affected by illumination gradients, local thresholding based on moving averages provides efficient, real-time segmentation. To prevent spatial bias, images are scanned row-by-row in a continuous zigzag pattern.

Recursive Moving Average Computation

Let $z_{k+1}$ denote the intensity of the pixel encountered at step $k + 1$ along the scan path. The moving average over a window of length $n$ is updated recursively:

$m(k + 1) = \frac{1}{n} \sum_{i=k+2-n}^{k+1} z_i = m(k) + \frac{1}{n}(z_{k+1} - z_k)$

Border pixels are padded with $n - 1$ zeros, and the initial value is $m(1) = z_1 / n$.

Segmentation Criterion

The local threshold at $(x, y)$ is scaled using a constant multiplier $b$:

$T_{xy} = b \cdot m_{xy}$

$g(x, y) = \begin{cases} 1 & \text{if } f(x, y) > T_{xy} \\ 0 & \text{if } f(x, y) \le T_{xy} \end{cases}$

In document restoration corrupted by spot or sinusoidal shading, parameters such as $n = 20$ and $b = 0.5$ consistently isolate dark characters from non-uniform backgrounds.

Region-Based Segmentation: Growing, Splitting, and Merging

Region Growing by Pixel Aggregation

Region growing groups pixels into larger spatial regions based on predefined criteria. Starting from a set of "seed" points, regions grow by iteratively appending adjacent pixels that satisfy a similarity predicate (e.g., intensity, texture, color, or shape).

Algorithm Based on 8-Connectivity

Let $f(x, y)$ be the input image array, $S(x, y)$ be a binary seed array ($1$ for seeds, $0$ otherwise), and $Q$ be the similarity predicate:

1. Find all connected components in $S(x, y)$ and erode each component down to a single pixel; label each resulting seed pixel as $1$, and all others as $0$.

2. Create a binary image $f_Q(x, y)$ evaluating the predicate $Q$ for every pixel:

$f_Q(x, y) = \begin{cases} 1 & \text{if } Q(x, y) \text{ is satisfied} \\ 0 & \text{otherwise} \end{cases}$

3. Form the segmented image $g$ by appending to each seed point in $S$ all 1-valued points in $f_Q$ that are 8-connected to that seed point.

4. Label each disjoint connected component in $g$ with a distinct region label.

A common similarity predicate uses an absolute intensity difference threshold $T$:

$Q = \begin{cases} \text{TRUE} & \text{if } |f(x, y) - f(\text{seed})| \le T \\ \text{FALSE} & \text{otherwise} \end{cases}$

Worked Numerical Example: 4-Connectivity vs. 8-Connectivity

Consider the following $7 \times 7$ digital image array, where the initial seed point is the center pixel at coordinates $(4, 4)$ (value $60$), and the similarity predicate is $|f(x, y) - f(4, 4)| \le 5$:

Row 1: [10, 10, 10, 10, 10, 10, 10]
Row 2: [10, 10, 10, 69, 70, 10, 10]
Row 3: [59, 10, 60, 64, 59, 56, 60]
Row 4: [10, 59, 10, 60, 70, 10, 62]
Row 5: [10, 60, 59, 65, 67, 10, 65]
Row 6: [10, 10, 10, 10, 10, 10, 10]
Row 7: [10, 10, 10, 10, 10, 10, 10]

The acceptable intensity range for region inclusion is $[60 - 5,\, 60 + 5] = [55,\, 65]$.

• Under 4-Connectivity:

  • From $(4, 4)$ (value $60$):

  • North neighbor $(3, 4)$ has value $64 \implies$ accepted.

  • South neighbor $(5, 4)$ has value $65 \implies$ accepted.

  • East neighbor $(4, 5)$ has value $70 \implies$ rejected.

  • West neighbor $(4, 3)$ has value $10 \implies$ rejected.

  • From $(3, 4)$ ($64$):

  • North neighbor $(2, 4)$ has value $69 \implies$ rejected.

  • East neighbor $(3, 5)$ has value $59 \implies$ accepted.

  • West neighbor $(3, 3)$ has value $60 \implies$ accepted.

  • From $(3, 5)$ ($59$):

  • East neighbor $(3, 6)$ has value $56 \implies$ accepted.

  • From $(3, 6)$ ($56$):

  • East neighbor $(3, 7)$ has value $60 \implies$ accepted.

  • From $(5, 4)$ ($65$):

  • West neighbor $(5, 3)$ has value $59 \implies$ accepted.

  • East neighbor $(5, 5)$ has value $67 \implies$ rejected.

  • From $(5, 3)$ ($59$):

  • West neighbor $(5, 2)$ has value $60 \implies$ accepted.

  • Pixels at $(3, 1)$ ($59$), $(4, 2)$ ($59$), and $(4, 7)$ ($62$) are isolated from 4-connected growth by surrounding values of $10$.

• Under 8-Connectivity:

  • The pixel at $(4, 2)$ ($59$) is diagonally 8-connected to $(3, 3)$ ($60$) and $(5, 3)$ ($59$), so it is included in the grown region.

  • Pixel $(3, 1)$ ($59$) connects to $(4, 2)$ ($59$).

  • Pixel $(4, 7)$ ($62$) connects diagonally to $(3, 6)$ ($56$) and $(5, 7)$ ($65$).

  • Thus, 8-connectivity bridges diagonal paths across background barriers that halt 4-connected aggregation.