Histogram Equalization

Histogram equalization is a contrast enhancement method that remaps the intensity levels of an image so that, after the mapping, the levels are used approximately equally often. The mapping is the image’s own cumulative histogram, scaled to the output range: crowded parts of the histogram are spread apart and sparse parts are compressed. It needs no parameters and runs in a single pass. Its adaptive variants, especially contrast-limited adaptive histogram equalization (CLAHE), equalize local regions separately and are widely used in medical imaging and as preprocessing.

Problem

A low-contrast image uses only a narrow band of the available intensity levels: a hazy scene might occupy levels 90–150 of an 8-bit range, so differences between objects are small and hard to see. A linear contrast stretch scales every part of the histogram by the same factor, so crowded peaks stay crowded. Histogram equalization instead chooses a nonlinear mapping from the data, giving more output levels to intensity ranges that contain more pixels.

Inputs and Outputs

• Input: a grayscale image with $L$ levels (256 for 8-bit data), or the luminance channel of a color image.
• Output: an image of the same size and bit depth, together with the lookup table $T$ that maps each input level to an output level.

Intuition

Line up all pixels from darkest to brightest and assign output values by rank: the darkest 1% of pixels get the lowest 1% of the output range, the next 1% the next slice, and so on. The fraction of pixels at or below a level is exactly the cumulative histogram, so mapping each level through the cumulative histogram is rank-based relabeling.

Algorithm

1. Compute the histogram $h(k)$ of the image, the number of pixels at each level $k = 0, \dots, L-1$.
2. Accumulate it into the cumulative distribution $\mathrm{CDF}(k)$, the fraction of pixels with level at most $k$.
3. Build the lookup table $T(k) = \operatorname{round}\big((L-1)\,\mathrm{CDF}(k)\big)$.
4. Replace every pixel value $I(\mathbf{x})$ by $T(I(\mathbf{x}))$.

A common variant subtracts $\mathrm{CDF}_{\min}$, the CDF at the lowest occupied level, and divides by $1 - \mathrm{CDF}_{\min}$, so that the darkest pixels map exactly to 0.

Mathematical Formulation

With $N$ pixels, the mapping is

$$T(k) = \operatorname{round}\!\left( (L-1)\, \mathrm{CDF}(k) \right), \qquad \mathrm{CDF}(k) = \frac{1}{N} \sum_{j=0}^{k} h(j).$$

Why it flattens the histogram. Treat intensity as a continuous random variable $r$ with density $p_r$ and cumulative distribution $F_r$, and let $s = (L-1)\,F_r(r)$. Because $F_r$ is continuous and non-decreasing, $P(s \le (L-1)u) = P(F_r(r) \le u) = u$ for $u \in [0, 1]$, so $s$ is uniformly distributed on $[0, L-1]$. This is the probability integral transform: any continuous variable passed through its own CDF becomes uniform.

Why it is only approximate for images. Digital images are discrete. All pixels that share an input level receive the same output level, so a level that holds 10% of the pixels still produces a single spike holding 10% of the output; the method can move spikes but never split them. Rounding also sends several sparsely populated input levels to the same output level, and leaves gaps between occupied output levels. The output histogram is therefore flat only on average, and the number of distinct levels can never increase.

Histogram matching. Equalization is a special case of histogram matching, also called histogram specification, which transforms an image so that its histogram approximates an arbitrary target. If $G$ is the target CDF, each level is mapped through the source CDF and then through the inverse of the target CDF, $k \mapsto G^{-1}\big(\mathrm{CDF}(k)\big)$; in the discrete case $G^{-1}(q)$ is the smallest level $j$ with $G(j) \ge q$. A uniform target gives equalization.

Parameters

Global equalization has no parameters besides $L$. A practical choice is which pixels enter the histogram: excluding a large uniform background keeps it from dominating the mapping. Adaptive variants add a region size and a clip limit (see Variants).

Complexity

Building the histogram takes $O(N)$ time, the cumulative sum and lookup table $O(L)$, and applying the table $O(N)$, so the total is $O(N + L)$ time and $O(L)$ extra memory.

