diff --git a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
index da306d546c..dee01c4868 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10028_jsb.py
@@ -88,12 +88,15 @@
 logger.info("Perform phase flip to input images.")
 src = src.phase_flip().cache()
 
+# Legacy MATLAB cropped the images to an odd resolution.
+src = src.crop_pad(src.L - 1).cache()
+
 # 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..43daaeb7bc 100644
--- a/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
+++ b/gallery/experiments/experimental_abinitio_pipeline_10073_jsb.py
@@ -82,12 +82,15 @@
 logger.info("Perform phase flip to input images.")
 src = src.phase_flip().cache()
 
+# Legacy MATLAB cropped the images to an odd resolution.
+src = src.crop_pad(src.L - 1).cache()
+
 # 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/src/aspire/image/image.py b/src/aspire/image/image.py
index ed061ecec7..185b9f798c 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:
@@ -44,19 +46,11 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, legacy=False):
             "`normalize_bg` is currently limited to 1D image stacks."
         )
     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
+    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  # Use doubles for accuracy and MATLAB repro
+    grid = grid_2d(L, shifted=shifted, indexing="yx", dtype=grid_dtype)
     mask = grid["r"] > bg_radius
 
     if do_ramp:
@@ -66,14 +60,14 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, legacy=False):
             (
                 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))
@@ -85,10 +79,14 @@ def normalize_bg(imgs, bg_radius=1.0, do_ramp=True, legacy=False):
         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)
+    # 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]
 
-    return (imgs - mean[:, None, None]) / std[:, None, None]
+    # Restore input dtype
+    return imgs.astype(input_dtype, copy=False)
 
 
 def load_mrc(filepath):
@@ -490,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)
@@ -515,12 +514,12 @@ 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)
 
-    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.
 
@@ -528,8 +527,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.
         """
 
@@ -540,25 +540,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/image/xform.py b/src/aspire/image/xform.py
index 10fddf1df2..316d8ca3c8 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__)
 
@@ -199,7 +200,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 +208,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 +228,30 @@ 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 CropPad(Xform):
+    """
+    A Xform that crops or pads an Image object to a specified size.
+    """
+
+    def __init__(self, L, fill_value=0):
+        """
+        Initialize Xform to crop Image to a specific size.
+
+        :param L: int - new size
+        :param fill_value: Optional value for padding, default 0.
+        """
+        self.L = L
+        self.fill_value = fill_value
+        super().__init__()
+
+    def _forward(self, im, indices):
+        return crop_pad_2d(im, self.L, self.fill_value)
+
+    def __str__(self):
+        return f"CropPad({self.L}, {self.fill_value}) Xform"
 
 
 class LegacyWhiten(Xform):
@@ -256,7 +282,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/src/aspire/source/image.py b/src/aspire/source/image.py
index f4f30bf481..b38c667748 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 (
+    CropPad,
     Downsample,
     FilterXform,
     IndexedXform,
@@ -768,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 workflow downsampling.
+
+        For uses other than MALTAB reproduction, prefer `downsample`.
+
+        :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}."
@@ -777,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
@@ -788,6 +815,41 @@ def downsample(self, L, zero_nyquist=True, legacy=False):
 
         self.L = L
 
+    @_as_copy
+    def crop_pad(self, L, fill_value=0):
+        """
+        Crop or pad images to size L.
+
+        Cropping and padding makes no adjustments for centering conventions,
+        but does maintain `pixel_size`.
+
+        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 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:
+            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}"
+            )
+        else:
+            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))
+
+        self.L = L
+
     @_as_copy
     def whiten(self, noise_estimate=None, epsilon=None):
         """
@@ -944,25 +1006,41 @@ 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):
+        """
+        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:
+        # 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
+        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.
-        :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.
+        :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(
@@ -971,7 +1049,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,
             )
         )
 
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 94c5c1dae4..ebd6e1c593 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.normalize_background(legacy=True)
+    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)
@@ -114,6 +123,37 @@ def test_norm_background_legacy(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.
@@ -269,3 +309,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 excercise `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 excercise `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)