3D Image Orientation

This script illustrates the estimation of orientation angles using a synchronization matrix and the voting method, based on simulated data projected from a 3D cryo-EM map.

import logging
import os

import numpy as np

from aspire.abinitio import CLSyncVoting
from aspire.operators import RadialCTFFilter
from aspire.source import OrientedSource, Simulation
from aspire.utils import mean_aligned_angular_distance
from aspire.volume import Volume

logger = logging.getLogger(__name__)

file_path = os.path.join(
    os.path.dirname(os.getcwd()), "data", "clean70SRibosome_vol_65p.mrc"
)

logger.info(
    "This script illustrates orientation estimation using "
    "synchronization matrix and voting method"
)

Initialize Simulation Object and CTF Filters

# Define a precision for this experiment
dtype = np.float32

# Set the sizes of images
img_size = 33

# Set the total number of images generated from the 3D map
num_imgs = 128

# Specify the CTF parameters not used for this example
# but necessary for initializing the simulation object
pixel_size = 5  # Pixel size of the images (in angstroms)
voltage = 200  # Voltage (in KV)
defocus_min = 1.5e4  # Minimum defocus value (in angstroms)
defocus_max = 2.5e4  # Maximum defocus value (in angstroms)
defocus_ct = 7  # Number of defocus groups.
Cs = 2.0  # Spherical aberration
alpha = 0.1  # Amplitude contrast

logger.info("Initialize simulation object and CTF filters.")
# Create CTF filters
filters = [
    RadialCTFFilter(pixel_size, voltage, defocus=d, Cs=2.0, alpha=0.1)
    for d in np.linspace(defocus_min, defocus_max, defocus_ct)
]

Downsampling

# Load the map file of a 70S Ribosome and downsample the 3D map to desired resolution.
# The downsampling can be done by the internal function of Volume object.
logger.info(
    f"Load 3D map and downsample 3D map to desired grids "
    f"of {img_size} x {img_size} x {img_size}."
)
vols = Volume.load(file_path, dtype=dtype)
vols = vols.downsample(img_size)

Create Simulation Object and Obtain True Rotation Angles

# Create a simulation object with specified filters and the downsampled 3D map
logger.info("Use downsampled map to creat simulation object.")
sim = Simulation(L=img_size, n=num_imgs, vols=vols, unique_filters=filters, dtype=dtype)

logger.info("Get true rotation angles generated randomly by the simulation object.")
rots_true = sim.rotations

Estimate Orientation and Rotation Angles

# Initialize an orientation estimation object and create an ``OrientedSource`` object
# to perform viewing angle estimation. Here, because of the small image size of the
# ``Simulation``, we customize the ``CLSyncVoting`` method to use fewer theta values
# when searching for common-lines between pairs of images. Additionally, since we are
# processing images with no noise, we opt not to use a ``fuzzy_mask``, an option that
# improves common-line detection in higher noise regimes.
logger.info("Estimate rotation angles using synchronization matrix and voting method.")
orient_est = CLSyncVoting(sim, n_theta=36, mask=False)
oriented_src = OrientedSource(sim, orient_est)
rots_est = oriented_src.rotations

Searching over common line pairs:   0%|          | 0/8128 [00:00<?, ?it/s]
Searching over common line pairs:   5%|▌         | 445/8128 [00:00<00:01, 4449.73it/s]
Searching over common line pairs:  11%|█         | 904/8128 [00:00<00:01, 4531.28it/s]
Searching over common line pairs:  17%|█▋        | 1365/8128 [00:00<00:01, 4566.62it/s]
Searching over common line pairs:  22%|██▏       | 1824/8128 [00:00<00:01, 4575.30it/s]
Searching over common line pairs:  28%|██▊       | 2286/8128 [00:00<00:01, 4588.94it/s]
Searching over common line pairs:  34%|███▍      | 2746/8128 [00:00<00:01, 4592.05it/s]
Searching over common line pairs:  39%|███▉      | 3206/8128 [00:00<00:01, 4591.07it/s]
Searching over common line pairs:  45%|████▌     | 3666/8128 [00:00<00:00, 4593.12it/s]
Searching over common line pairs:  51%|█████     | 4128/8128 [00:00<00:00, 4599.21it/s]
Searching over common line pairs:  56%|█████▋    | 4588/8128 [00:01<00:00, 4592.96it/s]
Searching over common line pairs:  62%|██████▏   | 5048/8128 [00:01<00:00, 4588.03it/s]
Searching over common line pairs:  68%|██████▊   | 5507/8128 [00:01<00:00, 4542.85it/s]
Searching over common line pairs:  73%|███████▎  | 5962/8128 [00:01<00:00, 4509.37it/s]
Searching over common line pairs:  79%|███████▉  | 6418/8128 [00:01<00:00, 4521.58it/s]
Searching over common line pairs:  85%|████████▍ | 6874/8128 [00:01<00:00, 4532.94it/s]
Searching over common line pairs:  90%|█████████ | 7332/8128 [00:01<00:00, 4546.13it/s]
Searching over common line pairs:  96%|█████████▌| 7788/8128 [00:01<00:00, 4549.08it/s]
Searching over common line pairs: 100%|██████████| 8128/8128 [00:01<00:00, 4557.74it/s]

Mean Angular Distance

# ``mean_aligned_angular_distance`` will perform global alignment of the estimated rotations
# to the ground truth and find the mean angular distance between them (in degrees).
mean_ang_dist = mean_aligned_angular_distance(rots_est, rots_true)
logger.info(
    f"Mean angular distance between estimates and ground truth: {mean_ang_dist} degrees"
)

# Basic Check
assert mean_ang_dist < 10