Thresholding

Thresholding is the simplest way to segment an image: pixels above a threshold go to one class and the rest to the other, giving a binary mask. It is the classical image segmentation method when objects and background differ in brightness, as with printed text, cells under a microscope, or parts on a backlit conveyor. The difficult part is choosing the threshold. Global methods, such as Otsu’s, pick one value from the image histogram; local methods compute a separate threshold for every pixel from its neighborhood, which copes with uneven illumination.

Problem

Given a grayscale image whose objects are consistently brighter (or darker) than their surroundings, assign each pixel to object or background. A fixed threshold rarely transfers between images, because exposure, illumination, staining and gain shift the levels, so the threshold is chosen from each image’s own data.

Inputs and Outputs

Input: a grayscale image $I$ with $L$ intensity levels (256 for 8-bit images), or one derived channel such as a gradient magnitude.

Output: a binary mask $B(\mathbf{x}) = [\,I(\mathbf{x}) > T\,]$ for a global threshold, or $[\,I(\mathbf{x}) > T(\mathbf{x})\,]$ for a local threshold surface. Polarity is a choice: dark text is foreground below the threshold. With $M - 1$ thresholds the output is a label map with $M$ classes.

Intuition

If objects and background each have a typical brightness, the histogram has two peaks, and any threshold in the valley between them works. Noise widens the peaks until the valley is shallow or gone, so Otsu instead scored each candidate threshold by how well it splits the pixels into two tight, well-separated groups, using the whole histogram.

No single number suits a page photographed under a desk lamp, where the paper on the dark side may be darker than the ink on the bright side. Local methods compare each pixel only with its neighborhood, where the illumination is roughly constant.

Algorithm

Global thresholding with Otsu’s criterion:

1. Compute the normalized histogram $p_i$, $i = 0, \dots, L-1$.
2. Compute the cumulative sums $\omega(t) = \sum_{i \le t} p_i$ and $\mu(t) = \sum_{i \le t} i\,p_i$, and the global mean $\mu_T = \mu(L-1)$.
3. For every $t$ with $0 < \omega(t) < 1$, evaluate the between-class variance $\sigma_B^2(t)$ (below).
4. Take $t^* = \arg\max_t \sigma_B^2(t)$ and label pixels with $I > t^*$ as class 1.

A local method replaces step 4 with a per-pixel rule: compute the mean $m(\mathbf{x})$ and standard deviation $s(\mathbf{x})$ in a $w \times w$ window around each pixel, form $T(\mathbf{x})$ from them, and compare.

Mathematical Formulation

Otsu’s method

A threshold $t$ splits the levels into $C_0 = \{0, \dots, t\}$ and $C_1 = \{t+1, \dots, L-1\}$, with probabilities $\omega_0 = \omega(t)$ and $\omega_1 = 1 - \omega(t)$ and means $\mu_0 = \mu(t)/\omega_0$ and $\mu_1 = (\mu_T - \mu(t))/\omega_1$. Otsu (1979) took three measures of class separability from discriminant analysis, built from the within-class variance $\sigma_W^2 = \omega_0\sigma_0^2 + \omega_1\sigma_1^2$, the between-class variance

$$\sigma_B^2(t) = \omega_0\,\omega_1\,(\mu_1 - \mu_0)^2 = \frac{\big(\mu_T\,\omega(t) - \mu(t)\big)^2}{\omega(t)\,\big(1 - \omega(t)\big)},$$

and the total variance $\sigma_T^2$. Because $\sigma_W^2 + \sigma_B^2 = \sigma_T^2$ for every $t$ and $\sigma_T^2$ does not depend on $t$, maximizing $\sigma_B^2$ is the same as minimizing $\sigma_W^2$. Otsu adopted the ratio $\eta = \sigma_B^2/\sigma_T^2$, the simplest of his three measures, because $\sigma_B^2$ needs only the zeroth- and first-order cumulative moments $\omega(t)$ and $\mu(t)$, while $\sigma_W^2$ needs second-order statistics.

The maximum $\eta^* = \sigma_B^2(t^*)/\sigma_T^2$ lies in $[0, 1]$ and measures how separable the two classes are. Otsu noted that it is invariant to affine changes of the gray scale, equals 0 only for a constant image, and equals 1 only for an image with two values. He also stated, from his experiments but without proof, that $\sigma_B^2(t)$ is smooth and unimodal; in our check, three equal Gaussian modes give two local maxima.

