Frequency Domain Filtering

Frequency domain filtering designs a linear, shift-invariant filter by its transfer function: a gain $H$ for every spatial frequency, by which the image’s spectrum is multiplied before it is transformed back. By the convolution theorem this is the same operation as convolution with the inverse transform of $H$, so the frequency domain is less a different kind of filtering than a different place to specify and compute it. It is the natural place when the property that matters is stated in frequencies (pass everything below a cutoff, remove one periodic pattern, keep a band of scales), when the kernel would be large in space, and when a product of two images must be split into a sum, as homomorphic filtering does for illumination and reflectance.

The Fourier transform article defines the transform and introduces the ideal, Gaussian, high-pass and notch filters in outline, and image filtering places them among the other filters. This article covers the design of the filters themselves: how a transfer function is sampled on the discrete Fourier grid and what that sampling assumes, how the ideal, Gaussian and Butterworth low-pass filters trade ringing, sharpness and passed noise against each other at the same cutoff, why “the Butterworth filter” has two definitions in practice, and how filtering the logarithm of an image separates a slowly varying illumination from the reflectance it multiplies.

Purpose

A spatial kernel says how each output pixel mixes its neighbours; a transfer function says what happens to each grating the image is made of. Many design goals are only simple in the second form. A low-pass filter for noise reduction should keep coarse structure and remove fine structure above some frequency; a filter against scan lines should remove two spikes in the spectrum and nothing else, as the notch example in image processing shows; a band-pass filter should select a range of scales. Stating the gain directly makes the trade-offs visible: how sharp the transition between pass and stop is, how much the passband is attenuated, and how much of the stopband leaks through.

The frequency domain is also the cheaper place to apply a large kernel, since the cost of the fast Fourier transforms does not depend on the kernel’s size (image processing gives the operation counts). Briand and Vacher (2016) note that a filter defined by a continuous function of frequency covers every image size with one formula, and call such definitions a classical approach to building multi-scale structures such as the steerable pyramid.

Homomorphic filtering extends the method to a nonlinear problem. An image formed as illumination times reflectance is not a sum of a smooth and a detailed part, so no linear filter can separate them directly. Its logarithm is such a sum, and a high-pass filter applied to the logarithm suppresses the slowly varying term.

When to Use It

Frequency domain filtering helps when:

• the filter is specified by its response: a cutoff, a band, a notch at known frequencies, or a measured spectrum as in the Wiener filter;
• the kernel is large, so that the FFT is cheaper than direct convolution;
• the disturbance is periodic and concentrated in a few spectral peaks;
• a multiplicative component varies slowly, as uneven illumination or the bias field of an MR image does, and can be moved into a sum by a logarithm.

It is the wrong tool when:

• the kernel is small: a $3 \times 3$ or separable Gaussian filter is cheaper in space;
• the filter is nonlinear (median, bilateral, morphological), since such filters have no transfer function (image filtering);
• the image’s borders do not match and cannot be padded, since the discrete transform treats the image as one period of a tiling (see Pitfalls);
• the filter must adapt to the content. A transfer function is the same everywhere in the image, so it smooths edges and flat regions alike; edge-preserving filters and the base–detail decompositions of tone mapping exist for that reason.

How It Works

The algorithm

For an $M \times N$ image $f$ and a transfer function $H$, the filtered image is

$$g = \mathcal{F}^{-1}\{ H \cdot \mathcal{F}\{f\} \},$$

where $\mathcal{F}$ is the discrete Fourier transform (DFT) and the product is taken element by element on the DFT grid. With NumPy’s ordering, the DFT bin $(k, l)$ holds the frequency $(u, v) = (k/M, l/N)$ cycles per pixel, with indices above half the size standing for negative frequencies (np.fft.fftfreq returns these values). A radial filter is a function of $D = \sqrt{u^2 + v^2}$, evaluated at every bin. If $H$ is real and symmetric, $H(-u, -v) = H(u, v)$, the output of a real image is real, up to rounding. The result is the circular convolution of $f$ with the kernel $h = \mathcal{F}^{-1}\{H\}$: content leaving one border re-enters at the opposite one.

Briand and Vacher (2016) analyse what this procedure computes when $H$ is given as a continuous function on the Nyquist domain. They interpret the filtering as the convolution of a distribution standing for the filter with a trigonometric polynomial that interpolates the image, and obtain three equally licit variants of the same algorithm: sampling $H$ at the DFT frequencies (their method 1, the classical DFT computation), a real convention that gives each coefficient on the Nyquist boundary the average of the filter’s values at $-\pi$ and $\pi$ (method 2), and a windowed variant that sets the boundary coefficients to zero (method 3). The variants differ only in the DFT coefficients on the boundary of the Nyquist domain, which exist only for even image sizes. For odd sizes, or for a filter that is zero on that boundary, all three coincide; for a radial filter, methods 1 and 2 coincide. In their experiments the differences were not visible to the eye, of the order of one grey level for 8-bit images in general, and possibly larger for filters with large values on the boundary.

Low-pass filters at one cutoff

A radial low-pass filter with cutoff $D_0$ passes frequencies with $D < D_0$ and attenuates the rest. Three families differ in how the transition happens.

