From 648c30c137f4fad923bfa7e85b1658915d1a945a Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 4 Aug 2025 07:54:03 -0400
Subject: [PATCH 01/23] Add crop preprocess method

---
 src/aspire/image/xform.py  | 22 ++++++++++++++++++++++
 src/aspire/source/image.py | 22 ++++++++++++++++++++++
 2 files changed, 44 insertions(+)

diff --git a/src/aspire/image/xform.py b/src/aspire/image/xform.py
index 10fddf1df2..7448f12c6e 100644
--- a/src/aspire/image/xform.py
+++ b/src/aspire/image/xform.py
@@ -228,6 +228,28 @@ def __str__(self):
         return f"Downsample (Resolution {self.resolution})"
 
 
+class Crop(Xform):
+    """
+    A Xform that crops an Image object to a size specified by this Xform's size.
+    """
+
+    def __init__(self, L):
+        """
+        Initialize Xform to crop Image to a specific size.
+
+        :param L: int - new size, should be <= the current size
+            of this Image
+        """
+        self.L = L
+        super().__init__()
+
+    def _forward(self, im, indices):
+        return im[..., : self.L, : self.L]
+
+    def __str__(self):
+        return f"Crop (Size {self.L})"
+
+
 class LegacyWhiten(Xform):
     """
     A Xform that implements MATLAB legacy whitening.
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index f4f30bf481..2b540eb415 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -12,6 +12,7 @@
 from aspire.abinitio import CLOrient3D, CLSync3N
 from aspire.image import Image, normalize_bg
 from aspire.image.xform import (
+    Crop,
     Downsample,
     FilterXform,
     IndexedXform,
@@ -788,6 +789,27 @@ def downsample(self, L, zero_nyquist=True, legacy=False):
 
         self.L = L
 
+    @_as_copy
+    def crop(self, L):
+        """
+        Crop images down to size L.
+
+        Used for reproducing legacy MATLAB workflows.
+        For other applications, `downsample` is preferred.
+
+        Note, cropping makes no adjustments for centering/offsets etc.
+        """
+
+        if L > self.L:
+            raise ValueError(
+                "Max desired resolution {L} should be less than the current resolution {self.L}."
+            )
+        logger.info(f"Cropping shape of source images = {L,L}")
+
+        self.generation_pipeline.add_xform(Crop(L=L))
+
+        self.L = L
+
     @_as_copy
     def whiten(self, noise_estimate=None, epsilon=None):
         """

From 629ac2d08bfaab83a6f3c92fae5c3895e3da78de Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 4 Aug 2025 08:12:42 -0400
Subject: [PATCH 02/23] Create legacy_normalize_background with workflow
 default radius

---
 src/aspire/image/image.py  | 27 ++++++++++-----------------
 src/aspire/source/image.py | 32 +++++++++++++++++++++++++-------
 2 files changed, 35 insertions(+), 24 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index ed061ecec7..64b2bd84dd 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -24,19 +24,21 @@
 logger = logging.getLogger(__name__)
 
 
-def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, legacy=False):
+def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, shifted=False, ddof=0):
     """
-    Normalize backgrounds and apply to a stack of images
+    Normalize backgrounds and apply to a stack of images.
+
+    To recreate legacy MATLAB workflow results, review parameters used in
+    `ImageSource.legacy_normalize_background`.
 
     :param imgs: A stack of images in N-by-L-by-L array
     :param bg_radius: Radius cutoff to be considered as background (in image size)
     :param do_ramp: When it is `True`, fit a ramping background to the data
-            and subtract. Namely perform normalization based on values from each image.
-            Otherwise, a constant background level from all images is used.
-    :param legacy: Option to match Matlab legacy normalize_background. Default, False,
-        uses ASPIRE-Python implementation. When True, ramping is disabled, a shifted
-        2d grid and alternative `bg_radius` is used to generate the background mask,
-        and standard deviation is computed using N - 1 degrees of freedom.
+        and subtract. Namely perform normalization based on values from each image.
+        Otherwise, a constant background level from all images is used.
+    :param shifted: Optionally shifts 2d grid by 1/2 pixel for even
+        resolution to replicate MATLAB.
+    :param ddof: Degrees of freedom for standard deviation.
     :return: The modified images
     """
     if imgs.ndim > 3:
@@ -45,15 +47,6 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, legacy=False):
         )
     L = imgs.shape[-1]
 
