---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.10.3
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
mystnb:
  execution_timeout: 180
---

(ch:fourier-transform)=
# Fourier Analysis in Imaging

The Fourier transform is the single most important mathematical tool in quantitative imaging. It provides a language for describing images in terms of their **spatial frequency content** — which periodicities are present, at what amplitude, and at what phase. Understanding Fourier analysis is essential for making sense of resolution, contrast, filtering, and the tomographic reconstruction methods described in {numref}`ch:tomography`.

(sec:ft-1d)=
## The 1D Fourier Transform

Any periodic signal can be decomposed into a sum of sinusoids — this is the central insight of Fourier analysis. For a non-periodic signal, the continuous **Fourier transform** generalises this to a continuous integral:

$$
\hat{f}(\nu) = \int_{-\infty}^{\infty} f(x)\, e^{-2\pi i \nu x}\, dx
$$ (eq:ft-1d)

where $f(x)$ is the signal in real space (e.g., position $x$), $\hat{f}(\nu)$ is its Fourier transform in frequency space, and $\nu$ is the **spatial frequency** (cycles per unit length). The inverse transform recovers the original signal:

$$
f(x) = \int_{-\infty}^{\infty} \hat{f}(\nu)\, e^{2\pi i \nu x}\, d\nu
$$ (eq:ift-1d)

The Fourier transform $\hat{f}(\nu)$ is in general complex. Its **magnitude** $|\hat{f}(\nu)|$ is the amplitude spectrum, and $\arg \hat{f}(\nu)$ is the phase spectrum. The **power spectrum** is $|\hat{f}(\nu)|^2$.

### Key Properties

| Property | Real space | Fourier space |
|---|---|---|
| Linearity | $a f + b g$ | $a\hat{f} + b\hat{g}$ |
| Shift | $f(x - x_0)$ | $e^{-2\pi i \nu x_0}\hat{f}(\nu)$ — phase shift only |
| Scale | $f(ax)$ | $\frac{1}{\|a\|}\hat{f}(\nu/a)$ — wide → narrow, narrow → wide |
| Convolution | $f * g$ | $\hat{f} \cdot \hat{g}$ — convolution becomes multiplication |
| Parseval | $\int \|f\|^2 dx$ | $\int \|\hat{f}\|^2 d\nu$ — energy is conserved |

The **convolution theorem** is particularly important for imaging: the action of any linear, shift-invariant imaging system (a lens, a detector, a filter) can be described by its **point spread function (PSF)** in real space, or equivalently by its **transfer function** (the Fourier transform of the PSF) in frequency space. Multiplying in frequency space is far more efficient than convolving in real space for large images.

The interactive below lets you build a 1D signal from sinusoidal components and see the resulting power spectrum. Adding components of different frequencies, amplitudes, and phases shows how complex waveforms arise from simple harmonic building blocks.

```{admonition} Interactive element
:class: tip
Click **Live Code** to activate, expand the **Show code** toggle, and click ▶ to run. Use the sliders to adjust the frequencies and amplitudes of three sine-wave components. The signal panel shows the resulting waveform; the spectrum panel shows its power spectrum.
```

