From 1ed4de6e7423fc118de49b231d29b0d6c24a585a Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 18 Jul 2024 11:16:28 -0400
Subject: [PATCH 01/24] add pixel_size attribute to Volume class

---
 src/aspire/volume/volume.py | 29 ++++++++++++++++++++---------
 1 file changed, 20 insertions(+), 9 deletions(-)

diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index 7e131a0190..4857e85bb9 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -58,7 +58,7 @@ class Volume:
     Volume is an (N1 x ...) x L x L x L array, along with associated utility methods.
     """
 
-    def __init__(self, data, dtype=None, symmetry_group=None):
+    def __init__(self, data, dtype=None, pixel_size=None, symmetry_group=None):
         """
         A stack of one or more volumes.
 
@@ -76,6 +76,10 @@ def __init__(self, data, dtype=None, symmetry_group=None):
             `(..., resolution, resolution, resolution)`.
         :param dtype: Optionally cast `data` to this dtype.
             Defaults to `data.dtype`.
+        :param pixel_size: Optional voxel_size in Angstroms.
+            When provided will be saved with `map`/`mrc` metadata.
+            Default of `None` will not write to file,
+            but will be considered unit pixels (1) for FSC.
         :param symmetry_group: A SymmetryGroup instance or string indicating symmetry of the Volume.
 
         :return: A Volume instance holding `data`.
@@ -107,6 +111,9 @@ def __init__(self, data, dtype=None, symmetry_group=None):
         self.n_vols = np.prod(self.stack_shape)
         self.resolution = self._data.shape[-1]
         self.size = self._data.size
+        self.pixel_size = None
+        if pixel_size is not None:
+            self.pixel_size = float(pixel_size)
 
         # Set symmetry_group. If None, default to 'C1'.
         self._set_symmetry_group(symmetry_group)
@@ -497,10 +504,14 @@ def downsample(self, ds_res, mask=None, zero_nyquist=True):
         out = fft.centered_ifftn(fx)
         out = out.real * (ds_res**3 / self.resolution**3)
 
+        # Optionally scale pixel size
+        ds_pixel_size = self.pixel_size
+        if ds_pixel_size is not None:
+            ds_pixel_size *= self.resolution / ds_res
+
         # returns a new Volume object
         return self.__class__(
-            xp.asnumpy(out),
-            symmetry_group=self.symmetry_group,
+            xp.asnumpy(out), pixel_size=ds_pixel_size, symmetry_group=self.symmetry_group
         ).stack_reshape(original_stack_shape)
 
     def shift(self):
@@ -592,6 +603,8 @@ def save(self, filename, overwrite=False):
             )
 
         with mrcfile.new(filename, overwrite=overwrite) as mrc:
+            if self.pixel_size is not None:
+                mrc.voxel_size(self.pixel_size)
             mrc.set_data(self._data.astype(np.float32))
 
         if self.dtype != np.float32:
@@ -624,7 +637,7 @@ def load(cls, filename, permissive=True, dtype=None, symmetry_group=None):
 
         return cls(loaded_data, symmetry_group=symmetry_group, dtype=dtype)
 
-    def fsc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
+    def fsc(self, other, cutoff=None, method="fft", plot=False):
         r"""
         Compute the Fourier shell correlation between two volumes.
 
@@ -641,8 +654,6 @@ def fsc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
         :param cutoff: Cutoff value, traditionally `.143`.
             Default `None` implies `cutoff=1` and excludes
             plotting cutoff line.
-        :param pixel_size: Pixel size in angstrom.  Default `None`
-            implies unit in pixels, equivalent to pixel_size=1.
         :param method: Selects either 'fft' (on cartesian grid),
             or 'nufft' (on polar grid). Defaults to 'fft'.
         :param plot: Optionally plot to screen or file.
@@ -662,7 +673,7 @@ def fsc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
         fsc = FourierShellCorrelation(
             a=self.asnumpy(),
             b=other.asnumpy(),
-            pixel_size=pixel_size,
+            pixel_size=self.pixel_size,
             method=method,
         )
 
@@ -681,7 +692,7 @@ def empty_like(v):
         :param v: Volume instance
         :return: Volume instance
         """
-        return Volume(np.empty(v.shape, dtype=v.dtype))
+        return Volume(np.empty(v.shape, dtype=v.dtype), pixel_size=v.pixel_size)
 
     @staticmethod
     def zeros_like(v):
@@ -691,7 +702,7 @@ def zeros_like(v):
         :param v: Volume instance
         :return: Volume instance
         """
-        return Volume(np.zeros(v.shape, dtype=v.dtype))
+        return Volume(np.zeros(v.shape, dtype=v.dtype), pixel_size=v.pixel_size)
 
 
 class CartesianVolume(Volume):

From f4ff4db109d8d9988d6a970e593fb5561d96b218 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 18 Jul 2024 11:25:05 -0400
Subject: [PATCH 02/24] volume load pixel size

---
 src/aspire/volume/volume.py | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index 4857e85bb9..5265fef9b7 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -625,6 +625,12 @@ def load(cls, filename, permissive=True, dtype=None, symmetry_group=None):
         """
         with mrcfile.open(filename, permissive=permissive) as mrc:
             loaded_data = mrc.data
+            pixel_size = mrc.voxel_size
+
+        # Handle default pixel_size.
+        # `mrcfile` defaults to zero, where we use `None`.
+        if pixel_size == 0:
+            pixel_size = None
 
         # FINUFFT work around
         if loaded_data.dtype == np.float32:
@@ -635,7 +641,12 @@ def load(cls, filename, permissive=True, dtype=None, symmetry_group=None):
         if loaded_data.dtype != dtype:
             logger.info(f"{filename} with dtype {loaded_data.dtype} loaded as {dtype}")
 
-        return cls(loaded_data, symmetry_group=symmetry_group, dtype=dtype)
+        return cls(
+            loaded_data,
+            pixel_size=pixel_size,
+            symmetry_group=symmetry_group,
+            dtype=dtype,
+        )
 
     def fsc(self, other, cutoff=None, method="fft", plot=False):
         r"""

From aa9310693f4954afa7d490992bb9692db13baf73 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 18 Jul 2024 12:57:13 -0400
Subject: [PATCH 03/24] volume load pixel size

---
 src/aspire/volume/volume.py | 103 ++++++++++++++++++++++++++++--------
 1 file changed, 81 insertions(+), 22 deletions(-)

diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index 5265fef9b7..315e3be2a0 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -147,7 +147,9 @@ def astype(self, dtype, copy=True):
         :return: Volume instance
         """
         return self.__class__(
-            self.asnumpy().astype(dtype, copy=copy), symmetry_group=self.symmetry_group
+            self.asnumpy().astype(dtype, copy=copy),
+            pixel_size=self.pixel_size,
+            symmetry_group=self.symmetry_group,
         )
 
     def _check_key_dims(self, key):
@@ -158,7 +160,11 @@ def _check_key_dims(self, key):
 
     def __getitem__(self, key):
         self._check_key_dims(key)
-        return self.__class__(self._data[key], symmetry_group=self.symmetry_group)
+        return self.__class__(
+            self._data[key],
+            pixel_size=self.pixel_size,
+            symmetry_group=self.symmetry_group,
+        )
 
     def __setitem__(self, key, value):
         self._check_key_dims(key)
@@ -249,6 +255,7 @@ def stack_reshape(self, *args):
 
         return self.__class__(
             self._data.reshape(*shape, *self._data.shape[-3:]),
+            pixel_size=self.pixel_size,
             symmetry_group=self.symmetry_group,
         )
 
@@ -265,9 +272,15 @@ def __len__(self):
     def __add__(self, other):
         symmetry = self._result_symmetry(other)
         if isinstance(other, Volume):
-            res = self.__class__(self._data + other.asnumpy(), symmetry_group=symmetry)
+            res = self.__class__(
+                self._data + other.asnumpy(),
+                pixel_size=self.pixel_size,
+                symmetry_group=symmetry,
+            )
         else:
-            res = self.__class__(self._data + other, symmetry_group=symmetry)
+            res = self.__class__(
+                self._data + other, pixel_size=self.pixel_size, symmetry_group=symmetry
+            )
 
         return res
 
@@ -277,21 +290,37 @@ def __radd__(self, otherL):
     def __sub__(self, other):
         symmetry = self._result_symmetry(other)
         if isinstance(other, Volume):
-            res = self.__class__(self._data - other.asnumpy(), symmetry_group=symmetry)
+            res = self.__class__(
+                self._data - other.asnumpy(),
+                pixel_size=self.pixel_size,
+                symmetry_group=symmetry,
+            )
         else:
-            res = self.__class__(self._data - other, symmetry_group=symmetry)
+            res = self.__class__(
+                self._data - other, pixel_size=self.pixel_size, symmetry_group=symmetry
+            )
 
         return res
 
     def __rsub__(self, otherL):
-        return self.__class__(otherL - self._data)
+        return self.__class__(
+            otherL - self._data,
+            pixel_size=self.pixel_size,
+            symmetry_group=self.symmetry_group,
+        )
 
     def __mul__(self, other):
         symmetry = self._result_symmetry(other)
         if isinstance(other, Volume):
-            res = self.__class__(self._data * other.asnumpy(), symmetry_group=symmetry)
+            res = self.__class__(
+                self._data * other.asnumpy(),
+                pixel_size=self.pixel_size,
+                symmetry_group=symmetry,
+            )
         else:
-            res = self.__class__(self._data * other, symmetry_group=symmetry)
+            res = self.__class__(
+                self._data * other, pixel_size=self.pixel_size, symmetry_group=symmetry
+            )
 
         return res
 
@@ -304,9 +333,15 @@ def __truediv__(self, other):
         """
         symmetry = self._result_symmetry(other)
         if isinstance(other, Volume):
-            res = self.__class__(self._data / other.asnumpy(), symmetry_group=symmetry)
+            res = self.__class__(
+                self._data / other.asnumpy(),
+                pixel_size=self.pixel_size,
+                symmetry_group=symmetry,
+            )
         else:
-            res = self.__class__(self._data / other, symmetry_group=symmetry)
+            res = self.__class__(
+                self._data / other, pixel_size=self.pixel_size, symmetry_group=symmetry
+            )
 
         return res
 
@@ -314,7 +349,10 @@ def __rtruediv__(self, otherL):
         """
         Right scalar division, follows numpy semantics.
         """
-        return otherL * Volume(1.0 / self._data)
+        return otherL * Volume(
+            1.0 / self._data,
+            pixel_size=self.pixel_size,
+        )
 
     def project(self, rot_matrices, zero_nyquist=True):
         """
@@ -381,7 +419,7 @@ def project(self, rot_matrices, zero_nyquist=True):
             im_f[:, :, 0] = 0
 
         im_f = fft.centered_ifft2(im_f)
-
+        # todo add pixel_size to Image
         return aspire.image.Image(xp.asnumpy(im_f.real))
 
     def to_vec(self):
@@ -426,7 +464,7 @@ def transpose(self):
         v = self._data.reshape(-1, *self._data.shape[-3:])
         vt = np.transpose(v, (0, -1, -2, -3))
         vt = vt.reshape(*original_stack_shape, *self._data.shape[-3:])
-        return self.__class__(vt, symmetry_group=symmetry)
+        return self.__class__(vt, pixel_size=self.pixel_size, symmetry_group=symmetry)
 
     @property
     def T(self):
@@ -469,7 +507,11 @@ def flip(self, axis=-3):
                     f"Cannot flip axis {ax}: stack axis. Did you mean {ax-4}?"
                 )
 
-        return self.__class__(np.flip(self._data, axis), symmetry_group=symmetry)
+        return self.__class__(
+            np.flip(self._data, axis),
+            pixel_size=self.pixel_size,
+            symmetry_group=symmetry,
+        )
 
     def downsample(self, ds_res, mask=None, zero_nyquist=True):
         """
@@ -583,7 +625,7 @@ def rotate(self, rot_matrices, zero_nyquist=True):
             np.real(fft.centered_ifftn(xp.asarray(vol_f), axes=(-3, -2, -1)))
         )
 
-        return self.__class__(vol, symmetry_group=symmetry)
+        return self.__class__(vol, pixel_size=self.pixel_size, symmetry_group=symmetry)
 
     def denoise(self):
         raise NotImplementedError
@@ -625,12 +667,7 @@ def load(cls, filename, permissive=True, dtype=None, symmetry_group=None):
         """
         with mrcfile.open(filename, permissive=permissive) as mrc:
             loaded_data = mrc.data
-            pixel_size = mrc.voxel_size
-
-        # Handle default pixel_size.
-        # `mrcfile` defaults to zero, where we use `None`.
-        if pixel_size == 0:
-            pixel_size = None
+            pixel_size = Volume._vx_array_to_size(mrc.voxel_size)
 
         # FINUFFT work around
         if loaded_data.dtype == np.float32:
