From 2746c50822d8c7f11aee666dee44ebf1745b4bb7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sat, 26 Sep 2026 22:54:37 +0200 Subject: [PATCH] BeamCenterFromBackground: skip pixels clearly outside the band before sqrt and atan2 60-70% of a large detector lies outside the resolution band the walk bins, and every one paid a square root and two arctangents on every iteration. A test on tan(2theta) = rho / lz with a 0.1% margin, far above float rounding, drops them first; every pixel the exact test keeps still reaches it. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- .../geom_refinement/BeamCenterFromBackground.cpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp index 9129589bc..b2c306b9e 100644 --- a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp @@ -127,6 +127,12 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe std::vector block_jy(static_cast(BLOCKS) * n_cells); std::vector block_count(static_cast(BLOCKS) * n_cells); + // Most pixels lie outside the band. Those clearly outside it in tan(2theta) = rho / lz - by a + // margin far above float rounding - skip the square root and both atan2 below; every pixel the exact + // test would keep still reaches it. + const double tan_lo = std::tan(static_cast(tt_lo)) * (1.0 - 1e-3); + const double tan_hi = tt_hi < PI / 2 ? std::tan(static_cast(tt_hi)) * (1.0 + 1e-3) : INFINITY; + float step_x = 0.0f, step_y = 0.0f, sigma_x = 0.0f, sigma_y = 0.0f; float previous_x = 0.0f, previous_y = 0.0f; int reversals = 0; @@ -153,7 +159,13 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe const float lx = rot[0] * u + rot[1] * v + rot[2] * distance; const float ly = rot[3] * u + rot[4] * v + rot[5] * distance; const float lz = rot[6] * u + rot[7] * v + rot[8] * distance; - const float rho = std::sqrt(lx * lx + ly * ly); + const float rho_sq = lx * lx + ly * ly; + if (lz > 0.0f) { + const double lz_sq = static_cast(lz) * lz; + if (rho_sq < tan_lo * tan_lo * lz_sq || rho_sq > tan_hi * tan_hi * lz_sq) + continue; + } + const float rho = std::sqrt(rho_sq); const float two_theta = std::atan2(rho, lz); if (two_theta < tt_lo || two_theta >= tt_hi || rho == 0.0f) continue;