From f7099292452007fdfe95747bbefb27ecfbf48b49 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 16 Oct 2025 16:00:37 -0400
Subject: [PATCH 01/25] initial ImageSource.save with optics block

---
 src/aspire/source/image.py      | 89 +++++++++++++++++++++++++++++++--
 src/aspire/storage/starfile.py  |  1 +
 tests/test_coordinate_source.py |  2 +-
 3 files changed, 87 insertions(+), 5 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index f1dff74404..3395dd854d 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1289,6 +1289,79 @@ def _populate_local_metadata(self):
         """
         return []
 
+    @staticmethod
+    def _prepare_relion_optics_blocks(metadata):
+        """
+        Split metadata into RELION>=3.1 style `data_optics` and `data_particles` blocks.
+
+        The optics block has one row per optics group with:
+        `_rlnOpticsGroup`, `_rlnOpticsGroupName`, and optics metadata columns.
+        The particle block keeps the remaining columns and includes a per-particle
+        `_rlnOpticsGroup` that references the optics block.
+        """
+        # All possible optics group fields
+        # TODO: Check we have all of them
+        all_optics_fields = [
+            "_rlnImagePixelSize",
+            "_rlnMicrographPixelSize",
+            "_rlnSphericalAberration",
+            "_rlnVoltage",
+            "_rlnImageSize",
+            "_rlnAmplitudeContrast",
+        ]
+
+        # TODO: Need to add _rlnImageSize here and above
+        required_fields = ["_rlnSphericalAberration", "_rlnVoltage"]
+        pixel_fields = ["_rlnImagePixelSize", "_rlnMicrographPixelSize"]
+
+        has_required = all(field in metadata for field in required_fields)
+        has_pixel_field = any(field in metadata for field in pixel_fields)
+
+        # TODO: Add warning?
+        if not (has_required and has_pixel_field):
+            # Optics metadata incomplete, fall back to legacy single block.
+            return None, metadata
+
+        # Collect all optics group fields present in metadata and determine unique optics groups.
+        optics_value_fields = [
+            field for field in all_optics_fields if field in metadata
+        ]
+        n_rows = len(metadata["_rlnImagePixelSize"])
+
+        group_lookup = OrderedDict()  # Stores distinct optics groups
+        optics_groups = np.empty(n_rows, dtype=int)
+
+        for idx in range(n_rows):
+            signature = tuple(metadata[field][idx] for field in optics_value_fields)
+            if signature not in group_lookup:
+                group_lookup[signature] = len(group_lookup) + 1
+            optics_groups[idx] = group_lookup[signature]
+
+        metadata["_rlnOpticsGroup"] = optics_groups
+
+        optics_block = OrderedDict()
+        optics_block["_rlnOpticsGroup"] = []
+        optics_block["_rlnOpticsGroupName"] = []
+        for field in optics_value_fields:
+            optics_block[field] = []
+
+        for signature, group_id in group_lookup.items():
+            optics_block["_rlnOpticsGroup"].append(group_id)
+            optics_block["_rlnOpticsGroupName"].append(f"opticsGroup{group_id}")
+            for field, value in zip(optics_value_fields, signature):
+                optics_block[field].append(value)
+
+        # Collect particle_block metadata
+        particle_block = OrderedDict()
+        if "_rlnOpticsGroup" in metadata:
+            particle_block["_rlnOpticsGroup"] = metadata["_rlnOpticsGroup"]
+        for key, value in metadata.items():
+            if key in optics_value_fields or key == "_rlnOpticsGroup":
+                continue
+            particle_block[key] = value
+
+        return optics_block, particle_block
+
     def save_metadata(self, starfile_filepath, batch_size=512, save_mode=None):
         """
         Save updated metadata to a STAR file
@@ -1324,12 +1397,20 @@ def save_metadata(self, starfile_filepath, batch_size=512, save_mode=None):
             for x in np.char.split(metadata["_rlnImageName"].astype(np.str_), sep="@")
         ]
 
+        # Separate metadata into optics and particle blocks
+        optics_block, particle_block = self._prepare_relion_optics_blocks(metadata)
+
         # initialize the star file object and save it
         odict = OrderedDict()
-        # since our StarFile only has one block, the convention is to save it with the header "data_", i.e. its name is blank
-        # if we had a block called "XYZ" it would be saved as "XYZ"
-        # thus we index the metadata block with ""
-        odict[""] = metadata
+
+        # StarFile uses the `odict` keys to label the starfile block headers "data_(key)". Following RELION>=3.1
+        # convention we label the blocks "data_optics" and "data_particles".
+        if optics_block is None:
+            odict["particles"] = particle_block
+        else:
+            odict["optics"] = optics_block
+            odict["particles"] = particle_block
+
         out_star = StarFile(blocks=odict)
         out_star.write(starfile_filepath)
         return filename_indices
diff --git a/src/aspire/storage/starfile.py b/src/aspire/storage/starfile.py
index 98b3219607..9b55610935 100644
--- a/src/aspire/storage/starfile.py
+++ b/src/aspire/storage/starfile.py
@@ -135,6 +135,7 @@ def write(self, filepath):
         # create an empty Document
         _doc = cif.Document()
         filepath = str(filepath)
+
         for name, block in self.blocks.items():
             # construct new empty block
             _block = _doc.add_new_block(name)
diff --git a/tests/test_coordinate_source.py b/tests/test_coordinate_source.py
index d2f93d5f65..783790f145 100644
--- a/tests/test_coordinate_source.py
+++ b/tests/test_coordinate_source.py
@@ -535,7 +535,7 @@ def testSave(self):
             self.assertTrue(np.array_equal(imgs.asnumpy()[i], saved_mrcs_stack[i]))
         # assert that the star file has the correct metadata
         self.assertEqual(
-            list(saved_star[""].keys()),
+            list(saved_star["particles"].keys()),
             [
                 "_rlnImagePixelSize",
                 "_rlnSymmetryGroup",

From 133a3d7b2dfc47b8ed70f82203ee07fdb7cac76a Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 23 Oct 2025 16:02:32 -0400
Subject: [PATCH 02/25] add _rlnImageSize column to metadata when saving.

---
 src/aspire/source/image.py     | 7 +++++--
 src/aspire/storage/starfile.py | 1 -
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 3395dd854d..b9d38d67cc 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1310,8 +1310,7 @@ def _prepare_relion_optics_blocks(metadata):
             "_rlnAmplitudeContrast",
         ]
 
-        # TODO: Need to add _rlnImageSize here and above
-        required_fields = ["_rlnSphericalAberration", "_rlnVoltage"]
+        required_fields = ["_rlnSphericalAberration", "_rlnVoltage", "_rlnImageSize"]
         pixel_fields = ["_rlnImagePixelSize", "_rlnMicrographPixelSize"]
 
         has_required = all(field in metadata for field in required_fields)
@@ -1397,6 +1396,10 @@ def save_metadata(self, starfile_filepath, batch_size=512, save_mode=None):
             for x in np.char.split(metadata["_rlnImageName"].astype(np.str_), sep="@")
         ]
 
+        # Populate _rlnImageSize column, required for optics_block below
+        if "_rlnImageSize" not in metadata:
+            metadata["_rlnImageSize"] = np.full(self.n, self.L, dtype=int)
+
         # Separate metadata into optics and particle blocks
         optics_block, particle_block = self._prepare_relion_optics_blocks(metadata)
 
