Image Offsets Swap

src/aspire/image/image.py

@@ -510,7 +510,9 @@ def _im_translate(self, shifts):
             xp.ceil(xp.arange(-L / 2, L / 2, dtype=self.dtype))
         )
         grid_1d = grid_shifted * 2 * xp.pi / L
-        om_x, om_y = xp.meshgrid(grid_1d, grid_1d, indexing="ij")
+
+        # Grid indexing changed to "xy" to match Relion shift conventions.
+        om_x, om_y = xp.meshgrid(grid_1d, grid_1d, indexing="xy")
 
         phase_shifts_x = -shifts[:, 0].reshape((n_shifts, 1, 1))
         phase_shifts_y = -shifts[:, 1].reshape((n_shifts, 1, 1))

tests/saved_test_data/rln_proj_64_centered.star

@@ -0,0 +1,32 @@
+
+# version 30001
+
+data_optics
+
+loop_ 
+_rlnOpticsGroup #1 
+_rlnOpticsGroupName #2 
+_rlnVoltage #3 
+_rlnSphericalAberration #4 
+_rlnImagePixelSize #5 
+_rlnImageSize #6 
+_rlnImageDimensionality #7 
+           1    optics1   300.000000     2.700000     1.000000           64            2 
+
+
+# version 30001
+
+data_particles
+
+loop_ 
+_rlnAngleRot #1 
+_rlnAngleTilt #2 
+_rlnAnglePsi #3 
+_rlnOriginXAngst #4 
+_rlnOriginYAngst #5 
+_rlnOpticsGroup #6 
+_rlnImageName #7 
+  235.820138   113.086030    50.468981     0.000000     0.000000            1 000001@rln_proj_64_centered.mrcs
+   86.698555    31.958115   139.545228     0.000000     0.000000            1 000002@rln_proj_64_centered.mrcs
+   48.456166    71.176316   185.304830     0.000000     0.000000            1 000003@rln_proj_64_centered.mrcs
+  215.714386   105.017323   154.043384     0.000000     0.000000            1 000004@rln_proj_64_centered.mrcs

tests/saved_test_data/rln_proj_64_shifted.star

@@ -0,0 +1,32 @@
+
+# version 30001
+
+data_optics
+
+loop_ 
+_rlnOpticsGroup #1 
+_rlnOpticsGroupName #2 
+_rlnVoltage #3 
+_rlnSphericalAberration #4 
+_rlnImagePixelSize #5 
+_rlnImageSize #6 
+_rlnImageDimensionality #7 
+           1    optics1   300.000000     2.700000     1.000000           64            2
+
+
+# version 30001
+
+data_particles
+
+loop_ 
+_rlnAngleRot #1 
+_rlnAngleTilt #2 
+_rlnAnglePsi #3 
+_rlnOriginX #4 
+_rlnOriginY #5 
+_rlnOpticsGroup #6 
+_rlnImageName #7
+  235.820138   113.086030    50.468981     6.000000     10.000000            1 000001@rln_proj_64_shifted.mrcs
+   86.698555    31.958115   139.545228     10.000000    -5.000000            1 000002@rln_proj_64_shifted.mrcs
+   48.456166    71.176316   185.304830     -8.000000    11.000000            1 000003@rln_proj_64_shifted.mrcs
+  215.714386   105.017323   154.043384     -13.000000   -3.000000            1 000004@rln_proj_64_shifted.mrcs

tests/saved_test_data/rln_proj_65_centered.star

@@ -0,0 +1,35 @@
+
+# version 30001
+
+data_optics
+
+loop_ 
+_rlnOpticsGroup #1 
+_rlnOpticsGroupName #2 
+_rlnVoltage #3 
+_rlnSphericalAberration #4 
+_rlnImagePixelSize #5 
+_rlnImageSize #6 
+_rlnImageDimensionality #7 
+           1    optics1   300.000000     2.700000     1.000000           65            2 
+
+
+# version 30001
+
+data_particles
+
+loop_ 
+_rlnAngleRot #1 
+_rlnAngleTilt #2 
+_rlnAnglePsi #3 
+_rlnOriginXAngst #4 
+_rlnOriginYAngst #5 
+_rlnOpticsGroup #6 
+_rlnImageName #7 
+_rlnOriginX #8 
+_rlnOriginY #9 
+  235.820138   113.086030    50.468981     0.000000     0.000000       1 000001@rln_proj_65_centered.mrcs     0.00000     0.000000 
+   86.698555    31.958115   139.545228     0.000000     0.000000       1 000002@rln_proj_65_centered.mrcs     0.00000     0.000000 
+   48.456166    71.176316   185.304830     0.000000     0.000000       1 000003@rln_proj_65_centered.mrcs     0.00000     0.000000 
+  215.714386   105.017323   154.043384     0.000000     0.000000       1 000004@rln_proj_65_centered.mrcs     0.00000     0.000000 
+

