From d70399a290b592f1cd9ff2471ce67b4f9c72290b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 20 Sep 2026 16:18:15 +0200 Subject: [PATCH] Fit the rotation step to the whole sweep, not to one recorded pair A 2700 frame, 540 degree sweep was refused at load as "screening images taken at scattered angles". Nothing about it is irregular: every consecutive difference is a clean step, and the headers simply write the start angle modulo 360, so the second and third revolutions repeat angles the first already had. The unwrap handled that correctly. What failed was the STEP. It was taken to be the smallest difference between two adjacent frames - a single recorded value, and of all the samples the one biased furthest low by read-back noise. These headers are written from 32-bit floats, so that one sample was 1.3e-4 low, and placing a frame by dividing its distance from the FIRST frame by that step multiplied the error by the frame number: by frame 1924 the drift had reached a quarter of a step and the on-grid check refused the rest of the sweep. The fix separates the two jobs the step was doing. The smallest difference is now only a guess, used to count how many steps apart each pair of NEIGHBOURS is; a local difference carries the noise of one reading and nothing accumulates, however long the sweep. The step itself is then a straight line fitted through (steps, angle) over every frame, so the noise averages out instead of one unlucky pair setting the scale for thousands of frames. Over the whole corpus this changes the step for 4 of 82 numbered series and leaves the other 78 bit-identical; of the 4, three were drifting by 0.17 to 0.35 of a step across the sweep and one by 1e-8. Counting steps between neighbours is also what carries a sweep past a full turn: each fold is one step forward, so the count keeps climbing through 360 and a frame taken on the second revolution lands beyond the first rather than on top of it. A series that genuinely goes BACK - two sweeps of one crystal concatenated - folds to a large negative count in a single pair and still collides, so that refusal keeps working; its message now says what it means. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- reader/SweepLayout.cpp | 80 ++++++++++++++++++++++++++++++--------- tests/SweepLayoutTest.cpp | 59 +++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 18 deletions(-) diff --git a/reader/SweepLayout.cpp b/reader/SweepLayout.cpp index 8127f7741..981fa0137 100644 --- a/reader/SweepLayout.cpp +++ b/reader/SweepLayout.cpp @@ -117,29 +117,33 @@ Layout Place(const std::vector &frames, const std::string &logger_name) { for (size_t i = 1; i < frames.size(); i++) angle[i] = angle[i - 1] + Fold(frames[i].angle_deg - angle[i - 1]); - // ---- one rotation step + // ---- a first guess at the rotation step // - // The step is the smallest move between two frames that ARE adjacent in the series, signed with - // the way it went. Taking the FIRST pair instead - which is what this code used to do - reads a - // gap as the step and compresses the whole sweep by however much is missing. The header's own + // The smallest move between two frames that ARE adjacent in the series, signed with the way it + // went. Taking the FIRST pair instead - which is what this code used to do - reads a gap as the + // step and compresses the whole sweep by however much is missing. The header's own // Angle_increment is not used for this: it is the oscillation WIDTH, which a series with // overlapping or spaced wedges does not step by. // // A difference smaller than half the oscillation width is not a step but jitter in the recorded // angle: no instrument slices finer than it exposes, so wedges overlapping twofold would be a // read-back wobble, and taking one as the step would spread the sweep over millions of slots. + // + // It is only a guess, and deliberately the crudest one: a single recorded difference, and the + // one biased furthest low by the read-back noise. It is used below to count STEPS between + // neighbours, never to place a frame outright. const double too_fine = 0.5 * std::abs(f0.increment_deg); - double step = 0; + double guess = 0; for (size_t i = 1; i < frames.size(); i++) { const double d = Fold(angle[i] - angle[i - 1]); - if (std::abs(d) > std::max(too_fine, 1e-6) && (step == 0 || std::abs(d) < std::abs(step))) - step = d; + if (std::abs(d) > std::max(too_fine, 1e-6) && (guess == 0 || std::abs(d) < std::abs(guess))) + guess = d; } Layout out; // A series that never turns: a grid scan, a set of stills, or a single image. There is no sweep // to place anything on, so the files are the slots and the header's nominal increment stands. - if (step == 0) { + if (guess == 0) { out.files.reserve(frames.size()); for (const auto &f : frames) out.files.push_back(f.path); @@ -149,15 +153,51 @@ Layout Place(const std::vector &frames, const std::string &logger_name) { return out; } - // ---- every frame on that step, or this is not a rotation series + // ---- how many steps apart each pair of NEIGHBOURS is + // + // Counted from one local difference at a time, so nothing accumulates: a recorded angle carries + // a few parts in 100000 of read-back noise (these headers are written from 32-bit floats), which + // is nowhere near half a step for one pair however long the sweep is. Measuring the frame's + // position against the FIRST frame instead multiplies that noise by the frame number - on a 2700 + // frame sweep the guess above was 1.3e-4 low and by frame 1924 the drift had reached a quarter + // of a step, so a perfectly regular series was refused as scattered. + // + // This is also what carries a sweep past a full turn. Each fold is one step forward, so the + // count keeps climbing through 360 and a frame taken on the second revolution lands beyond the + // first, not on top of it. A series that genuinely jumps BACK - two sweeps of one crystal + // concatenated - folds to a large negative count in a single pair and collides below. std::vector slot(frames.size()); - std::vector off_grid; + slot[0] = 0; + for (size_t i = 1; i < frames.size(); i++) + slot[i] = slot[i - 1] + std::llround(Fold(angle[i] - angle[i - 1]) / guess); + + // ---- the step the whole series agrees on + // + // A straight line through (steps, angle): every frame's own header has a say, so the read-back + // noise averages out instead of one unlucky pair setting the scale for thousands of frames. + double mean_slot = 0, mean_angle = 0; for (size_t i = 0; i < frames.size(); i++) { - const double k = (angle[i] - angle[0]) / step; - slot[i] = std::llround(k); - if (std::abs(k - static_cast(slot[i])) > 0.25) - off_grid.push_back(frames[i].path); + mean_slot += static_cast(slot[i]); + mean_angle += angle[i]; } + mean_slot /= static_cast(frames.size()); + mean_angle /= static_cast(frames.size()); + double cov = 0, var = 0; + for (size_t i = 0; i < frames.size(); i++) { + const double ds = static_cast(slot[i]) - mean_slot; + cov += ds * (angle[i] - mean_angle); + var += ds * ds; + } + // Two frames a step apart give var > 0; var == 0 only if every frame landed on one step, which + // the guess above has already ruled out. + const double step = var > 0 ? cov / var : guess; + const double start = mean_angle - step * mean_slot; + + // ---- every frame on that step, or this is not a rotation series + std::vector off_grid; + for (size_t i = 0; i < frames.size(); i++) + if (std::abs(angle[i] - (start + step * static_cast(slot[i]))) > 0.25 * std::abs(step)) + off_grid.push_back(frames[i].path); if (!off_grid.empty()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, fmt::format("The images named here do not lie on one rotation series: " @@ -180,14 +220,18 @@ Layout Place(const std::vector &frames, const std::string &logger_name) { duplicates.push_back(frames[i].path); at = frames[i].path; } + // A sweep of more than one full turn reaches the same recorded angle again and is NOT this: it + // gets there by stepping forward the whole way, so its second revolution occupies fresh slots. + // What collides here is a series that jumped back - two sweeps of one crystal concatenated. if (!duplicates.empty()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, - fmt::format("The images named here repeat an angle already taken by " - "another image of the series: {}. Two sweeps of the same " - "crystal have to be processed separately.", + fmt::format("The images named here go back to a point of the rotation " + "an earlier image of the series already took: {}. Two " + "sweeps of the same crystal have to be processed " + "separately.", NameSome(duplicates))); - out.start_deg = angle[0] + step * static_cast(first); + out.start_deg = start + step * static_cast(first); out.increment_deg = step; out.present = frames.size(); diff --git a/tests/SweepLayoutTest.cpp b/tests/SweepLayoutTest.cpp index 9655afe01..ee472c4d2 100644 --- a/tests/SweepLayoutTest.cpp +++ b/tests/SweepLayoutTest.cpp @@ -3,6 +3,8 @@ #include +#include + #include "../common/JFJochException.h" #include "../reader/SweepLayout.h" @@ -55,6 +57,63 @@ TEST_CASE("SweepLayout_PastFullTurn") { CHECK(l.increment_deg == Catch::Approx(0.1)); } +TEST_CASE("SweepLayout_MultiTurn") { + // One continuous 540 degree sweep, 2700 frames of 0.2 deg, whose headers write the start angle + // modulo 360: from frame 1476 on, every angle repeats one the first revolution already had. The + // second and third revolutions must occupy fresh slots, not land on top of the first. + std::vector angles; + for (int i = 0; i < 2700; i++) + angles.push_back(std::fmod(65.0 + 0.2 * i, 360.0)); + const auto l = sweep::Place(Series(angles, 0.2), "test"); + REQUIRE(l.files.size() == 2700); + CHECK(l.present == 2700); + CHECK(l.start_deg == Catch::Approx(65.0)); + CHECK(l.increment_deg == Catch::Approx(0.2)); + // The sweep really is 540 degrees wide, and no frame was dropped on top of another. + CHECK(l.increment_deg * static_cast(l.files.size() - 1) == Catch::Approx(539.8)); + CHECK(l.files.back() == "f2699.cbf"); +} + +TEST_CASE("SweepLayout_MultiTurnWithReadbackNoise") { + // The same sweep as a real header writes it: the recorded angles come from 32-bit floats, so + // each is a few parts in 100000 off the grid. Estimating the step from ONE recorded difference + // scales that noise by the frame number - the series it was measured on drifted a quarter of a + // step by frame 1924 and was refused as scattered. + std::vector angles; + for (int i = 0; i < 2700; i++) + angles.push_back(std::fmod(static_cast(static_cast(65.0 + 0.2 * i)), 360.0)); + const auto l = sweep::Place(Series(angles, 0.2), "test"); + REQUIRE(l.files.size() == 2700); + CHECK(l.present == 2700); + CHECK(l.increment_deg == Catch::Approx(0.2).epsilon(1e-6)); + CHECK(l.start_deg == Catch::Approx(65.0).epsilon(1e-6)); +} + +TEST_CASE("SweepLayout_MultiTurnGapped") { + // Past 360 AND missing frames: the two have to work together, because a gap is what the step + // guess has to survive and a full turn is what the step count has to carry through. + std::vector angles; + for (int i = 0; i < 2000; i++) + if (i % 7 != 3) + angles.push_back(std::fmod(10.0 + 0.25 * i, 360.0)); + const auto l = sweep::Place(Series(angles, 0.25), "test"); + CHECK(l.files.size() == 2000); + CHECK(l.present == angles.size()); + CHECK(l.increment_deg == Catch::Approx(0.25)); + CHECK(l.files[3].empty()); + CHECK(l.files[1999] == "f" + std::to_string(angles.size() - 1) + ".cbf"); +} + +TEST_CASE("SweepLayout_ConcatenatedSweepsStillRefused") { + // Two 90 degree sweeps of one crystal in one directory. Unlike a sweep that runs past 360, this + // one goes BACK: the second series returns to angles the first already covered without having + // stepped forward through a whole turn to get there. + std::vector angles; + for (int i = 0; i < 90; i++) angles.push_back(0.5 * i); + for (int i = 0; i < 90; i++) angles.push_back(0.5 * i); + CHECK_THROWS_AS(sweep::Place(Series(angles, 0.5), "test"), JFJochException); +} + TEST_CASE("SweepLayout_Reversed") { const auto l = sweep::Place(Series({10.0, 9.5, 9.0, 8.0}, 0.5), "test"); REQUIRE(l.files.size() == 5);