-    # Make adjustments for legacy mode
-    shifted = False
-    ddof = 0  # Degrees of freedom for standard deviation
-    if legacy:
-        do_ramp = False
-        shifted = True  # Shifts 2d grid by 1/2 pixel for even resolution
-        bg_radius = 2 * (L // 2) / L
-        ddof = 1
-
     # Generate background mask
     input_dtype = imgs.dtype
     grid = grid_2d(L, shifted=shifted, indexing="yx", dtype=input_dtype)
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 2b540eb415..4fea725b00 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -966,8 +966,25 @@ def invert_contrast(self, batch_size=512):
         logger.info("Adding Scaling Xform to end of generation pipeline")
         self.generation_pipeline.add_xform(Multiply(scale_factor))
 
+    def legacy_normalize_background(self):
+        """
+        Match MATLAB's `normalize_background` workflow method.
+
+        Ramping is disabled.
+        A shifted 2d grid and alternative `bg_radius` is used to generate the background mask.
+        Standard deviation is computed using N - 1 degrees of freedom.
+        """
+
+        # Radius definition is here:
+        # https://github.com/PrincetonUniversity/aspire/blob/760a43b35453e55ff2d9354339e9ffa109a25371/workflow/cryo_workflow_preprocess_execute.m#L166
+        bg_radius = 2 * np.floor(self.L * 0.45) / self.L
+
+        return self.normalize_background(
+            bg_radius=bg_radius, do_ramp=False, shifted=True, ddof=1
+        )
+
     @_as_copy
-    def normalize_background(self, bg_radius=1.0, do_ramp=True, legacy=False):
+    def normalize_background(self, bg_radius=1.0, do_ramp=True, shifted=False, ddof=0):
         """
         Normalize the images by the noise background
 
@@ -980,11 +997,8 @@ def normalize_background(self, bg_radius=1.0, do_ramp=True, legacy=False):
         :param do_ramp: When it is `True`, fit a ramping background to the data
             and subtract. Namely perform normalization based on values from each image.
             Otherwise, a constant background level from all images is used.
-        :param legacy: Option to match Matlab legacy normalize_background. Default, False,
-            uses ASPIRE-Python implementation. When True, ramping is disable, a shifted
-            2d grid and alternative `bg_radius` is used to generate the background mask,
-            and standard deviation is computed using N - 1 degrees of freedom.
-        :return: On return, the `ImageSource` object has been modified in place.
+
+        :return: Returns`ImageSource` object.
         """
 
         logger.info(
@@ -993,7 +1007,11 @@ def normalize_background(self, bg_radius=1.0, do_ramp=True, legacy=False):
         )
         self.generation_pipeline.add_xform(
             LambdaXform(
-                normalize_bg, bg_radius=bg_radius, do_ramp=do_ramp, legacy=legacy
+                normalize_bg,
+                bg_radius=bg_radius,
+                do_ramp=do_ramp,
+                shifted=shifted,
+                ddof=ddof,
             )
         )
 

From 65a5fe4f5a85e51454e995ad4f629b19030fb131 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 4 Aug 2025 08:49:49 -0400
Subject: [PATCH 03/23] Create legacy_downsample and lightly refactor
 downsample/Downsample methods

---
 src/aspire/image/image.py  | 21 +++++++++++----------
 src/aspire/source/image.py | 32 +++++++++++++++++++++++++++++---
 2 files changed, 40 insertions(+), 13 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 64b2bd84dd..a72738c06a 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -513,7 +513,7 @@ def legacy_whiten(self, psd, delta):
 
         return Image(res)
 
-    def downsample(self, ds_res, zero_nyquist=True, legacy=False):
+    def downsample(self, ds_res, zero_nyquist=True, centered_fft=True):
         """
         Downsample Image to a specific resolution. This method returns a new Image.
 
@@ -521,8 +521,9 @@ def downsample(self, ds_res, zero_nyquist=True, legacy=False):
             of this Image
         :param zero_nyquist: Option to keep or remove Nyquist frequency for even
             resolution (boolean). Defaults to zero_nyquist=True, removing the Nyquist frequency.
-        :param legacy: Option to match legacy Matlab downsample method (boolean).
-            Default of False uses `centered_fft` to maintain ASPIRE-Python centering conventions.
+        :param centered_fft: Default of True uses `centered_fft` to
+            maintain ASPIRE-Python centering conventions.
+
         :return: The downsampled Image object.
         """
 
@@ -533,25 +534,25 @@ def downsample(self, ds_res, zero_nyquist=True, legacy=False):
         # because all of the subsequent calls until `asnumpy` are GPU
         # when xp and fft in `cupy` mode.
 
-        if legacy:
-            fx = fft.fftshift(fft.fft2(xp.asarray(im._data)))
-        else:
+        if centered_fft:
             # compute FT with centered 0-frequency
             fx = fft.centered_fft2(xp.asarray(im._data))
+        else:
+            fx = fft.fftshift(fft.fft2(xp.asarray(im._data)))
 
         # crop 2D Fourier transform for each image
         crop_fx = crop_pad_2d(fx, ds_res)
 
         # If downsampled resolution is even, optionally zero out the nyquist frequency.
-        if ds_res % 2 == 0 and zero_nyquist and not legacy:
+        if ds_res % 2 == 0 and zero_nyquist:
             crop_fx[:, 0, :] = 0
             crop_fx[:, :, 0] = 0
 
         # take back to real space, discard complex part, and scale
-        if legacy:
-            out = fft.ifft2(fft.ifftshift(crop_fx))
-        else:
+        if centered_fft:
             out = fft.centered_ifft2(crop_fx)
+        else:
+            out = fft.ifft2(fft.ifftshift(crop_fx))
 
         # The parenths are required because dtype casting semantics
         # differs between Numpy 1, 2, and CuPy.
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 4fea725b00..26b20b34a2 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -769,8 +769,32 @@ def _images(self, indices):
         Subclasses handle cached image check as well as applying transforms in the generation pipeline.
         """
 
+    def legacy_downsample(self, L):
+        """
+        Reproduce MATLAB's `downsample` workflow method.
+
+        Downsample Image to `L-by-L` pixels.
+
+        :param L: int - new resolution, should be <= the current resolution
+            of this Image
+        :return: The downsampled `ImageSource` object.
+        """
+        return self.downsample(L=L, zero_nyquist=False, centered_fft=False)
+
     @_as_copy
-    def downsample(self, L, zero_nyquist=True, legacy=False):
+    def downsample(self, L, zero_nyquist=True, centered_fft=True):
+        """
+        Downsample Image to `L-by-L` pixels.
+
+        :param L: int - new resolution, should be <= the current resolution
+            of this Image
+        :param zero_nyquist: Option to keep or remove Nyquist frequency for even
+            resolution (boolean). Defaults to zero_nyquist=True, removing the Nyquist frequency.
+        :param centered_fft: Default of True uses `centered_fft` to
+            maintain ASPIRE-Python centering conventions.
+        :return: The downsampled `ImageSource` object.
+        """
+
         if L > self.L:
             raise ValueError(
                 "Max desired resolution {L} should be less than the current resolution {self.L}."
@@ -778,7 +802,9 @@ def downsample(self, L, zero_nyquist=True, legacy=False):
         logger.info(f"Setting max. resolution of source = {L}")
 
         self.generation_pipeline.add_xform(
-            Downsample(resolution=L, zero_nyquist=zero_nyquist, legacy=legacy)
+            Downsample(
+                resolution=L, zero_nyquist=zero_nyquist, centered_fft=centered_fft
+            )
         )
 
         ds_factor = self.L / L
@@ -968,7 +994,7 @@ def invert_contrast(self, batch_size=512):
 
     def legacy_normalize_background(self):
         """
-        Match MATLAB's `normalize_background` workflow method.
+        Reproduce MATLAB's `normalize_background` workflow method.
 
         Ramping is disabled.
         A shifted 2d grid and alternative `bg_radius` is used to generate the background mask.

From f3117e29f55c584c607e0ac55f3002ea5c64e4a3 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 4 Aug 2025 09:01:42 -0400
Subject: [PATCH 04/23] Adjust/cleanup Xforms

---
 src/aspire/image/xform.py         | 18 ++++++++++--------
 tests/test_preprocess_pipeline.py |  2 +-
 2 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/src/aspire/image/xform.py b/src/aspire/image/xform.py
index 7448f12c6e..7bfdccf55f 100644
--- a/src/aspire/image/xform.py
+++ b/src/aspire/image/xform.py
@@ -199,7 +199,7 @@ class Downsample(LinearXform):
     A Xform that downsamples an Image object to a resolution specified by this Xform's resolution.
     """
 
-    def __init__(self, resolution, zero_nyquist=True, legacy=False):
+    def __init__(self, resolution, zero_nyquist=True, centered_fft=True):
         """
         Initialize Xform to downsample Image to a specific resolution.
 
@@ -207,17 +207,19 @@ def __init__(self, resolution, zero_nyquist=True, legacy=False):
             of this Image
         :param zero_nyquist: Option to keep or remove Nyquist frequency for even
             resolution (boolean). Defaults to zero_nyquist=True, removing the Nyquist frequency.
-        :param legacy: Option to match legacy Matlab downsample method (boolean).
-            Default of False uses `centered_fft` to maintain ASPIRE-Python centering conventions.
+        :param centered_fft: Default of True uses `centered_fft` to
+            maintain ASPIRE-Python centering conventions.
         """
         self.resolution = resolution
         self.zero_nyquist = zero_nyquist
-        self.legacy = legacy
+        self.centered_fft = centered_fft
         super().__init__()
 
     def _forward(self, im, indices):
         return im.downsample(
-            self.resolution, zero_nyquist=self.zero_nyquist, legacy=self.legacy
+            self.resolution,
+            zero_nyquist=self.zero_nyquist,
+            centered_fft=self.centered_fft,
         )
 
     def _adjoint(self, im, indices):
@@ -225,7 +227,7 @@ def _adjoint(self, im, indices):
         raise NotImplementedError("Adjoint of downsampling not implemented yet.")
 
     def __str__(self):
-        return f"Downsample (Resolution {self.resolution})"
+        return f"Downsample (resolution={self.resolution}, zero_nyquist={self.zero_nyquist}, centered_fft={self.centered_fft}) Xform"
 
 
 class Crop(Xform):
@@ -247,7 +249,7 @@ def _forward(self, im, indices):
         return im[..., : self.L, : self.L]
 
     def __str__(self):
-        return f"Crop (Size {self.L})"
+        return f"Crop({self.L}) Xform"
 
 
 class LegacyWhiten(Xform):
@@ -278,7 +280,7 @@ def _forward(self, im, indices):
         return im.legacy_whiten(self.psd, self.delta)
 
     def __str__(self):
-        return "Legacy Whitening Xform."
+        return "LegacyWhiten() Xform"
 
 
 class FilterXform(SymmetricXform):
diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index 94c5c1dae4..928cfb4c7f 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -101,7 +101,7 @@ def test_norm_background_legacy(L, dtype):
     bg_radius = 2 * (L // 2) / L
     grid = grid_2d(sim.L, shifted=True, indexing="yx", dtype=dtype)
     mask = grid["r"] > bg_radius
-    sim = sim.normalize_background(legacy=True)
+    sim = sim.legacy_normalize_background()
     imgs_nb = sim.images[:].asnumpy()
     new_mean = np.mean(imgs_nb[:, mask])
     new_variance = np.var(imgs_nb[:, mask], ddof=1)

From 19970ae6db47a40a2ce966c2f1e2cfe4ab6567b7 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 4 Aug 2025 09:37:08 -0400
Subject: [PATCH 05/23] update unit tests

---
 tests/test_downsample.py          | 14 ++++++++++++--
 tests/test_preprocess_pipeline.py | 19 ++++++++++++++-----
 2 files changed, 26 insertions(+), 7 deletions(-)

diff --git a/tests/test_downsample.py b/tests/test_downsample.py
index 525dba1de6..2879fc3dd7 100644
--- a/tests/test_downsample.py
+++ b/tests/test_downsample.py
@@ -126,6 +126,12 @@ def test_integer_offsets():
 RES_DS = [32, 33]
 LEGACY = [True, False]
 
+# These parameters are defined in the top level `ImageSource.legacy_downsample` wrapper.
+im_ds_legacy_flags = {
+    "zero_nyquist": False,
+    "centered_fft": False,
+}
+
 
 @pytest.fixture(params=DTYPES, ids=lambda x: f"dtype={x}", scope="module")
 def dtype(request):
@@ -165,7 +171,10 @@ def test_downsample_project(volume, res_ds, legacy):
     """
     rot = np.eye(3, dtype=volume.dtype)  # project along z-axis
     im_ds_proj = volume.downsample(res_ds, legacy=legacy).project(rot)
-    im_proj_ds = volume.project(rot).downsample(res_ds, legacy=legacy)
+    if legacy:
+        im_proj_ds = volume.project(rot).downsample(res_ds, **im_ds_legacy_flags)
+    else:
+        im_proj_ds = volume.project(rot).downsample(res_ds)
 
     tol = 1e-09
     if volume.dtype == np.float32:
@@ -199,7 +208,8 @@ def test_downsample_legacy(volume, res_ds):
     ims = src.images[:]
 
     # Legacy downsampled images.
-    ims_ds_legacy = ims.downsample(res_ds, legacy=True)
+    #   Params are defined in `ImageSource.legacy_downsample`
+    ims_ds_legacy = ims.downsample(res_ds, **im_ds_legacy_flags)
 
     # ASPIRE-Python downsample with centering adjustments for odd resolution images.
     shifts = 0.5 * np.ones((n_img, 2), dtype=dtype)
diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index 928cfb4c7f..f34f844619 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -93,15 +93,24 @@ def testNormBackground(L, dtype):
 
 
 @pytest.mark.parametrize("L, dtype", params)
-def test_norm_background_legacy(L, dtype):
-    # Legacy normalize_background uses a shifted grid, a different
+def test_norm_background_legacy_outofcore(L, dtype):
+    """
+    This executes normalize_background with the parameters found to reproduce MATLAB's "outofcore" method.
+    """
+    # Legacy "outofcore" normalize_background defaults to a shifted grid, a different
     # mask radius, disabled ramping, and N - 1 degrees of freedom
     # when computing standard deviation.
+    norm_bg_outofcore_flags = {
+        "bg_radius": 2 * (L // 2) / L,
+        "do_ramp": False,
+        "shifted": True,
+        "ddof": 1,
+    }
+
     sim = get_sim_object(L, dtype)
-    bg_radius = 2 * (L // 2) / L
     grid = grid_2d(sim.L, shifted=True, indexing="yx", dtype=dtype)
-    mask = grid["r"] > bg_radius
-    sim = sim.legacy_normalize_background()
+    mask = grid["r"] > norm_bg_outofcore_flags["bg_radius"]
+    sim = sim.normalize_background(**norm_bg_outofcore_flags)
     imgs_nb = sim.images[:].asnumpy()
     new_mean = np.mean(imgs_nb[:, mask])
     new_variance = np.var(imgs_nb[:, mask], ddof=1)

From 69be0804d764a85cd98fb7ab9397ec9be8becec2 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 4 Aug 2025 09:40:09 -0400
Subject: [PATCH 06/23] update gallery jsb examples

---
 .../experiments/experimental_abinitio_pipeline_10028_jsb.py   | 4 ++--
 .../experiments/experimental_abinitio_pipeline_10073_jsb.py   | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index da306d546c..70219ee646 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -90,10 +90,10 @@
 
 # Downsample the images.
 logger.info(f"Set the resolution to {img_size} X {img_size}")
-src = src.downsample(img_size).cache()
+src = src.legacy_downsample(img_size).cache()
 
 # Normalize the background of the images.
-src = src.normalize_background().cache()
+src = src.legacy_normalize_background().cache()
 
 # Estimate the noise and whiten based on the estimated noise.
 src = src.legacy_whiten().cache()
diff --git a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
index c1e9d673b0..232870a511 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -84,10 +84,10 @@
 
 # Downsample the images.
 logger.info(f"Set the resolution to {img_size} X {img_size}")
-src = src.downsample(img_size).cache()
+src = src.legacy_downsample(img_size).cache()
 
 # Normalize the background of the images.
-src = src.normalize_background().cache()
+src = src.legacy_normalize_background().cache()
 
 # Estimate the noise and whiten based on the estimated noise.
 src = src.legacy_whiten().cache()

From e577713ea146f38623a39661cfeeb41d23731498 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 11 Aug 2025 07:45:47 -0400
Subject: [PATCH 07/23] Add legacy crop to JSB examples

---
 .../experiments/experimental_abinitio_pipeline_10028_jsb.py    | 3 +++
 .../experiments/experimental_abinitio_pipeline_10073_jsb.py    | 3 +++
 2 files changed, 6 insertions(+)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index 70219ee646..32744d1875 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -88,6 +88,9 @@
 logger.info("Perform phase flip to input images.")
 src = src.phase_flip().cache()
 
+# Legacy MATLAB right cropped the images to an odd resolution.
+src = src.crop(src.L - 1).cache()
+
 # Downsample the images.
 logger.info(f"Set the resolution to {img_size} X {img_size}")
 src = src.legacy_downsample(img_size).cache()
diff --git a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
index 232870a511..a516524c06 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -82,6 +82,9 @@
 logger.info("Perform phase flip to input images.")
 src = src.phase_flip().cache()
 
+# Legacy MATLAB right cropped the images to an odd resolution.
+src = src.crop(src.L - 1).cache()
+
 # Downsample the images.
 logger.info(f"Set the resolution to {img_size} X {img_size}")
 src = src.legacy_downsample(img_size).cache()

From fcf52ecf7fea02ed5de5d5eb5a9f22037928fae0 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 11 Aug 2025 07:58:33 -0400
Subject: [PATCH 08/23] Use DiracBasis3D for jsb examples

---
 .../experiments/experimental_abinitio_pipeline_10028_jsb.py  | 5 ++++-
 .../experiments/experimental_abinitio_pipeline_10073_jsb.py  | 5 ++++-
 src/aspire/source/image.py                                   | 2 +-
 3 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index 32744d1875..c601e75fce 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -29,6 +29,7 @@
 import logging
 from pathlib import Path
 
+from aspire.basis import DiracBasis3D
 from aspire.denoising import LegacyClassAvgSource
 from aspire.reconstruction import MeanEstimator
 from aspire.source import OrientedSource, RelionSource
@@ -156,7 +157,9 @@
 logger.info("Begin Volume reconstruction")
 
 # Set up an estimator to perform the backprojection.
-estimator = MeanEstimator(oriented_src)
+# Legacy MATLAB FIRM used Dirac basis.
+basis3d = DiracBasis3D(oriented_src.L)
+estimator = MeanEstimator(oriented_src, basis=basis3d)
 
 # Perform the estimation and save the volume.
 estimated_volume = estimator.estimate()
diff --git a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
index a516524c06..56753752f8 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -32,6 +32,7 @@
 import numpy as np
 
 from aspire.abinitio import CLSync3N
+from aspire.basis import DiracBasis3D
 from aspire.denoising import LegacyClassAvgSource
 from aspire.reconstruction import MeanEstimator
 from aspire.source import ArrayImageSource, OrientedSource, RelionSource
@@ -160,7 +161,9 @@
 logger.info("Begin Volume reconstruction")
 
 # Set up an estimator to perform the backprojection.
-estimator = MeanEstimator(oriented_src)
+# Legacy MATLAB FIRM used Dirac basis.
+basis3d = DiracBasis3D(oriented_src.L)
+estimator = MeanEstimator(oriented_src, basis=basis3d)
 
 # Perform the estimation and save the volume.
 estimated_volume = estimator.estimate()
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 26b20b34a2..bb88f9ff1d 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -830,7 +830,7 @@ def crop(self, L):
             raise ValueError(
                 "Max desired resolution {L} should be less than the current resolution {self.L}."
             )
-        logger.info(f"Cropping shape of source images = {L,L}")
+        logger.info(f"Cropping shape of source images = {L, L}")
 
         self.generation_pipeline.add_xform(Crop(L=L))
 

From e5a0f5683e82715a6fa88d7ec8689748f7e245db Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Tue, 12 Aug 2025 07:16:11 -0400
Subject: [PATCH 09/23] forgot basis dtypes

---
 gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py | 2 +-
 gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index c601e75fce..9787074856 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -158,7 +158,7 @@
 
 # Set up an estimator to perform the backprojection.
 # Legacy MATLAB FIRM used Dirac basis.
-basis3d = DiracBasis3D(oriented_src.L)
+basis3d = DiracBasis3D(oriented_src.L, dtype=oriented_src.dtype)
 estimator = MeanEstimator(oriented_src, basis=basis3d)
 
 # Perform the estimation and save the volume.
diff --git a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
index 56753752f8..070290c2d0 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -162,7 +162,7 @@
 
 # Set up an estimator to perform the backprojection.
 # Legacy MATLAB FIRM used Dirac basis.
-basis3d = DiracBasis3D(oriented_src.L)
+basis3d = DiracBasis3D(oriented_src.L, dtype=oriented_src.dtype)
 estimator = MeanEstimator(oriented_src, basis=basis3d)
 
 # Perform the estimation and save the volume.

From da27022f417cffc93fed09d3641783a7504062d5 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Wed, 20 Aug 2025 14:21:40 -0400
Subject: [PATCH 10/23] use our existing crop_pad_2d

---
 src/aspire/image/xform.py  | 12 +++++++-----
 src/aspire/source/image.py |  6 +++---
 2 files changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/aspire/image/xform.py b/src/aspire/image/xform.py
index 7bfdccf55f..57740e18c9 100644
--- a/src/aspire/image/xform.py
+++ b/src/aspire/image/xform.py
@@ -5,6 +5,7 @@
 from joblib import Memory
 
 from aspire.image import Image
+from aspire.utils import crop_pad_2d
 
 logger = logging.getLogger(__name__)
 
@@ -230,26 +231,27 @@ def __str__(self):
         return f"Downsample (resolution={self.resolution}, zero_nyquist={self.zero_nyquist}, centered_fft={self.centered_fft}) Xform"
 
 
-class Crop(Xform):
+class CropPad(Xform):
     """
     A Xform that crops an Image object to a size specified by this Xform's size.
     """
 
-    def __init__(self, L):
+    def __init__(self, L, fill_value=0):
         """
         Initialize Xform to crop Image to a specific size.
 
         :param L: int - new size, should be <= the current size
-            of this Image
+            of this Image.
         """
         self.L = L
+        self.fill_value = fill_value
         super().__init__()
 
     def _forward(self, im, indices):
-        return im[..., : self.L, : self.L]
+        return crop_pad_2d(im, self.L, self.fill_value)
 
     def __str__(self):
-        return f"Crop({self.L}) Xform"
+        return f"CropPad({self.L}, {self.fill_value}) Xform"
 
 
 class LegacyWhiten(Xform):
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index bb88f9ff1d..bfa9c96673 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -12,7 +12,7 @@
 from aspire.abinitio import CLOrient3D, CLSync3N
 from aspire.image import Image, normalize_bg
 from aspire.image.xform import (
-    Crop,
+    CropPad,
     Downsample,
     FilterXform,
     IndexedXform,
@@ -826,13 +826,13 @@ def crop(self, L):
         Note, cropping makes no adjustments for centering/offsets etc.
         """
 
-        if L > self.L:
+        if L >= self.L:
             raise ValueError(
                 "Max desired resolution {L} should be less than the current resolution {self.L}."
             )
         logger.info(f"Cropping shape of source images = {L, L}")
 
-        self.generation_pipeline.add_xform(Crop(L=L))
+        self.generation_pipeline.add_xform(CropPad(L=L))
 
         self.L = L
 

From 7ea823490ec867c2a93b3d565e03324763b0bddc Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 21 Aug 2025 07:32:24 -0400
Subject: [PATCH 11/23] Normalize bg should be computed in double precision

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

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index a72738c06a..f75a5d0e6e 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -46,10 +46,11 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, shifted=False, ddof=0):
             "`normalize_bg` is currently limited to 1D image stacks."
         )
     L = imgs.shape[-1]
+    input_dtype = imgs.dtype
 
     # Generate background mask
-    input_dtype = imgs.dtype
-    grid = grid_2d(L, shifted=shifted, indexing="yx", dtype=input_dtype)
+    grid_dtype = np.float64
+    grid = grid_2d(L, shifted=shifted, indexing="yx", dtype=grid_dtype)
     mask = grid["r"] > bg_radius
 
     if do_ramp:
@@ -59,14 +60,14 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, shifted=False, ddof=0):
             (
                 grid["x"][mask].flatten(),
                 grid["y"][mask].flatten(),
-                np.ones(grid["y"][mask].flatten().size, dtype=input_dtype),
+                np.ones(grid["y"][mask].flatten().size, dtype=grid_dtype),
             )
         ).T
         ramp_all = np.vstack(
             (
                 grid["x"].flatten(),
                 grid["y"].flatten(),
-                np.ones(L * L, dtype=input_dtype),
+                np.ones(L * L, dtype=grid_dtype),
             )
         ).T
         mask_reshape = mask.reshape((L * L))
@@ -78,10 +79,14 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, shifted=False, ddof=0):
         imgs = imgs.reshape((-1, L, L))
 
     # Apply mask images and calculate mean and std values of background
-    mean = np.mean(imgs[:, mask], axis=1)
-    std = np.std(imgs[:, mask], ddof=ddof, axis=1)
-
-    return (imgs - mean[:, None, None]) / std[:, None, None]
+    # These should be computed and normalized as doubles
+    bg_pixels = imgs[:, mask].astype(np.float64, copy=False)
+    mean = np.mean(bg_pixels, axis=1)
+    std = np.std(bg_pixels, ddof=ddof, axis=1)
+    imgs = (imgs - mean[:, None, None]) / std[:, None, None]
+
+    # Restore input dtype
+    return imgs.astype(input_dtype, copy=False)
 
 
 def load_mrc(filepath):

From 663c1413303811c82f60f7337782d0b27311799a Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 21 Aug 2025 07:37:12 -0400
Subject: [PATCH 12/23] Whitening should also be in double precision

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

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index f75a5d0e6e..4b57a2c4df 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -488,10 +488,11 @@ def legacy_whiten(self, psd, delta):
         else:
             slc = slice(k - L_half - 1, k + L_half - 1)
 
+        # Note these computations should be in double precision
         for i, proj in enumerate(self.asnumpy()):
 
             # Zero pad the image to twice the size
-            padded_proj[slc, slc] = xp.asarray(proj)
+            padded_proj[slc, slc] = xp.asarray(proj, dtype=np.float64)
 
             # Take the Fourier Transform of the padded image.
             fpadded_proj = fft.centered_fft2(padded_proj)
@@ -513,8 +514,8 @@ def legacy_whiten(self, psd, delta):
 
             filtered_proj = filtered_proj[slc, slc].real
 
-            # Assign the resulting image.
-            res[i] = xp.asnumpy(filtered_proj)
+            # Assign the resulting image, cast if required.
+            res[i] = xp.asnumpy(filtered_proj).astype(res.dtype, copy=False)
 
         return Image(res)
 

From cdb906630a9ceccecb680da0643a2aafbe93d318 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Thu, 21 Aug 2025 10:44:03 -0400
Subject: [PATCH 13/23] Might as well extend crop to crop_pad

---
 ...experimental_abinitio_pipeline_10028_jsb.py |  2 +-
 ...experimental_abinitio_pipeline_10073_jsb.py |  2 +-
 src/aspire/source/image.py                     | 18 +++++++++++-------
 3 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index 9787074856..f156b6107d 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -90,7 +90,7 @@
 src = src.phase_flip().cache()
 
 # Legacy MATLAB right cropped the images to an odd resolution.
-src = src.crop(src.L - 1).cache()
+src = src.crop_pad(src.L - 1).cache()
 
 # Downsample the images.
 logger.info(f"Set the resolution to {img_size} X {img_size}")
diff --git a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
index 070290c2d0..bcf606c56f 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -84,7 +84,7 @@
 src = src.phase_flip().cache()
 
 # Legacy MATLAB right cropped the images to an odd resolution.
-src = src.crop(src.L - 1).cache()
+src = src.crop_pad(src.L - 1).cache()
 
 # Downsample the images.
 logger.info(f"Set the resolution to {img_size} X {img_size}")
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index bfa9c96673..57018e6a16 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -816,9 +816,9 @@ def downsample(self, L, zero_nyquist=True, centered_fft=True):
         self.L = L
 
     @_as_copy
-    def crop(self, L):
+    def crop_pad(self, L, fill_value=0):
         """
-        Crop images down to size L.
+        Crop or pad images to size L.
 
         Used for reproducing legacy MATLAB workflows.
         For other applications, `downsample` is preferred.
@@ -826,13 +826,17 @@ def crop(self, L):
         Note, cropping makes no adjustments for centering/offsets etc.
         """
 
-        if L >= self.L:
-            raise ValueError(
-                "Max desired resolution {L} should be less than the current resolution {self.L}."
+        if L < self.L:
+            logger.info(f"Cropping shape of source images = {L, L}")
+        elif L > self.L:
+            logger.info(
+                f"Padding shape of source images = {L, L} with fill_value={fill_value}"
             )
-        logger.info(f"Cropping shape of source images = {L, L}")
+        else:
+            logger.info(f"Shape of source images already {L, L}, skipping.")
+            return
 
-        self.generation_pipeline.add_xform(CropPad(L=L))
+        self.generation_pipeline.add_xform(CropPad(L=L, fill_value=fill_value))
 
         self.L = L
 

From 77ef944497561e628c2ff7125d20f54c95c5007f Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 22 Aug 2025 08:35:14 -0400
Subject: [PATCH 14/23] revert/postpose DiracBasis

---
 .../experiments/experimental_abinitio_pipeline_10028_jsb.py  | 5 +----
 .../experiments/experimental_abinitio_pipeline_10073_jsb.py  | 5 +----
 2 files changed, 2 insertions(+), 8 deletions(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index f156b6107d..7718ddaff3 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -29,7 +29,6 @@
 import logging
 from pathlib import Path
 
-from aspire.basis import DiracBasis3D
 from aspire.denoising import LegacyClassAvgSource
 from aspire.reconstruction import MeanEstimator
 from aspire.source import OrientedSource, RelionSource
@@ -157,9 +156,7 @@
 logger.info("Begin Volume reconstruction")
 
 # Set up an estimator to perform the backprojection.
-# Legacy MATLAB FIRM used Dirac basis.
-basis3d = DiracBasis3D(oriented_src.L, dtype=oriented_src.dtype)
-estimator = MeanEstimator(oriented_src, basis=basis3d)
+estimator = MeanEstimator(oriented_src)
 
 # Perform the estimation and save the volume.
 estimated_volume = estimator.estimate()
diff --git a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
index bcf606c56f..935f92d033 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -32,7 +32,6 @@
 import numpy as np
 
 from aspire.abinitio import CLSync3N
-from aspire.basis import DiracBasis3D
 from aspire.denoising import LegacyClassAvgSource
 from aspire.reconstruction import MeanEstimator
 from aspire.source import ArrayImageSource, OrientedSource, RelionSource
@@ -161,9 +160,7 @@
 logger.info("Begin Volume reconstruction")
 
 # Set up an estimator to perform the backprojection.
-# Legacy MATLAB FIRM used Dirac basis.
-basis3d = DiracBasis3D(oriented_src.L, dtype=oriented_src.dtype)
-estimator = MeanEstimator(oriented_src, basis=basis3d)
+estimator = MeanEstimator(oriented_src)
 
 # Perform the estimation and save the volume.
 estimated_volume = estimator.estimate()

From d017fb355cdeb7d928a2ecd01472c68a91be202a Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 22 Aug 2025 10:19:57 -0400
Subject: [PATCH 15/23] cleanup docstrings

---
 src/aspire/source/image.py | 37 ++++++++++++++++++++++++++-----------
 1 file changed, 26 insertions(+), 11 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 57018e6a16..733a0dfce3 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -771,9 +771,9 @@ def _images(self, indices):
 
     def legacy_downsample(self, L):
         """
-        Reproduce MATLAB's `downsample` workflow method.
+        Reproduce MATLAB's workflow downsampling.
 
-        Downsample Image to `L-by-L` pixels.
+        For uses other than MALTAB reproduction, prefer `downsample`.
 
         :param L: int - new resolution, should be <= the current resolution
             of this Image
@@ -821,9 +821,20 @@ def crop_pad(self, L, fill_value=0):
         Crop or pad images to size L.
 
         Used for reproducing legacy MATLAB workflows.
-        For other applications, `downsample` is preferred.
+        For most applications, `downsample` is preferred.
+
+        Cropping and padding makes no adjustments for centering conventions,
+        but does maintain `pixel_size`.
 
-        Note, cropping makes no adjustments for centering/offsets etc.
+        Take care regarding the cropping convention.
+        Cropping a single pixel from even down to odd left crops.
+        Cropping a single pixel from odd down to even right crops.
+        Calling this crop method for multiple pixels will crop equally from both
+        sides with any single remainder pixel following applied as above.
+
+        :param L: int - new image size in pixels.
+        :param fill_value: Value used in padding, defaults to 0.
+        :return: Cropped or padded `ImageSource`.
         """
 
         if L < self.L:
@@ -833,7 +844,9 @@ def crop_pad(self, L, fill_value=0):
                 f"Padding shape of source images = {L, L} with fill_value={fill_value}"
             )
         else:
-            logger.info(f"Shape of source images already {L, L}, skipping.")
+            logger.warning(
+                f"Shape of source images already {L, L}, skipping `CropPad`."
+            )
             return
 
         self.generation_pipeline.add_xform(CropPad(L=L, fill_value=fill_value))
@@ -998,11 +1011,13 @@ def invert_contrast(self, batch_size=512):
 
     def legacy_normalize_background(self):
         """
-        Reproduce MATLAB's `normalize_background` workflow method.
+        Reproduce MATLAB's Normalize Background workflow method.
 
         Ramping is disabled.
         A shifted 2d grid and alternative `bg_radius` is used to generate the background mask.
         Standard deviation is computed using N - 1 degrees of freedom.
+
+        :return: `ImageSource` with normalized background.
         """
 
         # Radius definition is here:
@@ -1016,19 +1031,19 @@ def legacy_normalize_background(self):
     @_as_copy
     def normalize_background(self, bg_radius=1.0, do_ramp=True, shifted=False, ddof=0):
         """
-        Normalize the images by the noise background
+        Normalize the images by the noise background.
 
         This is done by shifting the image density by the mean value of background
         and scaling the image density by the standard deviation of background.
-        From the implementation level, we modify the `ImageSource` in-place by
-        appending the `Add` and `Multiple` filters to the generation pipeline.
 
         :param bg_radius: Radius cutoff to be considered as background (in image size)
         :param do_ramp: When it is `True`, fit a ramping background to the data
             and subtract. Namely perform normalization based on values from each image.
             Otherwise, a constant background level from all images is used.
-
-        :return: Returns`ImageSource` object.
+        :param shifted: Optionally shifts 2d grid by 1/2 pixel for even
+            resolution to replicate MATLAB.
+        :param ddof: Degrees of freedom for standard deviation.
+        :return: `ImageSource` object with normalized background.
         """
 
         logger.info(

From f27f7a9e70a8e8d5128b04d15d98a9b1d32bb0fa Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 22 Aug 2025 10:20:11 -0400
Subject: [PATCH 16/23] cleanup docstrings

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

diff --git a/src/aspire/image/xform.py b/src/aspire/image/xform.py
index 57740e18c9..838b77c28a 100644
--- a/src/aspire/image/xform.py
+++ b/src/aspire/image/xform.py
@@ -233,7 +233,7 @@ def __str__(self):
 
 class CropPad(Xform):
     """
-    A Xform that crops an Image object to a size specified by this Xform's size.
+    A Xform that crops or pads an Image object to a specified size.
     """
 
     def __init__(self, L, fill_value=0):
@@ -242,6 +242,7 @@ def __init__(self, L, fill_value=0):
 
         :param L: int - new size, should be <= the current size
             of this Image.
+        :param fill_value: Optional value for padding, default 0.
         """
         self.L = L
         self.fill_value = fill_value

From c4bf0b131789ba998183902da89a41f2810386cd Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 22 Aug 2025 10:20:25 -0400
Subject: [PATCH 17/23] add crop_pad unit tests

---
 tests/test_preprocess_pipeline.py | 139 ++++++++++++++++++++++++++++++
 1 file changed, 139 insertions(+)

diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index f34f844619..95d8960097 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -278,3 +278,142 @@ def testInvertContrast(L, dtype):
     # dtype of returned images should be the same
     assert dtype == imgs1_rc.dtype
     assert dtype == imgs2_rc.dtype
+
+
+@pytest.mark.parametrize("L, dtype", params)
+def test_crop(L, dtype):
+    """
+    Test cropping and document convention via code.
+    """
+
+    sim1 = get_sim_object(L, dtype)
+    ref_images = sim1.images[:].asnumpy()
+
+    if L % 2:  # L odd
+        # Cropping odd by one should remove last row and last col.
+        crop_odd_to_even_one = sim1.crop_pad(L - 1)
+        np.testing.assert_allclose(
+            crop_odd_to_even_one.images[:], ref_images[..., :-1, :-1]
+        )
+
+        # Cropping by two should remove first+last row and first+last col.
+        crop_odd_to_odd_two = sim1.crop_pad(L - 2)
+        np.testing.assert_allclose(
+            crop_odd_to_odd_two.images[:], ref_images[..., 1:-1, 1:-1]
+        )
+
+        # Cropping by many even should remove equal first+last rows and first+last cols.
+        many_even = 32
+        k = many_even // 2
+        crop_odd_to_odd_many_even = sim1.crop_pad(L - many_even)
+        np.testing.assert_allclose(
+            crop_odd_to_odd_many_even.images[:], ref_images[..., k:-k, k:-k]
+        )
+
+        # Cropping by many odd should remove remove equal first+(last+1) rows and first+(last+1) cols.
+        many_odd = 33
+        k = many_odd // 2
+        crop_odd_to_even_many_odd = sim1.crop_pad(L - many_odd)
+        np.testing.assert_allclose(
+            crop_odd_to_even_many_odd.images[:], ref_images[..., k : -k - 1, k : -k - 1]
+        )
+    else:  # L even
+        # Cropping even by one should remove first row and first col.
+        crop_even_to_odd_one = sim1.crop_pad(L - 1)
+        np.testing.assert_allclose(
+            crop_even_to_odd_one.images[:], ref_images[..., 1:, 1:]
+        )
+
+        # Cropping by two should remove first+last row and first+last col.
+        crop_even_to_even_two = sim1.crop_pad(L - 2)
+        np.testing.assert_allclose(
+            crop_even_to_even_two.images[:], ref_images[..., 1:-1, 1:-1]
+        )
+
+        # Cropping by many even should remove equal first+last rows and first+last cols.
+        many_even = 32
+        k = many_even // 2
+        crop_even_to_even_many_even = sim1.crop_pad(L - many_even)
+        np.testing.assert_allclose(
+            crop_even_to_even_many_even.images[:], ref_images[..., k:-k, k:-k]
+        )
+
+        # Cropping by many odd should remove remove equal (first+1)+last) rows and (first+1)+last cols.
+        many_odd = 33
+        k = many_odd // 2
+        crop_even_to_even_many_odd = sim1.crop_pad(L - many_odd)
+        np.testing.assert_allclose(
+            crop_even_to_even_many_odd.images[:],
+            ref_images[..., k + 1 : -k, k + 1 : -k],
+        )
+
+
+@pytest.mark.parametrize("L, dtype", params)
+def test_pad(L, dtype):
+    """
+    Test pad and document convention via code.
+    """
+
+    sim1 = get_sim_object(L, dtype)
+    ref_images = sim1.images[:].asnumpy()
+
+    if L % 2:  # L odd
+        # Padding odd by one should zero pad first row and first col.
+        pad_odd_to_even_one = sim1.crop_pad(L + 1)
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (1, 0), (1, 0)))
+        np.testing.assert_allclose(pad_odd_to_even_one.images[:], ref)
+
+        # Padding odd by two should zero pad first row and first col.
+        pad_odd_to_odd_two = sim1.crop_pad(L + 2)
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (1, 1), (1, 1)))
+        np.testing.assert_allclose(pad_odd_to_odd_two.images[:], ref)
+
+        # Padding odd to even by many should zero pad the first+1 and last cols equally.
+        many_odd = 33
+        pad_odd_to_even_many = sim1.crop_pad(L + many_odd)
+        k = many_odd // 2
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (k + 1, k), (k + 1, k)))
+        np.testing.assert_allclose(pad_odd_to_even_many.images[:], ref)
+
+        # Padding odd to odd by many should pad the first and last cols equally.
+        # This test will also excecise `fill_value`
+        many_even = 32
+        fill = -1
+        pad_odd_to_odd_many = sim1.crop_pad(L + many_even, fill_value=fill)
+        k = many_even // 2
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (k, k), (k, k)), constant_values=fill)
+        np.testing.assert_allclose(pad_odd_to_odd_many.images[:], ref)
+    else:  # L even
+        # Padding even by one should zero pad last row and last col.
+        pad_even_to_odd_one = sim1.crop_pad(L + 1)
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (0, 1), (0, 1)))
+        np.testing.assert_allclose(pad_even_to_odd_one.images[:], ref)
+
+        # Padding even by two should zero pad first row and first col.
+        pad_even_to_even_two = sim1.crop_pad(L + 2)
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (1, 1), (1, 1)))
+        np.testing.assert_allclose(pad_even_to_even_two.images[:], ref)
+
+        # Padding even to even by many should zero pad the first and last cols equally.
+        many_even = 32
+        pad_even_to_even_many = sim1.crop_pad(L + many_even)
+        k = many_even // 2
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (k, k), (k, k)))
+        np.testing.assert_allclose(pad_even_to_even_many.images[:], ref)
+
+        # Padding even to odd by many should pad the first and last+1 cols equally.
+        # This test will also excecise `fill_value`
+        many_odd = 33
+        fill = -1
+        pad_even_to_odd_many = sim1.crop_pad(L + many_odd, fill_value=fill)
+        k = many_odd // 2
+        # Test image content
+        ref = np.pad(ref_images, ((0, 0), (k, k + 1), (k, k + 1)), constant_values=fill)
+        np.testing.assert_allclose(pad_even_to_odd_many.images[:], ref)

