Practical: Fourier Optics and Image Processing

12. Practical: Fourier Optics and Image Processing#

12.1. Introduction#

Digital images are two-dimensional intensity functions $f(x,y)$. The information in an image is encoded in a pattern of variation — bright and dark regions, sharp edges, smooth gradients. Fourier analysis gives us a way to decompose any such pattern into a sum of simple sinusoidal waves, each characterised by a frequency, an amplitude, and a phase. This decomposition is not merely mathematical convenience: in electron microscopy it is the physical basis of image formation (through the contrast transfer function) and of reconstruction algorithms (through the Fourier slice theorem).

This practical builds up Fourier intuition from the ground up:

What is a wave?
Fourier series — decomposing arbitrary 1D signals
Frequency spectra — looking at signals through the lens of their frequency content
2D Fourier analysis — extending to images
The convolution theorem — efficient filtering
The contrast transfer function — how the microscope shapes the signal

Interactive elements

Click Live Code in the top toolbar to activate the kernel, then expand each Show code toggle and click ▶ to run the cell.

12.1.1. Setup#

12.2. Part 1 – What is a wave?#

A sinusoidal wave has three fundamental parameters:

\[ f(t) = A \cdot \cos(2\pi f \cdot t - \varphi) \]

Parameter	Symbol	Effect
Amplitude	$A$	Height of the wave
Frequency	$f$	Number of cycles per unit time
Phase	$\varphi$	Horizontal shift

12.3. Part 2 – Fourier series#

Any periodic signal can be expressed as a sum of sinusoidal waves (Fourier series):

\[ b(x) = A_0 + \sum_{k=1}^{F} A_k \cdot \cos\!\left(\frac{2\pi k x}{P} - \varphi_k\right) \]

where $F$ is the maximum frequency, $P$ the period, and $A_k$, $\varphi_k$ the amplitude and phase at frequency $k$.

The amplitudes and phases are computed from the signal via the Fourier coefficients:

\[ a_k = \frac{2}{P}\sum_{i=0}^{P} b(i)\cos\!\frac{2\pi k i}{P}, \qquad b_k = \frac{2}{P}\sum_{i=0}^{P} b(i)\sin\!\frac{2\pi k i}{P}, \qquad A_k = \sqrt{a_k^2 + b_k^2} \]

12.3.1. The box function#

A classic example is the box (rectangular) function:

\[\begin{split} b(x) = \begin{cases} 1 & -a/2 < x < a/2 \\ 0 & \text{elsewhere} \end{cases} \end{split}\]

Its Fourier series contains only odd harmonics, and the coefficients decay as $1/k$. As we add more terms, the reconstruction improves but Gibbs ringing remains at the edges.

Show code cell source

Hide code cell source