```{code-cell} ipython3
:tags: [hide-input]
import io
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from ipywidgets import FloatSlider, IntSlider, VBox, HBox, Output, Label
from IPython.display import display, Image, clear_output

_x = np.linspace(0, 1, 512)
_freqs = np.fft.rfftfreq(512, d=1/512)

def draw_1d_ft(f1=1, A1=1.0, f2=3, A2=0.5, f3=7, A3=0.25):
    sig = (A1 * np.sin(2*np.pi*f1*_x) +
           A2 * np.sin(2*np.pi*f2*_x) +
           A3 * np.sin(2*np.pi*f3*_x))
    ps = np.abs(np.fft.rfft(sig))**2 / 512

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))
    ax1.plot(_x, sig, color='steelblue', linewidth=2)
    ax1.set_xlabel("Position x", fontsize=10)
    ax1.set_ylabel("Signal amplitude", fontsize=10)
    ax1.set_title("Signal (sum of three sinusoids)", fontsize=10)
    ax1.grid(alpha=0.3)

    ax2.stem(_freqs[:30], ps[:30], linefmt='C1-', markerfmt='C1o', basefmt='k-')
    ax2.set_xlabel("Spatial frequency (cycles)", fontsize=10)
    ax2.set_ylabel("Power", fontsize=10)
    ax2.set_title("Power spectrum", fontsize=10)
    ax2.grid(alpha=0.3)

    fig.tight_layout()
    buf = io.BytesIO(); fig.savefig(buf, format='png', dpi=96); buf.seek(0)
    display(Image(data=buf.read())); plt.close(fig)

style = {"description_width": "80px"}
sl_kw = {"layout": {"width": "320px"}}
f1_sl = IntSlider(value=1,  min=1, max=15, step=1,  description="Freq 1", style=style, **sl_kw)
A1_sl = FloatSlider(value=1.0, min=0.0, max=2.0, step=0.1, description="Amp 1",  style=style, **sl_kw)
f2_sl = IntSlider(value=3,  min=1, max=20, step=1,  description="Freq 2", style=style, **sl_kw)
A2_sl = FloatSlider(value=0.5, min=0.0, max=2.0, step=0.1, description="Amp 2",  style=style, **sl_kw)
f3_sl = IntSlider(value=7,  min=1, max=25, step=1,  description="Freq 3", style=style, **sl_kw)
A3_sl = FloatSlider(value=0.25, min=0.0, max=2.0, step=0.1, description="Amp 3", style=style, **sl_kw)

out_1d = Output()
def update_1d(_=None):
    with out_1d:
        clear_output(wait=True)
        draw_1d_ft(f1_sl.value, A1_sl.value, f2_sl.value, A2_sl.value, f3_sl.value, A3_sl.value)
for s in [f1_sl, A1_sl, f2_sl, A2_sl, f3_sl, A3_sl]:
    s.observe(update_1d, names='value')

display(VBox([HBox([VBox([f1_sl, A1_sl]), VBox([f2_sl, A2_sl]), VBox([f3_sl, A3_sl])]), out_1d]))
update_1d()
```

(sec:ft-2d)=
## The 2D Fourier Transform in Images

For a 2D image $f(x, y)$, the Fourier transform extends naturally to two dimensions:

$$
\hat{f}(k_x, k_y) = \iint f(x, y)\, e^{-2\pi i (k_x x + k_y y)}\, dx\, dy
$$ (eq:ft-2d)

Here $(k_x, k_y)$ are the two components of the **spatial frequency vector** $\mathbf{k}$. The magnitude $|\mathbf{k}| = \sqrt{k_x^2 + k_y^2}$ is the radial spatial frequency (1/resolution), and the direction of $\mathbf{k}$ indicates the orientation of the corresponding sinusoidal pattern in the image.

### Interpreting the 2D FT

After applying `fftshift` to move the DC component (zero frequency) to the centre, the 2D FT has a characteristic structure:

- **Centre** (DC component): the average value of the image
- **Near centre** (low frequencies): coarse structure, global intensity variation
- **Away from centre** (high frequencies): fine detail, sharp edges
- **Bright lines or arcs**: features with a strong periodic structure at that frequency and orientation

Several important image features have recognisable signatures in the 2D FT:
- A **horizontal line** in the image appears as a **vertical line** in the FT (perpendicular to the feature)
- A **dot** (isotropic blob) appears as a **ring** in the FT
- A **ring** in the image appears as a **ring** in the FT (at the reciprocal radius)