Implementation

The example equalizes a synthetic low-contrast image with two regions, then measures flatness as the largest gap between the image’s cumulative histogram and the uniform one (0 would be perfectly flat).

import numpy as np

rng = np.random.default_rng(0)

# Synthetic low-contrast 8-bit image: a dark left half and a brighter right
# half, each with a gentle ramp and noise, all within levels 75-135.
h, w = 128, 256
x = np.linspace(0, 1, w)
base = np.where(x < 0.5, 75 + 20 * x, 115 + 20 * x)       # two regions
img = np.tile(base, (h, 1)) + rng.normal(0, 3, (h, w))
img = np.clip(np.round(img), 0, 255).astype(np.uint8)

def equalize(im, L=256):
    hist = np.bincount(im.ravel(), minlength=L)
    cdf = np.cumsum(hist) / im.size                       # CDF(k)
    T = np.round((L - 1) * cdf).astype(np.uint8)          # T(k) = round((L-1) CDF(k))
    return T[im], T

def report(name, im, L=256):
    hist = np.bincount(im.ravel(), minlength=L)
    occupied = np.count_nonzero(hist)
    # Flatness: how far the CDF is from the uniform CDF over 0..L-1
    cdf = np.cumsum(hist) / im.size
    uniform = np.arange(1, L + 1) / L
    ks = np.abs(cdf - uniform).max()
    print(f"{name:<10} range {im.min():>3}-{im.max():>3}  occupied levels {occupied:>3}"
          f"  std {im.std():5.1f}  max |CDF - uniform| {ks:.3f}")

eq, T = equalize(img)
report("input", img)
report("equalized", eq)

# T is monotonic, so pixel order is preserved; merged input levels share an output.
print("T monotonic:", bool(np.all(np.diff(T.astype(int)) >= 0)))
used_in = np.unique(img)
print("distinct input levels:", used_in.size, "-> distinct outputs:", np.unique(T[used_in]).size)

Output:

input      range  64-145  occupied levels  64  std  25.3  max |CDF - uniform| 0.448
equalized  range   0-255  occupied levels  44  std  73.5  max |CDF - uniform| 0.046
T monotonic: True
distinct input levels: 64 -> distinct outputs: 44

The output spans the full range, its standard deviation roughly triples, and its cumulative histogram is close to uniform (largest deviation 0.046 against 0.448 before). Yet the number of occupied levels drops from 64 to 44: the sparsely populated levels in the noise tails were merged by rounding.

Properties and Behavior

• Monotonic. $T$ is non-decreasing, so the order of intensities is preserved; levels may merge but never swap.
• Global. One lookup table is applied to every pixel, regardless of where it is.
• Data-dependent. Adding or cropping a large background region changes the mapping everywhere.
• Can over-enhance. The slope of $T$ at level $k$ is proportional to $h(k)$, so a dominant peak, such as a uniform sky or a black background, is stretched over a large part of the output range while the subject is compressed.
• Can amplify noise. For the same reason, a nearly uniform region whose values cluster in a tall, narrow peak has its small noise fluctuations stretched across many output levels, making noise in flat regions more visible.
• Idempotent (plain form). With the mapping above, equalizing an already equalized image returns it unchanged. If $k^\ast$ is the highest input level mapped to an output level $m$, the fraction of output pixels at or below $m$ is $\mathrm{CDF}(k^\ast)$, so the second pass sends $m$ to $\operatorname{round}\big((L-1)\,\mathrm{CDF}(k^\ast)\big) = T(k^\ast) = m$. The $\mathrm{CDF}_{\min}$ variant is not idempotent in general: when several input levels merge into output 0, the second pass sees a larger $\mathrm{CDF}_{\min}$ and shifts the levels again.

Color images

Equalizing the red, green, and blue channels independently changes the ratios between channels and therefore shifts hues: a gray pixel can become tinted. The usual approach is to convert to a color space that separates brightness from chromaticity, such as CIELAB, YCbCr, or HSV, equalize only the lightness or luma channel, and convert back. Some colors then fall outside the RGB gamut and must be clipped.

