From e7ef72ba407578cca7fa1c2178f243c2bc2c0b81 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 24 Sep 2026 08:48:31 +0200 Subject: [PATCH] Rotation prediction: place each partial at its slice of the rocking curve A rotation frame records the part of a reflection's rocking curve inside its own oscillation, and while the crystal turns through the curve the spot walks along its Debye ring. The predictor put every partial at the exact diffracting condition, so a partial recorded on the curve's flank was integrated pixels away from where its flux landed. The walk is largest where the reflection moves nearly tangent to the Ewald sphere (low |zeta|), whose curves are widest. The prediction now turns S about the beam by the rotation's component along the ring times the flux-weighted centre of the frame's slice of the curve (a truncated-normal mean, RockingSlice.h), on the CPU and the GPU predictor alike. 2theta is unchanged; a frame that straddles the condition symmetrically gets no shift. Measured (observed r1 centroid minus prediction, tangential, I/sig > 10): against the slice model correlation 0.93-0.98 before, 0.0 after; residual rms 1.0-1.8 px -> 0.35-0.48 px at |zeta| < 0.3. Emulated capture of low-|zeta| partials at high angle 0.50 -> 0.94 of a 12 px aperture. CPU runs on four rotation sets (400-800 frames): R_meas -0.03..-0.09 pp, ISa +1.2..+1.8, rejected observations roughly halved, space groups unchanged. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- .../bragg_prediction/BraggPredictionRot.cpp | 12 +++ .../bragg_prediction/BraggPredictionRotGPU.cu | 18 +++++ .../bragg_prediction/RockingSlice.h | 42 ++++++++++ tests/CalcBraggPredictionTest.cpp | 76 +++++++++++++++++++ 4 files changed, 148 insertions(+) create mode 100644 image_analysis/bragg_prediction/RockingSlice.h diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.cpp b/image_analysis/bragg_prediction/BraggPredictionRot.cpp index 1d4321d05..495a40a6c 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRot.cpp @@ -5,6 +5,7 @@ #include "BraggPredictionRot.h" #include "../SensorAbsorption.h" #include "../bragg_integration/SystematicAbsence.h" +#include "RockingSlice.h" int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const CrystalLattice &lattice, @@ -47,6 +48,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys const float m3_S0 = m3 * S0; int i = 0; + const Coord beam_dir = S0.Normalize(); const float mos_angle_rad = settings.mosaicity_deg * static_cast(PI) / 180.f; const float half_wedge_angle_rad = settings.wedge_deg * static_cast(PI) / 180.f / 2.0f ; @@ -163,6 +165,16 @@ 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; + // Place the partial where this frame recorded it: at the centre of the frame's slice of + // the rocking curve, which is the exact-condition position walked along its Debye ring by + // the rotation between the two (RockingSlice.h). Walking along the ring is turning S + // about the beam, by the rotation's component along the ring over the ring's radius. + const float c_slice = RockingSliceCentroid_rad(phi, half_wedge_angle_rad, c1, partiality); + const Coord S_par = beam_dir * (beam_dir * S), S_perp = S - S_par; + const Coord S_turn = beam_dir % S; // |S_turn| = |S_perp| + const float psi = c_slice * ((m2 % p) * S_turn) / (S_perp * S_perp); + S = S_par + S_perp * std::cos(psi) + S_turn * std::sin(psi); + // 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; diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index 59855cf14..07ce14209 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -4,6 +4,7 @@ #include #include "../../common/JFJochMath.h" #include "BraggPredictionRotGPU.h" +#include "RockingSlice.h" #include "../SensorAbsorption.h" #ifdef JFJOCH_USE_CUDA @@ -156,6 +157,23 @@ namespace { - erff((phi - half_wedge) * c1)) / 2.0f; + // The partial sits at the centre of the frame's slice of its rocking curve, walked along + // the Debye ring from the exact condition: S turned about the beam (see the CPU predictor). + float kx = C.S0.x, ky = C.S0.y, kz = C.S0.z; + normalize3(kx, ky, kz); + const float kS = dot3(kx, ky, kz, Sx, Sy, Sz); + const float qx = Sx - kS * kx, qy = Sy - kS * ky, qz = Sz - kS * kz; // S_perp + float tx, ty, tz; // S_turn = k x S + cross3(kx, ky, kz, Sx, Sy, Sz, tx, ty, tz); + float wx, wy, wz; // m2 x p + cross3(C.m2.x, C.m2.y, C.m2.z, px, py, pz, wx, wy, wz); + const float c_slice = RockingSliceCentroid_rad(phi, half_wedge, c1, partiality); + const float psi = c_slice * dot3(wx, wy, wz, tx, ty, tz) / dot3(qx, qy, qz, qx, qy, qz); + const float cos_psi = cosf(psi), sin_psi = sinf(psi); + Sx = kS * kx + qx * cos_psi + tx * sin_psi; + Sy = kS * ky + qy * cos_psi + ty * sin_psi; + Sz = kS * kz + qz * cos_psi + tz * sin_psi; + // Use S (rotated) for projection float Srx = C.rot[0] * Sx + C.rot[1] * Sy + C.rot[2] * Sz; float Sry = C.rot[3] * Sx + C.rot[4] * Sy + C.rot[5] * Sz; diff --git a/image_analysis/bragg_prediction/RockingSlice.h b/image_analysis/bragg_prediction/RockingSlice.h new file mode 100644 index 000000000..6178f6bac --- /dev/null +++ b/image_analysis/bragg_prediction/RockingSlice.h @@ -0,0 +1,42 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// Where on its rocking curve a rotation frame records a reflection. +// +// The rocking curve is a Gaussian in the rotation angle, centred on the exact diffracting condition, +// and a frame records the part of it inside its own oscillation: [phi - w/2, phi + w/2] in the +// convention of the rotation predictor, where phi is the frame's offset from that condition and +// c1 = zeta / (sqrt(2) sigma_M) is the curve's reciprocal width. The flux-weighted centre of that slice +// is the mean of a truncated normal, +// +// c = (exp(-a^2) - exp(-b^2)) / (2 sqrt(pi) c1 p), a = (phi - w/2) c1, b = (phi + w/2) c1, +// +// with p the slice's mass - the partiality. A frame whose oscillation straddles the condition +// symmetrically gets c = 0; one that only catches the curve's tail gets a c of up to several curve +// widths. +// +// It matters because the recorded spot is not fixed on the detector while the crystal rotates through +// the curve: it walks along its Debye ring, one detector velocity per radian of rotation, and a partial +// recorded at c sits c times that velocity from the exact-condition position. Measured on rotation +// data the walk follows this c with slope ~1 (0.9 at high angle, where the r1 centroid it was measured +// with is truncated least), and at |zeta| < 0.2 it is pixels - the reflection moves nearly tangent to +// the Ewald sphere, so its curve is wide and the spot travels far along the ring while it lasts. + +#include + +#ifdef __CUDACC__ +#define ROCKING_SLICE_HD __host__ __device__ inline +#else +#define ROCKING_SLICE_HD inline +#endif + +// Flux-weighted centre of the frame's slice of the rocking curve, relative to the exact diffracting +// condition, in radians of rotation. partiality is the slice's own mass (the caller has it already). +ROCKING_SLICE_HD float RockingSliceCentroid_rad(float phi, float half_wedge, float c1, float partiality) { + if (!(partiality > 0.0f) || !(c1 > 0.0f)) + return 0.0f; + const float a = (phi - half_wedge) * c1, b = (phi + half_wedge) * c1; + return (expf(-a * a) - expf(-b * b)) / (2.0f * sqrtf(3.14159265f) * c1 * partiality); +} diff --git a/tests/CalcBraggPredictionTest.cpp b/tests/CalcBraggPredictionTest.cpp index 752823cf6..9d94c9a90 100644 --- a/tests/CalcBraggPredictionTest.cpp +++ b/tests/CalcBraggPredictionTest.cpp @@ -6,6 +6,9 @@ #include "../image_analysis/bragg_prediction/BraggPrediction.h" #include #include "../image_analysis/SensorAbsorption.h" +#include "../image_analysis/bragg_prediction/BraggPredictionRot.h" +#include "../image_analysis/bragg_prediction/RockingSlice.h" +#include // The flight-path term is a zero-parameter prediction: a NIST attenuation coefficient, the stated // sample-to-detector distance, and Beer-Lambert. Nothing about it is fitted, so it can be checked @@ -282,6 +285,79 @@ TEST_CASE("BraggPrediction_systematic_absences") { } } +TEST_CASE("RockingSliceCentroid_TruncatedNormalMean", "[rocking_slice]") { + // The closed form against a direct quadrature of the slice, on a frame that holds the centre, one on + // the curve's flank and one on its far tail. + const float sigma = 0.004f, half_wedge = 0.0015f, c1 = 1.0f / (std::sqrt(2.0f) * sigma); + for (float phi : {0.0f, 0.001f, -0.006f, 0.012f}) { + double sw = 0.0, stw = 0.0; + for (int i = 0; i <= 20000; ++i) { + const double t = phi - half_wedge + 2.0 * half_wedge * i / 20000.0; + const double w = std::exp(-t * t / (2.0 * sigma * sigma)); + sw += w; stw += t * w; + } + const float partiality = (std::erf((phi + half_wedge) * c1) - std::erf((phi - half_wedge) * c1)) / 2.0f; + INFO("phi " << phi); + CHECK(RockingSliceCentroid_rad(phi, half_wedge, c1, partiality) == Catch::Approx(stw / sw).margin(2e-6)); + } + CHECK(RockingSliceCentroid_rad(0.0f, half_wedge, c1, 0.3f) == 0.0f); +} + +TEST_CASE("BraggPredictionRot_PartialWalksAlongItsRing", "[rocking_slice]") { + // Each rotation frame predicts a partial where the frame's slice of its rocking curve puts it: the + // exact-condition position walked along the Debye ring. Over the frames of one reflection the + // predicted positions therefore stay at one distance from the beam and move monotonically along + // the ring; a prediction at the exact condition would be the same point on every frame. + DiffractionExperiment experiment(DetJF4M()); + experiment.DetectorDistance_mm(100.0).BeamX_pxl(1100.0).BeamY_pxl(1000.0).IncidentEnergy_keV(12.4); + const GoniometerAxis axis("omega", 0.0f, 0.1f, Coord(-1, 0, 0), {}); + experiment.Goniometer(axis); + const CrystalLattice lattice(Coord{40, 0, 0}, Coord{0, 50, 0}, Coord{0, 0, 60}); + BraggPredictionSettings settings{.high_res_A = 2.5, .ewald_dist_cutoff = 0.0015, + .max_h = 20, .max_k = 25, .max_l = 30, + .wedge_deg = 0.1f, .mosaicity_deg = 0.1f}; + const float bx = 1100.0f, by = 1000.0f; + + std::map, std::vector>> track; // (radius, azimuth) + BraggPredictionRot pred; + for (int frame = 0; frame < 300; ++frame) { + const auto latt = lattice.Multiply(axis.GetTransformationAngle(frame * 0.1f)); + const int n = pred.Calc(experiment, latt, settings); + for (int i = 0; i < n; ++i) { + const auto &r = pred.GetReflections().at(i); + if (r.zeta > 0.3f) continue; + const float dx = r.predicted_x - bx, dy = r.predicted_y - by; + track[{r.h, r.k, r.l}].emplace_back(std::hypot(dx, dy), std::atan2(dy, dx)); + } + } + int tested = 0; + for (const auto &[hkl, pts] : track) { + if (pts.size() < 20) continue; + float rmin = pts[0].first, rmax = pts[0].first; + float along_first = 0.0f, along_last = 0.0f; + int sign_changes = 0; + float prev_step = 0.0f; + for (size_t j = 0; j < pts.size(); ++j) { + rmin = std::min(rmin, pts[j].first); + rmax = std::max(rmax, pts[j].first); + const float along = (pts[j].second - pts[0].second) * pts[0].first; // px along the ring + if (j == 0) along_first = along; + along_last = along; + if (j > 0) { + const float step = (pts[j].second - pts[j - 1].second); + if (step * prev_step < 0.0f) ++sign_changes; + if (step != 0.0f) prev_step = step; + } + } + INFO("hkl " << std::get<0>(hkl) << " " << std::get<1>(hkl) << " " << std::get<2>(hkl) << " frames " << pts.size()); + CHECK(rmax - rmin < 0.05f); // stays on its ring + CHECK(std::fabs(along_last - along_first) > 0.1f); // but walks along it + CHECK(sign_changes == 0); // in one direction + ++tested; + } + CHECK(tested > 20); +} + #ifdef JFJOCH_USE_CUDA #include "../image_analysis/bragg_prediction/BraggPredictionGPU.h" #include "../image_analysis/bragg_prediction/BraggPredictionRotGPU.h"