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 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

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

print(
    "This script illustrates orientation estimation using "
    "synchronization matrix and voting method"
)

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

print("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)
]

Initialize simulation object and CTF filters.

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.
print(
    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)

Load 3D map and downsample 3D map to desired grids of 33 x 33 x 33.

Create Simulation Object and Obtain True Rotation Angles

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

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

Use downsampled map to creat simulation object.
Get true rotation angles generated randomly by the simulation object.

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.
print(
    "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

Estimate rotation angles and offsets using synchronization matrix and voting method.

Searching over common line pairs:   0%|          | 0/8128 [00:00<?, ?it/s]
Searching over common line pairs:   7%|▋         | 584/8128 [00:00<00:01, 5833.47it/s]
Searching over common line pairs:  14%|█▍        | 1169/8128 [00:00<00:01, 5838.38it/s]
Searching over common line pairs:  22%|██▏       | 1753/8128 [00:00<00:01, 5821.78it/s]
Searching over common line pairs:  29%|██▊       | 2336/8128 [00:00<00:01, 5791.74it/s]
Searching over common line pairs:  36%|███▌      | 2918/8128 [00:00<00:00, 5799.01it/s]
Searching over common line pairs:  43%|████▎     | 3499/8128 [00:00<00:00, 5801.58it/s]
Searching over common line pairs:  50%|█████     | 4081/8128 [00:00<00:00, 5806.75it/s]
Searching over common line pairs:  57%|█████▋    | 4664/8128 [00:00<00:00, 5811.55it/s]
Searching over common line pairs:  65%|██████▍   | 5247/8128 [00:00<00:00, 5816.71it/s]
Searching over common line pairs:  72%|███████▏  | 5829/8128 [00:01<00:00, 5816.03it/s]
Searching over common line pairs:  79%|███████▉  | 6411/8128 [00:01<00:00, 5812.40it/s]
Searching over common line pairs:  86%|████████▌ | 6993/8128 [00:01<00:00, 5795.66it/s]
Searching over common line pairs:  93%|█████████▎| 7573/8128 [00:01<00:00, 5780.72it/s]
Searching over common line pairs: 100%|██████████| 8128/8128 [00:01<00:00, 5797.86it/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.60901e-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.33993e-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)
print(
    f"Mean angular distance between estimates and ground truth: {mean_ang_dist} degrees"
)

# Basic Check
assert mean_ang_dist < 10

Mean angular distance between estimates and ground truth: 8.615618705749512 degrees

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()
print(f"Mean offset error in pixels {offs_err}, approx {offs_err/img_size*100:.1f}%")

Mean offset error in pixels 0.724621057510376, approx 2.2%