Multilevel Otsu. With $M$ classes, $\sigma_B^2$ becomes a function of $M - 1$ thresholds, maximized jointly. Otsu cautioned that the thresholds become less credible as the number of classes grows, and considered two and three classes enough for most applications; his cell example separated background, cytoplasm and nucleus.

Other global criteria

• Iterative selection (isodata). Ridler and Calvard (1978) proposed an iteration that, as Sezgin and Sankur (2004) summarize it, sets each new threshold to the average of the two class means, $t_{n+1} = \tfrac{1}{2}\big(\mu_0(t_n) + \mu_1(t_n)\big)$, until it stops changing. In our reading, this is two-cluster k-means on intensities, whose objective is Otsu’s within-class variance: Otsu’s exhaustive search finds its global minimum, the iteration a fixed point that depends on the start.
• Minimum-error thresholding. Kittler and Illingworth (1986) model the histogram as a mixture of two Gaussians with unequal variances and choose the threshold that best fits it. In Sezgin and Sankur’s notation, with $P = \omega(t)$ and $\sigma_f, \sigma_b$ the standard deviations of the classes below and above $t$, the threshold minimizes $P \log \sigma_f + (1-P)\log\sigma_b - P\log P - (1-P)\log(1-P)$.
• Entropy-based thresholding. Kapur, Sahoo and Wong (1985) treat object and background as two sources and choose the threshold that maximizes the sum of the entropies of the two class distributions $p_i/\omega_0$ and $p_i/\omega_1$.

OpenCV’s threshold function also offers a triangle method next to Otsu’s.

Local thresholds

Niblack’s method (1986), as Sezgin and Sankur summarize it, sets $T(\mathbf{x}) = m(\mathbf{x}) + k\,s(\mathbf{x})$; Trier and Jain, cited there, found a $15 \times 15$ window and $k = -0.2$ satisfactory. Sauvola and Pietikäinen (2000) modified the rule; as given by Shafait et al. (2008), it is

$$T(\mathbf{x}) = m(\mathbf{x}) \left[ 1 + k \left( \frac{s(\mathbf{x})}{R} - 1 \right) \right],$$

where $R$ is the dynamic range of the standard deviation (128 for 8-bit images) and $k$ is positive. Shafait et al. explain the effect: in high-contrast windows $s \approx R$ and $T \approx m$, as with Niblack, while in low-contrast background $T$ drops below the mean, so dark background areas without text are not marked as foreground. OpenCV’s adaptiveThreshold uses the window mean, or a Gaussian-weighted mean, minus a constant $C$, normally positive.

Parameters

• Bins. Floating-point images must be binned, which quantizes the threshold to a bin edge.
• Window size $w$ (local methods). Large enough that most windows contain both foreground and background, and small enough that illumination is nearly constant inside each.
• $k$ and $R$ (Sauvola). Sezgin and Sankur used $k = 0.5$ and $R = 128$; Shafait et al. report that Sauvola and Pietikäinen also used $k = 0.5$, give a range of 0.2 to 0.5, and find results not very sensitive to $k$. scikit-image’s threshold_sauvola defaults to $k = 0.2$.

Complexity

A global threshold costs one pass to build the histogram, $O(N)$ for $N$ pixels, plus an $O(L)$ scan of the cumulative sums for Otsu’s method. Exhaustive multilevel Otsu grows as $O(L^{M-1})$. A direct local method costs $O(Nw^2)$; with integral images of $I$ and $I^2$, any window sum takes four lookups and the cost is $O(N)$, independent of $w$. Shafait et al. note that the squared integral image can overflow 32-bit integers.

Implementation

The NumPy code builds two 8-bit images with known masks: bright disks and a rectangle on an evenly lit background (A), and a page of dark strokes under illumination falling from 1.0 on the right to 0.25 on the left (B). It computes Otsu’s threshold from cumulative sums, checks it by brute-force minimization of the within-class variance, runs the Ridler–Calvard iteration, and applies Sauvola’s and Niblack’s rules to image B with integral-image window statistics. Scores are IoU of the foreground and the fraction of misclassified pixels.

import numpy as np

rng = np.random.default_rng(0)
yy, xx = np.mgrid[0:128, 0:192]

# Image A: bright blobs on a dark, evenly lit background (a bimodal histogram).
gt_a = ((yy - 40)**2 + (xx - 50)**2 < 22**2) | ((yy - 85)**2 + (xx - 135)**2 < 30**2) \
       | ((yy > 95) & (yy < 115) & (xx > 20) & (xx < 80))