• Ideal: $H(D) = 1$ for $D < D_0$ and $0$ otherwise. In one dimension its kernel is the sinc $\sin(2\pi D_0 x)/(\pi x)$, whose slowly decaying lobes make a filtered step overshoot and ripple: the Gibbs phenomenon. A unit step becomes $\tfrac12 + \operatorname{Si}(2\pi D_0 x)/\pi$, with $\operatorname{Si}$ the sine integral, whose first maximum, at $x = 1/(2D_0)$, overshoots by $\operatorname{Si}(\pi)/\pi - \tfrac12 \approx 0.0895$ of the step, whatever the cutoff, so a lower cutoff makes the ripples wider, not weaker, as MRI shows for truncated k-space.
• Gaussian: $H(D) = e^{-D^2/(2 D_\sigma^2)}$, the transform of a spatial Gaussian of standard deviation $1/(2\pi D_\sigma)$. Its kernel is positive, so a step stays monotonic (Gaussian filter).
• Butterworth: a family indexed by its order $n$ that runs between the two. scikit-image 0.25.2 (skimage.filters.butterworth) uses by default

$$H_{\text{sq}}(D) = \frac{1}{1 + (D/D_0)^{2n}},$$

with gain $\tfrac12$ at the cutoff, and with squared_butterworth=False its square root,

$$H_{\text{plain}}(D) = \frac{1}{\sqrt{1 + (D/D_0)^{2n}}},$$

with gain $1/\sqrt 2$ at the cutoff. The function’s documentation says that the “Butterworth filter” of image processing textbooks is often the square of the traditional Butterworth filter, and it cites Butterworth’s 1930 paper on filter amplifiers for the traditional form; it defines the cutoff as a fraction of the sampling frequency, between 0 and 0.5 (scikit-image 0.25.2 source). The squared form is the traditional form’s squared magnitude, so a squared Butterworth filter of order $n$ is the plain one of the same order applied twice. Both tend to the ideal filter as $n \to \infty$, and both fall off as $D^{-2n}$ or $D^{-n}$ far above the cutoff.

The two conventions give different filters for the same pair $(D_0, n)$, and the difference is not small: at $D = 2 D_0$ and $n = 2$ the squared gain is $1/17 \approx 0.059$ and the plain gain $1/\sqrt{17} \approx 0.24$. A cutoff and an order therefore do not specify a filter until the gain at the cutoff is stated.

The noise a low-pass filter passes follows from Parseval’s theorem: white noise of variance $\sigma^2$ leaves the filter with variance $\sigma^2 \cdot \frac{1}{MN} \sum_{k,l} |H_{kl}|^2$. In the continuous plane, our derivation (substituting $t = (D/D_0)^2$ into $\int_0^\infty 2\pi D\, |H(D)|^2\, dD$ and using the Beta integral) gives, relative to the ideal filter’s $\pi D_0^2$,

$$\frac{\int |H_{\text{sq}}|^2}{\pi D_0^2} = \frac{1}{n}\left(1 - \frac{1}{n}\right) \frac{\pi}{\sin(\pi/n)}, \qquad \frac{\int |H_{\text{plain}}|^2}{\pi D_0^2} = \frac{\pi/n}{\sin(\pi/n)} ,$$

with the limit 1 for the squared form at $n = 1$ and a divergent integral for the plain form at $n = 1$, whose $|H|^2$ decays only as $D^{-2}$. The squared form passes least noise at $n = 2$ ($\pi/4 \approx 0.785$ of the ideal’s) and approaches the ideal’s from below as $n$ grows ($0.833$ at $n = 4$); the plain form passes more than the ideal at every order ($\pi/2$ at $n = 2$) and approaches it from above. For a Gaussian, $\int |H|^2 = \pi D_\sigma^2$; gain $\tfrac12$ at $D_0$ means $D_0^2 = 2 \ln 2\, D_\sigma^2$, so it passes $1/(2 \ln 2) \approx 0.721$ of the ideal’s, and gain $1/\sqrt2$ ($D_0^2 = \ln 2\, D_\sigma^2$) passes $1/\ln 2 \approx 1.443$. On a finite grid the integrals stop at the Nyquist frequency, which matters only for slowly decaying responses.

High-pass, band and notch filters

As in Fourier transform, every low-pass filter $H$ defines a high-pass filter $1 - H$. scikit-image computes its high-pass Butterworth as $(D/D_0)^{2n}/(1 + (D/D_0)^{2n})$, which is $1 - H_{\text{sq}}$, and in the plain convention takes its square root, which is not $1 - H_{\text{plain}}$ (scikit-image 0.25.2 source). Adding a multiple of the high-pass output back to the image sharpens it, the mechanism of image sharpening. A band-pass filter is a product or difference of low-pass filters at two cutoffs, a band-reject filter one minus a band-pass, and a notch filter a band-reject filter centred on isolated frequencies and their mirror images. A transfer function of the form $\gamma_L + (\gamma_H - \gamma_L)(1 - H)$, with $\gamma_L < 1 < \gamma_H$, attenuates low frequencies without removing them and amplifies high ones; a homomorphic filter can take this form.

Homomorphic filtering

Under the model $I(x, y) = L(x, y)\, r(x, y)$, with illumination $L$ and reflectance $r$, the logarithm turns the product into a sum,

$$\log I = \log L + \log r .$$