```{code-cell} ipython3
:tags: [remove-input]

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

def show_ft(img, ax_img, ax_ft, title):
    ax_img.imshow(img, cmap='gray'); ax_img.axis('off'); ax_img.set_title(title, fontsize=9)
    FT = np.log1p(np.abs(np.fft.fftshift(np.fft.fft2(img))))
    ax_ft.imshow(FT, cmap='inferno'); ax_ft.axis('off')
    ax_ft.set_title(f"2D FT (log power)", fontsize=9)

n = 128
y, x = np.mgrid[0:n, 0:n].astype(float)
cx, cy = n/2, n/2
examples = []

# Horizontal line
img = np.zeros((n, n)); img[n//2-2:n//2+2, :] = 1; examples.append((img, "Horizontal line"))
# Vertical line
img = np.zeros((n, n)); img[:, n//2-2:n//2+2] = 1; examples.append((img, "Vertical line"))
# Ring
R = np.sqrt((x-cx)**2 + (y-cy)**2)
img = ((R > 0.25*n) & (R < 0.30*n)).astype(float); examples.append((img, "Ring"))
# Two blobs
img = np.zeros((n, n))
img += np.exp(-((x-cx-20)**2 + (y-cy)**2)/(2*8**2))
img += np.exp(-((x-cx+20)**2 + (y-cy)**2)/(2*8**2))
examples.append((img, "Two blobs (spacing 40 px)"))
# Random noise
img = np.random.default_rng(0).standard_normal((n, n)); examples.append((img, "White noise"))
# Low-pass filtered noise
ft = np.fft.fft2(img)
kx = np.fft.fftfreq(n); ky = np.fft.fftfreq(n)
KX, KY = np.meshgrid(kx, ky)
filt = np.exp(-(KX**2+KY**2)/(2*0.05**2))
img = np.real(np.fft.ifft2(ft * filt)); examples.append((img, "Low-pass filtered noise"))

fig, axes = plt.subplots(len(examples), 2, figsize=(8, 16))
for row, (img, title) in enumerate(examples):
    show_ft(img, axes[row, 0], axes[row, 1], title)

plt.suptitle("2D Fourier transform of different image features", fontsize=10, y=1.01)
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')
```

(sec:ft-filtering)=
## Filtering in Fourier Space

Because convolution in real space equals multiplication in Fourier space (the convolution theorem), **spatial filtering** — suppressing or enhancing certain frequencies — is most easily implemented by multiplying the 2D FT by a **filter function** $H(k_x, k_y)$ and then inverse-transforming:

$$
g(x, y) = \mathcal{F}^{-1}\!\left[ H(k_x, k_y) \cdot \hat{f}(k_x, k_y) \right]
$$ (eq:filtering)

Common filter types:

**Low-pass filter**: $H = 1$ for $|\mathbf{k}| < k_\text{cut}$, $H = 0$ otherwise. Attenuates high frequencies → smooths the image, removes noise but also blurs detail. In cryo-EM, low-pass filtering is used to remove high-frequency noise from noisy micrographs.

**High-pass filter**: $H = 0$ for $|\mathbf{k}| < k_\text{cut}$, $H = 1$ otherwise. Attenuates low frequencies → enhances edges and fine detail, but amplifies high-frequency noise. The **ramp filter** used in filtered backprojection (see {numref}`ch:tomography`) is a radial high-pass filter $H = |\mathbf{k}|$.

**Band-pass filter**: passes a range of frequencies. Used for specific feature enhancement (e.g., suppressing both very low and very high frequencies while retaining intermediate-scale structures).

A sharp cutoff in Fourier space (a "top-hat" filter) causes **ringing artifacts** (Gibbs phenomenon) in real space. In practice, filters are smoothed with a Gaussian or cosine taper. The choice of filter shape involves a trade-off between frequency selectivity and real-space artifacts.

```{admonition} Interactive element
:class: tip
Click **Live Code** to activate, expand the **Show code** toggle, and click ▶ to run. Choose a filter type and adjust the cutoff frequency. Watch how the image changes and what the Fourier-space filter looks like.
```

