Median Filter

The median filter replaces each pixel of an image by the median of the pixel values in a small window around it. It is the standard remedy for impulse (salt-and-pepper) noise: isolated pixels with extreme values are discarded outright instead of being averaged into their neighbors, and a clean step edge passes through unchanged. It is the best-known nonlinear image filter and the simplest rank, or order-statistic, filter.

Running medians are also a smoothing tool for data series in statistics, where John Tukey popularized them; his book Exploratory Data Analysis (1977) describes smoothers built from repeated running medians of three. In image processing, much of the work on the filter has concerned its cost, which grows with the window size unless a specialized algorithm is used.

Problem

Linear smoothing, such as a box or Gaussian filter, spreads an impulse of value 255 among pixels near 100 into a visible blob, and blurs edges. The median filter discards such outliers and creates no intermediate values across a step.

Inputs and Outputs

The input is a grayscale image and a window, usually a $k \times k$ square with odd $k$ (3, 5, 7), sometimes a cross or a disk. The output is an image of the same size in which each pixel is the median of the input values in the window centered on it. Every output value is one of the input values in its window, so no new gray levels appear.

Intuition

The median of a set of numbers is the middle value after sorting, so making the largest value much larger does not change it. The median stays within the range of the uncorrupted values as long as fewer than half of the values are corrupted: its breakdown point is 50%, while that of the mean is zero, since one arbitrarily large value can move the mean arbitrarily far.

In a $3 \times 3$ window, up to four of the nine pixels can be impulses and the output is still one of the clean values. Near a straight step edge the pixels on the center’s side of the edge are in the majority, so the median picks a value from that side and the edge stays sharp.

Algorithm

For each output pixel $[m, n]$ :

Collect the $N$ input values in the window centered on $[m, n]$ , extending the image at the border.
Find the median of these values: the $\lceil N/2 \rceil$ -th smallest for odd $N$ .
Write it to $g[m, n]$ .

Every window sees the original input values. Feeding already-filtered values into later windows gives a different filter, the recursive median filter.

Mathematical Formulation

With a window $W$ of offsets,

g[m, n] = \operatorname{median}\{\, f[m + i,\, n + j] : (i, j) \in W \,\}.

The median is also the value that minimizes the sum of absolute deviations,

g[m, n] = \arg\min_{v} \sum_{(i, j) \in W} \big| f[m + i,\, n + j] - v \big| ,

whereas the local mean minimizes the sum of squared deviations, which lets large deviations dominate.

Parameters

Window size. A $3 \times 3$ window removes isolated impulses and changes the image least. Larger windows tolerate higher impulse densities but remove larger details and flatten texture.
Window shape. A square window rounds the corners of shapes; a cross-shaped (plus) window preserves corners and horizontal and vertical lines better but removes less noise.
Iterations. Repeated application removes more noise and drives the image toward a root signal (below).

Complexity

The naive method sorts the $k^2$ values of each window, which costs $O(k^2 \log k)$ per pixel; a linear-time selection algorithm reduces this to $O(k^2)$ . A box filter, by contrast, costs $O(1)$ per pixel with running sums.

Huang, Yang, and Tang (1979) exploited the overlap between neighboring windows in an 8-bit image by keeping a 256-bin histogram of the window. Moving one pixel to the right, the algorithm removes the $k$ values of the leaving column, adds the $k$ values of the entering one, and finds the new median by stepping through the histogram from the previous one. The authors give the running time as $O(n)$ per pixel for an $m \times n$ window.

Perreault and Hébert (2007) removed the dependence on window size. They keep one histogram per image column, covering the $k$ rows of the current window; moving down a row updates each column histogram by removing one pixel and adding one. Along a row, the window histogram gains the histogram of the entering column and loses that of the leaving column. Both steps cost a number of operations set by the number of histogram bins rather than by $k$ , so the cost per pixel is $O(1)$ . The constant is large, since each histogram addition touches every bin, so the authors use vectorized (SIMD) arithmetic and multilevel histograms; they also extend the method to higher-precision data.

Implementation

The example adds salt-and-pepper noise at two densities, and Gaussian noise, to a synthetic 8-bit image with flat regions, edges, and a one-pixel-wide bright line. It compares median filters with a $3 \times 3$ box filter and a Gaussian, reporting PSNR against the clean image, the standard deviation left inside a flat square, and the contrast of the thin line (105 gray levels when clean).

import numpy as np
from scipy import ndimage

rng = np.random.default_rng(0)

# Synthetic 8-bit test image: a horizontal ramp, a bright square, a dark disk
# and a one-pixel-wide bright vertical line at column 230.
yy, xx = np.mgrid[0:256, 0:256]
img = 60.0 + 100.0 * xx / 255
img[40:120, 40:120] = 220.0
img[(yy - 180) ** 2 + (xx - 170) ** 2 < 40 ** 2] = 20.0
img[:, 230] = 255.0

def psnr(a, b):
    return 10 * np.log10(255.0 ** 2 / np.mean((a - b) ** 2))

def salt_and_pepper(x, fraction):
    out = x.copy()
    hit = rng.random(x.shape) < fraction
    out[hit] = rng.choice([0.0, 255.0], size=hit.sum())
    return out