If $\log L$ varies slowly and $\log r$ carries the fine structure, a high-pass transfer function $H_{\text{hp}}$ applied to $\log I$ suppresses the first and keeps most of the second, and exponentiation returns to intensities:

$$\hat r = \exp\!\big( \mathcal{F}^{-1}\{ H_{\text{hp}} \cdot \mathcal{F}\{\log I\} \} \big).$$

The filter acts linearly on the logarithm, so the residual of each component can be computed separately: the output is $H_{\text{hp}} * \log L + H_{\text{hp}} * \log r$. With a Gaussian low-pass $G$ and $H_{\text{hp}} = 1 - G$, the log-domain output is $\log I - G * \log I$.

Retinex algorithms reach a similar operator from a model of lightness perception. Petro, Sbert and Morel (2014) write the single-scale retinex as $R = \log I - \log(I * F)$, with a normalized surround $F$ such as a Gaussian, and show that under the assumption $I = S r$ with a smoothly varying illumination $S$ it is approximately $\log(r/\bar r)$, the logarithm of reflectance relative to its local weighted average, and hence independent of the illumination. They note that moving the logarithm inside the convolution, $\log I - (\log I) * G$, gives a classical high-pass linear filter applied to $\log I$ rather than $I$. That second form is the Gaussian homomorphic filter above. The two differ in the mean they subtract, arithmetic for the retinex and geometric for the homomorphic filter, and Jobson, Rahman and Woodell’s 1997 study of the centre/surround retinex, as Petro et al. report, found experimentally that placing the logarithm after the surround was better. Morel, Petro and Sbert (2014, §4.1), describing a classification of methods for correcting intensity inhomogeneity in MRI, list homomorphic filtering and unsharp masking among the filtering methods, which treat the non-uniform illumination as a low-frequency artefact and remove it by subtracting a low-pass version of the image. The N4 bias-field correction described in medical image preprocessing also works on log intensities, where the field is additive, but estimates the field iteratively instead of with a fixed filter.

How to Apply It

Example: Butterworth order as a continuum

The first script compares low-pass filters of one cutoff, $D_0 = c = 0.05$ cycles per pixel, on a $256 \times 256$ image of a disc of radius 64 pixels with value 1 on a background of 0. Each pixel holds the fraction of its area inside the disc ($8 \times 8$ subsamples), so the edge is sampled without staircase artefacts. Every filter multiplies the full complex DFT by a real radial transfer function sampled at the DFT frequencies, Briand and Vacher’s method 1, which for these radial filters coincides with their method 2; the filtering is circular, and the disc lies 64 pixels from every border. The filters are the ideal low-pass filter, a Gaussian with gain $\tfrac12$ at $c$, a Gaussian with gain $1/\sqrt2$ at $c$, and Butterworth filters of orders 1, 2, 4, 8, 16 and 32 in both conventions, squared (scikit-image’s default, gain $\tfrac12$ at $c$) and plain (gain $1/\sqrt2$). Since the sources leave open which gain defines “the” cutoff, the script prints both conventions, each with a Gaussian matched to the same gain.

For each filter it prints:

• the overshoot, the largest value minus 1, and the undershoot, minus the smallest value, both in percent of the step and measured in the ring of pixels within 16 pixels of the edge (a value of zero or below means none);
• the 10–90% width of the edge on the centre row, searched from 16 pixels inside the edge outward, with linear interpolation between pixels;
• the variance of filtered white noise relative to that of the ideal filter, $\sum |H|^2 / \sum |H_{\text{ideal}}|^2$, the expected value from Parseval’s theorem;
• for one field of white Gaussian noise of standard deviation 0.1, the measured root-mean-square of the filtered noise divided by its expected value $0.1 \sqrt{\operatorname{mean}|H|^2}$, with a tolerance of four standard deviations of that ratio, $4 \cdot \tfrac12 \sqrt{2 \sum |H|^4} / \sum |H|^2$ (our derivation for real white noise, whose DFT coefficients come in conjugate pairs).

It also checks the data (the disc’s area against $\pi R^2$, the noise’s standard deviation), the agreement of the formulas with skimage.filters.butterworth in both conventions, and that every filter preserves the image’s sum, since $H(0) = 1$.

The measures are restricted to the edge ring for a reason of construction: the ripples a filter makes along a circular edge are circles concentric with the disc, so the ripples from every part of the edge reach its centre at the same distance and add up there. A search starting from the centre therefore measures that focusing rather than the edge itself, and the first version of the planning script, which measured overshoot over the whole image and searched for the edge from the centre, failed for the ideal filter for this reason; the ring was adopted before the example was first run.

