Intensity Transforms and Frequency-Domain Filtering

A thermal camera on a furnace line produces images where the hot zone reads 4000 and the cool background reads 3. Displayed normally, the background is a solid black wall and the entire scene looks like one bright blob on nothing. Everything is technically there. None of it is visible.

This post is about the operations that fix that — and about a second way of looking at an image that makes blurring, sharpening and denoising all the same operation with a different setting.

Prerequisites: image fundamentals for histograms and convolution, and math for CV for the idea of a kernel.

Where this came from

Two problems created this field, and both are worth knowing because both still describe why you would use it.

In 1964 the Jet Propulsion Laboratory received the first close-up photographs of the Moon from Ranger 7. The images arrived corrupted by transmission noise and geometric distortion from the vidicon camera. Computers were used to correct them, and that work is generally taken as the start of digital image processing.

In 1979 Godfrey Hounsfield and Allan Cormack shared the Nobel Prize in Physiology or Medicine for computed tomography. A CT scan is not photographed — it is computed from X-ray absorption measurements. Without image processing there is no image at all.

The pattern is the same in both: the sensor gives you numbers, not a picture. Turning those numbers into something a human or an algorithm can use is the job.

Two kinds of operation

flowchart TD
  A["Image operation"] --> B["Point operation"]
  A --> C["Neighbourhood operation"]
  B --> B1["Output at (u,v) depends<br/>only on input at (u,v)"]
  B --> B2["Negative, log, gamma,<br/>thresholding, slicing"]
  C --> C1["Output at (u,v) depends<br/>on a window around (u,v)"]
  C --> C2["Blur, sharpen, median,<br/>edge detection"]
  style B fill:#eef
  style C fill:#efe

A point operation is a lookup table. There are only 256 possible input values in an 8-bit image, so you can precompute all 256 outputs once and apply them with a single array index. That makes point operations essentially free, no matter how complicated the maths looks.

Neighbourhood operations cost real time because each output pixel requires reading a whole window.

Point operations, worked

Let us use a small strip of real intensities to make each transform concrete:

The negative

The simplest one. For an image with L levels:

g = (L - 1) - f

This looks like a novelty until you work with medical images. Radiologists routinely invert X-rays and mammograms, because small bright details on a dark background are harder for the human eye to pick out than small dark details on a light background. Nothing is added — the same information is simply easier to see.

The log transform

Here is the fix for the furnace problem.

s = c · log(1 + |r|)

The 1 + exists so that r = 0 maps to s = 0 rather than to negative infinity. The constant c scales the result back into display range. To find it, force the maximum input to map to the maximum output:

c = (L - 1) / log(1 + R_max)

With L = 256 and a sensor whose peak reading is R_max = 4000:

c = 255 / log(1 + 4000)
  = 255 / log(4001)
  = 255 / 8.294
  = 30.75

Now apply it to a range of readings:

Look at the last column. Linear scaling maps a reading of 47 to 3 — indistinguishable from black. The log transform maps it to 119, comfortably mid-grey.

Dynamic range and contrast are different things, and confusing them causes real mistakes. Dynamic range is the span of values present. Contrast is how much difference there is between neighbouring regions. An image can have huge dynamic range and terrible local contrast — the furnace image is exactly that. The log transform does not add information; it reallocates display range from where the data is sparse to where it is dense.

Gamma correction

s = c · r^γ

Working on normalised values in [0, 1] with c = 1, the exponent decides everything:

Take r = 0.25:

γ = 0.5:  s = 0.25^0.5 = 0.500   → doubled
γ = 1.0:  s = 0.25^1.0 = 0.250   → unchanged
γ = 2.0:  s = 0.25^2.0 = 0.0625  → quartered

The property that makes gamma useful is that it never clips. Both endpoints stay put, so no matter how aggressive the exponent, you never lose highlights to white or shadows to black. Brightness offsets clip; gamma does not.

Gray-level slicing