@@ -715,6 +752,28 @@ def zeros_like(v):
         """
         return Volume(np.zeros(v.shape, dtype=v.dtype), pixel_size=v.pixel_size)
 
+    @staticmethod
+    def _vx_array_to_size(vx):
+        """
+        Utility to convert from several possible `mrcfile.voxel_size` representations to a single (float) value or None.
+        """
+        # Convert from recarray to single values,
+        #   checks uniformity.
+        if isinstance(vx, np.recarray):
+            if vx.x != vx.y or vx.x != vx.z:
+                raise ValueError(f"Voxel sizes are not uniform: {vx}")
+            vx = vx.x
+
+        # Convert `0` to `None`
+        if (isinstance(vx, int) or isinstance(vx, float)) and vx == 0:
+            vx = None
+
+        # Consistently return a `float` when not None
+        if vx is not None:
+            vx = float(vx)
+
+        return vx
+
 
 class CartesianVolume(Volume):
     def expand(self, basis):

From f0ad627625199bf6d759d0e79d330451490d3a96 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 18 Jul 2024 13:55:27 -0400
Subject: [PATCH 04/24] volume load pixel size

---
 src/aspire/volume/volume.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index 315e3be2a0..2dc0c9a309 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -645,9 +645,10 @@ def save(self, filename, overwrite=False):
             )
 
         with mrcfile.new(filename, overwrite=overwrite) as mrc:
-            if self.pixel_size is not None:
-                mrc.voxel_size(self.pixel_size)
             mrc.set_data(self._data.astype(np.float32))
+            # Note assigning voxel_size must come after `set_data`
+            if self.pixel_size is not None:
+                mrc.voxel_size = self.pixel_size
 
         if self.dtype != np.float32:
             logger.info(f"Volume with dtype {self.dtype} saved with dtype float32")
@@ -757,6 +758,7 @@ def _vx_array_to_size(vx):
         """
         Utility to convert from several possible `mrcfile.voxel_size` representations to a single (float) value or None.
         """
+
         # Convert from recarray to single values,
         #   checks uniformity.
         if isinstance(vx, np.recarray):
@@ -765,7 +767,9 @@ def _vx_array_to_size(vx):
             vx = vx.x
 
         # Convert `0` to `None`
-        if (isinstance(vx, int) or isinstance(vx, float)) and vx == 0:
+        if (
+            isinstance(vx, int) or isinstance(vx, float) or isinstance(vx, np.ndarray)
+        ) and vx == 0:
             vx = None
 
         # Consistently return a `float` when not None

From 5195b2de799f6dc72defc8c7264b3bd8d36bb982 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 18 Jul 2024 14:14:44 -0400
Subject: [PATCH 05/24] rm old pixel_size attr from FSC calls

---
 tests/test_fourier_correlation.py     | 18 ++++++-------
 tests/test_mean_estimator_boosting.py |  2 +-
 tests/test_volume.py                  | 39 +++++++++++++++++++++++++--
 3 files changed, 46 insertions(+), 13 deletions(-)

diff --git a/tests/test_fourier_correlation.py b/tests/test_fourier_correlation.py
index 282089ce12..c53ba40ae1 100644
--- a/tests/test_fourier_correlation.py
+++ b/tests/test_fourier_correlation.py
@@ -178,7 +178,7 @@ def test_frc_plot(image_fixture, method):
 def test_fsc_id(volume_fixture, method):
     vol, _ = volume_fixture
 
-    fsc_resolution, fsc = vol.fsc(vol, pixel_size=1, cutoff=0.143, method=method)
+    fsc_resolution, fsc = vol.fsc(vol, cutoff=0.143, method=method)
     assert np.isclose(fsc_resolution[0], 2, rtol=0.02)
     assert np.allclose(fsc, 1, rtol=0.01)
 
@@ -186,11 +186,11 @@ def test_fsc_id(volume_fixture, method):
 def test_fsc_trunc(volume_fixture, method):
     vol_a, vol_b = volume_fixture
 
-    fsc_resolution, fsc = vol_a.fsc(vol_b, pixel_size=1, cutoff=0.143, method=method)
+    fsc_resolution, fsc = vol_a.fsc(vol_b, cutoff=0.143, method=method)
     assert fsc_resolution[0] > 3.0
 
     # The follow should correspond to the test_fsc_plot below.
-    fsc_resolution, fsc = vol_a.fsc(vol_b, pixel_size=1, cutoff=0.5, method=method)
+    fsc_resolution, fsc = vol_a.fsc(vol_b, cutoff=0.5, method=method)
     assert fsc_resolution[0] > 3.9
 
 
@@ -202,13 +202,13 @@ def test_fsc_vol_plot(volume_fixture):
 
     # Plot to screen
     with matplotlib_no_gui():
-        _ = vol_a.fsc(vol_b, pixel_size=1, cutoff=0.5, plot=True)
+        _ = vol_a.fsc(vol_b, cutoff=0.5, plot=True)
 
     # Plot to file
     # Also tests `cutoff=None`
     with tempfile.TemporaryDirectory() as tmp_input_dir:
         file_path = os.path.join(tmp_input_dir, "vol_fsc_curve.png")
-        vol_a.fsc(vol_b, pixel_size=1, cutoff=None, plot=file_path)
+        vol_a.fsc(vol_b, cutoff=None, plot=file_path)
         assert os.path.exists(file_path)
 
 
@@ -218,9 +218,7 @@ def test_fsc_plot(volume_fixture, method):
     """
     vol_a, vol_b = volume_fixture
 
-    fsc = FourierShellCorrelation(
-        vol_a.asnumpy(), vol_b.asnumpy(), pixel_size=1, method=method
-    )
+    fsc = FourierShellCorrelation(vol_a.asnumpy(), vol_b.asnumpy(), method=method)
 
     with matplotlib_no_gui():
         fsc.plot(cutoff=0.5)
@@ -314,7 +312,7 @@ def test_vol_type_mismatch():
     b = a.asnumpy()
 
     with pytest.raises(TypeError, match=r"`other` volume must be an `Volume` instance"):
-        _ = a.fsc(b, pixel_size=1, cutoff=0.143)
+        _ = a.fsc(b, cutoff=0.143)
 
 
 # Broadcasting
@@ -378,7 +376,7 @@ def test_fsc_id_bcast(volume_fixture, method):
     k = 3
     vol_b = Volume(np.tile(vol.asnumpy(), (3, 1, 1, 1)))
 
