From baa98d70b63bd39c9820a6828c68cb7cc862f50e Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Mon, 22 Jul 2024 13:51:54 -0400 Subject: [PATCH 01/25] Added backproject script and stub in the image folder --- src/aspire/image/line.py | 96 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) create mode 100644 src/aspire/image/line.py diff --git a/src/aspire/image/line.py b/src/aspire/image/line.py new file mode 100644 index 0000000000..078e6b37c3 --- /dev/null +++ b/src/aspire/image/line.py @@ -0,0 +1,96 @@ +import aspire +import numpy as np + + +class Line: + def __init__(self, data, dtype = np.float64): + """ + Initialize a Line Object. Change later (similar intuition from Image class) + Question: Is it a line or collection of line object? + + :param data: Numpy array containing image data with shape + `(..., resolution, resolution)`. + :param dtype: Optionally cast `data` to this dtype. + Defaults to `data.dtype`. + """ + self.dtype = np.dtype(dtype) + if data.ndim == 2: + data = data[np.newaxis, :, :] + if data.ndim < 3: + raise('Projection Dimensions should be more than Three-Dimensions') + self._data = data.astype(self.dtype, copy=False) + self.stack_shape = self._data.shape[:-2] + self.n_images = self._data.shape[0] #broken for higher dimensional stacks + self.n_lines = self._data.shape[-2] + self.n_points = self._data.shape[-1] + # self.n_dim = (self._data.shape[1], self._data.shape[2]) + + def __str__(self): + return f"Line(n_images = {self.n_images}, n_points = {self.n_points})" + + @property + def stack(self): + return self.n_images + + # talk about angles and if they're supposed to be a certain input/output + # why this method and not another (explain design choices) + def back_project(self, angles): + """ + Back Projection Method for a single stack of lines. + + :param filter_name: string, optional + Filter used in frequency domain filtering. Assign None to use no filter. + :param angles: array + assuming not perfectly radial angles + :return: stack of reconstructed + """ + assert len(angles) == self.n_lines, "Angles must match the number of lines." + original_stack = self.stack_shape + n_angles = len(angles) + + ## our implementation + n_img, n_angles, n_rad = sinogram.shape + assert n_angles == len(angles), "gonna have a bad time" + L = n_rad + sinogram = np.fft.ifftshift(self.n_images, axes= -1) + sinogram_ft = np.fft.rfft(self.n_images, axis= -1) + + #grid generation + y_idx = np.fft.rfftfreq(n_rad) * np.pi * 2 + n_real_points = len(y_idx) + pts = np.empty((2, len(angles), n_real_points), dtype=self.dtype) + pts[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] + pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] + + imgs = aspire.nufft.anufft( + sinogram_ft.reshape(n_img, -1), + pts.reshape(2, n_real_points * len(angles)), + sz=(L, L), + real=True + ).reshape(n_img, L, L) + + return aspire.image.Image(imgs) + + def image_filter(self, filter_name, projections): + """ + Filter Method for projections. Will apply filter to line projection to get collection of projections (ramp, cosine, ... , etc.) + :param projections: Collection of line projections that need to be filtered. + :return: Filtered Projections. + """ + if projections is None: + raise ValueError('The input projections must not be None') + + filter_types = ('ramp', 'shepp-i logan', 'cosine', 'hamming', 'hann', None) + if filter_name is not filter_type: + raise ValueError(f"Unknown filter: {filter_name}") + + # skimage filter + fourier_filter = _get_fourier_filter(projection_size_padded, filter_name) + projection = fft(img, axis=0) * fourier_filter + radon_filtered = np.real(ifft(projection, axis=0)[:img_shape, :]) + + """ + step 0: Look more into filter function from skimage + thoughts: apply filter to each point + """ + pass From 170a2194a626d56c8cdefd3610b0002a601ad097 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 23 Jul 2024 03:42:36 -0400 Subject: [PATCH 02/25] Added One-Dimension Test for Backproject --- src/aspire/image/line.py | 71 +++++++++++++++------------------------- tests/test_sinogram.py | 32 ++++++++++++++++++ 2 files changed, 59 insertions(+), 44 deletions(-) diff --git a/src/aspire/image/line.py b/src/aspire/image/line.py index 078e6b37c3..9f2b30c800 100644 --- a/src/aspire/image/line.py +++ b/src/aspire/image/line.py @@ -1,9 +1,19 @@ -import aspire +import logging +import os +from pathlib import Path + import numpy as np +import aspire +from aspire.image import Image +from aspire.nufft import anufft, nufft + +# noticed a lot of classes had these already, might be helpful for pathing, logging info, etc. (inc. os, logging) +logger = logging.getLogger(__name__) + class Line: - def __init__(self, data, dtype = np.float64): + def __init__(self, data, dtype=np.float64): """ Initialize a Line Object. Change later (similar intuition from Image class) Question: Is it a line or collection of line object? @@ -17,11 +27,11 @@ def __init__(self, data, dtype = np.float64): if data.ndim == 2: data = data[np.newaxis, :, :] if data.ndim < 3: - raise('Projection Dimensions should be more than Three-Dimensions') + raise ("Projection Dimensions should be more than Three-Dimensions") self._data = data.astype(self.dtype, copy=False) self.stack_shape = self._data.shape[:-2] - self.n_images = self._data.shape[0] #broken for higher dimensional stacks - self.n_lines = self._data.shape[-2] + self.n_images = self._data.shape[0] # broken for higher dimensional stacks + self.n_lines = self._data.shape[-2] self.n_points = self._data.shape[-1] # self.n_dim = (self._data.shape[1], self._data.shape[2]) @@ -32,65 +42,38 @@ def __str__(self): def stack(self): return self.n_images - # talk about angles and if they're supposed to be a certain input/output - # why this method and not another (explain design choices) def back_project(self, angles): """ Back Projection Method for a single stack of lines. - + :param filter_name: string, optional Filter used in frequency domain filtering. Assign None to use no filter. :param angles: array assuming not perfectly radial angles :return: stack of reconstructed """ - assert len(angles) == self.n_lines, "Angles must match the number of lines." - original_stack = self.stack_shape - n_angles = len(angles) - - ## our implementation + sinogram = self._data n_img, n_angles, n_rad = sinogram.shape - assert n_angles == len(angles), "gonna have a bad time" + assert n_angles == len( + angles + ), "Number of angles must match the number of projections" + L = n_rad - sinogram = np.fft.ifftshift(self.n_images, axes= -1) - sinogram_ft = np.fft.rfft(self.n_images, axis= -1) + sinogram = np.fft.ifftshift(sinogram, axes=-1) + sinogram_ft = np.fft.rfft(sinogram, axis=-1) - #grid generation + # grid generation with real points y_idx = np.fft.rfftfreq(n_rad) * np.pi * 2 n_real_points = len(y_idx) pts = np.empty((2, len(angles), n_real_points), dtype=self.dtype) pts[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] - + imgs = aspire.nufft.anufft( sinogram_ft.reshape(n_img, -1), pts.reshape(2, n_real_points * len(angles)), sz=(L, L), - real=True + real=True, ).reshape(n_img, L, L) - - return aspire.image.Image(imgs) - - def image_filter(self, filter_name, projections): - """ - Filter Method for projections. Will apply filter to line projection to get collection of projections (ramp, cosine, ... , etc.) - :param projections: Collection of line projections that need to be filtered. - :return: Filtered Projections. - """ - if projections is None: - raise ValueError('The input projections must not be None') - - filter_types = ('ramp', 'shepp-i logan', 'cosine', 'hamming', 'hann', None) - if filter_name is not filter_type: - raise ValueError(f"Unknown filter: {filter_name}") - # skimage filter - fourier_filter = _get_fourier_filter(projection_size_padded, filter_name) - projection = fft(img, axis=0) * fourier_filter - radon_filtered = np.real(ifft(projection, axis=0)[:img_shape, :]) - - """ - step 0: Look more into filter function from skimage - thoughts: apply filter to each point - """ - pass + return aspire.image.Image(imgs) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 56aa6776e0..6459145f52 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -4,6 +4,7 @@ from skimage.transform import radon from aspire.image import Image +from aspire.image.line import Line from aspire.utils import grid_2d # Relative tolerance comparing line projections to scikit @@ -136,3 +137,34 @@ def test_multidim(num_ang): reference_sinograms, axis=-1 ) np.testing.assert_array_less(_nrms, SK_TOL, "Error in image projections.") + + +def test_back_project_single(masked_image, num_ang): + """ + Test Line.backproject on a single stack of line projections or sinogram. Compares the reconstructed image to original image. + """ + # I'll be creating a sinogram representation of camera man + # reusing our grid fixture + # and testing the skimage's backwards project without a filter (note there currently is no filter so blurry + angles = np.linspace(0, 360, num_ang, endpoint=False) + rads = angles / 180 * np.pi + sinogram = Line(masked_image.project(rads)) + back_project = sinogram.back_project(num_ang) + + assert masked_img.shape == back_project.shape, "Shape must be the same." + + # no filter for now + sk_image_iradon = iradon(masked_image, theta=np.degrees(angles), filter_name=None) + + nrms = np.sqrt( + np.mean((sk_image_iradon - back_project) ** 2, axis=-1) + ) / np.linalg.norm(back_project, axis=-1) + np.testing.assert_array_less(nrms, SK_TOL, "Error in image reconstruction.") + + +def test_back_project_multidim(num_ang): + """ + Test Line.backproject on a stack of images. Extension of back_project_single but for multi-dimensional stacks. + """ + # assume once we can get this functioning for single stack + # we can get this working for multiple From 2f563d99250314899872026df6fc0d6fbf2cd519 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 23 Jul 2024 11:12:00 -0400 Subject: [PATCH 03/25] Stashing --- src/aspire/image/line.py | 33 +++++++++++++++------------------ 1 file changed, 15 insertions(+), 18 deletions(-) diff --git a/src/aspire/image/line.py b/src/aspire/image/line.py index 9f2b30c800..435c0aff5e 100644 --- a/src/aspire/image/line.py +++ b/src/aspire/image/line.py @@ -2,18 +2,16 @@ import os from pathlib import Path -import numpy as np - -import aspire from aspire.image import Image from aspire.nufft import anufft, nufft +import aspire +import numpy as np # noticed a lot of classes had these already, might be helpful for pathing, logging info, etc. (inc. os, logging) logger = logging.getLogger(__name__) - class Line: - def __init__(self, data, dtype=np.float64): + def __init__(self, data, dtype = np.float64): """ Initialize a Line Object. Change later (similar intuition from Image class) Question: Is it a line or collection of line object? @@ -27,11 +25,11 @@ def __init__(self, data, dtype=np.float64): if data.ndim == 2: data = data[np.newaxis, :, :] if data.ndim < 3: - raise ("Projection Dimensions should be more than Three-Dimensions") + raise('Projection Dimensions should be more than Three-Dimensions') self._data = data.astype(self.dtype, copy=False) self.stack_shape = self._data.shape[:-2] - self.n_images = self._data.shape[0] # broken for higher dimensional stacks - self.n_lines = self._data.shape[-2] + self.n_images = self._data.shape[0] #broken for higher dimensional stacks + self.n_lines = self._data.shape[-2] self.n_points = self._data.shape[-1] # self.n_dim = (self._data.shape[1], self._data.shape[2]) @@ -45,7 +43,7 @@ def stack(self): def back_project(self, angles): """ Back Projection Method for a single stack of lines. - + :param filter_name: string, optional Filter used in frequency domain filtering. Assign None to use no filter. :param angles: array @@ -54,26 +52,25 @@ def back_project(self, angles): """ sinogram = self._data n_img, n_angles, n_rad = sinogram.shape - assert n_angles == len( - angles - ), "Number of angles must match the number of projections" + assert n_angles == len(angles), "Number of angles must match the number of projections" L = n_rad - sinogram = np.fft.ifftshift(sinogram, axes=-1) - sinogram_ft = np.fft.rfft(sinogram, axis=-1) + sinogram = np.fft.ifftshift(sinogram, axes= -1) + sinogram_ft = np.fft.rfft(sinogram, axis= -1) - # grid generation with real points + #grid generation with real points y_idx = np.fft.rfftfreq(n_rad) * np.pi * 2 n_real_points = len(y_idx) pts = np.empty((2, len(angles), n_real_points), dtype=self.dtype) pts[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] - + imgs = aspire.nufft.anufft( sinogram_ft.reshape(n_img, -1), pts.reshape(2, n_real_points * len(angles)), sz=(L, L), - real=True, + real=True ).reshape(n_img, L, L) - + return aspire.image.Image(imgs) + From 09087868ad6a21fe05f723943b5b38d427684b00 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Thu, 25 Jul 2024 17:16:10 -0400 Subject: [PATCH 04/25] Fixed Scaling Issue with BackProject and Integrated NRMSE to One Stack Test --- src/aspire/image/line.py | 45 +++++++++++++------------ tests/test_sinogram.py | 72 +++++++++++++++++++++++++++++----------- 2 files changed, 77 insertions(+), 40 deletions(-) diff --git a/src/aspire/image/line.py b/src/aspire/image/line.py index 435c0aff5e..b68a84383d 100644 --- a/src/aspire/image/line.py +++ b/src/aspire/image/line.py @@ -1,17 +1,15 @@ import logging -import os -from pathlib import Path -from aspire.image import Image -from aspire.nufft import anufft, nufft -import aspire import numpy as np +import aspire + # noticed a lot of classes had these already, might be helpful for pathing, logging info, etc. (inc. os, logging) logger = logging.getLogger(__name__) + class Line: - def __init__(self, data, dtype = np.float64): + def __init__(self, data, dtype=np.float64): """ Initialize a Line Object. Change later (similar intuition from Image class) Question: Is it a line or collection of line object? @@ -25,13 +23,14 @@ def __init__(self, data, dtype = np.float64): if data.ndim == 2: data = data[np.newaxis, :, :] if data.ndim < 3: - raise('Projection Dimensions should be more than Three-Dimensions') + assert "Projection Dimensions should be more than Three-Dimensions" self._data = data.astype(self.dtype, copy=False) - self.stack_shape = self._data.shape[:-2] - self.n_images = self._data.shape[0] #broken for higher dimensional stacks - self.n_lines = self._data.shape[-2] - self.n_points = self._data.shape[-1] - # self.n_dim = (self._data.shape[1], self._data.shape[2]) + + # self.stack_shape = self._data.shape[:-2] + # self.n_images = self._data.shape[0] #broken for higher dimensional stacks + # self.n_lines = self._data.shape[-2] + # self.n_points = self._data.shape[-1] + # self.n_dim = (self._data.shape[1], self._data.shape[2]) def __str__(self): return f"Line(n_images = {self.n_images}, n_points = {self.n_points})" @@ -43,7 +42,7 @@ def stack(self): def back_project(self, angles): """ Back Projection Method for a single stack of lines. - + :param filter_name: string, optional Filter used in frequency domain filtering. Assign None to use no filter. :param angles: array @@ -52,25 +51,29 @@ def back_project(self, angles): """ sinogram = self._data n_img, n_angles, n_rad = sinogram.shape - assert n_angles == len(angles), "Number of angles must match the number of projections" + assert n_angles == len( + angles + ), "Number of angles must match the number of projections" L = n_rad - sinogram = np.fft.ifftshift(sinogram, axes= -1) - sinogram_ft = np.fft.rfft(sinogram, axis= -1) + sinogram = np.fft.ifftshift(sinogram, axes=-1) + sinogram_ft = np.fft.rfft(sinogram, axis=-1) - #grid generation with real points + # grid generation with real points y_idx = np.fft.rfftfreq(n_rad) * np.pi * 2 n_real_points = len(y_idx) pts = np.empty((2, len(angles), n_real_points), dtype=self.dtype) pts[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] - + imgs = aspire.nufft.anufft( sinogram_ft.reshape(n_img, -1), pts.reshape(2, n_real_points * len(angles)), sz=(L, L), - real=True + real=True, ).reshape(n_img, L, L) - + + # normalization which gives us roughly the same error regardless of angles + imgs = imgs / (n_real_points * len(angles)) + return aspire.image.Image(imgs) - diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 6459145f52..189a8266ae 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -1,15 +1,17 @@ import numpy as np import pytest from skimage import data -from skimage.transform import radon +from skimage.transform import iradon, radon from aspire.image import Image from aspire.image.line import Line from aspire.utils import grid_2d # Relative tolerance comparing line projections to scikit -# The same tolerance will be used in all scikit comparisons -SK_TOL = 0.005 +# The same tolerance will be used in all scikit forward and backward comparisons +SK_TOL_FORWARDPROJECT = 0.005 + +SK_TOL_BACKPROJECT = 0.2 IMG_SIZES = [ 511, @@ -95,7 +97,9 @@ def test_image_project(masked_image, num_ang): reference_sinogram, axis=-1 ) - np.testing.assert_array_less(nrms, SK_TOL, "Error in image projections.") + np.testing.assert_array_less( + nrms, SK_TOL_FORWARDPROJECT, "Error in image projections." + ) def test_multidim(num_ang): @@ -136,35 +140,65 @@ def test_multidim(num_ang): _nrms = np.sqrt(np.mean((s - reference_sinograms) ** 2, axis=-1)) / np.linalg.norm( reference_sinograms, axis=-1 ) - np.testing.assert_array_less(_nrms, SK_TOL, "Error in image projections.") + np.testing.assert_array_less( + _nrms, SK_TOL_FORWARDPROJECT, "Error in image projections." + ) def test_back_project_single(masked_image, num_ang): """ Test Line.backproject on a single stack of line projections or sinogram. Compares the reconstructed image to original image. """ - # I'll be creating a sinogram representation of camera man - # reusing our grid fixture - # and testing the skimage's backwards project without a filter (note there currently is no filter so blurry angles = np.linspace(0, 360, num_ang, endpoint=False) rads = angles / 180 * np.pi - sinogram = Line(masked_image.project(rads)) - back_project = sinogram.back_project(num_ang) + sinogram_np = masked_image.project(rads) + sinogram = Line(sinogram_np) + back_project = sinogram.back_project(rads) - assert masked_img.shape == back_project.shape, "Shape must be the same." + assert masked_image.shape == back_project.shape, "The shape must be the same." - # no filter for now - sk_image_iradon = iradon(masked_image, theta=np.degrees(angles), filter_name=None) + # generate circular mask w/ radius 1 to reconstructed image + # aim to remove discrepencies for the edges of the image + g = grid_2d(sinogram_np.shape[2], normalized=True, shifted=True) + mask = g["r"] < 1 + our_back_project = back_project.asnumpy()[0] * mask + + # generating sci-kit image backproject method w/ no filter + sk_image_iradon = iradon(sinogram_np[0].T, theta=angles[::-1], filter_name=None) - nrms = np.sqrt( - np.mean((sk_image_iradon - back_project) ** 2, axis=-1) - ) / np.linalg.norm(back_project, axis=-1) - np.testing.assert_array_less(nrms, SK_TOL, "Error in image reconstruction.") + # we apply a normalized root mean square error on the images to find relative error to range of ref. image + # Note: toleranc is typically < 0.2 regardless of angles, pixels, etc. + nrmse = np.sqrt(np.mean((our_back_project - sk_image_iradon) ** 2)) / ( + np.max(sk_image_iradon - np.min(sk_image_iradon)) + ) + assert ( + nrmse < SK_TOL_BACKPROJECT + ), f"NRMSE is too high: {nrmse}, expected less than {SK_TOL_BACKPROJECT}" def test_back_project_multidim(num_ang): """ Test Line.backproject on a stack of images. Extension of back_project_single but for multi-dimensional stacks. """ - # assume once we can get this functioning for single stack - # we can get this working for multiple + L = 512 # pixels + n = 3 + m = 2 + + # Generate a mask + g = grid_2d(L, normalized=True, shifted=True) + mask = g["r"] < 1 + + # Generate images + imgs = Image(np.random.random((m, n, L, L))) * mask + + # Generate line project angles + angles = np.linspace(0, 360, num_ang, endpoint=False) + rads = angles / 180.0 * np.pi + s = imgs.project(rads) + sinogram = Line(s) + return sinogram + # back project + grid + + # sci-kit back project + + # compare nrmse for all images in the stack From d374b81e07b89c7ffc85fb156d8e0f30575827f9 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Fri, 26 Jul 2024 11:11:13 -0400 Subject: [PATCH 05/25] fixed single back_project test --- tests/test_sinogram.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 189a8266ae..13b77c4a90 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -180,23 +180,11 @@ def test_back_project_multidim(num_ang): """ Test Line.backproject on a stack of images. Extension of back_project_single but for multi-dimensional stacks. """ - L = 512 # pixels - n = 3 - m = 2 - # Generate a mask - g = grid_2d(L, normalized=True, shifted=True) - mask = g["r"] < 1 # Generate images - imgs = Image(np.random.random((m, n, L, L))) * mask # Generate line project angles - angles = np.linspace(0, 360, num_ang, endpoint=False) - rads = angles / 180.0 * np.pi - s = imgs.project(rads) - sinogram = Line(s) - return sinogram # back project + grid # sci-kit back project From 974b0392934080881b427bba35f9fe5c166828a7 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 26 Jul 2024 11:54:31 -0400 Subject: [PATCH 06/25] reorg Line to avoid circ import. Interop Image/Line classes --- src/aspire/image/image.py | 3 ++- src/aspire/line/__init__.py | 1 + src/aspire/{image => line}/line.py | 7 +++---- tests/test_sinogram.py | 2 +- 4 files changed, 7 insertions(+), 6 deletions(-) create mode 100644 src/aspire/line/__init__.py rename src/aspire/{image => line}/line.py (94%) diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py index 1bbe2ba160..263eaa1780 100644 --- a/src/aspire/image/image.py +++ b/src/aspire/image/image.py @@ -8,6 +8,7 @@ from PIL import Image as PILImage from scipy.linalg import lstsq +import aspire.line import aspire.volume from aspire.nufft import anufft, nufft from aspire.numeric import fft, xp @@ -238,7 +239,7 @@ def project(self, angles): # Radon transform, output: (stack size, angles, points) image_rt = np.fft.fftshift(np.fft.irfft(image_ft, n=n_points, axis=-1), axes=-1) image_rt = image_rt.reshape(*original_stack, n_angles, n_points) - return image_rt + return aspire.line.Line(image_rt) @property def res(self): diff --git a/src/aspire/line/__init__.py b/src/aspire/line/__init__.py new file mode 100644 index 0000000000..d3856d67ad --- /dev/null +++ b/src/aspire/line/__init__.py @@ -0,0 +1 @@ +from .line import Line diff --git a/src/aspire/image/line.py b/src/aspire/line/line.py similarity index 94% rename from src/aspire/image/line.py rename to src/aspire/line/line.py index b68a84383d..4b947592ed 100644 --- a/src/aspire/image/line.py +++ b/src/aspire/line/line.py @@ -2,9 +2,9 @@ import numpy as np -import aspire +import aspire.image +from aspire.nufft import anufft -# noticed a lot of classes had these already, might be helpful for pathing, logging info, etc. (inc. os, logging) logger = logging.getLogger(__name__) @@ -66,7 +66,7 @@ def back_project(self, angles): pts[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] - imgs = aspire.nufft.anufft( + imgs = anufft( sinogram_ft.reshape(n_img, -1), pts.reshape(2, n_real_points * len(angles)), sz=(L, L), @@ -75,5 +75,4 @@ def back_project(self, angles): # normalization which gives us roughly the same error regardless of angles imgs = imgs / (n_real_points * len(angles)) - return aspire.image.Image(imgs) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 13b77c4a90..41259ea563 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -4,7 +4,7 @@ from skimage.transform import iradon, radon from aspire.image import Image -from aspire.image.line import Line +from aspire.line import Line from aspire.utils import grid_2d # Relative tolerance comparing line projections to scikit From f3fe20c0dcd57bf1f1cbaa84c46bf401beb49dca Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 26 Jul 2024 12:10:50 -0400 Subject: [PATCH 07/25] adjust tests towards Line/Image interop [skip ci] Co-authored-by: Marc Karimi --- tests/test_sinogram.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 41259ea563..f410940180 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -151,15 +151,15 @@ def test_back_project_single(masked_image, num_ang): """ angles = np.linspace(0, 360, num_ang, endpoint=False) rads = angles / 180 * np.pi - sinogram_np = masked_image.project(rads) - sinogram = Line(sinogram_np) + sinogram = masked_image.project(rads) + sinogram_np = sinogram.asnumpy() back_project = sinogram.back_project(rads) assert masked_image.shape == back_project.shape, "The shape must be the same." # generate circular mask w/ radius 1 to reconstructed image # aim to remove discrepencies for the edges of the image - g = grid_2d(sinogram_np.shape[2], normalized=True, shifted=True) + g = grid_2d(back_project.resolution, normalized=True, shifted=True) mask = g["r"] < 1 our_back_project = back_project.asnumpy()[0] * mask From 46af4a166110238f6f6523be757b02cab488b0a3 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Fri, 26 Jul 2024 15:50:31 -0400 Subject: [PATCH 08/25] passing the 20/20 test cases and added Attributes + Methods to the Line class [need fix] --- src/aspire/line/line.py | 118 +++++++++++++++++++++++++++++++++++++--- tests/test_sinogram.py | 6 +- 2 files changed, 114 insertions(+), 10 deletions(-) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index 4b947592ed..6a19f03b94 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -1,5 +1,8 @@ import logging +import os +from warnings import catch_warnings, filterwarnings, simplefilter, warn +import mrcfile import numpy as np import aspire.image @@ -12,7 +15,6 @@ class Line: def __init__(self, data, dtype=np.float64): """ Initialize a Line Object. Change later (similar intuition from Image class) - Question: Is it a line or collection of line object? :param data: Numpy array containing image data with shape `(..., resolution, resolution)`. @@ -25,15 +27,117 @@ def __init__(self, data, dtype=np.float64): if data.ndim < 3: assert "Projection Dimensions should be more than Three-Dimensions" self._data = data.astype(self.dtype, copy=False) + self.ndim = self._data.ndim + self.shape = self._data.shape + self.stack_shape = self._data.shape[:-2] + self.stack_n_dim = self._data.ndim - 2 # fix + self.n = np.product(self.stack_shape) # stack number + self.n_angles = self._data.shape[-1] # fix + self.n_radial_points = self._data.shape[-1] - # self.stack_shape = self._data.shape[:-2] - # self.n_images = self._data.shape[0] #broken for higher dimensional stacks - # self.n_lines = self._data.shape[-2] - # self.n_points = self._data.shape[-1] - # self.n_dim = (self._data.shape[1], self._data.shape[2]) + # Numpy interop + # https://numpy.org/devdocs/user/basics.interoperability.html#the-array-interface-protocol + self.__array_interface__ = self._data.__array_interface__ + self.__array__ = self._data + + def _check_key_dims(self, key): + if isinstance(key, tuple) and (len(key) > self._data.ndim): + raise ValueError( + f"Line stack_dim is {self.stack_n_dim}, slice length must be =< {self.n_dim}" + ) + + def __getitem__(self, key): + self._check_key_dims(key) + return self.__class__(self._data[key]) + + def __setitem__(self, key, value): + self._check_key_dims(key) + self._data[key] = value + + def stack_reshape(self, *args): + """ + Reshape the stack axis. + + :*args: Integer(s) or tuple describing the intended shape. + + :returns: Line instance + """ + + # If we're passed a tuple, use that + if len(args) == 1 and isinstance(args[0], tuple): + shape = args[0] + else: + # Otherwise use the variadic args + shape = args + + # Sanity check the size + if shape != (-1,) and np.prod(shape) != self.n: + raise ValueError( + f"Number of sinogram images {self.n_images} cannot be reshaped to {shape}." + ) + + return self.__class__(self._data.reshape(*shape, *self._data.shape[-2:])) + + def asnumpy(self): + """ + Return image data as a (, angles, radians) + read-only array view. + + :return: read-only ndarray view + """ + + view = self._data.view() + view.flags.writeable = False + return view + + def copy(self): + return self.__class__(self._data.copy()) + + # fix later + def save(self, mrcs_filepath, overwrite=False): + if self.stack_ndim > 1: + raise NotImplementedError("`save` is currently limited to 1D image stacks.") + + with mrcfile.new(mrcs_filepath, overwrite=overwrite) as mrc: + # original input format (the image index first) + mrc.set_data(self._data.astype(np.float32)) + + # fix later + @staticmethod + def load(filepath, dtype=None): + """ + Load raw data from supported files. + + Currently MRC and TIFF are supported. + + :param filepath: File path (string). + :param dtype: Optionally force cast to `dtype`. + Default dtype is inferred from the file contents. + :return: numpy array of image data. + """ + + # Get the file extension + ext = os.path.splitext(filepath)[1] + + # On unsupported extension, raise with suggested file types + if ext not in Image.extensions: + raise RuntimeError( + f"Attempting to open unsupported file extension '{ext}', try {list(Image.extensions.keys())}." + ) + + # Call the appropriate file reader + im = Image.extensions[ext](filepath) + + # Attempt casting when user provides dtype + if dtype is not None: + im = im.astype(dtype, copy=False) + + # Return as Image instance + return Image(im) def __str__(self): - return f"Line(n_images = {self.n_images}, n_points = {self.n_points})" + # fix later + return f"Line(n_images = {self.n}, n_angles = {self.n_points}, n_radial_points = {self.n_radial_points})" @property def stack(self): diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index f410940180..e8dbf797df 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -93,9 +93,9 @@ def test_image_project(masked_image, num_ang): assert reference_sinogram.shape == (len(angles), ny), "Incorrect Shape" # compare project method on ski-image reference - nrms = np.sqrt(np.mean((s[0] - reference_sinogram) ** 2, axis=-1)) / np.linalg.norm( - reference_sinogram, axis=-1 - ) + nrms = np.sqrt( + np.mean((s[0]._data - reference_sinogram) ** 2, axis=-1) + ) / np.linalg.norm(reference_sinogram, axis=-1) np.testing.assert_array_less( nrms, SK_TOL_FORWARDPROJECT, "Error in image projections." From 9af2d836b12d25dff7225d39e85f2e32e1138401 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 30 Jul 2024 05:04:39 -0400 Subject: [PATCH 09/25] finished multidim test --- src/aspire/line/line.py | 67 ++++++++--------------------------------- tests/test_sinogram.py | 49 ++++++++++++++++++++++++------ 2 files changed, 52 insertions(+), 64 deletions(-) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index 6a19f03b94..58cdad88a3 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -30,9 +30,9 @@ def __init__(self, data, dtype=np.float64): self.ndim = self._data.ndim self.shape = self._data.shape self.stack_shape = self._data.shape[:-2] - self.stack_n_dim = self._data.ndim - 2 # fix - self.n = np.product(self.stack_shape) # stack number - self.n_angles = self._data.shape[-1] # fix + self.stack_n_dim = self._data.ndim - 2 + self.n = np.product(self.stack_shape) + self.n_angles = self._data.shape[-2] self.n_radial_points = self._data.shape[-1] # Numpy interop @@ -93,50 +93,7 @@ def asnumpy(self): def copy(self): return self.__class__(self._data.copy()) - # fix later - def save(self, mrcs_filepath, overwrite=False): - if self.stack_ndim > 1: - raise NotImplementedError("`save` is currently limited to 1D image stacks.") - - with mrcfile.new(mrcs_filepath, overwrite=overwrite) as mrc: - # original input format (the image index first) - mrc.set_data(self._data.astype(np.float32)) - - # fix later - @staticmethod - def load(filepath, dtype=None): - """ - Load raw data from supported files. - - Currently MRC and TIFF are supported. - - :param filepath: File path (string). - :param dtype: Optionally force cast to `dtype`. - Default dtype is inferred from the file contents. - :return: numpy array of image data. - """ - - # Get the file extension - ext = os.path.splitext(filepath)[1] - - # On unsupported extension, raise with suggested file types - if ext not in Image.extensions: - raise RuntimeError( - f"Attempting to open unsupported file extension '{ext}', try {list(Image.extensions.keys())}." - ) - - # Call the appropriate file reader - im = Image.extensions[ext](filepath) - - # Attempt casting when user provides dtype - if dtype is not None: - im = im.astype(dtype, copy=False) - - # Return as Image instance - return Image(im) - def __str__(self): - # fix later return f"Line(n_images = {self.n}, n_angles = {self.n_points}, n_radial_points = {self.n_radial_points})" @property @@ -153,30 +110,30 @@ def back_project(self, angles): assuming not perfectly radial angles :return: stack of reconstructed """ - sinogram = self._data - n_img, n_angles, n_rad = sinogram.shape - assert n_angles == len( - angles + assert ( + len(angles) == self.n_angles ), "Number of angles must match the number of projections" - L = n_rad + original_stack_shape = self.stack_shape + sinogram = self.stack_reshape(-1) + L = self.n_radial_points sinogram = np.fft.ifftshift(sinogram, axes=-1) sinogram_ft = np.fft.rfft(sinogram, axis=-1) # grid generation with real points - y_idx = np.fft.rfftfreq(n_rad) * np.pi * 2 + y_idx = np.fft.rfftfreq(self.n_radial_points) * np.pi * 2 n_real_points = len(y_idx) pts = np.empty((2, len(angles), n_real_points), dtype=self.dtype) pts[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] imgs = anufft( - sinogram_ft.reshape(n_img, -1), + sinogram_ft.reshape(self.n, -1), pts.reshape(2, n_real_points * len(angles)), sz=(L, L), real=True, - ).reshape(n_img, L, L) + ).reshape(self.n, L, L) # normalization which gives us roughly the same error regardless of angles imgs = imgs / (n_real_points * len(angles)) - return aspire.image.Image(imgs) + return aspire.image.Image(imgs).stack_reshape(original_stack_shape) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index e8dbf797df..f9a24b1590 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -71,7 +71,7 @@ def masked_image(dtype, img_size): # Image.project and compare results to skimage.radon -def test_image_project(masked_image, num_ang): +def test_project_single(masked_image, num_ang): """ Test Image.project on a single stack of images. Compares project method output with skimage project. """ @@ -102,7 +102,7 @@ def test_image_project(masked_image, num_ang): ) -def test_multidim(num_ang): +def test_project_multidim(num_ang): """ Test Image.project on stacks of images. Extension of test_image_project but for multi-dimensional stacks. """ @@ -167,7 +167,7 @@ def test_back_project_single(masked_image, num_ang): sk_image_iradon = iradon(sinogram_np[0].T, theta=angles[::-1], filter_name=None) # we apply a normalized root mean square error on the images to find relative error to range of ref. image - # Note: toleranc is typically < 0.2 regardless of angles, pixels, etc. + # Note: tolerance is typically < 0.2 regardless of angles, pixels, etc. nrmse = np.sqrt(np.mean((our_back_project - sk_image_iradon) ** 2)) / ( np.max(sk_image_iradon - np.min(sk_image_iradon)) ) @@ -178,15 +178,46 @@ def test_back_project_single(masked_image, num_ang): def test_back_project_multidim(num_ang): """ - Test Line.backproject on a stack of images. Extension of back_project_single but for multi-dimensional stacks. + Test Line.backproject on a stack of images. Extension of back_project_single but for multi-dimensional stacks. Similar to forward_multidim test. """ - # Generate a mask + L = 512 # pixels + n = 3 + m = 2 + + g = grid_2d(L, normalized=True, shifted=True) + mask = g["r"] < 1 # Generate images + imgs = Image(np.random.random((m, n, L, L))) * mask + angles = np.linspace(0, 360, num_ang, endpoint=False) + rads = angles / 180 * np.pi - # Generate line project angles - # back project + grid + # apply a forward project on the image, then backwards + ours_forward = imgs.project(rads) + ours_backward = ours_forward.back_project(rads) - # sci-kit back project + # Compare + reference_back_projects = np.empty((m, n, L, L)) + for i in range(m): + for j in range(n): + img = imgs[i, j] + # Compute the singleton case, and compare with stack. + single_sinogram = img.project(rads) + back_project = single_sinogram.back_project(rads) + + # These should be allclose up to determinism. + np.testing.assert_allclose(ours_backward[i, j : j + 1], back_project[0]) + + # Next individually compute sk's iradon transform for each image. + reference_back_projects[i, j] = iradon( + single_sinogram.asnumpy()[0].T, theta=angles[::-1], filter_name=None + ) - # compare nrmse for all images in the stack + # apply a mask, then find the NRMSE on the collection of images + # similar tolerance level to single project test + nrmse = np.sqrt( + np.mean((ours_backward.asnumpy() * mask - reference_back_projects) ** 2) + ) / (np.max(reference_back_projects) - np.min(reference_back_projects)) + assert ( + nrmse < SK_TOL_BACKPROJECT + ), f"NRMSE is too high for image ({i},{j}): {nrmse}, expected less than {SK_TOL_BACKPROJECT}" From 286b8fe66a4c341e3b458c0210824a4d4a02cf31 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 30 Jul 2024 11:54:56 -0400 Subject: [PATCH 10/25] removed unused statements --- src/aspire/line/line.py | 3 --- tests/test_sinogram.py | 1 - 2 files changed, 4 deletions(-) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index 58cdad88a3..f9c688f71d 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -1,8 +1,5 @@ import logging -import os -from warnings import catch_warnings, filterwarnings, simplefilter, warn -import mrcfile import numpy as np import aspire.image diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index f9a24b1590..9b27b10dca 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -4,7 +4,6 @@ from skimage.transform import iradon, radon from aspire.image import Image -from aspire.line import Line from aspire.utils import grid_2d # Relative tolerance comparing line projections to scikit From 3ac2b28db9237834d8198ce882ffd122e7e3689c Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 30 Jul 2024 12:18:31 -0400 Subject: [PATCH 11/25] initial fft changes --- src/aspire/image/image.py | 14 +++++++------- src/aspire/numeric/cupy_fft.py | 12 ++++++++++++ src/aspire/numeric/scipy_fft.py | 9 +++++++++ 3 files changed, 28 insertions(+), 7 deletions(-) diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py index 263eaa1780..b589471577 100644 --- a/src/aspire/image/image.py +++ b/src/aspire/image/image.py @@ -222,24 +222,24 @@ def project(self, angles): original_stack = self.stack_shape # 2-D grid - radial_idx = np.fft.rfftfreq(n_points) * np.pi * 2 + radial_idx = fft.rfftfreq(n_points) * xp.pi * 2 n_real_points = len(radial_idx) n_angles = len(angles) - pts = np.empty((2, n_angles, n_real_points), dtype=self.dtype) - pts[0] = radial_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] - pts[1] = radial_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] + pts = xp.empty((2, n_angles, n_real_points), dtype=self.dtype) + pts[0] = radial_idx[xp.newaxis, :] * xp.sin(angles)[:, xp.newaxis] + pts[1] = radial_idx[xp.newaxis, :] * xp.cos(angles)[:, xp.newaxis] pts = pts.reshape(2, n_real_points * n_angles) # compute the polar nufft (NUFFT) - image_ft = nufft(self.stack_reshape(-1)._data, pts).reshape( + image_ft = nufft(xp.asarray(self.stack_reshape(-1)._data), pts).reshape( self.n_images, n_angles, n_real_points ) # Radon transform, output: (stack size, angles, points) - image_rt = np.fft.fftshift(np.fft.irfft(image_ft, n=n_points, axis=-1), axes=-1) + image_rt = fft.fftshift(fft.irfft(image_ft, n=n_points, axis=-1), axes=-1) image_rt = image_rt.reshape(*original_stack, n_angles, n_points) - return aspire.line.Line(image_rt) + return aspire.line.Line(xp.asnumpy(image_rt)) @property def res(self): diff --git a/src/aspire/numeric/cupy_fft.py b/src/aspire/numeric/cupy_fft.py index 6ad6a4e9da..ce537a1cba 100644 --- a/src/aspire/numeric/cupy_fft.py +++ b/src/aspire/numeric/cupy_fft.py @@ -100,3 +100,15 @@ def dct(self, x, **kwargs): @_preserve_host def idct(self, x, **kwargs): return cufft.idct(x, **kwargs) + + @_preserve_host + def rfftfreq(self, x, **kwargs): + return cufft.rfftfreq(x, **kwargs) + + @_preserve_host + def irfft(self, x, **kwargs): + return cufft.irfft(x, **kwargs) + + @_preserve_host + def rfft(self, x, **kwargs): + return cufft.rfft(x, **kwargs) diff --git a/src/aspire/numeric/scipy_fft.py b/src/aspire/numeric/scipy_fft.py index 3891d45671..0ef5c95f16 100644 --- a/src/aspire/numeric/scipy_fft.py +++ b/src/aspire/numeric/scipy_fft.py @@ -39,3 +39,12 @@ def dct(self, x, **kwargs): def idct(self, x, **kwargs): return sp.fft.idct(x, **kwargs) + + def rfftfreq(self, x, **kwargs): + return sp.fft.rfftfreq(x, **kwargs) + + def irfft(self, x, **kwargs): + return sp.fft.irfft(x, **kwargs) + + def rfft(self, x, **kwargs): + return sp.fft.rfft(x, **kwargs) From ac6d747b15a1f2fa62bbebc5e231e637a4f7fd24 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 30 Jul 2024 12:38:37 -0400 Subject: [PATCH 12/25] stashing gpu fixes --- src/aspire/image/image.py | 1 + src/aspire/numeric/cupy_fft.py | 5 ++--- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py index b589471577..09a7dc56cd 100644 --- a/src/aspire/image/image.py +++ b/src/aspire/image/image.py @@ -225,6 +225,7 @@ def project(self, angles): radial_idx = fft.rfftfreq(n_points) * xp.pi * 2 n_real_points = len(radial_idx) n_angles = len(angles) + angles = xp.asarray(angles) pts = xp.empty((2, n_angles, n_real_points), dtype=self.dtype) pts[0] = radial_idx[xp.newaxis, :] * xp.sin(angles)[:, xp.newaxis] diff --git a/src/aspire/numeric/cupy_fft.py b/src/aspire/numeric/cupy_fft.py index ce537a1cba..b491a0dcd1 100644 --- a/src/aspire/numeric/cupy_fft.py +++ b/src/aspire/numeric/cupy_fft.py @@ -101,9 +101,8 @@ def dct(self, x, **kwargs): def idct(self, x, **kwargs): return cufft.idct(x, **kwargs) - @_preserve_host - def rfftfreq(self, x, **kwargs): - return cufft.rfftfreq(x, **kwargs) + def rfftfreq(self, n, **kwargs): + return cufft.rfftfreq(n, **kwargs) @_preserve_host def irfft(self, x, **kwargs): From 3d2a3a92179f03e38f7e820f28943a84b47789e3 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 30 Jul 2024 13:41:05 -0400 Subject: [PATCH 13/25] forward gpu --- src/aspire/config_default.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aspire/config_default.yaml b/src/aspire/config_default.yaml index def78983c0..fed4cea50a 100644 --- a/src/aspire/config_default.yaml +++ b/src/aspire/config_default.yaml @@ -1,9 +1,9 @@ version: 0.12.3 common: # numeric module to use - one of numpy/cupy - numeric: numpy + numeric: cupy # fft backend to use - one of pyfftw/scipy/cupy/mkl - fft: scipy + fft: cupy # Set cache directory for ASPIRE example data. # By default the cache location will be set by pooch.os_cache(), From ad1208104e35407e80e518c840fb6becd30c4cd5 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Fri, 2 Aug 2024 00:52:32 -0400 Subject: [PATCH 14/25] changed backproject to run on gpu (cupy) --- src/aspire/line/line.py | 13 +++++++------ tests/test_sinogram.py | 15 ++++++++++----- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index f9c688f71d..fb1b180a6b 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -4,6 +4,7 @@ import aspire.image from aspire.nufft import anufft +from aspire.numeric import fft, xp logger = logging.getLogger(__name__) @@ -114,15 +115,15 @@ def back_project(self, angles): original_stack_shape = self.stack_shape sinogram = self.stack_reshape(-1) L = self.n_radial_points - sinogram = np.fft.ifftshift(sinogram, axes=-1) - sinogram_ft = np.fft.rfft(sinogram, axis=-1) + sinogram = fft.ifftshift(sinogram, axes=-1) + sinogram_ft = fft.rfft(sinogram, axis=-1) # grid generation with real points - y_idx = np.fft.rfftfreq(self.n_radial_points) * np.pi * 2 + y_idx = xp.fft.rfftfreq(self.n_radial_points) * xp.pi * 2 n_real_points = len(y_idx) - pts = np.empty((2, len(angles), n_real_points), dtype=self.dtype) - pts[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis] - pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis] + pts = xp.empty((2, len(angles), n_real_points), dtype=self.dtype) + pts[0] = y_idx[xp.newaxis, :] * xp.sin(angles)[:, xp.newaxis] + pts[1] = y_idx[xp.newaxis, :] * xp.cos(angles)[:, xp.newaxis] imgs = anufft( sinogram_ft.reshape(self.n, -1), diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 9b27b10dca..d186866012 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -4,6 +4,7 @@ from skimage.transform import iradon, radon from aspire.image import Image +from aspire.numeric import xp from aspire.utils import grid_2d # Relative tolerance comparing line projections to scikit @@ -146,9 +147,12 @@ def test_project_multidim(num_ang): def test_back_project_single(masked_image, num_ang): """ - Test Line.backproject on a single stack of line projections or sinogram. Compares the reconstructed image to original image. + Test Line.backproject on a single stack of line projections or sinogram. Compares the reconstructed image to original. """ - angles = np.linspace(0, 360, num_ang, endpoint=False) + angles = xp.asarray(np.linspace(0, 360, num_ang, endpoint=False)) + angles_np = xp.asnumpy( + angles + ) # skimage requires numpy array while we're using cupy arrays rads = angles / 180 * np.pi sinogram = masked_image.project(rads) sinogram_np = sinogram.asnumpy() @@ -163,7 +167,7 @@ def test_back_project_single(masked_image, num_ang): our_back_project = back_project.asnumpy()[0] * mask # generating sci-kit image backproject method w/ no filter - sk_image_iradon = iradon(sinogram_np[0].T, theta=angles[::-1], filter_name=None) + sk_image_iradon = iradon(sinogram_np[0].T, theta=angles_np[::-1], filter_name=None) # we apply a normalized root mean square error on the images to find relative error to range of ref. image # Note: tolerance is typically < 0.2 regardless of angles, pixels, etc. @@ -188,7 +192,8 @@ def test_back_project_multidim(num_ang): # Generate images imgs = Image(np.random.random((m, n, L, L))) * mask - angles = np.linspace(0, 360, num_ang, endpoint=False) + angles = xp.asarray(np.linspace(0, 360, num_ang, endpoint=False)) + angles_np = xp.asnumpy(angles) # need this for the skimage transformations rads = angles / 180 * np.pi # apply a forward project on the image, then backwards @@ -209,7 +214,7 @@ def test_back_project_multidim(num_ang): # Next individually compute sk's iradon transform for each image. reference_back_projects[i, j] = iradon( - single_sinogram.asnumpy()[0].T, theta=angles[::-1], filter_name=None + single_sinogram.asnumpy()[0].T, theta=angles_np[::-1], filter_name=None ) # apply a mask, then find the NRMSE on the collection of images From d5464a52e13e94d45282ce50502dfa21bcb86891 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Fri, 2 Aug 2024 11:20:04 -0400 Subject: [PATCH 15/25] revert config --- src/aspire/config_default.yaml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/aspire/config_default.yaml b/src/aspire/config_default.yaml index fed4cea50a..def78983c0 100644 --- a/src/aspire/config_default.yaml +++ b/src/aspire/config_default.yaml @@ -1,9 +1,9 @@ version: 0.12.3 common: # numeric module to use - one of numpy/cupy - numeric: cupy + numeric: numpy # fft backend to use - one of pyfftw/scipy/cupy/mkl - fft: cupy + fft: scipy # Set cache directory for ASPIRE example data. # By default the cache location will be set by pooch.os_cache(), From ac201e7d8464f04f39a9a529c167b0742339f057 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Fri, 2 Aug 2024 12:47:05 -0400 Subject: [PATCH 16/25] fixed gpu issues --- src/aspire/line/line.py | 6 +++--- tests/test_sinogram.py | 13 ++++--------- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index fb1b180a6b..7787367c3c 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -113,13 +113,13 @@ def back_project(self, angles): ), "Number of angles must match the number of projections" original_stack_shape = self.stack_shape - sinogram = self.stack_reshape(-1) + sinogram = xp.asarray(self.stack_reshape(-1)._data) L = self.n_radial_points sinogram = fft.ifftshift(sinogram, axes=-1) sinogram_ft = fft.rfft(sinogram, axis=-1) # grid generation with real points - y_idx = xp.fft.rfftfreq(self.n_radial_points) * xp.pi * 2 + y_idx = fft.rfftfreq(self.n_radial_points) * xp.pi * 2 n_real_points = len(y_idx) pts = xp.empty((2, len(angles), n_real_points), dtype=self.dtype) pts[0] = y_idx[xp.newaxis, :] * xp.sin(angles)[:, xp.newaxis] @@ -134,4 +134,4 @@ def back_project(self, angles): # normalization which gives us roughly the same error regardless of angles imgs = imgs / (n_real_points * len(angles)) - return aspire.image.Image(imgs).stack_reshape(original_stack_shape) + return aspire.image.Image(xp.asnumpy(imgs)).stack_reshape(original_stack_shape) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index d186866012..ddef7c480f 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -4,7 +4,6 @@ from skimage.transform import iradon, radon from aspire.image import Image -from aspire.numeric import xp from aspire.utils import grid_2d # Relative tolerance comparing line projections to scikit @@ -149,10 +148,7 @@ def test_back_project_single(masked_image, num_ang): """ Test Line.backproject on a single stack of line projections or sinogram. Compares the reconstructed image to original. """ - angles = xp.asarray(np.linspace(0, 360, num_ang, endpoint=False)) - angles_np = xp.asnumpy( - angles - ) # skimage requires numpy array while we're using cupy arrays + angles = np.linspace(0, 360, num_ang, endpoint=False) rads = angles / 180 * np.pi sinogram = masked_image.project(rads) sinogram_np = sinogram.asnumpy() @@ -167,7 +163,7 @@ def test_back_project_single(masked_image, num_ang): our_back_project = back_project.asnumpy()[0] * mask # generating sci-kit image backproject method w/ no filter - sk_image_iradon = iradon(sinogram_np[0].T, theta=angles_np[::-1], filter_name=None) + sk_image_iradon = iradon(sinogram_np[0].T, theta=angles[::-1], filter_name=None) # we apply a normalized root mean square error on the images to find relative error to range of ref. image # Note: tolerance is typically < 0.2 regardless of angles, pixels, etc. @@ -192,8 +188,7 @@ def test_back_project_multidim(num_ang): # Generate images imgs = Image(np.random.random((m, n, L, L))) * mask - angles = xp.asarray(np.linspace(0, 360, num_ang, endpoint=False)) - angles_np = xp.asnumpy(angles) # need this for the skimage transformations + angles = np.linspace(0, 360, num_ang, endpoint=False) rads = angles / 180 * np.pi # apply a forward project on the image, then backwards @@ -214,7 +209,7 @@ def test_back_project_multidim(num_ang): # Next individually compute sk's iradon transform for each image. reference_back_projects[i, j] = iradon( - single_sinogram.asnumpy()[0].T, theta=angles_np[::-1], filter_name=None + single_sinogram.asnumpy()[0].T, theta=angles[::-1], filter_name=None ) # apply a mask, then find the NRMSE on the collection of images From a15df2fe765696021babd7f1d875d84f0e738e06 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Fri, 2 Aug 2024 17:17:39 -0400 Subject: [PATCH 17/25] added angles array --- src/aspire/line/line.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index 7787367c3c..4b57867a04 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -117,6 +117,7 @@ def back_project(self, angles): L = self.n_radial_points sinogram = fft.ifftshift(sinogram, axes=-1) sinogram_ft = fft.rfft(sinogram, axis=-1) + angles = xp.asarray(angles) # grid generation with real points y_idx = fft.rfftfreq(self.n_radial_points) * xp.pi * 2 From e9429366fef8b472459ccd87a66a6f626f58fccc Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Tue, 6 Aug 2024 15:48:53 -0400 Subject: [PATCH 18/25] Pr and naming changes --- src/aspire/line/line.py | 34 +++++++++++++++++----------------- tests/test_sinogram.py | 32 +++++++++++++++++++------------- 2 files changed, 36 insertions(+), 30 deletions(-) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index 4b57867a04..8efee5b679 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -12,10 +12,14 @@ class Line: def __init__(self, data, dtype=np.float64): """ - Initialize a Line Object. Change later (similar intuition from Image class) + Initialize a Line Object. This is a stack of one or more line projections or sinograms. + + The stack can be multidimensional with 'n' equal to the product + of the stack dimensions. Singletons will be expanded into a stack + with one entry. :param data: Numpy array containing image data with shape - `(..., resolution, resolution)`. + `(..., angles, radial points)`. :param dtype: Optionally cast `data` to this dtype. Defaults to `data.dtype`. """ @@ -23,7 +27,9 @@ def __init__(self, data, dtype=np.float64): if data.ndim == 2: data = data[np.newaxis, :, :] if data.ndim < 3: - assert "Projection Dimensions should be more than Three-Dimensions" + raise ValueError( + f"Invalid data shape: {data.shape}. Expected shape: (..., angles, radial_points), where '...' is the stack number." + ) self._data = data.astype(self.dtype, copy=False) self.ndim = self._data.ndim self.shape = self._data.shape @@ -71,7 +77,7 @@ def stack_reshape(self, *args): # Sanity check the size if shape != (-1,) and np.prod(shape) != self.n: raise ValueError( - f"Number of sinogram images {self.n_images} cannot be reshaped to {shape}." + f"Number of sinogram images {self.n} cannot be reshaped to {shape}." ) return self.__class__(self._data.reshape(*shape, *self._data.shape[-2:])) @@ -94,23 +100,17 @@ def copy(self): def __str__(self): return f"Line(n_images = {self.n}, n_angles = {self.n_points}, n_radial_points = {self.n_radial_points})" - @property - def stack(self): - return self.n_images - - def back_project(self, angles): + def backproject(self, angles): """ Back Projection Method for a single stack of lines. - :param filter_name: string, optional - Filter used in frequency domain filtering. Assign None to use no filter. - :param angles: array - assuming not perfectly radial angles - :return: stack of reconstructed + :param angles: np.ndarray + 1D array of angles in radians. Each entry in the array corresponds to a different number of angles which are used to reconstruct the image. + :return: Aspire Image + An Image object containing the original stack size with a newly reconstructed numpy array of the images. Expected return shape should be (n_images, n_angles, n_radial_points) """ - assert ( - len(angles) == self.n_angles - ), "Number of angles must match the number of projections" + if len(angles) != self.n_angles: + raise ValueError("Number of angles must match the number of projections.") original_stack_shape = self.stack_shape sinogram = xp.asarray(self.stack_reshape(-1)._data) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index ddef7c480f..ed729d0b6f 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -144,15 +144,18 @@ def test_project_multidim(num_ang): ) -def test_back_project_single(masked_image, num_ang): +def test_backproject_single(masked_image, num_ang): """ - Test Line.backproject on a single stack of line projections or sinogram. Compares the reconstructed image to original. + Test Line.backproject on a single stack of line projections (sinograms). + + This test compares the reconstructed image from the `backproject` method to + the skimage method `iradon.` """ angles = np.linspace(0, 360, num_ang, endpoint=False) rads = angles / 180 * np.pi sinogram = masked_image.project(rads) sinogram_np = sinogram.asnumpy() - back_project = sinogram.back_project(rads) + back_project = sinogram.backproject(rads) assert masked_image.shape == back_project.shape, "The shape must be the same." @@ -168,16 +171,18 @@ def test_back_project_single(masked_image, num_ang): # we apply a normalized root mean square error on the images to find relative error to range of ref. image # Note: tolerance is typically < 0.2 regardless of angles, pixels, etc. nrmse = np.sqrt(np.mean((our_back_project - sk_image_iradon) ** 2)) / ( - np.max(sk_image_iradon - np.min(sk_image_iradon)) + np.max(sk_image_iradon) - np.min(sk_image_iradon) ) assert ( nrmse < SK_TOL_BACKPROJECT ), f"NRMSE is too high: {nrmse}, expected less than {SK_TOL_BACKPROJECT}" -def test_back_project_multidim(num_ang): +def test_backproject_multidim(num_ang): """ - Test Line.backproject on a stack of images. Extension of back_project_single but for multi-dimensional stacks. Similar to forward_multidim test. + Test Line.backproject on a stack of line projections. + + Extension of the `backproject_single` test but checks for multi-dimensional stacks. """ L = 512 # pixels n = 3 @@ -193,7 +198,7 @@ def test_back_project_multidim(num_ang): # apply a forward project on the image, then backwards ours_forward = imgs.project(rads) - ours_backward = ours_forward.back_project(rads) + ours_backward = ours_forward.backproject(rads) # Compare reference_back_projects = np.empty((m, n, L, L)) @@ -202,7 +207,7 @@ def test_back_project_multidim(num_ang): img = imgs[i, j] # Compute the singleton case, and compare with stack. single_sinogram = img.project(rads) - back_project = single_sinogram.back_project(rads) + back_project = single_sinogram.backproject(rads) # These should be allclose up to determinism. np.testing.assert_allclose(ours_backward[i, j : j + 1], back_project[0]) @@ -212,11 +217,12 @@ def test_back_project_multidim(num_ang): single_sinogram.asnumpy()[0].T, theta=angles[::-1], filter_name=None ) - # apply a mask, then find the NRMSE on the collection of images - # similar tolerance level to single project test + # apply a mask, then find the NRMSE on the collection of images + # similar tolerance level to single project test nrmse = np.sqrt( np.mean((ours_backward.asnumpy() * mask - reference_back_projects) ** 2) ) / (np.max(reference_back_projects) - np.min(reference_back_projects)) - assert ( - nrmse < SK_TOL_BACKPROJECT - ), f"NRMSE is too high for image ({i},{j}): {nrmse}, expected less than {SK_TOL_BACKPROJECT}" + + np.testing.assert_array_less( + nrmse, SK_TOL_BACKPROJECT, "Error with the reconstructed images." + ) From 5046aa0c688a68a7fd5fb44338f4662b5872e8ec Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Thu, 8 Aug 2024 15:42:46 -0400 Subject: [PATCH 19/25] fixup backproject tranform and tests --- src/aspire/line/line.py | 4 +++- tests/test_sinogram.py | 30 +++++++++++++++++++----------- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/src/aspire/line/line.py b/src/aspire/line/line.py index 8efee5b679..44202e04d8 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/line/line.py @@ -117,6 +117,8 @@ def backproject(self, angles): L = self.n_radial_points sinogram = fft.ifftshift(sinogram, axes=-1) sinogram_ft = fft.rfft(sinogram, axis=-1) + sinogram_ft *= xp.pi # Fix scale to match + sinogram_ft[..., 0] /= 2 # Fix DC angles = xp.asarray(angles) # grid generation with real points @@ -134,5 +136,5 @@ def backproject(self, angles): ).reshape(self.n, L, L) # normalization which gives us roughly the same error regardless of angles - imgs = imgs / (n_real_points * len(angles)) + imgs = imgs / (self.n_radial_points * len(angles)) return aspire.image.Image(xp.asnumpy(imgs)).stack_reshape(original_stack_shape) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index ed729d0b6f..6da1ee8211 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -10,7 +10,7 @@ # The same tolerance will be used in all scikit forward and backward comparisons SK_TOL_FORWARDPROJECT = 0.005 -SK_TOL_BACKPROJECT = 0.2 +SK_TOL_BACKPROJECT = 0.0025 IMG_SIZES = [ 511, @@ -162,19 +162,18 @@ def test_backproject_single(masked_image, num_ang): # generate circular mask w/ radius 1 to reconstructed image # aim to remove discrepencies for the edges of the image g = grid_2d(back_project.resolution, normalized=True, shifted=True) - mask = g["r"] < 1 + mask = g["r"] < 0.99 our_back_project = back_project.asnumpy()[0] * mask # generating sci-kit image backproject method w/ no filter - sk_image_iradon = iradon(sinogram_np[0].T, theta=angles[::-1], filter_name=None) + sk_image_iradon = iradon(sinogram_np[0].T, theta=-angles, filter_name=None) * mask # we apply a normalized root mean square error on the images to find relative error to range of ref. image - # Note: tolerance is typically < 0.2 regardless of angles, pixels, etc. nrmse = np.sqrt(np.mean((our_back_project - sk_image_iradon) ** 2)) / ( np.max(sk_image_iradon) - np.min(sk_image_iradon) ) - assert ( - nrmse < SK_TOL_BACKPROJECT + np.testing.assert_array_less( + nrmse, SK_TOL_BACKPROJECT ), f"NRMSE is too high: {nrmse}, expected less than {SK_TOL_BACKPROJECT}" @@ -189,7 +188,7 @@ def test_backproject_multidim(num_ang): m = 2 g = grid_2d(L, normalized=True, shifted=True) - mask = g["r"] < 1 + mask = g["r"] < 0.99 # Generate images imgs = Image(np.random.random((m, n, L, L))) * mask @@ -213,15 +212,24 @@ def test_backproject_multidim(num_ang): np.testing.assert_allclose(ours_backward[i, j : j + 1], back_project[0]) # Next individually compute sk's iradon transform for each image. - reference_back_projects[i, j] = iradon( - single_sinogram.asnumpy()[0].T, theta=angles[::-1], filter_name=None + reference_back_projects[i, j] = ( + iradon( + single_sinogram.asnumpy()[0].T, theta=-1 * angles, filter_name=None + ) + * mask ) # apply a mask, then find the NRMSE on the collection of images # similar tolerance level to single project test nrmse = np.sqrt( - np.mean((ours_backward.asnumpy() * mask - reference_back_projects) ** 2) - ) / (np.max(reference_back_projects) - np.min(reference_back_projects)) + np.mean( + (ours_backward.asnumpy() * mask - reference_back_projects), axis=(-2, -1) + ) + ** 2 + ) / ( + np.max(reference_back_projects, axis=(-2, -1)) + - np.min(reference_back_projects, axis=(-2, -1)) + ) np.testing.assert_array_less( nrmse, SK_TOL_BACKPROJECT, "Error with the reconstructed images." From c2a3502bfa3445b0a8439f2bc48ad7716ac9b6b9 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 9 Aug 2024 11:26:45 -0400 Subject: [PATCH 20/25] Changed name Line to Sinogram in dir/code/imports --- src/aspire/image/image.py | 4 ++-- src/aspire/line/__init__.py | 1 - src/aspire/sinogram/__init__.py | 1 + src/aspire/{line/line.py => sinogram/sinogram.py} | 10 +++++----- tests/test_sinogram.py | 6 +++--- 5 files changed, 11 insertions(+), 11 deletions(-) delete mode 100644 src/aspire/line/__init__.py create mode 100644 src/aspire/sinogram/__init__.py rename src/aspire/{line/line.py => sinogram/sinogram.py} (92%) diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py index 09a7dc56cd..57b32041cd 100644 --- a/src/aspire/image/image.py +++ b/src/aspire/image/image.py @@ -8,7 +8,7 @@ from PIL import Image as PILImage from scipy.linalg import lstsq -import aspire.line +import aspire.sinogram import aspire.volume from aspire.nufft import anufft, nufft from aspire.numeric import fft, xp @@ -240,7 +240,7 @@ def project(self, angles): # Radon transform, output: (stack size, angles, points) image_rt = fft.fftshift(fft.irfft(image_ft, n=n_points, axis=-1), axes=-1) image_rt = image_rt.reshape(*original_stack, n_angles, n_points) - return aspire.line.Line(xp.asnumpy(image_rt)) + return aspire.sinogram.Sinogram(xp.asnumpy(image_rt)) @property def res(self): diff --git a/src/aspire/line/__init__.py b/src/aspire/line/__init__.py deleted file mode 100644 index d3856d67ad..0000000000 --- a/src/aspire/line/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from .line import Line diff --git a/src/aspire/sinogram/__init__.py b/src/aspire/sinogram/__init__.py new file mode 100644 index 0000000000..98e489eedf --- /dev/null +++ b/src/aspire/sinogram/__init__.py @@ -0,0 +1 @@ +from .sinogram import Sinogram diff --git a/src/aspire/line/line.py b/src/aspire/sinogram/sinogram.py similarity index 92% rename from src/aspire/line/line.py rename to src/aspire/sinogram/sinogram.py index 44202e04d8..2eea4f0c1a 100644 --- a/src/aspire/line/line.py +++ b/src/aspire/sinogram/sinogram.py @@ -9,10 +9,10 @@ logger = logging.getLogger(__name__) -class Line: +class Sinogram: def __init__(self, data, dtype=np.float64): """ - Initialize a Line Object. This is a stack of one or more line projections or sinograms. + Initialize a Sinogram Object. This is a stack of one or more line projections or sinograms. The stack can be multidimensional with 'n' equal to the product of the stack dimensions. Singletons will be expanded into a stack @@ -47,7 +47,7 @@ def __init__(self, data, dtype=np.float64): def _check_key_dims(self, key): if isinstance(key, tuple) and (len(key) > self._data.ndim): raise ValueError( - f"Line stack_dim is {self.stack_n_dim}, slice length must be =< {self.n_dim}" + f"Sinogram stack_dim is {self.stack_n_dim}, slice length must be =< {self.n_dim}" ) def __getitem__(self, key): @@ -64,7 +64,7 @@ def stack_reshape(self, *args): :*args: Integer(s) or tuple describing the intended shape. - :returns: Line instance + :returns: Sinogram instance """ # If we're passed a tuple, use that @@ -98,7 +98,7 @@ def copy(self): return self.__class__(self._data.copy()) def __str__(self): - return f"Line(n_images = {self.n}, n_angles = {self.n_points}, n_radial_points = {self.n_radial_points})" + return f"Sinogram(n_images = {self.n}, n_angles = {self.n_points}, n_radial_points = {self.n_radial_points})" def backproject(self, angles): """ diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 6da1ee8211..a553034cdf 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -6,7 +6,7 @@ from aspire.image import Image from aspire.utils import grid_2d -# Relative tolerance comparing line projections to scikit +# Relative tolerance comparing sinogram projections to scikit # The same tolerance will be used in all scikit forward and backward comparisons SK_TOL_FORWARDPROJECT = 0.005 @@ -146,7 +146,7 @@ def test_project_multidim(num_ang): def test_backproject_single(masked_image, num_ang): """ - Test Line.backproject on a single stack of line projections (sinograms). + Test Sinogram.backproject on a single stack of line projections (sinograms). This test compares the reconstructed image from the `backproject` method to the skimage method `iradon.` @@ -179,7 +179,7 @@ def test_backproject_single(masked_image, num_ang): def test_backproject_multidim(num_ang): """ - Test Line.backproject on a stack of line projections. + Test Sinogram.backproject on a stack of line projections. Extension of the `backproject_single` test but checks for multi-dimensional stacks. """ From 9cd418032e8af5a05242f6386ccdb5b6f1696234 Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Fri, 9 Aug 2024 11:31:12 -0400 Subject: [PATCH 21/25] cleaned up tox syntax remarks --- tests/test_sinogram.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index a553034cdf..16241b5a14 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -97,7 +97,7 @@ def test_project_single(masked_image, num_ang): ) / np.linalg.norm(reference_sinogram, axis=-1) np.testing.assert_array_less( - nrms, SK_TOL_FORWARDPROJECT, "Error in image projections." + nrms, SK_TOL_FORWARDPROJECT, err_msg="Error in image projections." ) @@ -140,7 +140,7 @@ def test_project_multidim(num_ang): reference_sinograms, axis=-1 ) np.testing.assert_array_less( - _nrms, SK_TOL_FORWARDPROJECT, "Error in image projections." + _nrms, SK_TOL_FORWARDPROJECT, err_msg="Error in image projections." ) @@ -173,8 +173,10 @@ def test_backproject_single(masked_image, num_ang): np.max(sk_image_iradon) - np.min(sk_image_iradon) ) np.testing.assert_array_less( - nrmse, SK_TOL_BACKPROJECT - ), f"NRMSE is too high: {nrmse}, expected less than {SK_TOL_BACKPROJECT}" + nrmse, + SK_TOL_BACKPROJECT, + err_msg=f"NRMSE is too high: {nrmse}, expected less than {SK_TOL_BACKPROJECT}", + ) def test_backproject_multidim(num_ang): @@ -232,5 +234,5 @@ def test_backproject_multidim(num_ang): ) np.testing.assert_array_less( - nrmse, SK_TOL_BACKPROJECT, "Error with the reconstructed images." + nrmse, SK_TOL_BACKPROJECT, err_msg="Error with the reconstructed images." ) From a6720dee47bf38a66c841943bac1db6ccaf4f3d8 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Thu, 15 Aug 2024 13:42:55 -0400 Subject: [PATCH 22/25] fixed docs --- src/aspire/sinogram/sinogram.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/aspire/sinogram/sinogram.py b/src/aspire/sinogram/sinogram.py index 2eea4f0c1a..404fcdd790 100644 --- a/src/aspire/sinogram/sinogram.py +++ b/src/aspire/sinogram/sinogram.py @@ -105,9 +105,12 @@ def backproject(self, angles): Back Projection Method for a single stack of lines. :param angles: np.ndarray - 1D array of angles in radians. Each entry in the array corresponds to a different number of angles which are used to reconstruct the image. - :return: Aspire Image - An Image object containing the original stack size with a newly reconstructed numpy array of the images. Expected return shape should be (n_images, n_angles, n_radial_points) + 1D array of angles in radians. Each entry in the array + corresponds to a different number of angles which are used to + reconstruct the image. + :return: An Image object containing the original stack size + with a newly reconstructed numpy array of the images. + Expected return shape should be (..., n_radial_points, n_radial_points) """ if len(angles) != self.n_angles: raise ValueError("Number of angles must match the number of projections.") From ff02e768608576c5fd6a8158e2cb7525d3868f41 Mon Sep 17 00:00:00 2001 From: Marc Karimi Date: Mon, 19 Aug 2024 13:30:43 -0400 Subject: [PATCH 23/25] fixed minor errors and created tests for str, repr methods --- src/aspire/sinogram/sinogram.py | 26 +++++++++++++++++++------- tests/test_sinogram.py | 32 +++++++++++++++++++++++++++++++- 2 files changed, 50 insertions(+), 8 deletions(-) diff --git a/src/aspire/sinogram/sinogram.py b/src/aspire/sinogram/sinogram.py index 404fcdd790..e4d1f23784 100644 --- a/src/aspire/sinogram/sinogram.py +++ b/src/aspire/sinogram/sinogram.py @@ -10,11 +10,11 @@ class Sinogram: - def __init__(self, data, dtype=np.float64): + def __init__(self, data, dtype=None): """ Initialize a Sinogram Object. This is a stack of one or more line projections or sinograms. - The stack can be multidimensional with 'n' equal to the product + The stack can be multidimensional with 'self.n' equal to the product of the stack dimensions. Singletons will be expanded into a stack with one entry. @@ -22,14 +22,21 @@ def __init__(self, data, dtype=np.float64): `(..., angles, radial points)`. :param dtype: Optionally cast `data` to this dtype. Defaults to `data.dtype`. + + :return: Sinogram instance holding `data`. """ - self.dtype = np.dtype(dtype) + if dtype is None: + self.dtype = data.dtype + else: + self.dtype = np.dtype(dtype) + if data.ndim == 2: data = data[np.newaxis, :, :] if data.ndim < 3: raise ValueError( f"Invalid data shape: {data.shape}. Expected shape: (..., angles, radial_points), where '...' is the stack number." ) + self._data = data.astype(self.dtype, copy=False) self.ndim = self._data.ndim self.shape = self._data.shape @@ -64,7 +71,7 @@ def stack_reshape(self, *args): :*args: Integer(s) or tuple describing the intended shape. - :returns: Sinogram instance + :return: Sinogram instance """ # If we're passed a tuple, use that @@ -98,7 +105,13 @@ def copy(self): return self.__class__(self._data.copy()) def __str__(self): - return f"Sinogram(n_images = {self.n}, n_angles = {self.n_points}, n_radial_points = {self.n_radial_points})" + return f"Sinogram(n_images = {self.n}, n_angles = {self.n_angles}, n_radial_points = {self.n_radial_points})" + + def __repr__(self): + msg = f"Sinogram: {self.n} images of dtype {self.dtype}, " + msg += f"arranged as a stack with shape {self.stack_shape}. " + msg += f"Each image has {self.n_angles} angles and {self.n_radial_points} radial points." + return msg def backproject(self, angles): """ @@ -106,7 +119,7 @@ def backproject(self, angles): :param angles: np.ndarray 1D array of angles in radians. Each entry in the array - corresponds to a different number of angles which are used to + corresponds to different angles which are used to reconstruct the image. :return: An Image object containing the original stack size with a newly reconstructed numpy array of the images. @@ -138,6 +151,5 @@ def backproject(self, angles): real=True, ).reshape(self.n, L, L) - # normalization which gives us roughly the same error regardless of angles imgs = imgs / (self.n_radial_points * len(angles)) return aspire.image.Image(xp.asnumpy(imgs)).stack_reshape(original_stack_shape) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 16241b5a14..66b6ac830a 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -62,7 +62,7 @@ def masked_image(dtype, img_size): Creates a masked image fixture using camera data from Scikit-Image. """ g = grid_2d(img_size, normalized=True, shifted=True) - mask = g["r"] < 1 + mask = g["r"] < 0.99 image = data.camera().astype(dtype) image = image[:img_size, :img_size] @@ -236,3 +236,33 @@ def test_backproject_multidim(num_ang): np.testing.assert_array_less( nrmse, SK_TOL_BACKPROJECT, err_msg="Error with the reconstructed images." ) + + +# testing the str method +def test_sinogram_str_method(masked_image, num_ang): + angles = np.linspace(0, 360, num_ang, endpoint=False) + rads = angles / 180 * np.pi + sinogram = masked_image.project(rads) + n_images = sinogram.n + n_angles = sinogram.n_angles + n_radial_points = sinogram.n_radial_points + expected_str = f"Sinogram(n_images = {n_images}, n_angles = {n_angles}, n_radial_points = {n_radial_points})" + assert str(sinogram) == expected_str + + +# testing the repr method +def test_sinogram_repr_method(masked_image, num_ang): + angles = np.linspace(0, 360, num_ang, endpoint=False) + rads = angles / 180 * np.pi + sinogram = masked_image.project(rads) + n_images = sinogram.n + dtype = sinogram.dtype + stack_shape = sinogram.stack_shape + n_angles = sinogram.n_angles + n_radial_points = sinogram.n_radial_points + expected_repr = ( + f"Sinogram: {n_images} images of dtype {dtype}, " + f"arranged as a stack with shape {stack_shape}. " + f"Each image has {n_angles} angles and {n_radial_points} radial points." + ) + assert repr(sinogram) == expected_repr From 6028b7ec7a4ff024e6108f9ca5fd9addaf34338f Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 26 Aug 2024 09:11:53 -0400 Subject: [PATCH 24/25] Back Project ~> Backproject --- src/aspire/sinogram/sinogram.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aspire/sinogram/sinogram.py b/src/aspire/sinogram/sinogram.py index e4d1f23784..34451a396d 100644 --- a/src/aspire/sinogram/sinogram.py +++ b/src/aspire/sinogram/sinogram.py @@ -115,7 +115,7 @@ def __repr__(self): def backproject(self, angles): """ - Back Projection Method for a single stack of lines. + Backprojection method for a single stack of lines. :param angles: np.ndarray 1D array of angles in radians. Each entry in the array From 4d86b1ef9b17600c98ef07995e23290aef2c5c4c Mon Sep 17 00:00:00 2001 From: Garrett Wright Date: Mon, 26 Aug 2024 09:12:29 -0400 Subject: [PATCH 25/25] replace `_data` in test_sinagram --- tests/test_sinogram.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py index 66b6ac830a..bbf448db97 100644 --- a/tests/test_sinogram.py +++ b/tests/test_sinogram.py @@ -88,12 +88,13 @@ def test_project_single(masked_image, num_ang): # the original author of this method. # # Note, transpose sk output to match (angles, points) - reference_sinogram = radon(masked_image._data[0], theta=angles[::-1]).T + # Note, `radon` does not admit read only views, so the slice is copied. + reference_sinogram = radon(masked_image.asnumpy()[0].copy(), theta=angles[::-1]).T assert reference_sinogram.shape == (len(angles), ny), "Incorrect Shape" # compare project method on ski-image reference nrms = np.sqrt( - np.mean((s[0]._data - reference_sinogram) ** 2, axis=-1) + np.mean((s[0].asnumpy() - reference_sinogram) ** 2, axis=-1) ) / np.linalg.norm(reference_sinogram, axis=-1) np.testing.assert_array_less( @@ -134,7 +135,10 @@ def test_project_multidim(num_ang): np.testing.assert_allclose(s[i, j : j + 1], single_sinogram) # Next individually compute sk's radon transform for each image. - reference_sinograms[i, j] = radon(img._data[0], theta=angles[::-1]).T + # Note, `radon` does not admit read only views, so the slice is copied. + reference_sinograms[i, j] = radon( + img.asnumpy()[0].copy(), theta=angles[::-1] + ).T _nrms = np.sqrt(np.mean((s - reference_sinograms) ** 2, axis=-1)) / np.linalg.norm( reference_sinograms, axis=-1