From af3b9f5cc0ab2456e7a9cc1a40113b20209cdd8f Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Thu, 13 Jun 2024 13:41:07 -0400
Subject: [PATCH 01/24] added 2D projection stub

---
 src/aspire/image/image.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index fd160d1644..392b5fd15d 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -186,6 +186,9 @@ def __init__(self, data, dtype=None):
         self.__array_interface__ = self._data.__array_interface__
         self.__array__ = self._data
 
+    def project(self, angles):
+        """docstring"""
+
     @property
     def res(self):
         warn(

From 6db02df309e82afd47c55d559341d61b36a09ecb Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Mon, 17 Jun 2024 11:57:06 -0400
Subject: [PATCH 02/24] initial test file add

---
 tests/test_sinogram.py | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 tests/test_sinogram.py

diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
new file mode 100644
index 0000000000..5871ed8eef
--- /dev/null
+++ b/tests/test_sinogram.py
@@ -0,0 +1 @@
+import pytest

From 34cd127515e7838eaff9647a619f714272dd3d46 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Tue, 18 Jun 2024 12:33:31 -0400
Subject: [PATCH 03/24] Stashing initial project with test placeholder

---
 src/aspire/image/image.py | 37 ++++++++++++++++++++++++++++++++++++-
 tests/test_sinogram.py    | 22 ++++++++++++++++++++++
 2 files changed, 58 insertions(+), 1 deletion(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 392b5fd15d..dfa526082b 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -9,6 +9,7 @@
 from scipy.linalg import lstsq
 
 import aspire.volume
+import finufft
 from aspire.nufft import anufft
 from aspire.numeric import fft, xp
 from aspire.utils import FourierRingCorrelation, anorm, crop_pad_2d, grid_2d
@@ -187,7 +188,41 @@ def __init__(self, data, dtype=None):
         self.__array__ = self._data
 
     def project(self, angles):
-        """docstring"""
+        """docstring
+        angles: radians
+        """
+        n_points = self.resolution  # number of points to sample on radial line in polar grid
+
+        nufft_type=2
+        eps=1e-8
+        
+        n_trans = self.n_images
+        assert n_trans == 1
+
+        # 2-D grid
+        
+        y_idx = np.arange(-n_points / 2, n_points / 2) / n_points * 2
+
+        x_theta = y_idx[:, np.newaxis] * np.sin(angles)[np.newaxis, :]
+        x_theta = np.pi * x_theta.flatten()
+
+        y_theta = y_idx[:, np.newaxis] * np.cos(angles)[np.newaxis, :]
+        y_theta = np.pi * y_theta.flatten()
+
+        # NUFFT
+        plan = finufft.Plan(nufft_type, (self.resolution, self.resolution), n_trans, eps)
+        plan.setpts(x_theta, y_theta)
+
+        freqs = np.abs(np.pi * y_idx)
+        n_lines = len(angles)
+
+        # compute the polar nufft
+        image_ft = plan.execute(self._data.astype(np.complex128)).reshape(n_points, n_lines)
+
+        # compute the Radon transform (sinogram)
+
+        image_rt = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(image_ft, axes=0), axis=0), axes=0).real
+        return image_rt
 
     @property
     def res(self):
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 5871ed8eef..8179e83788 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -1 +1,23 @@
 import pytest
+from skimage import data
+from skimage.transform import radon
+from aspire.image import Image
+import numpy as np
+
+#Image.project and compare results to skimage.radon
+
+def test_image_project():
+    image = Image(data.camera().astype(np.float64))
+    ny = image.resolution
+    angles = np.linspace(0, 360, ny, endpoint=False)
+    rads = angles / 180 * np.pi
+    s = image.project(rads)
+
+    # add reference skimage radon here
+
+    
+    #compare s with reference
+    
+    
+
+

From ec00ddebfd837d3e69d80e231bbf338b104d6cd7 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Tue, 18 Jun 2024 12:35:17 -0400
Subject: [PATCH 04/24] Style Updates

---
 src/aspire/image/image.py | 26 +++++++++++++++++---------
 tests/test_sinogram.py    | 13 +++++--------
 2 files changed, 22 insertions(+), 17 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index dfa526082b..02185ff0a2 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -2,6 +2,7 @@
 import os
 from warnings import catch_warnings, filterwarnings, simplefilter, warn
 
+import finufft
 import matplotlib.pyplot as plt
 import mrcfile
 import numpy as np
@@ -9,7 +10,6 @@
 from scipy.linalg import lstsq
 
 import aspire.volume
-import finufft
 from aspire.nufft import anufft
 from aspire.numeric import fft, xp
 from aspire.utils import FourierRingCorrelation, anorm, crop_pad_2d, grid_2d
@@ -191,16 +191,18 @@ def project(self, angles):
         """docstring
         angles: radians
         """
-        n_points = self.resolution  # number of points to sample on radial line in polar grid
+        n_points = (
+            self.resolution
+        )  # number of points to sample on radial line in polar grid
+
+        nufft_type = 2
+        eps = 1e-8
 
-        nufft_type=2
-        eps=1e-8
-        
         n_trans = self.n_images
         assert n_trans == 1
 
         # 2-D grid
-        
+
         y_idx = np.arange(-n_points / 2, n_points / 2) / n_points * 2
 
         x_theta = y_idx[:, np.newaxis] * np.sin(angles)[np.newaxis, :]
@@ -210,18 +212,24 @@ def project(self, angles):
         y_theta = np.pi * y_theta.flatten()
 
         # NUFFT
-        plan = finufft.Plan(nufft_type, (self.resolution, self.resolution), n_trans, eps)
+        plan = finufft.Plan(
+            nufft_type, (self.resolution, self.resolution), n_trans, eps
+        )
         plan.setpts(x_theta, y_theta)
 
         freqs = np.abs(np.pi * y_idx)
         n_lines = len(angles)
 
         # compute the polar nufft
-        image_ft = plan.execute(self._data.astype(np.complex128)).reshape(n_points, n_lines)
+        image_ft = plan.execute(self._data.astype(np.complex128)).reshape(
+            n_points, n_lines
+        )
 
         # compute the Radon transform (sinogram)
 
-        image_rt = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(image_ft, axes=0), axis=0), axes=0).real
+        image_rt = np.fft.fftshift(
+            np.fft.ifft(np.fft.ifftshift(image_ft, axes=0), axis=0), axes=0
+        ).real
         return image_rt
 
     @property
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 8179e83788..d3b7d6fb3d 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -1,10 +1,12 @@
+import numpy as np
 import pytest
 from skimage import data
 from skimage.transform import radon
+
 from aspire.image import Image
-import numpy as np
 