def getbox(width=200, N=1000):
    """Box function centred at N//2 with given pixel width."""
    sig = np.zeros(N)
    half = width // 2
    sig[N//2 - half : N//2 + half + 1] = 1
    return sig

def fourier_series_reconstruct(signal, tot_freq):
    """Reconstruct signal from its first tot_freq Fourier components."""
    N    = len(signal)
    x    = np.arange(N)
    recon = np.zeros(N)
    for k in range(tot_freq):
        ak = 2/N * np.sum(signal * np.cos(2*np.pi*k*x/N))
        bk = 2/N * np.sum(signal * np.sin(2*np.pi*k*x/N))
        Ak = np.sqrt(ak**2 + bk**2)
        pk = np.arctan2(bk, ak)
        recon += Ak * np.cos(2*np.pi*k*x/N - pk)
    recon -= recon.mean() - signal.mean()
    return recon

width_sl  = IntSlider(value=200, min=50, max=500, step=50,
                       description="Box width (px)", style=_sty, layout=_sl)
nwave_sl  = IntSlider(value=15,  min=1,  max=80,  step=1,
                       description="# waves F",      style=_sty, layout=_sl)
out_box = Output()

def _update_box(_=None):
    sig   = getbox(width_sl.value)
    recon = fourier_series_reconstruct(sig, nwave_sl.value)
    xax   = np.linspace(-500, 500, len(sig))
    fig, axes = plt.subplots(1, 2, figsize=(13, 3.8))
    axes[0].plot(xax, sig,   color='gray',     lw=1.2, label='Box function')
    axes[0].plot(xax, recon, color='steelblue', lw=1.6, label=f'Reconstruction (F={nwave_sl.value})')
    axes[0].set_xlabel("x (pixels)"); axes[0].set_ylabel("Amplitude")
    axes[0].set_title(f"Fourier series reconstruction: {nwave_sl.value} waves")
    axes[0].set_ylim(-0.4, 1.5); axes[0].legend(fontsize=9); axes[0].grid(alpha=0.3)

    # Amplitude spectrum
    N    = len(sig)
    ks   = np.arange(1, nwave_sl.value+1)
    x    = np.arange(N)
    amps = []
    for k in ks:
        ak = 2/N * np.sum(sig * np.cos(2*np.pi*k*x/N))
        bk = 2/N * np.sum(sig * np.sin(2*np.pi*k*x/N))
        amps.append(np.sqrt(ak**2 + bk**2))
    axes[1].bar(ks, amps, color='steelblue', edgecolor='k', linewidth=0.4)
    axes[1].set_xlabel("Frequency k"); axes[1].set_ylabel("Amplitude $A_k$")
    axes[1].set_title("Amplitude spectrum of the box function")
    axes[1].grid(alpha=0.3, axis='y')
    plt.tight_layout()
    with out_box:
        clear_output(wait=True); display(fig2img(fig))

for w in [width_sl, nwave_sl]:
    w.observe(_update_box, names='value')
display(VBox([width_sl, nwave_sl, out_box]))
_update_box()

12.4. Part 3 – Frequency spectrum (1D FFT)#

Instead of computing Fourier coefficients term by term, the discrete Fourier transform (DFT) computes all of them at once:

\[ F(k) = \frac{1}{P} \sum_{m=0}^{P-1} b(m) \cdot e^{-i 2\pi k m / P} \]

The fast Fourier transform (FFT) computes this in $O(N \log N)$ rather than $O(N^2)$. The frequency spectrum (plot of $|F(k)|$ vs $k$) reveals which frequencies are present in the signal.

12.4.1. Nyquist sampling theorem#

If a signal contains no frequencies higher than $W$ Hz, it is fully determined by samples taken every $1/(2W)$ seconds:

\[ f_\text{max} \leq \frac{1}{2 \Delta t} \quad \Leftrightarrow \quad \Delta t \leq \frac{1}{2 f_\text{max}} \]

The Nyquist frequency $f_N = 1/(2\Delta t)$ is the highest frequency measurable at a given sampling rate. In EM images, the pixel size $\Delta x$ sets $f_N = 1/(2\Delta x)$, and only features larger than $2\Delta x$ can be resolved.

Show code cell source

Hide code cell source

a1_sl = FloatSlider(value=1.0, min=0.0, max=5.0, step=0.5, description="Amp₁", style=_sty, layout=_sl)
f1_sl = IntSlider(  value=10,  min=1,   max=100, step=1,   description="Freq₁", style=_sty, layout=_sl)
a2_sl = FloatSlider(value=1.0, min=0.0, max=5.0, step=0.5, description="Amp₂", style=_sty, layout=_sl)
f2_sl = IntSlider(  value=20,  min=1,   max=100, step=1,   description="Freq₂", style=_sty, layout=_sl)
a3_sl = FloatSlider(value=1.0, min=0.0, max=5.0, step=0.5, description="Amp₃", style=_sty, layout=_sl)
f3_sl = IntSlider(  value=40,  min=1,   max=100, step=1,   description="Freq₃", style=_sty, layout=_sl)
dc_sl = FloatSlider(value=5.0, min=-5.0, max=10.0, step=0.5, description="DC offset", style=_sty, layout=_sl)
out_fft1 = Output()

def _update_fft1(_=None):
    tot_time = 1.0; N = 1000
    t = np.linspace(0, tot_time, N)
    sig = (dc_sl.value
           + a1_sl.value * np.cos(2*np.pi*f1_sl.value*t)
           + a2_sl.value * np.cos(2*np.pi*f2_sl.value*t)
           + a3_sl.value * np.cos(2*np.pi*f3_sl.value*t))
    dft   = np.fft.fft(sig)
    freqs = np.fft.fftfreq(N, tot_time/N)
    amp   = np.abs(dft) / N

    fig, axes = plt.subplots(1, 2, figsize=(13, 3.8))
    axes[0].plot(t, sig, color='steelblue', lw=1.4)
    axes[0].set_xlabel("Time (s)"); axes[0].set_ylabel("Amplitude")
    axes[0].set_title("Signal"); axes[0].grid(alpha=0.3)
    flim = min(100, max(f1_sl.value, f2_sl.value, f3_sl.value) + 20)
    axes[1].bar(freqs, amp, width=freqs[1]-freqs[0] if len(freqs) > 1 else 1,
                color='steelblue', edgecolor='none')
    axes[1].set_xlabel("Frequency (Hz)"); axes[1].set_ylabel("Amplitude")
    axes[1].set_xlim(-flim, flim)
    axes[1].set_title("Frequency spectrum (FFT)")
    axes[1].grid(alpha=0.3)
    plt.tight_layout()
    with out_fft1:
        clear_output(wait=True); display(fig2img(fig))

for w in [a1_sl, f1_sl, a2_sl, f2_sl, a3_sl, f3_sl, dc_sl]:
    w.observe(_update_fft1, names='value')
display(VBox([HBox([a1_sl, f1_sl]), HBox([a2_sl, f2_sl]),
              HBox([a3_sl, f3_sl]), dc_sl, out_fft1]))
_update_fft1()

Question 1

Given a detector that records 15 frames per 100 s. What is the maximum frequency component that can be recovered? What does this imply for the pixel size required to resolve a feature of width 2 Å in a cryo-EM image?

12.5. Part 4 – 2D Fourier analysis#

Images are 2D signals. A 2D sinusoidal wave has the form:

\[ f(x,y) = A \cdot \sin(2\pi f_x x + 2\pi f_y y + \varphi) \]

with separate spatial frequencies $f_x$ (cycles per pixel along $x$) and $f_y$ (cycles per pixel along $y$). The 2D Fourier transform of an image gives the amplitude and phase at every spatial frequency $(f_x, f_y)$.

12.5.1. Properties#

Translation property: Shifting an image in real space does not change the amplitude spectrum — only the phase spectrum.

Rotation property: Rotating an image rotates its Fourier transform by the same angle.

The interactive below shows 2D sine waves at different spatial frequencies, and the effects of translation and rotation on the Fourier transform of a real image.

12.5.2. 2D FT of real images#

Different image features produce characteristic Fourier signatures. A ring in real space produces a ring in Fourier space; a line produces a perpendicular line through the origin.

12.5.3. Translation and rotation properties#

Question 2

Play with the translation and rotation of the image and observe the amplitude spectrum. What changes and what stays the same in each case? Can you explain why?

12.6. Part 5 – The convolution theorem#

Convolution describes how one function modifies another by sweeping over all positions:

\[ (f * g)(y) = \int_{-\infty}^{\infty} f(x)\,g(y-x)\,dx \]

In image processing, convolving an image with a kernel (a small matrix) applies a local operation: smoothing, sharpening, or edge detection. The convolution theorem states:

\[ \mathcal{F}\{f * g\} = \mathcal{F}\{f\} \cdot \mathcal{F}\{g\} \]

Convolution in real space equals multiplication in Fourier space. This makes filtering in Fourier space much faster for large kernels than direct spatial convolution.

12.6.1. Common kernels#

Show code cell source

Hide code cell source

import io, numpy as np, matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from scipy.signal import convolve2d
from IPython.display import display, Image

def sobel_x(n=3):
    return np.array([[-1,0,1],[-2,0,2],[-1,0,1]], float)

def sobel_y(n=3):
    return np.array([[1,2,1],[0,0,0],[-1,-2,-1]], float)

def laplacian():
    return np.array([[0,1,0],[1,-4,1],[0,1,0]], float)

def gaussian_kernel(size=9, sigma=2.0):
    ax = np.arange(-(size//2), size//2+1)
    k  = np.exp(-ax**2/(2*sigma**2))
    k  = np.outer(k, k); return k / k.sum()

def box_blur(size=5):
    return np.ones((size,size)) / size**2

def sharpen():
    return np.array([[0,-1,0],[-1,5,-1],[0,-1,0]], float)

Q64 = make_letter('Q', 64)

kernels = {
    'Gaussian blur (σ=2)': gaussian_kernel(9, 2.0),
    'Box blur (5×5)':       box_blur(5),
    'Laplacian':            laplacian(),
    'Sobel X':              sobel_x(),
    'Sobel Y':              sobel_y(),
    'Sharpen':              sharpen(),
}

fig, axes = plt.subplots(3, len(kernels), figsize=(3.5*len(kernels), 9))
for col, (name, kernel) in enumerate(kernels.items()):
    filtered = convolve2d(Q64, kernel, mode='same', boundary='wrap')
    kshow    = kernel.copy()
    axes[0, col].imshow(Q64,      cmap='gray',    origin='lower'); axes[0, col].axis('off')
    axes[0, col].set_title(name,  fontsize=7)
    n_k      = kshow.shape[0]
    axes[1, col].imshow(kshow,    cmap='RdBu',    origin='lower',
                        vmin=-abs(kshow).max(), vmax=abs(kshow).max())
    axes[1, col].set_title(f"Kernel ({n_k}×{n_k})", fontsize=7); axes[1, col].axis('off')
    axes[2, col].imshow(filtered, cmap='gray',    origin='lower'); axes[2, col].axis('off')
    axes[2, col].set_title("Result", fontsize=7)

axes[0, 0].set_ylabel("Input", fontsize=9)
axes[1, 0].set_ylabel("Kernel", fontsize=9)
axes[2, 0].set_ylabel("Filtered output", fontsize=9)
plt.tight_layout()
_buf = io.BytesIO(); fig.savefig(_buf, format='png', bbox_inches='tight', dpi=96); _buf.seek(0)
display(Image(_buf.read()))
plt.close('all')

12.6.2. Interactive kernel#

Show code cell source

Hide code cell source

ker_dd  = Dropdown(options=['Gaussian blur','Box blur','Laplacian','Sobel X','Sobel Y','Sharpen'],
                   value='Gaussian blur', description="Kernel:", style=_sty)
ksz_sl  = IntSlider(value=9, min=3, max=21, step=2,
                     description="Kernel size", style=_sty, layout=_sl)
out_ker = Output()

def _update_ker(_=None):
    name = ker_dd.value; ksz = ksz_sl.value
    if   name == 'Gaussian blur': kernel = gaussian_kernel(ksz, ksz/6)
    elif name == 'Box blur':      kernel = box_blur(ksz)
    elif name == 'Laplacian':     kernel = laplacian()
    elif name == 'Sobel X':       kernel = sobel_x()
    elif name == 'Sobel Y':       kernel = sobel_y()
    else:                         kernel = sharpen()

    filtered  = convolve2d(Q64, kernel, mode='same', boundary='wrap')
    ps_input  = power_spectrum_2d(Q64)
    ps_kernel = power_spectrum_2d(np.pad(kernel, (32-kernel.shape[0]//2,
                                                  32-kernel.shape[0]//2+1)))
    ps_out    = power_spectrum_2d(filtered)

    fig, axes = plt.subplots(2, 3, figsize=(13, 8))
    axes[0,0].imshow(Q64,       cmap='gray',    origin='lower'); axes[0,0].axis('off')
    axes[0,0].set_title("Input image")
    axes[0,1].imshow(kernel, cmap='RdBu', origin='lower',
                     vmin=-abs(kernel).max(), vmax=abs(kernel).max())
    axes[0,1].set_title(f"Kernel: {name}\n({kernel.shape[0]}×{kernel.shape[1]})")
    axes[0,1].axis('off')
    axes[0,2].imshow(filtered,  cmap='gray',    origin='lower'); axes[0,2].axis('off')
    axes[0,2].set_title("Filtered output (real space)")
    axes[1,0].imshow(ps_input,  cmap='inferno', origin='lower'); axes[1,0].axis('off')
    axes[1,0].set_title("FT of input")
    axes[1,1].imshow(ps_kernel, cmap='inferno', origin='lower'); axes[1,1].axis('off')
    axes[1,1].set_title("FT of kernel (transfer function)")
    axes[1,2].imshow(ps_out,    cmap='inferno', origin='lower'); axes[1,2].axis('off')
    axes[1,2].set_title("FT of output = FT(in) × FT(kernel)")
    plt.suptitle("Convolution theorem: multiplication in Fourier space", fontsize=10)
    plt.tight_layout()
    with out_ker:
        clear_output(wait=True); display(fig2img(fig))

for w in [ker_dd, ksz_sl]:
    w.observe(_update_ker, names='value')
display(VBox([ker_dd, ksz_sl, out_ker]))
_update_ker()

Question 3

What does the FT of a Gaussian blur kernel look like? What does this tell you about which spatial frequencies are attenuated? How does this compare to the Laplacian kernel?

12.7. Part 6 – Low-pass and high-pass filtering in Fourier space#

We can design filters by directly modifying the Fourier transform: set Fourier coefficients outside (or inside) a radius to zero before inverse-transforming. This is equivalent to convolving with the Fourier transform of the mask function.

Show code cell source

Hide code cell source

ftype_dd = Dropdown(options=['Low-pass', 'High-pass', 'Band-pass'],
                    value='Low-pass', description="Filter type:", style=_sty)
fcut_sl  = FloatSlider(value=0.2, min=0.02, max=0.5, step=0.02,
                        description="Cutoff f₁",    style=_sty, layout=_sl)
fcut2_sl = FloatSlider(value=0.4, min=0.02, max=0.5, step=0.02,
                        description="Cutoff f₂ (band)", style=_sty, layout=_sl)
out_lp   = Output()

def freq_mask(n, f1, f2=None, ftype='Low-pass'):
    """Build a circular frequency mask."""
    y, x  = np.mgrid[0:n, 0:n].astype(float)
    y -= n//2; x -= n//2
    R = np.sqrt(x**2 + y**2) / n
    if ftype == 'Low-pass':
        return (R <= f1).astype(float)
    elif ftype == 'High-pass':
        return (R >= f1).astype(float)
    else:
        return ((R >= f1) & (R <= (f2 or f1))).astype(float)

def _update_lp(_=None):
    n      = Q64.shape[0]
    ftype  = ftype_dd.value
    f1     = fcut_sl.value
    f2     = fcut2_sl.value
    mask   = freq_mask(n, f1, f2, ftype)
    F      = np.fft.fftshift(np.fft.fft2(Q64))
    F_filt = F * mask
    filtered = np.real(np.fft.ifft2(np.fft.ifftshift(F_filt)))

    fig, axes = plt.subplots(1, 4, figsize=(16, 4.2))
    axes[0].imshow(Q64,      cmap='gray',    origin='lower'); axes[0].axis('off')
    axes[0].set_title("Original")
    axes[1].imshow(mask,     cmap='Blues',   origin='lower'); axes[1].axis('off')
    axes[1].set_title(f"Fourier mask\n({ftype}, f₁={f1:.2f})")
    axes[2].imshow(np.log(np.abs(F_filt)+1), cmap='inferno', origin='lower')
    axes[2].axis('off'); axes[2].set_title("Filtered FT")
    axes[3].imshow(filtered, cmap='gray',    origin='lower'); axes[3].axis('off')
    axes[3].set_title("Filtered image")
    plt.tight_layout()
    with out_lp:
        clear_output(wait=True); display(fig2img(fig))

for w in [ftype_dd, fcut_sl, fcut2_sl]:
    w.observe(_update_lp, names='value')
display(VBox([ftype_dd, fcut_sl, fcut2_sl, out_lp]))
_update_lp()

12.8. Part 7 – The contrast transfer function (CTF)#

In phase-contrast electron microscopy, the image intensity is related to the projected potential of the specimen through the contrast transfer function (CTF). The CTF modulates the Fourier amplitudes as a function of spatial frequency $k$:

\[ \text{CTF}(k) = -\sin\!\left[\Delta\varphi + \frac{-\pi}{2}C_s\lambda^3 k^4 + \pi\lambda\Delta_f k^2\right] \]

where:

Symbol	Quantity
$k$	spatial frequency (Å⁻¹)
$C_s$	spherical aberration coefficient
$\lambda$	relativistic electron wavelength
$\Delta_f$	defocus (positive = underfocus)
$\Delta\varphi$	additional phase shift (e.g. from phase plate)

The electron wavelength depends on the accelerating voltage $V$ (in eV):

\[ \lambda = \frac{12.264}{\sqrt{V(1 + V \cdot 0.98 \times 10^{-6})}} \; \text{Å} \]

Defocusing causes phase reversals at specific spatial frequencies (the CTF zeros). Frequencies near these zeros are not faithfully transferred to the image.

Show code cell source

Hide code cell source

def relativistic_lambda(voltage_eV):
    """Relativistic electron wavelength in Angstrom."""
    return 12.264 / np.sqrt(voltage_eV * (1 + voltage_eV * 0.98e-6))

def ctf_1d(defocus_um, voltage_kV=300, Cs_mm=2.7, delta_phi=0, B_factor=0,
           n_pts=256, max_freq=0.5):
    """Compute 1D CTF curve."""
    k       = np.linspace(0, max_freq, n_pts)
    V_eV    = voltage_kV * 1e3
    lam     = relativistic_lambda(V_eV)
    df_A    = defocus_um * 1e4        # µm → Å
    Cs_A    = Cs_mm * 1e7             # mm → Å
    gamma   = (-np.pi/2)*Cs_A*lam**3*k**4 + np.pi*lam*df_A*k**2
    ctf     = -np.sin(delta_phi + gamma)
    if B_factor > 0:
        ctf *= np.exp(-B_factor * k**2)
    return k, ctf

df_sl   = FloatSlider(value=2.0,  min=0.5,  max=10.0, step=0.5,
                       description="Defocus (µm)",   style=_sty, layout=_sl)
volt_sl = FloatSlider(value=300,  min=100,  max=300,  step=100,
                       description="Voltage (kV)",   style=_sty, layout=_sl)
cs_sl   = FloatSlider(value=2.7,  min=0.0,  max=5.0,  step=0.1,
                       description="Cs (mm)",         style=_sty, layout=_sl)
phi_sl  = FloatSlider(value=0.0,  min=-np.pi/2, max=np.pi/2, step=0.1,
                       description="Phase shift Δφ", style=_sty, layout=_sl)
bf_sl   = FloatSlider(value=0.0,  min=0.0,  max=200.0, step=10.0,
                       description="B-factor (Å²)",  style=_sty, layout=_sl)
out_ctf = Output()

def _update_ctf(_=None):
    k, ctf_curve = ctf_1d(df_sl.value, volt_sl.value, cs_sl.value,
                           phi_sl.value, bf_sl.value)
    lam = relativistic_lambda(volt_sl.value * 1e3)

    # Apply CTF to Q image
    n     = Q64.shape[0]
    apix  = 1.0
    ky    = np.fft.fftfreq(n, apix)
    kx    = np.fft.rfftfreq(n, apix)
    KX, KY = np.meshgrid(kx, ky)
    K2D   = np.sqrt(KX**2 + KY**2) * 0.5 / (n * apix)  # scale to ~0–0.5 Å⁻¹ range
    df_A  = df_sl.value * 1e4
    Cs_A  = cs_sl.value * 1e7
    gamma2 = (-np.pi/2)*Cs_A*lam**3*K2D**4 + np.pi*lam*df_A*K2D**2
    ctf2d = -np.sin(phi_sl.value + gamma2)
    if bf_sl.value > 0:
        ctf2d *= np.exp(-bf_sl.value * K2D**2)

    Fq    = np.fft.rfftn(Q64)
    Q_ctf = np.fft.irfftn(Fq * ctf2d, Q64.shape)
    Q_noisy = Q_ctf + np.random.default_rng(1).standard_normal(Q64.shape)*0.3

    fig, axes = plt.subplots(1, 4, figsize=(16, 4.5))
    axes[0].plot(k, ctf_curve, color='steelblue', lw=1.8)
    axes[0].axhline(0, color='k', lw=0.8, ls='--')
    axes[0].set_xlabel("Spatial frequency (Å⁻¹)"); axes[0].set_ylabel("CTF")
    axes[0].set_title(f"CTF (Δf={df_sl.value} µm, V={volt_sl.value:.0f} kV)")
    axes[0].set_ylim(-1.2, 1.2); axes[0].grid(alpha=0.3)

    ps_ctf2d = np.log(np.abs(np.fft.fftshift(np.fft.fft2(Q_ctf)))**2 + 1)
    axes[1].imshow(Q64,     cmap='gray', origin='lower'); axes[1].axis('off')
    axes[1].set_title("True image")
    axes[2].imshow(Q_noisy, cmap='gray', origin='lower'); axes[2].axis('off')
    axes[2].set_title("CTF-modulated image")
    axes[3].imshow(ps_ctf2d, cmap='inferno', origin='lower'); axes[3].axis('off')
    axes[3].set_title("Power spectrum (Thon rings)")
    plt.tight_layout()
    with out_ctf:
        clear_output(wait=True); display(fig2img(fig))

for w in [df_sl, volt_sl, cs_sl, phi_sl, bf_sl]:
    w.observe(_update_ctf, names='value')
display(VBox([df_sl, volt_sl, cs_sl, phi_sl, bf_sl, out_ctf]))
_update_ctf()

12.8.1. CTF correction#

To recover the true image from a CTF-modulated observation, we need to correct for the CTF. Three standard methods:

Method 1 – Phase flipping: Multiply Fourier amplitudes by the sign of the CTF. Brings all amplitudes to positive, but does not correct their magnitudes.

Method 2 – Full CTF correction with threshold: Divide by the CTF, but ignore frequencies where $|\text{CTF}| < \epsilon$ (near zeros) to avoid division instability.

Method 3 – Wiener filter: Divide by CTF with regularisation: $$ \hat{F}(k) = \frac{F_\text{obs}(k) \cdot \text{CTF}(k)}{\text{CTF}(k)^2 + \text{SNR}^{-1}} $$

Show code cell source

Hide code cell source

corr_dd = Dropdown(options=['Phase flip', 'Full CTF (threshold)', 'Wiener filter'],
                   value='Wiener filter', description="Correction:", style=_sty)
snr_sl  = FloatSlider(value=1.0,   min=0.01, max=10.0, step=0.1,
                       description="SNR (Wiener)", style=_sty, layout=_sl)
thr_sl  = FloatSlider(value=0.05,  min=0.01, max=0.5,  step=0.01,
                       description="Threshold ε",  style=_sty, layout=_sl)
df2_sl  = FloatSlider(value=2.0,   min=0.5,  max=10.0, step=0.5,
                       description="Defocus (µm)", style=_sty, layout=_sl)
out_ctfc = Output()

def _update_ctfc(_=None):
    lam   = relativistic_lambda(300e3)
    n     = Q64.shape[0]; apix = 1.0
    ky    = np.fft.fftfreq(n, apix); kx = np.fft.rfftfreq(n, apix)
    KX, KY = np.meshgrid(kx, ky)
    K2D   = np.sqrt(KX**2 + KY**2) * 0.5 / (n * apix)
    df_A  = df2_sl.value * 1e4; Cs_A = 2.7e7
    gamma2 = (-np.pi/2)*Cs_A*lam**3*K2D**4 + np.pi*lam*df_A*K2D**2
    ctf2d = -np.sin(gamma2)

    rng   = np.random.default_rng(1)
    Fq    = np.fft.rfftn(Q64)
    F_obs = Fq * ctf2d + rng.standard_normal(Fq.shape)*0.3

    method = corr_dd.value
    if method == 'Phase flip':
        F_corr = F_obs * np.sign(ctf2d)
    elif method == 'Full CTF (threshold)':
        eps    = thr_sl.value
        F_corr = np.where(np.abs(ctf2d) >= eps, F_obs / ctf2d, 0.0)
    else:
        snr    = snr_sl.value
        F_corr = F_obs * ctf2d / (ctf2d**2 + 1.0/snr)

    Q_corr  = np.fft.irfftn(F_corr, Q64.shape)
    Q_ctf_i = np.fft.irfftn(F_obs,  Q64.shape)

    fig, axes = plt.subplots(1, 3, figsize=(12, 4))
    axes[0].imshow(Q64,     cmap='gray', origin='lower'); axes[0].axis('off')
    axes[0].set_title("True image")
    axes[1].imshow(Q_ctf_i, cmap='gray', origin='lower'); axes[1].axis('off')
    axes[1].set_title(f"CTF-modulated (Δf={df2_sl.value} µm)")
    axes[2].imshow(Q_corr,  cmap='gray', origin='lower'); axes[2].axis('off')
    axes[2].set_title(f"Corrected: {method}")
    plt.tight_layout()
    with out_ctfc:
        clear_output(wait=True); display(fig2img(fig))

for w in [corr_dd, snr_sl, thr_sl, df2_sl]:
    w.observe(_update_ctfc, names='value')
display(VBox([df2_sl, corr_dd, thr_sl, snr_sl, out_ctfc]))
_update_ctfc()

Question 4

Compare the three CTF correction methods. What are the trade-offs?

Phase flipping: what information is still lost?
Full CTF correction with threshold: what happens to frequencies near the CTF zeros?
Wiener filter: what is the optimal SNR for your simulated data?

12.9. Summary#

Topic	Key concept
Sinusoidal waves	$f(t) = A\cos(2\pi f t - \varphi)$; three parameters: $A$, $f$, $\varphi$
Fourier series	Any periodic signal = sum of sinusoids; coefficients via $a_k$, $b_k$ integrals
FFT	Fast algorithm ($O(N\log N)$) for discrete FT; frequency spectrum = $
Nyquist	Sampling at $\Delta t$ can recover frequencies up to $1/(2\Delta t)$
2D FT	Generalises 1D; amplitude spectrum = $
Convolution theorem	$\mathcal{F}\{f*g\} = \mathcal{F}\{f\}\cdot\mathcal{F}\{g\}$; filter in Fourier space
CTF	$-\sin(\gamma(k))$; phase reversals at zeros; correct by phase flip, division, or Wiener filter

Parameter	Symbol	Effect
Amplitude	\(A\)	Height of the wave
Frequency	\(f\)	Number of cycles per unit time
Phase	\(\varphi\)	Horizontal shift

Symbol	Quantity
\(k\)	spatial frequency (Å⁻¹)
\(C_s\)	spherical aberration coefficient
\(\lambda\)	relativistic electron wavelength
\(\Delta_f\)	defocus (positive = underfocus)
\(\Delta\varphi\)	additional phase shift (e.g. from phase plate)

Topic	Key concept
Sinusoidal waves	\(f(t) = A\cos(2\pi f t - \varphi)\); three parameters: \(A\), \(f\), \(\varphi\)
Fourier series	Any periodic signal = sum of sinusoids; coefficients via \(a_k\), \(b_k\) integrals
FFT	Fast algorithm (\(O(N\log N)\)) for discrete FT; frequency spectrum = $
Nyquist	Sampling at \(\Delta t\) can recover frequencies up to \(1/(2\Delta t)\)
2D FT	Generalises 1D; amplitude spectrum = $
Convolution theorem	\(\mathcal{F}\{f*g\} = \mathcal{F}\{f\}\cdot\mathcal{F}\{g\}\); filter in Fourier space
CTF	\(-\sin(\gamma(k))\); phase reversals at zeros; correct by phase flip, division, or Wiener filter

Practical: Fourier Optics and Image Processing

Contents

12. Practical: Fourier Optics and Image Processing#

12.1. Introduction#

12.1.1. Setup#

12.2. Part 1 – What is a wave?#

12.3. Part 2 – Fourier series#

12.3.1. The box function#

12.4. Part 3 – Frequency spectrum (1D FFT)#

12.4.1. Nyquist sampling theorem#

12.5. Part 4 – 2D Fourier analysis#

12.5.1. Properties#

12.5.2. 2D FT of real images#

12.5.3. Translation and rotation properties#

12.6. Part 5 – The convolution theorem#

12.6.1. Common kernels#

12.6.2. Interactive kernel#

12.7. Part 6 – Low-pass and high-pass filtering in Fourier space#

12.8. Part 7 – The contrast transfer function (CTF)#

12.8.1. CTF correction#

12.9. Summary#