diff --git a/src/aspire/source/__init__.py b/src/aspire/source/__init__.py index 7a9a67cdd1..87f96482de 100644 --- a/src/aspire/source/__init__.py +++ b/src/aspire/source/__init__.py @@ -15,7 +15,12 @@ Multiply, Pipeline, ) -from aspire.operators import LambdaFilter, MultiplicativeFilter, PowerFilter +from aspire.operators import ( + IdentityFilter, + LambdaFilter, + MultiplicativeFilter, + PowerFilter, +) from aspire.storage import MrcStats, StarFile, StarFileBlock from aspire.utils import ensure from aspire.utils.coor_trans import grid_2d @@ -113,7 +118,7 @@ def __init__(self, L, n, dtype="double", metadata=None, memory=None): degrees=True, ) - self.unique_filters = None + self.unique_filters = [IdentityFilter()] self.generation_pipeline = Pipeline(xforms=None, memory=memory) self._metadata_out = None diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py new file mode 100644 index 0000000000..ff5f7035ab --- /dev/null +++ b/tutorials/examples/basic_image_array.py @@ -0,0 +1,173 @@ +import matplotlib.pyplot as plt +import numpy as np +from scipy import misc + +from aspire.image import Image +from aspire.image.xform import NoiseAdder +from aspire.noise import AnisotropicNoiseEstimator +from aspire.operators import FunctionFilter, ScalarFilter +from aspire.source import ArrayImageSource + +# ------------------------------------------------------------------------------ +# Lets get some basic image data as a numpy array. + +# Scipy ships with a portrait. +# We'll take the grayscale representation as floating point data. +stock_img = misc.face(gray=True).astype(np.float32) + +# Crop to a square +n_pixels = min(stock_img.shape) +stock_img = stock_img[0:n_pixels, 0:n_pixels] +# Normalize to [0,1] +stock_img /= np.max(stock_img) + + +# ------------------------------------------------------------------------------ +# Now that we have an example array, we'll begin using the ASPIRE toolkit. + +# First we'll make an ASPIRE Image instance out of our data. +# This is a light wrapper over the numpy array. Many ASPIRE internals +# are built around an Image class. + +# Construct the Image class by passing it an array of data. +img = Image(stock_img) +# Downsample (just to speeds things up) +new_resolution = img.res // 4 +img = img.downsample(new_resolution) + + +# We'll begin processing by adding some noise. +# We'd like to create uniform noise for a 2d image with prescibed variance, +noise_var = np.var(img.asnumpy()) * 5 +noise_filter = ScalarFilter(dim=2, value=noise_var) + +# Then create a NoiseAdder. +noise = NoiseAdder(seed=123, noise_filter=noise_filter) + +# We can apply the NoiseAdder to our image data. +img_with_noise = noise.forward(img) + +# We'll plot the original and first noisy image, +# because we only have one image in our Image stack right now. +fig, axs = plt.subplots(1, 2) +axs[0].imshow(img[0], cmap=plt.cm.gray) +axs[0].set_title("Starting Image") +axs[1].imshow(img_with_noise[0], cmap=plt.cm.gray) +axs[1].set_title("Noisy Image") +plt.show() + + +# ------------------------------------------------------------------------------ +# Great, now we have enough to try an experiment. +# This time we'll use a stack of images. +# +# In real use, you would probably bring your own array of images, +# or use a `Simulation` object. For now we'll create some arrays as before. +# For demonstration we'll setup a stack of n_imgs, +# with each image just being a copy of the data from `img`. +n_imgs = 128 +imgs_data = np.empty((n_imgs, img.res, img.res), dtype=np.float64) +for i in range(n_imgs): + imgs_data[i] = img[0] +imgs = Image(imgs_data) + + +# Lets say we want to add different kind of noise to the images. +# We can create our own function. Here we want to apply in two dimensions. +def noise_function(x, y): + return np.exp(-(x * x + y * y) / (2 * 0.3 ** 2)) + + +# We can create a custom filter from that function. +f = FunctionFilter(noise_function) +# And use the filter to add the noise to our stack of images. +noise_adder = NoiseAdder(seed=123, noise_filter=f) +imgs_with_noise = noise_adder.forward(imgs) + +# Let's see the first two noisy images. +# They should each display slightly different noise. +fig, axs = plt.subplots(2, 2) +for i, img in enumerate(imgs_with_noise[0:2]): + axs[0, i].imshow(img, cmap=plt.cm.gray) + axs[0, i].set_title(f"Custom Noisy Image {i}") + img_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(img))) + axs[1, i].imshow(np.log(1 + img_with_noise_f), cmap=plt.cm.gray) + axs[1, i].set_title(f"Custom Noisy Spectrum Image {i}") +plt.show() + + +# ------------------------------------------------------------------------------ +# Now we''l use an ASPIRE pipeline to Whiten the image stack. +# Here we will introduce our `Source` class and demonstrate applying a `xform`. + +# "Source" classes are what we use in processing pipelines. +# They provide a consistent interface to a variety of underlying data sources. +# In this case, we'll just use our Image in an ArrayImageSource to run a small experiment. +# +# If you were happy with the experiment design on an array of test data, +# the source is easily swapped out to something like RelionSource, +# which might point at a stack of images too large to fit in memory at once. +# The implementation of batching for memory management +# would be managed behind the scenes for you. + +imgs_src = ArrayImageSource(imgs_with_noise) + +# We'll copy the orginals for comparison later, before we process them further. +noisy_imgs_copy = imgs_src.images(0, n_imgs) + +# One of the tools we can use is a NoiseEstimator, +# which consumes from a Source. +noise_estimator = AnisotropicNoiseEstimator(imgs_src) + +# Once we have the estimator instance, +# we can use it in a transform applied to our Source. +imgs_src.whiten(noise_estimator.filter) + + +# Peek at two whitened images and their corresponding spectrum. +fig, axs = plt.subplots(2, 2) +for i, img in enumerate(imgs_src.images(0, 2)): + axs[0, i].imshow(img, cmap=plt.cm.gray) + axs[0, i].set_title(f"Whitened Noisy Image {i}") + img_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(img))) + axs[1, i].imshow(np.log(1 + img_with_noise_f), cmap=plt.cm.gray) + axs[1, i].set_title(f"Whitened Noisy Image Spectrum {i}") +plt.show() + + +# We'll also want to take a look at the spectrum power distribution. +# Since we just want to see the character of what is happening, +# I'll assume each pixel's contribution is placed at their lower left corner, +# and compute a crude radial profile. +# Code from a discussion at https://stackoverflow.com/questions/21242011/most-efficient-way-to-calculate-radial-profile. +def radial_profile(data): + y, x = np.indices((data.shape)) + # Distance from origin to lower left corner + r = np.sqrt(x ** 2 + y ** 2).astype(np.int) + binsum = np.bincount(r.ravel(), np.log(1 + data.ravel())) + bincount = np.bincount(r.ravel()) + # Return the mean per bin + return binsum / bincount + + +# Lets pickout several images and plot their the radial profile of their noise. +colors = ["r", "g", "b", "k", "c"] +for i, img in enumerate(imgs_src.images(0, len(colors))): + img_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(noisy_imgs_copy[i]))) + plt.plot( + radial_profile(img_with_noise_f), + color=colors[i], + label=f"Noisy Image Radial Profile {i}", + ) + + whitened_img_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(img))) + plt.plot( + radial_profile(whitened_img_with_noise_f), + color=colors[i], + linestyle="--", + label=f"Whitened Noisy Image Radial Profile {i}", + ) + +plt.title("Spectrum Profiles") +plt.legend() +plt.show()