The plan, fixed before the first run, recorded four checks and four hypotheses. The hypotheses: (HA1) the order-1 filters of both conventions do not overshoot (at most 0.01%), overshoot grows strictly with the order in both conventions, the order-32 filters stay below the ideal filter’s overshoot and within 1.5 percentage points of it, and neither Gaussian overshoots (below 0.01%); (HA2) the noise ratios follow the derivation above: the squared form is lowest at $n = 2$ and rises after it, the plain form falls strictly with the order, the Gaussian with gain $\tfrac12$ is below every squared Butterworth, and the Gaussian with gain $1/\sqrt2$ lies between the plain filters of orders 2 and 4; (HA3) the squared Butterworth widths for $n \ge 2$ are within 0.25 pixel of the ideal filter’s and the order-1 width is more than 1 pixel wider, the plain widths increase strictly with the order and all stay below the ideal’s, the Gaussian with gain $\tfrac12$ is wider than every squared Butterworth with $n \ge 2$ and narrower than order 1, and the Gaussian with gain $1/\sqrt2$ is within 0.05 pixel of the plain order-2 width and does not overshoot; (HA4) for every filter that overshoots by more than 0.5%, the overshoot inside the disc exceeds the undershoot outside it. All four hypotheses were computed with a planning script on the same set-up before the pilot, so the example illustrates them rather than tests them; only the noise realization is random. Two of them were set in view of those computations. A first version of HA3, written from theory, predicted that the edge width would shrink with the order in both conventions; the planning computation contradicted it, and HA3 was replaced by the computed predictions above. The 1.5-point bar of HA1 was set with the plain order-32 filter’s computed gap to the ideal filter in view, after a first bar of 1 point, which that filter could not meet. The script uses NumPy 2.3.5, SciPy 1.16.3 (through scikit-image) and scikit-image 0.25.2 on one BLAS thread, runs in under a second, and printed identical output on two runs.

import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"          # one BLAS thread, before NumPy
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
from skimage.filters import butterworth

N, R, c, sigma = 256, 64.0, 0.05, 0.1             # size, disc radius, cutoff (cycles/px), noise
k = np.fft.fftfreq(N)                             # DFT frequencies, cycles per pixel
F = np.hypot(k[:, None], k[None, :])              # radial frequency of every DFT bin

# A disc of value 1 on 0, each pixel holding the fraction of its area inside (8 x 8 subsamples).
yy, xx = np.mgrid[0:N, 0:N].astype(float)
sub = (np.arange(8) + 0.5) / 8 - 0.5
disc = sum(((yy + dy - 128) ** 2 + (xx + dx - 128) ** 2 < R ** 2).astype(float)
           for dy in sub for dx in sub) / 64
noise = sigma * np.random.default_rng(0).standard_normal((N, N))

def apply(H, img):                                # filter = multiply the DFT by H
    return np.real(np.fft.ifft2(np.fft.fft2(img) * H))

def butter(n, squared):                           # squared: gain 1/2 at c; else 1/sqrt(2)
    H = 1 / (1 + (F / c) ** (2 * n))
    return H if squared else np.sqrt(H)

filters = {"ideal": (F < c).astype(float),
           "Gaussian, H(c)=0.5": np.exp(-np.log(2) * (F / c) ** 2),
           "Gaussian, H(c)=0.707": np.exp(-0.5 * np.log(2) * (F / c) ** 2)}
orders = [1, 2, 4, 8, 16, 32]
for sq, tag in [(True, "squared"), (False, "plain")]:
    for n in orders:
        filters[f"Butterworth {tag} n={n}"] = butter(n, sq)

rad = np.hypot(yy - 128, xx - 128)
ring = np.abs(rad - R) <= 16                      # the edge zone where ringing is measured

def width(out):                                   # 10-90% width on the centre row, left edge
    p = out[128, 80::-1]                          # from x = 80 (inside the disc) outward
    i = np.argmax(p < 0.9)
    j = i + np.argmax(p[i:] < 0.1)
    x90 = 80 - (i - 1) - (p[i - 1] - 0.9) / (p[i - 1] - p[i])
    x10 = 80 - (j - 1) - (p[j - 1] - 0.1) / (p[j - 1] - p[j])
    return x90 - x10

print(f"data: disc area / (pi R^2) - 1 = {disc.sum() / (np.pi * R ** 2) - 1:+.2e} (|.| < 1e-3), "
      f"noise std / sigma - 1 = {noise.std() / sigma - 1:+.2e} (|.| < 0.011)")
lib = max(np.abs(butterworth(disc, c, high_pass=False, order=n, squared_butterworth=sq)
                 - apply(butter(n, sq), disc)).max() for n in orders for sq in (True, False))
print(f"check: max |skimage butterworth - formula| = {lib:.1e} (< 1e-12)")

ideal_var = np.mean(filters["ideal"] ** 2)
print(f"\n{'filter':26s} {'over %':>8s} {'under %':>8s} {'10-90 px':>9s}"
      f" {'noise var/ideal':>15s} {'rms meas/exp':>12s} {'tol':>7s}")
dsum = 0.0
for name, H in filters.items():
    out = apply(H, disc)
    dsum = max(dsum, abs(out.sum() - disc.sum()))
    var = np.mean(H ** 2)                         # expected noise variance / sigma^2 (Parseval)
    rms = np.sqrt(np.mean(apply(H, noise) ** 2)) / (sigma * np.sqrt(var))
    tol = 4 * 0.5 * np.sqrt(2 * np.sum(H ** 4)) / np.sum(H ** 2)  # 4 sd of the rms ratio
    print(f"{name:26s} {100 * (out[ring].max() - 1):+8.4f} {-100 * out[ring].min():+8.4f}"
          f" {width(out):9.4f} {var / ideal_var:15.4f} {rms:12.4f} {tol:7.4f}")
print(f"\ncheck: max |sum(filtered) - sum(disc)| = {dsum:.1e} (< 1e-8)")