-#Image.project and compare results to skimage.radon
+# Image.project and compare results to skimage.radon
+
 
 def test_image_project():
     image = Image(data.camera().astype(np.float64))
@@ -15,9 +17,4 @@ def test_image_project():
 
     # add reference skimage radon here
 
-    
-    #compare s with reference
-    
-    
-
-
+    # compare s with reference

From 53791dbdb2569c5f3efc798fcf2cbae15cc3199a Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 21 Jun 2024 12:00:34 -0400
Subject: [PATCH 05/24] Pytest fixtures

---
 tests/test_sinogram.py | 62 ++++++++++++++++++++++++++++++++++++++----
 1 file changed, 56 insertions(+), 6 deletions(-)

diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index d3b7d6fb3d..cf095f649a 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -1,20 +1,70 @@
+import itertools
+
 import numpy as np
 import pytest
 from skimage import data
-from skimage.transform import radon
+from skimage.transform import radon, resize
 
 from aspire.image import Image
+from aspire.utils import grid_2d
+
+# parameter img_sizes: 511, 512
+IMG_SIZES = [
+    511,
+    512,
+]
+
+# parameter dtype: float32, float64
+DTYPES = [
+    np.float32,
+    np.float64,
+]
+
+
+@pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
+def dtype(request):
+    """
+    Dtypes for image.
+    """
+    return request.param
 
-# Image.project and compare results to skimage.radon
 
+@pytest.fixture(params=IMG_SIZES, ids=lambda x: f"px={x}", scope="module")
+def img_size(request):
+    """
+    Image size.
+    """
+    return request.param
 
-def test_image_project():
-    image = Image(data.camera().astype(np.float64))
-    ny = image.resolution
+
+@pytest.fixture
+def masked_image(dtype, img_size):
+    """
+    Construct a masked image fixture that takes paramters
+    """
+    g = grid_2d(img_size, normalized=True, shifted=True)
+    mask = g["r"] < 1
+
+    # add more logic to check the sizes and readjust accordingly
+    image = data.camera().astype(dtype)
+    image = image[:img_size, :img_size]
+    return Image(image * mask)
+
+
+# Image.project and compare results to skimage.radon
+def test_image_project(masked_image):
+    ny = masked_image.resolution
     angles = np.linspace(0, 360, ny, endpoint=False)
     rads = angles / 180 * np.pi
-    s = image.project(rads)
+    s = masked_image.project(rads)
 
     # add reference skimage radon here
+    n = masked_image._data[0]
+    print(s.shape)
+    print(n.shape)
+    reference_sinogram = radon(n, theta=angles)
 
     # compare s with reference
+    np.testing.assert_allclose(s, reference_sinogram, rtol=11, atol=1e-8)
+
+    # create fixture called masked_image(img_size) -> return: masked image of size (grid generation goes in fixture)

From 705739219415d7393cd061256cdfe07d3b842516 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 21 Jun 2024 12:21:14 -0400
Subject: [PATCH 06/24] Cleanup

---
 src/aspire/image/image.py | 6 ++----
 tests/test_sinogram.py    | 8 ++------
 2 files changed, 4 insertions(+), 10 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 02185ff0a2..40546aff1a 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -191,9 +191,8 @@ def project(self, angles):
         """docstring
         angles: radians
         """
-        n_points = (
-            self.resolution
-        )  # number of points to sample on radial line in polar grid
+        # number of points to sample on radial line in polar grid
+        n_points = self.resolution
 
         nufft_type = 2
         eps = 1e-8
@@ -217,7 +216,6 @@ def project(self, angles):
         )
         plan.setpts(x_theta, y_theta)
 
-        freqs = np.abs(np.pi * y_idx)
         n_lines = len(angles)
 
         # compute the polar nufft
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index cf095f649a..de65ce252c 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -1,9 +1,7 @@
-import itertools
-
 import numpy as np
 import pytest
 from skimage import data
-from skimage.transform import radon, resize
+from skimage.transform import radon
 
 from aspire.image import Image
 from aspire.utils import grid_2d
@@ -60,9 +58,7 @@ def test_image_project(masked_image):
 
     # add reference skimage radon here
     n = masked_image._data[0]
-    print(s.shape)
-    print(n.shape)
-    reference_sinogram = radon(n, theta=angles)
+    reference_sinogram = radon(n, theta=angles[::-1])
 
     # compare s with reference
     np.testing.assert_allclose(s, reference_sinogram, rtol=11, atol=1e-8)

From 8822a6c49308e0d4ad97b55a919f2456a03f3bc2 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 21 Jun 2024 13:01:21 -0400
Subject: [PATCH 07/24] changed nufft call

---
 src/aspire/image/image.py | 24 +++++-------------------
 1 file changed, 5 insertions(+), 19 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 40546aff1a..2d771766b0 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -2,7 +2,6 @@
 import os
 from warnings import catch_warnings, filterwarnings, simplefilter, warn
 
-import finufft
 import matplotlib.pyplot as plt
 import mrcfile
 import numpy as np
@@ -10,7 +9,7 @@
 from scipy.linalg import lstsq
 
 import aspire.volume
-from aspire.nufft import anufft
+from aspire.nufft import anufft, nufft
 from aspire.numeric import fft, xp
 from aspire.utils import FourierRingCorrelation, anorm, crop_pad_2d, grid_2d
 from aspire.volume import SymmetryGroup
@@ -194,37 +193,24 @@ def project(self, angles):
         # number of points to sample on radial line in polar grid
         n_points = self.resolution
 
-        nufft_type = 2
-        eps = 1e-8
-
         n_trans = self.n_images
         assert n_trans == 1
 
         # 2-D grid
-
         y_idx = np.arange(-n_points / 2, n_points / 2) / n_points * 2
+        pts = np.empty((2, n_points * len(angles)), dtype=self.dtype)
 
         x_theta = y_idx[:, np.newaxis] * np.sin(angles)[np.newaxis, :]
-        x_theta = np.pi * x_theta.flatten()
+        pts[0] = np.pi * x_theta.flatten()
 
         y_theta = y_idx[:, np.newaxis] * np.cos(angles)[np.newaxis, :]
-        y_theta = np.pi * y_theta.flatten()
+        pts[1] = np.pi * y_theta.flatten()
 
         # NUFFT
-        plan = finufft.Plan(
-            nufft_type, (self.resolution, self.resolution), n_trans, eps
-        )
-        plan.setpts(x_theta, y_theta)
-
-        n_lines = len(angles)
-
         # compute the polar nufft
-        image_ft = plan.execute(self._data.astype(np.complex128)).reshape(
-            n_points, n_lines
-        )
+        image_ft = nufft(self._data, pts).reshape(n_points, n_points)
 
         # compute the Radon transform (sinogram)