Sometimes you do not want a smooth curve at all — you want one band of intensities and nothing else.

Suppose a thermal inspection cares only about surfaces between 180 and 220 in the mapped image. Two variants:

Binary slicing — the band becomes white, everything else black:

Preserving slicing — the band is brightened, everything else keeps its original value:

Binary slicing gives you a clean mask for measuring area. Preserving slicing keeps context so a human can see where the highlighted region sits. Choose by whether the output is for a machine or a person.

import numpy as np

def slice_levels(img, low, high, preserve=False):
    mask = (img >= low) & (img <= high)
    out = img.copy() if preserve else np.zeros_like(img)
    out[mask] = 255
    return out

The frequency view

Now the second half, and the more powerful idea.

Any image can be written as a sum of 2D sine waves of different frequencies, orientations and amplitudes. The Fourier transform tells you the amplitude of each. Nothing is lost — you can transform back exactly.

flowchart LR
  A["Image f(u,v)<br/>spatial domain"] -->|"Fourier transform"| B["Spectrum F(u,v)<br/>frequency domain"]
  B -->|"multiply by H(u,v)"| C["Filtered G = F·H"]
  C -->|"inverse Fourier"| D["Output g(u,v)"]
  A -->|"convolve with h"| D
  style B fill:#eef
  style C fill:#eef

Follow the two paths in that diagram. You can go the long way round — transform, multiply, transform back — or you can convolve directly in the spatial domain. They give the same answer. That is the convolution theorem:

g = f * h    ⟺    G = F · H

Convolution in space is multiplication in frequency.

What this actually tells you

This equivalence explains things that otherwise have to be memorised separately.

Take the two kernels you already know:

Sum the mean kernel’s coefficients: 9 × (1/9) = 1. A kernel whose coefficients sum to 1 passes a constant region unchanged — constant means frequency zero, so it passes DC. And it averages away rapid variation, so it kills high frequency. That is a low-pass filter, derived from nothing but the numbers in the box.

Now sum the sharpen kernel: 8/9 − 8×(1/9) = 0. A kernel summing to zero outputs zero on any flat region. It responds only where things change. That is a high-pass filter.

You can read a kernel’s frequency behaviour straight off its coefficient sum. Sum to 1, low-pass. Sum to 0, high-pass.

Filtering in the frequency domain directly

An ideal low-pass filter keeps everything within a radius D₀ of the centre and deletes the rest:

H(u,v) = 1   if D(u,v) ≤ D₀
         0   if D(u,v) >  D₀

where D(u,v) is the distance from the frequency-domain origin.

This is intuitive and it produces a visible artifact. A perfectly sharp cutoff in frequency corresponds to a spatial kernel that oscillates — it rings. The result has faint ripples radiating from every strong edge, like ripples in water.

The Butterworth filter softens the transition:

H(u,v) = 1 / (1 + [D(u,v)/D₀]^(2n))

The order n controls sharpness. Work through the values at the cutoff and around it, with D₀ = 30:

Higher order means a sharper transition, better frequency separation, and more ringing. Lower order means gentler separation and clean edges. There is no setting that gives both, and that is a property of the mathematics rather than a limitation of the implementation.

The honest summary: in day-to-day work you will filter spatially. The frequency view earns its place by explaining why the spatial filters behave as they do, and by handling periodic noise, which nothing spatial handles well.

Doing it in code

import cv2
import numpy as np

# --- Point operations as lookup tables (256 entries, applied instantly) ---

def log_lut(r_max=255.0):
    c = 255.0 / np.log(1.0 + r_max)
    table = c * np.log(1.0 + np.arange(256, dtype=np.float64))
    return np.clip(table, 0, 255).astype(np.uint8)

def gamma_lut(gamma):
    table = ((np.arange(256) / 255.0) ** gamma) * 255.0
    return np.clip(table, 0, 255).astype(np.uint8)

