Practical: Electron Tomography

14. Practical: Electron Tomography#

14.1. Introduction#

Unlike single-particle analysis, electron tomography (cryo-ET) reconstructs a 3D volume from images of a single specimen tilted to many angles by physically rotating the sample stage. We work through the problem in 2D: the unknown object is a 2D image, and we acquire 1D projections at various tilt angles, then reconstruct the 2D image from those projections.

Axis convention: The electron beam travels along $z$. The tilt axis is $y$ (pointing through the imaging plane). Tilting the sample about $y$ rotates the $xz$-plane, so each tilt angle $\alpha$ produces a 1D projection along $x$ of the rotated image.

A typical cryo-ET experiment collects tilt images from $-60°$ to $+60°$ in steps of $1°$–$3°$, limited by grid-bar occlusion and increased sample thickness at high tilt.

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.

14.1.1. Helper functions#

Show code cell source

Hide code cell source

import io, warnings
import numpy as np
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from scipy.ndimage import rotate as nd_rotate, gaussian_filter
from ipywidgets import IntSlider, FloatSlider, VBox, HBox, Layout, Output
from IPython.display import display, Image, clear_output
warnings.filterwarnings('ignore')

_sl  = Layout(width="420px")
_sty = {"description_width": "160px"}

def fig2img(fig, dpi=100):
    buf = io.BytesIO()
    fig.savefig(buf, format='png', dpi=dpi, bbox_inches='tight')
    buf.seek(0); plt.close(fig)
    return Image(data=buf.read())

# ── test image ────────────────────────────────────────────────────────────────

def make_2D_test_image(size=100):
    """
    2D phantom representing a cross-section of a biological specimen:
      - ring  = microtubule cross-section
      - lines = membrane sheets
      - dot   = protein complex
    Pure-numpy implementation (no skimage needed).
    """
    y, x = np.mgrid[0:size, 0:size].astype(float)
    img  = np.zeros((size, size))

    # Ring (microtubule): hollow circle, top-left quadrant
    cx1, cy1 = size*0.33, size*0.33
    R1 = np.sqrt((x - cx1)**2 + (y - cy1)**2)
    img += 0.9 * ((R1 > size*0.12) & (R1 < size*0.18))

    # Horizontal membrane sheet
    cy2 = size*0.68
    img += 0.85 * (np.abs(y - cy2) < size*0.018)

    # Vertical membrane sheet
    cx3 = size*0.67
    img += 0.85 * (np.abs(x - cx3) < size*0.018) * ((y > size*0.4) & (y < size*0.9))

    # Dot (protein complex), top-right quadrant
    cx4, cy4 = size*0.72, size*0.35
    R4 = np.sqrt((x - cx4)**2 + (y - cy4)**2)
    img += 0.75 * (R4 < size*0.06)

    return np.clip(img, 0, 1)[::-1]

# ── Radon transform ───────────────────────────────────────────────────────────

def radon(im, angles):
    """
    Compute 1D projections of im at each tilt angle.
    Rotation is applied to the image; projection is the column sum.
    Returns sinogram of shape (len(angles), im.shape[1]).
    """
    sino = np.zeros((len(angles), im.shape[1]))
    for k, ang in enumerate(angles):
        rotated    = nd_rotate(im, ang, reshape=False, cval=0)
        sino[k, :] = rotated.sum(axis=0)
    return sino

# ── backprojection ────────────────────────────────────────────────────────────

def repeat_line(line, N):
    """Repeat 1D line along a new axis N times → shape (N, len(line))."""
    return np.tile(line[np.newaxis, :], (N, 1))

def backproject(sinogram, angles):
    """
    Simple (unfiltered) backprojection.
    Each 1D projection is repeated into a 2D slab, rotated, and summed.
    """
    n = sinogram.shape[1]
    recon = np.zeros((n, n))
    for proj, ang in zip(sinogram, angles):
        slab    = repeat_line(proj / n, n)
        recon  += nd_rotate(slab, -ang, reshape=False, cval=0)
    return recon / max(len(angles), 1)