diff --git a/src/aspire/storage/starfile.py b/src/aspire/storage/starfile.py
index 9b55610935..98b3219607 100644
--- a/src/aspire/storage/starfile.py
+++ b/src/aspire/storage/starfile.py
@@ -135,7 +135,6 @@ def write(self, filepath):
         # create an empty Document
         _doc = cif.Document()
         filepath = str(filepath)
-
         for name, block in self.blocks.items():
             # construct new empty block
             _block = _doc.add_new_block(name)

From bc01f4ced21e483b73ad984cd6585d838566600e Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 24 Oct 2025 08:50:09 -0400
Subject: [PATCH 03/25] Add _rlnImageSize to coordinate source test.

---
 tests/test_coordinate_source.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/test_coordinate_source.py b/tests/test_coordinate_source.py
index 783790f145..c8fb20ca8e 100644
--- a/tests/test_coordinate_source.py
+++ b/tests/test_coordinate_source.py
@@ -542,6 +542,7 @@ def testSave(self):
                 "_rlnImageName",
                 "_rlnCoordinateX",
                 "_rlnCoordinateY",
+                "_rlnImageSize",
             ],
         )
         # assert that all the correct coordinates were saved

From 0b1204f9ad171dc8d0b1eca994d25626427a4eba Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 24 Oct 2025 09:53:37 -0400
Subject: [PATCH 04/25] use rlnImageName to get n_rows

---
 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 b9d38d67cc..cdfa4205be 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1325,7 +1325,7 @@ def _prepare_relion_optics_blocks(metadata):
         optics_value_fields = [
             field for field in all_optics_fields if field in metadata
         ]
-        n_rows = len(metadata["_rlnImagePixelSize"])
+        n_rows = len(metadata["_rlnImageName"])
 
         group_lookup = OrderedDict()  # Stores distinct optics groups
         optics_groups = np.empty(n_rows, dtype=int)

From e4e24b0ad6b1044bc7ac200b308a20f81e022685 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 24 Oct 2025 16:02:54 -0400
Subject: [PATCH 05/25] test for sim save with optics block

---
 src/aspire/utils/relion_interop.py |  1 +
 tests/test_simulation.py           | 55 +++++++++++++++++++++++++++++-
 2 files changed, 55 insertions(+), 1 deletion(-)

diff --git a/src/aspire/utils/relion_interop.py b/src/aspire/utils/relion_interop.py
index 996d321366..90e9a555b2 100644
--- a/src/aspire/utils/relion_interop.py
+++ b/src/aspire/utils/relion_interop.py
@@ -21,6 +21,7 @@
     "_rlnCtfFigureOfMerit": float,
     "_rlnMagnification": float,
     "_rlnImagePixelSize": float,
+    "_rlnImageSize": int,
     "_rlnAmplitudeContrast": float,
     "_rlnImageName": str,
     "_rlnOriginalName": str,
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 5859aae934..082ae8e59a 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -9,7 +9,7 @@
 from aspire.noise import WhiteNoiseAdder
 from aspire.operators import RadialCTFFilter
 from aspire.source import RelionSource, Simulation, _LegacySimulation
-from aspire.utils import utest_tolerance
+from aspire.utils import RelionStarFile, utest_tolerance
 from aspire.volume import LegacyVolume, SymmetryGroup, Volume
 
 from .test_utils import matplotlib_dry_run
@@ -627,6 +627,59 @@ def testSimulationSaveFile(self):
             )
 
 
+def test_simulation_save_optics_block(tmp_path):
+    res = 32
+
+    # Radial CTF Filters. Should make 3 distinct optics blocks
+    kv_min, kv_max, kv_ct = 200, 300, 3
+    voltages = np.linspace(kv_min, kv_max, kv_ct)
+    ctf_filters = [RadialCTFFilter(voltage=kv) for kv in voltages]
+
+    # Generate and save Simulation
+    sim = Simulation(
+        n=9, L=res, C=1, unique_filters=ctf_filters, pixel_size=1.34
+    ).cache()
+    starpath = tmp_path / "sim.star"
+    sim.save(starpath, overwrite=True)
+
+    star = RelionStarFile(str(starpath))
+    assert star.relion_version == "3.1"
+    assert set(star.blocks.keys()) == {"optics", "particles"}
+
+    optics = star["optics"]
+    expected_optics_fields = [
+        "_rlnOpticsGroup",
+        "_rlnOpticsGroupName",
+        "_rlnImagePixelSize",
+        "_rlnSphericalAberration",
+        "_rlnVoltage",
+        "_rlnImageSize",
+        "_rlnAmplitudeContrast",
+    ]
+    for field in expected_optics_fields:
+        assert field in optics
+
+    np.testing.assert_array_equal(
+        optics["_rlnOpticsGroup"], np.arange(1, kv_ct + 1, dtype=int)
+    )
+    np.testing.assert_array_equal(
+        optics["_rlnOpticsGroupName"],
+        np.array([f"opticsGroup{i}" for i in range(1, kv_ct + 1)], dtype=object),
+    )
+    np.testing.assert_array_equal(optics["_rlnImageSize"], np.full(kv_ct, res))
+
+    # Depending on Simulation random indexing, voltages will be unordered
+    np.testing.assert_allclose(np.sort(optics["_rlnVoltage"]), voltages)
+
+    particles = star["particles"]
+    assert "_rlnOpticsGroup" in particles
+    assert len(particles["_rlnOpticsGroup"]) == sim.n
+    np.testing.assert_array_equal(
+        np.sort(np.unique(particles["_rlnOpticsGroup"])),
+        np.arange(1, kv_ct + 1, dtype=int),
+    )
+
+
 def test_default_symmetry_group():
     # Check that default is "C1".
     sim = Simulation()

From 7c88e0a90b5f3f8c8a371e0afb1948ef04fe508f Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 27 Oct 2025 10:48:42 -0400
Subject: [PATCH 06/25] add _rlnImageDimensionality

---
 src/aspire/source/image.py         | 16 +++++++++++++---
 src/aspire/utils/relion_interop.py |  1 +
 tests/test_simulation.py           | 12 +++++++++++-
 3 files changed, 25 insertions(+), 4 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index cdfa4205be..cf369c584b 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1306,11 +1306,18 @@ def _prepare_relion_optics_blocks(metadata):
             "_rlnMicrographPixelSize",
             "_rlnSphericalAberration",
             "_rlnVoltage",
-            "_rlnImageSize",
             "_rlnAmplitudeContrast",
+            "_rlnImageSize",
+            "_rlnImageDimensionality",
         ]
 
-        required_fields = ["_rlnSphericalAberration", "_rlnVoltage", "_rlnImageSize"]
+        required_fields = [
+            "_rlnSphericalAberration",
+            "_rlnVoltage",
+            "_rlnAmplitudeContrast",
+            "_rlnImageSize",
+            "_rlnImageDimensionality",
+        ]
         pixel_fields = ["_rlnImagePixelSize", "_rlnMicrographPixelSize"]
 
         has_required = all(field in metadata for field in required_fields)
@@ -1396,10 +1403,13 @@ def save_metadata(self, starfile_filepath, batch_size=512, save_mode=None):
             for x in np.char.split(metadata["_rlnImageName"].astype(np.str_), sep="@")
         ]
 
-        # Populate _rlnImageSize column, required for optics_block below
+        # Populate _rlnImageSize, _rlnImageDimensionality columns, required for optics_block below
         if "_rlnImageSize" not in metadata:
             metadata["_rlnImageSize"] = np.full(self.n, self.L, dtype=int)
 