gray = cv2.imread("input.png", cv2.IMREAD_GRAYSCALE)
logged    = cv2.LUT(gray, log_lut())
brightened = cv2.LUT(gray, gamma_lut(0.5))
negative   = 255 - gray            # no LUT needed

cv2.LUT is the right tool here. It replaces a per-pixel pow call with a single array index, and on a 4K image that is the difference between milliseconds and hundreds of milliseconds.

For real high-dynamic-range data, do not force it into uint8 first:

raw = load_thermal()                        # float32, range 3 .. 4000
c = 255.0 / np.log(1.0 + raw.max())
disp = np.clip(c * np.log(1.0 + raw), 0, 255).astype(np.uint8)

Frequency-domain filtering:

def butterworth_lowpass(shape, d0, n=2):
    h, w = shape
    v = np.arange(h)[:, None] - h / 2
    u = np.arange(w)[None, :] - w / 2
    d = np.sqrt(u**2 + v**2)
    return 1.0 / (1.0 + (d / d0) ** (2 * n))

f  = np.fft.fftshift(np.fft.fft2(gray.astype(np.float32)))
H  = butterworth_lowpass(gray.shape, d0=30, n=2)
out = np.real(np.fft.ifft2(np.fft.ifftshift(f * H)))
out = np.clip(out, 0, 255).astype(np.uint8)

The two fftshift calls are the step people get wrong. fft2 puts the zero frequency at the array corner; fftshift moves it to the centre so a centred radial mask works. You must undo it with ifftshift before the inverse transform, or the output comes back scrambled.

To see a spectrum:

spectrum = 20 * np.log(np.abs(f) + 1)       # log again — spectra have huge dynamic range
spectrum = cv2.normalize(spectrum, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)

Note that the log transform reappears here for exactly the reason from the top of this post: the DC term dwarfs everything else, and without a log you see one bright dot and nothing more.

Practice task

About forty-five minutes.

1. Find or make a high-dynamic-range image — a photo of a lit window from inside a dark room works well. Display it with plain linear scaling, then with the log transform. Compute c yourself from the actual maximum.
2. Apply gamma at 0.4, 0.7, 1.5 and 2.5. Plot the histogram after each. Notice that nothing ever piles up at 0 or 255.
3. Compare a +60 brightness shift against γ = 0.6. Both brighten. Look at the top of the histogram in each case and identify which one clipped.
4. Take a photo and add periodic noise: img + 30*np.sin(np.arange(w)[None,:] * 0.4). Try to remove it with a Gaussian blur. Then look at the FFT magnitude, find the two bright points off-centre, zero them, and invert.
5. Run an ideal low-pass at D₀ = 20, then Butterworth n = 2 at the same cutoff. Zoom into a strong edge in each and find the ringing.

Step 4 is the one that makes the frequency domain click. The noise that no spatial filter could touch becomes two pixels you delete by hand.

Summary

Point operations change each pixel using only its own value, which makes them lookup tables and therefore free. The negative aids human reading, the log transform compresses dynamic range so faint detail survives the trip to an 8-bit display, gamma adjusts perceived brightness without ever clipping, and gray-level slicing isolates a known band when physics has already given you the threshold.

The frequency view says any image is a sum of sine waves, with smooth regions at low frequency and edges and noise at high frequency. The convolution theorem — g = f * h ⟺ G = F · H — means every spatial filter you already use is a frequency filter in disguise, and you can read which kind off the kernel’s coefficient sum.

Filtering directly in frequency is worth it for periodic interference and for very large kernels. Its cost is ringing whenever the cutoff is sharp, which the Butterworth family trades against frequency selectivity. You cannot have both.

What comes next

Edge detection and thresholding takes the high-pass idea and turns it into a working edge pipeline, including why Canny blurs before it differentiates — which is now a statement about frequency rather than an arbitrary step.

If you want to see the frequency domain used for something other than enhancement, how images and video are stored shows JPEG using the closely related DCT to compress rather than filter, keeping the low frequencies and rounding the high ones to zero.