# ── filtered backprojection ───────────────────────────────────────────────────

def ramp_filter(sinogram):
    """Apply ramp filter |ν| to each projection in the sinogram."""
    n       = sinogram.shape[1]
    freqs   = np.fft.rfftfreq(n)
    ramp    = np.abs(freqs)
    sino_f  = np.zeros_like(sinogram)
    for i, proj in enumerate(sinogram):
        fft_proj    = np.fft.rfft(proj)
        sino_f[i]   = np.fft.irfft(fft_proj * ramp, n)
    return sino_f

def filtered_backproject(sinogram, angles):
    """Filtered backprojection with ramp filter."""
    sino_f = ramp_filter(sinogram)
    return backproject(sino_f, angles)

# ── Fourier power spectrum ────────────────────────────────────────────────────

def power_spectrum(im):
    F   = np.fft.fft2(im)
    F   = np.fft.fftshift(F)
    return np.log(np.abs(F)**2 + 1)

print("Helper functions loaded.")

14.2. The 2D test image#

We work with a 2D phantom that mimics a biological cross-section: a hollow ring (microtubule), two line segments (membrane sheets), and a small dot (protein complex).

14.3. Tilting the sample and acquiring projections#

To acquire a tilt series, the specimen is physically rotated about the $y$-axis. At each tilt angle $\alpha$, the transmitted electrons form a 1D projection (the column sum of the rotated image). The interactive below shows the specimen at a chosen tilt angle together with the resulting 1D projection profile.

14.4. The sinogram#

Collecting projections over the full tilt range and stacking them row by row produces the sinogram — so named because a point object traces a sinusoidal curve across it as the tilt angle varies.

14.5. Task 1 – Backprojection#

The simplest reconstruction method is backprojection: for each projection at angle $\alpha$, repeat the 1D line into a 2D slab, rotate it back by $-\alpha$, and sum over all angles.

\[ A_\text{BP}(x,z) = \sum_\alpha p_\alpha\!\left(x\cos\alpha + z\sin\alpha\right) \]

This is correct in principle but overweights low-spatial-frequency components (every point in the image receives contributions from every projection), producing a blurry, star-shaped reconstruction artifact.

Task 1

Take a 1D projection projection.
Repeat it $N$ times along a new axis using repeat_line to get a 2D slab of shape $(N, N)$.
Rotate the slab back by $-\alpha$ using nd_rotate.
Sum over all angles. Divide by the number of angles.

14.6. Task 2 – Filtered backprojection#

The blurring artifact of plain backprojection is caused by the over-sampling of low frequencies: every projection adds a constant along its direction, contributing to the low-frequency region of Fourier space more than to high frequencies. The fix is to pre-filter each projection with a ramp filter $|\nu|$ before backprojecting. This exactly counteracts the $1/|\nu|$ weighting introduced by backprojection.

\[ A_\text{FBP}(x,z) = \sum_\alpha \bigl[p_\alpha * h\bigr]\!\left(x\cos\alpha + z\sin\alpha\right), \qquad \hat h(\nu) = |\nu| \]

Task 2

Take a 1D projection.
Compute its Fourier transform using np.fft.rfft.
Multiply by the ramp filter $|\nu|$ (use np.fft.rfftfreq for the frequency axis).
Inverse-Fourier-transform back.
Use the filtered projections for backprojection.

14.7. The Fourier Slice Theorem#

The projection-slice theorem (also called the Fourier slice theorem) states that the 1D Fourier transform of a projection at angle $\alpha$ equals a central slice through the 2D Fourier transform of the object at the same angle:

\[ \mathcal{F}_1\{p_\alpha\}(\nu) = \mathcal{F}_2\{A\}(\nu\cos\alpha,\, \nu\sin\alpha) \]