-
         image_rt = np.fft.fftshift(
             np.fft.ifft(np.fft.ifftshift(image_ft, axes=0), axis=0), axes=0
         ).real

From cf157cbec2f186de264ecd07549ce4b089133ef1 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 21 Jun 2024 13:18:49 -0400
Subject: [PATCH 08/24] added stub for image stack line project

marc to continue
---
 tests/test_sinogram.py | 30 ++++++++++++++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index de65ce252c..020e8b58bf 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.source import Simulation
 from aspire.utils import grid_2d
 
 # parameter img_sizes: 511, 512
@@ -64,3 +65,32 @@ def test_image_project(masked_image):
     np.testing.assert_allclose(s, reference_sinogram, rtol=11, atol=1e-8)
 
     # create fixture called masked_image(img_size) -> return: masked image of size (grid generation goes in fixture)
+
+
+def test_multidim():
+    """
+    Test Image.project on stacks of images.
+    """
+
+    L = 32  # pixels
+    n = 3
+
+    # Generate a mask
+    g = grid_2d(L, normalized=True, shifted=True)
+    mask = g["r"] < 1
+
+    # Generate a simulation
+    src = Simulation(n=n, L=L, C=1, dtype=np.float64)
+    imgs = src.images[:]
+
+    # Generate line project angles
+    ang_degrees = np.linspace(0, 180, L)
+    ang_rads = ang_degrees * np.pi / 180.0
+
+    # Call the line projection method
+    s = imgs.project(ang_rads)
+
+    # # Compare with sk
+    # res = np.empty((n,L,L))
+    # for i,img in enumerate(imgs.asnumpy()):
+    #     #res[i] = radon(img ...)

From 6d0fa474c25ffc4f9c1267a072307cfb4cdd6716 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Mon, 24 Jun 2024 12:03:41 -0400
Subject: [PATCH 09/24] Dimensional Test Fix

---
 src/aspire/image/image.py | 14 ++++++++------
 tests/test_sinogram.py    | 16 +++++++++-------
 2 files changed, 17 insertions(+), 13 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 2d771766b0..ee0b46e3f0 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -194,7 +194,6 @@ def project(self, angles):
         n_points = self.resolution
 
         n_trans = self.n_images
-        assert n_trans == 1
 
         # 2-D grid
         y_idx = np.arange(-n_points / 2, n_points / 2) / n_points * 2
@@ -207,13 +206,16 @@ def project(self, angles):
         pts[1] = np.pi * y_theta.flatten()
 
         # NUFFT
-        # compute the polar nufft
-        image_ft = nufft(self._data, pts).reshape(n_points, n_points)
+        # compute the polar nufft, create a
+        image_ft = nufft(self._data, pts).reshape(self.n_images, n_points, n_points)
 
         # compute the Radon transform (sinogram)
-        image_rt = np.fft.fftshift(
-            np.fft.ifft(np.fft.ifftshift(image_ft, axes=0), axis=0), axes=0
-        ).real
+        image_rt = np.empty((self.n_images, n_points, n_points))
+        for i in range(n_trans):
+            image_rt[i] = np.fft.fftshift(
+                np.fft.ifft(np.fft.ifftshift(image_ft[i], axes=0), axis=0), axes=0
+            ).real
+        # previous code: image_rt = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(image_ft, axes=0), axis=0), axes=0).real
         return image_rt
 
     @property
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 020e8b58bf..b9b2f3912f 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -62,7 +62,7 @@ def test_image_project(masked_image):
     reference_sinogram = radon(n, theta=angles[::-1])
 
     # compare s with reference
-    np.testing.assert_allclose(s, reference_sinogram, rtol=11, atol=1e-8)
+    np.testing.assert_allclose(s[0], reference_sinogram, rtol=11, atol=1e-8)
 
     # create fixture called masked_image(img_size) -> return: masked image of size (grid generation goes in fixture)
 
@@ -72,7 +72,7 @@ def test_multidim():
     Test Image.project on stacks of images.
     """
 
-    L = 32  # pixels
+    L = 64  # pixels
     n = 3
 
     # Generate a mask
@@ -81,16 +81,18 @@ def test_multidim():
 
     # Generate a simulation
     src = Simulation(n=n, L=L, C=1, dtype=np.float64)
-    imgs = src.images[:]
+    imgs = src.images[:] * mask
 
     # Generate line project angles
-    ang_degrees = np.linspace(0, 180, L)
+    ang_degrees = np.linspace(0, 180, L, endpoint=False)
     ang_rads = ang_degrees * np.pi / 180.0
 
     # Call the line projection method
     s = imgs.project(ang_rads)
 
     # # Compare with sk
-    # res = np.empty((n,L,L))
-    # for i,img in enumerate(imgs.asnumpy()):
-    #     #res[i] = radon(img ...)
+    res = np.empty((n, L, L))
+    for i, img in enumerate(imgs._data):
+        res[i] = radon(img, theta=ang_rads[::-1])
+
+    np.testing.assert_allclose(s, res, rtol=12, atol=1e-8)

From eccb18240206154f95757a1571d81f4b219b8474 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Mon, 24 Jun 2024 12:26:06 -0400
Subject: [PATCH 10/24] Multidim FFT

---
 src/aspire/image/image.py | 12 +++++-------
 1 file changed, 5 insertions(+), 7 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index ee0b46e3f0..c88a17e6a3 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -196,7 +196,7 @@ def project(self, angles):
         n_trans = self.n_images
 
         # 2-D grid
-        y_idx = np.arange(-n_points / 2, n_points / 2) / n_points * 2
+        y_idx = np.arange(-n_points / 2, n_points / 2, dtype=self.dtype) / n_points * 2
         pts = np.empty((2, n_points * len(angles)), dtype=self.dtype)
 
         x_theta = y_idx[:, np.newaxis] * np.sin(angles)[np.newaxis, :]
@@ -210,12 +210,10 @@ def project(self, angles):
         image_ft = nufft(self._data, pts).reshape(self.n_images, n_points, n_points)
 
         # compute the Radon transform (sinogram)
-        image_rt = np.empty((self.n_images, n_points, n_points))
-        for i in range(n_trans):
-            image_rt[i] = np.fft.fftshift(
-                np.fft.ifft(np.fft.ifftshift(image_ft[i], axes=0), axis=0), axes=0
-            ).real
-        # previous code: image_rt = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(image_ft, axes=0), axis=0), axes=0).real
+        image_rt = np.fft.fftshift(
+            np.fft.ifftn(np.fft.ifftshift(image_ft, axes=(0, 1)), axes=(0, 1)),
+            axes=(0, 1),
+        ).real
         return image_rt
 
     @property

From b7de220611e27f7217db2e893036aff4e0bd5a74 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Tue, 25 Jun 2024 11:43:25 -0400
Subject: [PATCH 11/24] Integrated stack reshape to project

---
 src/aspire/image/image.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index c88a17e6a3..f30c82ca3b 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -192,8 +192,7 @@ def project(self, angles):
         """
         # number of points to sample on radial line in polar grid
         n_points = self.resolution