data: disc area / (pi R^2) - 1 = +5.14e-05 (|.| < 1e-3), noise std / sigma - 1 = -5.60e-04 (|.| < 0.011)
check: max |skimage butterworth - formula| = 8.9e-16 (< 1e-12)

filter                       over %  under %  10-90 px noise var/ideal rms meas/exp     tol
ideal                      +11.0007  +8.1976    9.1864          1.0000       0.9426  0.1254
Gaussian, H(c)=0.5          -0.0012  -0.0009    9.7047          0.7295       0.9561  0.1038
Gaussian, H(c)=0.707        -0.0000  -0.0000    6.9908          1.4589       0.9760  0.0734
Butterworth squared n=1     -0.3907  -0.2877   10.3558          1.0030       0.9733  0.0726
Butterworth squared n=2     +3.7910  +2.9870    9.0873          0.7942       0.9508  0.1112
Butterworth squared n=4     +7.5054  +6.2085    8.9605          0.8424       0.9427  0.1223
Butterworth squared n=8     +9.1537  +7.6323    8.9718          0.9080       0.9393  0.1247
Butterworth squared n=16    +9.8047  +8.0355    9.0157          0.9542       0.9404  0.1250
Butterworth squared n=32   +10.3617  +8.2299    9.1222          0.9804       0.9428  0.1251
Butterworth plain n=1       -0.1254  -0.0925    6.7619          4.8877       0.9955  0.0257
Butterworth plain n=2       +2.6626  +2.1540    7.0117          1.5802       0.9760  0.0707
Butterworth plain n=4       +6.1879  +5.2106    7.6338          1.1232       0.9597  0.1024
Butterworth plain n=8       +8.4912  +7.1729    8.2140          1.0377       0.9533  0.1151
Butterworth plain n=16      +9.4411  +7.9432    8.6125          1.0177       0.9497  0.1203
Butterworth plain n=32      +9.9337  +8.3395    8.8774          1.0116       0.9472  0.1227

check: max |sum(filtered) - sum(disc)| = 1.8e-12 (< 1e-8)

All four checks held. The disc’s area is within $10^{-4}$ of $\pi R^2$, the hand-written transfer functions agree with scikit-image’s butterworth to rounding in both conventions, and every filter preserves the sum. The measured noise is within its tolerance for all fifteen filters; every ratio is below 1, from 0.94 to 0.996, which in our reading reflects one noise field shared by all filters that happens to carry slightly less power than expected at the frequencies they pass.

Ringing (HA1 held). The order-1 Butterworth filters do not overshoot at all; in the edge ring their output stays below 1 and above 0; in our reading their long tails spread the edge so widely that the plateau is not reached within 16 pixels. From order 2 on, the overshoot grows strictly with the order in both conventions, from 3.8% to 10.4% for the squared form and from 2.7% to 9.9% for the plain form, and the order-32 filters are 0.64 and 1.07 percentage points below the ideal filter’s 11.0%. Neither Gaussian overshoots. A Butterworth filter is thus a dial between a filter with a positive kernel and the ideal one, and at order 8 it already overshoots by 77% (plain) to 83% (squared) of the ideal filter’s amount.

Ringing on a curved edge (HA4 held). For all eleven filters that ring, the overshoot inside the disc exceeds the undershoot outside it; for the ideal filter, 11.0% against 8.2%. The one-dimensional Gibbs value, 8.95%, lies between the two. In our reading this is the curvature of the edge: inside the disc the ripples from the edge converge, outside they spread, so a straight edge would give the one-dimensional value on both sides.

Edge width (HA3 held). With gain $\tfrac12$ at the cutoff, every squared Butterworth from order 2 on gives practically the same edge, 8.96 to 9.12 pixels against the ideal filter’s 9.19: beyond order 2 a higher order narrows the edge by at most 0.13 pixel, while the overshoot rises from 3.8% to as much as 10.4%. Order 1 is more than a pixel wider. With gain $1/\sqrt2$ the plain filters are sharper the lower their order, from 6.76 pixels at order 1 to 8.88 at order 32; in our reading a low-order plain filter passes much of the band above the cutoff, which narrows the edge. The Gaussians follow the same convention: 9.70 pixels with gain $\tfrac12$, wider than the squared Butterworths of order 2 and up, and 6.99 pixels with gain $1/\sqrt2$, 0.021 pixel narrower than the plain order-2 filter and without its 2.7% overshoot.

Noise (HA2 held). The expected noise ratios follow the continuous derivation, slightly modified by the finite grid. The squared Butterworth passes least noise at order 2, 0.79 of the ideal’s, and more at every other order; the plain form passes 4.9 times the ideal’s at order 1, where the integral is cut off only by the grid, and approaches 1 from above. The Gaussian with gain $\tfrac12$ passes 0.73, less than any squared Butterworth, and the Gaussian with gain $1/\sqrt2$ passes 1.46, between the plain orders 2 and 4.

Taken together, at matched cutoff and gain convention, the Gaussian in this test does at least as well as the Butterworth filter of comparable edge width on both ringing and noise: against squared order 2 it gives up 0.6 pixel of edge width for no overshoot and less noise, and against plain order 2 it matches the width with no overshoot and less noise. In our reading, a Butterworth filter of moderate order is chosen for what the example does not measure: a passband closer to 1 below the cutoff, so that frequencies the application needs keep their amplitude, with a sharper transition than a Gaussian. Order 2 of the squared form is the compromise among the Butterworth filters here: almost the full sharpness of the family, the least noise, and less than half the ideal filter’s overshoot.