img_a = np.where(gt_a, 170, 80) + rng.normal(0, 15, gt_a.shape)
img_a = np.clip(np.round(img_a), 0, 255).astype(np.uint8)

# Image B: a "page" of dark strokes (short bars and dotted underlines) whose
# illumination falls from 1.0 on the right to 0.25 on the left.
gt_b = np.zeros(gt_a.shape, bool)
for r in range(12, 100, 16):                       # six lines of "text"
    for c in range(8, 184, 9):                     # short vertical strokes
        h = rng.integers(6, 11)
        gt_b[r:r + h, c:c + 3] = True
    gt_b[r + 9:r + 11, 8:184:2] = True             # dotted underline per line
light = 0.25 + 0.75 * xx / xx.max()
img_b = 255 * light * np.where(gt_b, 0.25, 0.90) + rng.normal(0, 8, gt_b.shape)
img_b = np.clip(np.round(img_b), 0, 255).astype(np.uint8)

def otsu(img):
    """Otsu threshold t (class 0 = levels <= t) from cumulative sums of the histogram."""
    p = np.bincount(img.ravel(), minlength=256) / img.size
    levels = np.arange(256)
    w0 = np.cumsum(p)                              # omega_0(t)
    mu = np.cumsum(levels * p)                     # first cumulative moment mu(t)
    mu_T = mu[-1]
    with np.errstate(divide="ignore", invalid="ignore"):
        sb2 = (mu_T * w0 - mu)**2 / (w0 * (1 - w0))  # between-class variance
    sb2[~np.isfinite(sb2)] = 0
    return int(np.argmax(sb2)), sb2

def brute_force(img):
    """Threshold minimizing the within-class variance, computed from the pixels."""
    x = img.ravel().astype(float)
    best, best_t = np.inf, None
    for t in range(int(x.min()), int(x.max())):
        c0, c1 = x[x <= t], x[x > t]
        sw2 = (c0.size * c0.var() + c1.size * c1.var()) / x.size
        if sw2 < best:
            best, best_t = sw2, t
    return best_t, best

def isodata(img):
    """Ridler-Calvard iteration: move t to the midpoint of the two class means."""
    x, t = img.ravel().astype(float), img.mean()
    while True:
        t_new = (x[x <= t].mean() + x[x > t].mean()) / 2
        if abs(t_new - t) < 0.5:
            return t_new
        t = t_new

def box_mean(a, w):
    """Mean over a w x w window (w odd) from an integral image, edge-padded."""
    r = w // 2
    s = np.pad(np.pad(a, r, mode="edge").cumsum(0).cumsum(1), ((1, 0), (1, 0)))
    H, W = a.shape
    return (s[w:w + H, w:w + W] - s[:H, w:w + W] - s[w:w + H, :W] + s[:H, :W]) / w**2

def local_stats(img, w):
    f = img.astype(float)
    m = box_mean(f, w)
    s = np.sqrt(np.maximum(box_mean(f * f, w) - m * m, 0))
    return m, s

def report(name, pred, gt):
    iou = (pred & gt).sum() / (pred | gt).sum()
    err = (pred != gt).mean()
    print(f"  {name:<24} IoU {iou:.3f}  pixel error {err:.2%}")

for name, img, gt, dark_fg in [("A (even light)", img_a, gt_a, False),
                               ("B (light gradient)", img_b, gt_b, True)]:
    t, sb2 = otsu(img)
    tb, sw2 = brute_force(img)
    sT2 = img.astype(float).var()
    print(f"{name}: Otsu t = {t}, brute force t = {tb}; "
          f"sigma_W^2 + sigma_B^2 = {sw2 + sb2[t]:.2f}, sigma_T^2 = {sT2:.2f}, "
          f"eta = {sb2[t] / sT2:.3f}")
    print("  sigma_B^2 at t = 50, 100, 125, 150, 200:",
          np.round(sb2[[50, 100, 125, 150, 200]], 1))
    print(f"  isodata (Ridler-Calvard) t = {isodata(img):.1f}")
    report("Otsu", (img <= t) if dark_fg else (img > t), gt)

paper_as_ink = (img_b <= otsu(img_b)[0]) & ~gt_b       # Otsu errors on image B
print("  paper marked as ink, left/middle/right thirds:",
      [f"{paper_as_ink[:, i:i + 64].sum() / (~gt_b[:, i:i + 64]).sum():.1%}" for i in (0, 64, 128)])
m, s = local_stats(img_b, 31)
for k in (0.2, 0.5):
    report(f"Sauvola w=31 k={k}", img_b <= m * (1 + k * (s / 128 - 1)), gt_b)
