From fed077e683feffc70a6f91b5cee03a0119f1eadb Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 29 Aug 2026 22:11:09 +0200 Subject: [PATCH] geometry: hold the detector plane as axis vectors, and give the discrete part its own home The detector plane was three PONI angles and nothing else, so the two things it cannot express - an image mirrored in Y, and one mounted at a multiple of 90 degrees - had no home at all. They are now the DetectorOrientation carried by the detector setup, composed with the PONI rotation into one orthogonal matrix whose columns ARE the fast axis, the slow axis and the sample->PONI normal: lab = R(rot1, rot2, rot3) * Delta * ( (x-bx)*p , (y-by)*p , distance ) GetFastAxis/GetSlowAxis/GetNormalAxis read those columns and DetectorAxes() sets the plane from them, decomposing back to the angles; PoniRotMatrix and PoniAnglesFromMatrix are the conversion in both directions, exact on the canonical branch (rot2 in [-pi/2, pi/2]) and with a stated convention at gimbal lock. The angles stay stored rather than re-derived, so a geometry given as angles is written back as the same angles, to the bit. Delta is never inferred. In particular an arbitrary rot3 is NOT decomposed into a quarter turn plus a residual: rot3 is a fitted quantity, and a least-squares step must not be able to turn the stored image. It is set only where something states it - the detector setup, --detector-mirror-y / --detector-quarter-turns, or the value a file this system wrote records - and defaults to the identity, which makes the whole change a no-op for every existing detector and every existing file. It is a different setting from DetectorSetup::mirror_y, which flips the MODULE LAYOUT while an image is assembled and so decides what the stored pixels are. Merging the two would apply the mirror twice for every modular detector, or change the pixel content of every file written; both are ruled out. The new one earns its keep exactly where the old one is a no-op: a detector whose image arrives already assembled has no layout to flip. Both generators are signed permutations of the in-plane offset, so they preserve the distance from the PONI. That is why almost nothing downstream changes: everything needing an azimuth already goes through LabCoord, and everything that does not needs only a radius. The two hand-written copies of the rotation - XtalResidual and RingOptimizer - take the discrete part as four constants next to cos_rot3/sin_rot3, since it acts in the detector frame where rot3 acts in the laboratory and cannot be folded into it. RingOptimizer needs it despite being a radial fit: it fits the tilt, and the discrete part changes which way the tilt tips a ring. Carried as two optional CBOR keys and two detectorSpecific datasets, both back-compatible; the NXmx module axis vectors and the translation direction stop being hardcoded and are computed from it, reproducing today's values exactly at the identity. GetPoniRotMatrix is renamed GetDetectorMatrix, because it is no longer only the PONI rotation. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Lc5JG6kJqZoCWaoZ43JGTW --- common/CMakeLists.txt | 2 + common/Coord.cpp | 12 ++ common/Coord.h | 4 + common/DetectorOrientation.cpp | 51 +++++ common/DetectorOrientation.h | 42 ++++ common/DetectorSetup.cpp | 9 + common/DetectorSetup.h | 8 + common/DiffractionExperiment.cpp | 3 + common/DiffractionGeometry.cpp | 78 ++++++- common/DiffractionGeometry.h | 31 ++- common/JFJochMessages.h | 6 + docs/CBOR.md | 2 + docs/CHANGELOG.md | 1 + docs/DETECTOR_GEOMETRY.md | 39 ++++ docs/HDF5.md | 7 +- docs/RUGNUX.md | 2 + frame_serialize/CBORStream2Deserializer.cpp | 4 + frame_serialize/CBORStream2Serializer.cpp | 5 + .../bragg_prediction/BraggPrediction.cpp | 4 +- .../bragg_prediction/BraggPredictionGPU.cu | 2 +- .../bragg_prediction/BraggPredictionRot.cpp | 4 +- .../bragg_prediction/BraggPredictionRotGPU.cu | 2 +- .../BeamCenterFromBackground.cpp | 2 +- .../geom_refinement/BeamCenterFromSpots.cpp | 2 +- .../geom_refinement/GeometryRefiner.cpp | 3 +- image_analysis/geom_refinement/PostRefine.cpp | 5 +- .../geom_refinement/RingOptimizer.cpp | 25 ++- .../geom_refinement/XtalOptimizer.cpp | 3 +- image_analysis/geom_refinement/XtalResidual.h | 23 +- reader/HDF5MetadataSource.cpp | 10 + rugnux/rugnux_cli.cpp | 22 ++ tests/DiffractionGeometryTest.cpp | 202 +++++++++++++++++- viewer/JFJochHttpReader.cpp | 2 + writer/HDF5NXmx.cpp | 30 ++- 34 files changed, 605 insertions(+), 42 deletions(-) create mode 100644 common/DetectorOrientation.cpp create mode 100644 common/DetectorOrientation.h diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index 3b7b35082..59aaf46c3 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -104,6 +104,8 @@ ADD_LIBRARY(JFJochCommon STATIC GoniometerAxis.h DetectorTransformation.cpp DetectorTransformation.h + DetectorOrientation.cpp + DetectorOrientation.h CompressedImage.cpp CompressedImage.h Reflection.h diff --git a/common/Coord.cpp b/common/Coord.cpp index 8595c40ad..c43c67a1e 100644 --- a/common/Coord.cpp +++ b/common/Coord.cpp @@ -177,6 +177,18 @@ RotMatrix::RotMatrix(float alpha, const Coord &dir) { v[2][2] = t * n.z * n.z + c; } +RotMatrix::RotMatrix(const Coord &col0, const Coord &col1, const Coord &col2) { + for (int i = 0; i < 3; i++) { + v[i][0] = col0[i]; + v[i][1] = col1[i]; + v[i][2] = col2[i]; + } +} + +Coord RotMatrix::Column(int64_t i) const { + return {v[0][i], v[1][i], v[2][i]}; +} + Coord RotMatrix::operator*(const Coord &in) const { return { v[0][0] * in.x + v[0][1] * in.y + v[0][2] * in.z, diff --git a/common/Coord.h b/common/Coord.h index fd2f06d6a..ab0367405 100644 --- a/common/Coord.h +++ b/common/Coord.h @@ -49,6 +49,10 @@ class RotMatrix { public: RotMatrix(); RotMatrix(float alpha, const Coord &dir); + // From three columns. Only the arithmetic is enforced, not orthogonality - this is how a + // detector orientation is built from its fast, slow and normal axes. + RotMatrix(const Coord &col0, const Coord &col1, const Coord &col2); + [[nodiscard]] Coord Column(int64_t i) const; Coord operator*(const Coord &in) const; RotMatrix operator*(const RotMatrix &other) const; diff --git a/common/DetectorOrientation.cpp b/common/DetectorOrientation.cpp new file mode 100644 index 000000000..4e7095c4c --- /dev/null +++ b/common/DetectorOrientation.cpp @@ -0,0 +1,51 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "DetectorOrientation.h" +#include "JFJochException.h" + +DetectorOrientation::DetectorOrientation(bool mirror_y, int64_t quarter_turns) { + MirrorY(mirror_y); + QuarterTurns(quarter_turns); +} + +DetectorOrientation &DetectorOrientation::MirrorY(bool input) { + mirror_y = input; + return *this; +} + +DetectorOrientation &DetectorOrientation::QuarterTurns(int64_t input) { + if ((input < 0) || (input > 3)) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Quarter turns must be 0, 1, 2 or 3"); + quarter_turns = input; + return *this; +} + +bool DetectorOrientation::IsMirrorY() const { + return mirror_y; +} + +int64_t DetectorOrientation::GetQuarterTurns() const { + return quarter_turns; +} + +bool DetectorOrientation::IsIdentity() const { + return !mirror_y && (quarter_turns == 0); +} + +RotMatrix DetectorOrientation::Matrix() const { + // Columns of Rz(k*90 deg), in the internal frame (x = column, y = row downward, z = beam). y + // points down, so a positive right-handed turn about +z takes +x to +y - clockwise on screen. + const float c[4] = {1, 0, -1, 0}; + const float s[4] = {0, 1, 0, -1}; + const Coord rz_x = {c[quarter_turns], s[quarter_turns], 0}; + const Coord rz_y = {-s[quarter_turns], c[quarter_turns], 0}; + + // Rz * diag(1,-1,1): the mirror negates the second column. + return {rz_x, mirror_y ? -rz_y : rz_y, {0, 0, 1}}; +} + +bool DetectorOrientation::operator==(const DetectorOrientation &other) const { + return (mirror_y == other.mirror_y) && (quarter_turns == other.quarter_turns); +} diff --git a/common/DetectorOrientation.h b/common/DetectorOrientation.h new file mode 100644 index 000000000..e9a48a2ce --- /dev/null +++ b/common/DetectorOrientation.h @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "Coord.h" + +// How the stored image is laid out in the detector plane: mirrored in Y, and/or turned by a multiple +// of 90 degrees about the beam. Both are exact pixel remappings, which an arbitrary in-plane rotation +// (PONI rot3) is not - a viewer can show the image the right way up from these two without resampling +// anything. +// +// It is applied to the offset from the PONI, in stored-image millimetres, BEFORE the continuous PONI +// tilt takes that offset to the laboratory: +// +// lab = R(rot1, rot2, rot3) * Matrix() * ( (x-beam_x)*pixel, (y-beam_y)*pixel, distance ) +// +// Mirror first, then the quarter turns. Every element of the group the two generate can be written +// that way, so the order is a convention rather than a derivation, and this is the one. +// +// This is NOT the same thing as DetectorSetup::mirror_y, which describes the raw readout -> assembled +// image module layout and is spent before the geometry sees anything. This one describes the assembled +// image -> detector canonical frame, changes no pixel, and defaults to the identity. +class DetectorOrientation { + bool mirror_y = false; + int64_t quarter_turns = 0; // 0..3, right-handed about the beam = clockwise on the displayed image +public: + DetectorOrientation() = default; + DetectorOrientation(bool mirror_y, int64_t quarter_turns); + + DetectorOrientation& MirrorY(bool input); + DetectorOrientation& QuarterTurns(int64_t input); + + [[nodiscard]] bool IsMirrorY() const; + [[nodiscard]] int64_t GetQuarterTurns() const; + [[nodiscard]] bool IsIdentity() const; + + // Rz(quarter_turns * 90 deg) * diag(1,-1,1)^mirror_y. Entries are exactly 0 and +-1. + [[nodiscard]] RotMatrix Matrix() const; + + bool operator==(const DetectorOrientation &other) const; +}; diff --git a/common/DetectorSetup.cpp b/common/DetectorSetup.cpp index 7de513437..04d2edbea 100644 --- a/common/DetectorSetup.cpp +++ b/common/DetectorSetup.cpp @@ -277,6 +277,15 @@ bool DetectorSetup::IsMirrorY() const { return mirror_y; } +DetectorSetup &DetectorSetup::ImageOrientation(const DetectorOrientation &input) { + image_orientation = input; + return *this; +} + +DetectorOrientation DetectorSetup::GetImageOrientation() const { + return image_orientation; +} + DetectorSetup & DetectorSetup::ReadOutTime(std::chrono::nanoseconds input) { if (input.count() < 0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, diff --git a/common/DetectorSetup.h b/common/DetectorSetup.h index 92a36b853..a7095ffc4 100644 --- a/common/DetectorSetup.h +++ b/common/DetectorSetup.h @@ -8,6 +8,7 @@ #include #include "DetectorGeometry.h" +#include "DetectorOrientation.h" #include "../jungfrau/JFModuleGainCalibration.h" #include "DetectorGeometryFixed.h" #include "DetectorGeometryModular.h" @@ -40,6 +41,11 @@ class DetectorSetup { // A property of the configuration, not something derivable from the assembled image: it says // how the modules were laid out to reach the MX convention of row 0 at the top. bool mirror_y = true; + // How the ASSEMBLED image sits in the detector plane - mirrored in Y, turned by a multiple of 90 + // degrees. A different thing from mirror_y above, which is spent on the module layout before the + // geometry sees anything; this one changes no pixel, only how a pixel coordinate is taken to the + // laboratory. Identity by default, which is every detector assembled by this system. + DetectorOrientation image_orientation; std::chrono::nanoseconds read_out_time; std::chrono::nanoseconds min_count_time; std::chrono::nanoseconds min_frame_time; @@ -79,6 +85,7 @@ public: DetectorSetup& BaseIPv4Addr(const std::string &input); DetectorSetup& ModuleSync(bool input); DetectorSetup& MirrorY(bool input); + DetectorSetup& ImageOrientation(const DetectorOrientation &input); DetectorSetup& ReadOutTime(std::chrono::nanoseconds input); DetectorSetup& Geometry(const DetectorGeometryFixed& input); DetectorSetup& BitDepthImage(int64_t input); @@ -112,6 +119,7 @@ public: [[nodiscard]] std::string GetBaseIPv4Addr() const; [[nodiscard]] bool IsModuleSync() const; [[nodiscard]] bool IsMirrorY() const; + [[nodiscard]] DetectorOrientation GetImageOrientation() const; [[nodiscard]] std::chrono::nanoseconds GetReadOutTime() const; [[nodiscard]] std::chrono::nanoseconds GetMinFrameTime() const; [[nodiscard]] std::chrono::nanoseconds GetMinCountTime() const; diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index e257d871b..719f28e31 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -660,6 +660,8 @@ void DiffractionExperiment::FillMessage(StartMessage &message) const { message.image_size_x = GetXPixelsNum(); message.image_size_y = GetYPixelsNum(); message.mirror_y = IsDetectorMirroredY(); + message.detector_orientation_mirror_y = detector.GetImageOrientation().IsMirrorY(); + message.detector_orientation_quarter_turns = detector.GetImageOrientation().GetQuarterTurns(); message.saturation_value = SaturationValueFromLimit(GetSaturationLimit()); // The marker actually stored in the pixels: UINTx_MAX unsigned, INTx_MIN signed. GetUnderflow() // was written here, which is -1 for an unsigned image and so matches no pixel it can contain. @@ -1598,6 +1600,7 @@ DiffractionGeometry DiffractionExperiment::GetDiffractionGeometry() const { .PoniRot1_rad(dataset.GetPoniRot1_rad()) .PoniRot2_rad(dataset.GetPoniRot2_rad()) .PoniRot3_rad(dataset.GetPoniRot3_rad()) + .Orientation(detector.GetImageOrientation()) .Rotation(dataset.GetGoniometer()); return g; } diff --git a/common/DiffractionGeometry.cpp b/common/DiffractionGeometry.cpp index 9f3ed32c5..085cc5997 100644 --- a/common/DiffractionGeometry.cpp +++ b/common/DiffractionGeometry.cpp @@ -2,16 +2,39 @@ // SPDX-License-Identifier: GPL-3.0-only #include "JFJochMath.h" +#include #include #include "DiffractionGeometry.h" #include "RawToConvertedGeometry.h" +RotMatrix PoniRotMatrix(float rot1, float rot2, float rot3) { + return RotMatrix(-rot3, {0,0,1}) + * RotMatrix(-rot2, {1,0,0}) + * RotMatrix(rot1, {0,1,0}); +} + +void PoniAnglesFromMatrix(const RotMatrix &rot_matrix, float &rot1, float &rot2, float &rot3) { + const Coord fast = rot_matrix.Column(0); + const Coord slow = rot_matrix.Column(1); + const Coord normal = rot_matrix.Column(2); + + rot2 = asinf(std::clamp(-slow.z, -1.0f, 1.0f)); + if (fabsf(cosf(rot2)) < 1e-6f) { + // Gimbal lock: only rot1 +- rot3 is determined, so put it all into rot1. + rot1 = atan2f(normal.x, fast.x); + rot3 = 0.0f; + } else { + rot1 = atan2f(-fast.z, normal.z); + rot3 = atan2f(slow.x, slow.y); + } +} + Coord DiffractionGeometry::LabCoord(float x, float y) const { Coord detectorCoord = {(x - beam_x_pxl) * pixel_size_mm , (y - beam_y_pxl) * pixel_size_mm , det_distance_mm}; - return poni_rot * detectorCoord; + return det_matrix * detectorCoord; } std::pair DiffractionGeometry::GetDirectBeam_pxl() const { @@ -28,7 +51,7 @@ Coord DiffractionGeometry::DetectorToRecip(float x, float y) const { std::pair DiffractionGeometry::RecipToDetector(const Coord &recip) const { auto S_unrotated = recip + GetScatteringVector(); - auto S = poni_rot.transpose() * S_unrotated; + auto S = det_matrix.transpose() * S_unrotated; if (S.z <= 0) return {NAN, NAN}; @@ -180,27 +203,25 @@ float DiffractionGeometry::AngleFromEwaldSphere_deg(const Coord &p0) const { return angle_deg(p_star, p0); } -void DiffractionGeometry::UpdatePoniRotMatrix() { - poni_rot = RotMatrix(-poni_rot_3, {0,0,1}) - * RotMatrix(-poni_rot_2, {1,0,0}) - * RotMatrix(poni_rot_1, {0,1,0}); +void DiffractionGeometry::UpdateDetectorMatrix() { + det_matrix = PoniRotMatrix(poni_rot_1, poni_rot_2, poni_rot_3) * orientation.Matrix(); } DiffractionGeometry &DiffractionGeometry::PoniRot1_rad(float input) { poni_rot_1 = input; - UpdatePoniRotMatrix(); + UpdateDetectorMatrix(); return *this; } DiffractionGeometry &DiffractionGeometry::PoniRot2_rad(float input) { poni_rot_2 = input; - UpdatePoniRotMatrix(); + UpdateDetectorMatrix(); return *this; } DiffractionGeometry &DiffractionGeometry::PoniRot3_rad(float input) { poni_rot_3 = input; - UpdatePoniRotMatrix(); + UpdateDetectorMatrix(); return *this; } @@ -216,6 +237,41 @@ float DiffractionGeometry::GetPoniRot3_rad() const { return poni_rot_3; } +DiffractionGeometry &DiffractionGeometry::Orientation(const DetectorOrientation &input) { + orientation = input; + UpdateDetectorMatrix(); + return *this; +} + +DetectorOrientation DiffractionGeometry::GetOrientation() const { + return orientation; +} + +DiffractionGeometry &DiffractionGeometry::DetectorAxes(const Coord &fast, const Coord &slow) { + const Coord f = fast.Normalize(); + const Coord s = slow.Normalize(); + // The normal is not free: it is the sample->PONI direction, and whether it is +fast x slow or + // -fast x slow is exactly whether the stored image is mirrored, which the orientation already says. + const Coord n = orientation.IsMirrorY() ? -(f % s) : (f % s); + + PoniAnglesFromMatrix(RotMatrix(f, s, n) * orientation.Matrix().transpose(), + poni_rot_1, poni_rot_2, poni_rot_3); + UpdateDetectorMatrix(); + return *this; +} + +Coord DiffractionGeometry::GetFastAxis() const { + return det_matrix.Column(0); +} + +Coord DiffractionGeometry::GetSlowAxis() const { + return det_matrix.Column(1); +} + +Coord DiffractionGeometry::GetNormalAxis() const { + return det_matrix.Column(2); +} + std::pair DiffractionGeometry::ResPhiToPxl(float d_A, float phi_rad) const { // Guard invalid inputs if (wavelength_A <= 0.0f || d_A <= wavelength_A / 2.0f) @@ -241,8 +297,8 @@ Coord DiffractionGeometry::ProjectToEwaldSphere(const Coord &p0) const { return S - S0; } -const RotMatrix &DiffractionGeometry::GetPoniRotMatrix() const { - return poni_rot; +const RotMatrix &DiffractionGeometry::GetDetectorMatrix() const { + return det_matrix; } std::optional DiffractionGeometry::GetRotation() const { diff --git a/common/DiffractionGeometry.h b/common/DiffractionGeometry.h index bf1924f20..92cd46cbd 100644 --- a/common/DiffractionGeometry.h +++ b/common/DiffractionGeometry.h @@ -5,8 +5,20 @@ #include "JFJochException.h" #include "Coord.h" +#include "DetectorOrientation.h" #include "GoniometerAxis.h" +// The two directions of the PONI convention, as pure functions, so the conversion can be exercised +// on its own. rot_matrix = Rz(-rot3) * Rx(-rot2) * Ry(+rot1) in the internal frame (x = column, +// y = row downward, z = beam); its columns are the lab directions of a +1 column step, a +1 row step +// and the sample->PONI vector. +RotMatrix PoniRotMatrix(float rot1, float rot2, float rot3); + +// The inverse. rot2 comes back in [-pi/2, pi/2] and rot1, rot3 in (-pi, pi], which is the canonical +// branch: on it the round trip is the identity. At rot2 = +-pi/2 only rot1 +- rot3 is determined, and +// the convention is to put it all into rot1 and leave rot3 at zero. +void PoniAnglesFromMatrix(const RotMatrix &rot_matrix, float &rot1, float &rot2, float &rot3); + class DiffractionGeometry { float beam_x_pxl = 0.0; float beam_y_pxl = 0.0; @@ -16,10 +28,14 @@ class DiffractionGeometry { float poni_rot_1 = 0.0f; float poni_rot_2 = 0.0f; float poni_rot_3 = 0.0f; - RotMatrix poni_rot; + DetectorOrientation orientation; + // The full detector orientation: the PONI rotation composed with the discrete image orientation. + // Its columns are the fast, slow and normal axes. Orthogonal, but improper when the image is + // mirrored, so transpose() is still its inverse. + RotMatrix det_matrix; std::optional axis; - void UpdatePoniRotMatrix(); + void UpdateDetectorMatrix(); public: DiffractionGeometry &BeamX_pxl(float input); DiffractionGeometry &BeamY_pxl(float input); @@ -29,6 +45,11 @@ public: DiffractionGeometry &PoniRot1_rad(float input); DiffractionGeometry &PoniRot2_rad(float input); DiffractionGeometry &PoniRot3_rad(float input); + DiffractionGeometry &Orientation(const DetectorOrientation &input); + // Sets the detector plane from its two axis vectors (unit, orthogonal). The discrete orientation + // is left as it is - it says how the image is stored, which two vectors cannot - and the PONI + // angles are re-derived so that the two views stay in step. + DiffractionGeometry &DetectorAxes(const Coord &fast, const Coord &slow); DiffractionGeometry &Rotation(const std::optional &input); [[nodiscard]] float GetBeamX_pxl() const; @@ -40,6 +61,10 @@ public: [[nodiscard]] float GetPoniRot1_rad() const; [[nodiscard]] float GetPoniRot2_rad() const; [[nodiscard]] float GetPoniRot3_rad() const; + [[nodiscard]] DetectorOrientation GetOrientation() const; + [[nodiscard]] Coord GetFastAxis() const; // lab direction of a +1 column step + [[nodiscard]] Coord GetSlowAxis() const; // lab direction of a +1 row step + [[nodiscard]] Coord GetNormalAxis() const; // sample -> PONI direction [[nodiscard]] std::pair GetDirectBeam_pxl() const; [[nodiscard]] std::optional GetRotation() const; @@ -62,5 +87,5 @@ public: // eq. 18 in https://journals.iucr.org/d/issues/2014/08/00/dz5332/index.html [[nodiscard]] float AngleFromEwaldSphere_deg(const Coord &p0) const; - [[nodiscard]] const RotMatrix& GetPoniRotMatrix() const; + [[nodiscard]] const RotMatrix& GetDetectorMatrix() const; }; diff --git a/common/JFJochMessages.h b/common/JFJochMessages.h index 59fb283e6..b8e0fd0d8 100644 --- a/common/JFJochMessages.h +++ b/common/JFJochMessages.h @@ -221,6 +221,12 @@ struct StartMessage { // the MX convention - row 0 at the top of the detector, seen from the sample - and is what // Jungfraujoch has always produced, so it is also what absence of the field means. bool mirror_y = true; + // How the assembled image sits in the detector plane, relative to the frame the PONI angles below + // are stated in: mirrored in Y, and/or turned by this many quarter turns about the beam. A + // different thing from mirror_y above - see DetectorSetup. Absence means the identity, which is + // what every stream written before these fields existed carries. + bool detector_orientation_mirror_y = false; + int64_t detector_orientation_quarter_turns = 0; uint64_t bit_depth_image; // user data std::optional bit_depth_readout; bool pixel_signed; // user data diff --git a/docs/CBOR.md b/docs/CBOR.md index c887b6f3d..9887d83e2 100644 --- a/docs/CBOR.md +++ b/docs/CBOR.md @@ -26,6 +26,8 @@ There are minor differences at the moment: | image_size_x | uint64 | Image width \[pixels\] | X | | image_size_y | uint64 | Image height \[pixels\] | X | | mirror_y | bool | Whether the assembled image is mirrored in Y relative to the detector's raw readout order. True is the MX convention - row 0 at the top of the detector seen from the sample - and is what absence of the key means | | +| detector_orientation_mirror_y | bool | Whether the assembled image is mirrored in Y relative to the frame the PONI angles are stated in. A different setting from `mirror_y` above, which is about the module layout; this one changes no pixel. Absence means false | | +| detector_orientation_quarter_turns | int | How many multiples of 90 degrees about the beam the assembled image is turned by, relative to the frame the PONI angles are stated in (0-3). Absence means 0 | | | incident_energy | float | X-ray energy \[eV\] | X | | incident_wavelength | float | X-ray wavelength \[Angstrom\] | X | | incident_wavelength_spread | float (optional) | FWHM of the X-ray wavelength distribution \[Angstrom\] (NXmx incident_wavelength_spread); omitted when the beam is monochromatic | | diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 50ecdf0f8..e2c0875ec 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -2,6 +2,7 @@ ## 1.0.0 ### 1.0.0-rc.166 * `jfjoch_viewer` opens PILATUS miniCBF sweeps - naming any frame opens the whole sweep - and can run a processing job on one. +* A detector whose stored image is mirrored in Y or mounted at a multiple of 90 degrees can be described as such, in the detector configuration or with `--detector-mirror-y` / `--detector-quarter-turns`, rather than having to be expressed as a detector rotation. * The rotation first pass refines twelve candidate lattices rather than four, so a correct cell that the pre-refinement ranking put fifth is still reached. * A candidate cell whose three rows are coplanar is rejected before refinement, instead of producing a not-a-number Jacobian and several hundred lines of solver output. * When the directions an FFT search shortlists all lie in one plane, one further transform is spent along the plane normal, which is where the missing row must be - so a crystal with a very long axis can index. diff --git a/docs/DETECTOR_GEOMETRY.md b/docs/DETECTOR_GEOMETRY.md index 63cb48b0c..c59b41e77 100644 --- a/docs/DETECTOR_GEOMETRY.md +++ b/docs/DETECTOR_GEOMETRY.md @@ -42,6 +42,45 @@ that half pixel; the pixel values the same run reports are ours. `Rot3` in that negated and turned by 180°: the half turn sets the azimuthal reference, because pyFAI's in-plane axes are the negatives of ours. It leaves 2θ untouched, so it moves only the azimuth. +## Inside: two axis vectors; outside: rot1/rot2/rot3 + +Internally the detector plane is one orthogonal matrix whose columns are the **fast axis** (the +laboratory direction of a +1 column step), the **slow axis** (+1 row step) and the **normal** (the +sample→PONI direction). Every geometry calculation — resolution, azimuth, polarization, prediction, +refinement — is that matrix applied to the offset of a pixel from the PONI. + +`rot1`/`rot2`/`rot3` remain the way the tilt is stated from outside, and the two views convert both +ways: `R = Rz(-rot3)·Rx(-rot2)·Ry(+rot1)` in the internal frame, and back from the columns as + +``` +rot2 = asin(-slow.z) rot1 = atan2(-fast.z, normal.z) rot3 = atan2(slow.x, slow.y) +``` + +with `rot2` in [-90°, 90°]. The angles are what is stored and what is written out, so a geometry +given as angles comes back exactly as it was given. + +## Mirrored and quarter-turned detectors + +On top of the continuous tilt the detector setup carries a **discrete image orientation**: whether the +stored image is mirrored in Y, and how many multiples of 90° about the beam it is turned by. It is +applied to the offset from the PONI before the tilt. + +The distinction matters because these two operations are exact pixel remappings — an image can be +shown the right way up without resampling anything — while an arbitrary in-plane rotation cannot. +`rot3` is therefore reserved for the genuinely arbitrary part: an in-plane angle is **never** +decomposed into a quarter turn plus a residual, and the discrete part is set only where something +states it (the detector configuration, `--detector-mirror-y` / `--detector-quarter-turns`, or the +value a Jungfraujoch-written file records). + +Both operations leave the distance from the PONI unchanged, so resolution, the solid-angle correction +and anything else that needs only a radius are unaffected by them. Polarization *is* affected, and +correctly so: it is computed from the azimuth in the **laboratory**, and what these operations change +is which pixel index lands at which laboratory azimuth. + +This is a different setting from the `mirror_y` in the JSON configuration below, which flips the +**module layout** while the image is being assembled and so decides what the stored pixels are. The +discrete image orientation changes no pixel at all. + ## Macromolecular crystallography convention for the vertical direction One place of confusion is the convention to have point (0,0) of the detector in the top left corner of the detector, with Y values increasing downwards. This is also consistent with computer image formats. diff --git a/docs/HDF5.md b/docs/HDF5.md index 8e02dd85c..db3a8f6b7 100644 --- a/docs/HDF5.md +++ b/docs/HDF5.md @@ -228,7 +228,10 @@ refined tilt into `rot1`/`rot2`/`rot3`; the broker writes the user-provided geom ### `/entry/instrument/detector/module` (NXdetector_module) `data_origin`, `data_size`, `fast_pixel_direction`, `slow_pixel_direction`, `module_offset` — all -NXmx (`fast/slow_pixel_direction` and `module_offset` carry transformation attributes). +NXmx (`fast/slow_pixel_direction` and `module_offset` carry transformation attributes). The two +pixel-direction vectors carry the discrete image orientation (mirror in Y, multiples of 90° about the +beam); for a detector this system assembled itself they are the McStas form of the internal +x and ++y, i.e. `(-1, 0, 0)` and `(0, -1, 0)`. ### `/entry/sample` (NXsample) @@ -500,6 +503,8 @@ group for compatibility with existing tooling: | `detector_distance` | m | duplicate of `distance` (Dectris/Neggia compatibility) | | `detector_number` | | detector identifier (Dectris convention) | | `mirror_y` (in `detectorSpecific`) | | whether the stored image is mirrored in Y relative to the raw readout; true is the MX convention (row 0 at the top) | +| `detector_orientation_mirror_y` (in `detectorSpecific`) | | whether the stored image is mirrored in Y relative to the frame `rot1`/`rot2`/`rot3` are stated in — a different setting from `mirror_y`, and one that changes no pixel | +| `detector_orientation_quarter_turns` (in `detectorSpecific`) | | multiples of 90° about the beam the stored image is turned by, relative to that same frame (0-3) | | `error_value` | | masked/error pixel sentinel: `UINTx_MAX` unsigned, `INTx_MIN` signed (NXmx has no equivalent). NXmx `underload_value` is written too: `INTx_MIN + 1` for signed, `0` for unsigned | | `bit_depth_image` | | stored image bit depth (DECTRIS convention, not NXmx). Equal to `bit_depth_readout` where that is written, i.e. for unsigned images | | `acquisition_type` | | always `triggered` (Dectris convention) | diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index e91a44de2..f8d7ccc55 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -1058,5 +1058,7 @@ Geometry overrides (defaults are taken from the input file; override them to rep | `--wavelength ` | Wavelength (Å) | | `--rot1 ` | PONI detector rotation 1 (rad) | | `--rot2 ` | PONI detector rotation 2 (rad) | +| `--detector-mirror-y` | The stored image is mirrored in Y relative to the frame the PONI angles are stated in | +| `--detector-quarter-turns <0-3>` | The stored image is turned by this many multiples of 90° about the beam relative to that same frame | | `--polarization ` | Polarization factor | | `--rotation-scale ` | Goniometer rotation scale: the stage turned `k` times the angle stored in the file (the commanded one). Applied to both passes, and overrides the scale rugnux fits for itself | diff --git a/frame_serialize/CBORStream2Deserializer.cpp b/frame_serialize/CBORStream2Deserializer.cpp index 458f53410..5b1605b9f 100644 --- a/frame_serialize/CBORStream2Deserializer.cpp +++ b/frame_serialize/CBORStream2Deserializer.cpp @@ -1318,6 +1318,10 @@ namespace { message.geometry_transformation_enabled = GetCBORBool(value); else if (key == "mirror_y") message.mirror_y = GetCBORBool(value); + else if (key == "detector_orientation_mirror_y") + message.detector_orientation_mirror_y = GetCBORBool(value); + else if (key == "detector_orientation_quarter_turns") + message.detector_orientation_quarter_turns = GetCBORInt(value); else if (key == "jungfrau_conversion_factor") message.jungfrau_conversion_factor = GetCBORFloat(value); else if (key == "arm_date") diff --git a/frame_serialize/CBORStream2Serializer.cpp b/frame_serialize/CBORStream2Serializer.cpp index 8aa7f343c..a80c97937 100644 --- a/frame_serialize/CBORStream2Serializer.cpp +++ b/frame_serialize/CBORStream2Serializer.cpp @@ -727,6 +727,11 @@ void CBORStream2Serializer::SerializeSequenceStart(const StartMessage& message) // Not a DECTRIS field - stream2 has nothing for the row direction, so a consumer that does not // know this key simply skips it and gets today's behaviour, which is what absence means. CBOR_ENC(mapEncoder, "mirror_y", message.mirror_y); + // Also not a DECTRIS field, and also skipped by a consumer that does not know it - absence means + // the identity, which is what every stream written before these keys existed carries. + CBOR_ENC(mapEncoder, "detector_orientation_mirror_y", message.detector_orientation_mirror_y); + CBOR_ENC(mapEncoder, "detector_orientation_quarter_turns", + message.detector_orientation_quarter_turns); CBOR_ENC_PIXEL_MASK(mapEncoder, message); CBOR_ENC_AZINT_MAP(mapEncoder, message); diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 7f939a1d5..de48b63a3 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -83,7 +83,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal const Coord Cstar = lattice.Cstar(); const Coord S0 = geom.GetScatteringVector(); - std::vector rot = geom.GetPoniRotMatrix().transpose().arr(); + std::vector rot = geom.GetDetectorMatrix().transpose().arr(); // Precompute detector geometry constants float beam_x = geom.GetBeamX_pxl(); @@ -168,7 +168,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal } } - // Inlined RecipToDector with rot1 and rot2 (rot3 = 0) + // Inlined RecipToDetector: the full transposed detector matrix, tilt and discrete orientation // Apply rotation matrix transpose float S_rot_x = rot[0] * S_x + rot[1] * S_y + rot[2] * S_z; float S_rot_y = rot[3] * S_x + rot[4] * S_y + rot[5] * S_z; diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu index 91a45ac74..aa9e90673 100644 --- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu @@ -176,7 +176,7 @@ namespace { kc.Cstar = lattice.Cstar(); kc.S0 = geom.GetScatteringVector(); kc.centering = centering; - auto rotT = geom.GetPoniRotMatrix().transpose().arr(); + auto rotT = geom.GetDetectorMatrix().transpose().arr(); for (int i = 0; i < 9; ++i) kc.rot[i] = rotT[i]; return kc; } diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.cpp b/image_analysis/bragg_prediction/BraggPredictionRot.cpp index 5ca4eab90..59fb41c81 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRot.cpp @@ -23,7 +23,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys const Coord Cstar = lattice.Cstar(); const Coord S0 = geom.GetScatteringVector(); - std::vector rot = geom.GetPoniRotMatrix().transpose().arr(); + std::vector rot = geom.GetDetectorMatrix().transpose().arr(); // Precompute detector geometry constants float beam_x = geom.GetBeamX_pxl(); @@ -148,7 +148,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys const float partiality = (std::erf((phi + half_wedge_angle_rad) * c1) - std::erf((phi - half_wedge_angle_rad) * c1)) / 2.0f; - // Inlined RecipToDector with rot1 and rot2 (rot3 = 0) + // Inlined RecipToDetector: the full transposed detector matrix, tilt and discrete orientation // Apply rotation matrix transpose float S_rot_x = rot[0] * S.x + rot[1] * S.y + rot[2] * S.z; float S_rot_y = rot[3] * S.x + rot[4] * S.y + rot[5] * S.z; diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index c67e44136..63126406b 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -239,7 +239,7 @@ namespace { kc.Cstar = lattice.Cstar(); kc.S0 = geom.GetScatteringVector(); - auto rotT = geom.GetPoniRotMatrix().transpose().arr(); + auto rotT = geom.GetDetectorMatrix().transpose().arr(); for (int i = 0; i < 9; ++i) kc.rot[i] = rotT[i]; kc.centering = settings.centering; diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp index bba29ad9a..c8e4ba95a 100644 --- a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp @@ -69,7 +69,7 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe const float tt_hi = 2.0f * std::asin(sin_high); const float d_tt = (tt_hi - tt_lo) / RADIAL_BINS; - const auto rot = geom.GetPoniRotMatrix().arr(); // row major + const auto rot = geom.GetDetectorMatrix().arr(); // row major const float pixel_size = geom.GetPixelSize_mm(); const float distance = geom.GetDetectorDistance_mm(); diff --git a/image_analysis/geom_refinement/BeamCenterFromSpots.cpp b/image_analysis/geom_refinement/BeamCenterFromSpots.cpp index 6f67cc844..baf96e118 100644 --- a/image_analysis/geom_refinement/BeamCenterFromSpots.cpp +++ b/image_analysis/geom_refinement/BeamCenterFromSpots.cpp @@ -172,7 +172,7 @@ class LabMirror { const int axis; public: LabMirror(const DiffractionGeometry &geometry, int mirror_axis) - : geom(geometry), inverse_rotation(geometry.GetPoniRotMatrix().transpose()), axis(mirror_axis) {} + : geom(geometry), inverse_rotation(geometry.GetDetectorMatrix().transpose()), axis(mirror_axis) {} // The image of a detector point: negate the mirrored lab component of the ray to it and project // the result back onto the detector. The reflecting plane contains the beam, so a point's image diff --git a/image_analysis/geom_refinement/GeometryRefiner.cpp b/image_analysis/geom_refinement/GeometryRefiner.cpp index edc370f77..45630b291 100644 --- a/image_analysis/geom_refinement/GeometryRefiner.cpp +++ b/image_analysis/geom_refinement/GeometryRefiner.cpp @@ -148,6 +148,7 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge const double rot3 = nominal_geom.GetPoniRot3_rad(); // Same for every spot of every frame, so taken once here rather than per residual. const double cos_rot3 = std::cos(rot3), sin_rot3 = std::sin(rot3); + const DetectorOrientation orientation = nominal_geom.GetOrientation(); const double lambda = nominal_geom.GetWavelength_A(); const double pixel = nominal_geom.GetPixelSize_mm(); @@ -210,7 +211,7 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge problem.AddResidualBlock( new ceres::AutoDiffCostFunction( new XtalResidual(s.x, s.y, lambda, pixel, cos_rot3, sin_rot3, 0.0, - h, k, l, system)), + h, k, l, system, orientation)), new ceres::CauchyLoss(loss_scale), beam, &distance_mm, detector_rot, rot_vec, orient[i].data(), cell_len, cell_ang); ++used_in_frame; diff --git a/image_analysis/geom_refinement/PostRefine.cpp b/image_analysis/geom_refinement/PostRefine.cpp index 945da29e8..c16755fa0 100644 --- a/image_analysis/geom_refinement/PostRefine.cpp +++ b/image_analysis/geom_refinement/PostRefine.cpp @@ -329,6 +329,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vectorh, pp->k, pp->l, s)) continue; XtalResidual r(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3, angle_rad(pp->img), - pp->h, pp->k, pp->l, sys); + pp->h, pp->k, pp->l, sys, orientation); double resid[3] = {0, 0, 0}; r(beam, dist, det_rot, rot_vec, p0, p1, p2, resid); c += resid[0]*resid[0] + resid[1]*resid[1] + resid[2]*resid[2]; ++n; @@ -683,7 +684,7 @@ PostRefineResult PostRefineRotationGeometry(const std::vector( new XtalResidualBeamDistance( XtalResidual(pp->obs_x, pp->obs_y, lambda_l, pixel_mm, cos_rot3, sin_rot3, - angle_rad(pp->img), pp->h, pp->k, pp->l, sys), + angle_rad(pp->img), pp->h, pp->k, pp->l, sys, orientation), fc, p0)), new ceres::CauchyLoss(0.02), beam, dist); } diff --git a/image_analysis/geom_refinement/RingOptimizer.cpp b/image_analysis/geom_refinement/RingOptimizer.cpp index 1e845d245..f29d32ed9 100644 --- a/image_analysis/geom_refinement/RingOptimizer.cpp +++ b/image_analysis/geom_refinement/RingOptimizer.cpp @@ -3,6 +3,7 @@ #include +#include "../../common/DetectorOrientation.h" #include "../../common/JFJochMath.h" #include "RingOptimizer.h" #include "ceres/ceres.h" @@ -10,11 +11,18 @@ struct RingResidual { RingResidual(double x, double y, double lambda, double pixel_size, - double expected_q) + double expected_q, + const DetectorOrientation &orientation) : obs_x(x), obs_y(y), lambda(lambda), pixel_size(pixel_size), - expected_len_recip_sq(expected_q * expected_q / (4.0 * PI * PI)) {} + expected_len_recip_sq(expected_q * expected_q / (4.0 * PI * PI)) { + const RotMatrix delta = orientation.Matrix(); + det_m00 = delta.Column(0).x; + det_m01 = delta.Column(1).x; + det_m10 = delta.Column(0).y; + det_m11 = delta.Column(1).y; + } template bool operator()(const T* const center_x, const T* const center_y, @@ -22,8 +30,13 @@ struct RingResidual { const T* const rot2, T* residual) const { // Calculate lab coordinates from observed pixel coordinates - T x_lab = (T(obs_x) - center_x[0]) * T(pixel_size); // convert to mm - T y_lab = (T(obs_y) - center_y[0]) * T(pixel_size); + T u_lab = (T(obs_x) - center_x[0]) * T(pixel_size); // convert to mm + T v_lab = (T(obs_y) - center_y[0]) * T(pixel_size); + // The discrete image orientation, which turns the offset from the PONI before the tilt acts. + // Identity unless a detector says otherwise. It cannot change a ring's radius, but it does + // change which way the tilt tips the ring, which is exactly what this fits. + T x_lab = det_m00 * u_lab + det_m01 * v_lab; + T y_lab = det_m10 * u_lab + det_m11 * v_lab; T z_lab = distance[0]; // Apply rotations around y and x axes @@ -53,6 +66,7 @@ struct RingResidual { const double lambda; const double pixel_size; const double expected_len_recip_sq; + double det_m00, det_m01, det_m10, det_m11; }; RingOptimizer::RingOptimizer(const DiffractionGeometry& geom) : reference(geom) {} @@ -74,7 +88,8 @@ DiffractionGeometry RingOptimizer::Run(const std::vector &in new RingResidual(pt.x, pt.y, reference.GetWavelength_A(), reference.GetPixelSize_mm(), - pt.q_expected)), + pt.q_expected, + reference.GetOrientation())), nullptr, ¢er_x, ¢er_y, diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 53e7d15eb..5397a5afa 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -282,7 +282,8 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, cos_rot3, sin_rot3, angle_rad, h, k, l, - data.crystal_system); + data.crystal_system, + data.geom.GetOrientation()); // Ceres has no per-residual weight; ScaledLoss(nullptr, a) multiplies the squared // residual by the constant a, i.e. it applies a weight of sqrt(a) to the residual. diff --git a/image_analysis/geom_refinement/XtalResidual.h b/image_analysis/geom_refinement/XtalResidual.h index 5e49aa1df..3d206210f 100644 --- a/image_analysis/geom_refinement/XtalResidual.h +++ b/image_analysis/geom_refinement/XtalResidual.h @@ -13,6 +13,7 @@ #include "../../common/JFJochException.h" #include "../../common/CrystalLattice.h" +#include "../../common/DetectorOrientation.h" // Rodrigues rotation with everything that depends only on the angle-axis worked out once. This is // ceres::AngleAxisRotatePoint term for term - so the rotated point is bit-identical to it - but that @@ -92,7 +93,8 @@ struct XtalResidual { double angle_rad, double exp_h, double exp_k, double exp_l, - gemmi::CrystalSystem symmetry) + gemmi::CrystalSystem symmetry, + const DetectorOrientation &orientation = {}) : obs_x(x), obs_y(y), inv_lambda(1.0/lambda), pixel_size(pixel_size), @@ -106,6 +108,12 @@ struct XtalResidual { if (std::fabs(lambda) < 1e-6) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Lambda cannot be close to zero"); + + const RotMatrix delta = orientation.Matrix(); + det_m00 = delta.Column(0).x; + det_m01 = delta.Column(1).x; + det_m10 = delta.Column(0).y; + det_m11 = delta.Column(1).y; } // The observed reciprocal vector: the spot's detector position taken through the current beam, @@ -159,9 +167,14 @@ struct XtalResidual { const double c3 = cos_rot3; const double s3 = sin_rot3; - // Detector coordinates in mm - const T det_x = (T(obs_x) - beam[0]) * T(pixel_size); - const T det_y = (T(obs_y) - beam[1]) * T(pixel_size); + // Detector coordinates in mm, then the discrete image orientation, which turns the offset + // from the PONI before the tilt acts on it. It cannot be folded into rot3: rot3 turns the + // whole detector in the laboratory, this turns the image within the detector plane, and + // Ry(rot1) and Rx(-rot2) sit between them. Identity unless a detector says otherwise. + const T det_u = (T(obs_x) - beam[0]) * T(pixel_size); + const T det_v = (T(obs_y) - beam[1]) * T(pixel_size); + const T det_x = det_m00 * det_u + det_m01 * det_v; + const T det_y = det_m10 * det_u + det_m11 * det_v; const D &det_z = distance_mm; // Apply Ry(rot1) first: rotate around Y @@ -338,6 +351,8 @@ struct XtalResidual { const double inv_lambda; const double pixel_size; const double cos_rot3, sin_rot3; + // The discrete image orientation's action on (u, v). Identity by default. + double det_m00, det_m01, det_m10, det_m11; const double exp_h; const double exp_k; const double exp_l; diff --git a/reader/HDF5MetadataSource.cpp b/reader/HDF5MetadataSource.cpp index fe4889994..a99e92f5c 100644 --- a/reader/HDF5MetadataSource.cpp +++ b/reader/HDF5MetadataSource.cpp @@ -818,6 +818,16 @@ HDF5MetadataSource::OpenResult HDF5MetadataSource::Open(const std::string &filen detector.MirrorY(master_file ->GetOptBool("/entry/instrument/detector/detectorSpecific/mirror_y") .value_or(true)); + // How the stored image sits in the detector plane. A different setting from mirror_y above, + // recorded separately by the writer; absence means the identity, which is what a file written + // before it existed - or by anything else - describes. + detector.ImageOrientation(DetectorOrientation( + master_file->GetOptBool( + "/entry/instrument/detector/detectorSpecific/detector_orientation_mirror_y") + .value_or(false), + master_file->GetOptInt( + "/entry/instrument/detector/detectorSpecific/detector_orientation_quarter_turns") + .value_or(0))); // Sensor thickness/material drive the parallax/absorption model, so take them from the file // rather than the DetectorSetup default (NXmx stores thickness in metres). if (master_file->Exists("/entry/instrument/detector/sensor_thickness")) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 20a916b99..2a99832a6 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -202,6 +202,8 @@ void print_usage() { std::cout << " --rot1 PONI rotation 1 (rad)" << std::endl; std::cout << " --rot2 PONI rotation 2 (rad)" << std::endl; std::cout << " --rot3 PONI rotation 3, about the beam (rad)" << std::endl; + std::cout << " --detector-mirror-y Stored image is mirrored in Y vs the detector" << std::endl; + std::cout << " --detector-quarter-turns <0-3> Stored image is turned by this many 90 deg about the beam" << std::endl; std::cout << " --polarization Polarization factor" << std::endl; } @@ -286,6 +288,8 @@ enum { OPT_ROT1, OPT_ROT2, OPT_ROT3, + OPT_DETECTOR_MIRROR_Y, + OPT_DETECTOR_QUARTER_TURNS, OPT_FFT_MIN_UNIT_CELL, OPT_POLARIZATION }; @@ -338,6 +342,8 @@ static option long_options[] = { {"rot1", required_argument, nullptr, OPT_ROT1}, {"rot2", required_argument, nullptr, OPT_ROT2}, {"rot3", required_argument, nullptr, OPT_ROT3}, + {"detector-mirror-y", no_argument, nullptr, OPT_DETECTOR_MIRROR_Y}, + {"detector-quarter-turns", required_argument, nullptr, OPT_DETECTOR_QUARTER_TURNS}, {"fft-min-unit-cell", required_argument, nullptr, OPT_FFT_MIN_UNIT_CELL}, {"polarization", required_argument, nullptr, OPT_POLARIZATION}, {"force-rotation-lattice", required_argument, nullptr, OPT_FORCE_ROTATION_LATTICE}, @@ -671,6 +677,8 @@ static int RunRugnux(int argc, char **argv) { // Geometry overrides (default: keep the value stored in the input file) std::optional beam_x, beam_y, detector_distance_mm, wavelength_A, rot1_rad, rot2_rad, rot3_rad, polarization_factor; + bool detector_mirror_y = false; + int64_t detector_quarter_turns = 0; std::optional smooth_g_deg_arg; // --smooth-g[=deg]; default 5 deg for rot3d, 0 (off) otherwise std::optional relative_b_deg_arg; // --relative-b[=deg]; per-batch relative-B width, 0 (off) unless given bool no_scaling_corrections = false; // --no-scaling-corrections: disable rot3d decay+absorption+modulation surfaces @@ -1239,6 +1247,14 @@ static int RunRugnux(int argc, char **argv) { case OPT_ROT1: rot1_rad = parse_float_arg(optarg, "--rot1", logger); break; case OPT_ROT2: rot2_rad = parse_float_arg(optarg, "--rot2", logger); break; case OPT_ROT3: rot3_rad = parse_float_arg(optarg, "--rot3", logger); break; + case OPT_DETECTOR_MIRROR_Y: detector_mirror_y = true; break; + case OPT_DETECTOR_QUARTER_TURNS: + detector_quarter_turns = atoi(optarg); + if (detector_quarter_turns < 0 || detector_quarter_turns > 3) { + logger.Error("--detector-quarter-turns must be 0, 1, 2 or 3"); + return 1; + } + break; case OPT_FFT_MIN_UNIT_CELL: fft_min_unit_cell_A = parse_float_arg(optarg, "--fft-min-unit-cell", logger); break; case OPT_POLARIZATION: polarization_factor = parse_float_arg(optarg, "--polarization", logger); break; @@ -1830,6 +1846,12 @@ static int RunRugnux(int argc, char **argv) { if (rot1_rad) experiment.PoniRot1_rad(rot1_rad.value()); if (rot2_rad) experiment.PoniRot2_rad(rot2_rad.value()); if (rot3_rad) experiment.PoniRot3_rad(rot3_rad.value()); + // The discrete part of the detector orientation, which the file does not state for a detector this + // system did not assemble. Never derived from rot3 - it says how the image is stored, and a fitted + // angle must not change that. + if (detector_mirror_y || detector_quarter_turns != 0) + experiment.Detector().ImageOrientation( + DetectorOrientation(detector_mirror_y, detector_quarter_turns)); // --polarization is applied after configure_offline_output below, which sets the rugnux default. // Calibrating from the run-summed profile needs the profile to be binned in azimuth - see diff --git a/tests/DiffractionGeometryTest.cpp b/tests/DiffractionGeometryTest.cpp index cb9da7d67..81f7651f5 100644 --- a/tests/DiffractionGeometryTest.cpp +++ b/tests/DiffractionGeometryTest.cpp @@ -6,6 +6,7 @@ #include #include "../common/DiffractionGeometry.h" #include "../common/DiffractionExperiment.h" +#include "../common/JFJochMath.h" TEST_CASE("RecipToDetector_1", "[LinearAlgebra][Coord]") { DiffractionExperiment x(DetJF(8, 2)); @@ -527,7 +528,7 @@ TEST_CASE("DiffractionGeometry_PONI_matrix_consistency") { .PixelSize_mm(0.075).Wavelength_A(1.0) .PoniRot1_rad(0.04).PoniRot2_rad(-0.025); - const auto& poni_rot = geom.GetPoniRotMatrix(); + const auto& poni_rot = geom.GetDetectorMatrix(); const auto poni_rot_T = poni_rot.transpose(); // Test: poni_rot * poni_rot^T should be identity (orthogonal matrix) @@ -631,3 +632,202 @@ TEST_CASE("DiffractionGeometry_Tilted_vs_PyFAI_and_DIALS", "[DiffractionGeometry CHECK(lab.z == Catch::Approx(-r.dials[2]).margin(margin)); } } + +// --------------------------------------------------------------------------------------------- +// PONI angles <-> detector axis vectors, and the discrete image orientation +// --------------------------------------------------------------------------------------------- + +namespace { + void CheckSameMatrix(const RotMatrix &a, const RotMatrix &b, float margin = 1e-6f) { + for (int i = 0; i < 3; i++) { + const Coord ca = a.Column(i), cb = b.Column(i); + CHECK(ca.x == Catch::Approx(cb.x).margin(margin)); + CHECK(ca.y == Catch::Approx(cb.y).margin(margin)); + CHECK(ca.z == Catch::Approx(cb.z).margin(margin)); + } + } +} + +TEST_CASE("PoniAngles_matrix_roundtrip") { + const float half_pi = static_cast(PI) / 2.0f; + // rot1 and rot3 are recovered by atan2, so the branch cut at +-pi makes an angle comparison there + // meaningless (+pi and -pi are the same rotation). The matrix comparison below covers it; the + // angle comparison uses everything else, including the exact multiples of 90 degrees that are not + // on the cut. + const std::vector angles = {0.0f, 0.01f, -0.03f, 0.7f, -1.2f, half_pi, -half_pi}; + + for (float rot1: angles) { + for (float rot3: angles) { + for (float rot2: {0.0f, 0.02f, -0.4f, 1.0f, -1.4f}) { + float r1, r2, r3; + PoniAnglesFromMatrix(PoniRotMatrix(rot1, rot2, rot3), r1, r2, r3); + CHECK(r1 == Catch::Approx(rot1).margin(1e-5)); + CHECK(r2 == Catch::Approx(rot2).margin(1e-5)); + CHECK(r3 == Catch::Approx(rot3).margin(1e-5)); + CheckSameMatrix(PoniRotMatrix(r1, r2, r3), PoniRotMatrix(rot1, rot2, rot3)); + } + } + } + + // A half turn is on the atan2 branch cut, so only the matrix can be required to come back. + for (float rot1: {static_cast(PI), -static_cast(PI)}) { + float r1, r2, r3; + PoniAnglesFromMatrix(PoniRotMatrix(rot1, 0.1f, 0.2f), r1, r2, r3); + CheckSameMatrix(PoniRotMatrix(r1, r2, r3), PoniRotMatrix(rot1, 0.1f, 0.2f)); + } + + // Gimbal lock: at rot2 = +-90 degrees only rot1 +- rot3 is determined, and the convention is to + // put it all into rot1. A triple that already has rot3 = 0 therefore comes back unchanged, and + // the matrix comes back whatever rot3 was. + for (float rot2: {half_pi, -half_pi}) { + for (float rot1: {0.0f, 0.3f, -1.2f}) { + float r1, r2, r3; + PoniAnglesFromMatrix(PoniRotMatrix(rot1, rot2, 0.0f), r1, r2, r3); + CHECK(r1 == Catch::Approx(rot1).margin(1e-5)); + CHECK(r2 == Catch::Approx(rot2).margin(1e-5)); + CHECK(r3 == 0.0f); + + PoniAnglesFromMatrix(PoniRotMatrix(rot1, rot2, 0.4f), r1, r2, r3); + CheckSameMatrix(PoniRotMatrix(r1, r2, r3), PoniRotMatrix(rot1, rot2, 0.4f)); + } + } +} + +TEST_CASE("DetectorAxes_roundtrip") { + for (int64_t quarter_turns = 0; quarter_turns < 4; quarter_turns++) { + for (bool mirror: {false, true}) { + for (float rot1: {0.0f, 0.05f, -0.9f}) { + for (float rot2: {0.0f, -0.03f, 1.1f}) { + for (float rot3: {0.0f, 0.2f, -1.5f}) { + DiffractionGeometry geom; + geom.Orientation(DetectorOrientation(mirror, quarter_turns)) + .PoniRot1_rad(rot1).PoniRot2_rad(rot2).PoniRot3_rad(rot3); + + const Coord fast = geom.GetFastAxis(); + const Coord slow = geom.GetSlowAxis(); + const RotMatrix before = geom.GetDetectorMatrix(); + + // Feeding the two axes straight back must not move anything. + DiffractionGeometry from_axes; + from_axes.Orientation(DetectorOrientation(mirror, quarter_turns)) + .DetectorAxes(fast, slow); + CheckSameMatrix(from_axes.GetDetectorMatrix(), before, 1e-5f); + CHECK(from_axes.GetPoniRot1_rad() == Catch::Approx(rot1).margin(1e-5)); + CHECK(from_axes.GetPoniRot2_rad() == Catch::Approx(rot2).margin(1e-5)); + CHECK(from_axes.GetPoniRot3_rad() == Catch::Approx(rot3).margin(1e-5)); + } + } + } + } + } +} + +TEST_CASE("DetectorOrientation_identity_is_todays_geometry") { + DiffractionGeometry with_default; + with_default.PoniRot1_rad(0.04f).PoniRot2_rad(-0.02f).PoniRot3_rad(0.11f); + + DiffractionGeometry with_identity; + with_identity.Orientation(DetectorOrientation(false, 0)) + .PoniRot1_rad(0.04f).PoniRot2_rad(-0.02f).PoniRot3_rad(0.11f); + + // Bit for bit: the discrete part must cost existing data nothing. + CheckSameMatrix(with_identity.GetDetectorMatrix(), with_default.GetDetectorMatrix(), 0.0f); + CheckSameMatrix(with_default.GetDetectorMatrix(), PoniRotMatrix(0.04f, -0.02f, 0.11f), 0.0f); + CHECK(with_default.GetOrientation().IsIdentity()); +} + +TEST_CASE("DetectorOrientation_maps_the_detector_plane") { + // Untilted, so the lab coordinate of a pixel is the discrete orientation applied to its offset + // from the PONI, in mm. + auto make = [](bool mirror, int64_t quarter_turns) { + DiffractionGeometry g; + g.BeamX_pxl(100).BeamY_pxl(200).DetectorDistance_mm(100).PixelSize_mm(0.1f) + .Orientation(DetectorOrientation(mirror, quarter_turns)); + return g; + }; + + // One pixel along the fast direction is 0.1 mm from the PONI. + const Coord fast_step = make(false, 0).LabCoord(101, 200) - make(false, 0).LabCoord(100, 200); + CHECK(fast_step.x == Catch::Approx(0.1).margin(1e-6)); + CHECK(fast_step.y == Catch::Approx(0.0).margin(1e-6)); + + // A quarter turn about the beam takes the fast direction to lab +y ... + CHECK(make(false, 1).GetFastAxis().y == Catch::Approx(1.0).margin(1e-6)); + // ... and the slow direction to lab -x. + CHECK(make(false, 1).GetSlowAxis().x == Catch::Approx(-1.0).margin(1e-6)); + // A mirror in Y leaves the fast direction alone and reverses the slow one. + CHECK(make(true, 0).GetFastAxis().x == Catch::Approx(1.0).margin(1e-6)); + CHECK(make(true, 0).GetSlowAxis().y == Catch::Approx(-1.0).margin(1e-6)); + // Two quarter turns is a half turn. + CHECK(make(false, 2).GetFastAxis().x == Catch::Approx(-1.0).margin(1e-6)); + CHECK(make(false, 2).GetSlowAxis().y == Catch::Approx(-1.0).margin(1e-6)); + + // Every orientation is orthogonal, and improper exactly when it mirrors. + for (int64_t k = 0; k < 4; k++) + for (bool mirror: {false, true}) { + const DetectorOrientation o(mirror, k); + CheckSameMatrix(o.Matrix() * o.Matrix().transpose(), RotMatrix(), 1e-6f); + const Coord expected_normal = mirror ? -(o.Matrix().Column(0) % o.Matrix().Column(1)) + : (o.Matrix().Column(0) % o.Matrix().Column(1)); + CheckSameMatrix(RotMatrix(o.Matrix().Column(0), o.Matrix().Column(1), expected_normal), + o.Matrix(), 1e-6f); + } +} + +TEST_CASE("DetectorOrientation_preserves_radius_and_solid_angle") { + // Both generators are signed permutations of (u, v), so the distance from the PONI - and with it + // the solid-angle correction, the resolution of a ring and every radius-only consumer - cannot + // move. This is why most of the pipeline needs no change. + const float ref = [] { + DiffractionGeometry g; + g.BeamX_pxl(500).BeamY_pxl(700).DetectorDistance_mm(120).PixelSize_mm(0.075f); + return g.CalcAzIntSolidAngleCorr(823, 311); + }(); + + for (int64_t k = 0; k < 4; k++) + for (bool mirror: {false, true}) { + DiffractionGeometry g; + g.BeamX_pxl(500).BeamY_pxl(700).DetectorDistance_mm(120).PixelSize_mm(0.075f) + .PoniRot1_rad(0.03f).PoniRot2_rad(-0.02f) + .Orientation(DetectorOrientation(mirror, k)); + CHECK(g.CalcAzIntSolidAngleCorr(823, 311) == Catch::Approx(ref).margin(1e-7)); + } +} + +TEST_CASE("DetectorOrientation_and_polarization") { + // Polarization depends on the azimuth in the LABORATORY, so what the discrete orientation changes + // is which pixel lands where. A quarter turn moves a pixel from the polarization plane to across + // it; a mirror in Y sends phi to -phi and so cannot move it at all. + auto corr = [](bool mirror, int64_t quarter_turns, float x, float y) { + DiffractionGeometry g; + g.BeamX_pxl(500).BeamY_pxl(500).DetectorDistance_mm(100).PixelSize_mm(0.075f) + .Orientation(DetectorOrientation(mirror, quarter_turns)); + return g.CalcAzIntPolarizationCorr(x, y, 0.99f); + }; + + const float along_x = corr(false, 0, 700, 500); + const float along_y = corr(false, 0, 500, 700); + CHECK(along_x != Catch::Approx(along_y)); + + CHECK(corr(false, 1, 700, 500) == Catch::Approx(along_y)); + CHECK(corr(true, 0, 700, 500) == Catch::Approx(along_x)); + CHECK(corr(true, 0, 500, 700) == Catch::Approx(along_y)); +} + +TEST_CASE("DetectorOrientation_recip_roundtrip") { + for (int64_t k = 0; k < 4; k++) + for (bool mirror: {false, true}) { + DiffractionGeometry geom; + geom.BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150) + .PixelSize_mm(0.075f).Wavelength_A(1.0f) + .PoniRot1_rad(0.05f).PoniRot2_rad(-0.03f).PoniRot3_rad(0.2f) + .Orientation(DetectorOrientation(mirror, k)); + + for (const auto &[x, y]: std::vector>{ + {500, 500}, {1500, 500}, {500, 1500}, {1200, 800}}) { + const auto [px, py] = geom.RecipToDetector(geom.DetectorToRecip(x, y)); + CHECK(px == Catch::Approx(x).margin(0.001)); + CHECK(py == Catch::Approx(y).margin(0.001)); + } + } +} diff --git a/viewer/JFJochHttpReader.cpp b/viewer/JFJochHttpReader.cpp index f0111062b..1c0bf8c81 100644 --- a/viewer/JFJochHttpReader.cpp +++ b/viewer/JFJochHttpReader.cpp @@ -232,6 +232,8 @@ std::shared_ptr JFJochHttpReader::UpdateDataset_i() { DetectorSetup detector = DetDECTRIS(msg->start_message->image_size_x, msg->start_message->image_size_y, msg->start_message->detector_description, {}); detector.PixelSize_um(msg->start_message->pixel_size_x * 1e6); + detector.ImageOrientation(DetectorOrientation(msg->start_message->detector_orientation_mirror_y, + msg->start_message->detector_orientation_quarter_turns)); detector.SaturationLimit(SaturationLimitFromValue(msg->start_message->saturation_value)); detector.MinFrameTime(std::chrono::microseconds(0)); detector.MinCountTime(std::chrono::microseconds(0)); diff --git a/writer/HDF5NXmx.cpp b/writer/HDF5NXmx.cpp index 52a2ae81f..89a5be7ac 100644 --- a/writer/HDF5NXmx.cpp +++ b/writer/HDF5NXmx.cpp @@ -6,6 +6,7 @@ #include "HDF5NXmx.h" +#include "../common/DetectorOrientation.h" #include "../common/GitInfo.h" #include "../include/spdlog/fmt/fmt.h" #include "MakeDirectory.h" @@ -427,6 +428,11 @@ void NXmx::Detector(const StartMessage &start) { // the geometry is but not how it was arrived at. This records the assembly setting itself, so a // re-opened file knows whether the stored image was mirrored rather than having to infer it. SaveScalar(det_specific, "mirror_y", start.mirror_y); + // Likewise for the discrete image orientation: the module axis vectors below carry its effect, + // these two carry the setting. + SaveScalar(det_specific, "detector_orientation_mirror_y", start.detector_orientation_mirror_y); + SaveScalar(det_specific, "detector_orientation_quarter_turns", + start.detector_orientation_quarter_turns); det_specific.NXClass("NXcollection"); if (!start.jfjoch_release.empty()) @@ -651,9 +657,22 @@ void NXmx::Metrology(const StartMessage &start, const EndMessage &end) { HDF5Group transformations(*hdf5_file, "/entry/instrument/detector/transformations"); transformations.NXClass("NXtransformations"); - std::vector vector{beam_center_x * start.pixel_size_x, - beam_center_y * start.pixel_size_y, - start.detector_distance}; + // Internal frame (x = column, y = row downward, z = beam) -> McStas, a half turn about z. Written + // as a subtraction from zero rather than a negation so that a zero component stays a positive + // zero, and an untilted, unturned detector writes the same axis vectors it always has. + auto to_mcstas = [](const Coord &v) { + return std::vector{0.0 - v.x, 0.0 - v.y, 0.0 + v.z}; + }; + + // The discrete image orientation turns the offset from the PONI before the rot1/rot2/rot3 chain + // acts, so it belongs to the module axes and to the translation, not to the arm rotations. + const RotMatrix delta = DetectorOrientation(start.detector_orientation_mirror_y, + start.detector_orientation_quarter_turns).Matrix(); + + // Sample -> module origin (pixel 0, 0), which is where LabCoord(0, 0) puts it. + std::vector vector = to_mcstas(delta * Coord(-beam_center_x * start.pixel_size_x, + -beam_center_y * start.pixel_size_y, + start.detector_distance)); double vector_length = sqrt(vector[0] * vector[0] + vector[1] * vector[1] + vector[2] * vector[2]); std::vector vector_norm{vector[0] / vector_length, vector[1]/vector_length, vector[2]/vector_length}; @@ -705,8 +724,9 @@ void NXmx::Metrology(const StartMessage &start, const EndMessage &end) { "rotation", std::vector{0.0, 0.0, -1.0}); - DetectorModule("module", origin, size, {-1,0,0}, {0,-1,0}, "translation", - start.pixel_size_x); + DetectorModule("module", origin, size, + to_mcstas(delta * Coord(1, 0, 0)), to_mcstas(delta * Coord(0, 1, 0)), + "translation", start.pixel_size_x); } void SaveUnitCell( HDF5Group& group, const std::string& name, const UnitCell& unit_cell) {