Example: homomorphic filtering

The second script builds an image whose illumination and reflectance are known. The illumination varies smoothly by a factor of 50, from 0.02 to 1: $\log L = \mu + \rho\,\big(\cos(2\pi(x - 32)/N) + \cos(2\pi(y - 32)/N)\big)/2$ with $N = 256$, brightest at pixel (32, 32) and darkest at (160, 160). The reflectance is a $64 \times 64$ tile of $8 \times 8$-pixel blocks with values drawn uniformly from 0.2 to 0.9, repeated $4 \times 4$ times, so that the brightest and the darkest $64 \times 64$ tiles hold the same reflectance pattern. The image is $I = L r$; a noisy copy adds Gaussian noise of standard deviation 0.0005, small enough that the darkest pixel is 8.8 standard deviations above zero, so the logarithm is always defined.

The low-pass filter is a Gaussian transfer function $G(D) = e^{-2\pi^2 \sigma^2 D^2}$ with $\sigma = 16$ pixels, applied through the DFT. The illumination and the filter were designed at planning so that the separation shows, and the example illustrates the mechanism rather than testing it. The illumination is made of cosines with one cycle across the image, which are eigenfunctions of every radial filter: the high-pass $1 - G$ multiplies them by exactly $1 - e^{-2\pi^2\sigma^2/N^2} = 0.0742$. The reflectance repeats every 64 pixels, so its spectrum lies on multiples of $1/64$ cycles per pixel, where the high-pass gain is $1 - e^{-2\pi^2\sigma^2/64^2} \approx 0.71$ and higher. The two components therefore share only the zero frequency, which the high-pass removes, and the output of the homomorphic filter is exactly $0.0742\,(\log L - \mu)$ plus the high-passed log reflectance, which the script checks. The two parts are orthogonal, so $\beta$ below is exactly 0.0742 for the homomorphic filter and 1 for $\log I$.

Five outputs are compared: the image $I$ and its logarithm, unfiltered, as references; the high-pass on $I$ itself, $I - G * I$; the homomorphic filter, $\log I - G * \log I$; and the single-scale retinex of Petro et al., $\log I - \log(G * I)$. For each it prints:

• the correlation of the output with the true $\log r$ over the whole image;
• $\beta$, the regression coefficient of the output on $\log L$, the share of the illumination’s log variation left in a log-domain output (1 for $\log I$; not defined for the linear outputs);
• the detail ratio, the standard deviation of the noise-free output in the bright tile divided by that in the dark tile, 1 when detail no longer depends on the illumination;
• for the noisy image, the noise ratio, the standard deviation of what the noise adds to the output in the dark tile divided by that in the bright tile;
• the detail-to-noise ratio (DNR) in each tile, the standard deviation of the noise-free output divided by that of the added noise, and the dark tile’s DNR relative to that of the high-pass on $I$.

The plan recorded two checks and five hypotheses before the first run, with values computed by a planning script on the same set-up: (HB1) $\beta$ of the homomorphic filter equals 0.0742 to four decimals, the retinex’s $\beta$ is within 0.01 of it, and $\beta$ of $\log I$ is 1; (HB2) the detail ratios of the homomorphic filter and the retinex lie between 0.95 and 1.05, while those of the high-pass on $I$ and of $I$ itself exceed 10; (HB3) the homomorphic filter and the retinex correlate with $\log r$ above 0.95, the high-pass on $I$ below 0.8, and $I$ and $\log I$ below 0.5; (HB4) the noise ratio of the high-pass on $I$ lies between 0.93 and 1.06, and those of the homomorphic filter and the retinex between 30.0 and 36.8; (HB5) the dark-tile DNR of the homomorphic filter and the retinex lies between 0.75 and 1.05 times that of the high-pass on $I$. The reasoning behind HB4 and HB5: to first order the logarithm turns additive noise $n$ into $n/I$, so noise in the log domain is about as many times larger in the dark tile as the illumination is lower, while within a tile, where the illumination is nearly constant, the logarithm scales detail and noise by about the same factor, so to first order it should not raise the dark tile’s ratio of detail to noise above that of the same high-pass applied to $I$. The script uses NumPy 2.3.5, runs in well under a second, and printed identical output on two runs.

import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"          # one BLAS thread, before NumPy
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np

N, s_g, sigma_n = 256, 16.0, 0.0005               # size, Gaussian sigma (px), sensor noise
k = np.fft.fftfreq(N)
F = np.hypot(k[:, None], k[None, :])
G = np.exp(-2 * np.pi ** 2 * s_g ** 2 * F ** 2)   # Gaussian low-pass transfer function
apply = lambda H, img: np.real(np.fft.ifft2(np.fft.fft2(img) * H))

# Illumination: smooth, from 0.02 to 1 (a factor of 50), brightest at (32, 32), darkest at (160, 160).
x = np.arange(N)
mu, rho = np.log(0.02) / 2, -np.log(0.02) / 2
logL = mu + rho * (np.cos(2 * np.pi * (x[None, :] - 32) / N)
                   + np.cos(2 * np.pi * (x[:, None] - 32) / N)) / 2
