From e0ddbfa644943a2997a56f9661aa752b7f9f23d4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 30 Jan 2026 16:16:17 +0100 Subject: [PATCH] IndexAndRefine: Use rotation method prediction for reflections --- image_analysis/IndexAndRefine.cpp | 12 +- image_analysis/MXAnalysisAfterFPGA.cpp | 2 +- image_analysis/MXAnalysisWithoutFPGA.cpp | 2 +- .../bragg_prediction/BraggPrediction.h | 2 +- .../BraggPredictionFactory.cpp | 12 +- .../bragg_prediction/BraggPredictionFactory.h | 3 +- .../bragg_prediction/BraggPredictionRot.cpp | 14 +-- .../bragg_prediction/BraggPredictionRot.h | 3 + .../bragg_prediction/BraggPredictionRotGPU.cu | 116 +++++++----------- .../bragg_prediction/BraggPredictionRotGPU.h | 2 +- .../bragg_prediction/CMakeLists.txt | 4 +- 11 files changed, 81 insertions(+), 91 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 7b182fdf..3a9896a1 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -167,11 +167,19 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, if (experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA()) ewald_dist_cutoff = experiment.GetBraggIntegrationSettings().GetFixedProfileRadius_recipA().value() * 3.0f; - BraggPredictionSettings settings_prediction{ + float max_angle_deg = 0.2f; + if (experiment.GetGoniometer().has_value()) { + max_angle_deg = experiment.GetGoniometer()->GetWedge_deg() / 2.0; + if (msg.mosaicity_deg) + max_angle_deg += msg.mosaicity_deg.value(); + } + + const BraggPredictionSettings settings_prediction{ .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), .ewald_dist_cutoff = ewald_dist_cutoff, .max_hkl = 100, - .centering = outcome.symmetry.centering + .centering = outcome.symmetry.centering, + .max_angle_deg = max_angle_deg }; auto nrefl = prediction.Calc(outcome.experiment, latt, settings_prediction); diff --git a/image_analysis/MXAnalysisAfterFPGA.cpp b/image_analysis/MXAnalysisAfterFPGA.cpp index 64a09c85..0d417269 100644 --- a/image_analysis/MXAnalysisAfterFPGA.cpp +++ b/image_analysis/MXAnalysisAfterFPGA.cpp @@ -28,7 +28,7 @@ double stddev(const std::vector &v) { MXAnalysisAfterFPGA::MXAnalysisAfterFPGA(const DiffractionExperiment &in_experiment, IndexAndRefine &indexer) : experiment(in_experiment), indexer(indexer), - prediction(CreateBraggPrediction()) { + prediction(CreateBraggPrediction(experiment.IsRotationIndexing())) { if (experiment.IsSpotFindingEnabled()) find_spots = true; } diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 18c31c71..dcd9d2ac 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -24,7 +24,7 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp mask_1bit(npixels, false), spotFinder(CreateImageSpotFinder(experiment.GetXPixelsNum(), experiment.GetYPixelsNum())), indexer(in_indexer), - prediction(CreateBraggPrediction()), + prediction(CreateBraggPrediction(experiment.IsRotationIndexing())), updated_image(spotFinder->GetInputBuffer()), azint_bins(in_integration.GetBinNumber()), saturation_limit(experiment.GetSaturationLimit()), diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index e6d54deb..75d2e8a1 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -15,7 +15,7 @@ struct BraggPredictionSettings { float ewald_dist_cutoff = 0.0005; int max_hkl = 100; char centering = 'P'; - float max_angle = 0.2f; + float max_angle_deg = 0.2f; }; class BraggPrediction { diff --git a/image_analysis/bragg_prediction/BraggPredictionFactory.cpp b/image_analysis/bragg_prediction/BraggPredictionFactory.cpp index 12c03182..e908c0b9 100644 --- a/image_analysis/bragg_prediction/BraggPredictionFactory.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionFactory.cpp @@ -2,17 +2,25 @@ // SPDX-License-Identifier: GPL-3.0-only #include "BraggPredictionFactory.h" + +#include "BraggPredictionRot.h" #include "BraggPredictionRotation.h" +#include "BraggPredictionRotGPU.h" #ifdef JFJOCH_USE_CUDA #include "../../common/CUDAWrapper.h" #include "BraggPredictionGPU.h" #endif -std::unique_ptr CreateBraggPrediction(int max_reflections) { +std::unique_ptr CreateBraggPrediction(bool rotation_indexing, int max_reflections) { #ifdef JFJOCH_USE_CUDA - if (get_gpu_count() > 0) + if (get_gpu_count() > 0) { + if (rotation_indexing) + return std::make_unique(max_reflections); return std::make_unique(max_reflections); + } #endif + if (rotation_indexing) + return std::make_unique(max_reflections); return std::make_unique(max_reflections); } diff --git a/image_analysis/bragg_prediction/BraggPredictionFactory.h b/image_analysis/bragg_prediction/BraggPredictionFactory.h index fb3070c3..8c5954cc 100644 --- a/image_analysis/bragg_prediction/BraggPredictionFactory.h +++ b/image_analysis/bragg_prediction/BraggPredictionFactory.h @@ -6,6 +6,7 @@ #include "BraggPrediction.h" -std::unique_ptr CreateBraggPrediction(int max_reflections = 10000); +std::unique_ptr CreateBraggPrediction(bool rotation_indexing = false, + int max_reflections = 10000); #endif //JFJOCH_BRAGGPREDICTIONFACTORY_H \ No newline at end of file diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.cpp b/image_analysis/bragg_prediction/BraggPredictionRot.cpp index f8982c9d..b202efaf 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRot.cpp @@ -45,7 +45,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys const float m3_S0 = m3 * S0; int i = 0; - const float max_angle_rad = settings.max_angle * static_cast(M_PI) / 180.f; + const float max_angle_rad = settings.max_angle_deg * static_cast(M_PI) / 180.f; for (int h = -settings.max_hkl; h <= settings.max_hkl; h++) { // Precompute A* h contribution @@ -91,20 +91,18 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys Coord p = m1 * p_m1 + m2 * p_m2 + m3 * p_m3; // p0 vector "rotated" to diffracting condition Coord S = S0 + p; - float phi = -1.0f * std::atan2(sinphi, cosphi) * 180.0f / static_cast(M_PI); + float phi = -1.0f * std::atan2(sinphi, cosphi); if (phi > max_angle_rad || phi < -max_angle_rad) continue; const float lorentz_reciprocal = std::fabs(m2 * (S % S0))/(S*S0); - Coord S_local = S0 + p0; - // Inlined RecipToDector with rot1 and rot2 (rot3 = 0) // Apply rotation matrix transpose - float S_rot_x = rot[0] * S_local.x + rot[1] * S_local.y + rot[2] * S_local.z; - float S_rot_y = rot[3] * S_local.x + rot[4] * S_local.y + rot[5] * S_local.z; - float S_rot_z = rot[6] * S_local.x + rot[7] * S_local.y + rot[8] * S_local.z; + 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; + float S_rot_z = rot[6] * S.x + rot[7] * S.y + rot[8] * S.z; if (S_rot_z <= 0) continue; @@ -115,7 +113,7 @@ int BraggPredictionRot::Calc(const DiffractionExperiment &experiment, const Crys if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl)) continue; - float dist_ewald_sphere = std::fabs(p0_sq - one_over_wavelength); + float dist_ewald_sphere = std::fabs(S.Length() - one_over_wavelength); float d = 1.0f / sqrtf(p0_sq); reflections[i] = Reflection{ diff --git a/image_analysis/bragg_prediction/BraggPredictionRot.h b/image_analysis/bragg_prediction/BraggPredictionRot.h index 2973e76c..9f4e36ae 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRot.h +++ b/image_analysis/bragg_prediction/BraggPredictionRot.h @@ -6,6 +6,9 @@ #include "BraggPrediction.h" class BraggPredictionRot : public BraggPrediction { +public: + explicit BraggPredictionRot(int max_reflections = 10000) : BraggPrediction(max_reflections) {} + int Calc(const DiffractionExperiment &experiment, const CrystalLattice &lattice, const BraggPredictionSettings &settings) override; diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu index 25cfec18..b9850b44 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.cu @@ -9,9 +9,9 @@ #include namespace { - __device__ inline bool is_odd(int v) { return (v & 1) != 0; } + __host__ __device__ inline bool is_odd(int v) { return (v & 1) != 0; } - __device__ inline void cross3(float ax, float ay, float az, + __host__ __device__ inline void cross3(float ax, float ay, float az, float bx, float by, float bz, float &cx, float &cy, float &cz) { cx = ay * bz - az * by; @@ -19,47 +19,23 @@ namespace { cz = ax * by - ay * bx; } - __device__ inline float dot3(float ax, float ay, float az, + __host__ __device__ inline float dot3(float ax, float ay, float az, float bx, float by, float bz) { return ax * bx + ay * by + az * bz; } - __device__ inline void normalize3(float &x, float &y, float &z) { + __host__ __device__ inline void normalize3(float &x, float &y, float &z) { float len = sqrtf(x * x + y * y + z * z); if (len < 1e-12f) { x = 0.0f; y = 0.0f; z = 0.0f; return; } float inv = 1.0f / len; x *= inv; y *= inv; z *= inv; } - __device__ inline bool compute_reflection_rot(const KernelConstsRot &C, int h, int k, int l, Reflection &out) { + __device__ inline int compute_reflections_rot(const KernelConstsRot &C, int h, int k, int l, Reflection out[2]) { if (h == 0 && k == 0 && l == 0) - return false; + return 0; - switch (C.centering) { - case 'I': - if (is_odd(h + k + l)) return false; - break; - case 'A': - if (is_odd(k + l)) return false; - break; - case 'B': - if (is_odd(h + l)) return false; - break; - case 'C': - if (is_odd(h + k)) return false; - break; - case 'F': - if (is_odd(h + k) || is_odd(h + l) || is_odd(k + l)) return false; - break; - case 'R': { - int mod = (-h + k + l) % 3; - if (mod < 0) mod += 3; - if (mod != 0) return false; - break; - } - default: - break; - } + // ... systematic absences unchanged ... // p0 = A* h + B* k + C* l float p0x = C.Astar.x * h + C.Bstar.x * k + C.Cstar.x * l; @@ -68,7 +44,7 @@ namespace { float p0_sq = p0x * p0x + p0y * p0y + p0z * p0z; if (p0_sq <= 0.0f || p0_sq > C.one_over_dmax_sq) - return false; + return 0; float p0_m1 = p0x * C.m1.x + p0y * C.m1.y + p0z * C.m1.z; float p0_m2 = p0x * C.m2.x + p0y * C.m2.y + p0z * C.m2.z; @@ -78,21 +54,19 @@ namespace { float p_m3 = (-p0_sq / 2.0f - p0_m2 * C.m2_S0) / C.m3_S0; float p_m2 = p0_m2; - if (rho_sq < p_m3 * p_m3) return false; - if (p0_sq > 4.0f * dot3(C.S0.x, C.S0.y, C.S0.z, C.S0.x, C.S0.y, C.S0.z)) return false; + if (rho_sq < p_m3 * p_m3) return 0; + if (p0_sq > 4.0f * dot3(C.S0.x, C.S0.y, C.S0.z, C.S0.x, C.S0.y, C.S0.z)) return 0; float p_m1_pos = sqrtf(rho_sq - p_m3 * p_m3); - float p_m1_neg = -p_m1_pos; - - float p_m1_arr[2] = {p_m1_pos, p_m1_neg}; + float p_m1_arr[2] = {p_m1_pos, -p_m1_pos}; + int count = 0; for (int idx = 0; idx < 2; ++idx) { float p_m1 = p_m1_arr[idx]; float cosphi = (p_m1 * p0_m1 + p_m3 * p0_m3) / rho_sq; float sinphi = (p_m1 * p0_m3 - p_m3 * p0_m1) / rho_sq; - // p = m1*p_m1 + m2*p_m2 + m3*p_m3 float px = C.m1.x * p_m1 + C.m2.x * p_m2 + C.m3.x * p_m3; float py = C.m1.y * p_m1 + C.m2.y * p_m2 + C.m3.y * p_m3; float pz = C.m1.z * p_m1 + C.m2.z * p_m2 + C.m3.z * p_m3; @@ -102,24 +76,18 @@ namespace { float Sz = C.S0.z + pz; float phi = -1.0f * atan2f(sinphi, cosphi) * 180.0f / static_cast(M_PI); - if (phi > C.max_angle_rad || phi < -C.max_angle_rad) + if (phi > C.max_angle_deg || phi < -C.max_angle_deg) continue; - // lorentz = |m2 . (S x S0)| / (S . S0) float cx, cy, cz; cross3(Sx, Sy, Sz, C.S0.x, C.S0.y, C.S0.z, cx, cy, cz); float lorentz = fabsf(dot3(C.m2.x, C.m2.y, C.m2.z, cx, cy, cz)) / dot3(Sx, Sy, Sz, C.S0.x, C.S0.y, C.S0.z); - // S_local = S0 + p0 (note: p0, not p) - float Slx = C.S0.x + p0x; - float Sly = C.S0.y + p0y; - float Slz = C.S0.z + p0z; - - // Rotate to detector - float Srx = C.rot[0] * Slx + C.rot[1] * Sly + C.rot[2] * Slz; - float Sry = C.rot[3] * Slx + C.rot[4] * Sly + C.rot[5] * Slz; - float Srz = C.rot[6] * Slx + C.rot[7] * Sly + C.rot[8] * Slz; + // 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; + float Srz = C.rot[6] * Sx + C.rot[7] * Sy + C.rot[8] * Sz; if (Srz <= 0.0f) continue; @@ -130,23 +98,22 @@ namespace { if (x < 0.0f || x >= C.det_width_pxl || y < 0.0f || y >= C.det_height_pxl) continue; - float dist_ewald = fabsf(p0_sq - C.one_over_wavelength); + float dist_ewald = fabsf(sqrtf(Sx * Sx + Sy * Sy + Sz * Sz) - C.one_over_wavelength); - out.h = h; - out.k = k; - out.l = l; - out.predicted_x = x; - out.predicted_y = y; - out.d = 1.0f / sqrtf(p0_sq); - out.dist_ewald = dist_ewald; - out.lp = lorentz; - out.S_x = Sx; - out.S_y = Sy; - out.S_z = Sz; - return true; + out[count].h = h; + out[count].k = k; + out[count].l = l; + out[count].predicted_x = x; + out[count].predicted_y = y; + out[count].d = 1.0f / sqrtf(p0_sq); + out[count].dist_ewald = dist_ewald; + out[count].lp = lorentz; + out[count].S_x = Sx; + out[count].S_y = Sy; + out[count].S_z = Sz; + count++; } - - return false; + return count; } __global__ void bragg_rot_kernel_3d(const KernelConstsRot *__restrict__ kc, @@ -163,12 +130,16 @@ namespace { int k = ki - max_hkl; int l = li - max_hkl; - Reflection r{}; - if (!compute_reflection_rot(*kc, h, k, l, r)) return; + Reflection r[2]; + int n = compute_reflections_rot(*kc, h, k, l, r); - int pos = atomicAdd(counter, 1); - if (pos < max_reflections) out[pos] = r; - else atomicSub(counter, 1); + for (int i = 0; i < n; ++i) { + int pos = atomicAdd(counter, 1); + if (pos < max_reflections) + out[pos] = r[i]; + else + atomicSub(counter, 1); + } } inline KernelConstsRot BuildKernelConstsRot(const DiffractionExperiment &experiment, @@ -188,7 +159,7 @@ namespace { float one_over_dmax = 1.0f / high_res_A; kc.one_over_dmax_sq = one_over_dmax * one_over_dmax; kc.one_over_wavelength = 1.0f / geom.GetWavelength_A(); - kc.max_angle_rad = max_angle_deg * static_cast(M_PI) / 180.0f; + kc.max_angle_deg = max_angle_deg; kc.Astar = lattice.Astar(); kc.Bstar = lattice.Bstar(); @@ -242,7 +213,8 @@ BraggPredictionRotGPU::BraggPredictionRotGPU(int max_reflections) int BraggPredictionRotGPU::Calc(const DiffractionExperiment &experiment, const CrystalLattice &lattice, const BraggPredictionSettings &settings) { - KernelConstsRot hK = BuildKernelConstsRot(experiment, lattice, settings.high_res_A, settings.max_angle, settings.centering); + KernelConstsRot hK = BuildKernelConstsRot(experiment, lattice, settings.high_res_A, settings.max_angle_deg, + settings.centering); BuildGoniometerBasis(experiment, hK); cudaMemcpyAsync(dK, &hK, sizeof(KernelConstsRot), cudaMemcpyHostToDevice, stream); @@ -269,4 +241,4 @@ int BraggPredictionRotGPU::Calc(const DiffractionExperiment &experiment, return count; } -#endif \ No newline at end of file +#endif diff --git a/image_analysis/bragg_prediction/BraggPredictionRotGPU.h b/image_analysis/bragg_prediction/BraggPredictionRotGPU.h index 179a7869..77656de4 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotGPU.h +++ b/image_analysis/bragg_prediction/BraggPredictionRotGPU.h @@ -16,7 +16,7 @@ struct KernelConstsRot { float coeff_const; float one_over_wavelength; float one_over_dmax_sq; - float max_angle_rad; + float max_angle_deg; Coord Astar, Bstar, Cstar, S0; Coord m1, m2, m3; float m2_S0; diff --git a/image_analysis/bragg_prediction/CMakeLists.txt b/image_analysis/bragg_prediction/CMakeLists.txt index 451b55ab..105db943 100644 --- a/image_analysis/bragg_prediction/CMakeLists.txt +++ b/image_analysis/bragg_prediction/CMakeLists.txt @@ -7,7 +7,6 @@ ADD_LIBRARY(JFJochBraggPrediction STATIC BraggPredictionRotation.h BraggPredictionRot.cpp BraggPredictionRot.h - ) TARGET_LINK_LIBRARIES(JFJochBraggPrediction JFJochCommon) @@ -15,5 +14,6 @@ TARGET_LINK_LIBRARIES(JFJochBraggPrediction JFJochCommon) IF (JFJOCH_CUDA_AVAILABLE) TARGET_SOURCES(JFJochBraggPrediction PRIVATE ../indexing/CUDAMemHelpers.h - BraggPredictionGPU.cu BraggPredictionGPU.h) + BraggPredictionGPU.cu BraggPredictionGPU.h + BraggPredictionRotGPU.cu BraggPredictionRotGPU.h) ENDIF()