report("Niblack w=31 k=-0.2", img_b <= m - 0.2 * s, gt_b)
blank = (yy >= 112) & (xx < 64)                    # empty margin, dark third
print(f"  in empty paper (bottom-left), marked as ink: Sauvola k=0.5 "
      f"{(img_b <= m * (1 + 0.5 * (s / 128 - 1)))[blank].mean():.1%}, "
      f"Niblack {(img_b <= m - 0.2 * s)[blank].mean():.1%}")

Output:

A (even light): Otsu t = 124, brute force t = 124; sigma_W^2 + sigma_B^2 = 1619.98, sigma_T^2 = 1619.98, eta = 0.863
  sigma_B^2 at t = 50, 100, 125, 150, 200: [  59.9 1250.6 1397.5 1330.8   46.9]
  isodata (Ridler-Calvard) t = 124.9
  Otsu                     IoU 0.993  pixel error 0.15%
B (light gradient): Otsu t = 127, brute force t = 127; sigma_W^2 + sigma_B^2 = 3687.69, sigma_T^2 = 3687.69, eta = 0.745
  sigma_B^2 at t = 50, 100, 125, 150, 200: [1216.5 2572.5 2745.4 2610.3 1390. ]
  isodata (Ridler-Calvard) t = 127.3
  Otsu                     IoU 0.316  pixel error 34.35%
  paper marked as ink, left/middle/right thirds: ['99.7%', '22.7%', '0.0%']
  Sauvola w=31 k=0.2       IoU 0.970  pixel error 0.49%
  Sauvola w=31 k=0.5       IoU 0.998  pixel error 0.03%
  Niblack w=31 k=-0.2      IoU 0.736  pixel error 5.69%
  in empty paper (bottom-left), marked as ink: Sauvola k=0.5 0.0%, Niblack 33.9%

On both images the cumulative-sum and brute-force thresholds agree, and $\sigma_W^2 + \sigma_B^2 = \sigma_T^2$. On image A (classes at 80 and 170) the threshold of 124 lies at the midpoint, under a single broad maximum of $\sigma_B^2$, and the mask has IoU 0.993; the Ridler–Calvard iteration agrees.

On image B Otsu’s method still finds a clear maximum ($\eta^* = 0.745$), but it splits the bright half of the page from the dim half, not ink from paper: almost all the paper in the left third, and 22.7% in the middle third, is labeled ink, and IoU drops to 0.316. A high $\eta^*$ says the histogram splits cleanly, not that the split is the wanted one. Sauvola’s rule with $k = 0.5$ recovers the strokes (IoU 0.998). With $k = 0.2$ the threshold sits closer to the local mean, and almost all of its 0.49% errors are paper pixels in the dimmest third, where near the left edge 20% of the paper level is under two noise standard deviations. Niblack’s rule marks a third of the empty margin as ink, because in a window with no text its threshold sits just below the local mean and noise crosses it. Over noise seeds 0–9 every score stayed within about 0.02 IoU of these; window sizes from 15 to 61 gave similar Sauvola results.

Properties and Behavior

• No spatial information. Shuffling the pixels leaves a global threshold unchanged.
• Benchmarks. In Sezgin and Sankur’s comparison of 40 global and local methods on 40 nondestructive-testing images and 40 degraded documents, Kittler and Illingworth’s method ranked first on both sets, and Sauvola’s ranked second on the documents; every method did poorly on at least one or two images, and the document results apply only to text degraded with noise and blur.

Limitations

• Uneven illumination. A global threshold fails when the background level varies across the image, as on image B. The usual remedies are local methods or flattening the background first, by dividing by a heavily smoothed copy or with a top-hat transform from mathematical morphology.
• Unimodal or skewed histograms. Small objects contribute a small mode or none. Sezgin and Sankur note that Otsu’s method gives satisfactory results when the two classes have similar pixel counts, which a small foreground violates.
• Overlapping classes. When noise or texture makes the class distributions overlap, every threshold misclassifies pixels; smoothing first with a Gaussian filter or a median filter narrows the peaks (see image noise).

Variants

• Hysteresis thresholding uses a high threshold to start regions and a low one to extend them through connected pixels, as in the Canny edge detector.
• Color thresholding applies thresholds per channel (a box in color space) or to one derived channel, such as hue or the distance to a reference color, in a suitable color space.
• Post-processing. Masks are usually cleaned with morphological opening and closing, then split into connected components, dropping those below a minimum size; touching objects stay merged unless split, for example by the watershed transform.