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

# Initialize an orientation estimation object and create an
# ``OrientedSource`` object to perform viewing angle and image offset
# 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 and offsets 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:   6%|▌         | 493/8128 [00:00<00:01, 4921.53it/s]
Searching over common line pairs:  12%|█▏        | 998/8128 [00:00<00:01, 4994.07it/s]
Searching over common line pairs:  18%|█▊        | 1498/8128 [00:00<00:01, 4987.40it/s]
Searching over common line pairs:  25%|██▍       | 1997/8128 [00:00<00:01, 4670.65it/s]
Searching over common line pairs:  31%|███       | 2502/8128 [00:00<00:01, 4801.39it/s]
Searching over common line pairs:  37%|███▋      | 3007/8128 [00:00<00:01, 4880.63it/s]
Searching over common line pairs:  43%|████▎     | 3500/8128 [00:00<00:00, 4893.97it/s]
Searching over common line pairs:  49%|████▉     | 4008/8128 [00:00<00:00, 4950.09it/s]
Searching over common line pairs:  56%|█████▌    | 4513/8128 [00:00<00:00, 4980.42it/s]
Searching over common line pairs:  62%|██████▏   | 5017/8128 [00:01<00:00, 4995.81it/s]
Searching over common line pairs:  68%|██████▊   | 5521/8128 [00:01<00:00, 5008.59it/s]
Searching over common line pairs:  74%|███████▍  | 6023/8128 [00:01<00:00, 5006.12it/s]
Searching over common line pairs:  80%|████████  | 6525/8128 [00:01<00:00, 5008.88it/s]
Searching over common line pairs:  86%|████████▋ | 7027/8128 [00:01<00:00, 5006.69it/s]
Searching over common line pairs:  93%|█████████▎| 7528/8128 [00:01<00:00, 5000.98it/s]
Searching over common line pairs:  99%|█████████▉| 8029/8128 [00:01<00:00, 4993.66it/s]
Searching over common line pairs: 100%|██████████| 8128/8128 [00:01<00:00, 4950.96it/s]

LSQR            Least-squares solution of  Ax = b
The matrix A has 8128 rows and 256 columns
damp = 0.00000000000000e+00   calc_var =        0
atol = 1.00e-06                 conlim = 1.00e+08
btol = 1.00e-06               iter_lim =      512

   Itn      x[0]       r1norm     r2norm   Compatible    LS      Norm A   Cond A
     0  0.00000e+00   2.322e+02  2.322e+02    1.0e+00  2.2e-02
     1  1.06962e+00   1.656e+02  1.656e+02    7.1e-01  2.1e-01   7.4e+00  1.0e+00
     2  8.60900e-01   1.614e+02  1.614e+02    7.0e-01  4.6e-02   1.0e+01  2.1e+00
     3  8.97186e-01   1.609e+02  1.609e+02    6.9e-01  1.4e-02   1.2e+01  3.2e+00
     4  9.52016e-01   1.609e+02  1.609e+02    6.9e-01  5.8e-03   1.4e+01  4.4e+00
     5  9.33992e-01   1.609e+02  1.609e+02    6.9e-01  2.2e-03   1.5e+01  5.6e+00
     6  9.37907e-01   1.609e+02  1.609e+02    6.9e-01  9.5e-04   1.7e+01  6.8e+00
     7  9.34790e-01   1.609e+02  1.609e+02    6.9e-01  4.5e-04   1.8e+01  8.1e+00
     8  9.32634e-01   1.609e+02  1.609e+02    6.9e-01  1.9e-04   1.9e+01  9.4e+00
     9  9.33399e-01   1.609e+02  1.609e+02    6.9e-01  8.3e-05   2.0e+01  1.1e+01
    10  9.32935e-01   1.609e+02  1.609e+02    6.9e-01  3.8e-05   2.1e+01  1.2e+01
    12  9.32869e-01   1.609e+02  1.609e+02    6.9e-01  8.4e-06   2.3e+01  1.4e+01
    13  9.32768e-01   1.609e+02  1.609e+02    6.9e-01  2.8e-06   2.4e+01  1.6e+01
    14  9.32769e-01   1.609e+02  1.609e+02    6.9e-01  8.1e-07   2.5e+01  1.7e+01

LSQR finished
The least-squares solution is good enough, given atol

istop =       2   r1norm = 1.6e+02   anorm = 2.5e+01   arnorm = 3.3e-03
itn   =      14   r2norm = 1.6e+02   acond = 1.7e+01   xnorm  = 2.4e+01

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

Offsets Estimation

# The ground truth offsets from the simulation can be compared to
# those estimated by the commonlines method above.

# Calculate Estimation error in pixels for each image.
offs_diff = np.sqrt(np.sum((oriented_src.offsets - sim.offsets) ** 2, axis=1))

# Calculate the mean error in pixels across all images.
offs_err = offs_diff.mean()
logger.info(
    f"Mean offset error in pixels {offs_err}, approx {offs_err/img_size*100:.1f}%"
)