From cb9ad5cf2cd75d348484995379f13661d3029c9b Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 8 Feb 2021 15:03:08 -0500 Subject: [PATCH 01/11] Stashing array image source end of day --- src/aspire/source/__init__.py | 8 +- tutorials/examples/basic_image_array.py | 128 ++++++++++++++++++++++++ 2 files changed, 135 insertions(+), 1 deletion(-) create mode 100644 tutorials/examples/basic_image_array.py diff --git a/src/aspire/source/__init__.py b/src/aspire/source/__init__.py index 7a9a67cdd1..ca4bead1e4 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 ( + LambdaFilter, + MultiplicativeFilter, + PowerFilter, + IdentityFilter, +) from aspire.storage import MrcStats, StarFile, StarFileBlock from aspire.utils import ensure from aspire.utils.coor_trans import grid_2d @@ -742,3 +747,4 @@ def __init__(self, im, metadata=None): L=im.res, n=im.n_images, dtype=im.dtype, metadata=metadata, memory=None ) self._cached_im = im + self.unique_filters = [IdentityFilter(dim=2)] diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py new file mode 100644 index 0000000000..eac51650ac --- /dev/null +++ b/tutorials/examples/basic_image_array.py @@ -0,0 +1,128 @@ +import numpy as np +from scipy import misc +import matplotlib.pyplot as plt +import colorednoise + +#------------------------------------------------------------------------------ +# Lets get some image data as a numpy array. + +# Grab some image, here we have a Procyon (common north american trash panda). +img_data = misc.face(gray=True) + +# Crop to a square +l = min(img_data.shape) +img_data = img_data[0:l,0:l] +print(f'Shape {img_data.shape}') +plt.imshow(img_data, cmap=plt.cm.gray) +plt.title('Original Image') +plt.show() + +# Initially, we will add some white noise to the image. +# Later we will try other types of noise in combination with whitening. + +mu = 0. +sigma = 256. +white = np.random.normal(mu, sigma, img_data.shape) + +img_data_withWhiteNoise = img_data + white +# plt.imshow(img_data_withWhiteNoise, cmap=plt.cm.gray) +# plt.title('With White Noise') +# plt.show() + + +#------------------------------------------------------------------------------ +# Now that we have an example array, we''ll begin using the ASPIRE toolkit. + +from aspire.image import Image +from aspire.source import ArrayImageSource +from aspire.noise import WhiteNoiseEstimator + +# 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 classes. + +# Construct the Image class by passing it the array data. +img = Image(img_data_withWhiteNoise) + +# "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 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. +# This would be managed behind the scenes for you. +img_src = ArrayImageSource(img) + +# ASPIRE's WhiteNoiseEstimator consumes from a Source +#noise_estimator = WhiteNoiseEstimator(img_src) + +# We can use that estimator to whiten all the images in the Source. +#img_src.whiten(noise_estimator.filter) + +# We can get a copy as numpy array instead of an ASPIRE source object. +#img_data_whitened = img_src.images(0,img_src.n).asnumpy() +# plt.imshow(img_data_whitened[0], cmap=plt.cm.gray) +# plt.title('Whitened') +# plt.show() + + +#------------------------------------------------------------------------------ +# Lets try adding a different sort of noise + +pink = (colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T) * sigma +pink2 = (colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T) * sigma +pink3 = (colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T) * sigma +brown = (colorednoise.powerlaw_psd_gaussian(2, img_data.shape) + colorednoise.powerlaw_psd_gaussian(2, img_data.shape).T) * sigma + +#noises = {'White': white, 'Pink': pink, 'Brown': brown} +noises = {'Pink1': pink, 'Pink2': pink2, 'Pink3': pink3} + +# We'll go through the same steps as before, +# but adding different noise to our original image data, +# and keeping some fourier space arrays to plot later. +myarray = np.zeros((3,l,l)) +myarray_f = np.zeros_like(myarray) +myarray_f_whitened = np.zeros_like(myarray) +for i, name in enumerate(sorted(noises)): + + myarray[i] = img_data + noises[name] + # Lets get the fourier space + myarray_f[i] = np.abs( + np.fft.fftshift( + np.fft.fft2(myarray[i]))) + +# Construct our Image and Source +images = Image(myarray) +imgs_src = ArrayImageSource(images) + +# Use ASPIRE pipeline to Whiten +noise_estimator = WhiteNoiseEstimator(imgs_src) +imgs_src.whiten(noise_estimator.filter) +img_data_whitened = imgs_src.images(0,imgs_src.n).asnumpy() + + +# Plot before and after whitening. +fig, axs = plt.subplots(3,4) +for i, name in enumerate(sorted(noises)): + + myarray[i] = img_data + noises[name] + # Lets get a picture of the whitened fourier space + myarray_f_whitened[i] = np.abs( + np.fft.fftshift( + np.fft.fft2(img_data_whitened[i]))) + + # and we can make the plots now + axs[i,0].imshow(myarray[i], cmap=plt.cm.gray) + axs[i,0].set_title(f"Image with {name} Noise") + + axs[i,1].imshow(np.log(1+myarray_f[i]), cmap=plt.cm.gray) + axs[i,1].set_title(f"{name} Noise Spectrum") + + axs[i,2].imshow(np.log(1+myarray_f_whitened[i]), cmap=plt.cm.gray) + axs[i,2].set_title(f"Whitened {name} Noise Spection") + + axs[i,3].imshow(img_data_whitened[i], cmap=plt.cm.gray) + axs[i,3].set_title(f"Image with {name} Noise Whitened") + +plt.show() + From 234b7aeb6ebf5bd6298e9b13a9205da3b53e34b0 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 8 Feb 2021 15:56:17 -0500 Subject: [PATCH 02/11] fixup flake8 and code style --- setup.py | 1 + src/aspire/source/__init__.py | 2 +- tutorials/examples/basic_image_array.py | 103 +++++++++++++----------- 3 files changed, 57 insertions(+), 49 deletions(-) diff --git a/setup.py b/setup.py index 45dc1ec37b..228c420c3a 100644 --- a/setup.py +++ b/setup.py @@ -24,6 +24,7 @@ def read(fname): author_email="devs.aspire@gmail.com", install_requires=[ "click", + "colorednoise", "finufft", "importlib_resources>=1.0.2", "joblib", diff --git a/src/aspire/source/__init__.py b/src/aspire/source/__init__.py index ca4bead1e4..152b49d161 100644 --- a/src/aspire/source/__init__.py +++ b/src/aspire/source/__init__.py @@ -16,10 +16,10 @@ Pipeline, ) from aspire.operators import ( + IdentityFilter, LambdaFilter, MultiplicativeFilter, PowerFilter, - IdentityFilter, ) from aspire.storage import MrcStats, StarFile, StarFileBlock from aspire.utils import ensure diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py index eac51650ac..5012fefebd 100644 --- a/tutorials/examples/basic_image_array.py +++ b/tutorials/examples/basic_image_array.py @@ -1,27 +1,31 @@ +import colorednoise +import matplotlib.pyplot as plt import numpy as np from scipy import misc -import matplotlib.pyplot as plt -import colorednoise -#------------------------------------------------------------------------------ +from aspire.image import Image +from aspire.noise import WhiteNoiseEstimator +from aspire.source import ArrayImageSource + +# ------------------------------------------------------------------------------ # Lets get some image data as a numpy array. # Grab some image, here we have a Procyon (common north american trash panda). img_data = misc.face(gray=True) # Crop to a square -l = min(img_data.shape) -img_data = img_data[0:l,0:l] -print(f'Shape {img_data.shape}') +n_pixels = min(img_data.shape) +img_data = img_data[0:n_pixels, 0:n_pixels] + plt.imshow(img_data, cmap=plt.cm.gray) -plt.title('Original Image') +plt.title("Original Image") plt.show() # Initially, we will add some white noise to the image. # Later we will try other types of noise in combination with whitening. -mu = 0. -sigma = 256. +mu = 0.0 +sigma = 256.0 white = np.random.normal(mu, sigma, img_data.shape) img_data_withWhiteNoise = img_data + white @@ -30,12 +34,8 @@ # plt.show() -#------------------------------------------------------------------------------ -# Now that we have an example array, we''ll begin using the ASPIRE toolkit. - -from aspire.image import Image -from aspire.source import ArrayImageSource -from aspire.noise import WhiteNoiseEstimator +# ------------------------------------------------------------------------------ +# 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 @@ -54,75 +54,82 @@ img_src = ArrayImageSource(img) # ASPIRE's WhiteNoiseEstimator consumes from a Source -#noise_estimator = WhiteNoiseEstimator(img_src) +# noise_estimator = WhiteNoiseEstimator(img_src) # We can use that estimator to whiten all the images in the Source. -#img_src.whiten(noise_estimator.filter) +# img_src.whiten(noise_estimator.filter) # We can get a copy as numpy array instead of an ASPIRE source object. -#img_data_whitened = img_src.images(0,img_src.n).asnumpy() +# img_data_whitened = img_src.images(0,img_src.n).asnumpy() # plt.imshow(img_data_whitened[0], cmap=plt.cm.gray) # plt.title('Whitened') # plt.show() -#------------------------------------------------------------------------------ +# ------------------------------------------------------------------------------ # Lets try adding a different sort of noise -pink = (colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T) * sigma -pink2 = (colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T) * sigma -pink3 = (colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T) * sigma -brown = (colorednoise.powerlaw_psd_gaussian(2, img_data.shape) + colorednoise.powerlaw_psd_gaussian(2, img_data.shape).T) * sigma - -#noises = {'White': white, 'Pink': pink, 'Brown': brown} -noises = {'Pink1': pink, 'Pink2': pink2, 'Pink3': pink3} +pink = ( + colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T +) * sigma +pink2 = ( + colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T +) * sigma +pink3 = ( + colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T +) * sigma +brown = ( + colorednoise.powerlaw_psd_gaussian(2, img_data.shape) + + colorednoise.powerlaw_psd_gaussian(2, img_data.shape).T +) * sigma + +# noises = {'White': white, 'Pink': pink, 'Brown': brown} +noises = {"Pink1": pink, "Pink2": pink2, "Pink3": pink3} # We'll go through the same steps as before, # but adding different noise to our original image data, # and keeping some fourier space arrays to plot later. -myarray = np.zeros((3,l,l)) +myarray = np.zeros((3, n_pixels, n_pixels)) myarray_f = np.zeros_like(myarray) myarray_f_whitened = np.zeros_like(myarray) for i, name in enumerate(sorted(noises)): myarray[i] = img_data + noises[name] # Lets get the fourier space - myarray_f[i] = np.abs( - np.fft.fftshift( - np.fft.fft2(myarray[i]))) + myarray_f[i] = np.abs(np.fft.fftshift(np.fft.fft2(myarray[i]))) # Construct our Image and Source images = Image(myarray) imgs_src = ArrayImageSource(images) -# Use ASPIRE pipeline to Whiten +# Use ASPIRE pipeline to Whiten noise_estimator = WhiteNoiseEstimator(imgs_src) imgs_src.whiten(noise_estimator.filter) -img_data_whitened = imgs_src.images(0,imgs_src.n).asnumpy() +img_data_whitened = imgs_src.images(0, imgs_src.n).asnumpy() # Plot before and after whitening. -fig, axs = plt.subplots(3,4) +fig, axs = plt.subplots(3, 4) for i, name in enumerate(sorted(noises)): myarray[i] = img_data + noises[name] # Lets get a picture of the whitened fourier space - myarray_f_whitened[i] = np.abs( - np.fft.fftshift( - np.fft.fft2(img_data_whitened[i]))) + myarray_f_whitened[i] = np.abs(np.fft.fftshift(np.fft.fft2(img_data_whitened[i]))) # and we can make the plots now - axs[i,0].imshow(myarray[i], cmap=plt.cm.gray) - axs[i,0].set_title(f"Image with {name} Noise") - - axs[i,1].imshow(np.log(1+myarray_f[i]), cmap=plt.cm.gray) - axs[i,1].set_title(f"{name} Noise Spectrum") - - axs[i,2].imshow(np.log(1+myarray_f_whitened[i]), cmap=plt.cm.gray) - axs[i,2].set_title(f"Whitened {name} Noise Spection") - - axs[i,3].imshow(img_data_whitened[i], cmap=plt.cm.gray) - axs[i,3].set_title(f"Image with {name} Noise Whitened") + axs[i, 0].imshow(myarray[i], cmap=plt.cm.gray) + axs[i, 0].set_title(f"Image with {name} Noise") + + axs[i, 1].imshow(np.log(1 + myarray_f[i]), cmap=plt.cm.gray) + axs[i, 1].set_title(f"{name} Noise Spectrum") + + axs[i, 2].imshow(np.log(1 + myarray_f_whitened[i]), cmap=plt.cm.gray) + axs[i, 2].set_title(f"Whitened {name} Noise Spection") + + axs[i, 3].imshow(img_data_whitened[i], cmap=plt.cm.gray) + axs[i, 3].set_title(f"Image with {name} Noise Whitened") plt.show() - From 6165ddbb130fef920769799cfa942f1bbbcf7eb8 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Wed, 10 Feb 2021 09:55:30 -0500 Subject: [PATCH 03/11] Plots spectrum profiles before and after whitening --- tutorials/examples/basic_image_array.py | 149 ++++++++++++++++-------- 1 file changed, 100 insertions(+), 49 deletions(-) diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py index 5012fefebd..a92b7764cd 100644 --- a/tutorials/examples/basic_image_array.py +++ b/tutorials/examples/basic_image_array.py @@ -1,24 +1,28 @@ import colorednoise +import matplotlib._color_data as mcd import matplotlib.pyplot as plt import numpy as np from scipy import misc +from skimage.measure import block_reduce from aspire.image import Image from aspire.noise import WhiteNoiseEstimator from aspire.source import ArrayImageSource # ------------------------------------------------------------------------------ -# Lets get some image data as a numpy array. +# Lets get some basic image data as a numpy array. -# Grab some image, here we have a Procyon (common north american trash panda). +# Scipy ships with a portrait of a Procyon (common north american trash panda). img_data = misc.face(gray=True) # Crop to a square n_pixels = min(img_data.shape) img_data = img_data[0:n_pixels, 0:n_pixels] +# Block (down)sample the image. +img_data = block_reduce(img_data, (2, 2)) plt.imshow(img_data, cmap=plt.cm.gray) -plt.title("Original Image") +plt.title("Starting Image") plt.show() # Initially, we will add some white noise to the image. @@ -28,10 +32,12 @@ sigma = 256.0 white = np.random.normal(mu, sigma, img_data.shape) -img_data_withWhiteNoise = img_data + white -# plt.imshow(img_data_withWhiteNoise, cmap=plt.cm.gray) -# plt.title('With White Noise') -# plt.show() +snr = np.var(img_data) / np.var(white) +print(f"Rough SNR {snr}") + +# We'll also compute the spectrum of the original image and white noise sample for later. +img_data_f = np.abs(np.fft.fftshift(np.fft.fft2(img_data))) +white_f = np.abs(np.fft.fftshift(np.fft.fft2(white))) # ------------------------------------------------------------------------------ @@ -41,8 +47,8 @@ # This is a light wrapper over the numpy array. Many ASPIRE internals # are built around an Image classes. -# Construct the Image class by passing it the array data. -img = Image(img_data_withWhiteNoise) +# Construct the Image class by passing it an array of data. +img = Image(img_data + white) # "Source" classes are what we use in processing pipelines. # They provide a consistent interface to a variety of underlying data sources. @@ -54,82 +60,127 @@ img_src = ArrayImageSource(img) # ASPIRE's WhiteNoiseEstimator consumes from a Source -# noise_estimator = WhiteNoiseEstimator(img_src) +noise_estimator = WhiteNoiseEstimator(img_src) # We can use that estimator to whiten all the images in the Source. -# img_src.whiten(noise_estimator.filter) +img_src.whiten(noise_estimator.filter) # We can get a copy as numpy array instead of an ASPIRE source object. -# img_data_whitened = img_src.images(0,img_src.n).asnumpy() -# plt.imshow(img_data_whitened[0], cmap=plt.cm.gray) -# plt.title('Whitened') -# plt.show() +img_data_whitened = img_src.images(0, img_src.n).asnumpy() # ------------------------------------------------------------------------------ -# Lets try adding a different sort of noise +# Lets try a small experiment, this time on a stack of thee images in an array. +# We'll go through the same steps as before. +# The following will generate additional distributions of noise. pink = ( colorednoise.powerlaw_psd_gaussian(1, img_data.shape) + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T ) * sigma -pink2 = ( - colorednoise.powerlaw_psd_gaussian(1, img_data.shape) - + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T -) * sigma -pink3 = ( - colorednoise.powerlaw_psd_gaussian(1, img_data.shape) - + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T -) * sigma +pink_f = np.abs(np.fft.fftshift(np.fft.fft2(pink))) + brown = ( colorednoise.powerlaw_psd_gaussian(2, img_data.shape) + colorednoise.powerlaw_psd_gaussian(2, img_data.shape).T ) * sigma +brown_f = np.abs(np.fft.fftshift(np.fft.fft2(brown))) + +# Storing noises in a dictionary for reference later. +noises = {"White": white, "Pink": pink, "Brown": brown} +noises_f = {"White": white_f, "Pink": pink_f, "Brown": brown_f} -# noises = {'White': white, 'Pink': pink, 'Brown': brown} -noises = {"Pink1": pink, "Pink2": pink2, "Pink3": pink3} -# We'll go through the same steps as before, -# but adding different noise to our original image data, -# and keeping some fourier space arrays to plot later. -myarray = np.zeros((3, n_pixels, n_pixels)) -myarray_f = np.zeros_like(myarray) -myarray_f_whitened = np.zeros_like(myarray) +# Setup some arrays to hold our data. +# In real use, you would probably bring your own stack of images, +# but we'll create some here as before. +stack = np.zeros((3, img_data.shape[-2], img_data.shape[-1])) +stack_f = np.zeros_like(stack) +stack_whitened_f = np.zeros_like(stack) +whitened_noises_f = dict() + for i, name in enumerate(sorted(noises)): + # Adding the different noises to our original image data. + stack[i] = img_data + noises[name] + # Compute and keep some fourier space arrays to plot later. + stack_f[i] = np.abs(np.fft.fftshift(np.fft.fft2(stack[i]))) + - myarray[i] = img_data + noises[name] - # Lets get the fourier space - myarray_f[i] = np.abs(np.fft.fftshift(np.fft.fft2(myarray[i]))) +# Construct our Image and Source. +images = Image(stack) +img_src = ArrayImageSource(images) -# Construct our Image and Source -images = Image(myarray) -imgs_src = ArrayImageSource(images) # Use ASPIRE pipeline to Whiten -noise_estimator = WhiteNoiseEstimator(imgs_src) -imgs_src.whiten(noise_estimator.filter) -img_data_whitened = imgs_src.images(0, imgs_src.n).asnumpy() +noise_estimator = WhiteNoiseEstimator(img_src) +img_src.whiten(noise_estimator.filter) +img_data_whitened = img_src.images(0, img_src.n).asnumpy() # Plot before and after whitening. fig, axs = plt.subplots(3, 4) for i, name in enumerate(sorted(noises)): - myarray[i] = img_data + noises[name] - # Lets get a picture of the whitened fourier space - myarray_f_whitened[i] = np.abs(np.fft.fftshift(np.fft.fft2(img_data_whitened[i]))) + stack[i] = img_data + noises[name] + # Lets save the whitened fourier space + stack_whitened_f[i] = np.abs(np.fft.fftshift(np.fft.fft2(img_data_whitened[i]))) - # and we can make the plots now - axs[i, 0].imshow(myarray[i], cmap=plt.cm.gray) + # and retrieve the whitened noise profile by subracting from the original signal + whitened_noises_f[name] = img_data_f - stack_whitened_f[i] + + # and we can make some plots now + axs[i, 0].imshow(stack[i], cmap=plt.cm.gray) axs[i, 0].set_title(f"Image with {name} Noise") - axs[i, 1].imshow(np.log(1 + myarray_f[i]), cmap=plt.cm.gray) + axs[i, 1].imshow(np.log(1 + stack_f[i]), cmap=plt.cm.gray) axs[i, 1].set_title(f"{name} Noise Spectrum") - axs[i, 2].imshow(np.log(1 + myarray_f_whitened[i]), cmap=plt.cm.gray) + axs[i, 2].imshow(np.log(1 + stack_whitened_f[i]), cmap=plt.cm.gray) axs[i, 2].set_title(f"Whitened {name} Noise Spection") axs[i, 3].imshow(img_data_whitened[i], cmap=plt.cm.gray) - axs[i, 3].set_title(f"Image with {name} Noise Whitened") + axs[i, 3].set_title(f"Image with Whitened {name} Noise") + +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 + + +# Setup some plot colors +legend_colors = { + "White": mcd.XKCD_COLORS["xkcd:black"], + "Pink": mcd.XKCD_COLORS["xkcd:pink"], + "Brown": mcd.XKCD_COLORS["xkcd:sienna"], +} + +# Loop through the sprectral profiles and plot. +for i, name in enumerate(sorted(noises)): + plt.plot(radial_profile(noises_f[name]), legend_colors[name], label=name) + plt.plot( + radial_profile(whitened_noises_f[name]), + color=f"{legend_colors[name]}", + linestyle="--", + label=f"Whitened {name}", + ) + +plt.title(f"Spectrum Profiles") +plt.legend() plt.show() + +# At this point we should see that ASPIRE's whitening procedure has +# effected the distribution of noise. +# In other tutorials a Simulation source is generally used +# in place of constructing your own image stacks. From b1b2ac8f288462046e5af9ecc193cebac65c1a51 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Wed, 10 Feb 2021 10:20:16 -0500 Subject: [PATCH 04/11] dim=2 for Identity filter was not needed I must have been confused --- src/aspire/source/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aspire/source/__init__.py b/src/aspire/source/__init__.py index 152b49d161..4512875b47 100644 --- a/src/aspire/source/__init__.py +++ b/src/aspire/source/__init__.py @@ -747,4 +747,4 @@ def __init__(self, im, metadata=None): L=im.res, n=im.n_images, dtype=im.dtype, metadata=metadata, memory=None ) self._cached_im = im - self.unique_filters = [IdentityFilter(dim=2)] + self.unique_filters = [IdentityFilter()] From c2a0cd9cedf82d6d5da921d82c3f2234fe06211a Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Wed, 10 Feb 2021 10:25:44 -0500 Subject: [PATCH 05/11] code cleanup --- tutorials/examples/basic_image_array.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py index a92b7764cd..a97b536a01 100644 --- a/tutorials/examples/basic_image_array.py +++ b/tutorials/examples/basic_image_array.py @@ -32,8 +32,6 @@ sigma = 256.0 white = np.random.normal(mu, sigma, img_data.shape) -snr = np.var(img_data) / np.var(white) -print(f"Rough SNR {snr}") # We'll also compute the spectrum of the original image and white noise sample for later. img_data_f = np.abs(np.fft.fftshift(np.fft.fft2(img_data))) @@ -167,7 +165,7 @@ def radial_profile(data): } # Loop through the sprectral profiles and plot. -for i, name in enumerate(sorted(noises)): +for name in sorted(noises): plt.plot(radial_profile(noises_f[name]), legend_colors[name], label=name) plt.plot( radial_profile(whitened_noises_f[name]), @@ -176,7 +174,7 @@ def radial_profile(data): label=f"Whitened {name}", ) -plt.title(f"Spectrum Profiles") +plt.title("Spectrum Profiles") plt.legend() plt.show() From 986526f9fd7532b148ea1f8bcde9cc2ab9637739 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 11 Feb 2021 07:31:17 -0500 Subject: [PATCH 06/11] Move uinique_filter init to base class --- src/aspire/source/__init__.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/aspire/source/__init__.py b/src/aspire/source/__init__.py index 4512875b47..87f96482de 100644 --- a/src/aspire/source/__init__.py +++ b/src/aspire/source/__init__.py @@ -118,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 @@ -747,4 +747,3 @@ def __init__(self, im, metadata=None): L=im.res, n=im.n_images, dtype=im.dtype, metadata=metadata, memory=None ) self._cached_im = im - self.unique_filters = [IdentityFilter()] From 7c57dfe00d9c5fcaba303e03c3694924f0a94155 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 11 Feb 2021 15:54:33 -0500 Subject: [PATCH 07/11] remove colorednoise, replaced with NoiseAdder --- setup.py | 1 - 1 file changed, 1 deletion(-) diff --git a/setup.py b/setup.py index 228c420c3a..45dc1ec37b 100644 --- a/setup.py +++ b/setup.py @@ -24,7 +24,6 @@ def read(fname): author_email="devs.aspire@gmail.com", install_requires=[ "click", - "colorednoise", "finufft", "importlib_resources>=1.0.2", "joblib", From e3a72cccfde47fb03b1a8e40e7de503294f55634 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 11 Feb 2021 15:55:07 -0500 Subject: [PATCH 08/11] Migrate image array example to ASPIRE NoiseAdder --- tutorials/examples/basic_image_array.py | 222 ++++++++++++------------ 1 file changed, 112 insertions(+), 110 deletions(-) diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py index a97b536a01..a53bec1580 100644 --- a/tutorials/examples/basic_image_array.py +++ b/tutorials/examples/basic_image_array.py @@ -1,41 +1,33 @@ -import colorednoise -import matplotlib._color_data as mcd import matplotlib.pyplot as plt import numpy as np from scipy import misc from skimage.measure import block_reduce from aspire.image import Image +from aspire.image.xform import NoiseAdder from aspire.noise import WhiteNoiseEstimator +from aspire.operators import ScalarFilter from aspire.source import ArrayImageSource # ------------------------------------------------------------------------------ # Lets get some basic image data as a numpy array. -# Scipy ships with a portrait of a Procyon (common north american trash panda). -img_data = misc.face(gray=True) +# Scipy ships with a portrait. +# We'll take the grayscale representation as floating point data. +stock_img_data = misc.face(gray=True).astype(np.float32) # Crop to a square -n_pixels = min(img_data.shape) -img_data = img_data[0:n_pixels, 0:n_pixels] -# Block (down)sample the image. -img_data = block_reduce(img_data, (2, 2)) +n_pixels = min(stock_img_data.shape) +stock_img_data = stock_img_data[0:n_pixels, 0:n_pixels] +# downsample (speeds things up) +stock_img_data = block_reduce(stock_img_data, (4, 4)) -plt.imshow(img_data, cmap=plt.cm.gray) +plt.imshow(stock_img_data, cmap=plt.cm.gray) plt.title("Starting Image") plt.show() -# Initially, we will add some white noise to the image. -# Later we will try other types of noise in combination with whitening. - -mu = 0.0 -sigma = 256.0 -white = np.random.normal(mu, sigma, img_data.shape) - - -# We'll also compute the spectrum of the original image and white noise sample for later. -img_data_f = np.abs(np.fft.fftshift(np.fft.fft2(img_data))) -white_f = np.abs(np.fft.fftshift(np.fft.fft2(white))) +# We'll also compute the spectrum of the original image sample for later. +stock_img_data_f = np.abs(np.fft.fftshift(np.fft.fft2(stock_img_data))) # ------------------------------------------------------------------------------ @@ -46,102 +38,119 @@ # are built around an Image classes. # Construct the Image class by passing it an array of data. -img = Image(img_data + white) +img = Image(stock_img_data) -# "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 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. -# This would be managed behind the scenes for you. -img_src = ArrayImageSource(img) +# We'll begin processing by adding some noise. +# We'd like to create uniform noise for a 2d image with prescibed variance, +# say yielding an SNR around 10. +noise_var = np.var(stock_img_data) * 10.0 +noise_filter = ScalarFilter(dim=2, value=noise_var) -# ASPIRE's WhiteNoiseEstimator consumes from a Source -noise_estimator = WhiteNoiseEstimator(img_src) +# Then create a NoiseAdder, +noise = NoiseAdder(seed=123, noise_filter=noise_filter) -# We can use that estimator to whiten all the images in the Source. -img_src.whiten(noise_estimator.filter) +# which we can apply to our image data. +img_with_noise = noise.forward(img) -# We can get a copy as numpy array instead of an ASPIRE source object. -img_data_whitened = img_src.images(0, img_src.n).asnumpy() +# We'll plot the first image (we only have one in our stack here). +plt.imshow(img_with_noise[0], cmap=plt.cm.gray) +plt.title("Noisy Image") +plt.show() # ------------------------------------------------------------------------------ -# Lets try a small experiment, this time on a stack of thee images in an array. -# We'll go through the same steps as before. - -# The following will generate additional distributions of noise. -pink = ( - colorednoise.powerlaw_psd_gaussian(1, img_data.shape) - + colorednoise.powerlaw_psd_gaussian(1, img_data.shape).T -) * sigma -pink_f = np.abs(np.fft.fftshift(np.fft.fft2(pink))) - -brown = ( - colorednoise.powerlaw_psd_gaussian(2, img_data.shape) - + colorednoise.powerlaw_psd_gaussian(2, img_data.shape).T -) * sigma -brown_f = np.abs(np.fft.fftshift(np.fft.fft2(brown))) - -# Storing noises in a dictionary for reference later. -noises = {"White": white, "Pink": pink, "Brown": brown} -noises_f = {"White": white_f, "Pink": pink_f, "Brown": brown_f} - - -# Setup some arrays to hold our data. +# Lets try an experiment, this time on a stack of images in an array. # In real use, you would probably bring your own stack of images, # but we'll create some here as before. -stack = np.zeros((3, img_data.shape[-2], img_data.shape[-1])) -stack_f = np.zeros_like(stack) -stack_whitened_f = np.zeros_like(stack) -whitened_noises_f = dict() - -for i, name in enumerate(sorted(noises)): - # Adding the different noises to our original image data. - stack[i] = img_data + noises[name] - # Compute and keep some fourier space arrays to plot later. - stack_f[i] = np.abs(np.fft.fftshift(np.fft.fft2(stack[i]))) - +n_imgs = 128 +imgs_data = np.empty( + (n_imgs, stock_img_data.shape[-2], stock_img_data.shape[-1]), dtype=np.float64 +) +for i in range(n_imgs): + imgs_data[i] = stock_img_data +imgs = Image(imgs_data) + +# Similar to before, we'll construct noise with a constant variance, +# but now for the whole stack. +noise_var = np.var(stock_img_data) * 10.0 +noise_adder = NoiseAdder(seed=123, noise_filter=ScalarFilter(dim=2, value=noise_var)) +imgs_with_noise = noise_adder.forward(imgs) + +# We can check some images, which should be the same up to noise. +n_check = 2 +fig, axs = plt.subplots(1, n_check) +for i in range(n_check): + axs[i].imshow(imgs_with_noise[i], cmap=plt.cm.gray) + axs[i].set_title(f"Noisy Image {i}") +plt.show() -# Construct our Image and Source. -images = Image(stack) -img_src = ArrayImageSource(images) +# Lets say we want to add additonal noise to half the images, every other image. +indices = range(1, n_imgs, 2) +noise_adder = NoiseAdder( + seed=123, noise_filter=ScalarFilter(dim=2, value=4 * noise_var) +) +# We can use the "indices" to selectively apply our xform. +# Note, that the xform will return a dense Image, so we need to match dimensions +# There is an IndexdXform that provides an alternative to this. +imgs_with_noise[indices] = noise_adder.forward(imgs_with_noise, indices=indices)[ + : len(indices) +] + +# We can check now that the second image has a different noise profile. +n_check = 2 +fig, axs = plt.subplots(1, n_check) +for i in range(n_check): + axs[i].imshow(imgs_with_noise[i], cmap=plt.cm.gray) + axs[i].set_title(f"Noisy Image {i}") +plt.show() +# ------------------------------------------------------------------------------ # Use ASPIRE pipeline to Whiten -noise_estimator = WhiteNoiseEstimator(img_src) -img_src.whiten(noise_estimator.filter) -img_data_whitened = img_src.images(0, img_src.n).asnumpy() - -# Plot before and after whitening. -fig, axs = plt.subplots(3, 4) -for i, name in enumerate(sorted(noises)): +# "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 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. +# This would be managed behind the scenes for you. +imgs_src = ArrayImageSource(imgs_with_noise) - stack[i] = img_data + noises[name] - # Lets save the whitened fourier space - stack_whitened_f[i] = np.abs(np.fft.fftshift(np.fft.fft2(img_data_whitened[i]))) - # and retrieve the whitened noise profile by subracting from the original signal - whitened_noises_f[name] = img_data_f - stack_whitened_f[i] +# One of the tools we can use is a NoiseEstimator, +# which consumes from a Source. +noise_estimator = WhiteNoiseEstimator(imgs_src) - # and we can make some plots now - axs[i, 0].imshow(stack[i], cmap=plt.cm.gray) - axs[i, 0].set_title(f"Image with {name} Noise") +# Once we have the estimator instance, +# we can use it in a transform applied to our Source. +imgs_src.whiten(noise_estimator.filter) - axs[i, 1].imshow(np.log(1 + stack_f[i]), cmap=plt.cm.gray) - axs[i, 1].set_title(f"{name} Noise Spectrum") +# We can get numpy arrays from our source if we want them. +whitened_imgs_data = imgs_src.images(0, imgs_src.n).asnumpy() - axs[i, 2].imshow(np.log(1 + stack_whitened_f[i]), cmap=plt.cm.gray) - axs[i, 2].set_title(f"Whitened {name} Noise Spection") +fig, axs = plt.subplots(2, 2) +for i in range(axs.shape[1]): + axs[0, i].imshow(imgs_with_noise[i], cmap=plt.cm.gray) + axs[0, i].set_title(f"Noisy Image {i}") + axs[1, i].imshow(whitened_imgs_data[i], cmap=plt.cm.gray) + axs[1, i].set_title(f"Whitened Noisy Image {i}") +plt.show() - axs[i, 3].imshow(img_data_whitened[i], cmap=plt.cm.gray) - axs[i, 3].set_title(f"Image with Whitened {name} Noise") +# Okay, so lets look at the spectrum +fig, axs = plt.subplots(2, 2) +for i in range(axs.shape[1]): + imgs_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(imgs_with_noise[i]))) + axs[0, i].imshow(np.log(1 + imgs_with_noise_f), cmap=plt.cm.gray) + axs[0, i].set_title(f"Spectrum of Noisy Image {i}") + whitened_imgs_f = np.abs(np.fft.fftshift(np.fft.fft2(whitened_imgs_data[i]))) + axs[1, i].imshow(np.log(1 + whitened_imgs_f), cmap=plt.cm.gray) + axs[1, i].set_title(f"Whitened Noisy Image {i}") plt.show() +# Stil hard to tell. # 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, @@ -157,28 +166,21 @@ def radial_profile(data): return binsum / bincount -# Setup some plot colors -legend_colors = { - "White": mcd.XKCD_COLORS["xkcd:black"], - "Pink": mcd.XKCD_COLORS["xkcd:pink"], - "Brown": mcd.XKCD_COLORS["xkcd:sienna"], -} +# Still not happy with this part +colors = ["r", "g"] +for i in range(axs.shape[1]): + imgs_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(imgs_with_noise[i]))) + noise_f = imgs_with_noise_f - stock_img_data_f + plt.plot(radial_profile(noise_f), color=colors[i], label=f"PSD of Noisy Image {i}") -# Loop through the sprectral profiles and plot. -for name in sorted(noises): - plt.plot(radial_profile(noises_f[name]), legend_colors[name], label=name) + whitened_imgs_f = np.abs(np.fft.fftshift(np.fft.fft2(whitened_imgs_data[i]))) + whitened_noise_f = stock_img_data_f - whitened_imgs_f plt.plot( - radial_profile(whitened_noises_f[name]), - color=f"{legend_colors[name]}", + radial_profile(whitened_noise_f), + color=colors[i], linestyle="--", - label=f"Whitened {name}", + label=f"PSD Whitened Noisy Image {i}", ) - plt.title("Spectrum Profiles") plt.legend() plt.show() - -# At this point we should see that ASPIRE's whitening procedure has -# effected the distribution of noise. -# In other tutorials a Simulation source is generally used -# in place of constructing your own image stacks. From 5849f428a4dab139def6333cd789e5eb8c5608b8 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 12 Feb 2021 15:43:28 -0500 Subject: [PATCH 09/11] Use suggested noise function and estimator in example --- tutorials/examples/basic_image_array.py | 165 +++++++++++------------- 1 file changed, 75 insertions(+), 90 deletions(-) diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py index a53bec1580..22c1709632 100644 --- a/tutorials/examples/basic_image_array.py +++ b/tutorials/examples/basic_image_array.py @@ -5,8 +5,8 @@ from aspire.image import Image from aspire.image.xform import NoiseAdder -from aspire.noise import WhiteNoiseEstimator -from aspire.operators import ScalarFilter +from aspire.noise import AnisotropicNoiseEstimator +from aspire.operators import FunctionFilter, ScalarFilter from aspire.source import ArrayImageSource # ------------------------------------------------------------------------------ @@ -14,20 +14,15 @@ # Scipy ships with a portrait. # We'll take the grayscale representation as floating point data. -stock_img_data = misc.face(gray=True).astype(np.float32) +stock_img = misc.face(gray=True).astype(np.float32) # Crop to a square -n_pixels = min(stock_img_data.shape) -stock_img_data = stock_img_data[0:n_pixels, 0:n_pixels] -# downsample (speeds things up) -stock_img_data = block_reduce(stock_img_data, (4, 4)) - -plt.imshow(stock_img_data, cmap=plt.cm.gray) -plt.title("Starting Image") -plt.show() - -# We'll also compute the spectrum of the original image sample for later. -stock_img_data_f = np.abs(np.fft.fftshift(np.fft.fft2(stock_img_data))) +n_pixels = min(stock_img.shape) +stock_img = stock_img[0:n_pixels, 0:n_pixels] +# Downsample (just to speeds things up) +stock_img = block_reduce(stock_img, (4, 4)) +# Normalize to [0,1] +stock_img /= np.max(stock_img) # ------------------------------------------------------------------------------ @@ -38,119 +33,106 @@ # are built around an Image classes. # Construct the Image class by passing it an array of data. -img = Image(stock_img_data) +img = Image(stock_img) # We'll begin processing by adding some noise. # We'd like to create uniform noise for a 2d image with prescibed variance, -# say yielding an SNR around 10. -noise_var = np.var(stock_img_data) * 10.0 +noise_var = np.var(stock_img) * 5 noise_filter = ScalarFilter(dim=2, value=noise_var) -# Then create a NoiseAdder, +# Then create a NoiseAdder. noise = NoiseAdder(seed=123, noise_filter=noise_filter) -# which we can apply to our image data. +# We can apply the NoiseAdder to our image data. img_with_noise = noise.forward(img) -# We'll plot the first image (we only have one in our stack here). -plt.imshow(img_with_noise[0], cmap=plt.cm.gray) -plt.title("Noisy Image") +# We'll plot the original and first noisy image (we only have one in our stack here). +fig, axs = plt.subplots(1, 2) +axs[0].imshow(stock_img, 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() # ------------------------------------------------------------------------------ -# Lets try an experiment, this time on a stack of images in an array. -# In real use, you would probably bring your own stack of images, -# but we'll create some here as before. +# Great, now we have enough to try an experiment, +# but this time on 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. n_imgs = 128 imgs_data = np.empty( - (n_imgs, stock_img_data.shape[-2], stock_img_data.shape[-1]), dtype=np.float64 + (n_imgs, stock_img.shape[-2], stock_img.shape[-1]), dtype=np.float64 ) for i in range(n_imgs): - imgs_data[i] = stock_img_data + imgs_data[i] = stock_img imgs = Image(imgs_data) -# Similar to before, we'll construct noise with a constant variance, -# but now for the whole stack. -noise_var = np.var(stock_img_data) * 10.0 -noise_adder = NoiseAdder(seed=123, noise_filter=ScalarFilter(dim=2, value=noise_var)) -imgs_with_noise = noise_adder.forward(imgs) -# We can check some images, which should be the same up to noise. -n_check = 2 -fig, axs = plt.subplots(1, n_check) -for i in range(n_check): - axs[i].imshow(imgs_with_noise[i], cmap=plt.cm.gray) - axs[i].set_title(f"Noisy Image {i}") -plt.show() +# 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)) -# Lets say we want to add additonal noise to half the images, every other image. -indices = range(1, n_imgs, 2) -noise_adder = NoiseAdder( - seed=123, noise_filter=ScalarFilter(dim=2, value=4 * noise_var) -) -# We can use the "indices" to selectively apply our xform. -# Note, that the xform will return a dense Image, so we need to match dimensions -# There is an IndexdXform that provides an alternative to this. -imgs_with_noise[indices] = noise_adder.forward(imgs_with_noise, indices=indices)[ - : len(indices) -] - -# We can check now that the second image has a different noise profile. -n_check = 2 -fig, axs = plt.subplots(1, n_check) -for i in range(n_check): - axs[i].imshow(imgs_with_noise[i], cmap=plt.cm.gray) - axs[i].set_title(f"Noisy Image {i}") + +# 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() # ------------------------------------------------------------------------------ # Use ASPIRE pipeline to Whiten +# 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 to run a small experiment. +# 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. -# This would be managed behind the scenes for you. +# 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 = WhiteNoiseEstimator(imgs_src) +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) -# We can get numpy arrays from our source if we want them. -whitened_imgs_data = imgs_src.images(0, imgs_src.n).asnumpy() - -fig, axs = plt.subplots(2, 2) -for i in range(axs.shape[1]): - axs[0, i].imshow(imgs_with_noise[i], cmap=plt.cm.gray) - axs[0, i].set_title(f"Noisy Image {i}") - axs[1, i].imshow(whitened_imgs_data[i], cmap=plt.cm.gray) - axs[1, i].set_title(f"Whitened Noisy Image {i}") -plt.show() - -# Okay, so lets look at the spectrum +# Peek at two whitened images and their corresponding spectrum. fig, axs = plt.subplots(2, 2) -for i in range(axs.shape[1]): - imgs_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(imgs_with_noise[i]))) - axs[0, i].imshow(np.log(1 + imgs_with_noise_f), cmap=plt.cm.gray) - axs[0, i].set_title(f"Spectrum of Noisy Image {i}") - whitened_imgs_f = np.abs(np.fft.fftshift(np.fft.fft2(whitened_imgs_data[i]))) - axs[1, i].imshow(np.log(1 + whitened_imgs_f), cmap=plt.cm.gray) - axs[1, i].set_title(f"Whitened Noisy Image {i}") +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() -# Stil hard to tell. # 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, @@ -166,21 +148,24 @@ def radial_profile(data): return binsum / bincount -# Still not happy with this part -colors = ["r", "g"] -for i in range(axs.shape[1]): - imgs_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(imgs_with_noise[i]))) - noise_f = imgs_with_noise_f - stock_img_data_f - plt.plot(radial_profile(noise_f), color=colors[i], label=f"PSD of Noisy Image {i}") +# 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_imgs_f = np.abs(np.fft.fftshift(np.fft.fft2(whitened_imgs_data[i]))) - whitened_noise_f = stock_img_data_f - whitened_imgs_f + whitened_img_with_noise_f = np.abs(np.fft.fftshift(np.fft.fft2(img))) plt.plot( - radial_profile(whitened_noise_f), + radial_profile(whitened_img_with_noise_f), color=colors[i], linestyle="--", - label=f"PSD Whitened Noisy Image {i}", + label=f"Whitened Noisy Image Radial Profile {i}", ) + plt.title("Spectrum Profiles") plt.legend() plt.show() From a8258f0e806924ffd79ed4a093353d757f78a179 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 18 Feb 2021 10:56:23 -0500 Subject: [PATCH 10/11] slight changes to tutorial comments --- tutorials/examples/basic_image_array.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py index 22c1709632..a983ed9f5b 100644 --- a/tutorials/examples/basic_image_array.py +++ b/tutorials/examples/basic_image_array.py @@ -30,7 +30,7 @@ # 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 classes. +# are built around an Image class. # Construct the Image class by passing it an array of data. img = Image(stock_img) @@ -46,7 +46,8 @@ # We can apply the NoiseAdder to our image data. img_with_noise = noise.forward(img) -# We'll plot the original and first noisy image (we only have one in our stack here). +# 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(stock_img, cmap=plt.cm.gray) axs[0].set_title("Starting Image") @@ -56,11 +57,13 @@ # ------------------------------------------------------------------------------ -# Great, now we have enough to try an experiment, -# but this time on a stack of images. +# 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 stock_img data. n_imgs = 128 imgs_data = np.empty( (n_imgs, stock_img.shape[-2], stock_img.shape[-1]), dtype=np.float64 @@ -95,7 +98,7 @@ def noise_function(x, y): # ------------------------------------------------------------------------------ -# Use ASPIRE pipeline to Whiten +# 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. From c34f4243bd8d05a14974024439560de8acf8565c Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 19 Feb 2021 08:10:23 -0500 Subject: [PATCH 11/11] replace use of skimage with ASPIRE --- tutorials/examples/basic_image_array.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/tutorials/examples/basic_image_array.py b/tutorials/examples/basic_image_array.py index a983ed9f5b..ff5f7035ab 100644 --- a/tutorials/examples/basic_image_array.py +++ b/tutorials/examples/basic_image_array.py @@ -1,7 +1,6 @@ import matplotlib.pyplot as plt import numpy as np from scipy import misc -from skimage.measure import block_reduce from aspire.image import Image from aspire.image.xform import NoiseAdder @@ -19,8 +18,6 @@ # Crop to a square n_pixels = min(stock_img.shape) stock_img = stock_img[0:n_pixels, 0:n_pixels] -# Downsample (just to speeds things up) -stock_img = block_reduce(stock_img, (4, 4)) # Normalize to [0,1] stock_img /= np.max(stock_img) @@ -34,10 +31,14 @@ # 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(stock_img) * 5 +noise_var = np.var(img.asnumpy()) * 5 noise_filter = ScalarFilter(dim=2, value=noise_var) # Then create a NoiseAdder. @@ -49,7 +50,7 @@ # 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(stock_img, cmap=plt.cm.gray) +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") @@ -63,13 +64,11 @@ # 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 stock_img data. +# with each image just being a copy of the data from `img`. n_imgs = 128 -imgs_data = np.empty( - (n_imgs, stock_img.shape[-2], stock_img.shape[-1]), dtype=np.float64 -) +imgs_data = np.empty((n_imgs, img.res, img.res), dtype=np.float64) for i in range(n_imgs): - imgs_data[i] = stock_img + imgs_data[i] = img[0] imgs = Image(imgs_data)