From 7caf1649432b7505dec9735a9124c4d4ee1fb347 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 22 Aug 2025 10:24:23 -0400
Subject: [PATCH 18/23] update crop comment

---
 gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py | 2 +-
 gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index 7718ddaff3..dee01c4868 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -88,7 +88,7 @@
 logger.info("Perform phase flip to input images.")
 src = src.phase_flip().cache()
 
-# Legacy MATLAB right cropped the images to an odd resolution.
+# Legacy MATLAB cropped the images to an odd resolution.
 src = src.crop_pad(src.L - 1).cache()
 
 # Downsample the images.
diff --git a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
index 935f92d033..43daaeb7bc 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -82,7 +82,7 @@
 logger.info("Perform phase flip to input images.")
 src = src.phase_flip().cache()
 
-# Legacy MATLAB right cropped the images to an odd resolution.
+# Legacy MATLAB cropped the images to an odd resolution.
 src = src.crop_pad(src.L - 1).cache()
 
 # Downsample the images.

From 36394743944f7f2ef27fc31040276caf1d5bf356 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 22 Aug 2025 10:25:44 -0400
Subject: [PATCH 19/23] words

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

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 733a0dfce3..48da7d259d 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -830,7 +830,7 @@ def crop_pad(self, L, fill_value=0):
         Cropping a single pixel from even down to odd left crops.
         Cropping a single pixel from odd down to even right crops.
         Calling this crop method for multiple pixels will crop equally from both
