Files
Jungfraujoch/tests/SweepLayoutTest.cpp
T
leonarski_fandClaude Opus 5 d70399a290 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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
2026-09-20 18:45:18 +02:00

156 lines
6.5 KiB
C++

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include <cmath>
#include "../common/JFJochException.h"
#include "../reader/SweepLayout.h"
namespace {
// A series of frames all taken at the same instrument setting, at the angles given.
std::vector<sweep::Frame> Series(const std::vector<double> &angles, double increment) {
std::vector<sweep::Frame> out;
out.reserve(angles.size());
for (size_t i = 0; i < angles.size(); i++)
out.push_back({"f" + std::to_string(i) + ".cbf", angles[i], increment, 0.2, 1000, 1000, 1.0});
return out;
}
} // namespace
TEST_CASE("SweepLayout_Contiguous") {
std::vector<double> angles;
for (int i = 0; i < 100; i++)
angles.push_back(20.0 + 0.1 * i);
const auto l = sweep::Place(Series(angles, 0.1), "test");
CHECK(l.files.size() == 100);
CHECK(l.present == 100);
CHECK(l.start_deg == Catch::Approx(20.0));
CHECK(l.increment_deg == Catch::Approx(0.1));
for (size_t i = 0; i < l.files.size(); i++)
CHECK(l.files[i] == "f" + std::to_string(i) + ".cbf");
}
TEST_CASE("SweepLayout_Gapped") {
// Frames 3, 4 and 7 of a ten-frame sweep never made it into the archive. The sweep is still ten
// steps wide and every frame keeps its own angle.
const auto l = sweep::Place(Series({0.0, 0.1, 0.2, 0.5, 0.6, 0.8, 0.9}, 0.1), "test");
REQUIRE(l.files.size() == 10);
CHECK(l.present == 7);
CHECK(l.increment_deg == Catch::Approx(0.1));
CHECK(l.files[2] == "f2.cbf");
CHECK(l.files[3].empty());
CHECK(l.files[4].empty());
CHECK(l.files[5] == "f3.cbf");
CHECK(l.files[7].empty());
CHECK(l.files[9] == "f6.cbf");
}
TEST_CASE("SweepLayout_PastFullTurn") {
// A writer that starts over at 0 rather than counting on past 360.
const auto l = sweep::Place(Series({359.7, 359.8, 359.9, 0.0, 0.1}, 0.1), "test");
CHECK(l.files.size() == 5);
CHECK(l.present == 5);
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<double> 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<double>(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<double> angles;
for (int i = 0; i < 2700; i++)
angles.push_back(std::fmod(static_cast<double>(static_cast<float>(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<double> 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<double> 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);
CHECK(l.increment_deg == Catch::Approx(-0.5));
CHECK(l.files[3].empty());
CHECK(l.files[4] == "f3.cbf");
}
TEST_CASE("SweepLayout_Stills") {
// Nothing turns: the files are the images, and the header's nominal increment stands.
const auto l = sweep::Place(Series({45.0, 45.0, 45.0}, 0.1), "test");
CHECK(l.files.size() == 3);
CHECK(l.present == 3);
CHECK(l.start_deg == Catch::Approx(45.0));
CHECK(l.increment_deg == Catch::Approx(0.1));
}
TEST_CASE("SweepLayout_ScreeningImagesRefused") {
// Five shots at scattered angles are not a sweep, and must not be laid out as one.
CHECK_THROWS_AS(sweep::Place(Series({0.0, 90.0, 45.0, 300.0, 270.0}, 0.5), "test"),
JFJochException);
}
TEST_CASE("SweepLayout_MovedDetectorRefused") {
auto frames = Series({0.0, 0.1, 0.2, 0.3}, 0.1);
frames[2].distance_m = 0.3;
CHECK_THROWS_AS(sweep::Place(frames, "test"), JFJochException);
}
TEST_CASE("SweepLayout_SecondWavelengthRefused") {
auto frames = Series({0.0, 0.1, 0.2, 0.3}, 0.1);
frames[3].wavelength_A = 1.9;
CHECK_THROWS_AS(sweep::Place(frames, "test"), JFJochException);
}
TEST_CASE("SweepLayout_RepeatedAngleRefused") {
// Two sweeps of the same crystal concatenated: the second covers angles the first already has.
CHECK_THROWS_AS(sweep::Place(Series({0.0, 0.1, 0.2, 0.1, 0.2}, 0.1), "test"), JFJochException);
}