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) {