tests/saved_test_data/rln_proj_65_shifted.star

@@ -0,0 +1,32 @@
+
+# version 30001
+
+data_optics
+
+loop_ 
+_rlnOpticsGroup #1 
+_rlnOpticsGroupName #2 
+_rlnVoltage #3 
+_rlnSphericalAberration #4 
+_rlnImagePixelSize #5 
+_rlnImageSize #6 
+_rlnImageDimensionality #7 
+           1    optics1   300.000000     2.700000     1.000000           65            2
+
+
+# version 30001
+
+data_particles
+
+loop_ 
+_rlnAngleRot #1 
+_rlnAngleTilt #2 
+_rlnAnglePsi #3 
+_rlnOriginX #4 
+_rlnOriginY #5 
+_rlnOpticsGroup #6 
+_rlnImageName #7
+  235.820138   113.086030    50.468981     6.000000     10.000000            1 000001@rln_proj_65_shifted.mrcs
+   86.698555    31.958115   139.545228     10.000000    -5.000000            1 000002@rln_proj_65_shifted.mrcs
+   48.456166    71.176316   185.304830     -8.000000    11.000000            1 000003@rln_proj_65_shifted.mrcs
+  215.714386   105.017323   154.043384     -13.000000   -3.000000            1 000004@rln_proj_65_shifted.mrcs

tests/test_anisotropic_noise.py

@@ -25,6 +25,10 @@ def setUp(self):
             dtype=self.dtype,
         )
 
+        # Keep hardcoded tests passing after fixing swapped offsets.
+        # See github issue #1146.
+        self.sim = self.sim.update(offsets=self.sim.offsets[:, [1, 0]])
+
     def tearDown(self):
         pass
 

tests/test_image.py

@@ -87,19 +87,24 @@ def testImShift(parity, dtype):
     im1 = im._im_translate(shifts)
     # test that float input returns the same thing
     im2 = im.shift(shifts.astype(dtype))
-    # ground truth numpy roll
-    im3 = np.roll(im_np[0, :, :], -shifts, axis=(0, 1))
+    # ground truth numpy roll.
+    # Note: NumPy axes 0 and 1 correspond to the row and column of an array,
+    # respectively, which corresponds to the y-axis and x-axis when that array
+    # represents an image. Since our shifts are (x-shifts, y-shifts), the axis
+    # parameter for np.roll() must be set to (1, 0) to accomodate.
+    im3 = np.roll(im_np[0, :, :], -shifts, axis=(1, 0))
 
     atol = utest_tolerance(dtype)
 
-    assert np.allclose(im0.asnumpy(), im1.asnumpy(), atol=atol)
-    assert np.allclose(im1.asnumpy(), im2.asnumpy(), atol=atol)
-    assert np.allclose(im0.asnumpy()[0, :, :], im3, atol=atol)
+    np.testing.assert_allclose(im0.asnumpy(), im1.asnumpy(), atol=atol)
+    np.testing.assert_allclose(im1.asnumpy(), im2.asnumpy(), atol=atol)
+    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)
+
     # test stack of shifts (same number as Image.num_img)
     # mix of odd and even
     shifts = np.array([[100, 200], [203, 150], [55, 307]])
@@ -111,15 +116,19 @@ def testImShiftStack(parity, dtype):
     # test that float input returns the same thing
     im2 = ims.shift(shifts.astype(dtype))
     # ground truth numpy roll
+    # Note: NumPy axes 0 and 1 correspond to the row and column of an array,
+    # respectively, which corresponds to the y-axis and x-axis when that array
+    # represents an image. Since our shifts are (x-shifts, y-shifts), the axis
+    # parameter for np.roll() must be set to (1, 0) to accomodate.
     im3 = np.array(
-        [np.roll(ims_np[i, :, :], -shifts[i], axis=(0, 1)) for i in range(n)]
+        [np.roll(ims_np[i, :, :], -shifts[i], axis=(1, 0)) for i in range(n)]
     )
 
     atol = utest_tolerance(dtype)
 
-    assert np.allclose(im0.asnumpy(), im1.asnumpy(), atol=atol)
-    assert np.allclose(im1.asnumpy(), im2.asnumpy(), atol=atol)
-    assert np.allclose(im0.asnumpy(), im3, atol=atol)
+    np.testing.assert_allclose(im0.asnumpy(), im1.asnumpy(), atol=atol)
+    np.testing.assert_allclose(im1.asnumpy(), im2.asnumpy(), atol=atol)
+    np.testing.assert_allclose(im0.asnumpy(), im3, atol=atol)
 
 
 def testImageShiftErrors():

