mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases 73 Wiki Activity
Files
rc166
Jungfraujoch/tests/DetGeomCalibTest.cpp
T
leonarski_fandClaude Opus 5 03456dcddd calibration: an option to hold the detector tilt fixed 
--mode calibration fits five parameters - beam centre, distance and the two
PONI tilts - and a program that cannot express a tilted detector has nowhere to
put the last two. Dropping them after the fact is worse than never fitting
them: the centre and the distance of a tilted fit have already absorbed the
tilt, so the flattened geometry is right nowhere.

rugnux --no-refine-tilt, the "Refine detector tilt" tick box on the viewer's
Calib page and RingOptimizer's refine_tilt argument hold rot1/rot2 at the value
the geometry came in with and fit the remaining three. That is the best
flat-detector answer, and the one such a program would refine to itself.

Measured on a five-distance calibrant series. At short distance the tilt is
real and reproducible - three independent fits agreeing to 0.01 deg, radial rms
1.4 -> 0.4 px - and its direct beam agrees with the background beam-centre
estimator to 0.05 px, so the tilted model is the physically right one. The
pinned fit then displaces the centre 2.6 px to absorb the tilt and lands within
0.03 px of the same place at every distance.

Past ~300 mm, where only two rings reach the detector, the tilt is instead
under-determined: it comes out with the opposite sign to every short-distance
fit and drags the PONI 28 px while the rms does not move (0.960 against 0.962).
The existing degeneracy guard only fires on a strictly single ring, so it does
not catch that; declining a tilt that does not pay for itself in rms is left
for a separate change.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MxrrPcxodNiXzhNiECCVp5
2026-08-30 20:44:12 +02:00

305 lines
10 KiB
C++
Raw Permalink Blame History

// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#include <catch2/catch_all.hpp>
#include "../image_analysis/geom_refinement/AssignSpotsToRings.h"
#include "../common/Definitions.h"
TEST_CASE("DetGeomCalib_FindCircleCenter") {
std::vector<SpotToSave> spots;
// Make two colinear circles
for (int i = 0; i < 7; i++) {
float angle = i * 2 * M_PI / 10;
float x = 200.0f + 100.0f * cosf(angle);
float y = 200.0f + 100.0f * sinf(angle);
spots.push_back(SpotToSave(x, y, 1000));
}
for (int i = 0; i < 7; i++) {
float angle = i * 2 * M_PI / 10;
float x = 200.0f + 50.0f * cosf(angle);
float y = 200.0f + 50.0f * sinf(angle);
spots.push_back(SpotToSave(x, y, 1000));
}
// Add some outliers
spots.push_back(SpotToSave(1000.0f, 1000.0f, 1000));
spots.push_back(SpotToSave(0.0f, 0.0f, 1000));
auto ret = FindCircleCenter(spots);
REQUIRE(ret.x == Catch::Approx(200.0f));
REQUIRE(ret.y == Catch::Approx(200.0f));
}
TEST_CASE("DetGeomCalib_FindCircleCenter_250") {
std::vector<SpotToSave> spots;
for (int i = 0; i < 250; i++) {
float angle = i * 2 * M_PI / 250;
float x = 200.0f + 100.0f * cosf(angle);
float y = 200.0f + 100.0f * sinf(angle);
spots.push_back(SpotToSave(x, y, 1000));
}
auto ret = FindCircleCenter(spots);
REQUIRE(ret.x == Catch::Approx(200.0f));
REQUIRE(ret.y == Catch::Approx(200.0f));
REQUIRE(ret.total_votes == 250 * 249 * 248 / 6); // N * (N-1) * (N-2) / 6
}
TEST_CASE("DetGeomCalib_dbscan") {
std::vector<float> spots_r;
for (int i = 0; i < 15; i++)
spots_r.push_back(100.0);
for (int i = 0; i < 15; i++)
spots_r.push_back(50.0);
for (int i = 0; i < 9; i++)
spots_r.push_back(30.0);
spots_r.push_back(70.0);
spots_r.push_back(11.0);
auto ret = ClusterSpotsIntoRings(spots_r, 0.1, 10);
REQUIRE(ret.size() == 2);
REQUIRE(ret[0][0] == 0);
REQUIRE(ret[0][1] == 1);
REQUIRE(ret[1][0] == 15);
}
TEST_CASE("DetGeomCalib_AnalyzeClusters") {
std::vector<std::vector<int>> dbscan_result;
std::vector<float> spots_r;
std::vector<int> tmp_1;
for (int i = 0; i < 15; i++) {
spots_r.push_back(100.0);
tmp_1.push_back(i);
}
std::vector<int> tmp_2;
for (int i = 0; i < 15; i++) {
spots_r.push_back(50.0);
tmp_2.push_back(15+i);
}
dbscan_result.push_back(tmp_1);
dbscan_result.push_back(tmp_2);
auto ret = AnalyzeClusters(spots_r, dbscan_result);
REQUIRE(ret.size() == 2);
REQUIRE(ret[0].R_obs == Catch::Approx(50.0f));
REQUIRE(ret[1].R_obs == Catch::Approx(100.0f));
}
TEST_CASE("DetGeomCalib_CalculateXtalRings_cubic") {
auto ret = CalculateXtalRings(UnitCell(2.0, 2.0, 2.0, 90, 90, 90));
CHECK(ret[0] == Catch::Approx(2.0 * M_PI * 1.0 / 2.0));
CHECK(ret[1] == Catch::Approx(2.0 * M_PI * sqrt( 2.0 )/ 2.0));
CHECK(ret[2] == Catch::Approx(2.0 * M_PI * sqrt( 3.0 )/ 2.0));
CHECK(ret[3] == Catch::Approx(2.0 * M_PI * 2.0/ 2.0));
// 7 cannot be obtained by h^2 + k^2 + l^2, while 8 can
CHECK(ret[6] == Catch::Approx(2.0 * M_PI * sqrt( 8.0 )/ 2.0));
}
TEST_CASE("DetGeomCalib_CalculateXtalRings_one_long_axis") {
auto ret = CalculateXtalRings(UnitCell(50.0, 2.0, 2.0, 90, 90, 90));
CHECK(ret[0] == Catch::Approx(2.0 * M_PI * 1.0 / 50.0));
CHECK(ret[1] == Catch::Approx(2.0 * M_PI * 2.0 / 50.0));
CHECK(ret[2] == Catch::Approx(2.0 * M_PI * 3.0 / 50.0));
CHECK(ret[3] == Catch::Approx(2.0 * M_PI * 4.0 / 50.0));
CHECK(ret[4] == Catch::Approx(2.0 * M_PI * 5.0 / 50.0));
}
TEST_CASE("DetGeomCalib_GuessDetectorDistance") {
std::vector<SpotToSave> spots;
DiffractionGeometry geom;
geom.Wavelength_A(1.0).BeamX_pxl(100.0).BeamY_pxl(200.0)
.DetectorDistance_mm(1000);
float lab6_a = 4.156468;
float ring_radius_pxl = geom.ResToPxl(lab6_a);
REQUIRE(GuessDetectorDistance(geom, ring_radius_pxl, lab6_a) == Catch::Approx(1000.0));
}
TEST_CASE("DetGeomCalib_GuessInitialGeometry") {
std::vector<SpotToSave> spots;
DiffractionGeometry geom;
geom.Wavelength_A(1.0).BeamX_pxl(100.0).BeamY_pxl(200.0)
.DetectorDistance_mm(1000);
float lab6_a = LAB6_CELL_A;
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a, i * M_PI * 2.0 / 30.0);
spots.push_back(SpotToSave{.x = x, .y = y});
}
DiffractionGeometry geom_out;
geom_out.Wavelength_A(1.0);
GuessInitialGeometry(geom_out, spots, LAB6_CELL_A);
REQUIRE(geom_out.GetBeamX_pxl() == geom.GetBeamX_pxl());
REQUIRE(geom_out.GetBeamY_pxl() == geom.GetBeamY_pxl());
REQUIRE(geom_out.GetDetectorDistance_mm() == Catch::Approx(geom.GetDetectorDistance_mm()));
}
TEST_CASE("DetGeomCalib_GuessGeometry") {
std::vector<SpotToSave> spots;
DiffractionGeometry geom;
geom.Wavelength_A(1.0).BeamX_pxl(100.0).BeamY_pxl(200.0)
.DetectorDistance_mm(100);
float lab6_a = LAB6_CELL_A;
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a, i * M_PI * 2.0 / 30.0);
spots.push_back(SpotToSave{.x = x, .y = y});
}
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a / sqrt(2), i * M_PI * 2.0 / 30.0);
spots.push_back(SpotToSave{.x = x, .y = y});
}
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a / sqrt(3.0f), i * M_PI * 2.0 / 30.0);
spots.push_back(SpotToSave{.x = x, .y = y});
}
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a / 2.0f, i * M_PI * 2.0 / 30.0);
spots.push_back(SpotToSave{.x = x, .y = y});
}
DiffractionGeometry geom_out;
geom_out.Wavelength_A(1.0).DetectorDistance_mm(200.0);
GuessGeometry(geom_out, spots, CalculateXtalRings(UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90)));
CHECK(fabsf(geom_out.GetBeamX_pxl() - geom.GetBeamX_pxl()) < 0.001f);
CHECK(fabsf(geom_out.GetBeamY_pxl() - geom.GetBeamY_pxl()) < 0.001f);
// This is wrong!!!!
CHECK(fabsf(geom_out.GetDetectorDistance_mm() - geom.GetDetectorDistance_mm()) < 0.01f);
}
TEST_CASE("DetGeomCalib_RingOptimizer") {
std::vector<RingOptimizerInput> spots;
DiffractionGeometry geom;
geom.Wavelength_A(1.0).BeamX_pxl(100.0).BeamY_pxl(200.0)
.DetectorDistance_mm(100).PoniRot1_rad(0.1).PoniRot2_rad(0.05);
float lab6_a = LAB6_CELL_A;
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a, i * M_PI * 2.0 / 30.0);
spots.push_back(RingOptimizerInput{
.x = x, .y = y,
.q_expected = M_PI * 2.0 / lab6_a
});
}
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a / sqrt(2), i * M_PI * 2.0 / 30.0);
spots.push_back(RingOptimizerInput{
.x = x, .y = y,
.q_expected = sqrtf(2.0f) * M_PI * 2.0 / lab6_a
});
}
DiffractionGeometry geom_i;
geom_i.Wavelength_A(1.0).BeamX_pxl(105).BeamY_pxl(195).DetectorDistance_mm(110);
RingOptimizer optimizer(geom_i);
DiffractionGeometry geom_o = optimizer.Run(spots);
CHECK(fabs(geom_o.GetBeamX_pxl() - geom.GetBeamX_pxl()) < 0.001f);
CHECK(fabs(geom_o.GetBeamY_pxl() - geom.GetBeamY_pxl()) < 0.001f);
CHECK(fabs(geom_o.GetDetectorDistance_mm() - geom.GetDetectorDistance_mm()) < 0.001f);
CHECK(fabs(geom_o.GetPoniRot1_rad() - geom.GetPoniRot1_rad()) < 0.001f);
CHECK(fabs(geom_o.GetPoniRot2_rad() - geom.GetPoniRot2_rad()) < 0.001f);
}
TEST_CASE("DetGeomCalib_OptimizeGeometry") {
DiffractionGeometry geom;
geom.Wavelength_A(1.0).BeamX_pxl(1000.0).BeamY_pxl(1275.0)
.DetectorDistance_mm(100).PoniRot1_rad(0.1).PoniRot2_rad(0.05);
float lab6_a = LAB6_CELL_A;
std::vector<SpotToSave> spots;
for (int d = 1; d < 7; d++) {
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a / sqrt(d), i * M_PI * 2.0 / 30.0);
spots.push_back(SpotToSave{.x = x, .y = y});
}
}
DiffractionGeometry geom_i;
geom_i.Wavelength_A(1.0).BeamX_pxl(995.0).BeamY_pxl(1277.0)
.DetectorDistance_mm(98).PoniRot1_rad(0.0975).PoniRot2_rad(0.055);
OptimizeGeometry(geom_i, spots, CalculateXtalRings(UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90)));
CHECK(geom_i.GetBeamX_pxl() == Catch::Approx(geom.GetBeamX_pxl()));
CHECK(geom_i.GetBeamY_pxl() == Catch::Approx(geom.GetBeamY_pxl()));
CHECK(geom_i.GetDetectorDistance_mm() == Catch::Approx(geom.GetDetectorDistance_mm()));
CHECK(geom_i.GetPoniRot1_rad() == Catch::Approx(geom.GetPoniRot1_rad()));
CHECK(geom_i.GetPoniRot2_rad() == Catch::Approx(geom.GetPoniRot2_rad()));
}
// The tilt is held where it was given, and the beam centre and distance are still fitted. This is what
// a calibration destined for a program with no way to express a tilted detector asks for.
TEST_CASE("DetGeomCalib_RingOptimizerFixedTilt") {
DiffractionGeometry geom;
geom.Wavelength_A(1.0).BeamX_pxl(100.0).BeamY_pxl(200.0).DetectorDistance_mm(100);
const float lab6_a = LAB6_CELL_A;
std::vector<RingOptimizerInput> spots;
for (int d = 1; d < 4; d++) {
for (int i = 0; i < 30; i++) {
auto [x, y] = geom.ResPhiToPxl(lab6_a / sqrtf(d), i * M_PI * 2.0 / 30.0);
spots.push_back(RingOptimizerInput{.x = x, .y = y,
.q_expected = sqrtf(d) * M_PI * 2.0 / lab6_a});
}
}
DiffractionGeometry geom_i;
geom_i.Wavelength_A(1.0).BeamX_pxl(105).BeamY_pxl(195).DetectorDistance_mm(110);
DiffractionGeometry geom_o = RingOptimizer(geom_i, false).Run(spots);
CHECK(geom_o.GetPoniRot1_rad() == 0.0f);
CHECK(geom_o.GetPoniRot2_rad() == 0.0f);
CHECK(geom_o.GetBeamX_pxl() == Catch::Approx(geom.GetBeamX_pxl()));
CHECK(geom_o.GetBeamY_pxl() == Catch::Approx(geom.GetBeamY_pxl()));
CHECK(geom_o.GetDetectorDistance_mm() == Catch::Approx(geom.GetDetectorDistance_mm()));
// A tilt already in the geometry stays exactly where it is rather than being refined away.
DiffractionGeometry geom_tilted(geom_i);
geom_tilted.PoniRot1_rad(0.1).PoniRot2_rad(0.05);
DiffractionGeometry geom_t = RingOptimizer(geom_tilted, false).Run(spots);
CHECK(geom_t.GetPoniRot1_rad() == Catch::Approx(0.1));
CHECK(geom_t.GetPoniRot2_rad() == Catch::Approx(0.05));
}
Reference in New Issue View Git Blame Copy Permalink