```{code-cell} ipython3
:tags: [hide-input]
import io
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from ipywidgets import FloatSlider, Dropdown, VBox, HBox, Output
from IPython.display import display, Image, clear_output

def _make_ti(n=128):
    y, x = np.mgrid[0:n, 0:n].astype(float)
    cx, cy = n/2, n/2
    img = np.zeros((n, n))
    R = np.sqrt((x-cx)**2 + (y-cy)**2)
    img += 0.8 * ((R > 0.22*n) & (R < 0.28*n)).astype(float)
    img += 0.9 * (np.abs(y - cy*1.5) < n*0.02).astype(float)
    img += 0.7 * (np.sqrt((x-cx*1.6)**2+(y-cy*0.5)**2) < n*0.06).astype(float)
    img += 0.05 * np.random.default_rng(1).standard_normal((n, n))
    return np.clip(img, 0, 1)

_fil_img = _make_ti(128)
_fil_n = 128
kx = np.fft.fftfreq(_fil_n); ky = np.fft.fftfreq(_fil_n)
_KX, _KY = np.meshgrid(kx, ky); _K = np.sqrt(_KX**2 + _KY**2)
_FT = np.fft.fft2(_fil_img)

def draw_filter(filter_type='Low-pass', cutoff=0.15):
    if filter_type == 'Low-pass':
        H = np.exp(-(_K/cutoff)**4)
        label = f"Low-pass (cutoff={cutoff:.2f})"
    elif filter_type == 'High-pass':
        H = 1 - np.exp(-(_K/cutoff)**4)
        label = f"High-pass (cutoff={cutoff:.2f})"
    else:
        H = _K / _K.max()
        label = "Ramp filter |k|"

    filtered = np.real(np.fft.ifft2(_FT * np.fft.ifftshift(H)))
    H_display = np.fft.fftshift(H)

    fig, axes = plt.subplots(1, 4, figsize=(16, 4.2))
    axes[0].imshow(_fil_img, cmap='gray'); axes[0].axis('off')
    axes[0].set_title("Original image", fontsize=9)
    axes[1].imshow(np.log1p(np.abs(np.fft.fftshift(_FT))), cmap='inferno'); axes[1].axis('off')
    axes[1].set_title("2D FT (log power)", fontsize=9)
    axes[2].imshow(H_display, cmap='RdBu_r', vmin=0, vmax=1); axes[2].axis('off')
    axes[2].set_title(f"Filter H(k)\n{label}", fontsize=9)
    axes[3].imshow(filtered, cmap='gray'); axes[3].axis('off')
    axes[3].set_title("Filtered image", fontsize=9)

    fig.suptitle(f"Fourier filtering: {label}", fontsize=10)
    fig.tight_layout()
    buf = io.BytesIO(); fig.savefig(buf, format='png', dpi=96); buf.seek(0)
    display(Image(data=buf.read())); plt.close(fig)

filter_dd = Dropdown(options=['Low-pass', 'High-pass', 'Ramp'],
                     value='Low-pass', description='Filter type:',
                     style={"description_width": "110px"}, layout={"width": "280px"})
cutoff_sl = FloatSlider(value=0.15, min=0.02, max=0.48, step=0.01,
                        description="Cutoff k",
                        style={"description_width": "90px"}, layout={"width": "340px"})
out_fil = Output()
def update_fil(_=None):
    with out_fil:
        clear_output(wait=True)
        draw_filter(filter_dd.value, cutoff_sl.value)
for w in [filter_dd, cutoff_sl]:
    w.observe(update_fil, names='value')
display(VBox([HBox([filter_dd, cutoff_sl]), out_fil]))
update_fil()
```

(sec:ft-resolution)=
## Resolution and the Frequency Limit

In any digital image, the maximum representable spatial frequency is set by the **pixel size** (the Nyquist limit): if the pixel size is $d$ (in Å/pixel), the highest frequency is $\nu_\text{max} = 1/(2d)$ Å⁻¹. This is the **Nyquist frequency**. Features finer than $2d$ cannot be represented in the image — they alias back to lower frequencies.

In cryo-EM, the nominal resolution is often defined as the spatial frequency at which the signal-to-noise ratio drops to 1, typically measured by the **Fourier Shell Correlation** (see {numref}`subsec:spa-fsc`). A map at "3 Å resolution" means the FSC drops below the 0.143 criterion at $1/3$ Å⁻¹.

(sec:ft-psf-ctf)=
## The Point Spread Function and the CTF

Every real imaging system has a **point spread function (PSF)**: the image of a perfect point source. By the convolution theorem, the image of any object $f$ through a system with PSF $h$ is:

$$
g = f * h \quad \Longleftrightarrow \quad \hat{g} = \hat{f} \cdot \hat{h}
$$

The Fourier transform of the PSF is the **optical transfer function (OTF)**. For an electron microscope operating in phase-contrast mode, the OTF reduces to the **contrast transfer function (CTF)**. The full derivation — starting from the weak phase object approximation and the wave-propagation picture of defocus-induced contrast — is given in {numref}`sec:tem-image-formation`. The result is:

$$
\mathrm{CTF}(\nu) = -\sin\!\left[\pi\lambda\Delta f\,\nu^2 - \frac{\pi}{2}C_s\lambda^3\nu^4\right] \cdot e^{-B\nu^2/4}
$$

where $\Delta f$ is the underfocus, $C_s$ the spherical aberration coefficient, $\lambda$ the electron wavelength, and $e^{-B\nu^2/4}$ the B-factor envelope from finite coherence. The CTF oscillates between $-1$ and $+1$; zeros — the **Thon rings** visible in the 2D power spectrum of a micrograph — mark frequencies where no information is transferred. CTF correction (phase flipping, Wiener filter, or multiplicity weighting across many defocus values) is required before coherent particle averaging in SPA.

```{code-cell} ipython3
:tags: [remove-input]

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

def _lam_A(kV=300):
    V = kV * 1e3
    m0, ec, c, h = 9.109e-31, 1.602e-19, 2.998e8, 6.626e-34
    return h / np.sqrt(2*m0*ec*V*(1+ec*V/(2*m0*c**2))) * 1e10

def _ctf(nu, df_um, Cs_mm, B=0, kV=300):
    lam = _lam_A(kV)
    df = df_um * 1e4; Cs = Cs_mm * 1e7
    return -np.sin(np.pi*lam*df*nu**2 - 0.5*np.pi*Cs*lam**3*nu**4) * np.exp(-B*nu**2/4)

nu_p = np.linspace(0.002, 0.44, 2000)
n2 = 280
nu_lin = np.linspace(-0.44, 0.44, n2)
NX, NY = np.meshgrid(nu_lin, nu_lin)
CTF_2d = _ctf(np.sqrt(NX**2+NY**2), 2.0, 2.7, B=150)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

defoci = [0.5, 1.0, 2.0, 3.0]
cols   = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
for df, col in zip(defoci, cols):
    ax1.plot(nu_p, _ctf(nu_p, df, 2.7, B=0), color=col, linewidth=1.5, label=f'Δf={df} μm')
ax1.plot(nu_p, _ctf(nu_p, 2.0, 2.7, B=150), color='#9467bd', linewidth=2,
         linestyle='--', label='Δf=2 μm, B=150 Ų')
ax1.axhline(0, color='k', linewidth=0.6); ax1.set_ylim(-1.25, 1.45)
ax1.set_xlabel('Spatial frequency ν (Å⁻¹)', fontsize=10)
ax1.set_ylabel('CTF(ν)', fontsize=10)
ax1.set_title('CTF(ν) at different defocus values  (Cs=2.7 mm, 300 kV)', fontsize=10)
ax1.legend(fontsize=9)
ax_top = ax1.twiny(); ax_top.set_xlim(ax1.get_xlim())
ticks = [0.05, 0.10, 0.20, 0.33]
ax_top.set_xticks(ticks); ax_top.set_xticklabels([f'{1/t:.0f}' for t in ticks], fontsize=8)
ax_top.set_xlabel('Resolution (Å)', fontsize=8)

im = ax2.imshow(CTF_2d, cmap='RdBu_r', vmin=-1, vmax=1,
                extent=[-0.44, 0.44, -0.44, 0.44], origin='lower')
ax2.set_xlabel('k_x (Å⁻¹)', fontsize=10); ax2.set_ylabel('k_y (Å⁻¹)', fontsize=10)
ax2.set_title('2D CTF — Thon rings  (Δf=2 μm, Cs=2.7 mm, B=150 Ų, 300 kV)', fontsize=10)
plt.colorbar(im, ax=ax2, label='CTF', shrink=0.85)

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')
```