Limitations

• No spatial awareness. A single global mapping cannot enhance a dark region and a bright region at the same time if their intensities overlap with other parts of the image.
• Not a quality criterion. A flat histogram is not an aim in itself; the right histogram depends on the scene (see image histogram).
• Quantization artifacts. Merging levels and leaving gaps produces banding in smooth gradients, especially in 8-bit images with narrow input histograms.
• Not for measurement. The mapping is nonlinear and image-specific, so equalized values lose any physical (radiometric) meaning.

Variants

Adaptive histogram equalization (AHE)

AHE replaces the global mapping by local ones: each pixel is mapped using the histogram of a contextual region around it, so contrast is enhanced relative to the local surroundings. Pizer et al. (1987) reviewed AHE and its variations for medical images and identified two problems with the basic method: it is slow, and it over-enhances noise in relatively homogeneous regions. (In a flat region the local histogram is a narrow peak, which the mapping stretches across the output range.) They described interpolated AHE to speed it up and concluded that clipped AHE should be a preferred method in medical imaging, fast enough for routine display.

Contrast-limited AHE (CLAHE)

CLAHE combines the interpolated and clipped forms. Zuiderveld (1994) published a widely used C implementation in Graphics Gems IV; in it:

• Tiles and interpolation. The image is divided into a grid of non-overlapping contextual regions (tiles), and one equalization mapping is computed per tile. Each pixel’s output is a bilinear interpolation of the mappings of the four nearest tile centers (two or one near the image border), which removes the block boundaries that independent per-tile equalization would produce. This costs $O(N)$ for the interpolation plus one histogram per tile, instead of one histogram per pixel.
• Clip limit. Before accumulating a tile’s histogram, every bin is clipped at a limit, and the clipped excess is redistributed evenly among the bins without pushing any bin above the limit. The limit is specified as a normalized multiple $c$ of the average bin count (tile pixels divided by the number of bins).

Clipping bounds the amplification. With $L$ bins, the slope of a tile’s mapping at level $k$ is $(L-1)\,h'(k)/n$, where $h'$ is the clipped histogram and $n$ the tile’s pixel count. In Zuiderveld’s scheme no bin exceeds the limit $c\,n/L$ after redistribution, so, before rounding,

$$T(k) - T(k-1) \le \frac{L-1}{L}\, c < c,$$

and no level is stretched by more than a factor $c$; noise in flat regions is no longer magnified without bound. A very large $c$ recovers AHE.

The example below implements this scheme in a simplified form (a single redistribution pass, which can leave bins up to $E/L$ above the limit, where $E$ is the clipped excess, so the slope bound loosens to $c + 1$; edge tiles extended to the border) and applies it to an image with a strong left-to-right illumination ramp and a faint stripe texture. Texture contrast is measured as the column-wise standard deviation, since the ramp is constant down each column.

import numpy as np

rng = np.random.default_rng(1)

def equalize(im, L=256):
    cdf = np.cumsum(np.bincount(im.ravel(), minlength=L)) / im.size
    return np.round((L - 1) * cdf).astype(np.uint8)[im]

def tile_maps(im, tiles, clip, L=256):
    """Per-tile equalization lookup tables, with optional histogram clipping."""
    th, tw = im.shape[0] // tiles, im.shape[1] // tiles
    maps = np.empty((tiles, tiles, L))
    for i in range(tiles):
        for j in range(tiles):
            t = im[i * th:(i + 1) * th, j * tw:(j + 1) * tw]
            hist = np.bincount(t.ravel(), minlength=L).astype(float)
            if clip is not None:                       # clip limit as a multiple of
                limit = clip * t.size / L              # the mean count per bin
                excess = np.maximum(hist - limit, 0).sum()
                hist = np.minimum(hist, limit) + excess / L   # redistribute evenly
            maps[i, j] = (L - 1) * np.cumsum(hist) / hist.sum()
    return maps, th, tw