# Reflectance: one 64 x 64 tile of 8 x 8 random blocks in [0.2, 0.9], repeated 4 x 4.
rng = np.random.default_rng(0)
tile = rng.uniform(0.2, 0.9, (8, 8)).repeat(8, 0).repeat(8, 1)
r = np.tile(tile, (4, 4))
I = np.exp(logL) * r
noisy = I + sigma_n * rng.standard_normal((N, N))

def methods(img):
    li = np.log(img)
    return {"I (no filter)": img, "log I (no filter)": li,
            "high-pass on I": img - apply(G, img),
            "homomorphic": li - apply(G, li),               # high-pass on log I
            "SSR log(I / G*I)": li - np.log(apply(G, img))}

bright, dark = np.s_[0:64, 0:64], np.s_[128:192, 128:192]
corr = lambda a, b: np.corrcoef(a.ravel(), b.ravel())[0, 1]
beta = lambda o: np.cov(o.ravel(), logL.ravel())[0, 1] / logL.var(ddof=1)
fac = 1 - np.exp(-2 * np.pi ** 2 * s_g ** 2 / N ** 2)   # high-pass gain at 1/N cycles/px

print(f"data: L {np.exp(logL).min():.4f}-{np.exp(logL).max():.4f}, r {r.min():.3f}-{r.max():.3f}, "
      f"min I / sigma_n = {I.min() / sigma_n:.2f}, noisy pixels <= 0: {(noisy <= 0).sum()}")
hpL, hpr = logL - apply(G, logL), np.log(r) - apply(G, np.log(r))
clean, noisy_out = methods(I), methods(noisy)
print(f"check: max |HP(log L) - {fac:.6f} (log L - mu)| = {np.abs(hpL - fac * (logL - mu)).max():.1e}, "
      f"max |homomorphic - HP(log L) - HP(log r)| = {np.abs(clean['homomorphic'] - hpL - hpr).max():.1e}"
      f" (both < 1e-12)")

print(f"\n{'output':18s} {'corr log r':>10s} {'beta log L':>10s} {'detail b/d':>10s}"
      f" {'noise d/b':>9s} {'DNR dark':>8s} {'DNR bright':>10s} {'DNR dark / HP on I':>18s}")
part = {m: noisy_out[m] - o for m, o in clean.items()}   # what the noise adds to each output
dnr = {m: o[dark].std() / part[m][dark].std() for m, o in clean.items()}
for m, o in clean.items():
    d = part[m]
    b = f"{beta(o):+10.4f}" if "log" in m or m == "homomorphic" else f"{'n/a':>10s}"
    print(f"{m:18s} {corr(o, np.log(r)):10.4f} {b} {o[bright].std() / o[dark].std():10.4f}"
          f" {d[dark].std() / d[bright].std():9.3f} {dnr[m]:8.3f} {o[bright].std() / d[bright].std():10.2f}"
          f" {dnr[m] / dnr['high-pass on I']:18.3f}")

data: L 0.0200-1.0000, r 0.202-0.898, min I / sigma_n = 8.79, noisy pixels <= 0: 0
check: max |HP(log L) - 0.074209 (log L - mu)| = 1.6e-15, max |homomorphic - HP(log L) - HP(log r)| = 2.3e-15 (both < 1e-12)

output             corr log r beta log L detail b/d noise d/b DNR dark DNR bright DNR dark / HP on I
I (no filter)          0.3301        n/a    36.2095     1.017   10.274     378.19              1.041
log I (no filter)      0.4040    +1.0000     1.0404    35.973    8.273     309.62              0.838
high-pass on I         0.6777        n/a    34.6397     1.016    9.874     347.60              1.000
homomorphic            0.9824    +0.0742     1.0024    35.961    7.922     285.56              0.802
SSR log(I / G*I)       0.9808    +0.0765     0.9972    35.966    8.013     287.38              0.812

Both checks held: the data have the planned ranges and no noisy pixel is zero or negative, and the high-passed illumination equals $0.074209\,(\log L - \mu)$ and the homomorphic output the sum of its two parts, both to $10^{-15}$.

What each output keeps (HB1, HB2 and HB3 held). The unfiltered image carries the reflectance detail in proportion to the illumination: its detail ratio is 36, and it correlates with $\log r$ at only 0.33. Taking the logarithm equalizes the detail (ratio 1.04) but leaves the illumination entirely in place ($\beta = 1$), so the correlation rises only to 0.40. The high-pass on $I$ removes the smooth component of $I$ but not the multiplication by $L$: the remaining detail is still 35 times stronger in the bright tile than in the dark one, and the correlation is 0.68. The homomorphic filter does both. It leaves 0.0742 of the illumination’s log variation, exactly the gain derived above, keeps detail independent of the illumination (ratio 1.0024), and correlates with $\log r$ at 0.98. The retinex behaves almost identically, with $\beta = 0.0765$ and a detail ratio of 0.9972; in our reading, the illumination is smooth enough over the Gaussian’s width for Petro et al.’s approximation to hold, so whether the logarithm comes before or after the surround matters little here. In our reading, the remaining gap to a correlation of 1 is reflectance removed with the illumination: a high-pass filter cannot tell a slow change of reflectance from a slow change of light, and the choice of $\sigma$ sets the balance, since a larger $\sigma$ keeps more of both.