Each projection therefore fills one line (passing through the origin) in 2D Fourier space. To reconstruct the full 3D (or 2D) object, we need enough projections to cover Fourier space densely.

The interactive below shows, for a chosen tilt angle:

The rotated specimen and its 1D projection
The 2D Fourier transform of the specimen with the corresponding central slice highlighted
A verification that the 1D FT of the projection matches the highlighted central slice

Show code cell source

Hide code cell source

from scipy.ndimage import map_coordinates

fst_sl  = IntSlider(value=0, min=-70, max=70, step=1,
                     description="Tilt angle (°)", style=_sty, layout=_sl)
out_fst = Output()

_n_fst  = im.shape[0]
_F2_obj = np.fft.fftshift(np.fft.fft2(im))

def _extract_central_slice(F2, angle_deg):
    """Extract central slice from 2D FT at given angle via bilinear interpolation."""
    n    = F2.shape[0]
    ang  = np.deg2rad(angle_deg)
    t    = np.linspace(-n//2, n//2, n)
    xs   = t * np.cos(ang) + n//2
    ys   = t * np.sin(ang) + n//2
    re   = map_coordinates(F2.real, [ys, xs], order=1, mode='nearest')
    im_p = map_coordinates(F2.imag, [ys, xs], order=1, mode='nearest')
    return re + 1j*im_p

def _update_fst(_=None):
    ang  = fst_sl.value
    rot  = nd_rotate(im, ang, reshape=False, cval=0)
    proj = rot.sum(axis=0) / _n_fst

    # 1D FT of projection
    fft1d      = np.fft.fftshift(np.fft.fft(proj))
    fft1d_amp  = np.log(np.abs(fft1d) + 1)

    # Central slice from 2D FT
    cslice     = _extract_central_slice(_F2_obj, ang)
    cslice_amp = np.log(np.abs(cslice) + 1)

    # Draw slice line on 2D FT image
    ps2 = np.log(np.abs(_F2_obj) + 1)
    n   = _n_fst
    t   = np.linspace(-n//2, n//2, n)
    ang_r = np.deg2rad(ang)
    xs_line = t * np.cos(ang_r) + n//2
    ys_line = t * np.sin(ang_r) + n//2

    fig, axes = plt.subplots(1, 4, figsize=(16, 4))
    axes[0].imshow(rot, cmap='gray', origin='lower'); axes[0].axis('off')
    axes[0].set_title(f"Specimen at α={ang}°")

    axes[1].plot(proj, np.arange(len(proj)), color='steelblue', lw=1.4)
    axes[1].set_xlabel("Projection"); axes[1].set_ylabel("x")
    axes[1].set_title("1D projection"); axes[1].grid(alpha=0.3)

    axes[2].imshow(ps2, cmap='inferno', origin='lower')
    axes[2].plot(xs_line, ys_line, 'w-', linewidth=1.5, alpha=0.8)
    axes[2].set_title(f"2D FT + central slice\nat {ang}°"); axes[2].axis('off')

    axes[3].plot(fft1d_amp,  label='FT of projection', color='steelblue', lw=1.4)
    axes[3].plot(cslice_amp, label='Central slice of 2D FT', color='red', lw=1.4, ls='--')
    axes[3].set_title("Fourier slice theorem\nverification")
    axes[3].set_xlabel("Frequency index"); axes[3].legend(fontsize=8); axes[3].grid(alpha=0.3)

    plt.tight_layout()
    with out_fst:
        clear_output(wait=True); display(fig2img(fig))

fst_sl.observe(_update_fst, names='value')
display(VBox([fst_sl, out_fst]))
_update_fst()

14.8. Missing wedge#

In practice the tilt range is limited to roughly $\pm 60°$ because:

Grid bars obstruct the beam at high tilt
Increased effective sample thickness at high tilt degrades image quality

This means the central region of Fourier space is never filled — the missing wedge. Features oriented along the tilt axis (perpendicular to the tilt direction) are most affected: they appear elongated in the reconstruction.

The interactive below shows the Fourier coverage at different tilt ranges. The missing wedge is the dark triangular gap.

Show code cell source

Hide code cell source

mt_mw  = IntSlider(value=60, min=10, max=90, step=5,
                    description="Max tilt (°)",       style=_sty, layout=_sl)
di_mw  = IntSlider(value=3,  min=1,  max=10, step=1,
                    description="Tilt increment (°)", style=_sty, layout=_sl)
out_mw = Output()

def _update_mw(_=None):
    n    = im.shape[0]
    angs = np.arange(-mt_mw.value, mt_mw.value+1, di_mw.value)
    sino = radon(im, angs)
    fbp  = filtered_backproject(sino, angs)

    # Coverage map: mark filled lines in Fourier space
    coverage = np.zeros((n, n))
    t = np.linspace(-n//2, n//2, n).astype(int)
    for ang in angs:
        ang_r = np.deg2rad(ang)
        xs = np.clip((t * np.cos(ang_r) + n//2).astype(int), 0, n-1)
        ys = np.clip((t * np.sin(ang_r) + n//2).astype(int), 0, n-1)
        coverage[ys, xs] = 1
    coverage_s = gaussian_filter(coverage, 1.5)

    fig, axes = plt.subplots(1, 4, figsize=(16, 4))
    axes[0].imshow(im,         cmap='gray',    origin='lower'); axes[0].axis('off')
    axes[0].set_title("True image")
    axes[1].imshow(fbp,        cmap='gray',    origin='lower'); axes[1].axis('off')
    axes[1].set_title(f"FBP (±{mt_mw.value}°, Δ={di_mw.value}°)")
    axes[2].imshow(coverage_s, cmap='Blues',   origin='lower'); axes[2].axis('off')
    axes[2].set_title(f"Fourier coverage\n({len(angs)} tilts)")

    # Draw missing-wedge wedge lines on FBP
    mid = n//2
    for sign in [+1, -1]:
        ang_r = np.deg2rad(sign * mt_mw.value + 90)
        t_line = np.linspace(-n//2, n//2, n)
        axes[3].plot(t_line*np.cos(ang_r)+mid, t_line*np.sin(ang_r)+mid,
                     'r--', lw=1.2)
    axes[3].imshow(power_spectrum(fbp), cmap='inferno', origin='lower')
    axes[3].set_title("FBP Fourier space\n(red = missing wedge boundary)")
    axes[3].axis('off')
    plt.suptitle(f"Missing wedge: tilt range ±{mt_mw.value}°", fontsize=10)
    plt.tight_layout()
    with out_mw:
        clear_output(wait=True); display(fig2img(fig))

for w in [mt_mw, di_mw]:
    w.observe(_update_mw, names='value')
display(VBox([mt_mw, di_mw, out_mw]))
_update_mw()

14.9. Task 3 – The Crowther criterion#

To fully reconstruct an object of diameter $D$ at resolution $r$, how many tilt angles $m$ do we need? The answer follows from the Fourier slice theorem and the sampling theorem.

Each projection fills one central line in 2D Fourier space. At resolution $r$, we need the outermost Fourier shell (radius $1/r$) to be densely sampled. Two adjacent central lines at the edge of the Fourier disk must be separated by at most $1/D$ (the Nyquist frequency across the object). The angular separation between adjacent projections at radius $1/r$ is therefore:

\[ \Delta\alpha \leq \frac{1/D}{1/r} = \frac{r}{D} \]

Since the total range is $\pi$ radians (projections at $\alpha$ and $\alpha+\pi$ are identical), the number of tilts needed is:

\[ \boxed{m \geq \frac{\pi}{\Delta\alpha} = \frac{\pi D}{r}} \]

This is the Crowther criterion. It tells us that to achieve resolution $r$ from an object of diameter $D$, we need at least $\pi D/r$ projections — independent of whether we are doing tomography or SPA (in SPA, particles in random orientations naturally fill the full angular range).

Show code cell source

Hide code cell source

# Interactive: show how many tilts are needed for different D and r
_sl_cr = Layout(width="440px")
D_sl = IntSlider(value=10, min=2, max=50, step=1,
                  description="Diameter D (nm)", style=_sty, layout=_sl_cr)
r_sl = FloatSlider(value=1.0, min=0.5, max=5.0, step=0.5,
                    description="Resolution r (nm)", style=_sty, layout=_sl_cr)
out_cr = Output()

def _update_cr(_=None):
    D     = D_sl.value; r = r_sl.value
    m_min = int(np.ceil(np.pi * D / r))
    D_arr = np.linspace(2, 50, 200)
    r_arr = np.array([0.5, 1.0, 2.0, 3.0, 5.0])

    fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
    for ri in r_arr:
        axes[0].plot(D_arr, np.pi * D_arr / ri, lw=1.6, label=f"r={ri} nm")
    axes[0].axvline(D, color='k', ls='--', lw=1.2)
    axes[0].axhline(m_min, color='red', ls=':', lw=1.2,
                    label=f"m={m_min} (D={D} nm, r={r} nm)")
    axes[0].set_xlabel("Object diameter D (nm)"); axes[0].set_ylabel("Min. tilts m")
    axes[0].set_title("Crowther criterion: m ≥ πD/r")
    axes[0].legend(fontsize=8); axes[0].grid(alpha=0.3)

    # Show what m tilts actually cover
    angs_cr = np.linspace(-90, 90, m_min, endpoint=False)
    n_cr    = 80
    cov_cr  = np.zeros((n_cr, n_cr))
    t_idx   = np.arange(n_cr) - n_cr//2
    for ang in angs_cr:
        ang_r = np.deg2rad(ang)
        xs = np.clip((t_idx*np.cos(ang_r)+n_cr//2).astype(int), 0, n_cr-1)
        ys = np.clip((t_idx*np.sin(ang_r)+n_cr//2).astype(int), 0, n_cr-1)
        cov_cr[ys, xs] = 1
    axes[1].imshow(gaussian_filter(cov_cr, 1), cmap='Blues', origin='lower')
    axes[1].set_title(f"Fourier coverage with m={m_min} tilts\n"
                      f"(D={D} nm, r={r} nm → full coverage at radius D/r={D//r:.0f})")
    axes[1].axis('off')
    plt.tight_layout()
    with out_cr:
        clear_output(wait=True); display(fig2img(fig))

for w in [D_sl, r_sl]:
    w.observe(_update_cr, names='value')
display(VBox([D_sl, r_sl, out_cr]))
_update_cr()

14.10. Summary#

Concept	Key idea
Tilt series	Physical rotation of specimen; each tilt gives one 1D projection
Sinogram	Stack of 1D projections; point objects trace sinusoids
Radon transform	Mathematical projection operator $p_\alpha(t) = \int A(t\cos\alpha - s\sin\alpha, t\sin\alpha + s\cos\alpha)\,ds$
Backprojection	Smear each projection back along its direction and sum
Ramp filter	$\hat h(\nu) =
Fourier slice theorem	$\mathcal{F}_1\{p_\alpha\}(\nu) = \mathcal{F}_2\{A\}(\nu\cos\alpha, \nu\sin\alpha)$
Missing wedge	Triangular gap in Fourier coverage from limited tilt range
Crowther criterion	$m \geq \pi D / r$ tilts for resolution $r$ from object of diameter $D$

Practical: Electron Tomography

Contents

14. Practical: Electron Tomography#

14.1. Introduction#

14.1.1. Helper functions#

14.2. The 2D test image#

14.3. Tilting the sample and acquiring projections#

14.4. The sinogram#

14.5. Task 1 – Backprojection#

14.6. Task 2 – Filtered backprojection#

14.7. The Fourier Slice Theorem#

14.8. Missing wedge#

14.9. Task 3 – The Crowther criterion#

14.10. Summary#