-
-        n_trans = self.n_images
+        original_stack = self.stack_shape
 
         # 2-D grid
         y_idx = np.arange(-n_points / 2, n_points / 2, dtype=self.dtype) / n_points * 2
@@ -206,14 +205,17 @@ def project(self, angles):
         pts[1] = np.pi * y_theta.flatten()
 
         # NUFFT
-        # compute the polar nufft, create a
-        image_ft = nufft(self._data, pts).reshape(self.n_images, n_points, n_points)
+        # compute the polar nufft
+        image_ft = nufft(self.stack_reshape(-1)._data, pts).reshape(
+            self.n_images, n_points, n_points
+        )
 
         # compute the Radon transform (sinogram)
         image_rt = np.fft.fftshift(
             np.fft.ifftn(np.fft.ifftshift(image_ft, axes=(0, 1)), axes=(0, 1)),
             axes=(0, 1),
         ).real
+        image_rt = image_rt.reshape(*original_stack, n_points, n_points)
         return image_rt
 
     @property

From 654cbbc941f13ef2d00f9a15774c42347b6363b8 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Wed, 26 Jun 2024 15:46:26 -0400
Subject: [PATCH 12/24] Fleshed out Image Project Single and Multidim Tests

---
 tests/test_sinogram.py | 37 ++++++++++++++++++++++++-------------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index b9b2f3912f..eb42b0a871 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -44,7 +44,6 @@ def masked_image(dtype, img_size):
     g = grid_2d(img_size, normalized=True, shifted=True)
     mask = g["r"] < 1
 
-    # add more logic to check the sizes and readjust accordingly
     image = data.camera().astype(dtype)
     image = image[:img_size, :img_size]
     return Image(image * mask)
@@ -52,6 +51,9 @@ def masked_image(dtype, img_size):
 
 # Image.project and compare results to skimage.radon
 def test_image_project(masked_image):
+    """
+    TestImage.project on a single stack of images. Compares project method with skimage.
+    """
     ny = masked_image.resolution
     angles = np.linspace(0, 360, ny, endpoint=False)
     rads = angles / 180 * np.pi
@@ -62,9 +64,15 @@ def test_image_project(masked_image):
     reference_sinogram = radon(n, theta=angles[::-1])
 
     # compare s with reference
-    np.testing.assert_allclose(s[0], reference_sinogram, rtol=11, atol=1e-8)
+    nrms = np.sqrt(np.mean((s[0] - reference_sinogram) ** 2, axis=0)) / np.linalg.norm(
+        reference_sinogram, axis=0
+    )
+    tol = 0.002
 
-    # create fixture called masked_image(img_size) -> return: masked image of size (grid generation goes in fixture)
+    # odd image tolerance (stink)
+    if masked_image.resolution % 2 == 1:
+        tol = 0.02
+    np.testing.assert_array_less(nrms, tol, "Error in test image")
 
 
 def test_multidim():
@@ -72,7 +80,7 @@ def test_multidim():
     Test Image.project on stacks of images.
     """
 
-    L = 64  # pixels
+    L = 512  # pixels
     n = 3
 
     # Generate a mask
@@ -84,15 +92,18 @@ def test_multidim():
     imgs = src.images[:] * mask
 
     # Generate line project angles
-    ang_degrees = np.linspace(0, 180, L, endpoint=False)
-    ang_rads = ang_degrees * np.pi / 180.0
-
-    # Call the line projection method
-    s = imgs.project(ang_rads)
+    angles = np.linspace(0, 180, L, endpoint=False)
+    rads = angles / 180.0 * np.pi
+    s = imgs.project(rads)
 
     # # Compare with sk
-    res = np.empty((n, L, L))
+    reference_sinograms = np.empty((n, L, L))
     for i, img in enumerate(imgs._data):
-        res[i] = radon(img, theta=ang_rads[::-1])
-
-    np.testing.assert_allclose(s, res, rtol=12, atol=1e-8)
+        reference_sinograms[i] = radon(img, theta=angles[::-1])
+
+    # decrease tolerance as L goes up
+    for i in range(n):
+        nrms = np.sqrt(
+            np.mean((s[i] - reference_sinograms[i]) ** 2, axis=0)
+        ) / np.linalg.norm(reference_sinograms[i], axis=0)
+        np.testing.assert_array_less(nrms, 0.05, err_msg=f"Error in image {i}")

From dd41f0f0429244e97ae4481f55255886bd7c88e0 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Thu, 27 Jun 2024 12:35:05 -0400
Subject: [PATCH 13/24] Fixed the grid issues yay

---
 src/aspire/image/image.py | 16 +++++++---------
 tests/test_sinogram.py    |  6 +-----
 2 files changed, 8 insertions(+), 14 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index f30c82ca3b..c38e57dcde 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -195,19 +195,17 @@ def project(self, angles):
         original_stack = self.stack_shape
 
         # 2-D grid
-        y_idx = np.arange(-n_points / 2, n_points / 2, dtype=self.dtype) / n_points * 2
-        pts = np.empty((2, n_points * len(angles)), dtype=self.dtype)
+        y_idx = np.fft.fftshift(np.fft.fftfreq(n_points))
+        pts = np.empty((2, n_points, len(angles)), dtype=self.dtype)
 
-        x_theta = y_idx[:, np.newaxis] * np.sin(angles)[np.newaxis, :]
-        pts[0] = np.pi * x_theta.flatten()
-
-        y_theta = y_idx[:, np.newaxis] * np.cos(angles)[np.newaxis, :]
-        pts[1] = np.pi * y_theta.flatten()
+        pts[0] = y_idx[:, np.newaxis] * np.sin(angles)[np.newaxis, :]
+        pts[1] = y_idx[:, np.newaxis] * np.cos(angles)[np.newaxis, :]
+        pts = pts.reshape(2, n_points * len(angles)) * 2 * np.pi
 
         # NUFFT
         # compute the polar nufft
         image_ft = nufft(self.stack_reshape(-1)._data, pts).reshape(
-            self.n_images, n_points, n_points
+            self.n_images, n_points, len(angles)
         )
 
         # compute the Radon transform (sinogram)
@@ -215,7 +213,7 @@ def project(self, angles):
             np.fft.ifftn(np.fft.ifftshift(image_ft, axes=(0, 1)), axes=(0, 1)),
             axes=(0, 1),
         ).real
-        image_rt = image_rt.reshape(*original_stack, n_points, n_points)
+        image_rt = image_rt.reshape(*original_stack, n_points, len(angles))
         return image_rt
 
     @property
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index eb42b0a871..efa2c4aa3d 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -55,7 +55,7 @@ def test_image_project(masked_image):
     TestImage.project on a single stack of images. Compares project method with skimage.
     """
     ny = masked_image.resolution