-    fsc_resolution, fsc = vol.fsc(vol_b, pixel_size=1, cutoff=0.143, method=method)
+    fsc_resolution, fsc = vol.fsc(vol_b, cutoff=0.143, method=method)
     assert np.allclose(
         fsc_resolution,
         [
diff --git a/tests/test_mean_estimator_boosting.py b/tests/test_mean_estimator_boosting.py
index 9251dee09e..6eac159115 100644
--- a/tests/test_mean_estimator_boosting.py
+++ b/tests/test_mean_estimator_boosting.py
@@ -122,7 +122,7 @@ def weighted_source(weighted_volume):
 def test_fsc(source, estimated_volume):
     """Compare estimated volume to source volume with FSC."""
     # Fourier Shell Correlation
-    fsc_resolution, fsc = source.vols.fsc(estimated_volume, pixel_size=1, cutoff=0.5)
+    fsc_resolution, fsc = source.vols.fsc(estimated_volume, cutoff=0.5)
 
     # Check that resolution is less than 2.1 pixels.
     np.testing.assert_array_less(fsc_resolution, 2.1)
diff --git a/tests/test_volume.py b/tests/test_volume.py
index 991cc3288e..bec31d6ed0 100644
--- a/tests/test_volume.py
+++ b/tests/test_volume.py
@@ -75,7 +75,7 @@ def vols_1(data_1):
 
 @pytest.fixture
 def vols_2(data_2):
-    return Volume(data_2)
+    return Volume(data_2, pixel_size=4.56)
 
 
 @pytest.fixture
@@ -291,6 +291,32 @@ def test_save_load(vols_1):
         assert np.allclose(vols_1, vols_loaded_single)
         assert isinstance(vols_loaded_double, Volume)
         assert np.allclose(vols_1, vols_loaded_double)
+        assert vols_loaded_single.pixel_size is None, "Pixel size should be None"
+        assert vols_loaded_double.pixel_size is None, "Pixel size should be None"
+
+
+def test_save_load_pixel_size(vols_2):
+    # Create a tmpdir in a context. It will be cleaned up on exit.
+    with tempfile.TemporaryDirectory() as tmpdir:
+        # Save the Volume object into an MRC files
+        mrcs_filepath = os.path.join(tmpdir, "test.mrc")
+        vols_2.save(mrcs_filepath)
+
+        # Load saved MRC file as a Volume of dtypes single and double.
+        vols_loaded_single = Volume.load(mrcs_filepath, dtype=np.float32)
+        vols_loaded_double = Volume.load(mrcs_filepath, dtype=np.float64)
+
+        # Confirm the pixel size is loaded
+        np.testing.assert_approx_equal(
+            vols_loaded_single.pixel_size,
+            vols_2.pixel_size,
+            err_msg="Incorrect pixel size in singles.",
+        )
+        np.testing.assert_approx_equal(
+            vols_loaded_double.pixel_size,
+            vols_2.pixel_size,
+            err_msg="Incorrect pixel size in doubles.",
+        )
 
 
 def test_project(vols_hot_cold):
@@ -545,11 +571,20 @@ def test_flip(vols_1, data_1):
 
 
 def test_downsample(res):
-    vols = Volume(np.load(os.path.join(DATA_DIR, "clean70SRibosome_vol.npy")))
+    vols = Volume(
+        np.load(os.path.join(DATA_DIR, "clean70SRibosome_vol.npy")), pixel_size=1.23
+    )
     result = vols.downsample(res)
     og_res = vols.resolution
     ds_res = result.resolution
 
+    # Confirm the pixel size is scaled
+    np.testing.assert_approx_equal(
+        result.pixel_size,
+        vols.pixel_size * res / ds_res,
+        err_msg="Incorrect pixel size.",
+    )
+
     # check signal energy
     np.testing.assert_allclose(
         anorm(vols.asnumpy(), axes=(1, 2, 3)) / og_res,

From 3b2a6a64d0b2195abb9c27c596fc20da6bffb5e2 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 18 Jul 2024 14:56:26 -0400
Subject: [PATCH 06/24] Add Image.voxel_size handling

---
 src/aspire/image/image.py         | 115 +++++++++++++++++++++++-------
 src/aspire/volume/volume.py       |   5 +-
 tests/test_fourier_correlation.py |  26 ++++---
 3 files changed, 104 insertions(+), 42 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index cbab268cf6..984b1e6540 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -80,7 +80,7 @@ def load_mrc(filepath):
     Load raw data from `.mrc` into an array.
 
     :param filepath: File path (string).
-    :return: numpy array of image data.
+    :return: (numpy array of image data, pixel_size)
     """
 
     # mrcfile tends to yield many warnings about EMPIAR datasets being corrupt
@@ -92,6 +92,7 @@ def load_mrc(filepath):
 
         with mrcfile.open(filepath, mode="r", permissive=True) as mrc:
             im = mrc.data
+            pixel_size = Image._vx_array_to_size(mrc.voxel_size)
 
         # Log each mrcfile warning to debug log, noting the associated file
         for w in ws:
@@ -110,19 +111,29 @@ def load_mrc(filepath):
                 f" Will attempt to continue processing {filepath}"
             )
 
-    return im
+    return im, pixel_size
 
 
 def load_tiff(filepath):
     """
     Load raw data from `.tiff` into an array.
 
+    Note, TIFF does not natively provide equivalent to pixel/voxel_size,
+    so users of TIFF files may need to manually assign `pixel_size` to
+    `Image` instances when required. Defaults to `pixel_size=None`.
+
     :param filepath: File path (string).
-    :return: numpy array of image data.
+    :return: (numpy array of image data, pixel_size=None)
     """
+    # Use PIL to open `filepath`
+    img = PILImage.open(filepath)
+
+    # Future todo, extract `voxel_size` if available in TIFF tags (custom tag?)
+    # For now, default to `None`.
+    voxel_size = None
 
-    # Use PIL to open `filepath` and cast to numpy array.
-    return np.array(PILImage.open(filepath))
+    # Cast image data as numpy array
+    return np.array(img), voxel_size
 
 
 class Image:
@@ -133,7 +144,7 @@ class Image:
         ".tiff": load_tiff,
     }
 
-    def __init__(self, data, dtype=None):
+    def __init__(self, data, pixel_size=None, dtype=None):
         """
         A stack of one or more images.
 
@@ -149,6 +160,10 @@ def __init__(self, data, dtype=None):
 
         :param data: Numpy array containing image data with shape
             `(..., resolution, resolution)`.
+        :param pixel_size: Optional pixel size in Angstroms.
+            When provided will be saved with `mrc` metadata.
+            Default of `None` will not write to file,
+            but will be considered unit pixels (1) for FSC.
         :param dtype: Optionally cast `data` to this dtype.
             Defaults to `data.dtype`.
 
@@ -180,6 +195,9 @@ def __init__(self, data, dtype=None):
         self.stack_shape = self._data.shape[:-2]
         self.n_images = np.prod(self.stack_shape)
         self.resolution = self._data.shape[-1]
+        self.pixel_size = None
+        if pixel_size is not None:
+            self.pixel_size = float(pixel_size)
 
         # Numpy interop
         # https://numpy.org/devdocs/user/basics.interoperability.html#the-array-interface-protocol
@@ -238,7 +256,7 @@ def _check_key_dims(self, key):
 
     def __getitem__(self, key):
         self._check_key_dims(key)
-        return self.__class__(self._data[key])
+        return self.__class__(self._data[key], pixel_size=self.pixel_size)
 
     def __setitem__(self, key, value):
         self._check_key_dims(key)
@@ -266,31 +284,34 @@ def stack_reshape(self, *args):
                 f"Number of images {self.n_images} cannot be reshaped to {shape}."
             )
 
-        return self.__class__(self._data.reshape(*shape, *self._data.shape[-2:]))
+        return self.__class__(
+            self._data.reshape(*shape, *self._data.shape[-2:]),
+            pixel_size=self.pixel_size,
+        )
 
     def __add__(self, other):
         if isinstance(other, Image):
             other = other._data
 
-        return self.__class__(self._data + other)
+        return self.__class__(self._data + other, pixel_size=self.pixel_size)
 
     def __sub__(self, other):
         if isinstance(other, Image):
             other = other._data
 
-        return self.__class__(self._data - other)
+        return self.__class__(self._data - other, pixel_size=self.pixel_size)
 
     def __mul__(self, other):
         if isinstance(other, Image):
             other = other._data
 
-        return self.__class__(self._data * other)
+        return self.__class__(self._data * other, pixel_size=self.pixel_size)
 
     def __neg__(self):
-        return self.__class__(-self._data)
+        return self.__class__(-self._data, pixel_size=self.pixel_size)
 
     def sqrt(self):
-        return self.__class__(np.sqrt(self._data))
+        return self.__class__(np.sqrt(self._data), pixel_size=self.pixel_size)
 
     @property
     def T(self):
@@ -312,7 +333,9 @@ def transpose(self):
         im = self.stack_reshape(-1)
         imt = np.transpose(im._data, (0, -1, -2))
 
-        return self.__class__(imt).stack_reshape(original_stack_shape)
+        return self.__class__(imt, pixel_size=self.pixel_size).stack_reshape(
+            original_stack_shape
+        )
 
     def flip(self, axis=-2):
         """
@@ -335,7 +358,7 @@ def flip(self, axis=-2):
                     f"Cannot flip axis {ax}: stack axis. Did you mean {ax-3}?"
                 )
 
-        return self.__class__(np.flip(self._data, axis))
+        return self.__class__(np.flip(self._data, axis), pixel_size=self.pixel_size)
 
     def __repr__(self):
         msg = f"{self.n_images} {self.dtype} images arranged as a {self.stack_shape} stack"
@@ -355,7 +378,7 @@ def asnumpy(self):
         return view
 
     def copy(self):
-        return self.__class__(self._data.copy())
+        return self.__class__(self._data.copy(), pixel_size=self.pixel_size)
 
     def shift(self, shifts):
         """
@@ -412,7 +435,14 @@ def downsample(self, ds_res, zero_nyquist=True):
         out = fft.centered_ifft2(crop_fx).real * (ds_res**2 / self.resolution**2)
         out = xp.asnumpy(out)
 
-        return self.__class__(out).stack_reshape(original_stack_shape)
+        # Optionally scale pixel size
+        ds_pixel_size = self.pixel_size
+        if ds_pixel_size is not None:
+            ds_pixel_size *= self.resolution / ds_res
+
+        return self.__class__(out, pixel_size=ds_pixel_size).stack_reshape(
+            original_stack_shape
+        )
 
     def filter(self, filter):
         """
@@ -441,7 +471,9 @@ def filter(self, filter):
 
         im = xp.asnumpy(im.real)
 
-        return self.__class__(im).stack_reshape(original_stack_shape)
+        return self.__class__(im, pixel_size=self.pixel_size).stack_reshape(
+            original_stack_shape
+        )
 
     def rotate(self):
         raise NotImplementedError
@@ -453,6 +485,9 @@ def save(self, mrcs_filepath, overwrite=False):
         with mrcfile.new(mrcs_filepath, overwrite=overwrite) as mrc:
             # original input format (the image index first)
             mrc.set_data(self._data.astype(np.float32))
+            # Note assigning voxel_size must come after `set_data`
+            if self.pixel_size is not None:
+                mrc.voxel_size = self.pixel_size
 
     @staticmethod
     def load(filepath, dtype=None):
@@ -477,14 +512,14 @@ def load(filepath, dtype=None):
             )
 
         # Call the appropriate file reader
-        im = Image.extensions[ext](filepath)
+        im, pixel_size = 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)
+        return Image(im, pixel_size=pixel_size)
 
     def _im_translate(self, shifts):
         """
@@ -535,7 +570,9 @@ def _im_translate(self, shifts):
         im_translated = xp.asnumpy(im_translated.real)
 
         # Reshape to stack shape
-        return self.__class__(im_translated).stack_reshape(stack_shape)
+        return self.__class__(im_translated, pixel_size=self.pixel_size).stack_reshape(
+            stack_shape
+        )
 
     def norm(self):
         return anorm(self._data)
@@ -602,7 +639,9 @@ def backproject(self, rot_matrices, symmetry_group=None, zero_nyquist=True):
 
         vol /= L
 
-        return aspire.volume.Volume(vol, symmetry_group=symmetry_group)
+        return aspire.volume.Volume(
+            vol, pixel_size=self.pixel_size, symmetry_group=symmetry_group
+        )
 
     def show(self, columns=5, figsize=(20, 10), colorbar=True):
         """
@@ -645,7 +684,7 @@ def show(self, columns=5, figsize=(20, 10), colorbar=True):
 
             plt.show()
 
-    def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
+    def frc(self, other, cutoff=None, method="fft", plot=False):
         r"""
         Compute the Fourier ring correlation between two images.
 
@@ -663,8 +702,6 @@ def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
             Default `None` implies `cutoff=1` and excludes
             plotting cutoff line.
 
-        :param pixel_size: Pixel size in angstrom.  Default `None`
-            implies unit in pixels, equivalent to pixel_size=1.
         :param method: Selects either 'fft' (on cartesian grid),
             or 'nufft' (on polar grid). Defaults to 'fft'.
         :param plot: Optionally plot to screen or file.
@@ -684,7 +721,7 @@ def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
         frc = FourierRingCorrelation(
             a=self.asnumpy(),
             b=other.asnumpy(),
-            pixel_size=pixel_size,
+            pixel_size=self.pixel_size,
             method=method,
         )
 
@@ -695,6 +732,32 @@ def frc(self, other, cutoff=None, pixel_size=None, method="fft", plot=False):
 
         return frc.analyze_correlations(cutoff), frc.correlations
 
+    @staticmethod
+    def _vx_array_to_size(vx):
+        """
+        Utility to convert from several possible `mrcfile.voxel_size`
+        representations to a single (float) value or None.
+        """
+
+        # Convert from recarray to single values,
+        #   checks uniformity.
+        if isinstance(vx, np.recarray):
+            if vx.x != vx.y:
+                raise ValueError(f"Voxel sizes are not uniform: {vx}")
+            vx = vx.x
+
+        # Convert `0` to `None`
+        if (
+            isinstance(vx, int) or isinstance(vx, float) or isinstance(vx, np.ndarray)
+        ) and vx == 0:
+            vx = None
+
+        # Consistently return a `float` when not None
+        if vx is not None:
+            vx = float(vx)
+
+        return vx
+
 
 class CartesianImage(Image):
     def expand(self, basis):
diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index 2dc0c9a309..1484a9c344 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -420,7 +420,7 @@ def project(self, rot_matrices, zero_nyquist=True):
 
         im_f = fft.centered_ifft2(im_f)
         # todo add pixel_size to Image
-        return aspire.image.Image(xp.asnumpy(im_f.real))
+        return aspire.image.Image(xp.asnumpy(im_f.real), pixel_size=self.pixel_size)
 
     def to_vec(self):
         """Returns an N x resolution ** 3 array."""
@@ -756,7 +756,8 @@ def zeros_like(v):
     @staticmethod
     def _vx_array_to_size(vx):
         """
-        Utility to convert from several possible `mrcfile.voxel_size` representations to a single (float) value or None.
+        Utility to convert from several possible `mrcfile.voxel_size`
+        representations to a single (float) value or None.
         """
 
         # Convert from recarray to single values,
diff --git a/tests/test_fourier_correlation.py b/tests/test_fourier_correlation.py
index c53ba40ae1..79240572f8 100644
--- a/tests/test_fourier_correlation.py
+++ b/tests/test_fourier_correlation.py
@@ -115,7 +115,7 @@ def volume_fixture(img_size, dtype):
 def test_frc_id(image_fixture, method):
     img, _, _ = image_fixture
 
-    frc_resolution, frc = img.frc(img, pixel_size=1, cutoff=0.143, method=method)
+    frc_resolution, frc = img.frc(img, cutoff=0.143, method=method)
     assert np.isclose(frc_resolution[0], 2, rtol=0.02)
     assert np.allclose(frc, 1, rtol=0.01)
 
@@ -123,14 +123,14 @@ def test_frc_id(image_fixture, method):
 def test_frc_trunc(image_fixture, method):
     img_a, img_b, _ = image_fixture
     assert img_a.dtype == img_b.dtype
-    frc_resolution, frc = img_a.frc(img_b, pixel_size=1, cutoff=0.143, method=method)
+    frc_resolution, frc = img_a.frc(img_b, cutoff=0.143, method=method)
     assert frc_resolution[0] > 3.0
 
 
 def test_frc_noise(image_fixture, method):
     img_a, _, img_n = image_fixture
 
-    frc_resolution, frc = img_a.frc(img_n, pixel_size=1, cutoff=0.143, method=method)
+    frc_resolution, frc = img_a.frc(img_n, cutoff=0.143, method=method)
     assert frc_resolution[0] > 3.5
 
 
@@ -142,13 +142,13 @@ def test_frc_img_plot(image_fixture):
 
     # Plot to screen
     with matplotlib_no_gui():
-        _ = img_a.frc(img_n, pixel_size=1, cutoff=0.143, plot=True)
+        _ = img_a.frc(img_n, cutoff=0.143, plot=True)
 
     # Plot to file
     # Also tests `cutoff=None`
     with tempfile.TemporaryDirectory() as tmp_input_dir:
         file_path = os.path.join(tmp_input_dir, "img_frc_curve.png")
-        img_a.frc(img_n, pixel_size=1, cutoff=None, plot=file_path)
+        img_a.frc(img_n, cutoff=None, plot=file_path)
         assert os.path.exists(file_path)
 
 
@@ -160,9 +160,7 @@ def test_frc_plot(image_fixture, method):
     """
     img_a, img_b, _ = image_fixture
 
-    frc = FourierRingCorrelation(
-        img_a.asnumpy(), img_b.asnumpy(), pixel_size=1, method=method
-    )
+    frc = FourierRingCorrelation(img_a.asnumpy(), img_b.asnumpy(), method=method)
 
     with matplotlib_no_gui():
         frc.plot(cutoff=0.5)
@@ -304,7 +302,7 @@ def test_img_type_mismatch():
     b = a.asnumpy()
 
     with pytest.raises(TypeError, match=r"`other` image must be an `Image` instance"):
-        _ = a.frc(b, pixel_size=1, cutoff=0.143)
+        _ = a.frc(b, cutoff=0.143)
 
 
 def test_vol_type_mismatch():
@@ -327,7 +325,7 @@ def test_frc_id_bcast(image_fixture, method):
     k = 3
     img_b = Image(np.tile(img, (3, 1, 1)))
 
-    frc_resolution, frc = img.frc(img_b, pixel_size=1, cutoff=0.143, method=method)
+    frc_resolution, frc = img.frc(img_b, cutoff=0.143, method=method)
     assert np.allclose(
         frc_resolution,
         [
@@ -342,7 +340,7 @@ def test_frc_id_bcast(image_fixture, method):
     # (1) x (1,3)
     img_b = img_b.stack_reshape(1, 3)
 
-    frc_resolution, frc = img.frc(img_b, pixel_size=1, cutoff=0.143, method=method)
+    frc_resolution, frc = img.frc(img_b, cutoff=0.143, method=method)
     assert np.allclose(
         frc_resolution,
         [
@@ -357,7 +355,7 @@ def test_frc_id_bcast(image_fixture, method):
     # (1) x (3,1)
     img_b = img_b.stack_reshape(3, 1)
 
-    frc_resolution, frc = img.frc(img_b, pixel_size=1, cutoff=0.143, method=method)
+    frc_resolution, frc = img.frc(img_b, cutoff=0.143, method=method)
     assert np.allclose(
         frc_resolution,
         [
@@ -398,12 +396,12 @@ def test_frc_img_plot_bcast(image_fixture):
 
     # Plot to screen, one:many
     with matplotlib_no_gui():
-        _ = img_a.frc(img_b, pixel_size=1, cutoff=0.143, plot=True)
+        _ = img_a.frc(img_b, cutoff=0.143, plot=True)
 
     # Plot to file, many elementwise
     with tempfile.TemporaryDirectory() as tmp_input_dir:
         file_path = os.path.join(tmp_input_dir, "img_frc_curve.png")
-        img_b.frc(img_b, pixel_size=1, cutoff=0.143, plot=file_path)
+        img_b.frc(img_b, cutoff=0.143, plot=file_path)
         assert os.path.exists(file_path)
 
 

From 78adb4ad99cc6bb618b14cf0ee0b1ad23326a662 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 22 Jul 2024 09:34:39 -0400
Subject: [PATCH 07/24] cleanup image tests round 1

---
 tests/test_image.py | 160 ++++++++++++++++++++++----------------------
 1 file changed, 81 insertions(+), 79 deletions(-)

diff --git a/tests/test_image.py b/tests/test_image.py
index 688d4169ec..6a9888479a 100644
--- a/tests/test_image.py
+++ b/tests/test_image.py
@@ -24,8 +24,22 @@
 n = 3
 mdim = 2
 
+PARITY = [0, 1]
+DTYPES = [np.float32, np.float64]
 
-def get_images(parity=0, dtype=np.float32):
+
+@pytest.fixture(params=PARITY, ids=lambda x: f"parity={x}", scope="module")
+def parity(request):
+    return request.param
+
+
+@pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
+def dtype(request):
+    return request.param
+
+
+@pytest.fixture(scope="module")
+def get_images(parity, dtype):
     size = 768 - parity
     # numpy array for top-level functions that directly expect it
     im_np = face(gray=True).astype(dtype)[np.newaxis, :size, :size]
@@ -33,36 +47,40 @@ def get_images(parity=0, dtype=np.float32):
     im_np /= denom  # Normalize test image data to 0,1
 
     # Independent Image object for testing Image methods
-    im = Image(im_np.copy())
+    im = Image(im_np.copy(), pixel_size=1.23)
     return im_np, im
 
 
-def get_stacks(parity=0, dtype=np.float32):
-    im_np, im = get_images(parity, dtype)
+@pytest.fixture(scope="module")
+def get_stacks(get_images, dtype):
+    im_np, im = get_images
 
     # Construct a simple stack of Images
-    ims_np = np.empty((n, *im_np.shape[1:]), dtype=dtype)
+    ims_np = np.empty((n, *im_np.shape[1:]), dtype=im_np.dtype)
     for i in range(n):
         ims_np[i] = im_np * (i + 1) / float(n)
 
     # Independent Image stack object for testing Image methods
-    ims = Image(ims_np)
+    ims = Image(ims_np.copy())
     return ims_np, ims
 
 
-def get_mdim_images(parity=0, dtype=np.float32):
-    ims_np, im = get_stacks(parity, dtype)
+# Note that `get_mdim_images` is mutated by some tests,
+# force per function scope.
+@pytest.fixture(scope="function")
+def get_mdim_images(get_stacks):
+    ims_np, im = get_stacks
     # Multi dimensional stack Image object
     mdim = 2
     mdim_ims_np = np.concatenate([ims_np] * mdim).reshape(mdim, *ims_np.shape)
 
     # Independent multidimensional Image stack object for testing Image methods
-    mdim_ims = Image(mdim_ims_np)
+    mdim_ims = Image(mdim_ims_np.copy())
     return mdim_ims_np, mdim_ims
 
 
-def testRepr():
-    _, mdim_ims = get_mdim_images()
+def testRepr(get_mdim_images):
+    _, mdim_ims = get_mdim_images
     r = repr(mdim_ims)
     logger.info(f"Image repr:\n{r}")
 
@@ -73,9 +91,8 @@ def testNonSquare():
         _ = Image(np.empty((4, 5)))
 
 
-@pytest.mark.parametrize("parity,dtype", params)
-def testImShift(parity, dtype):
-    im_np, im = get_images(parity, dtype)
+def testImShift(get_images, dtype):
+    im_np, im = get_images
     # Note that the _im_translate method can handle float input shifts, as it
     # computes the shifts in Fourier space, rather than performing a roll
     # However, NumPy's roll() only accepts integer inputs
@@ -101,10 +118,8 @@ def testImShift(parity, dtype):
     np.testing.assert_allclose(im0.asnumpy()[0, :, :], im3, atol=atol)
 
 
-@pytest.mark.parametrize("parity,dtype", params)
-def testImShiftStack(parity, dtype):
-    ims_np, ims = get_stacks(parity, dtype)
-
+def testImShiftStack(get_stacks, dtype):
+    ims_np, ims = get_stacks
     # test stack of shifts (same number as Image.num_img)
     # mix of odd and even
     shifts = np.array([[100, 200], [203, 150], [55, 307]])
@@ -131,8 +146,8 @@ def testImShiftStack(parity, dtype):
     np.testing.assert_allclose(im0.asnumpy(), im3, atol=atol)
 
 
-def testImageShiftErrors():
-    _, im = get_images(0, np.float32)
+def testImageShiftErrors(get_images):
+    _, im = get_images
     # test bad shift shape
     with pytest.raises(ValueError, match="Input shifts must be of shape"):
         _ = im.shift(np.array([100, 100, 100]))
@@ -141,18 +156,16 @@ def testImageShiftErrors():
         _ = im.shift(np.array([[100, 200], [100, 200]]))
 
 
-@pytest.mark.parametrize("parity,dtype", params)
-def testImageSqrt(parity, dtype):
-    im_np, im = get_images(parity, dtype)
-    ims_np, ims = get_stacks(parity, dtype)
+def testImageSqrt(get_images, get_stacks):
+    im_np, im = get_images
+    ims_np, ims = get_stacks
     assert np.allclose(im.sqrt().asnumpy(), np.sqrt(im_np))
     assert np.allclose(ims.sqrt().asnumpy(), np.sqrt(ims_np))
 
 
-@pytest.mark.parametrize("parity,dtype", params)
-def testImageTranspose(parity, dtype):
-    im_np, im = get_images(parity, dtype)
-    ims_np, ims = get_stacks(parity, dtype)
+def testImageTranspose(get_images, get_stacks):
+    im_np, im = get_images
+    ims_np, ims = get_stacks
     # test method and abbreviation
     assert np.allclose(im.T.asnumpy(), np.transpose(im_np, (0, 2, 1)))
     assert np.allclose(im.transpose().asnumpy(), np.transpose(im_np, (0, 2, 1)))
@@ -163,10 +176,9 @@ def testImageTranspose(parity, dtype):
         assert np.allclose(ims.transpose()[i], ims_np[i].T)
 
 
-@pytest.mark.parametrize("parity,dtype", params)
-def testImageFlip(parity, dtype):
-    im_np, im = get_images(parity, dtype)
-    ims_np, ims = get_stacks(parity, dtype)
+def testImageFlip(get_images, get_stacks):
+    im_np, im = get_images
+    ims_np, ims = get_stacks
     for axis in powerset(range(1, 3)):
         if not axis:
             # test default
@@ -188,31 +200,31 @@ def testImageFlip(parity, dtype):
             _ = im.flip(axis)
 
 
-def testShape():
-    ims_np, ims = get_stacks()
+def testShape(get_stacks):
+    ims_np, ims = get_stacks
     assert ims.shape == ims_np.shape
     assert ims.stack_shape == ims_np.shape[:-2]
     assert ims.stack_ndim == 1
 
 
-def testMultiDimShape():
-    ims_np, ims = get_stacks()
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimShape(get_stacks, get_mdim_images):
+    ims_np, ims = get_stacks
+    mdim_ims_np, mdim_ims = get_mdim_images
     assert mdim_ims.shape == mdim_ims_np.shape
     assert mdim_ims.stack_shape == mdim_ims_np.shape[:-2]
     assert mdim_ims.stack_ndim == mdim
     assert mdim_ims.n_images == mdim * ims.n_images
 
 
-def testBadKey():
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testBadKey(get_mdim_images):
+    mdim_ims_np, mdim_ims = get_mdim_images
     with pytest.raises(ValueError, match="slice length must be"):
         _ = mdim_ims[tuple(range(mdim_ims.ndim + 1))]
 
 
-def testMultiDimGets():
-    ims_np, ims = get_stacks()
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimGets(get_stacks, get_mdim_images):
+    ims_np, ims = get_stacks
+    mdim_ims_np, mdim_ims = get_mdim_images
     for X in mdim_ims:
         assert np.allclose(ims_np, X)
 
@@ -220,9 +232,9 @@ def testMultiDimGets():
     assert np.allclose(mdim_ims[:, 1:], ims[1:])
 
 
-def testMultiDimSets():
-    ims_np, ims = get_stacks()
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimSets(get_stacks, get_mdim_images):
+    ims_np, ims = get_stacks
+    mdim_ims_np, mdim_ims = get_mdim_images
     mdim_ims[0, 1] = 123
     # Check the values changed
     assert np.allclose(mdim_ims[0, 1], 123)
@@ -232,9 +244,9 @@ def testMultiDimSets():
     assert np.allclose(mdim_ims[1, :], ims_np)
 
 
-def testMultiDimSetsSlice():
-    ims_np, ims = get_stacks()
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimSetsSlice(get_stacks, get_mdim_images):
+    ims_np, ims = get_stacks
+    mdim_ims_np, mdim_ims = get_mdim_images
     # Test setting a slice
     mdim_ims[0, 1:] = 456
     # Check the values changed
@@ -244,9 +256,9 @@ def testMultiDimSetsSlice():
     assert np.allclose(mdim_ims[1, :], ims_np)
 
 
-def testMultiDimReshape():
+def testMultiDimReshape(get_mdim_images):
     # Try mdim reshape
-    mdim_ims_np, mdim_ims = get_mdim_images()
+    mdim_ims_np, mdim_ims = get_mdim_images
     X = mdim_ims.stack_reshape(*mdim_ims.stack_shape[::-1])
     assert X.stack_shape == mdim_ims.stack_shape[::-1]
     # Compare with direct np.reshape of axes of ndarray
@@ -254,22 +266,22 @@ def testMultiDimReshape():
     assert np.allclose(X.asnumpy(), mdim_ims_np.reshape(shape))
 
 
-def testMultiDimFlattens():
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimFlattens(get_mdim_images):
+    mdim_ims_np, mdim_ims = get_mdim_images
     # Try flattening
     X = mdim_ims.stack_reshape(mdim_ims.n_images)
     assert X.stack_shape, (mdim_ims.n_images,)
 
 
-def testMultiDimFlattensTrick():
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimFlattensTrick(get_mdim_images):
+    mdim_ims_np, mdim_ims = get_mdim_images
     # Try flattening with -1
     X = mdim_ims.stack_reshape(-1)
     assert X.stack_shape == (mdim_ims.n_images,)
 
 
-def testMultiDimReshapeTuples():
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimReshapeTuples(get_mdim_images):
+    mdim_ims_np, mdim_ims = get_mdim_images
     # Try flattening with (-1,)
     X = mdim_ims.stack_reshape((-1,))
     assert X.stack_shape, (mdim_ims.n_images,)
@@ -279,8 +291,8 @@ def testMultiDimReshapeTuples():
     assert X.stack_shape == mdim_ims.stack_shape[::-1]
 
 
-def testMultiDimBadReshape():
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimBadReshape(get_mdim_images):
+    mdim_ims_np, mdim_ims = get_mdim_images
     # Incorrect flat shape
     with pytest.raises(ValueError, match="Number of images"):
         _ = mdim_ims.stack_reshape(8675309)
@@ -290,11 +302,11 @@ def testMultiDimBadReshape():
         _ = mdim_ims.stack_reshape(42, 8675309)
 
 
-def testMultiDimBroadcast():
-    ims_np, ims = get_stacks()
-    mdim_ims_np, mdim_ims = get_mdim_images()
+def testMultiDimBroadcast(get_stacks, get_mdim_images):
+    ims_np, ims = get_stacks
+    mdim_ims_np, mdim_ims = get_mdim_images
     X = mdim_ims + ims
-    assert np.allclose(X[0], 2 * ims.asnumpy())
+    np.testing.assert_allclose(X[0], 2 * ims.asnumpy())
 
 
 @matplotlib_dry_run
@@ -306,12 +318,12 @@ def testShow():
     im.show()
 
 
-def test_backproject_symmetry_group():
+def test_backproject_symmetry_group(dtype):
     """
     Test backproject SymmetryGroup pass through and error message.
     """
     ary = np.random.random((5, 8, 8))
-    im = Image(ary)
+    im = Image(ary, dtype=dtype)
     rots = Rotation.generate_random_rotations(5).matrices
 
     # Attempt backproject with bad symmetry group.
@@ -324,9 +336,7 @@ def test_backproject_symmetry_group():
     assert isinstance(vol.symmetry_group, CnSymmetryGroup)
 
     # Symmetry from instance.
-    vol = im.backproject(
-        rots, symmetry_group=CnSymmetryGroup(order=3, dtype=np.float32)
-    )
+    vol = im.backproject(rots, symmetry_group=CnSymmetryGroup(order=3, dtype=dtype))
     assert isinstance(vol.symmetry_group, CnSymmetryGroup)
 
 
@@ -381,7 +391,7 @@ def test_load_bad_ext():
         _ = Image.load("bad.ext")
 
 
-def test_load_mrc():
+def test_load_mrc(dtype):
     """
     Test `Image.load` round-trip.
     """
@@ -390,27 +400,19 @@ def test_load_mrc():
     filepath = os.path.join(DATA_DIR, "sample.mrc")
 
     # Load data from file
-    im = Image.load(filepath)
-    im_64 = Image.load(filepath, dtype=np.float64)
+    im = Image.load(filepath, dtype=dtype)
 
     with tempfile.TemporaryDirectory() as tmpdir_name:
         # tmp filename
         test_filepath = os.path.join(tmpdir_name, "test.mrc")
-        test_filepath_64 = os.path.join(tmpdir_name, "test_64.mrc")
 
         im.save(test_filepath)
-        im_64.save(test_filepath_64)
 
-        im2 = Image.load(test_filepath)
-        im2_64 = Image.load(test_filepath_64, dtype=np.float64)
+        im2 = Image.load(test_filepath, dtype)
 
     # Check the single precision round-trip
     assert np.array_equal(im, im2)
-    assert im2.dtype == np.float32
-
-    # check the double precision round-trip
-    assert np.array_equal(im_64, im2_64)
-    assert im2_64.dtype == np.float64
+    assert im2.dtype == dtype
 
 
 def test_load_tiff():

From 6a75a5fa802921adb1f4c7b7f691f7dd15b08bab Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 22 Jul 2024 09:46:43 -0400
Subject: [PATCH 08/24] add image save-load test

---
 tests/test_image.py | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/tests/test_image.py b/tests/test_image.py
index 6a9888479a..887e726c0d 100644
--- a/tests/test_image.py
+++ b/tests/test_image.py
@@ -438,3 +438,30 @@ def test_load_tiff():
 
     # Check contents
     assert np.array_equal(im, im2)
+
+
+def test_save_load_pixel_size(get_images, dtype):
+    """
+    Test saving and loading an MRC with pixel size attribute
+    """
+
+    im_np, im = get_images
+
+    with tempfile.TemporaryDirectory() as tmpdir_name:
+        # tmp filename
+        test_filepath = os.path.join(tmpdir_name, "test.mrc")
+
+        # Save image to file
+        im.save(test_filepath)
+
+        # Load image from file
+        im2 = Image.load(test_filepath, dtype)
+
+    # Check we've loaded the image data
+    np.testing.assert_allclose(im2, im)
+    # Check we've loaded the image dtype
+    assert im2.dtype == im.dtype, "Image dtype mismatched on save-load"
+    # Check we've loaded the pixel size
+    np.testing.assert_almost_equal(
+        im2.pixel_size, im.pixel_size, err_msg="Image pixel_size incorrect save-load"
+    )

From 98d72d9a23271b20c956c9c5ab32f0decc87d01d Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 22 Jul 2024 10:06:19 -0400
Subject: [PATCH 09/24] add Image.downsample pixel_size test

---
 tests/test_downsample.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/tests/test_downsample.py b/tests/test_downsample.py
index 305f6e4ec4..276927d81c 100644
--- a/tests/test_downsample.py
+++ b/tests/test_downsample.py
@@ -90,6 +90,9 @@ def test_downsample_2d_case(L, L_ds):
     assert (N, L_ds, L_ds) == imgs_ds.shape
     # check center points for all images
     assert checkCenterPoint(imgs_org, imgs_ds)
+    # Confirm default `pixel_size`
+    assert imgs_org.pixel_size is None
+    assert imgs_ds.pixel_size is None
 
 
 @pytest.mark.parametrize("L", [65, 66])
@@ -103,6 +106,9 @@ def test_downsample_3d_case(L, L_ds):
     assert checkCenterPoint(vols_org, vols_ds)
     # check signal energy is conserved
     assert checkSignalEnergy(vols_org, vols_ds)
+    # Confirm default `pixel_size`
+    assert vols_org.pixel_size is None
+    assert vols_ds.pixel_size is None
 
 
 def test_integer_offsets():
@@ -155,3 +161,25 @@ def test_downsample_project(volume, res_ds):
     if volume.dtype == np.float64:
         tol = 1e-09
     np.testing.assert_allclose(im_ds_proj, im_proj_ds, atol=tol)
+
+def test_pixel_size():
+    """
+    Test downsampling is rescaling the `pixel_size` attribute.
+    """
+    # Image sizes in pixels
+    L = 8  # original
+    dsL = 5  # downsampled
+
+    # Construct a small test Image
+    img = Image(np.random.random((1, L, L)).astype(DTYPE, copy=False), pixel_size=1.23)
+
+    # Downsample the image
+    result = img.downsample(dsL)
+
+    # Confirm the pixel size is scaled
+    np.testing.assert_approx_equal(
+        result.pixel_size,
+        img.pixel_size * L / dsL,
+        err_msg="Incorrect pixel size.",
+    )
+

From da4befecd8c777d49e7ce3ec29323f4434f58326 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 22 Jul 2024 16:13:02 -0400
Subject: [PATCH 10/24] Minimally add pixel_size to sources

---
 src/aspire/source/coordinates.py |  4 +++-
 src/aspire/source/image.py       | 23 +++++++++++++++++++----
 src/aspire/source/micrograph.py  | 31 ++++++++++++++++++++++++-------
 src/aspire/source/relion.py      |  6 ++++--
 src/aspire/source/simulation.py  |  2 +-
 5 files changed, 51 insertions(+), 15 deletions(-)

diff --git a/src/aspire/source/coordinates.py b/src/aspire/source/coordinates.py
index dca7aaf873..299422df70 100644
--- a/src/aspire/source/coordinates.py
+++ b/src/aspire/source/coordinates.py
@@ -490,7 +490,9 @@ def _images(self, indices):
                     cropped = self._crop_micrograph(arr, next(coord))
                     im[i] = cropped
         # Finally, apply transforms to resulting Image
-        return self.generation_pipeline.forward(Image(im), indices)
+        return self.generation_pipeline.forward(
+            Image(im, pixel_size=self.pixel_size), indices
+        )
 
     @staticmethod
     def _is_number(text):
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index fa5be1f7f7..1c243a6f2e 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -150,7 +150,14 @@ class ImageSource(ABC):
     _mutable = True
 
     def __init__(
-        self, L, n, dtype="double", metadata=None, memory=None, symmetry_group=None
+        self,
+        L,
+        n,
+        dtype="double",
+        metadata=None,
+        memory=None,
+        symmetry_group=None,
+        pixel_size=None,
     ):
         """
         A cryo-EM ImageSource object that supplies images along with other parameters for image manipulation.
@@ -163,6 +170,7 @@ def __init__(
             The path of the base directory to use as a data store or None. If None is given, no caching is performed.
         :param symmetry_group: A SymmetryGroup instance or string indicating the underlying symmetry of the molecule.
             Defaults to the `IdentitySymmetryGroup`, which represents an asymmetric particle, if none provided.
+        :param pixel_size: Pixel size of the images in Angstroms, default `None`.
         """
 
         # Instantiate the accessor for the `images` property
@@ -172,6 +180,9 @@ def __init__(
         self._n = None
         self.n = n
         self.dtype = np.dtype(dtype)
+        if pixel_size is not None:
+            pixel_size = float(pixel_size)
+        self.pixel_size = pixel_size
 
         # The private attribute '_cached_im' can be populated by calling this object's cache() method explicitly
         self._cached_im = None
@@ -736,7 +747,7 @@ def _apply_filters(
                 f"_apply_filters() passed {type(im_orig)} instead of Image instance"
             )
             # for now just convert it
-            im_orig = Image(im_orig)
+            im_orig = Image(im_orig, pixel_size=self.pixel_size)
 
         im = im_orig.copy()
 
@@ -1481,6 +1492,7 @@ def __init__(self, src, indices, memory=None):
             dtype=src.dtype,
             metadata=metadata,
             memory=memory,
+            pixel_size=src.pixel_size,
         )
 
         # Create filter indices, these are required to pass unharmed through filter eval code
@@ -1650,7 +1662,9 @@ class ArrayImageSource(ImageSource):
     if available, is consulted directly by the parent class, bypassing `_images`.
     """
 
-    def __init__(self, im, metadata=None, angles=None, symmetry_group=None):
+    def __init__(
+        self, im, metadata=None, angles=None, symmetry_group=None, pixel_size=None
+    ):
         """
         Initialize from an `Image` object.
 
@@ -1664,7 +1678,7 @@ def __init__(self, im, metadata=None, angles=None, symmetry_group=None):
         if not isinstance(im, Image):
             logger.info("Attempting to create an Image object from Numpy array.")
             try:
-                im = Image(im)
+                im = Image(im, pixel_size=pixel_size)
             except Exception as e:
                 raise RuntimeError(
                     "Creating Image object from Numpy array failed."
@@ -1678,6 +1692,7 @@ def __init__(self, im, metadata=None, angles=None, symmetry_group=None):
             metadata=metadata,
             memory=None,
             symmetry_group=symmetry_group,
+            pixel_size=im.pixel_size,
         )
 
         self._cached_im = im
diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py
index 182133d982..dd4d9e497c 100644
--- a/src/aspire/source/micrograph.py
+++ b/src/aspire/source/micrograph.py
@@ -17,11 +17,14 @@
 
 
 class MicrographSource(ABC):
-    def __init__(self, micrograph_count, micrograph_size, dtype):
+    def __init__(self, micrograph_count, micrograph_size, dtype, pixel_size=None):
         """ """
         self.micrograph_count = int(micrograph_count)
         self.micrograph_size = int(micrograph_size)
         self.dtype = np.dtype(dtype)
+        if pixel_size is not None:
+            pixel_size = float(pixel_size)
+        self.pixel_size = pixel_size
 
         self._images_accessor = _ImageAccessor(self._images, self.micrograph_count)
 
@@ -85,7 +88,7 @@ def show(self, *args, **kwargs):
         """
         Helper function to display micrograph. See Image.show().
         """
-        Image(self.asnumpy()).show(*args, **kwargs)
+        Image(self.asnumpy(), pixel_size=self.pixel_size).show(*args, **kwargs)
 
     @property
     def images(self):
@@ -107,7 +110,7 @@ def _images(self, indices):
 
 
 class ArrayMicrographSource(MicrographSource):
-    def __init__(self, micrographs, dtype=None):
+    def __init__(self, micrographs, dtype=None, pixel_size=None):
         """
         Instantiate a `MicrographSource` with `micrographs`.
 
@@ -140,6 +143,7 @@ def __init__(self, micrographs, dtype=None):
             micrograph_count=micrographs.shape[0],
             micrograph_size=micrographs.shape[-1],
             dtype=dtype or micrographs.dtype,
+            pixel_size=pixel_size,
         )
 
         # We're already backed by an array, access it directly.
@@ -152,11 +156,11 @@ def _images(self, indices):
         :param indices: A 1-D Numpy array of integer indices.
         :return: An array backed `MicrographSource` object representing the micrographs for `indices`.
         """
-        return Image(self._data[indices])
+        return Image(self._data[indices], pixel_size=self.pixel_size)
 
 
 class DiskMicrographSource(MicrographSource):
-    def __init__(self, micrographs_path, dtype=None):
+    def __init__(self, micrographs_path, dtype=None, pixel_size=None):
         """
         Instantiate a `MicrographSource` with `micrographs_path`.
 
@@ -190,11 +194,16 @@ def __init__(self, micrographs_path, dtype=None):
         # Load the first micrograph to infer shape/type
         # Size will be checked during on-the-fly loading of subsequent micrographs.
         micrograph0 = Image.load(self.micrograph_files[0])
+        if micrograph0.pixel_size is not None and micrograph0.pixel_size != pixel_size:
+            raise NotImplementedError(
+                f"Mismatched pixel size. {micrograph0.pixel_size} defined in {self.micrograph_files[0]}, but provided {pixel_size}."
+            )
 
         super().__init__(
             micrograph_count=len(self.micrograph_files),
             micrograph_size=micrograph0.resolution,
             dtype=dtype or micrograph0.dtype,
+            pixel_size=pixel_size,
         )
 
         # Prepare accessor to load files from disk on the fly.
@@ -262,8 +271,16 @@ def _images(self, indices):
                 )
             # Assign to array, implicitly performs casting to dtype
             micrographs[i] = micrograph.asnumpy()
+            # Assert pixel_size
+            if (
+                micrograph.pixel_size is not None
+                and micrograph.pixel_size != pixel_size
+            ):
+                raise NotImplementedError(
+                    f"Mismatched pixel size. {micrograph.pixel_size} defined in {self.micrograph_files[ind]}, but provided {pixel_size}."
+                )
 
-        return Image(micrographs)
+        return Image(micrographs, pixel_size=self.pixel_size)
 
 
 class MicrographSimulation(MicrographSource):
@@ -557,7 +574,7 @@ def _clean_images(self, indices):
             self.pad : self.micrograph_size + self.pad,
             self.pad : self.micrograph_size + self.pad,
         ]