def run(noisy, label):
    print(label)
    print(f"  {'filter':15s} {'PSNR':>7s} {'flat std':>9s} {'line':>6s}")
    for name, out in [
        ("noisy input",    noisy),
        ("median 3x3",     ndimage.median_filter(noisy, size=3)),
        ("median 5x5",     ndimage.median_filter(noisy, size=5)),
        ("box 3x3",        ndimage.uniform_filter(noisy, size=3)),
        ("Gaussian s=1",   ndimage.gaussian_filter(noisy, sigma=1.0)),
    ]:
        flat = out[50:110, 50:110].std()          # inside the square
        line = (out[10:30, 230] - out[10:30, 226]).mean()  # true contrast 105
        print(f"  {name:15s} {psnr(out, img):5.1f}dB {flat:9.1f} {line:6.1f}")

run(salt_and_pepper(img, 0.10), "10% salt-and-pepper noise:")
run(salt_and_pepper(img, 0.40), "40% salt-and-pepper noise:")
run(img + rng.normal(0.0, 15.0, img.shape), "Gaussian noise, sigma = 15:")

Output:

10% salt-and-pepper noise:
  filter             PSNR  flat std   line
  noisy input      15.3dB      50.3  121.2
  median 3x3       30.4dB       5.2    7.2
  median 5x5       30.8dB       0.0    2.0
  box 3x3          23.3dB      16.6   49.9
  Gaussian s=1     24.5dB      14.1   54.6
40% salt-and-pepper noise:
  filter             PSNR  flat std   line
  noisy input       9.3dB      92.1   74.5
  median 3x3       18.7dB      26.6   22.9
  median 5x5       27.5dB       0.0    2.0
  box 3x3          17.0dB      29.8   46.5
  Gaussian s=1     17.9dB      25.0   44.5
Gaussian noise, sigma = 15:
  filter             PSNR  flat std   line
  noisy input      24.6dB      15.2  110.4
  median 3x3       29.3dB       6.0   13.7
  median 5x5       30.3dB       3.6    3.9
  box 3x3          29.0dB       5.0   40.3
  Gaussian s=1     30.1dB       4.3   46.1

At 10% impulse density, both median filters raise the PSNR by about 15 dB, while the linear filters gain 8 to 9 dB and leave the flat square speckled, because each impulse is spread over its neighbors rather than removed. The median fails only where impulses of one kind fill more than half of the window. At 40% density (20% of each kind) that happens in a few percent of $3 \times 3$ windows but in under 0.1% of $5 \times 5$ windows, and the $5 \times 5$ filter restores the flat region exactly. The thin line shows the price: the median filters cut its contrast from 105 to between 2 and 23, while the linear filters keep a blurred trace of 40 to 55. With Gaussian noise, the $3 \times 3$ median leaves more noise in the flat region than the $3 \times 3$ box filter (6.0 against 5.0), as expected from the two estimators.

Properties and Behavior

Nonlinear. The median of a sum is not the sum of medians, so the filter is not a convolution and has no frequency response. It is shift-invariant: the same rule applies everywhere.
Not separable. A 1D median along rows followed by one along columns is cheaper, but its result generally differs from the 2D median of the square window.
Preserves step edges. A noise-free straight edge, or any monotone ramp in 1D, passes through unchanged.
Removes small structures. Any feature that covers less than half of the window is replaced by its surroundings. In a $3 \times 3$ window a one-pixel line covers three of nine pixels and disappears, and the corner of a bright square covers four of nine and is rounded off.
Root signals. Applying the filter twice can change the image again, so it is not idempotent. For 1D signals and a specific form of the filter, Gallagher and Wise (1981) derived necessary and sufficient conditions for a signal to be left unchanged, a root: it must be locally monotone. They also proved that repeated filtering of a signal of length $L$ reaches a root after at most $(L - 2)/2$ passes.

Limitations

Gaussian noise. For independent Gaussian noise the median is a less efficient estimator than the mean: for large windows its noise variance is about $\pi/2 \approx 1.57$ times that of a box filter of the same size. A linear filter of the same size, or an edge-preserving filter such as the bilateral filter, is often preferable when the noise is not impulsive.
Loss of fine detail. Small structures and fine texture are removed with the noise, and large windows produce flat, patchy regions.
High impulse densities. Once impulses approach half of a window the filter fails, and enlarging the window costs more detail.
Filtering every pixel. Uncorrupted pixels are changed too, slightly degrading clean regions; switching and adaptive variants avoid this.
Color. Filtering each channel separately can output a color that appears nowhere in the window (see the vector median below).

Variants

Weighted and center-weighted median. Each sample is counted according to an integer weight, as if repeated; giving the center pixel extra weight preserves more detail at the cost of less noise removal.
Adaptive median. Hwang and Haddad (1995) proposed two adaptive median filters with variable window size; one tests whether the center pixel is an impulse, the other estimates the size of the impulse noise, and both aim to remove dense impulse noise while preserving sharpness.
Switching median. An impulse detector selects the corrupted pixels, and only those are replaced by the median, leaving clean pixels untouched.
Vector median. Astola, Haavisto, and Neuvo (1990) introduced vector median filters, which process vector-valued samples, such as color pixels, as vectors rather than component by component. The vector median outputs the sample in the window whose total distance to all the others is smallest, so the result is always a color present in the input.
Other rank filters. Other order statistics give the minimum and maximum filters (grayscale erosion and dilation) and percentile filters.