-    angles = np.linspace(0, 360, ny, endpoint=False)
+    angles = np.linspace(0, 360, ny + 1, endpoint=False)
     rads = angles / 180 * np.pi
     s = masked_image.project(rads)
 
@@ -68,10 +68,6 @@ def test_image_project(masked_image):
         reference_sinogram, axis=0
     )
     tol = 0.002
-
-    # odd image tolerance (stink)
-    if masked_image.resolution % 2 == 1:
-        tol = 0.02
     np.testing.assert_array_less(nrms, tol, "Error in test image")
 
 

From 8473a218df1fcf4ca0702a27a7b0d30a3f15a64f Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 28 Jun 2024 11:52:09 -0400
Subject: [PATCH 14/24] Angle slow moving axis

---
 src/aspire/image/image.py | 14 +++++++-------
 tests/test_sinogram.py    | 10 ++++++----
 2 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index c38e57dcde..92aa780823 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -196,24 +196,24 @@ def project(self, angles):
 
         # 2-D grid
         y_idx = np.fft.fftshift(np.fft.fftfreq(n_points))
-        pts = np.empty((2, n_points, len(angles)), dtype=self.dtype)
+        pts = np.empty((2, len(angles), n_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[0] = y_idx[np.newaxis, :] * np.sin(angles)[:, np.newaxis]
+        pts[1] = y_idx[np.newaxis, :] * np.cos(angles)[:, np.newaxis]
         pts = pts.reshape(2, n_points * len(angles)) * 2 * np.pi
 
         # NUFFT
         # compute the polar nufft
         image_ft = nufft(self.stack_reshape(-1)._data, pts).reshape(
-            self.n_images, n_points, len(angles)
+            self.n_images, len(angles), n_points
         )
 
         # compute the Radon transform (sinogram)
         image_rt = np.fft.fftshift(
-            np.fft.ifftn(np.fft.ifftshift(image_ft, axes=(0, 1)), axes=(0, 1)),
-            axes=(0, 1),
+            np.fft.ifftn(np.fft.ifftshift(image_ft, axes=(0, 2)), axes=(0, 2)),
+            axes=(0, 2),
         ).real
-        image_rt = image_rt.reshape(*original_stack, n_points, len(angles))
+        image_rt = image_rt.reshape(*original_stack, len(angles), n_points)
         return image_rt
 
     @property
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index efa2c4aa3d..2abe42f12e 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -58,14 +58,16 @@ def test_image_project(masked_image):
     angles = np.linspace(0, 360, ny + 1, endpoint=False)
     rads = angles / 180 * np.pi
     s = masked_image.project(rads)
+    assert s.shape == (1, len(angles), ny)
 
     # add reference skimage radon here
     n = masked_image._data[0]
-    reference_sinogram = radon(n, theta=angles[::-1])
-
+    reference_sinogram = radon(n, theta=angles[::-1]).T  # transpose angles, points
+    assert reference_sinogram.shape == (len(angles), ny)
     # compare s with reference
-    nrms = np.sqrt(np.mean((s[0] - reference_sinogram) ** 2, axis=0)) / np.linalg.norm(
-        reference_sinogram, axis=0
+
+    nrms = np.sqrt(np.mean((s[0] - reference_sinogram) ** 2, axis=1)) / np.linalg.norm(
+        reference_sinogram, axis=1
     )
     tol = 0.002
     np.testing.assert_array_less(nrms, tol, "Error in test image")

From b736730443758db4b537dd09927b869499b23824 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 28 Jun 2024 15:50:29 -0400
Subject: [PATCH 15/24] Replaced FFT with rfft

---
 src/aspire/image/image.py | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 92aa780823..e12cf29598 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -5,6 +5,7 @@
 import matplotlib.pyplot as plt
 import mrcfile
 import numpy as np
+from numpy.fft import irfft
 from PIL import Image as PILImage
 from scipy.linalg import lstsq
 
@@ -195,24 +196,27 @@ def project(self, angles):
         original_stack = self.stack_shape
 
         # 2-D grid
-        y_idx = np.fft.fftshift(np.fft.fftfreq(n_points))
-        pts = np.empty((2, len(angles), n_points), dtype=self.dtype)
+        y_idx = np.fft.rfftfreq(n_points) * np.pi * 2
+        n_real_points = len(y_idx)
+
+        # y_idx = np.fft.fftshift(np.fft.fftfreq(n_points))
+        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 = pts.reshape(2, n_points * len(angles)) * 2 * np.pi
+        pts = pts.reshape(2, n_real_points * len(angles))
 
         # NUFFT
         # compute the polar nufft
         image_ft = nufft(self.stack_reshape(-1)._data, pts).reshape(
-            self.n_images, len(angles), n_points
+            self.n_images, len(angles), n_real_points
         )
 
         # compute the Radon transform (sinogram)
         image_rt = np.fft.fftshift(
-            np.fft.ifftn(np.fft.ifftshift(image_ft, axes=(0, 2)), axes=(0, 2)),
+            np.fft.irfftn(image_ft, s=(self.n_images, n_points), axes=(0, 2)),
             axes=(0, 2),
-        ).real
+        )
         image_rt = image_rt.reshape(*original_stack, len(angles), n_points)
         return image_rt
 

From c3ecb8f1bf68d79b5c82dbcdae05c664b413287b Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 28 Jun 2024 16:00:17 -0400
Subject: [PATCH 16/24] Cleaned up other unit tests

---
 src/aspire/image/image.py |  1 -
 tests/test_sinogram.py    | 12 +++++-------
 2 files changed, 5 insertions(+), 8 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index e12cf29598..f0ef24dcc7 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -5,7 +5,6 @@
 import matplotlib.pyplot as plt
 import mrcfile
 import numpy as np
-from numpy.fft import irfft
 from PIL import Image as PILImage
 from scipy.linalg import lstsq
 
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 2abe42f12e..fc8e712408 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -4,7 +4,6 @@
 from skimage.transform import radon
 
 from aspire.image import Image
-from aspire.source import Simulation
 from aspire.utils import grid_2d
 
 # parameter img_sizes: 511, 512
@@ -85,9 +84,8 @@ def test_multidim():
     g = grid_2d(L, normalized=True, shifted=True)
     mask = g["r"] < 1
 
-    # Generate a simulation
-    src = Simulation(n=n, L=L, C=1, dtype=np.float64)
-    imgs = src.images[:] * mask
+    # Generate images
+    imgs = Image(np.random.random((n, L, L))) * mask
 
     # Generate line project angles
     angles = np.linspace(0, 180, L, endpoint=False)
@@ -97,11 +95,11 @@ def test_multidim():
     # # Compare with sk
     reference_sinograms = np.empty((n, L, L))
     for i, img in enumerate(imgs._data):
-        reference_sinograms[i] = radon(img, theta=angles[::-1])
+        reference_sinograms[i] = radon(img, theta=angles[::-1]).T
 
     # decrease tolerance as L goes up
     for i in range(n):
         nrms = np.sqrt(
-            np.mean((s[i] - reference_sinograms[i]) ** 2, axis=0)
-        ) / np.linalg.norm(reference_sinograms[i], axis=0)
+            np.mean((s[i] - reference_sinograms[i]) ** 2, axis=1)
+        ) / np.linalg.norm(reference_sinograms[i], axis=1)
         np.testing.assert_array_less(nrms, 0.05, err_msg=f"Error in image {i}")

From 0fab57f0d67817cc3a64771309a4e91d3f08b3d0 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 5 Jul 2024 20:48:03 -0400
Subject: [PATCH 17/24] Added Doc Test and Cleaned up Code

---
 src/aspire/image/image.py | 15 ++++++++-------
 tests/test_sinogram.py    | 17 ++++++++---------
 2 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index f0ef24dcc7..d5d1b9885e 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -187,8 +187,12 @@ def __init__(self, data, dtype=None):
         self.__array__ = self._data
 
     def project(self, angles):
-        """docstring
-        angles: radians
+        """
+        Computes the Radon Transform on an Image Stack using Non-Uniform Fast Fourier Transforms. This method projects the Image stack along different angles and returns the Radon Transform.
+
+        :param angles: A 1-D Numpy Array of angles in Radians. This is used to compute the Radon Transform at different angles.
+        :return: Radon transform of the Image Stack.
+        :rtype: Ndarray (stack size, number of angles, image resolution)
         """
         # number of points to sample on radial line in polar grid
         n_points = self.resolution
@@ -198,20 +202,17 @@ def project(self, angles):
         y_idx = np.fft.rfftfreq(n_points) * np.pi * 2
         n_real_points = len(y_idx)
 
-        # y_idx = np.fft.fftshift(np.fft.fftfreq(n_points))
         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 = pts.reshape(2, n_real_points * len(angles))
 
-        # NUFFT
-        # compute the polar nufft
+        # compute the polar nufft (NUFFT)
         image_ft = nufft(self.stack_reshape(-1)._data, pts).reshape(
             self.n_images, len(angles), n_real_points
         )
 
-        # compute the Radon transform (sinogram)
+        # Radon transform
         image_rt = np.fft.fftshift(
             np.fft.irfftn(image_ft, s=(self.n_images, n_points), axes=(0, 2)),
             axes=(0, 2),
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index fc8e712408..a935f68ccc 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -38,7 +38,7 @@ def img_size(request):
 @pytest.fixture
 def masked_image(dtype, img_size):
     """
-    Construct a masked image fixture that takes paramters
+    Creates a masked image fixture using camera data from Skikit-Image.
     """
     g = grid_2d(img_size, normalized=True, shifted=True)
     mask = g["r"] < 1
@@ -51,7 +51,7 @@ def masked_image(dtype, img_size):
 # Image.project and compare results to skimage.radon
 def test_image_project(masked_image):
     """
-    TestImage.project on a single stack of images. Compares project method with skimage.
+    Test Image.project on a single stack of images. Compares project method with skimage.
     """
     ny = masked_image.resolution
     angles = np.linspace(0, 360, ny + 1, endpoint=False)
@@ -59,22 +59,22 @@ def test_image_project(masked_image):
     s = masked_image.project(rads)
     assert s.shape == (1, len(angles), ny)
 
-    # add reference skimage radon here
+    # ski-kit image radon reference
     n = masked_image._data[0]
     reference_sinogram = radon(n, theta=angles[::-1]).T  # transpose angles, points
     assert reference_sinogram.shape == (len(angles), ny)
-    # compare s with reference
 
+    # 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
     )
     tol = 0.002
-    np.testing.assert_array_less(nrms, tol, "Error in test image")
+    np.testing.assert_array_less(nrms, tol, "Error in image projections.")
 
 
 def test_multidim():
     """
-    Test Image.project on stacks of images.
+    Test Image.project on stacks of images. Extension of test_image_project but for multi-dimensional stacks.
     """
 
     L = 512  # pixels
@@ -92,14 +92,13 @@ def test_multidim():
     rads = angles / 180.0 * np.pi
     s = imgs.project(rads)
 
-    # # Compare with sk
+    # Compare with ski-image
     reference_sinograms = np.empty((n, L, L))
     for i, img in enumerate(imgs._data):
         reference_sinograms[i] = radon(img, theta=angles[::-1]).T
 
-    # decrease tolerance as L goes up
     for i in range(n):
         nrms = np.sqrt(
             np.mean((s[i] - reference_sinograms[i]) ** 2, axis=1)
         ) / np.linalg.norm(reference_sinograms[i], axis=1)
-        np.testing.assert_array_less(nrms, 0.05, err_msg=f"Error in image {i}")
+        np.testing.assert_array_less(nrms, 0.05, err_msg=f"Error in image {i}.")

From c98c93eca88f222c20a0e6e978ecd0e2398a08e8 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Wed, 10 Jul 2024 19:59:12 -0400
Subject: [PATCH 18/24] fixup sinogram tests and simpler multi test

Co-authored-by: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
---
 tests/test_sinogram.py | 37 +++++++++++++++++++++++++++----------
 1 file changed, 27 insertions(+), 10 deletions(-)

diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index a935f68ccc..323c1f2eb0 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -6,6 +6,10 @@
 from aspire.image import Image
 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.002
+
 # parameter img_sizes: 511, 512
 IMG_SIZES = [
     511,
@@ -59,17 +63,23 @@ def test_image_project(masked_image):
     s = masked_image.project(rads)
     assert s.shape == (1, len(angles), ny)
 
-    # ski-kit image radon reference
-    n = masked_image._data[0]
-    reference_sinogram = radon(n, theta=angles[::-1]).T  # transpose angles, points
-    assert reference_sinogram.shape == (len(angles), ny)
+    # sci-kit image `radon` reference
+    #
+    # Note, Image.project's angles are wrt projection line (ie
+    # grid), while sk's radon are wrt the image. To correspond the
+    # rotations are inverted.  This was the convention prefered by
+    # 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
+    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
     )
-    tol = 0.002
-    np.testing.assert_array_less(nrms, tol, "Error in image projections.")
+
+    np.testing.assert_array_less(nrms, SK_TOL, "Error in image projections.")
 
 
 def test_multidim():
@@ -92,13 +102,20 @@ def test_multidim():
     rads = angles / 180.0 * np.pi
     s = imgs.project(rads)
 
-    # Compare with ski-image
+    # Compare
     reference_sinograms = np.empty((n, L, L))
-    for i, img in enumerate(imgs._data):
-        reference_sinograms[i] = radon(img, theta=angles[::-1]).T
+    for i, img in enumerate(imgs):
+        # Compute the singleton case, and compare with the stack
+        single_sinogram = img.project(rads)
+        # These should be allclose up to determinism in the FFT and NUFFT.
+        np.testing.assert_allclose(s[i : i + 1], single_sinogram)
+
+        # Next individually compute sk's radon transform for each image.
+        reference_sinograms[i] = radon(img._data[0], theta=angles[::-1]).T
 
+    # Compare all lines in each sinogram with sk-image
     for i in range(n):
         nrms = np.sqrt(
             np.mean((s[i] - reference_sinograms[i]) ** 2, axis=1)
         ) / np.linalg.norm(reference_sinograms[i], axis=1)
-        np.testing.assert_array_less(nrms, 0.05, err_msg=f"Error in image {i}.")
+        np.testing.assert_array_less(nrms, SK_TOL, err_msg=f"Error in image {i}.")

From 7740110137e76a83afe876c86bff6731be69e298 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Wed, 10 Jul 2024 20:07:57 -0400
Subject: [PATCH 19/24] fix irfft and shift

Co-authored-by: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
---
 src/aspire/image/image.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index d5d1b9885e..dd0adefd13 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -213,10 +213,7 @@ def project(self, angles):
         )
 
         # Radon transform
-        image_rt = np.fft.fftshift(
-            np.fft.irfftn(image_ft, s=(self.n_images, n_points), axes=(0, 2)),
-            axes=(0, 2),
-        )
+        image_rt = np.fft.fftshift(np.fft.irfft(image_ft, n=n_points, axis=-1), axes=-1)
         image_rt = image_rt.reshape(*original_stack, len(angles), n_points)
         return image_rt
 

From d9d498dc4809957ff5c9908f567a6ec421988b59 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Thu, 11 Jul 2024 12:49:52 -0400
Subject: [PATCH 20/24] added angles but need to change multidim

---
 src/aspire/image/image.py | 17 +++++++++--------
 tests/test_sinogram.py    | 14 ++++++++++----
 2 files changed, 19 insertions(+), 12 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index dd0adefd13..f6e07dcdc6 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -199,22 +199,23 @@ def project(self, angles):
         original_stack = self.stack_shape
 
         # 2-D grid
-        y_idx = np.fft.rfftfreq(n_points) * np.pi * 2
-        n_real_points = len(y_idx)
+        radial_idx = np.fft.rfftfreq(n_points) * np.pi * 2
+        n_real_points = len(radial_idx)
+        n_angles = len(angles)
 
-        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 = pts.reshape(2, n_real_points * 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 = pts.reshape(2, n_real_points * n_angles)
 
         # compute the polar nufft (NUFFT)
         image_ft = nufft(self.stack_reshape(-1)._data, pts).reshape(
-            self.n_images, len(angles), n_real_points
+            self.n_images, n_angles, n_real_points
         )
 
         # Radon transform
         image_rt = np.fft.fftshift(np.fft.irfft(image_ft, n=n_points, axis=-1), axes=-1)
-        image_rt = image_rt.reshape(*original_stack, len(angles), n_points)
+        image_rt = image_rt.reshape(*original_stack, n_angles, n_points)
         return image_rt
 
     @property
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 323c1f2eb0..94cc172189 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -10,18 +10,17 @@
 # The same tolerance will be used in all scikit comparisons
 SK_TOL = 0.002
 
-# parameter img_sizes: 511, 512
 IMG_SIZES = [
     511,
     512,
 ]
 
-# parameter dtype: float32, float64
 DTYPES = [
     np.float32,
     np.float64,
 ]
 
+ANGLES = [1, 50, 90, 117, 180, 360]
 
 @pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
 def dtype(request):
@@ -38,6 +37,13 @@ def img_size(request):
     """
     return request.param
 
+@pytest.fixture(params=ANGLES, ids=lambda x: f"angles={x}", scope="module")
+def num_ang(request):
+    """
+    Angles.
+    """
+    return request.param
+
 
 @pytest.fixture
 def masked_image(dtype, img_size):
@@ -53,12 +59,12 @@ def masked_image(dtype, img_size):
 
 
 # Image.project and compare results to skimage.radon
-def test_image_project(masked_image):
+def test_image_project(masked_image, num_ang):
     """
     Test Image.project on a single stack of images. Compares project method with skimage.
     """
     ny = masked_image.resolution
-    angles = np.linspace(0, 360, ny + 1, endpoint=False)
+    angles = np.linspace(0, 360, num_ang, endpoint=False)
     rads = angles / 180 * np.pi
     s = masked_image.project(rads)
     assert s.shape == (1, len(angles), ny)

From 50be3d3939c72a7131f7ce3ae0a311f6652b0d61 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Thu, 11 Jul 2024 15:37:37 -0400
Subject: [PATCH 21/24] Added Changes from PR: parameterized angles, adjusted
 tests accordingly, renamed variables

---
 tests/test_sinogram.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 94cc172189..69fdd1da37 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -8,7 +8,7 @@
 
 # Relative tolerance comparing line projections to scikit
 # The same tolerance will be used in all scikit comparisons
-SK_TOL = 0.002
+SK_TOL = 0.005
 
 IMG_SIZES = [
     511,
@@ -22,6 +22,7 @@
 
 ANGLES = [1, 50, 90, 117, 180, 360]
 
+
 @pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
 def dtype(request):
     """
@@ -37,6 +38,7 @@ def img_size(request):
     """
     return request.param
 
+
 @pytest.fixture(params=ANGLES, ids=lambda x: f"angles={x}", scope="module")
 def num_ang(request):
     """
@@ -88,7 +90,7 @@ def test_image_project(masked_image, num_ang):
     np.testing.assert_array_less(nrms, SK_TOL, "Error in image projections.")
 
 
-def test_multidim():
+def test_multidim(num_ang):
     """
     Test Image.project on stacks of images. Extension of test_image_project but for multi-dimensional stacks.
     """
@@ -104,12 +106,12 @@ def test_multidim():
     imgs = Image(np.random.random((n, L, L))) * mask
 
     # Generate line project angles
-    angles = np.linspace(0, 180, L, endpoint=False)
+    angles = np.linspace(0, 180, num_ang, endpoint=False)
     rads = angles / 180.0 * np.pi
     s = imgs.project(rads)
 
     # Compare
-    reference_sinograms = np.empty((n, L, L))
+    reference_sinograms = np.empty((n, num_ang, L))
     for i, img in enumerate(imgs):
         # Compute the singleton case, and compare with the stack
         single_sinogram = img.project(rads)

From a8adc0505c473ab5d117f8c7a6abbc3e87c12667 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Thu, 18 Jul 2024 00:46:02 -0400
Subject: [PATCH 22/24] Added extra comments + Integrated Changes from
 lineproject_dbg2 branch

---
 src/aspire/image/image.py |  2 +-
 tests/test_sinogram.py    | 57 ++++++++++++++++++++++-----------------
 2 files changed, 34 insertions(+), 25 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index f6e07dcdc6..613f393993 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -213,7 +213,7 @@ def project(self, angles):
             self.n_images, n_angles, n_real_points
         )
 
-        # Radon transform
+        # 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
diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 69fdd1da37..4dee693351 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -20,7 +20,14 @@
     np.float64,
 ]
 
-ANGLES = [1, 50, 90, 117, 180, 360]
+ANGLES = [
+    1,
+    50,
+    pytest.param(90, marks=pytest.mark.expensive),
+    pytest.param(117, marks=pytest.mark.expensive),
+    pytest.param(180, marks=pytest.mark.expensive),
+    pytest.param(360, marks=pytest.mark.expensive),
+]
 
 
 @pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
@@ -42,7 +49,7 @@ def img_size(request):
 @pytest.fixture(params=ANGLES, ids=lambda x: f"angles={x}", scope="module")
 def num_ang(request):
     """
-    Angles.
+    Angles (Degrees).
     """
     return request.param
 
@@ -50,7 +57,7 @@ def num_ang(request):
 @pytest.fixture
 def masked_image(dtype, img_size):
     """
-    Creates a masked image fixture using camera data from Skikit-Image.
+    Creates a masked image fixture using camera data from Scikit-Image.
     """
     g = grid_2d(img_size, normalized=True, shifted=True)
     mask = g["r"] < 1
@@ -63,7 +70,7 @@ def masked_image(dtype, img_size):
 # Image.project and compare results to skimage.radon
 def test_image_project(masked_image, num_ang):
     """
-    Test Image.project on a single stack of images. Compares project method with skimage.
+    Test Image.project on a single stack of images. Compares project method output with skimage project.
     """
     ny = masked_image.resolution
     angles = np.linspace(0, 360, num_ang, endpoint=False)
@@ -83,8 +90,8 @@ 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] - reference_sinogram) ** 2, axis=-1)) / np.linalg.norm(
+        reference_sinogram, axis=-1
     )
 
     np.testing.assert_array_less(nrms, SK_TOL, "Error in image projections.")
@@ -97,33 +104,35 @@ def test_multidim(num_ang):
 
     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((n, L, L))) * mask
+    imgs = Image(np.random.random((m, n, L, L))) * mask
 
     # Generate line project angles
-    angles = np.linspace(0, 180, num_ang, endpoint=False)
+    angles = np.linspace(0, 360, num_ang, endpoint=False)
     rads = angles / 180.0 * np.pi
     s = imgs.project(rads)
 
     # Compare
-    reference_sinograms = np.empty((n, num_ang, L))
-    for i, img in enumerate(imgs):
-        # Compute the singleton case, and compare with the stack
-        single_sinogram = img.project(rads)
-        # These should be allclose up to determinism in the FFT and NUFFT.
-        np.testing.assert_allclose(s[i : i + 1], single_sinogram)
-
-        # Next individually compute sk's radon transform for each image.
-        reference_sinograms[i] = radon(img._data[0], theta=angles[::-1]).T
-
-    # Compare all lines in each sinogram with sk-image
-    for i in range(n):
-        nrms = np.sqrt(
-            np.mean((s[i] - reference_sinograms[i]) ** 2, axis=1)
-        ) / np.linalg.norm(reference_sinograms[i], axis=1)
-        np.testing.assert_array_less(nrms, SK_TOL, err_msg=f"Error in image {i}.")
+    reference_sinograms = np.empty((m, n, num_ang, 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)
+
+            # These should be allclose up to determinism in the FFT and NUFFT.
+            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
+
+    _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.")

From adfa6e0c5d44ba89f9eb02a2f0eaa2bfdd812dc7 Mon Sep 17 00:00:00 2001
From: Marc Karimi <mk3652@decaf.math.princeton.edu.private>
Date: Fri, 19 Jul 2024 04:10:13 -0400
Subject: [PATCH 23/24] Changed angle fixture description + Id Name

---
 tests/test_sinogram.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/test_sinogram.py b/tests/test_sinogram.py
index 4dee693351..56aa6776e0 100644
--- a/tests/test_sinogram.py
+++ b/tests/test_sinogram.py
@@ -46,10 +46,10 @@ def img_size(request):
     return request.param
 
 
-@pytest.fixture(params=ANGLES, ids=lambda x: f"angles={x}", scope="module")
+@pytest.fixture(params=ANGLES, ids=lambda x: f"n_angles={x}", scope="module")
 def num_ang(request):
     """
-    Angles (Degrees).
+    Number of angles in radon transform.
     """
     return request.param
 

From 7664e9277c5f38b1333bd123e3411266d6fb584a Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 19 Jul 2024 08:01:36 -0400
Subject: [PATCH 24/24] Docstring len cleanup

---
 src/aspire/image/image.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 613f393993..20d998afe6 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -188,9 +188,13 @@ def __init__(self, data, dtype=None):
 
     def project(self, angles):
         """
-        Computes the Radon Transform on an Image Stack using Non-Uniform Fast Fourier Transforms. This method projects the Image stack along different angles and returns the Radon Transform.
+        Computes the Radon Transform on an Image Stack using
+        Non-Uniform Fast Fourier Transforms. This method projects the
+        Image stack along different angles and returns the Radon
+        Transform.
 
-        :param angles: A 1-D Numpy Array of angles in Radians. This is used to compute the Radon Transform at different angles.
+        :param angles: A 1-D Numpy Array of angles in Radians.
+            This is used to compute the Radon Transform at different angles.
         :return: Radon transform of the Image Stack.
         :rtype: Ndarray (stack size, number of angles, image resolution)
         """