What happens to noise (HB4 and HB5 held). The noise ratio of the high-pass on $I$ is 1.016: a linear filter treats the stationary noise the same in both tiles. In every log-domain output the noise in the dark tile is about 36 times that in the bright tile (35.961 for the homomorphic filter and 35.966 for the retinex, inside the planned range). That is about the factor by which the bright tile is better lit, which the detail ratio of $I$ (36) reflects since both tiles hold the same reflectance, as the first-order argument predicts. The homomorphic filter equalizes detail and noise alike: in the bright tile the detail-to-noise ratio is about 286, in the dark tile 7.9, and the dark tile’s DNR is 0.80 of what the high-pass on $I$ achieves (0.81 for the retinex; the logarithm alone gives 0.84). The logarithm makes dark regions as visible as bright ones without making them any cleaner, and here slightly less clean, so in a real photograph the shadows come out noisy. Denoising before the logarithm, or limiting the gain where $I$ is small, is one remedy, in our reading.

No setting of either example was changed after its first run.

Variants

• Other smooth transitions than the Gaussian and Butterworth ones, for example the steerable pyramid’s low-pass filter, which Briand and Vacher (2016) define as 1 below $r = \pi/4$, 0 above $\pi/2$ and $\cos\big(\tfrac{\pi}{2}\log_2(4r/\pi)\big)$ in between, with $r$ in radians per sample.
• Band-pass and band-reject filters keep or remove an annulus of frequencies; oriented band-pass filters, such as Gabor filters, add an angular selection.
• Homomorphic filtering applies any of these to the logarithm of the image. Single- and multiscale retinex methods are relatives that subtract the logarithm of a surround average instead; multiscale retinex takes a weighted sum of single-scale outputs at several Gaussian scales, which, as Petro et al. (2014) report, Jobson, Rahman and Woodell’s second 1997 paper, on the multiscale retinex, fixed experimentally at 15, 80 and 250.
• Filters estimated from data. The Wiener filter computes the gain at each frequency from the power spectra of the signal and the noise, and image deblurring regularizes the inverse of a blur’s transfer function. Both are frequency domain filters whose transfer function is derived rather than chosen.

Best Practices

• State the convention with the cutoff. Give the cutoff in cycles per pixel (or as a fraction of the sampling frequency, as scikit-image does) and the gain at the cutoff. In scikit-image 0.25.2, butterworth defaults to the squared form (squared_butterworth=True, gain $\tfrac12$ at the cutoff), to a high-pass filter (high_pass=True), and to cutoff_frequency_ratio=0.005 and order=2.0; a call that means a traditional low-pass Butterworth filter must set two of these explicitly. Compare filters at a matched gain, not a matched parameter.
• Pad or taper the borders. The DFT filters the periodic extension of the image. Padding by reflection or replication, then cropping, keeps content from one border out of the other; scikit-image’s butterworth has an npad argument that pads each edge by replicating the border pixels before filtering.
• Prefer smooth transitions where ringing matters; neither of the example’s Gaussians overshoots.
• Keep the transfer function symmetric for real images, $H(-u, -v) = H(u, v)$, so that the output is real; then discard the rounding-level imaginary part, or use real-input transforms such as rfft2, which scikit-image uses for real images.
• For homomorphic filtering, keep $I$ positive and the noise in mind. Add a small offset or clip before the logarithm, and denoise first when the dark regions matter.

Pitfalls

• Circular convolution. Without padding, content wraps around, and borders that do not match produce the seams that the Fourier transform article describes.
• Mistaking one Butterworth convention for the other. Filters from different libraries or textbooks with the same order and cutoff can differ by a factor of two in their attenuation in decibels, and a cutoff set for gain $1/\sqrt2$ is effectively lower in the squared convention.
• Off-centre frequency grids. A transfer function built on an fftshifted grid must be shifted back before it multiplies an unshifted spectrum; a mismatch puts the passband in the corners of the spectrum and produces a high-pass result. For even sizes, the Nyquist frequency is its own mirror image on the DFT grid, and Briand and Vacher (2016) show that how its coefficients are treated is a legitimate choice whose effect is small but not zero.
• Higher order is not sharper. At a fixed gain at the cutoff, raising a Butterworth filter’s order above 2 adds ringing without narrowing edges in the example, and in the plain convention it widens them.
• Separation that is not there. A homomorphic filter can only remove the part of the illumination below its cutoff and keep the part of the reflectance above it. In our reading, a hard shadow has an edge as sharp as any reflectance edge, so its boundary survives the filter, and a large region of uniform reflectance loses its level along with the illumination.

Evidence and Impact

The correctness of the basic algorithm is a theorem, the convolution theorem, and Briand and Vacher (2016) give it a precise meaning for filters defined by continuous functions. They use the first of their three methods in their implementation of the Heeger–Bergen texture synthesis algorithm because it guarantees exact reconstruction from the pyramid.

For illumination removal, Petro, Sbert and Morel (2014) analyse multiscale retinex, a log-domain centre–surround filter they describe as probably the most successful of its kind, and propose an implementation that keeps the original method and replaces its post-processing; they conclude that the method must be applied to each colour channel when colour balance is wanted, or to luminance only for local contrast enhancement. Morel, Petro and Sbert (2014) illustrate their own screened-Poisson method, which belongs to the same group of filtering methods as homomorphic filtering, on fundus photographs and MR images.