def adaptive_equalize(im, tiles=8, clip=None):
    """AHE (clip=None) or CLAHE with bilinear interpolation between tile maps."""
    maps, th, tw = tile_maps(im, tiles, clip)
    # Continuous tile coordinates of each pixel; tile centers sit at integers.
    fy = np.clip((np.arange(im.shape[0]) + 0.5) / th - 0.5, 0, tiles - 1)
    fx = np.clip((np.arange(im.shape[1]) + 0.5) / tw - 0.5, 0, tiles - 1)
    y0 = np.minimum(fy.astype(int), tiles - 2); wy = (fy - y0)[:, None]
    x0 = np.minimum(fx.astype(int), tiles - 2); wx = (fx - x0)[None, :]
    Y0, X0 = y0[:, None], x0[None, :]
    m = lambda dy, dx: maps[Y0 + dy, X0 + dx, im]      # map of a neighboring tile
    out = ((1 - wy) * ((1 - wx) * m(0, 0) + wx * m(0, 1))
           + wy * ((1 - wx) * m(1, 0) + wx * m(1, 1)))
    return np.round(out).astype(np.uint8)

# Test image: a strong left-to-right illumination ramp (40 -> 200) carrying a
# faint horizontal-stripe texture (amplitude 4) and mild noise.
h = w = 256
yy, xx = np.mgrid[0:h, 0:w]
img = 40 + 160 * xx / (w - 1) + 4 * np.sin(2 * np.pi * yy / 16)
img = np.clip(np.round(img + rng.normal(0, 1, (h, w))), 0, 255).astype(np.uint8)

def texture_contrast(im):
    # The ramp is constant down each column, so the column-wise standard
    # deviation measures the stripe texture (plus noise).
    return im.astype(float).std(axis=0).mean()

results = {
    "input": img,
    "global HE": equalize(img),
    "AHE (no clip)": adaptive_equalize(img, clip=None),
    "CLAHE clip=2": adaptive_equalize(img, clip=2.0),
    "CLAHE clip=4": adaptive_equalize(img, clip=4.0),
}
for name, im in results.items():
    print(f"{name:<14} texture contrast {texture_contrast(im):5.1f}"
          f"   mean of left/right quarters {im[:, :64].mean():5.1f} / {im[:, -64:].mean():5.1f}")

Output:

input          texture contrast   3.0   mean of left/right quarters  59.8 / 180.2
global HE      texture contrast   4.7   mean of left/right quarters  32.8 / 223.7
AHE (no clip)  texture contrast  27.9   mean of left/right quarters 120.1 / 144.4
CLAHE clip=2   texture contrast   7.8   mean of left/right quarters  73.8 / 172.0
CLAHE clip=4   texture contrast  12.4   mean of left/right quarters  86.0 / 164.7

The ramp already spreads the global histogram, so global equalization is nearly linear and barely changes the faint texture. Unclipped AHE magnifies the texture (and the noise) about ninefold but almost erases the large-scale ramp. CLAHE sits between the two, and the clip limit sets the trade-off: the texture gain stays below the bound $c + 1$, while more of the global brightness ordering survives.

Zuiderveld’s code uses 128 histogram bins when none is specified, requires at least 2 regions per direction, and advises against more than 256 contextual regions. OpenCV’s createCLAHE defaults to an 8 × 8 tile grid and a clip limit of 40, normalized the same way.

Where It Is Used

• Medical imaging. AHE and CLAHE were developed and evaluated largely for medical image display, where useful contrast is often subtle and local; Pizer et al. (1987) recommended clipped AHE for routine use there.
• Preprocessing. Equalization and CLAHE normalize contrast before thresholding, feature detection, or recognition, reducing sensitivity to lighting; CLAHE on the luminance channel is a common step for low-light images.
• Visualization. High-bit-depth data, such as thermal, microscopy, and satellite imagery, must be mapped to 8 bits for display; equalization provides an automatic mapping.
• Normalization across acquisitions. Histogram matching aligns intensity distributions of images from different sensors or sessions before comparison.