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);