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) <noreply@anthropic.com>

.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

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.

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));
-}
+}
+// 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 <input.h5> -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> 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<float>(r.x), static_cast<float>(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));
+    }
+}

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 <master.h5>
+
+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]))