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"