-        sides with any single remainder pixel following applied as above.
+        sides with any single remainder pixel applied as above.
 
         :param L: int - new image size in pixels.
         :param fill_value: Value used in padding, defaults to 0.

From 76d2c9083b4ac8f067a6d0a7d4ef05aa6e9770ed Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Fri, 22 Aug 2025 10:36:24 -0400
Subject: [PATCH 20/23] add test for legacy_normalize_background

---
 tests/test_preprocess_pipeline.py | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index 95d8960097..558d1eefcb 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -123,6 +123,37 @@ def test_norm_background_legacy_outofcore(L, dtype):
     np.testing.assert_equal(dtype, imgs_nb.dtype)
 
 
+@pytest.mark.parametrize("L, dtype", params)
+def test_legacy_normalize_background(L, dtype):
+    """
+    This executes legacy_normalize_background.
+    """
+    # Legacy normalize_background defaults to a shifted grid, a different
+    # 0.45 mask radius, disabled ramping, and N - 1 degrees of freedom
+    # when computing standard deviation.
+    norm_bg_legacy_flags = {
+        "bg_radius": 2 * np.floor(L * 0.45) / L,
+        "do_ramp": False,
+        "shifted": True,
+        "ddof": 1,
+    }
+
+    sim = get_sim_object(L, dtype)
+    grid = grid_2d(sim.L, shifted=True, indexing="yx", dtype=dtype)
+    mask = grid["r"] > norm_bg_legacy_flags["bg_radius"]
+    sim = sim.legacy_normalize_background()
+    imgs_nb = sim.images[:].asnumpy()
+    new_mean = np.mean(imgs_nb[:, mask])
+    new_variance = np.var(imgs_nb[:, mask], ddof=1)
+
+    # new mean of noise should be close to zero and variance should be close to 1
+    np.testing.assert_array_less(new_mean, 3e-4)
+    np.testing.assert_array_less(abs(new_variance - 1), 2e-3)
+
+    # dtype of returned images should be the same
+    np.testing.assert_equal(dtype, imgs_nb.dtype)
+
+
 @pytest.mark.parametrize("dtype", [np.float32, np.float64])
 def testWhiten(dtype):
     # Note this atol holds only for L even. Odd tested in testWhiten2.