```{admonition} Interactive element
:class: tip
Click **Live Code** to activate, expand the **Show code** toggle, and click ▶ to run. Drag the defocus slider to shift the CTF zeros; increase B-factor to see the high-frequency envelope collapse; compare 120/200/300 kV.
```

```{code-cell} ipython3
:tags: [hide-input]
import io
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from ipywidgets import FloatSlider, Dropdown, VBox, HBox, Output
from IPython.display import display, Image, clear_output

def _lam_A(kV=300):
    V = kV * 1e3
    m0, ec, c, h = 9.109e-31, 1.602e-19, 2.998e8, 6.626e-34
    return h / np.sqrt(2*m0*ec*V*(1+ec*V/(2*m0*c**2))) * 1e10

def _ctf(nu, df_um, Cs_mm, B=0, kV=300):
    lam = _lam_A(kV)
    df = df_um * 1e4; Cs = Cs_mm * 1e7
    return -np.sin(np.pi*lam*df*nu**2 - 0.5*np.pi*Cs*lam**3*nu**4) * np.exp(-B*nu**2/4)

_nu1d = np.linspace(0.002, 0.50, 2000)
_n2 = 220
_nu_lin = np.linspace(-0.50, 0.50, _n2)
_NX2d, _NY2d = np.meshgrid(_nu_lin, _nu_lin)
_NU2d = np.sqrt(_NX2d**2 + _NY2d**2)

def draw_ctf(df_um, Cs_mm, B, kV):
    ctf1d = _ctf(_nu1d, df_um, Cs_mm, B, kV)
    env1d = np.exp(-B * _nu1d**2 / 4)
    ctf2d = _ctf(_NU2d, df_um, Cs_mm, B, kV)
    lam = _lam_A(kV)
    nu_zero = np.sqrt(1/(lam*df_um*1e4)) if df_um > 0 else np.nan

    fig, axes = plt.subplots(1, 3, figsize=(16, 5))

    ax = axes[0]
    ax.plot(_nu1d, ctf1d, color='#2c7bb6', linewidth=2, label='CTF(ν)')
    ax.plot(_nu1d,  env1d, color='#d62728', linewidth=1.5, linestyle='--', label='envelope')
    ax.plot(_nu1d, -env1d, color='#d62728', linewidth=1.5, linestyle='--')
    ax.axhline(0, color='k', linewidth=0.6)
    if np.isfinite(nu_zero) and nu_zero < 0.5:
        ax.axvline(nu_zero, color='orange', lw=1.2, ls=':', label=f'1st zero ≈{1/nu_zero:.0f} Å')
    ax.set_xlabel('ν (Å⁻¹)', fontsize=10); ax.set_ylabel('CTF(ν)', fontsize=10)
    ax.set_title(f'1D CTF  Δf={df_um:.1f} μm, Cs={Cs_mm:.1f} mm, B={B:.0f} Ų, {kV:.0f} kV', fontsize=10)
    ax.legend(fontsize=9); ax.set_ylim(-1.3, 1.5)
    ax_top = ax.twiny(); ax_top.set_xlim(ax.get_xlim())
    tks = [0.05, 0.10, 0.20, 0.33]
    ax_top.set_xticks(tks); ax_top.set_xticklabels([f'{1/t:.0f}' for t in tks], fontsize=8)
    ax_top.set_xlabel('Resolution (Å)', fontsize=9)

    ax = axes[1]
    im = ax.imshow(ctf2d, cmap='RdBu_r', vmin=-1, vmax=1,
                   extent=[-0.5, 0.5, -0.5, 0.5], origin='lower')
    ax.set_xlabel('k_x (Å⁻¹)', fontsize=10); ax.set_ylabel('k_y (Å⁻¹)', fontsize=10)
    ax.set_title('2D CTF — Thon rings', fontsize=10)
    plt.colorbar(im, ax=ax, label='CTF', shrink=0.85)

    ax = axes[2]
    rng = np.random.default_rng(42)
    noise_ft = np.fft.fft2(rng.standard_normal((_n2, _n2)))
    img = np.real(np.fft.ifft2(noise_ft * np.fft.ifftshift(ctf2d)))
    ax.imshow(img, cmap='gray', vmin=-3*img.std(), vmax=3*img.std())
    ax.axis('off')
    ax.set_title('Simulated micrograph (CTF × white noise)', fontsize=10)

    fig.tight_layout()
    buf = io.BytesIO(); fig.savefig(buf, format='png', dpi=96); buf.seek(0); plt.close(fig)
    display(Image(data=buf.read()))

style = {'description_width': '110px'}
sl = {'layout': {'width': '320px'}}
df_sl  = FloatSlider(value=2.0, min=0.2, max=5.0, step=0.1, description='Defocus (μm)', style=style, **sl)
Cs_sl  = FloatSlider(value=2.7, min=0.0, max=4.0, step=0.1, description='Cs (mm)',      style=style, **sl)
B_sl   = FloatSlider(value=100, min=0,   max=400,  step=10,  description='B-factor (Ų)',style=style, **sl)
kV_dd  = Dropdown(options=[120, 200, 300], value=300, description='Voltage (kV)',
                  style=style, layout={'width': '260px'})

out_ctf = Output()
def update_ctf(_=None):
    with out_ctf:
        clear_output(wait=True)
        draw_ctf(df_sl.value, Cs_sl.value, B_sl.value, kV_dd.value)
for w in [df_sl, Cs_sl, B_sl, kV_dd]:
    w.observe(update_ctf, names='value')
display(VBox([HBox([VBox([df_sl, Cs_sl]), VBox([B_sl, kV_dd])]), out_ctf]))
update_ctf()
```
(sec:ft-slice-preview)=
## From Fourier Transforms to Tomographic Reconstruction

