Note
Go to the end to download the full example code.
Ab-initio Pipeline Demonstration¶
This tutorial demonstrates some key components of an ab-initio
reconstruction pipeline using synthetic data generated with ASPIRE’s
Simulation class of objects.
Download an Example Volume¶
We begin by downloading a high resolution volume map of the 80S Ribosome, sourced from EMDB: https://www.ebi.ac.uk/emdb/EMD-2660. This is one of several volume maps that can be downloaded with ASPIRE’s data downloading utility by using the following import. sphinx_gallery_start_ignore flake8: noqa sphinx_gallery_end_ignore
from aspire.downloader import emdb_2660
# Load 80s Ribosome as a ``Volume`` object.
original_vol = emdb_2660()
# Downsample the volume
res = 41
vol = original_vol.downsample(res)
Note
A Volume can be saved using the Volume.save() method as follows:
fn = f"downsampled_80s_ribosome_size{res}.mrc"
vol.save(fn, overwrite=True)
Create a Simulation Source¶
ASPIRE’s Simulation class can be used to generate a synthetic
dataset of projection images. A Simulation object produces
random projections of a supplied Volume and applies noise and CTF
filters. The resulting stack of 2D images is stored in an Image
object.
CTF Filters¶
Let’s start by creating CTF filters. The operators package
contains a collection of filter classes that can be supplied to a
Simulation. We use RadialCTFFilter to generate a set of CTF
filters with various defocus values.
# Create CTF filters
import numpy as np
from aspire.operators import RadialCTFFilter
# Radial CTF Filter
defocus_min = 15000 # unit is angstroms
defocus_max = 25000
defocus_ct = 7
ctf_filters = [
RadialCTFFilter(pixel_size=vol.pixel_size, defocus=d)
for d in np.linspace(defocus_min, defocus_max, defocus_ct)
]
Initialize Simulation Object¶
We feed our Volume and filters into Simulation to generate
the dataset of images. When controlled white Gaussian noise is
desired, WhiteNoiseAdder.from_snr() can be used to generate a
simulation data set around a specific SNR.
Alternatively, users can bring their own images using an
ArrayImageSource, or define their own custom noise functions via
Simulation(..., noise_adder=CustomNoiseAdder(...)). Examples
can be found in tutorials/class_averaging.py and
experiments/simulated_abinitio_pipeline.py.
from aspire.noise import WhiteNoiseAdder
from aspire.source import Simulation
# set parameters
n_imgs = 2500
# SNR target for white gaussian noise.
snr = 0.5
Note
Note, the SNR value was chosen based on other parameters for this quick tutorial, and can be changed to adjust the power of the additive noise.
# For this ``Simulation`` we set all 2D offset vectors to zero,
# but by default offset vectors will be randomly distributed.
src = Simulation(
n=n_imgs, # number of projections
vols=vol, # volume source
offsets=0, # Default: images are randomly shifted
unique_filters=ctf_filters,
noise_adder=WhiteNoiseAdder.from_snr(snr=snr), # desired SNR
)
0%| | 0/5 [00:00<?, ?it/s]
20%|██ | 1/5 [00:00<00:00, 4.97it/s]
40%|████ | 2/5 [00:00<00:00, 5.05it/s]
60%|██████ | 3/5 [00:00<00:00, 5.06it/s]
80%|████████ | 4/5 [00:00<00:00, 5.03it/s]
100%|██████████| 5/5 [00:00<00:00, 5.26it/s]
100%|██████████| 5/5 [00:00<00:00, 5.15it/s]
Several Views of the Projection Images¶
We can access several views of the projection images.
# with no corruption applied
src.projections[0:10].show()
# with no noise corruption
src.clean_images[0:10].show()
# with noise and CTF corruption
src.images[0:10].show()
CTF Correction¶
We apply phase_flip() to correct for CTF effects.
src = src.phase_flip()
Cache¶
We apply cache to store the results of the ImageSource
pipeline up to this point. This is optional, but can provide
benefit when used intently on machines with adequate memory.
src = src.cache()
src.images[0:10].show()
0%| | 0/5 [00:00<?, ?it/s]
20%|██ | 1/5 [00:01<00:04, 1.20s/it]
40%|████ | 2/5 [00:02<00:03, 1.20s/it]
60%|██████ | 3/5 [00:03<00:02, 1.22s/it]
80%|████████ | 4/5 [00:04<00:01, 1.21s/it]
100%|██████████| 5/5 [00:05<00:00, 1.16s/it]
100%|██████████| 5/5 [00:05<00:00, 1.18s/it]
Class Averaging¶
We use RIRClass2D object to classify the images via the
rotationally invariant representation (RIR) algorithm. Class
selection is customizable. The classification module also includes a
set of protocols for selecting a set of images to be used after
classification. Here we’re using the simplest
DebugClassAvgSource which internally uses the
TopClassSelector to select the first n_classes images from
the source. In practice, the selection is done by sorting class
averages based on some configurable notion of quality (contrast,
neighbor distance etc).
from aspire.classification import RIRClass2D
# set parameters
n_classes = 200
n_nbor = 6
# We will customize our class averaging source. Note that the
# ``fspca_components`` and ``bispectrum_components`` were selected for
# this small tutorial.
rir = RIRClass2D(
src,
fspca_components=40,
bispectrum_components=30,
n_nbor=n_nbor,
)
from aspire.denoising import DebugClassAvgSource
avgs = DebugClassAvgSource(
src=src,
classifier=rir,
)
# We'll continue our pipeline using only the first ``n_classes`` from
# ``avgs``. The ``cache()`` call is used here to precompute results
# for the ``:n_classes`` slice. This avoids recomputing the same
# images twice when peeking in the next cell then requesting them in
# the following ``CLSyncVoting`` algorithm. Outside of demonstration
# purposes, where we are repeatedly peeking at various stage results,
# such caching can be dropped allowing for more lazy evaluation.
avgs = avgs[:n_classes].cache()
0%| | 0/1 [00:00<?, ?it/s]
0%| | 0/5 [00:00<?, ?it/s]
100%|██████████| 5/5 [00:00<00:00, 820.77it/s]
0%| | 0/5 [00:00<?, ?it/s]
20%|██ | 1/5 [00:00<00:01, 3.55it/s]
40%|████ | 2/5 [00:00<00:00, 3.64it/s]
60%|██████ | 3/5 [00:00<00:00, 3.69it/s]
80%|████████ | 4/5 [00:01<00:00, 3.69it/s]
100%|██████████| 5/5 [00:01<00:00, 3.87it/s]
100%|██████████| 5/5 [00:01<00:00, 3.77it/s]
0%| | 0/200 [00:00<?, ?it/s]
1%| | 2/200 [00:00<00:15, 12.67it/s]
2%|▏ | 4/200 [00:00<00:15, 12.67it/s]
3%|▎ | 6/200 [00:00<00:15, 12.68it/s]
4%|▍ | 8/200 [00:00<00:15, 12.71it/s]
5%|▌ | 10/200 [00:00<00:15, 12.66it/s]
6%|▌ | 12/200 [00:00<00:14, 12.64it/s]
7%|▋ | 14/200 [00:01<00:14, 12.64it/s]
8%|▊ | 16/200 [00:01<00:14, 12.65it/s]
9%|▉ | 18/200 [00:01<00:14, 12.64it/s]
10%|█ | 20/200 [00:01<00:14, 12.56it/s]
11%|█ | 22/200 [00:01<00:14, 12.59it/s]
12%|█▏ | 24/200 [00:01<00:13, 12.63it/s]
13%|█▎ | 26/200 [00:02<00:13, 12.62it/s]
14%|█▍ | 28/200 [00:02<00:13, 12.61it/s]
15%|█▌ | 30/200 [00:02<00:13, 12.64it/s]
16%|█▌ | 32/200 [00:02<00:13, 12.60it/s]
17%|█▋ | 34/200 [00:02<00:13, 12.59it/s]
18%|█▊ | 36/200 [00:02<00:13, 12.60it/s]
19%|█▉ | 38/200 [00:03<00:12, 12.58it/s]
20%|██ | 40/200 [00:03<00:13, 12.26it/s]
21%|██ | 42/200 [00:03<00:12, 12.36it/s]
22%|██▏ | 44/200 [00:03<00:12, 12.39it/s]
23%|██▎ | 46/200 [00:03<00:12, 12.46it/s]
24%|██▍ | 48/200 [00:03<00:12, 12.54it/s]
25%|██▌ | 50/200 [00:03<00:11, 12.61it/s]
26%|██▌ | 52/200 [00:04<00:11, 12.63it/s]
27%|██▋ | 54/200 [00:04<00:11, 12.65it/s]
28%|██▊ | 56/200 [00:04<00:11, 12.62it/s]
29%|██▉ | 58/200 [00:04<00:11, 12.58it/s]
30%|███ | 60/200 [00:04<00:11, 12.22it/s]
31%|███ | 62/200 [00:04<00:11, 12.33it/s]
32%|███▏ | 64/200 [00:05<00:10, 12.36it/s]
33%|███▎ | 66/200 [00:05<00:10, 12.40it/s]
34%|███▍ | 68/200 [00:05<00:10, 12.42it/s]
35%|███▌ | 70/200 [00:05<00:10, 12.42it/s]
36%|███▌ | 72/200 [00:05<00:10, 12.35it/s]
37%|███▋ | 74/200 [00:05<00:10, 11.99it/s]
38%|███▊ | 76/200 [00:06<00:10, 12.13it/s]
39%|███▉ | 78/200 [00:06<00:10, 12.16it/s]
40%|████ | 80/200 [00:06<00:09, 12.26it/s]
41%|████ | 82/200 [00:06<00:09, 12.18it/s]
42%|████▏ | 84/200 [00:06<00:09, 12.34it/s]
43%|████▎ | 86/200 [00:06<00:09, 12.28it/s]
44%|████▍ | 88/200 [00:07<00:09, 12.41it/s]
45%|████▌ | 90/200 [00:07<00:08, 12.47it/s]
46%|████▌ | 92/200 [00:07<00:08, 12.50it/s]
47%|████▋ | 94/200 [00:07<00:08, 12.56it/s]
48%|████▊ | 96/200 [00:07<00:08, 12.56it/s]
49%|████▉ | 98/200 [00:07<00:08, 12.59it/s]
50%|█████ | 100/200 [00:08<00:07, 12.59it/s]
51%|█████ | 102/200 [00:08<00:07, 12.61it/s]
52%|█████▏ | 104/200 [00:08<00:07, 12.64it/s]
53%|█████▎ | 106/200 [00:08<00:07, 12.63it/s]
54%|█████▍ | 108/200 [00:08<00:07, 12.63it/s]
55%|█████▌ | 110/200 [00:08<00:07, 12.68it/s]
56%|█████▌ | 112/200 [00:08<00:06, 12.65it/s]
57%|█████▋ | 114/200 [00:09<00:06, 12.58it/s]
58%|█████▊ | 116/200 [00:09<00:06, 12.50it/s]
59%|█████▉ | 118/200 [00:09<00:06, 12.44it/s]
60%|██████ | 120/200 [00:09<00:06, 12.44it/s]
61%|██████ | 122/200 [00:09<00:06, 12.49it/s]
62%|██████▏ | 124/200 [00:09<00:06, 12.44it/s]
63%|██████▎ | 126/200 [00:10<00:05, 12.48it/s]
64%|██████▍ | 128/200 [00:10<00:05, 12.50it/s]
65%|██████▌ | 130/200 [00:10<00:05, 12.51it/s]
66%|██████▌ | 132/200 [00:10<00:05, 12.52it/s]
67%|██████▋ | 134/200 [00:10<00:05, 12.53it/s]
68%|██████▊ | 136/200 [00:10<00:05, 12.41it/s]
69%|██████▉ | 138/200 [00:11<00:04, 12.47it/s]
70%|███████ | 140/200 [00:11<00:04, 12.51it/s]
71%|███████ | 142/200 [00:11<00:04, 12.52it/s]
72%|███████▏ | 144/200 [00:11<00:04, 12.53it/s]
73%|███████▎ | 146/200 [00:11<00:04, 12.52it/s]
74%|███████▍ | 148/200 [00:11<00:04, 12.51it/s]
75%|███████▌ | 150/200 [00:12<00:03, 12.52it/s]
76%|███████▌ | 152/200 [00:12<00:03, 12.53it/s]
77%|███████▋ | 154/200 [00:12<00:03, 12.56it/s]
78%|███████▊ | 156/200 [00:12<00:03, 12.57it/s]
79%|███████▉ | 158/200 [00:12<00:03, 12.56it/s]
80%|████████ | 160/200 [00:12<00:03, 12.37it/s]
81%|████████ | 162/200 [00:12<00:03, 12.41it/s]
82%|████████▏ | 164/200 [00:13<00:02, 12.42it/s]
83%|████████▎ | 166/200 [00:13<00:02, 12.44it/s]
84%|████████▍ | 168/200 [00:13<00:02, 12.49it/s]
85%|████████▌ | 170/200 [00:13<00:02, 12.49it/s]
86%|████████▌ | 172/200 [00:13<00:02, 12.47it/s]
87%|████████▋ | 174/200 [00:13<00:02, 12.47it/s]
88%|████████▊ | 176/200 [00:14<00:01, 12.41it/s]
89%|████████▉ | 178/200 [00:14<00:01, 12.44it/s]
90%|█████████ | 180/200 [00:14<00:01, 12.48it/s]
91%|█████████ | 182/200 [00:14<00:01, 12.51it/s]
92%|█████████▏| 184/200 [00:14<00:01, 12.55it/s]
93%|█████████▎| 186/200 [00:14<00:01, 12.40it/s]
94%|█████████▍| 188/200 [00:15<00:00, 12.42it/s]
95%|█████████▌| 190/200 [00:15<00:00, 12.45it/s]
96%|█████████▌| 192/200 [00:15<00:00, 12.50it/s]
97%|█████████▋| 194/200 [00:15<00:00, 12.49it/s]
98%|█████████▊| 196/200 [00:15<00:00, 12.48it/s]
99%|█████████▉| 198/200 [00:15<00:00, 12.44it/s]
100%|██████████| 200/200 [00:16<00:00, 12.47it/s]
100%|██████████| 200/200 [00:16<00:00, 12.49it/s]
0%| | 0/200 [00:00<?, ?it/s]
10%|▉ | 19/200 [00:00<00:00, 185.98it/s]
19%|█▉ | 38/200 [00:00<00:00, 185.44it/s]
28%|██▊ | 57/200 [00:00<00:00, 186.73it/s]
38%|███▊ | 76/200 [00:00<00:00, 186.06it/s]
48%|████▊ | 95/200 [00:00<00:00, 185.74it/s]
57%|█████▋ | 114/200 [00:00<00:00, 186.50it/s]
66%|██████▋ | 133/200 [00:00<00:00, 186.93it/s]
76%|███████▌ | 152/200 [00:00<00:00, 182.80it/s]
86%|████████▌ | 172/200 [00:00<00:00, 185.15it/s]
96%|█████████▌| 191/200 [00:01<00:00, 185.72it/s]
100%|██████████| 200/200 [00:01<00:00, 185.45it/s]
100%|██████████| 1/1 [00:24<00:00, 24.05s/it]
100%|██████████| 1/1 [00:24<00:00, 24.05s/it]
View the Class Averages¶
# Show class averages
avgs.images[0:10].show()
# Show original images corresponding to those classes. This 1:1
# comparison is only expected to work because we used
# ``TopClassSelector`` to classify our images.
src.images[0:10].show()
Orientation Estimation¶
We create an OrientedSource, which consumes an ImageSource object, an
orientation estimator, and returns a new source which lazily estimates orientations.
In this case we supply avgs for our source and a CLSyncVoting
class instance for our orientation estimator. The CLSyncVoting algorithm employs
a common-lines method with synchronization and voting.
from aspire.abinitio import CLSyncVoting
from aspire.source import OrientedSource
# Stash true rotations for later comparison
true_rotations = src.rotations[:n_classes]
# For this low resolution example we will customize the ``CLSyncVoting``
# instance to use fewer theta points ``n_theta`` then the default value of 360.
orient_est = CLSyncVoting(avgs, n_theta=72)
# Instantiate an ``OrientedSource``.
oriented_src = OrientedSource(avgs, orient_est)
Mean Error of Estimated Rotations¶
ASPIRE has the built-in utility function, mean_aligned_angular_distance, which globally
aligns the estimated rotations to the true rotations and computes the mean
angular distance (in degrees).
from aspire.utils import mean_aligned_angular_distance
# Compare with known true rotations
mean_ang_dist = mean_aligned_angular_distance(oriented_src.rotations, true_rotations)
print(f"Mean aligned angular distance: {mean_ang_dist} degrees")
Searching over common line pairs: 0%| | 0/19900 [00:00<?, ?it/s]
Searching over common line pairs: 1%|▏ | 275/19900 [00:00<00:07, 2747.78it/s]
Searching over common line pairs: 3%|▎ | 568/19900 [00:00<00:06, 2853.69it/s]
Searching over common line pairs: 4%|▍ | 859/19900 [00:00<00:06, 2876.28it/s]
Searching over common line pairs: 6%|▌ | 1151/19900 [00:00<00:06, 2890.22it/s]
Searching over common line pairs: 7%|▋ | 1442/19900 [00:00<00:06, 2894.81it/s]
Searching over common line pairs: 9%|▊ | 1735/19900 [00:00<00:06, 2903.43it/s]
Searching over common line pairs: 10%|█ | 2027/19900 [00:00<00:06, 2908.66it/s]
Searching over common line pairs: 12%|█▏ | 2319/19900 [00:00<00:06, 2911.99it/s]
Searching over common line pairs: 13%|█▎ | 2611/19900 [00:00<00:05, 2912.90it/s]
Searching over common line pairs: 15%|█▍ | 2903/19900 [00:01<00:05, 2911.09it/s]
Searching over common line pairs: 16%|█▌ | 3195/19900 [00:01<00:05, 2911.27it/s]
Searching over common line pairs: 18%|█▊ | 3487/19900 [00:01<00:05, 2911.54it/s]
Searching over common line pairs: 19%|█▉ | 3779/19900 [00:01<00:05, 2913.67it/s]
Searching over common line pairs: 20%|██ | 4071/19900 [00:01<00:05, 2913.38it/s]
Searching over common line pairs: 22%|██▏ | 4363/19900 [00:01<00:05, 2910.55it/s]
Searching over common line pairs: 23%|██▎ | 4656/19900 [00:01<00:05, 2913.42it/s]
Searching over common line pairs: 25%|██▍ | 4948/19900 [00:01<00:05, 2907.12it/s]
Searching over common line pairs: 26%|██▋ | 5239/19900 [00:01<00:05, 2903.72it/s]
Searching over common line pairs: 28%|██▊ | 5530/19900 [00:01<00:04, 2902.09it/s]
Searching over common line pairs: 29%|██▉ | 5821/19900 [00:02<00:04, 2900.56it/s]
Searching over common line pairs: 31%|███ | 6112/19900 [00:02<00:04, 2889.52it/s]
Searching over common line pairs: 32%|███▏ | 6401/19900 [00:02<00:04, 2877.12it/s]
Searching over common line pairs: 34%|███▎ | 6692/19900 [00:02<00:04, 2886.62it/s]
Searching over common line pairs: 35%|███▌ | 6983/19900 [00:02<00:04, 2891.00it/s]
Searching over common line pairs: 37%|███▋ | 7273/19900 [00:02<00:04, 2891.71it/s]
Searching over common line pairs: 38%|███▊ | 7564/19900 [00:02<00:04, 2896.30it/s]
Searching over common line pairs: 39%|███▉ | 7856/19900 [00:02<00:04, 2901.97it/s]
Searching over common line pairs: 41%|████ | 8147/19900 [00:02<00:04, 2904.02it/s]
Searching over common line pairs: 42%|████▏ | 8439/19900 [00:02<00:03, 2908.08it/s]
Searching over common line pairs: 44%|████▍ | 8731/19900 [00:03<00:03, 2908.87it/s]
Searching over common line pairs: 45%|████▌ | 9023/19900 [00:03<00:03, 2910.33it/s]
Searching over common line pairs: 47%|████▋ | 9315/19900 [00:03<00:03, 2912.92it/s]
Searching over common line pairs: 48%|████▊ | 9607/19900 [00:03<00:03, 2912.84it/s]
Searching over common line pairs: 50%|████▉ | 9899/19900 [00:03<00:03, 2911.53it/s]
Searching over common line pairs: 51%|█████ | 10191/19900 [00:03<00:03, 2907.37it/s]
Searching over common line pairs: 53%|█████▎ | 10482/19900 [00:03<00:03, 2901.82it/s]
Searching over common line pairs: 54%|█████▍ | 10773/19900 [00:03<00:03, 2903.29it/s]
Searching over common line pairs: 56%|█████▌ | 11065/19900 [00:03<00:03, 2905.93it/s]
Searching over common line pairs: 57%|█████▋ | 11357/19900 [00:03<00:02, 2909.15it/s]
Searching over common line pairs: 59%|█████▊ | 11648/19900 [00:04<00:02, 2906.40it/s]
Searching over common line pairs: 60%|█████▉ | 11939/19900 [00:04<00:02, 2907.43it/s]
Searching over common line pairs: 61%|██████▏ | 12230/19900 [00:04<00:02, 2906.33it/s]
Searching over common line pairs: 63%|██████▎ | 12521/19900 [00:04<00:02, 2905.44it/s]
Searching over common line pairs: 64%|██████▍ | 12813/19900 [00:04<00:02, 2907.98it/s]
Searching over common line pairs: 66%|██████▌ | 13104/19900 [00:04<00:02, 2906.86it/s]
Searching over common line pairs: 67%|██████▋ | 13395/19900 [00:04<00:02, 2905.81it/s]
Searching over common line pairs: 69%|██████▉ | 13687/19900 [00:04<00:02, 2907.39it/s]
Searching over common line pairs: 70%|███████ | 13979/19900 [00:04<00:02, 2909.41it/s]
Searching over common line pairs: 72%|███████▏ | 14270/19900 [00:04<00:01, 2908.65it/s]
Searching over common line pairs: 73%|███████▎ | 14562/19900 [00:05<00:01, 2910.88it/s]
Searching over common line pairs: 75%|███████▍ | 14854/19900 [00:05<00:01, 2912.69it/s]
Searching over common line pairs: 76%|███████▌ | 15147/19900 [00:05<00:01, 2915.40it/s]
Searching over common line pairs: 78%|███████▊ | 15439/19900 [00:05<00:01, 2914.73it/s]
Searching over common line pairs: 79%|███████▉ | 15731/19900 [00:05<00:01, 2911.20it/s]
Searching over common line pairs: 81%|████████ | 16023/19900 [00:05<00:01, 2907.83it/s]
Searching over common line pairs: 82%|████████▏ | 16314/19900 [00:05<00:01, 2898.03it/s]
Searching over common line pairs: 83%|████████▎ | 16606/19900 [00:05<00:01, 2903.27it/s]
Searching over common line pairs: 85%|████████▍ | 16897/19900 [00:05<00:01, 2894.34it/s]
Searching over common line pairs: 86%|████████▋ | 17190/19900 [00:05<00:00, 2902.68it/s]
Searching over common line pairs: 88%|████████▊ | 17482/19900 [00:06<00:00, 2905.55it/s]
Searching over common line pairs: 89%|████████▉ | 17774/19900 [00:06<00:00, 2909.11it/s]
Searching over common line pairs: 91%|█████████ | 18066/19900 [00:06<00:00, 2911.10it/s]
Searching over common line pairs: 92%|█████████▏| 18358/19900 [00:06<00:00, 2911.26it/s]
Searching over common line pairs: 94%|█████████▎| 18650/19900 [00:06<00:00, 2910.51it/s]
Searching over common line pairs: 95%|█████████▌| 18942/19900 [00:06<00:00, 2908.43it/s]
Searching over common line pairs: 97%|█████████▋| 19234/19900 [00:06<00:00, 2909.88it/s]
Searching over common line pairs: 98%|█████████▊| 19525/19900 [00:06<00:00, 2908.74it/s]
Searching over common line pairs: 100%|█████████▉| 19816/19900 [00:06<00:00, 2907.44it/s]
Searching over common line pairs: 100%|██████████| 19900/19900 [00:06<00:00, 2904.21it/s]
LSQR Least-squares solution of Ax = b
The matrix A has 19900 rows and 400 columns
damp = 0.00000000000000e+00 calc_var = 0
atol = 1.00e-06 conlim = 1.00e+08
btol = 1.00e-06 iter_lim = 800
Itn x[0] r1norm r2norm Compatible LS Norm A Cond A
0 0.00000e+00 7.900e+01 7.900e+01 1.0e+00 4.1e-02
1 -9.09180e-02 7.372e+01 7.372e+01 9.3e-01 7.6e-02 9.1e+00 1.0e+00
2 -6.20053e-02 7.347e+01 7.347e+01 9.3e-01 1.5e-02 1.3e+01 2.0e+00
3 -7.09242e-02 7.345e+01 7.345e+01 9.3e-01 4.6e-03 1.5e+01 3.1e+00
4 -6.65593e-02 7.345e+01 7.345e+01 9.3e-01 2.2e-03 1.7e+01 4.3e+00
5 -6.65370e-02 7.345e+01 7.345e+01 9.3e-01 9.2e-04 1.9e+01 5.7e+00
6 -6.57091e-02 7.345e+01 7.345e+01 9.3e-01 3.0e-04 2.1e+01 6.9e+00
7 -6.61623e-02 7.345e+01 7.345e+01 9.3e-01 1.1e-04 2.2e+01 8.1e+00
8 -6.59082e-02 7.345e+01 7.345e+01 9.3e-01 4.6e-05 2.3e+01 9.3e+00
9 -6.60291e-02 7.345e+01 7.345e+01 9.3e-01 1.9e-05 2.5e+01 1.0e+01
10 -6.59594e-02 7.345e+01 7.345e+01 9.3e-01 7.7e-06 2.6e+01 1.2e+01
11 -6.59790e-02 7.345e+01 7.345e+01 9.3e-01 3.5e-06 2.7e+01 1.3e+01
12 -6.59683e-02 7.345e+01 7.345e+01 9.3e-01 1.3e-06 2.8e+01 1.4e+01
13 -6.59708e-02 7.345e+01 7.345e+01 9.3e-01 4.4e-07 3.0e+01 1.5e+01
LSQR finished
The least-squares solution is good enough, given atol
istop = 2 r1norm = 7.3e+01 anorm = 3.0e+01 arnorm = 9.6e-04
itn = 13 r2norm = 7.3e+01 acond = 1.5e+01 xnorm = 3.4e+00
Mean aligned angular distance: 13.057314471222043 degrees
Volume Reconstruction¶
Now that we have our class averages and rotation estimates, we can estimate the mean volume by supplying the class averages and basis for back projection.
from aspire.reconstruction import MeanEstimator
# Setup an estimator to perform the back projection.
estimator = MeanEstimator(oriented_src)
# Perform the estimation and save the volume.
estimated_volume = estimator.estimate()
Comparison of Estimated Volume with Source Volume¶
To get a visual confirmation that our results are sane, we rotate the estimated volume by the estimated rotations and project along the z-axis. These estimated projections should align with the original projection images.
# Get the first 10 projections from the estimated volume using the
# estimated orientations. Recall that ``project`` returns an
# ``Image`` instance, which we can peek at with ``show``.
projections_est = estimated_volume.project(oriented_src.rotations[0:10])
# We view the first 10 projections of the estimated volume.
projections_est.show()
# For comparison, we view the first 10 source projections.
src.projections[0:10].show()
Total running time of the script: (1 minutes 6.861 seconds)