+        if "_rlnImageDimensionality" not in metadata:
+            metadata["_rlnImageDimensionality"] = np.full(self.n, 2, dtype=int)
+
         # Separate metadata into optics and particle blocks
         optics_block, particle_block = self._prepare_relion_optics_blocks(metadata)
 
diff --git a/src/aspire/utils/relion_interop.py b/src/aspire/utils/relion_interop.py
index 90e9a555b2..c807c7d863 100644
--- a/src/aspire/utils/relion_interop.py
+++ b/src/aspire/utils/relion_interop.py
@@ -20,6 +20,7 @@
     "_rlnDetectorPixelSize": float,
     "_rlnCtfFigureOfMerit": float,
     "_rlnMagnification": float,
+    "_rlnImageDimensionality": int,
     "_rlnImagePixelSize": float,
     "_rlnImageSize": int,
     "_rlnAmplitudeContrast": float,
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 082ae8e59a..c1f273e082 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -653,12 +653,16 @@ def test_simulation_save_optics_block(tmp_path):
         "_rlnImagePixelSize",
         "_rlnSphericalAberration",
         "_rlnVoltage",
-        "_rlnImageSize",
         "_rlnAmplitudeContrast",
+        "_rlnImageSize",
+        "_rlnImageDimensionality",
     ]
+
+    # Check all required fields are present
     for field in expected_optics_fields:
         assert field in optics
 
+    # Optics group and group name should 1-indexed
     np.testing.assert_array_equal(
         optics["_rlnOpticsGroup"], np.arange(1, kv_ct + 1, dtype=int)
     )
@@ -666,11 +670,17 @@ def test_simulation_save_optics_block(tmp_path):
         optics["_rlnOpticsGroupName"],
         np.array([f"opticsGroup{i}" for i in range(1, kv_ct + 1)], dtype=object),
     )
+
+    # Check image size and image dimensionality
     np.testing.assert_array_equal(optics["_rlnImageSize"], np.full(kv_ct, res))
 
+    optics_dim = np.array(optics["_rlnImageDimensionality"], dtype=int)
+    np.testing.assert_array_equal(optics_dim, np.full(len(optics_dim), 2))
+
     # Depending on Simulation random indexing, voltages will be unordered
     np.testing.assert_allclose(np.sort(optics["_rlnVoltage"]), voltages)
 
+    # Check that each row of the data_particles block has an associated optics group
     particles = star["particles"]
     assert "_rlnOpticsGroup" in particles
     assert len(particles["_rlnOpticsGroup"]) == sim.n

From 3e78f9ac54a4585d80d53aded5dee01a55a62256 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 27 Oct 2025 11:39:52 -0400
Subject: [PATCH 07/25] test save/load roundtrip w/ phase_flip.

---
 tests/test_simulation.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index c1f273e082..46b9c5603d 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -689,6 +689,12 @@ def test_simulation_save_optics_block(tmp_path):
         np.arange(1, kv_ct + 1, dtype=int),
     )
 
+    # Test phase_flip after save/load round trip to ensure correct optics group mapping
+    rln_src = RelionSource(starpath)
+    np.testing.assert_allclose(
+        sim.phase_flip().images[:], rln_src.phase_flip().images[:]
+    )
+
 
 def test_default_symmetry_group():
     # Check that default is "C1".

From 816923c8e6171ff95e1f89edb21b121505eb2fa6 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 27 Oct 2025 13:14:35 -0400
Subject: [PATCH 08/25] update coord source test

---
 tests/test_coordinate_source.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/test_coordinate_source.py b/tests/test_coordinate_source.py
index c8fb20ca8e..2a63c7b0b9 100644
--- a/tests/test_coordinate_source.py
+++ b/tests/test_coordinate_source.py
@@ -543,6 +543,7 @@ def testSave(self):
                 "_rlnCoordinateX",
                 "_rlnCoordinateY",
                 "_rlnImageSize",
+                "_rlnImageDimensionality",
             ],
         )
         # assert that all the correct coordinates were saved

From 85b9fd4c7437251efc4a7f0b6d84a8c0b8647978 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 27 Oct 2025 15:48:37 -0400
Subject: [PATCH 09/25] cleanup

---
 tests/test_simulation.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 46b9c5603d..70fd75a339 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -644,7 +644,7 @@ def test_simulation_save_optics_block(tmp_path):
 
     star = RelionStarFile(str(starpath))
     assert star.relion_version == "3.1"
-    assert set(star.blocks.keys()) == {"optics", "particles"}
+    assert star.blocks.keys() == {"optics", "particles"}
 
     optics = star["optics"]
     expected_optics_fields = [
@@ -668,16 +668,16 @@ def test_simulation_save_optics_block(tmp_path):
     )
     np.testing.assert_array_equal(
         optics["_rlnOpticsGroupName"],
-        np.array([f"opticsGroup{i}" for i in range(1, kv_ct + 1)], dtype=object),
+        np.array([f"opticsGroup{i}" for i in range(1, kv_ct + 1)]),
     )
 
-    # Check image size and image dimensionality
+    # Check image size (res) and image dimensionality (2)
     np.testing.assert_array_equal(optics["_rlnImageSize"], np.full(kv_ct, res))
+    np.testing.assert_array_equal(
+        optics["_rlnImageDimensionality"], np.full(len(optics_dim), 2)
+    )
 
-    optics_dim = np.array(optics["_rlnImageDimensionality"], dtype=int)
-    np.testing.assert_array_equal(optics_dim, np.full(len(optics_dim), 2))
-
-    # Depending on Simulation random indexing, voltages will be unordered
+    # Due to Simulation random indexing, voltages will be unordered
     np.testing.assert_allclose(np.sort(optics["_rlnVoltage"]), voltages)
 
     # Check that each row of the data_particles block has an associated optics group

From 10056b1bfd7c86e8fdb39a58d6173b0604641117 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 28 Oct 2025 09:14:24 -0400
Subject: [PATCH 10/25] remove unused var

---
 tests/test_simulation.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 70fd75a339..944136c999 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -673,9 +673,7 @@ def test_simulation_save_optics_block(tmp_path):
 
     # Check image size (res) and image dimensionality (2)
     np.testing.assert_array_equal(optics["_rlnImageSize"], np.full(kv_ct, res))
-    np.testing.assert_array_equal(
-        optics["_rlnImageDimensionality"], np.full(len(optics_dim), 2)
-    )
+    np.testing.assert_array_equal(optics["_rlnImageDimensionality"], np.full(kv_ct, 2))
 
     # Due to Simulation random indexing, voltages will be unordered
     np.testing.assert_allclose(np.sort(optics["_rlnVoltage"]), voltages)

From 6a846443a270ce23dae3da83695ab29359120b0c Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 28 Oct 2025 10:27:40 -0400
Subject: [PATCH 11/25] missing optics fields warning. warning test. code
 comments.