-        return Image(clean_micrograph)
+        return Image(clean_micrograph, pixel_size=self.pixel_size)
 
     def get_micrograph_index(self, particle_index):
         """
diff --git a/src/aspire/source/relion.py b/src/aspire/source/relion.py
index 99907cbf6a..bd6d660dd3 100644
--- a/src/aspire/source/relion.py
+++ b/src/aspire/source/relion.py
@@ -59,7 +59,6 @@ def __init__(
 
         self.filepath = filepath
         self.data_folder = data_folder
-        self.pixel_size = pixel_size
         self.B = B
         self.n_workers = n_workers
         self.max_rows = max_rows
@@ -112,6 +111,7 @@ def __init__(
             metadata=metadata,
             symmetry_group=symmetry_group,
             memory=memory,
+            pixel_size=pixel_size,
         )
 
         # CTF estimation parameters coming from Relion
@@ -272,4 +272,6 @@ def load_single_mrcs(filepath, indices):
         logger.debug(f"Loading {len(indices)} images complete")
 
         # Finally, apply transforms to resulting Image
-        return self.generation_pipeline.forward(Image(im), indices)
+        return self.generation_pipeline.forward(
+            Image(im, pixel_size=self.pixel_size), indices
+        )
diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py
index 304d5be56d..c23d50df11 100644
--- a/src/aspire/source/simulation.py
+++ b/src/aspire/source/simulation.py
@@ -260,7 +260,7 @@ def _projections(self, indices):
             im_k = self.vols[k - 1].project(rot_matrices=rot)
             im[idx_k, :, :] = im_k.asnumpy()
 
-        return Image(im)
+        return Image(im, pixel_size=self.pixel_size)
 
     @property
     def clean_images(self):

From 3d8686191c8235979c617ea77741fcef912287c4 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 22 Jul 2024 16:14:59 -0400
Subject: [PATCH 11/24] tox caught incorrect var

---
 src/aspire/source/micrograph.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py
index dd4d9e497c..8a8b616f6e 100644
--- a/src/aspire/source/micrograph.py
+++ b/src/aspire/source/micrograph.py
@@ -274,10 +274,10 @@ def _images(self, indices):
             # Assert pixel_size
             if (
                 micrograph.pixel_size is not None
-                and micrograph.pixel_size != pixel_size
+                and micrograph.pixel_size != self.pixel_size
             ):
                 raise NotImplementedError(
-                    f"Mismatched pixel size. {micrograph.pixel_size} defined in {self.micrograph_files[ind]}, but provided {pixel_size}."
+                    f"Mismatched pixel size. {micrograph.pixel_size} defined in {self.micrograph_files[ind]}, but provided {self.pixel_size}."
                 )
 
         return Image(micrographs, pixel_size=self.pixel_size)

From 60d6b26e0253df98c7f684dd612072413b2740bc Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 26 Jul 2024 09:11:52 -0400
Subject: [PATCH 12/24] lint

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

diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index 1484a9c344..bb04309d33 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -553,7 +553,9 @@ def downsample(self, ds_res, mask=None, zero_nyquist=True):
 
         # returns a new Volume object
         return self.__class__(
-            xp.asnumpy(out), pixel_size=ds_pixel_size, symmetry_group=self.symmetry_group
+            xp.asnumpy(out),
+            pixel_size=ds_pixel_size,
+            symmetry_group=self.symmetry_group,
         ).stack_reshape(original_stack_shape)
 
     def shift(self):

From 8fbdc86a92e22e2a0449a8877d5fc888aa6a0a9b Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 26 Jul 2024 09:32:59 -0400
Subject: [PATCH 13/24] Change default CTFFilter pixel_size from 10 to 1

---
 gallery/tutorials/tutorials/cov3d_simulation.py | 4 +++-
 src/aspire/operators/filters.py                 | 8 ++++----
 tests/test_anisotropic_noise.py                 | 4 +++-
 tests/test_filters.py                           | 8 ++++++--
 tests/test_simulation.py                        | 8 ++++++--
 5 files changed, 22 insertions(+), 10 deletions(-)

diff --git a/gallery/tutorials/tutorials/cov3d_simulation.py b/gallery/tutorials/tutorials/cov3d_simulation.py
index 5fced70fbb..741a47de99 100644
--- a/gallery/tutorials/tutorials/cov3d_simulation.py
+++ b/gallery/tutorials/tutorials/cov3d_simulation.py
@@ -47,7 +47,9 @@
     L=img_size,
     n=num_imgs,
     vols=vols,
-    unique_filters=[RadialCTFFilter(defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)],
+    unique_filters=[
+        RadialCTFFilter(pixel_size=10, defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)
+    ],
     dtype=dtype,
 )
 
diff --git a/src/aspire/operators/filters.py b/src/aspire/operators/filters.py
index bb7491c780..e75187fb4a 100644
--- a/src/aspire/operators/filters.py
+++ b/src/aspire/operators/filters.py
@@ -403,7 +403,7 @@ def __init__(self, dim=None):
 class CTFFilter(Filter):
     def __init__(
         self,
-        pixel_size=10,
+        pixel_size=1,
         voltage=200,
         defocus_u=15000,
         defocus_v=15000,
@@ -415,7 +415,7 @@ def __init__(
         """
         A CTF (Contrast Transfer Function) Filter
 