From a6a74acb8ae862a4fb4497054f4812d920de1589 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 25 Aug 2025 13:42:41 -0400
Subject: [PATCH 21/23] remove CropPad docstring oversight

---
 src/aspire/image/xform.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/aspire/image/xform.py b/src/aspire/image/xform.py
index 838b77c28a..316d8ca3c8 100644
--- a/src/aspire/image/xform.py
+++ b/src/aspire/image/xform.py
@@ -240,8 +240,7 @@ def __init__(self, L, fill_value=0):
         """
         Initialize Xform to crop Image to a specific size.
 
-        :param L: int - new size, should be <= the current size
-            of this Image.
+        :param L: int - new size
         :param fill_value: Optional value for padding, default 0.
         """
         self.L = L

From 3b232966cfcd29bf1d3b88ff871f1eb6fe5a72fd Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Mon, 25 Aug 2025 13:44:06 -0400
Subject: [PATCH 22/23] fixed comment typo

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

diff --git a/tests/test_preprocess_pipeline.py b/tests/test_preprocess_pipeline.py
index 558d1eefcb..ebd6e1c593 100644
--- a/tests/test_preprocess_pipeline.py
+++ b/tests/test_preprocess_pipeline.py
@@ -410,7 +410,7 @@ def test_pad(L, dtype):
         np.testing.assert_allclose(pad_odd_to_even_many.images[:], ref)
 
         # Padding odd to odd by many should pad the first and last cols equally.
