Note
Go to the end to download the full example code.
Image Preprocessing¶
This script illustrates the preprocess steps implemented prior to starting the pipeline of reconstructing a 3D map using simulated 2D images.
import logging
import os
import matplotlib.pyplot as plt
import numpy as np
from aspire.noise import WhiteNoiseAdder
from aspire.operators import RadialCTFFilter
from aspire.source.simulation import Simulation
from aspire.volume import Volume
logger = logging.getLogger(__name__)
file_path = os.path.join(
os.path.dirname(os.getcwd()), "data", "clean70SRibosome_vol_65p.mrc"
)
Specify Parameters¶
# Set the initial simulated full size of images.
img_size = 33
# Set the demonstration downsampled image size.
ds_img_size = 15
# Set the total number of images generated from the 3D map
num_imgs = 512
# Set the noise variance and build the noise filter
noise_variance = 4e-1
noise_adder = WhiteNoiseAdder(var=noise_variance)
# Specify the CTF parameters not used for this example
# but necessary for initializing the simulation object
pixel_size = 5 * 65 / img_size # 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
Build Simulation Object and Apply Noise¶
logger.info("Initialize simulation object and CTF filters.")
# Create CTF filters
ctf_filters = [
RadialCTFFilter(pixel_size, voltage, defocus=d, Cs=2.0, alpha=0.1)
for d in np.linspace(defocus_min, defocus_max, defocus_ct)
]
# Load the map file of a 70S ribosome and downsample the 3D map to desired image size.
logger.info("Load 3D map from mrc file")
vols = Volume.load(file_path)
# Downsample the volume to a desired image size and increase density
# by 1.0e5 time for a better graph view
logger.info(f"Downsample map to a image size of {img_size} x {img_size} x {img_size}")
vols = vols.downsample(img_size) * 1.0e5
# Create a simulation object with specified filters and the downsampled 3D map
logger.info("Use downsampled map to create simulation object.")
source = Simulation(
L=img_size,
n=num_imgs,
vols=vols,
unique_filters=ctf_filters,
noise_adder=noise_adder,
)
Apply Independent Preprocessing Techniques¶
Now we’ll apply each technique sequentially. This is easily
accomplished because each preprocessing technique returns a new
ImageSource
object. In this case we assign each to a new
variable source_*
. That leaves the original source
object
untouched.
logger.info("Obtain original images.")
imgs_od = source.images[0].asnumpy()
logger.info("Perform phase flip to input images.")
source_pf = source.phase_flip()
logger.info(f"Downsample image size to {ds_img_size} X {ds_img_size}")
source_ds = source.downsample(ds_img_size)
logger.info("Normalize images to background noise.")
source_nb = source.normalize_background()
logger.info("Whiten noise of images")
source_wt = source.whiten()
logger.info("Invert the global density contrast if need")
source_rc = source.invert_contrast()
0%| | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:00<00:00, 1.02it/s]
100%|██████████| 1/1 [00:00<00:00, 1.02it/s]
0%| | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:00<00:00, 1.02it/s]
100%|██████████| 1/1 [00:00<00:00, 1.02it/s]
Plot First Image from Each Preprocess Step¶
Plot the first images.
logger.info("Plot first image from each preprocess steps")
idm = 0
plt.subplot(2, 3, 1)
plt.imshow(imgs_od[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("original image")
plt.subplot(2, 3, 2)
plt.imshow(source_pf.images[0].asnumpy()[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("phase flip")
plt.subplot(2, 3, 3)
plt.imshow(source_ds.images[0].asnumpy()[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("downsample")
plt.subplot(2, 3, 4)
plt.imshow(source_nb.images[0].asnumpy()[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("normalize background")
plt.subplot(2, 3, 5)
plt.imshow(source_wt.images[0].asnumpy()[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("noise whitening")
plt.subplot(2, 3, 6)
plt.imshow(source_rc.images[0].asnumpy()[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("invert contrast")
plt.tight_layout()
Apply Sequential Preprocessing Techniques¶
Now we’ll apply each technique sequentially.
This is accomplished by reassigning each new ImageSource
to the same variable.
In this case we reassign to source
.
Note, after each source
assignment we’ll manually save off the images for the plot below.
# We'll copy ``source`` so we can use it in a later section.
# Since ``source`` objects are designed to follow an immutable usage by default (like Numpy arrays),
# we can copy a source just by copying the object.
source_copy = source
logger.info("Obtain original images.")
imgs_seq_od = source.images[0].asnumpy()
logger.info("Perform phase flip to input images.")
source = source.phase_flip()
imgs_seq_pf = source.images[0].asnumpy()
logger.info(f"Downsample image size to {ds_img_size} X {ds_img_size}")
source = source.downsample(ds_img_size)
imgs_seq_ds = source.images[0].asnumpy()
logger.info("Normalize images to background noise.")
source = source.normalize_background()
imgs_seq_nb = source.images[0].asnumpy()
logger.info("Whiten noise of images")
source = source.whiten()
imgs_seq_wt = source.images[0].asnumpy()
logger.info("Invert the global density contrast if need")
source = source.invert_contrast()
imgs_seq_rc = source.images[0].asnumpy()
0%| | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:01<00:00, 1.01s/it]
100%|██████████| 1/1 [00:01<00:00, 1.01s/it]
0%| | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:01<00:00, 1.05s/it]
100%|██████████| 1/1 [00:01<00:00, 1.05s/it]
Plot First Image from Each Preprocess Step¶
Plot the first images.
logger.info("Plot first image from each preprocess steps")
idm = 0
plt.subplot(2, 3, 1)
plt.imshow(imgs_od[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("original image")
plt.subplot(2, 3, 2)
plt.imshow(imgs_seq_pf[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("phase flip")
plt.subplot(2, 3, 3)
plt.imshow(imgs_seq_ds[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("downsample")
plt.subplot(2, 3, 4)
plt.imshow(imgs_seq_nb[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("normalize background")
plt.subplot(2, 3, 5)
plt.imshow(imgs_seq_wt[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("noise whitening")
plt.subplot(2, 3, 6)
plt.imshow(imgs_seq_rc[idm], cmap="gray")
plt.colorbar(orientation="horizontal")
plt.title("invert contrast")
plt.tight_layout()
Apply Chained Preprocessing Techniques¶
Now we’ll apply the preprocessing in a chain syntax
# We'll reset our ``source`` to the reference copy we started with.
source = source_copy
logger.info("Perform phase flip to input images.")
logger.info(f"Downsample image size to {ds_img_size} X {ds_img_size}")
logger.info("Normalize images to background noise.")
logger.info("Whiten noise of images")
logger.info("Invert the global density contrast if need")
source = (
source.phase_flip()
.downsample(ds_img_size)
.normalize_background()
.whiten()
.invert_contrast()
)
# Assign the first image from the preprocessed chain.
imgs_chained = source.images[0].asnumpy()
# This preprocessing chain should correspond to applying each
# operation sequentially.
assert np.allclose(imgs_chained, imgs_seq_rc)
0%| | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:01<00:00, 1.02s/it]
100%|██████████| 1/1 [00:01<00:00, 1.02s/it]
0%| | 0/1 [00:00<?, ?it/s]
100%|██████████| 1/1 [00:01<00:00, 1.11s/it]
100%|██████████| 1/1 [00:01<00:00, 1.11s/it]
Total running time of the script: (0 minutes 7.971 seconds)