-        :param pixel_size:  Pixel size in angstrom
+        :param pixel_size:  Pixel size in angstrom, default 1.
         :param voltage:     Electron voltage in kV
         :param defocus_u:   Defocus depth along the u-axis in angstrom
         :param defocus_v:   Defocus depth along the v-axis in angstrom
@@ -425,7 +425,7 @@ def __init__(
         :param B:           Envelope decay in inverse square angstrom (default 0)
         """
         super().__init__(dim=2, radial=defocus_u == defocus_v)
-        self.pixel_size = pixel_size
+        self.pixel_size = float(pixel_size)
         self.voltage = voltage
         self.wavelength = voltage_to_wavelength(self.voltage)
         self.defocus_u = defocus_u
@@ -482,7 +482,7 @@ def scale(self, c=1):
 
 class RadialCTFFilter(CTFFilter):
     def __init__(
-        self, pixel_size=10, voltage=200, defocus=15000, Cs=2.26, alpha=0.07, B=0
+        self, pixel_size=1, voltage=200, defocus=15000, Cs=2.26, alpha=0.07, B=0
     ):
         super().__init__(
             pixel_size=pixel_size,
diff --git a/tests/test_anisotropic_noise.py b/tests/test_anisotropic_noise.py
index 2fd1d13ca8..caaedc4aff 100644
--- a/tests/test_anisotropic_noise.py
+++ b/tests/test_anisotropic_noise.py
@@ -20,7 +20,9 @@ def setUp(self):
             n=1024,
             vols=self.vol,
             unique_filters=[
-                RadialCTFFilter(defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)
+                # Set legacy pixel size
+                RadialCTFFilter(pixel_size=10, defocus=d)
+                for d in np.linspace(1.5e4, 2.5e4, 7)
             ],
             dtype=self.dtype,
         )
diff --git a/tests/test_filters.py b/tests/test_filters.py
index 911e3b347b..b0b23bb74f 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -118,7 +118,8 @@ def testRadialCTFFilter(self):
         self.assertEqual(result.shape, (256,))
 
     def testRadialCTFFilterGrid(self):
-        filter = RadialCTFFilter(defocus=2.5e4)
+        # Set legacy pixel size
+        filter = RadialCTFFilter(pixel_size=10, defocus=2.5e4)
         result = filter.evaluate_grid(8, dtype=self.dtype)
 
         self.assertEqual(result.shape, (8, 8))
@@ -218,7 +219,10 @@ def testRadialCTFFilterGrid(self):
         )
 
     def testRadialCTFFilterMultiplierGrid(self):
-        filter = RadialCTFFilter(defocus=2.5e4) * RadialCTFFilter(defocus=2.5e4)
+        # Set legacy pixel size
+        filter = RadialCTFFilter(pixel_size=10, defocus=2.5e4) * RadialCTFFilter(
+            pixel_size=10, defocus=2.5e4
+        )
         result = filter.evaluate_grid(8, dtype=self.dtype)
 
         self.assertEqual(result.shape, (8, 8))
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index ad8a7ff4e1..69b0fef5d3 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -116,8 +116,10 @@ def setUp(self):
             n=self.n,
             L=self.L,
             vols=self.vols,
+            # Set legacy pixel_size
             unique_filters=[
-                RadialCTFFilter(defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)
+                RadialCTFFilter(pixel_size=10, defocus=d)
+                for d in np.linspace(1.5e4, 2.5e4, 7)
             ],
             noise_adder=WhiteNoiseAdder(var=1),
             dtype=self.dtype,
@@ -168,7 +170,9 @@ def testSimulationCached(self):
             vols=self.vols,
             offsets=self.sim.offsets,
             unique_filters=[
-                RadialCTFFilter(defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)
+                # Set legacy pixel size
+                RadialCTFFilter(pixel_size=10, defocus=d)
+                for d in np.linspace(1.5e4, 2.5e4, 7)
             ],
             noise_adder=WhiteNoiseAdder(var=1),
             dtype=self.dtype,

From f3620a9bed67b1e7f70092d33c76bc126a64e3b1 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 29 Jul 2024 12:55:23 -0400
Subject: [PATCH 14/24] add pixel_size synth volume classes

---
 src/aspire/source/simulation.py       | 28 ++++++++++++++++++++
 src/aspire/volume/volume.py           |  6 ++++-
 src/aspire/volume/volume_synthesis.py | 37 +++++++++++++++++++--------
 tests/test_simulation.py              | 32 ++++++++++++++++++-----
 tests/test_synthetic_volume.py        | 15 +++++++++--
 5 files changed, 97 insertions(+), 21 deletions(-)

diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py
index c23d50df11..2fc47e9afa 100644
--- a/src/aspire/source/simulation.py
+++ b/src/aspire/source/simulation.py
@@ -50,6 +50,7 @@ def __init__(
         memory=None,
         noise_adder=None,
         symmetry_group=None,
+        pixel_size=None,
     ):
         """
         A `Simulation` object that supplies images along with other parameters for image manipulation.
@@ -79,6 +80,7 @@ def __init__(
         :param noise_adder: Optionally append instance of `NoiseAdder`
             to generation pipeline.
         :param symmetry_group: A SymmetryGroup instance or string indicating symmetry of the molecule.
+        :param pixel_size: Pixel size of the images in Angstroms, default `None`.
 
         :return: A Simulation object.
         """
@@ -91,6 +93,7 @@ def __init__(
             self.vols = AsymmetricVolume(
                 L=L or 8,
                 C=C,
+                pixel_size=pixel_size,
                 seed=self.seed,
                 dtype=dtype or np.float32,
             ).generate()
@@ -122,6 +125,7 @@ def __init__(
             dtype=self.vols.dtype,
             memory=memory,
             symmetry_group=symmetry_group,
+            pixel_size=self.vols.pixel_size,
         )
 
         # If a user provides both `L` and `vols`, resolution should match.
@@ -153,6 +157,7 @@ def __init__(
         if unique_filters is None:
             unique_filters = []
         self.unique_filters = unique_filters
+        self._check_filter_pixel_size(unique_filters)
         # sim_filters must be a deep copy so that it is not changed
         # when unique_filters is changed
         self.sim_filters = copy.deepcopy(unique_filters)
@@ -231,6 +236,29 @@ def _populate_ctf_metadata(self, filter_indices):
             filter_values,
         )
 
+    def _check_filter_pixel_size(self, unique_filters):
+        """
+        Private method to ensure user provided filters match `Simulation` pixel size.
+
+        When `Simulation.pixel_size` is not `None`, any
+        `unique_filters` having a non-matching `pixel_size` attribute
+        will raise.
+        """
+
+        # Skip when Simulation pixel_size is not explicitly provided.
+        if self.pixel_size is None:
+            return
+
+        for f in unique_filters:
+            f_pixel_size = getattr(f, "pixel_size", None)
+            if f_pixel_size is not None and not np.isclose(
+                f_pixel_size, self.pixel_size
+            ):
+                raise ValueError(
+                    f"`Simulation.pixel_size` {self.pixel_size} does not match filter {f} {f_pixel_size}."
+                    "Ensure provided `pixel_size` attributes match."
+                )
+
     @property
     def projections(self):
         """
diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index bb04309d33..05e08352c3 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -260,10 +260,14 @@ def stack_reshape(self, *args):
         )
 
     def __repr__(self):
+        px_msg = "."
+        if self.pixel_size is not None:
+            px_msg = f" with pixel_size={self.pixel_size} Angstroms."
+
         msg = (
             f"{self.n_vols} {self.dtype} volumes arranged as a {self.stack_shape} stack"
         )
-        msg += f" each of size {self.resolution}x{self.resolution}x{self.resolution}."
+        msg += f" each of size {self.resolution}x{self.resolution}x{self.resolution}{px_msg}"
         return msg
 
     def __len__(self):
diff --git a/src/aspire/volume/volume_synthesis.py b/src/aspire/volume/volume_synthesis.py
index b9514df5ea..e0f1ef84f7 100644
--- a/src/aspire/volume/volume_synthesis.py
+++ b/src/aspire/volume/volume_synthesis.py
@@ -16,11 +16,12 @@
 
 
 class SyntheticVolumeBase(abc.ABC):
-    def __init__(self, L, C, seed=None, dtype=np.float64):
+    def __init__(self, L, C, pixel_size=None, seed=None, dtype=np.float64):
         self.L = L
         self.C = C
         self.seed = seed
         self.dtype = dtype
+        self.pixel_size = pixel_size
 
     @abc.abstractmethod
     def generate(self):
@@ -39,18 +40,21 @@ class GaussianBlobsVolume(SyntheticVolumeBase):
     A base class for all volumes which are generated with randomized 3D Gaussians.
     """
 
-    def __init__(self, L, C, K=16, alpha=1, seed=None, dtype=np.float64):
+    def __init__(
+        self, L, C, K=16, alpha=1, pixel_size=None, seed=None, dtype=np.float64
+    ):
         """
         :param L: Resolution of the Volume(s) in pixels.
         :param C: Number of Volumes to generate.
         :param K: Number of Gaussian blobs used to construct the Volume(s).
         :param alpha: Scaling factor for variance of Gaussian blobs. Default=1.
+        :param pixel_size: Optional voxel_size in Angstroms. Default=1.
         :param seed: Random seed for generating random Gaussian blobs.
         :param dtype: dtype for Volume(s)
         """
         self.K = int(K)
         self.alpha = float(alpha)
-        super().__init__(L=L, C=C, seed=seed, dtype=dtype)
+        super().__init__(L=L, C=C, pixel_size=pixel_size, seed=seed, dtype=dtype)
         self._set_symmetry_group()
 
     @abc.abstractproperty
@@ -75,7 +79,11 @@ def generate(self):
         """
         vol = self._gaussian_blob_vols()
         bump_mask = bump_3d(self.L, spread=5, dtype=self.dtype)
-        return Volume(bump_mask * vol, symmetry_group=self.symmetry_group)
+        return Volume(
+            bump_mask * vol,
+            symmetry_group=self.symmetry_group,
+            pixel_size=self.pixel_size,
+        )
 
     def _gaussian_blob_vols(self):
         """
@@ -168,18 +176,23 @@ class CnSymmetricVolume(GaussianBlobsVolume):
     A Volume object with cyclically symmetric volumes constructed of random 3D Gaussian blobs.
     """
 
-    def __init__(self, L, C, order, K=16, alpha=1, seed=None, dtype=np.float64):
+    def __init__(
+        self, L, C, order, K=16, alpha=1, pixel_size=None, seed=None, dtype=np.float64
+    ):
         """
         :param L: Resolution of the Volume(s) in pixels.
         :param C: Number of Volumes to generate.
         :param order: An integer representing the cyclic order of the Volume(s).
         :param K: Number of Gaussian blobs used to construct the Volume(s).
+        :param pixel_size: Optional voxel_size in Angstroms. Default=1.
         :param seed: Random seed for generating random Gaussian blobs.
         :param dtype: dtype for Volume(s)
         """
         self.order = int(order)
         self._check_order()
-        super().__init__(L=L, C=C, K=K, alpha=alpha, seed=seed, dtype=dtype)
+        super().__init__(
+            L=L, C=C, K=K, alpha=alpha, pixel_size=pixel_size, seed=seed, dtype=dtype
+        )
 
     def _check_order(self):
         if self.order < 2:
@@ -239,8 +252,10 @@ class AsymmetricVolume(CnSymmetricVolume):
     An asymmetric Volume constructed of random 3D Gaussian blobs with compact support in the unit sphere.
     """
 
-    def __init__(self, L, C, K=64, seed=None, dtype=np.float64):
-        super().__init__(L=L, C=C, K=K, order=1, seed=seed, dtype=dtype)
+    def __init__(self, L, C, K=64, pixel_size=None, seed=None, dtype=np.float64):
+        super().__init__(
+            L=L, C=C, K=K, order=1, pixel_size=pixel_size, seed=seed, dtype=dtype
+        )
 
     def _check_order(self):
         if self.order != 1:
@@ -260,8 +275,8 @@ class LegacyVolume(AsymmetricVolume):
     An asymmetric Volume object used for testing of legacy code.
     """
 
-    def __init__(self, L, C=2, K=16, seed=0, dtype=np.float64):
-        super().__init__(L=L, C=C, K=K, seed=seed, dtype=dtype)
+    def __init__(self, L, C=2, K=16, pixel_size=None, seed=0, dtype=np.float64):
+        super().__init__(L=L, C=C, K=K, pixel_size=pixel_size, seed=seed, dtype=dtype)
 
     def generate(self):
         """
@@ -272,4 +287,4 @@ def generate(self):
         # Swap axes to retain Legacy xyz-indexing.
         vols = np.swapaxes(vols, 1, 3)
 
-        return Volume(vols)
+        return Volume(vols, pixel_size=self.pixel_size)
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 69b0fef5d3..92a29e225e 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -21,15 +21,19 @@ class SingleSimTestCase(TestCase):
     """Test we can construct a length 1 Sim."""
 
     def setUp(self):
-        self.sim = Simulation(
-            n=1,
-            L=8,
-        )
+        self._pixel_size = 1.23  # Test value
+
+        self.sim = Simulation(n=1, L=8, pixel_size=self._pixel_size)
 
     def testImage(self):
         """Test we can get an Image from a length 1 Sim."""
         _ = self.sim.images[0]
 
+    def testPixelSize(self):
+        """Test pixel_size is passing through Simulation."""
+        self.assertTrue(self.sim.pixel_size == self._pixel_size)
+        self.assertTrue(self.sim.pixel_size == self.sim.vols.pixel_size)
+
     @matplotlib_dry_run
     def testImageShow(self):
         self.sim.images[:].show()
@@ -106,9 +110,12 @@ def setUp(self):
         self.n = 1024
         self.L = 8
         self.dtype = np.float32
+        # Set legacy pixel_size
+        self._pixel_size = 10
 
         self.vols = LegacyVolume(
             L=self.L,
+            pixel_size=self._pixel_size,
             dtype=self.dtype,
         ).generate()
 
@@ -116,9 +123,8 @@ def setUp(self):
             n=self.n,
             L=self.L,
             vols=self.vols,
-            # Set legacy pixel_size
             unique_filters=[
-                RadialCTFFilter(pixel_size=10, defocus=d)
+                RadialCTFFilter(pixel_size=self._pixel_size, defocus=d)
                 for d in np.linspace(1.5e4, 2.5e4, 7)
             ],
             noise_adder=WhiteNoiseAdder(var=1),
@@ -171,7 +177,7 @@ def testSimulationCached(self):
             offsets=self.sim.offsets,
             unique_filters=[
                 # Set legacy pixel size
-                RadialCTFFilter(pixel_size=10, defocus=d)
+                RadialCTFFilter(pixel_size=self._pixel_size, defocus=d)
                 for d in np.linspace(1.5e4, 2.5e4, 7)
             ],
             noise_adder=WhiteNoiseAdder(var=1),
@@ -663,3 +669,15 @@ def test_cached_image_accessors():
     np.testing.assert_allclose(cached_src.projections[:], src.projections[:])
     np.testing.assert_allclose(cached_src.images[:], src.images[:])
     np.testing.assert_allclose(cached_src.clean_images[:], src.clean_images[:])
+
+
+def test_mismatched_pixel_size():
+    """
+    Confirm raises error when explicit Simulation and CTFFilter pixel sizes mismatch.
+    """
+    # Create a CTF with a pixel_size
+    filts = [RadialCTFFilter(pixel_size=5)]
+
+    # Try to create a Simulation with a different pixel_size
+    with raises(ValueError, match=r"pixel_size.*does not match filter.*"):
+        _ = Simulation(L=8, n=1, C=1, pixel_size=10, unique_filters=filts)
diff --git a/tests/test_synthetic_volume.py b/tests/test_synthetic_volume.py
index ddcdcbcab5..01fb74ec90 100644
--- a/tests/test_synthetic_volume.py
+++ b/tests/test_synthetic_volume.py
@@ -85,6 +85,10 @@ def vol_fixture(request, dtype_fixture):
     if len(params) > 2:
         vol_kwargs["order"] = params[2]
 
+    # Assign some volumes a pixel_size, leave others as default.
+    if res % 2:
+        vol_kwargs["pixel_size"] = 3.0
+
     return vol_class(**vol_kwargs)
 
 
@@ -96,8 +100,15 @@ def test_volume_repr(vol_fixture):
 
 
 def test_volume_generate(vol_fixture):
-    """Test that a volume is generated"""
-    _ = vol_fixture.generate()
+    """
+    Test that a volume is generated
+    and stores pixel_size when provided.
+    """
+    v = vol_fixture.generate()
+
+    # In vol_fixture, we assign pixel_size to volumes having odd voxel sizes.
+    if vol_fixture.L % 2:
+        assert v.pixel_size == 3
 
 
 def test_simulation_init(vol_fixture):

From 1e9543ff7f9ce531217714b7684aa7b19614e76b Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 29 Jul 2024 14:49:45 -0400
Subject: [PATCH 15/24] self review cleanup

---
 gallery/tutorials/aspire_introduction.py | 2 +-
 gallery/tutorials/pipeline_demo.py       | 2 +-
 src/aspire/image/image.py                | 6 +++---
 src/aspire/source/image.py               | 3 ++-
 src/aspire/source/micrograph.py          | 1 +
 src/aspire/source/simulation.py          | 2 +-
 src/aspire/volume/volume.py              | 5 ++---
 src/aspire/volume/volume_synthesis.py    | 4 ++--
 8 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/gallery/tutorials/aspire_introduction.py b/gallery/tutorials/aspire_introduction.py
index 648750ac01..cffe6d544e 100644
--- a/gallery/tutorials/aspire_introduction.py
+++ b/gallery/tutorials/aspire_introduction.py
@@ -571,7 +571,7 @@ def noise_function(x, y):
 
 # Generate several CTFs.
 ctf_filters = [
-    RadialCTFFilter(pixel_size=5, defocus=d)
+    RadialCTFFilter(pixel_size=vol_ds.pixel_size, defocus=d)
     for d in np.linspace(defocus_min, defocus_max, defocus_ct)
 ]
 
diff --git a/gallery/tutorials/pipeline_demo.py b/gallery/tutorials/pipeline_demo.py
index 8910436de2..77d304b156 100644
--- a/gallery/tutorials/pipeline_demo.py
+++ b/gallery/tutorials/pipeline_demo.py
@@ -63,7 +63,7 @@
 defocus_ct = 7
 
 ctf_filters = [
-    RadialCTFFilter(pixel_size=5, defocus=d)
+    RadialCTFFilter(pixel_size=vol.pixel_size, defocus=d)
     for d in np.linspace(defocus_min, defocus_max, defocus_ct)
 ]
 
diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 984b1e6540..bce9535665 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -130,10 +130,10 @@ def load_tiff(filepath):
 
     # Future todo, extract `voxel_size` if available in TIFF tags (custom tag?)
     # For now, default to `None`.
-    voxel_size = None
+    pixel_size = None
 
     # Cast image data as numpy array
-    return np.array(img), voxel_size
+    return np.array(img), pixel_size
 
 
 class Image:
@@ -160,7 +160,7 @@ def __init__(self, data, pixel_size=None, dtype=None):
 
         :param data: Numpy array containing image data with shape
             `(..., resolution, resolution)`.
-        :param pixel_size: Optional pixel size in Angstroms.
+        :param pixel_size: Optional pixel size in angstroms.
             When provided will be saved with `mrc` metadata.
             Default of `None` will not write to file,
             but will be considered unit pixels (1) for FSC.
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 1c243a6f2e..4d256bd414 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -170,7 +170,7 @@ def __init__(
             The path of the base directory to use as a data store or None. If None is given, no caching is performed.
         :param symmetry_group: A SymmetryGroup instance or string indicating the underlying symmetry of the molecule.
             Defaults to the `IdentitySymmetryGroup`, which represents an asymmetric particle, if none provided.
-        :param pixel_size: Pixel size of the images in Angstroms, default `None`.
+        :param pixel_size: Pixel size of the images in angstroms, default `None`.
         """
 
         # Instantiate the accessor for the `images` property
@@ -1673,6 +1673,7 @@ def __init__(
         :param metadata: A Dataframe of metadata information corresponding to this ImageSource's images
         :param angles: Optional n-by-3 array of rotation angles corresponding to `im`.
         :param symmetry_group: A SymmetryGroup instance or string indicating the underlying symmetry of the molecule.
+        :param pixel_size: Pixel size of the images in angstroms, default `None`.
         """
 
         if not isinstance(im, Image):
diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py
index 8a8b616f6e..ece811eb30 100644
--- a/src/aspire/source/micrograph.py
+++ b/src/aspire/source/micrograph.py
@@ -122,6 +122,7 @@ def __init__(self, micrographs, dtype=None, pixel_size=None):
             Currently only `float32` and `float64` are supported.
             Note, due to limitations of common MRC implementations,
             saving is limited to single precision.
+        :param pixel_size: Pixel size of the images in angstroms, default `None`.
         """
 
         # Check micrographs is an array
diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py
index 2fc47e9afa..7d80b570d2 100644
--- a/src/aspire/source/simulation.py
+++ b/src/aspire/source/simulation.py
@@ -80,7 +80,7 @@ def __init__(
         :param noise_adder: Optionally append instance of `NoiseAdder`
             to generation pipeline.
         :param symmetry_group: A SymmetryGroup instance or string indicating symmetry of the molecule.
-        :param pixel_size: Pixel size of the images in Angstroms, default `None`.
+        :param pixel_size: Pixel size of the images in angstroms, default `None`.
 
         :return: A Simulation object.
         """
diff --git a/src/aspire/volume/volume.py b/src/aspire/volume/volume.py
index 05e08352c3..5e2212e958 100644
--- a/src/aspire/volume/volume.py
+++ b/src/aspire/volume/volume.py
@@ -76,7 +76,7 @@ def __init__(self, data, dtype=None, pixel_size=None, symmetry_group=None):
             `(..., resolution, resolution, resolution)`.
         :param dtype: Optionally cast `data` to this dtype.
             Defaults to `data.dtype`.
-        :param pixel_size: Optional voxel_size in Angstroms.
+        :param pixel_size: Optional voxel_size in angstroms.
             When provided will be saved with `map`/`mrc` metadata.
             Default of `None` will not write to file,
             but will be considered unit pixels (1) for FSC.
@@ -262,7 +262,7 @@ def stack_reshape(self, *args):
     def __repr__(self):
         px_msg = "."
         if self.pixel_size is not None:
-            px_msg = f" with pixel_size={self.pixel_size} Angstroms."
+            px_msg = f" with pixel_size={self.pixel_size} angstroms."
 
         msg = (
             f"{self.n_vols} {self.dtype} volumes arranged as a {self.stack_shape} stack"
@@ -423,7 +423,6 @@ def project(self, rot_matrices, zero_nyquist=True):
             im_f[:, :, 0] = 0
 
         im_f = fft.centered_ifft2(im_f)
-        # todo add pixel_size to Image
         return aspire.image.Image(xp.asnumpy(im_f.real), pixel_size=self.pixel_size)
 
     def to_vec(self):
diff --git a/src/aspire/volume/volume_synthesis.py b/src/aspire/volume/volume_synthesis.py
index e0f1ef84f7..59947b79ea 100644
--- a/src/aspire/volume/volume_synthesis.py
+++ b/src/aspire/volume/volume_synthesis.py
@@ -48,7 +48,7 @@ def __init__(
         :param C: Number of Volumes to generate.
         :param K: Number of Gaussian blobs used to construct the Volume(s).
         :param alpha: Scaling factor for variance of Gaussian blobs. Default=1.
-        :param pixel_size: Optional voxel_size in Angstroms. Default=1.
+        :param pixel_size: Optional voxel_size in angstroms. Default=1.
         :param seed: Random seed for generating random Gaussian blobs.
         :param dtype: dtype for Volume(s)
         """
@@ -184,7 +184,7 @@ def __init__(
         :param C: Number of Volumes to generate.
         :param order: An integer representing the cyclic order of the Volume(s).
         :param K: Number of Gaussian blobs used to construct the Volume(s).
-        :param pixel_size: Optional voxel_size in Angstroms. Default=1.
+        :param pixel_size: Optional voxel_size in angstroms. Default=1.
         :param seed: Random seed for generating random Gaussian blobs.
         :param dtype: dtype for Volume(s)
         """

From dc60bedaa29de403c5f5c8dba9d2f4e2359656aa Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 1 Aug 2024 11:41:44 -0400
Subject: [PATCH 16/24] add units to error messages

---
 src/aspire/source/micrograph.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py
index ece811eb30..521c416b4c 100644
--- a/src/aspire/source/micrograph.py
+++ b/src/aspire/source/micrograph.py
@@ -197,7 +197,7 @@ def __init__(self, micrographs_path, dtype=None, pixel_size=None):
         micrograph0 = Image.load(self.micrograph_files[0])
         if micrograph0.pixel_size is not None and micrograph0.pixel_size != pixel_size:
             raise NotImplementedError(
-                f"Mismatched pixel size. {micrograph0.pixel_size} defined in {self.micrograph_files[0]}, but provided {pixel_size}."
+                f"Mismatched pixel size. {micrograph0.pixel_size} angstroms defined in {self.micrograph_files[0]}, but provided {pixel_size} angstroms."
             )
 
         super().__init__(
@@ -278,7 +278,7 @@ def _images(self, indices):
                 and micrograph.pixel_size != self.pixel_size
             ):
                 raise NotImplementedError(
-                    f"Mismatched pixel size. {micrograph.pixel_size} defined in {self.micrograph_files[ind]}, but provided {self.pixel_size}."
+                    f"Mismatched pixel size. {micrograph.pixel_size} angstroms defined in {self.micrograph_files[ind]}, but provided {self.pixel_size} angstroms."
                 )
 
         return Image(micrographs, pixel_size=self.pixel_size)

From 112f4ac7ee80409743ad239aa468cb0457ac3365 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 1 Aug 2024 11:46:28 -0400
Subject: [PATCH 17/24] add units to error messages

---
 src/aspire/source/simulation.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/aspire/source/simulation.py b/src/aspire/source/simulation.py
index 7d80b570d2..331b86e442 100644
--- a/src/aspire/source/simulation.py
+++ b/src/aspire/source/simulation.py
@@ -255,7 +255,7 @@ def _check_filter_pixel_size(self, unique_filters):
                 f_pixel_size, self.pixel_size
             ):
                 raise ValueError(
-                    f"`Simulation.pixel_size` {self.pixel_size} does not match filter {f} {f_pixel_size}."
+                    f"`Simulation.pixel_size` {self.pixel_size} does not match filter {f} pixel size {f_pixel_size}."
                     "Ensure provided `pixel_size` attributes match."
                 )
 

From 794b57add8c38706a46ece9e795829fdefa40717 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 1 Aug 2024 11:49:02 -0400
Subject: [PATCH 18/24] add pixel_size to image repr

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

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index bce9535665..1bbe2ba160 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -361,8 +361,12 @@ def flip(self, axis=-2):
         return self.__class__(np.flip(self._data, axis), pixel_size=self.pixel_size)
 
     def __repr__(self):
+        px_msg = "."
+        if self.pixel_size is not None:
+            px_msg = f" with pixel_size={self.pixel_size} angstroms."
+
         msg = f"{self.n_images} {self.dtype} images arranged as a {self.stack_shape} stack"
-        msg += f" each of size {self.resolution}x{self.resolution}."
+        msg += f" each of size {self.resolution}x{self.resolution}{px_msg}"
         return msg
 
     def asnumpy(self):

From 525b4077c80b4b3cb554d0bbbe6a83b034ed45d8 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 1 Aug 2024 11:51:32 -0400
Subject: [PATCH 19/24] correct synth vol doc strings

---
 src/aspire/volume/volume_synthesis.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/aspire/volume/volume_synthesis.py b/src/aspire/volume/volume_synthesis.py
index 59947b79ea..43f794bfaf 100644
--- a/src/aspire/volume/volume_synthesis.py
+++ b/src/aspire/volume/volume_synthesis.py
@@ -48,7 +48,10 @@ def __init__(
         :param C: Number of Volumes to generate.
         :param K: Number of Gaussian blobs used to construct the Volume(s).
         :param alpha: Scaling factor for variance of Gaussian blobs. Default=1.
-        :param pixel_size: Optional voxel_size in angstroms. Default=1.
+        :param pixel_size: Optional voxel_size in angstroms.
+            When provided will be saved with `map`/`mrc` metadata.
+            Default of `None` will not write to file,
+            but will be considered unit pixels (1) for FSC.
         :param seed: Random seed for generating random Gaussian blobs.
         :param dtype: dtype for Volume(s)
         """
@@ -184,7 +187,10 @@ def __init__(
         :param C: Number of Volumes to generate.
         :param order: An integer representing the cyclic order of the Volume(s).
         :param K: Number of Gaussian blobs used to construct the Volume(s).
-        :param pixel_size: Optional voxel_size in angstroms. Default=1.
+        :param pixel_size: Optional voxel_size in angstroms.
+            When provided will be saved with `map`/`mrc` metadata.
+            Default of `None` will not write to file,
+            but will be considered unit pixels (1) for FSC.
         :param seed: Random seed for generating random Gaussian blobs.
         :param dtype: dtype for Volume(s)
         """

From 96e271f83dde87073b8dff155863a174cd35f257 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 1 Aug 2024 11:57:59 -0400
Subject: [PATCH 20/24] add direct volume pixel_size attr test

---
 tests/test_volume.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/tests/test_volume.py b/tests/test_volume.py
index bec31d6ed0..9059a89a67 100644
--- a/tests/test_volume.py
+++ b/tests/test_volume.py
@@ -30,6 +30,7 @@ def res_id(params):
 
 
 RES = [42, 43]
+TEST_PX_SZ = 4.56
 
 
 @pytest.fixture(params=RES, ids=res_id, scope="module")
@@ -75,7 +76,7 @@ def vols_1(data_1):
 
 @pytest.fixture
 def vols_2(data_2):
-    return Volume(data_2, pixel_size=4.56)
+    return Volume(data_2, pixel_size=TEST_PX_SZ)
 
 
 @pytest.fixture
@@ -295,6 +296,13 @@ def test_save_load(vols_1):
         assert vols_loaded_double.pixel_size is None, "Pixel size should be None"
 
 
+def test_volume_pixel_size(vols_2):
+    """
+    Test volume is storing pixel_size attribute.
+    """
+    assert np.isclose(TEST_PX_SZ, vols_2.pixel_size), "Incorrect Volume pixel_size"
+
+
 def test_save_load_pixel_size(vols_2):
     # Create a tmpdir in a context. It will be cleaned up on exit.
     with tempfile.TemporaryDirectory() as tmpdir:

From 5f01347801fb575fd581b5c7aff1c273ab10ab31 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 26 Aug 2024 08:45:58 -0400
Subject: [PATCH 21/24] Use ValueError for px sz mismatch

---
 src/aspire/source/micrograph.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/aspire/source/micrograph.py b/src/aspire/source/micrograph.py
index 521c416b4c..2d654401b5 100644
--- a/src/aspire/source/micrograph.py
+++ b/src/aspire/source/micrograph.py
@@ -196,7 +196,7 @@ def __init__(self, micrographs_path, dtype=None, pixel_size=None):
         # Size will be checked during on-the-fly loading of subsequent micrographs.
         micrograph0 = Image.load(self.micrograph_files[0])
         if micrograph0.pixel_size is not None and micrograph0.pixel_size != pixel_size:
-            raise NotImplementedError(
+            raise ValueError(
                 f"Mismatched pixel size. {micrograph0.pixel_size} angstroms defined in {self.micrograph_files[0]}, but provided {pixel_size} angstroms."
             )
 
@@ -277,7 +277,7 @@ def _images(self, indices):
                 micrograph.pixel_size is not None
                 and micrograph.pixel_size != self.pixel_size
             ):
-                raise NotImplementedError(
+                raise ValueError(
                     f"Mismatched pixel size. {micrograph.pixel_size} angstroms defined in {self.micrograph_files[ind]}, but provided {self.pixel_size} angstroms."
                 )
 

From bfd4476812e6fae8550f95aca3420754dc9745c7 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 26 Aug 2024 08:59:40 -0400
Subject: [PATCH 22/24] Test px sz assignment with assert_approx_equal

---
 tests/test_synthetic_volume.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/tests/test_synthetic_volume.py b/tests/test_synthetic_volume.py
index 01fb74ec90..fec7591764 100644
--- a/tests/test_synthetic_volume.py
+++ b/tests/test_synthetic_volume.py
@@ -20,6 +20,9 @@
 # dtype fixture to pass into volume fixture.
 DTYPES = [np.float32, pytest.param(np.float64, marks=pytest.mark.expensive)]
 
+# Pixel sized used to test assignment
+PXSZ = 3.0
+
 
 @pytest.fixture(params=DTYPES)
 def dtype_fixture(request):
@@ -87,7 +90,7 @@ def vol_fixture(request, dtype_fixture):
 
     # Assign some volumes a pixel_size, leave others as default.
     if res % 2:
-        vol_kwargs["pixel_size"] = 3.0
+        vol_kwargs["pixel_size"] = PXSZ
 
     return vol_class(**vol_kwargs)
 
@@ -108,7 +111,7 @@ def test_volume_generate(vol_fixture):
 
     # In vol_fixture, we assign pixel_size to volumes having odd voxel sizes.
     if vol_fixture.L % 2:
-        assert v.pixel_size == 3
+        np.testing.assert_approx_equal(v.pixel_size, PXSZ)
 
 
 def test_simulation_init(vol_fixture):

From 76ac0b804fd766add8ce2b24161dd7ce90391be2 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 26 Aug 2024 14:22:32 -0400
Subject: [PATCH 23/24] merge conflict lint

---
 tests/test_downsample.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_downsample.py b/tests/test_downsample.py
index 276927d81c..4c990212ce 100644
--- a/tests/test_downsample.py
+++ b/tests/test_downsample.py
@@ -162,6 +162,7 @@ def test_downsample_project(volume, res_ds):
         tol = 1e-09
     np.testing.assert_allclose(im_ds_proj, im_proj_ds, atol=tol)
 
+
 def test_pixel_size():
     """
     Test downsampling is rescaling the `pixel_size` attribute.
@@ -182,4 +183,3 @@ def test_pixel_size():
         img.pixel_size * L / dsL,
         err_msg="Incorrect pixel size.",
     )
-

From 393b089b80a816c86f24a67db0c257842a43837f Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 26 Aug 2024 14:35:09 -0400
Subject: [PATCH 24/24] merge conflict test

---
 tests/test_volume.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_volume.py b/tests/test_volume.py
index 9059a89a67..ac86c4096b 100644
--- a/tests/test_volume.py
+++ b/tests/test_volume.py
@@ -589,7 +589,7 @@ def test_downsample(res):
     # Confirm the pixel size is scaled
     np.testing.assert_approx_equal(
         result.pixel_size,
-        vols.pixel_size * res / ds_res,
+        vols.pixel_size * og_res / ds_res,
         err_msg="Incorrect pixel size.",
     )