From 659ba0d5b94a80fa3fe5c4900c3a6908dc7ef09e Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sat, 22 Aug 2026 20:50:38 +0200 Subject: [PATCH] Cross-check a tilted detector against pyFAI and DIALS Nothing constrained the tilted geometry. Every existing test is either self-consistent or moves one angle at a time, and every file CI writes has zero tilt - where a swapped axis, a reversed composition order and a wrong pivot all give exactly the same answer. Two real bugs lived in that gap. Two checks, because there are two things to guard. DiffractionGeometry_Tilted_vs_PyFAI_and_DIALS pins the model itself: all three PONI angles non-zero and of mixed sign, compared per pixel against reference positions from pyFAI (an independent implementation of the convention) and from DIALS, to 2 um. Because the references are quoted in their own frames and the test applies the documented mappings - pyFAI (t1,t2,t3) -> our (x,y,z), and imgCIF = ours turned 180 degrees about x - it pins those relations too, not just the arithmetic. The comment gives the snippets to regenerate both sets. tests/nxmx_geometry_dials_test.py guards the writer, which is where the bugs actually were and which the unit test cannot reach. CI writes a master, patches the geometry in and asks DIALS where the panel is. Only the angle VALUES are patched; the axis vectors, the depends_on chain and the pivot stay as the writer emitted them, so they remain under test. Verified to fail on the pre-fix encoding: 7.3 mm, exit 1, naming the chain as the thing to look at. Recorded in both, because it cost an hour: compare via get_origin() and the fast/slow axes, NOT get_pixel_lab_coord(), which applies a parallax correction from the sensor thickness that Jungfraujoch does not model - about 0.1 mm at the detector edge, easily mistaken for a geometry error. Co-Authored-By: Claude Opus 5 (1M context) --- .gitea/workflows/build_and_test.yml | 10 +++ docs/CHANGELOG.md | 1 + tests/DiffractionGeometryTest.cpp | 85 ++++++++++++++++++++++- tests/nxmx_geometry_dials_test.py | 104 ++++++++++++++++++++++++++++ 4 files changed, 199 insertions(+), 1 deletion(-) create mode 100644 tests/nxmx_geometry_dials_test.py diff --git a/.gitea/workflows/build_and_test.yml b/.gitea/workflows/build_and_test.yml index fa954e31..233acf43 100644 --- a/.gitea/workflows/build_and_test.yml +++ b/.gitea/workflows/build_and_test.yml @@ -354,6 +354,16 @@ jobs: dials.import single_master.h5 | tee import.log grep -q "num images: 100" import.log || { echo "DIALS did not read all 100 images"; exit 1; } xia2.ssx image=single_master.h5 space_group=P43212 unit_cell=78.551,78.551,36.914,90.000,90.000,90.000 + # Everything above writes an UNTILTED detector, where every plausible encoding of the tilt + # agrees exactly - so none of it constrains the NXmx transformation chain. This does. + - name: Check tilted detector geometry against DIALS + shell: bash + run: | + source /opt/dials-v3-27-0/dials_env.sh + mkdir -p test04 + cd test04 + ../build/tools/jfjoch_hdf5_test ../tests/test_data/compression_benchmark.h5 -n1 -S -o geom + dials.python ../tests/nxmx_geometry_dials_test.py geom_master.h5 xds-durin-test: name: XDS test (durin plugin) runs-on: jfjoch_rocky9 diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index a69fc632..958bde74 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -11,6 +11,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: `_process.h5` now describes the pixel format of the images it links to, instead of the container rugnux processes in. * HDF5: the detector tilt (`rot1`/`rot2`/`rot3`) is now exported correctly in the NXmx transformation chain; a tilted geometry was previously written so that other programs placed the detector wrongly (tens of mm at a few degrees of tilt). Untilted geometries are unaffected. * rugnux: the `.poni` file written by `--mode calibration` no longer negates `Rot3`, which exported a detector rotation about the beam with the wrong sign. +* Tests: a tilted detector is now cross-checked against pyFAI and DIALS, in the unit tests and against a written file in CI. ### 1.0.0-rc.161 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. diff --git a/tests/DiffractionGeometryTest.cpp b/tests/DiffractionGeometryTest.cpp index 872f6963..cb9da7d6 100644 --- a/tests/DiffractionGeometryTest.cpp +++ b/tests/DiffractionGeometryTest.cpp @@ -547,4 +547,87 @@ TEST_CASE("DiffractionGeometry_PONI_matrix_consistency") { CHECK(S0.x == Catch::Approx(0.0f)); CHECK(S0.y == Catch::Approx(0.0f)); CHECK(S0.z == Catch::Approx(1.0f)); -} \ No newline at end of file +} +// Cross-check of a TILTED detector against two independent implementations, pyFAI and DIALS/dxtbx. +// +// Every other geometry test here is either self-consistent (round trips) or exercises one angle at a +// time. This one pins all three PONI angles at once, non-zero and of mixed sign, against reference +// positions computed outside Jungfraujoch. That matters because the errors this guards against are +// second order: a wrong composition order or a swapped axis is invisible unless two angles are +// non-zero simultaneously, and a wrong pivot is invisible to anything that only checks directions. +// +// HOW TO REGENERATE THE NUMBERS +// +// pyFAI (`pip install pyFAI`), which is an independent implementation of the PONI convention: +// +// from pyFAI.geometry import Geometry +// from pyFAI.detectors import Detector +// px = 75e-6 +// det = Detector(pixel1=px, pixel2=px, max_shape=(2164, 2030)) +// g = Geometry(dist=0.150, poni1=1275*px + px/2, poni2=1000*px + px/2, +// rot1=0.05, rot2=+0.03, rot3=0.02, # (+rot1, -rot2, +rot3); see WritePoniFile +// detector=det, wavelength=1e-10) +// t3, t1, t2 = g.calc_pos_zyx(d1=[y], d2=[x]) # metres, pyFAI's own axes +// lab_mm = (t2*1e3, t1*1e3, t3*1e3) # pyFAI (t1,t2,t3) -> our (x,y,z) +// +// The half pixel in poni1/poni2 is the origin convention (docs/DETECTOR_GEOMETRY.md): our beam +// centre is pixel-centred, pyFAI measures from the sensor edge. +// +// DIALS: write a master, patch this geometry into it, and read the panel back. +// +// source /opt/dials-v3-27-0/dials_env.sh +// build/tools/jfjoch_hdf5_test -n1 -S -o g # writes g_master.h5 +// # with h5py, set /entry/instrument/detector/{beam_center_x,beam_center_y,distance}, +// # transformations/{rot1,rot2,rot3}, and recompute transformations/translation - both its +// # magnitude and its @vector - as {bx*px, by*px, distance} normalised, since the writer +// # derives it from the beam centre and distance. +// p = ExperimentListFactory.from_filenames(['g_master.h5'])[0].detector[0] +// lab = p.get_origin() + x*px_mm*p.get_fast_axis() + y*px_mm*p.get_slow_axis() +// +// Use get_origin()/get_fast_axis()/get_slow_axis() as above, NOT get_pixel_lab_coord(): that applies +// a parallax correction from the sensor thickness and material which DiffractionGeometry does not +// model, and it costs ~0.1 mm at the detector edge - enough to look like a geometry error. +// +// DIALS reports in the imgCIF frame, which is ours turned 180 degrees about x (diag(1,-1,-1)) - a +// proper rotation, not a mirror. The test applies that mapping, so it pins the frame relation too. +TEST_CASE("DiffractionGeometry_Tilted_vs_PyFAI_and_DIALS", "[DiffractionGeometry]") { + DiffractionGeometry geom; + geom.BeamX_pxl(1000.0f).BeamY_pxl(1275.0f).DetectorDistance_mm(150.0f) + .PixelSize_mm(0.075f).Wavelength_A(1.0f) + .PoniRot1_rad(0.05f).PoniRot2_rad(-0.03f).PoniRot3_rad(0.02f); + + struct Reference { + int x, y; + float pyfai[3]; // our frame: x, y, z [mm] + float dials[3]; // imgCIF frame: x, y, z [mm] + }; + + const std::vector reference = { + { 0, 0, {-69.399541334f, -98.819975327f, 150.623559791f}, + {-69.399544982f, 98.819979800f, -150.623559832f}}, + { 2029, 2163, { 85.802269836f, 60.488193357f, 147.887421982f}, + { 85.802273560f, -60.488196450f, -147.887421893f}}, + { 1000, 1275, { 7.405508016f, -4.642730847f, 149.745128473f}, + { 7.405508016f, 4.642730847f, -149.745128473f}}, + { 300, 1800, {-44.232874953f, 35.676608756f, 153.548927011f}, + {-44.232877406f, -35.676610671f, -153.548927191f}}, + { 1700, 400, { 58.519162201f, -71.195011374f, 145.153948055f}, + { 58.519164629f, 71.195014535f, -145.153947836f}}, + }; + + // 2 um, i.e. 1/37 of a pixel. The references agree with each other to ~4e-6 mm; the margin is + // set by float32 rounding in DiffractionGeometry and in the HDF5 file the DIALS values came from. + const double margin = 2e-3; + + for (const auto &r: reference) { + const Coord lab = geom.LabCoord(static_cast(r.x), static_cast(r.y)); + + CHECK(lab.x == Catch::Approx(r.pyfai[0]).margin(margin)); + CHECK(lab.y == Catch::Approx(r.pyfai[1]).margin(margin)); + CHECK(lab.z == Catch::Approx(r.pyfai[2]).margin(margin)); + + CHECK(lab.x == Catch::Approx( r.dials[0]).margin(margin)); + CHECK(lab.y == Catch::Approx(-r.dials[1]).margin(margin)); + CHECK(lab.z == Catch::Approx(-r.dials[2]).margin(margin)); + } +} diff --git a/tests/nxmx_geometry_dials_test.py b/tests/nxmx_geometry_dials_test.py new file mode 100644 index 00000000..3a1b374a --- /dev/null +++ b/tests/nxmx_geometry_dials_test.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python +"""Check that DIALS places a TILTED detector where Jungfraujoch does. + +Run under the DIALS environment: + + source /opt/dials-v3-27-0/dials_env.sh + dials.python tests/nxmx_geometry_dials_test.py + +Why this exists: the NXmx transformation chain the writer produces is only exercised by the rest of +CI at zero tilt, where every plausible encoding of the tilt - right axes or swapped, right +composition order or reversed, right pivot or not - gives exactly the same answer. The errors are +second order and only appear with two angles non-zero at once, or with a non-zero distance being +rotated. A previous encoding placed the detector 40 mm out at a few degrees of tilt while passing +every zero-tilt check. + +The reference positions below are the same ones pinned in +tests/DiffractionGeometryTest.cpp ("DiffractionGeometry_Tilted_vs_PyFAI_and_DIALS"), which checks +them against pyFAI as well; that comment explains how to regenerate them. Keep the two in step. + +The script patches the geometry into an existing master rather than requiring the writing tool to +take tilt options: the angle VALUES are just scalars, while the axis vectors, the depends_on chain +and the pivot - the parts that were wrong - stay exactly as the writer emitted them, so they remain +under test. +""" + +import shutil +import sys + +import h5py +import numpy as np +from dxtbx.model.experiment_list import ExperimentListFactory + +# Geometry to impose. Mixed signs on purpose, and all three angles non-zero. +BEAM_X, BEAM_Y, DISTANCE_M = 1000.0, 1275.0, 0.150 +ROT1, ROT2, ROT3 = 0.05, -0.03, 0.02 + +# Expected lab position per pixel, in the imgCIF frame DIALS reports (mm). +EXPECTED = { + (0, 0): (-69.399544982, 98.819979800, -150.623559832), + (1000, 1275): (7.405508016, 4.642730847, -149.745128473), + (300, 1800): (-44.232877406, -35.676610671, -153.548927191), + (1700, 400): (58.519164629, 71.195014535, -145.153947836), +} + +# 2 um, matching the C++ test. A wrong axis, order or pivot is a millimetre-scale error. +TOLERANCE_MM = 2e-3 + + +def main(master): + patched = "geometry_check_master.h5" + shutil.copy(master, patched) + + with h5py.File(patched, "r+") as h: + det = h["/entry/instrument/detector"] + det["beam_center_x"][()] = BEAM_X + det["beam_center_y"][()] = BEAM_Y + det["distance"][()] = DISTANCE_M + transformations = det["transformations"] + transformations["rot1"][()] = ROT1 + transformations["rot2"][()] = ROT2 + transformations["rot3"][()] = ROT3 + # The writer derives the translation from the beam centre and the distance, so it has to be + # recomputed here - both the magnitude and the direction it points in. + pixel_m = float(det["x_pixel_size"][()]) + vector = np.array([BEAM_X * pixel_m, BEAM_Y * pixel_m, DISTANCE_M]) + length = np.linalg.norm(vector) + transformations["translation"][()] = length + transformations["translation"].attrs["vector"] = vector / length + + panel = ExperimentListFactory.from_filenames([patched])[0].detector[0] + origin = np.array(panel.get_origin()) + fast = np.array(panel.get_fast_axis()) + slow = np.array(panel.get_slow_axis()) + pixel_mm = panel.get_pixel_size()[0] + + # Deliberately NOT get_pixel_lab_coord(): it applies a parallax correction from the sensor + # thickness and material that Jungfraujoch does not model, worth ~0.1 mm at the detector edge. + failures = [] + for (x, y), expected in sorted(EXPECTED.items()): + got = origin + x * pixel_mm * fast + y * pixel_mm * slow + error = np.max(np.abs(got - np.array(expected))) + status = "ok" if error <= TOLERANCE_MM else "FAILED" + print(f" pixel ({x:5d},{y:5d}): max error {error:.2e} mm {status}") + if error > TOLERANCE_MM: + failures.append(f"pixel ({x},{y}): got {got}, expected {expected}") + + if failures: + print("\nDIALS does not place the tilted detector where Jungfraujoch does:") + for failure in failures: + print(" " + failure) + print("\nThe NXmx transformation chain in writer/HDF5NXmx.cpp is the thing to look at: the " + "rot1/rot2/rot3 axis vectors, the order of the depends_on chain, and whether " + "translation sits inside the rotations so the tilt pivots about the sample.") + return 1 + + print(f"\nAll {len(EXPECTED)} pixels agree within {TOLERANCE_MM} mm.") + return 0 + + +if __name__ == "__main__": + if len(sys.argv) != 2: + print(__doc__) + sys.exit(2) + sys.exit(main(sys.argv[1]))