---
 src/aspire/source/image.py  | 16 ++++++++++------
 tests/test_relion_source.py | 24 +++++++++++++++++++++++-
 tests/test_simulation.py    |  2 +-
 3 files changed, 34 insertions(+), 8 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index cf369c584b..02f4942d74 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1299,8 +1299,7 @@ def _prepare_relion_optics_blocks(metadata):
         The particle block keeps the remaining columns and includes a per-particle
         `_rlnOpticsGroup` that references the optics block.
         """
-        # All possible optics group fields
-        # TODO: Check we have all of them
+        # Columns that belong in RELION's optics table.
         all_optics_fields = [
             "_rlnImagePixelSize",
             "_rlnMicrographPixelSize",
@@ -1311,6 +1310,7 @@ def _prepare_relion_optics_blocks(metadata):
             "_rlnImageDimensionality",
         ]
 
+        # We only write a data_optics block when every required field is present.
         required_fields = [
             "_rlnSphericalAberration",
             "_rlnVoltage",
@@ -1323,18 +1323,21 @@ def _prepare_relion_optics_blocks(metadata):
         has_required = all(field in metadata for field in required_fields)
         has_pixel_field = any(field in metadata for field in pixel_fields)
 
-        # TODO: Add warning?
         if not (has_required and has_pixel_field):
             # Optics metadata incomplete, fall back to legacy single block.
+            logger.warning(
+                "Optics metadata incomplete, writing only data_particles block."
+            )
             return None, metadata
 
-        # Collect all optics group fields present in metadata and determine unique optics groups.
+        # Restrict to the optics columns that are actually present on this source.
         optics_value_fields = [
             field for field in all_optics_fields if field in metadata
         ]
         n_rows = len(metadata["_rlnImageName"])
 
-        group_lookup = OrderedDict()  # Stores distinct optics groups
+        # Map each unique optics tuple to a 1-based group ID in order encountered.
+        group_lookup = OrderedDict()
         optics_groups = np.empty(n_rows, dtype=int)
 
         for idx in range(n_rows):
@@ -1345,6 +1348,7 @@ def _prepare_relion_optics_blocks(metadata):
 
         metadata["_rlnOpticsGroup"] = optics_groups
 
+        # Build the optics block rows and assign group names.
         optics_block = OrderedDict()
         optics_block["_rlnOpticsGroup"] = []
         optics_block["_rlnOpticsGroupName"] = []
@@ -1357,7 +1361,7 @@ def _prepare_relion_optics_blocks(metadata):
             for field, value in zip(optics_value_fields, signature):
                 optics_block[field].append(value)
 
-        # Collect particle_block metadata
+        # Everything not lifted into the optics block stays with the particle metadata.
         particle_block = OrderedDict()
         if "_rlnOpticsGroup" in metadata:
             particle_block["_rlnOpticsGroup"] = metadata["_rlnOpticsGroup"]
diff --git a/tests/test_relion_source.py b/tests/test_relion_source.py
index 009ecd321d..d0b996c795 100644
--- a/tests/test_relion_source.py
+++ b/tests/test_relion_source.py
@@ -5,7 +5,7 @@
 import numpy as np
 import pytest
 
-from aspire.source import RelionSource, Simulation
+from aspire.source import ImageSource, RelionSource, Simulation
 from aspire.utils import RelionStarFile
 from aspire.volume import SymmetryGroup
 
@@ -61,6 +61,28 @@ def test_symmetry_group(caplog):
     assert str(src_override_sym.symmetry_group) == "C6"
 
 
+def test_prepare_relion_optics_blocks_warns(caplog):
+    """
+    Test we warn when optics group metadata is missing.
+    """
+    # metadata dict with no CTF values
+    metadata = {
+        "_rlnImagePixelSize": np.array([1.234]),
+        "_rlnImageSize": np.array([32]),
+        "_rlnImageDimensionality": np.array([2]),
+    }
+
+    caplog.clear()
+    with caplog.at_level(logging.WARNING):
+        optics_block, particle_block = ImageSource._prepare_relion_optics_blocks(
+            metadata.copy()
+        )
+
+    assert optics_block is None
+    assert particle_block == metadata
+    assert "Optics metadata incomplete" in caplog.text
+
+
 def test_pixel_size(caplog):
     """
     Instantiate RelionSource from starfiles containing the following pixel size
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 944136c999..2f323d3490 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -8,7 +8,7 @@
 
 from aspire.noise import WhiteNoiseAdder
 from aspire.operators import RadialCTFFilter
-from aspire.source import RelionSource, Simulation, _LegacySimulation
+from aspire.source import ImageSource, RelionSource, Simulation, _LegacySimulation
 from aspire.utils import RelionStarFile, utest_tolerance
 from aspire.volume import LegacyVolume, SymmetryGroup, Volume
 

From 35b053244be50add3432539f5fdfaefd817f9e9b Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Tue, 28 Oct 2025 10:28:34 -0400
Subject: [PATCH 12/25] removed unused import.

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

diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 2f323d3490..944136c999 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -8,7 +8,7 @@
 
 from aspire.noise import WhiteNoiseAdder
 from aspire.operators import RadialCTFFilter
-from aspire.source import ImageSource, RelionSource, Simulation, _LegacySimulation
+from aspire.source import RelionSource, Simulation, _LegacySimulation
 from aspire.utils import RelionStarFile, utest_tolerance
 from aspire.volume import LegacyVolume, SymmetryGroup, Volume
 

From 4da7ee831705e6de90572881b147f4088c42ce06 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 30 Oct 2025 11:24:10 -0400
Subject: [PATCH 13/25] deprecate _rlnAmplitude in favor of private __amplitude
 field.

---
 src/aspire/source/image.py         | 16 +++++++++++-----
 src/aspire/utils/relion_interop.py |  1 +
 tests/test_array_image_source.py   |  4 ++--
 tests/test_simulation.py           |  2 ++
 4 files changed, 16 insertions(+), 7 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 02f4942d74..fa1615c372 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -483,15 +483,21 @@ def offsets(self, values):
 
     @property
     def amplitudes(self):