The Fourier Slice Theorem ({numref}`sec:tomo-projection-theorem`) is a direct consequence of the Fourier transform's properties. Consider projecting a 2D image $f(x,y)$ along the $y$-axis to get a 1D profile $P(x) = \int f(x,y)\, dy$. Taking the 1D Fourier transform of $P$:

$$
\hat{P}(k_x) = \int P(x)\, e^{-2\pi i k_x x}\, dx
= \int\!\!\int f(x,y)\, e^{-2\pi i k_x x}\, dx\, dy
= \hat{F}(k_x, 0)
$$

The last equality shows that $\hat{P}(k_x)$ is simply the 2D Fourier transform of $f$ evaluated on the line $k_y = 0$ — the central horizontal slice. For a projection at angle $\theta$, the same argument shows that the FT of the projection fills in the central slice of the 2D FT at angle $\theta$.

This has a direct practical implication: **to reconstruct a 2D image from its projections, one can fill in the 2D Fourier space from the 1D FT of each projection, and then take the inverse 2D FT**. This is the basis of reconstruction by Fourier inversion — the foundation of all modern tomographic reconstruction algorithms.

The interactive in {numref}`sec:tomo-projection-theorem` demonstrates this theorem live: as you change the tilt angle, the red line in the 2D Fourier space rotates to show exactly which frequencies each new projection measures. Building up from many angles progressively fills Fourier space, and the missing wedge becomes visible as the angular range is limited.

(sec:ft-summary)=
## Summary

The key ideas in this chapter:

1. **Any image can be decomposed** into a sum of sinusoidal spatial-frequency components, characterised by frequency, amplitude, and phase.
2. **The 2D FT is the spatial-frequency representation** of an image. Low frequencies encode global structure; high frequencies encode fine detail and edges.
3. **Filtering** is multiplication in Fourier space — far more efficient than convolution in real space for large images.
4. **Resolution** is set by the highest reliably measured frequency. In cryo-EM, radiation dose and detector noise limit the usable resolution.
5. **The PSF/OTF (and CTF in EM)** fully characterise the imaging system in frequency space. Image formation is multiplication of the object's FT with the OTF.
6. **The Fourier Slice Theorem** — each projection fills one central line through Fourier space — is the mathematical foundation of tomographic reconstruction.