tests/test_relion_interop.py

@@ -4,16 +4,29 @@
 import pytest
 
 from aspire.source import RelionSource, Simulation
+from aspire.utils import utest_tolerance
 from aspire.volume import Volume
 
 DATA_DIR = os.path.join(os.path.dirname(__file__), "saved_test_data")
 
 
-STARFILE = ["rln_proj_65.star", "rln_proj_64.star"]
+STARFILE_ODD = [
+    "rln_proj_65_centered.star",
+    "rln_proj_65_shifted.star",
+]
 
+STARFILE_EVEN = [
+    "rln_proj_64_centered.star",
+    "rln_proj_64_shifted.star",
+]
 
-@pytest.fixture(params=STARFILE, scope="module")
+
+@pytest.fixture(params=STARFILE_ODD + STARFILE_EVEN, scope="module")
 def sources(request):
+    """
+    Initialize RelionSource from starfile and generate corresponding ASPIRE
+    Simulation source.
+    """
     starfile = os.path.join(DATA_DIR, request.param)
     rln_src = RelionSource(starfile)
 
@@ -24,9 +37,9 @@ def sources(request):
 
     # Create Simulation source using Volume and angles from Relion projections.
     # Note, for odd resolution Relion projections are shifted by 1 pixel in x and y.
-    offsets = 0
+    offsets = rln_src.offsets
     if rln_src.L % 2 == 1:
-        offsets = -np.ones((rln_src.n, 2), dtype=rln_src.dtype)
+        offsets -= np.ones((rln_src.n, 2), dtype=rln_src.dtype)
 
     sim_src = Simulation(
         n=rln_src.n,
@@ -39,6 +52,21 @@ def sources(request):
     return rln_src, sim_src
 
 
+@pytest.fixture(params=[STARFILE_ODD, STARFILE_EVEN], scope="module")
+def rln_sources(request):
+    """
+    Initialize centered and shifted RelionSource's generated using the
+    same viewing angles.
+    """
+    starfile_centered = os.path.join(DATA_DIR, request.param[0])
+    starfile_shifted = os.path.join(DATA_DIR, request.param[1])
+
+    rln_src_centered = RelionSource(starfile_centered)
+    rln_src_shifted = RelionSource(starfile_shifted)
+
+    return rln_src_centered, rln_src_shifted
+
+
 def test_projections_relative_error(sources):
     """Check the relative error between Relion and ASPIRE projection images."""
     rln_src, sim_src = sources
@@ -51,11 +79,11 @@ def test_projections_relative_error(sources):
     rln_np = (rln_np - np.mean(rln_np)) / np.std(rln_np)
     sim_np = (sim_np - np.mean(sim_np)) / np.std(sim_np)
 
-    # Check that relative error is less than 3%.
+    # Check that relative error is less than 4%.
     error = np.linalg.norm(rln_np - sim_np, axis=(1, 2)) / np.linalg.norm(
         rln_np, axis=(1, 2)
     )
-    np.testing.assert_array_less(error, 0.03)
+    np.testing.assert_array_less(error, 0.04)
 
 
 def test_projections_frc(sources):
@@ -67,4 +95,18 @@ def test_projections_frc(sources):
 
     # Check that estimated resolution is high (< 2.5 pixels) and correlation is close to 1.
     np.testing.assert_array_less(res, 2.5)
-    np.testing.assert_array_less(1 - corr[:, -2], 0.02)
+    np.testing.assert_array_less(1 - corr[:, -2], 0.025)
+
+
+def test_relion_source_centering(rln_sources):
+    """Test that centering by using provided Relion shifts works."""
+    rln_src_centered, rln_src_shifted = rln_sources
+    ims_centered = rln_src_centered.images[:]
+    ims_shifted = rln_src_shifted.images[:]
+
+    offsets = rln_src_shifted.offsets
+    np.testing.assert_allclose(
+        ims_centered.asnumpy(),
+        ims_shifted.shift(-offsets).asnumpy(),
+        atol=utest_tolerance(rln_src_centered.dtype),
+    )

tests/test_simulation.py

@@ -123,6 +123,10 @@ def setUp(self):
             dtype=self.dtype,
         )
 
+        # Keep hardcoded tests passing after fixing swapped offsets.
+        # See github issue #1146.
+        self.sim = self.sim.update(offsets=self.sim.offsets[:, [1, 0]])
+
     def tearDown(self):
         pass
 
@@ -162,6 +166,7 @@ def testSimulationCached(self):
             n=self.n,
             L=self.L,
             vols=self.vols,
+            offsets=self.sim.offsets,
             unique_filters=[
                 RadialCTFFilter(defocus=d) for d in np.linspace(1.5e4, 2.5e4, 7)
             ],