-        return np.atleast_1d(
-            self.get_metadata(
-                "_rlnAmplitude", default_value=np.array(1.0, dtype=self.dtype)
+        if self.has_metadata("__amplitude"):
+            values = self.get_metadata("__amplitude")
+        else:
+            values = self.get_metadata(
+                "_rlnAmplitude",
+                default_value=np.array(1.0, dtype=np.float64),
             )
-        )
+        return np.atleast_1d(np.asarray(values, dtype=np.float64))
 
     @amplitudes.setter
     def amplitudes(self, values):
-        return self.set_metadata("_rlnAmplitude", np.array(values, dtype=self.dtype))
+        values = np.asarray(values, dtype=np.float64)
+        self.set_metadata("__amplitude", values)
+        # Drop the legacy field if we encountered it while loading a STAR file.
+        self._metadata.pop("_rlnAmplitude", None)
 
     @property
     def angles(self):
diff --git a/src/aspire/utils/relion_interop.py b/src/aspire/utils/relion_interop.py
index c807c7d863..7f9916bdc2 100644
--- a/src/aspire/utils/relion_interop.py
+++ b/src/aspire/utils/relion_interop.py
@@ -12,6 +12,7 @@
 # of certain key fields used in the codebase,
 # which are originally read from Relion STAR files.
 relion_metadata_fields = {
+    "__amplitude": float,
     "_rlnVoltage": float,
     "_rlnDefocusU": float,
     "_rlnDefocusV": float,
diff --git a/tests/test_array_image_source.py b/tests/test_array_image_source.py
index 8dc2a28fb7..a2c3ee4f0c 100644
--- a/tests/test_array_image_source.py
+++ b/tests/test_array_image_source.py
@@ -323,10 +323,10 @@ def test_dtype_passthrough(dtype):
     # Check dtypes
     np.testing.assert_equal(src.dtype, dtype)
     np.testing.assert_equal(src.images[:].dtype, dtype)
-    np.testing.assert_equal(src.amplitudes.dtype, dtype)
 
-    # offsets are always stored as doubles
+    # offsets and amplitudes are always stored as doubles
     np.testing.assert_equal(src.offsets.dtype, np.float64)
+    np.testing.assert_equal(src.amplitudes.dtype, np.float64)
 
 
 def test_stack_1d_only():
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 944136c999..a066f6b67d 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -882,6 +882,8 @@ def check_metadata(sim_src, relion_src):
     those in a RelionSource.
     """
     for k, v in sim_src._metadata.items():
+        if k.startswith("__"):
+            continue
         try:
             np.testing.assert_array_equal(v, relion_src._metadata[k])
         except AssertionError:

From 1ef064c6deee8c7cd1d3098f062ffdee41a73574 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 31 Oct 2025 11:26:22 -0400
Subject: [PATCH 14/25] use _aspireAmplitude to retain ampitudes for saved
 files

---
 src/aspire/source/image.py         | 15 +++++----------
 src/aspire/source/relion.py        |  4 ++++
 src/aspire/utils/relion_interop.py |  3 ++-
 tests/test_simulation.py           |  2 --
 4 files changed, 11 insertions(+), 13 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index fa1615c372..280528367d 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -483,21 +483,16 @@ def offsets(self, values):
 
     @property
     def amplitudes(self):
-        if self.has_metadata("__amplitude"):
-            values = self.get_metadata("__amplitude")
-        else:
-            values = self.get_metadata(
-                "_rlnAmplitude",
-                default_value=np.array(1.0, dtype=np.float64),
-            )
+        values = self.get_metadata(
+            "_aspireAmplitude",
+            default_value=np.array(1.0, dtype=np.float64),
+        )
         return np.atleast_1d(np.asarray(values, dtype=np.float64))
 
     @amplitudes.setter
     def amplitudes(self, values):
         values = np.asarray(values, dtype=np.float64)
-        self.set_metadata("__amplitude", values)
-        # Drop the legacy field if we encountered it while loading a STAR file.
-        self._metadata.pop("_rlnAmplitude", None)
+        self.set_metadata("_aspireAmplitude", values)
 
     @property
     def angles(self):
diff --git a/src/aspire/source/relion.py b/src/aspire/source/relion.py
index c3620c9bdb..f9851eb69f 100644
--- a/src/aspire/source/relion.py
+++ b/src/aspire/source/relion.py
@@ -195,6 +195,10 @@ def populate_metadata(self):
 
         metadata = RelionStarFile(self.filepath).get_merged_data_block()
 
+        # Promote legacy _rlnAmplitude column to the ASPIRE-specific name.
+        if "_rlnAmplitude" in metadata and "_aspireAmplitude" not in metadata:
+            metadata["_aspireAmplitude"] = metadata.pop("_rlnAmplitude")
+
         # particle locations are stored as e.g. '000001@first_micrograph.mrcs'
         # in the _rlnImageName column. here, we're splitting this information
         # so we can get the particle's index in the .mrcs stack as an int
diff --git a/src/aspire/utils/relion_interop.py b/src/aspire/utils/relion_interop.py
index 7f9916bdc2..f7774689df 100644
--- a/src/aspire/utils/relion_interop.py
+++ b/src/aspire/utils/relion_interop.py
@@ -12,7 +12,8 @@
 # of certain key fields used in the codebase,
 # which are originally read from Relion STAR files.
 relion_metadata_fields = {
-    "__amplitude": float,
+    "_aspireAmplitude": float,
+    "_rlnAmplitude": float,
     "_rlnVoltage": float,
     "_rlnDefocusU": float,
     "_rlnDefocusV": float,
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index a066f6b67d..944136c999 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -882,8 +882,6 @@ def check_metadata(sim_src, relion_src):
     those in a RelionSource.
     """
     for k, v in sim_src._metadata.items():
-        if k.startswith("__"):
-            continue
         try:
             np.testing.assert_array_equal(v, relion_src._metadata[k])
         except AssertionError:

From 020e38cab82e55a3937666b0b8e0d4835a19f4ea Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 3 Nov 2025 10:29:34 -0500
Subject: [PATCH 15/25] ensure pixel size is added to mrc header

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

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 280528367d..afdbdbce16 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1499,6 +1499,8 @@ def save_images(
                 #   for large arrays.
                 stats.update_header(mrc)
 
+                # Add pixel size to header
+                mrc.voxel_size = self.pixel_size
         else:
             # save all images into multiple mrc files in batch size
             for i_start in np.arange(0, self.n, batch_size):
@@ -1512,6 +1514,7 @@ def save_images(
                     f"Saving ImageSource[{i_start}-{i_end-1}] to {mrcs_filepath}"
                 )
                 im = self.images[i_start:i_end]
+                im.pixel_size = self.pixel_size
                 im.save(mrcs_filepath, overwrite=overwrite)
 
     def estimate_signal_mean_energy(

From 65fbfcd57ad8ccd30b18ccd8888909aef18d4f93 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Mon, 3 Nov 2025 11:20:15 -0500
Subject: [PATCH 16/25] test mrc pixel_size in header

---
 tests/test_simulation.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 944136c999..f326612aa1 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -3,6 +3,7 @@
 import tempfile
 from unittest import TestCase
 
+import mrcfile
 import numpy as np
 import pytest
 
@@ -876,6 +877,25 @@ def test_save_overwrite(caplog):
         check_metadata(sim2, sim2_loaded_renamed)
 
 
+@pytest.mark.parametrize("batch_size", [1, 6])
+def test_simulation_save_sets_voxel_size(tmp_path, batch_size):
+    """
+    Test we save with pixel_size appended to the mrcfile header.
+    """
+    # Note, n=6 and batch_size=6 exercises save_mode=='single' branch.
+    sim = Simulation(n=6, L=24, pixel_size=1.37)
+    info = sim.save(tmp_path / "pixel_size.star", batch_size=batch_size, overwrite=True)
+
+    for stack_name in info["mrcs"]:
+        stack_path = tmp_path / stack_name
+        with mrcfile.open(stack_path, permissive=True) as f:
+            vs = f.voxel_size
+            header_vals = np.array(
+                [float(vs.x), float(vs.y), float(vs.z)], dtype=np.float64
+            )
+            np.testing.assert_allclose(header_vals, sim.pixel_size)
+
+
 def check_metadata(sim_src, relion_src):
     """
     Helper function to test if metadata fields in a Simulation match

From 5c288906c468699b13d37559d1e1b3a3865938d7 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Wed, 5 Nov 2025 10:11:15 -0500
Subject: [PATCH 17/25] Use defaultdict

---
 src/aspire/source/image.py | 8 ++------
 tests/test_simulation.py   | 2 +-
 2 files changed, 3 insertions(+), 7 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index afdbdbce16..d7af691abc 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -3,7 +3,7 @@
 import logging
 import os.path
 from abc import ABC, abstractmethod
-from collections import OrderedDict
+from collections import OrderedDict, defaultdict
 from collections.abc import Iterable
 
 import mrcfile
@@ -1350,11 +1350,7 @@ def _prepare_relion_optics_blocks(metadata):
         metadata["_rlnOpticsGroup"] = optics_groups
 
         # Build the optics block rows and assign group names.
-        optics_block = OrderedDict()
-        optics_block["_rlnOpticsGroup"] = []
-        optics_block["_rlnOpticsGroupName"] = []
-        for field in optics_value_fields:
-            optics_block[field] = []
+        optics_block = defaultdict(list)
 
         for signature, group_id in group_lookup.items():
             optics_block["_rlnOpticsGroup"].append(group_id)
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index f326612aa1..83c0932c17 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -663,7 +663,7 @@ def test_simulation_save_optics_block(tmp_path):
     for field in expected_optics_fields:
         assert field in optics
 
-    # Optics group and group name should 1-indexed
+    # Optics group and group name should be 1-indexed
     np.testing.assert_array_equal(
         optics["_rlnOpticsGroup"], np.arange(1, kv_ct + 1, dtype=int)
     )

From 8ab9ee7d9296d7299eba27274c97faba9448b4c1 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Wed, 5 Nov 2025 10:35:14 -0500
Subject: [PATCH 18/25] Remove old file compatibility

---
 src/aspire/source/relion.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/aspire/source/relion.py b/src/aspire/source/relion.py
index f9851eb69f..c3620c9bdb 100644
--- a/src/aspire/source/relion.py
+++ b/src/aspire/source/relion.py
@@ -195,10 +195,6 @@ def populate_metadata(self):
 
         metadata = RelionStarFile(self.filepath).get_merged_data_block()
 
-        # Promote legacy _rlnAmplitude column to the ASPIRE-specific name.
-        if "_rlnAmplitude" in metadata and "_aspireAmplitude" not in metadata:
-            metadata["_aspireAmplitude"] = metadata.pop("_rlnAmplitude")
-
         # particle locations are stored as e.g. '000001@first_micrograph.mrcs'
         # in the _rlnImageName column. here, we're splitting this information
         # so we can get the particle's index in the .mrcs stack as an int

From 60942fa4f2eece03ebe4895049d0c2d2b40411e8 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Wed, 5 Nov 2025 15:44:58 -0500
Subject: [PATCH 19/25] Remove _aspireAmplitude from relion_metadata_fields
 dict.

---
 src/aspire/utils/relion_interop.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/aspire/utils/relion_interop.py b/src/aspire/utils/relion_interop.py
index f7774689df..c807c7d863 100644
--- a/src/aspire/utils/relion_interop.py
+++ b/src/aspire/utils/relion_interop.py
@@ -12,8 +12,6 @@
 # of certain key fields used in the codebase,
 # which are originally read from Relion STAR files.
 relion_metadata_fields = {
-    "_aspireAmplitude": float,
-    "_rlnAmplitude": float,
     "_rlnVoltage": float,
     "_rlnDefocusU": float,
     "_rlnDefocusV": float,

From 7091078242ca5a9f71a0797f22ea3e79b2fe0884 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 7 Nov 2025 11:49:23 -0500
Subject: [PATCH 20/25] Test save on source slices

---
 tests/test_simulation.py | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 83c0932c17..191fca0810 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -695,6 +695,44 @@ def test_simulation_save_optics_block(tmp_path):
     )
 
 
+def test_simulation_slice_save_roundtrip(tmp_path):
+    # Radial CTF Filters
+    kv_min, kv_max, kv_ct = 200, 300, 3
+    voltages = np.linspace(kv_min, kv_max, kv_ct)
+    ctf_filters = [RadialCTFFilter(voltage=kv) for kv in voltages]
+
+    # Generate and save slice of Simulation
+    sim = Simulation(n=9, L=16, C=1, unique_filters=ctf_filters, pixel_size=1.34)
+    sliced_sim = sim[::2]
+    save_path = tmp_path / "sliced_sim.star"
+    sliced_sim.save(save_path, overwrite=True)
+
+    # Load saved slice and compare to original
+    reloaded = RelionSource(save_path)
+
+    # Check images
+    np.testing.assert_allclose(
+        reloaded.images[:].asnumpy(),
+        sliced_sim.images[:].asnumpy(),
+    )
+
+    # Check metadata related to optics block
+    metadata_fields = [
+        "_rlnVoltage",
+        "_rlnDefocusU",
+        "_rlnDefocusV",
+        "_rlnDefocusAngle",
+        "_rlnSphericalAberration",
+        "_rlnAmplitudeContrast",
+        "_rlnImagePixelSize",
+    ]
+    for field in metadata_fields:
+        np.testing.assert_allclose(
+            reloaded.get_metadata(field),
+            sliced_sim.get_metadata(field),
+        )
+
+
 def test_default_symmetry_group():
     # Check that default is "C1".
     sim = Simulation()

From 195f44c7f73deef27cfab259acf6076f16ff066d Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Thu, 13 Nov 2025 15:06:36 -0500
Subject: [PATCH 21/25] Always write an optics block

---
 src/aspire/source/image.py      | 31 +++++++++++++------------------
 src/aspire/source/relion.py     |  3 ++-
 tests/test_coordinate_source.py | 17 +++++++++++++++--
 tests/test_relion_source.py     | 19 ++++++++++++++++---
 4 files changed, 46 insertions(+), 24 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index d7af691abc..3b0db192cd 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1311,31 +1311,26 @@ def _prepare_relion_optics_blocks(metadata):
             "_rlnImageDimensionality",
         ]
 
-        # We only write a data_optics block when every required field is present.
-        required_fields = [
-            "_rlnSphericalAberration",
-            "_rlnVoltage",
-            "_rlnAmplitudeContrast",
-            "_rlnImageSize",
-            "_rlnImageDimensionality",
-        ]
-        pixel_fields = ["_rlnImagePixelSize", "_rlnMicrographPixelSize"]
+        # Some optics group fields might not always be present, but are necessary
+        # for reading the file in Relion. We ensure these fields exist and populate
+        # with a dummy value if not.
+        n_rows = len(metadata["_rlnImageName"])
 
-        has_required = all(field in metadata for field in required_fields)
-        has_pixel_field = any(field in metadata for field in pixel_fields)
+        def _ensure_column(field, value):
+            if field not in metadata:
+                logger.warning(
+                    f"Optics field {field} not found, populating with default value {value}"
+                )
+                metadata[field] = np.full(n_rows, value)
 
-        if not (has_required and has_pixel_field):
-            # Optics metadata incomplete, fall back to legacy single block.
-            logger.warning(
-                "Optics metadata incomplete, writing only data_particles block."
-            )
-            return None, metadata
+        _ensure_column("_rlnSphericalAberration", 2.0)
+        _ensure_column("_rlnVoltage", 300)
+        _ensure_column("_rlnAmplitudeContrast", 0.1)
 
         # Restrict to the optics columns that are actually present on this source.
         optics_value_fields = [
             field for field in all_optics_fields if field in metadata
         ]
-        n_rows = len(metadata["_rlnImageName"])
 
         # Map each unique optics tuple to a 1-based group ID in order encountered.
         group_lookup = OrderedDict()
diff --git a/src/aspire/source/relion.py b/src/aspire/source/relion.py
index c3620c9bdb..7d7cc3cd54 100644
--- a/src/aspire/source/relion.py
+++ b/src/aspire/source/relion.py
@@ -168,7 +168,8 @@ def __init__(
                 f"Found partially populated CTF Params."
                 f"  To automatically populate CTFFilters provide {CTF_params}"
             )
-
+            self.unique_filters = [IdentityFilter()]
+            self.filter_indices = np.zeros(self.n, dtype=int)
         # If no CTF info in STAR, we initialize the filter values of metadata with default values
         else:
             self.unique_filters = [IdentityFilter()]
diff --git a/tests/test_coordinate_source.py b/tests/test_coordinate_source.py
index 2a63c7b0b9..13219dcaa3 100644
--- a/tests/test_coordinate_source.py
+++ b/tests/test_coordinate_source.py
@@ -526,7 +526,7 @@ def testSave(self):
         # load saved particle stack
         saved_star = StarFile(star_path)
         # we want to read the saved mrcs file from the STAR file
-        image_name_column = saved_star.get_block_by_index(0)["_rlnImageName"]
+        image_name_column = saved_star.get_block_by_index(1)["_rlnImageName"]
         # we're reading a string of the form 0000X@mrcs_path.mrcs
         _particle, mrcs_path = image_name_column[0].split("@")
         saved_mrcs_stack = mrcfile.open(os.path.join(self.data_folder, mrcs_path)).data
@@ -537,15 +537,28 @@ def testSave(self):
         self.assertEqual(
             list(saved_star["particles"].keys()),
             [
-                "_rlnImagePixelSize",
+                "_rlnOpticsGroup",
                 "_rlnSymmetryGroup",
                 "_rlnImageName",
                 "_rlnCoordinateX",
                 "_rlnCoordinateY",
+            ],
+        )
+
+        self.assertEqual(
+            list(saved_star["optics"].keys()),
+            [
+                "_rlnOpticsGroup",
+                "_rlnOpticsGroupName",
+                "_rlnImagePixelSize",
+                "_rlnSphericalAberration",
+                "_rlnVoltage",
+                "_rlnAmplitudeContrast",
                 "_rlnImageSize",
                 "_rlnImageDimensionality",
             ],
         )
+
         # assert that all the correct coordinates were saved
         for i in range(10):
             self.assertEqual(
diff --git a/tests/test_relion_source.py b/tests/test_relion_source.py
index d0b996c795..31638bdb7e 100644
--- a/tests/test_relion_source.py
+++ b/tests/test_relion_source.py
@@ -70,6 +70,7 @@ def test_prepare_relion_optics_blocks_warns(caplog):
         "_rlnImagePixelSize": np.array([1.234]),
         "_rlnImageSize": np.array([32]),
         "_rlnImageDimensionality": np.array([2]),
+        "_rlnImageName": np.array(["000001@stack.mrcs"]),
     }
 
     caplog.clear()
@@ -78,9 +79,21 @@ def test_prepare_relion_optics_blocks_warns(caplog):
             metadata.copy()
         )
 
-    assert optics_block is None
-    assert particle_block == metadata
-    assert "Optics metadata incomplete" in caplog.text
+    # We should get and optics block
+    assert optics_block is not None
+
+    # Verify defaults were injected.
+    np.testing.assert_allclose(optics_block["_rlnImagePixelSize"], [1.234])
+    np.testing.assert_array_equal(optics_block["_rlnImageSize"], [32])
+    np.testing.assert_array_equal(optics_block["_rlnImageDimensionality"], [2])
+    np.testing.assert_allclose(optics_block["_rlnVoltage"], [300.0])
+    np.testing.assert_allclose(optics_block["_rlnSphericalAberration"], [2.0])
+    np.testing.assert_allclose(optics_block["_rlnAmplitudeContrast"], [0.1])
+
+    # Caplog should contain the warnings about the three missing fields.
+    assert "Optics field _rlnSphericalAberration not found" in caplog.text
+    assert "Optics field _rlnVoltage not found" in caplog.text
+    assert "Optics field _rlnAmplitudeContrast not found" in caplog.text
 
 
 def test_pixel_size(caplog):

From a3d6381f6fa431dca1537b5a5aa04bc5f16ec689 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 14 Nov 2025 08:10:13 -0500
Subject: [PATCH 22/25] Gallery example simulation -> relion

---
 .../save_simulation_relion_reconstruct.py     | 87 +++++++++++++++++++
 1 file changed, 87 insertions(+)
 create mode 100644 gallery/experiments/save_simulation_relion_reconstruct.py

diff --git a/gallery/experiments/save_simulation_relion_reconstruct.py b/gallery/experiments/save_simulation_relion_reconstruct.py
new file mode 100644
index 0000000000..babeffc351
--- /dev/null
+++ b/gallery/experiments/save_simulation_relion_reconstruct.py
@@ -0,0 +1,87 @@
+"""
+Simulated Stack → RELION Reconstruction
+=======================================
+
+This experiment shows how to:
+
+1. build a synthetic dataset with ASPIRE,
+2. write the stack via ``ImageSource.save`` so RELION can consume it, and
+3. call :code:`relion_reconstruct` on the saved STAR file.
+"""
+
+# %%
+# Imports
+# -------
+
+import logging
+from pathlib import Path
+
+import numpy as np
+
+from aspire.downloader import emdb_2660
+from aspire.noise import WhiteNoiseAdder
+from aspire.operators import RadialCTFFilter
+from aspire.source import RelionSource, Simulation
+
+logger = logging.getLogger(__name__)
+
+
+# %%
+# Configuration
+# -------------
+# We set a few parameters to initialize the Simulation.
+# You can safely alter ``n_particles`` (or change the voltages, etc.) when
+# trying this interactively; the defaults here are chosen for demonstrative purposes.
+
+output_dir = Path("relion_save_demo")
+output_dir.mkdir(exist_ok=True)
+
+n_particles = 512
+snr = 0.25
+voltages = np.linspace(200, 300, 3)  # kV settings for the radial CTF filters
+star_path = output_dir / f"sim_n{n_particles}.star"
+
+
+# %%
+# Volume and Filters
+# ------------------
+# Start from the EMDB-2660 ribosome map and build a small set of radial CTF filters
+# that RELION will recover as optics groups.
+
+vol = emdb_2660()
+ctf_filters = [RadialCTFFilter(voltage=kv) for kv in voltages]
+
+
+# %%
+# Simulate, Add Noise, Save
+# -------------------------
+# Initialize the Simulation:
+# mix the CTFs across the stack, add white noise at a target SNR,
+# and write the particles and metadata to a RELION-compatible STAR/MRC stack.
+
+sim = Simulation(
+    n=n_particles,
+    vols=vol,
+    unique_filters=ctf_filters,
+    noise_adder=WhiteNoiseAdder.from_snr(snr),
+)
+sim.save(star_path, overwrite=True)
+
+
+# %%
+# Running ``relion_reconstruct``
+# ------------------------------
+# ``relion_reconstruct`` is an external RELION command, so we just show the call.
+# Run this in a RELION-enabled shell after generating the STAR file above.
+
+relion_cmd = [
+    "relion_reconstruct",
+    "--i",
+    str(star_path),
+    "--o",
+    str(output_dir / "relion_recon.mrc"),
+    "--ctf",
+]
+
+print(" ".join(relion_cmd))
+

From 7324533b0065f0249744088bb80db2939c90e302 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 14 Nov 2025 08:25:15 -0500
Subject: [PATCH 23/25] clean up gallery

---
 gallery/experiments/save_simulation_relion_reconstruct.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/gallery/experiments/save_simulation_relion_reconstruct.py b/gallery/experiments/save_simulation_relion_reconstruct.py
index babeffc351..767abf6bfc 100644
--- a/gallery/experiments/save_simulation_relion_reconstruct.py
+++ b/gallery/experiments/save_simulation_relion_reconstruct.py
@@ -21,7 +21,7 @@
 from aspire.downloader import emdb_2660
 from aspire.noise import WhiteNoiseAdder
 from aspire.operators import RadialCTFFilter
-from aspire.source import RelionSource, Simulation
+from aspire.source import Simulation
 
 logger = logging.getLogger(__name__)
 
@@ -83,5 +83,4 @@
     "--ctf",
 ]
 
-print(" ".join(relion_cmd))
-
+logger.info(" ".join(relion_cmd))

From 665e038cab5ab319c57ec471f14d056b580feb82 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 14 Nov 2025 15:14:26 -0500
Subject: [PATCH 24/25] _aspireMetadata 'no_ctf' tag for missing optics fields.
 Detect ASPIRE-generated dummy values when loading. Unit test.

---
 src/aspire/source/image.py  | 12 +++++++++---
 src/aspire/source/relion.py | 17 +++++++++++++++--
 tests/test_simulation.py    | 28 ++++++++++++++++++++++++++++
 3 files changed, 52 insertions(+), 5 deletions(-)

diff --git a/src/aspire/source/image.py b/src/aspire/source/image.py
index 3b0db192cd..aa1b85a1b3 100644
--- a/src/aspire/source/image.py
+++ b/src/aspire/source/image.py
@@ -1316,16 +1316,22 @@ def _prepare_relion_optics_blocks(metadata):
         # with a dummy value if not.
         n_rows = len(metadata["_rlnImageName"])
 
+        missing_fields = []
+
         def _ensure_column(field, value):
             if field not in metadata:
+                missing_fields.append(field)
                 logger.warning(
                     f"Optics field {field} not found, populating with default value {value}"
                 )
                 metadata[field] = np.full(n_rows, value)
 
-        _ensure_column("_rlnSphericalAberration", 2.0)
-        _ensure_column("_rlnVoltage", 300)
-        _ensure_column("_rlnAmplitudeContrast", 0.1)
+        _ensure_column("_rlnSphericalAberration", 0)
+        _ensure_column("_rlnVoltage", 0)
+        _ensure_column("_rlnAmplitudeContrast", 0)
+
+        if missing_fields:
+            metadata["_aspireMetadata"] = np.full(n_rows, "no_ctf", dtype=object)
 
         # Restrict to the optics columns that are actually present on this source.
         optics_value_fields = [
diff --git a/src/aspire/source/relion.py b/src/aspire/source/relion.py
index 7d7cc3cd54..c063c3b7b2 100644
--- a/src/aspire/source/relion.py
+++ b/src/aspire/source/relion.py
@@ -125,6 +125,12 @@ def __init__(
             for key in offset_keys:
                 del self._metadata[key]
 
+        # Detect ASPIRE-generated dummy variables
+        aspire_metadata = metadata.get("_aspireMetadata")
+        dummy_ctf = isinstance(aspire_metadata, (list, np.ndarray)) and np.all(
+            np.asarray(aspire_metadata) == "no_ctf"
+        )
+
         # CTF estimation parameters coming from Relion
         CTF_params = [
             "_rlnVoltage",
@@ -162,14 +168,21 @@ def __init__(
             # self.unique_filters of the filter that should be applied
             self.filter_indices = filter_indices
 
+        # If we detect ASPIRE added dummy variables, log and initialize identity filter
+        elif dummy_ctf:
+            logger.info(
+                "Detected ASPIRE-generated dummy optics; initializing identity filters."
+            )
+            self.unique_filters = [IdentityFilter()]
+            self.filter_indices = np.zeros(self.n, dtype=int)
+
         # We have provided some, but not all the required params
         elif any(param in metadata for param in CTF_params):
             logger.warning(
                 f"Found partially populated CTF Params."
                 f"  To automatically populate CTFFilters provide {CTF_params}"
             )
-            self.unique_filters = [IdentityFilter()]
-            self.filter_indices = np.zeros(self.n, dtype=int)
+
         # If no CTF info in STAR, we initialize the filter values of metadata with default values
         else:
             self.unique_filters = [IdentityFilter()]
diff --git a/tests/test_simulation.py b/tests/test_simulation.py
index 191fca0810..c053b2eb0f 100644
--- a/tests/test_simulation.py
+++ b/tests/test_simulation.py
@@ -915,6 +915,34 @@ def test_save_overwrite(caplog):
         check_metadata(sim2, sim2_loaded_renamed)
 
 
+def test_save_load_dummy_ctf_values(tmp_path, caplog):
+    """
+    Test we populate optics group field with dummy values when none
+    are present. These values should be detected upon reloading the source.
+    """
+    star_path = tmp_path / "no_ctf.star"
+    sim = Simulation(n=8, L=16)  # no unique_filters, ie. no CTF info
+    sim.save(star_path, overwrite=True)
+
+    # STAR file should contain our fallback tag
+    star = RelionStarFile(star_path)
+    particles_block = star.get_block_by_index(1)
+    np.testing.assert_array_equal(
+        particles_block["_aspireMetadata"], np.full(sim.n, "no_ctf", dtype=object)
+    )
+
+    # Tag should survive round-trip
+    caplog.clear()
+    reloaded = RelionSource(star_path)
+    np.testing.assert_array_equal(
+        reloaded._metadata["_aspireMetadata"],
+        np.full(reloaded.n, "no_ctf", dtype=object),
+    )
+
+    # Check message is logged about detecting dummy variables
+    assert "Detected ASPIRE-generated dummy optics" in caplog.text
+
+
 @pytest.mark.parametrize("batch_size", [1, 6])
 def test_simulation_save_sets_voxel_size(tmp_path, batch_size):
     """

From 5cb9557d1d248c29634e482d25c6ed8a58b4ece1 Mon Sep 17 00:00:00 2001
From: Josh Carmichael <carmichael@princeton.edu>
Date: Fri, 14 Nov 2025 15:56:50 -0500
Subject: [PATCH 25/25] test cleanup

---
 tests/test_coordinate_source.py | 1 +
 tests/test_relion_source.py     | 6 +++---
 2 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/tests/test_coordinate_source.py b/tests/test_coordinate_source.py
index 13219dcaa3..8baf7d1338 100644
--- a/tests/test_coordinate_source.py
+++ b/tests/test_coordinate_source.py
@@ -542,6 +542,7 @@ def testSave(self):
                 "_rlnImageName",
                 "_rlnCoordinateX",
                 "_rlnCoordinateY",
+                "_aspireMetadata",
             ],
         )
 
diff --git a/tests/test_relion_source.py b/tests/test_relion_source.py
index 31638bdb7e..64703ed491 100644
--- a/tests/test_relion_source.py
+++ b/tests/test_relion_source.py
@@ -86,9 +86,9 @@ def test_prepare_relion_optics_blocks_warns(caplog):
     np.testing.assert_allclose(optics_block["_rlnImagePixelSize"], [1.234])
     np.testing.assert_array_equal(optics_block["_rlnImageSize"], [32])
     np.testing.assert_array_equal(optics_block["_rlnImageDimensionality"], [2])
-    np.testing.assert_allclose(optics_block["_rlnVoltage"], [300.0])
-    np.testing.assert_allclose(optics_block["_rlnSphericalAberration"], [2.0])
-    np.testing.assert_allclose(optics_block["_rlnAmplitudeContrast"], [0.1])
+    np.testing.assert_allclose(optics_block["_rlnVoltage"], [0])
+    np.testing.assert_allclose(optics_block["_rlnSphericalAberration"], [0])
+    np.testing.assert_allclose(optics_block["_rlnAmplitudeContrast"], [0])
 
     # Caplog should contain the warnings about the three missing fields.
     assert "Optics field _rlnSphericalAberration not found" in caplog.text