-        # This test will also excecise `fill_value`
+        # This test will also excercise `fill_value`
         many_even = 32
         fill = -1
         pad_odd_to_odd_many = sim1.crop_pad(L + many_even, fill_value=fill)
@@ -440,7 +440,7 @@ def test_pad(L, dtype):
         np.testing.assert_allclose(pad_even_to_even_many.images[:], ref)
 
         # Padding even to odd by many should pad the first and last+1 cols equally.
-        # This test will also excecise `fill_value`
+        # This test will also excercise `fill_value`
         many_odd = 33
         fill = -1
         pad_even_to_odd_many = sim1.crop_pad(L + many_odd, fill_value=fill)

From 74446177dd98caa2ff11ba5ee5b3af1a15f198a4 Mon Sep 17 00:00:00 2001
From: Garrett Wright <garrettwrong@gmail.com>
Date: Wed, 3 Sep 2025 09:29:34 -0400
Subject: [PATCH 23/23] Review string Updates

---
 src/aspire/image/image.py  | 2 +-
 src/aspire/source/image.py | 3 ---
 2 files changed, 1 insertion(+), 4 deletions(-)

diff --git a/src/aspire/image/image.py b/src/aspire/image/image.py
index 4b57a2c4df..185b9f798c 100644
--- a/src/aspire/image/image.py
+++ b/src/aspire/image/image.py
@@ -49,7 +49,7 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, shifted=False, ddof=0):
     input_dtype = imgs.dtype
 
     # Generate background mask
-    grid_dtype = np.float64
+    grid_dtype = np.float64  # Use doubles for accuracy and MATLAB repro
     grid = grid_2d(L, shifted=shifted, indexing="yx", dtype=grid_dtype)
     mask = grid["r"] > bg_radius
 
diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 48da7d259d..b38c667748 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -820,9 +820,6 @@ def crop_pad(self, L, fill_value=0):
         """
         Crop or pad images to size L.
 
-        Used for reproducing legacy MATLAB workflows.
-        For most applications, `downsample` is preferred.
-
         Cropping and padding makes no adjustments for centering conventions,
         but does maintain `pixel_size`.