diff --git a/image_analysis/geom_refinement/CMakeLists.txt b/image_analysis/geom_refinement/CMakeLists.txt index e751d11f..164a96d4 100644 --- a/image_analysis/geom_refinement/CMakeLists.txt +++ b/image_analysis/geom_refinement/CMakeLists.txt @@ -4,6 +4,8 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC RingOptimizer.h AssignSpotsToRings.cpp AssignSpotsToRings.h + RingsFromProfile.cpp + RingsFromProfile.h XtalOptimizer.cpp XtalOptimizer.h XtalResidual.h diff --git a/image_analysis/geom_refinement/RingsFromProfile.cpp b/image_analysis/geom_refinement/RingsFromProfile.cpp new file mode 100644 index 00000000..c9372845 --- /dev/null +++ b/image_analysis/geom_refinement/RingsFromProfile.cpp @@ -0,0 +1,125 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include "RingsFromProfile.h" +#include "AssignSpotsToRings.h" // CalculateXtalRings +#include "../../common/JFJochMath.h" + +namespace { + +// Peak position of one ring in one azimuthal sector, in q, or NaN if there is no peak worth using. +// +// The window is narrow and centred on where the ring is expected, so the background under it is close +// to a straight line: take it from the two bins at each end and interpolate. The position itself is the +// intensity-weighted centroid of everything above half the peak height, which is insensitive to the +// exact half-maximum crossing and needs no line-shape assumption - a powder ring is not Gaussian, it is +// the instrumental profile convolved with whatever strain and size broadening the standard has. +float SectorPeakQ(const std::vector &profile, int32_t q_bins, int phi_bin, + int lo_bin, int hi_bin, float low_q, float q_spacing, float min_peak_over_noise) { + const size_t row = static_cast(phi_bin) * static_cast(q_bins); + const auto value = [&](int i) { return profile[row + static_cast(i)]; }; + const auto q_of = [&](int i) { return low_q + (static_cast(i) + 0.5f) * q_spacing; }; + + const float bkg_lo = 0.5f * (value(lo_bin) + value(lo_bin + 1)); + const float bkg_hi = 0.5f * (value(hi_bin) + value(hi_bin - 1)); + const auto bkg_at = [&](int i) { + const float t = static_cast(i - lo_bin) / static_cast(hi_bin - lo_bin); + return bkg_lo + t * (bkg_hi - bkg_lo); + }; + + int peak = -1; + float peak_height = 0.0f; + for (int i = lo_bin + 2; i <= hi_bin - 2; ++i) { + const float h = value(i) - bkg_at(i); + if (h > peak_height) { peak_height = h; peak = i; } + } + if (peak < 0) + return NAN; + + // Scatter of the background shoulders, as the noise this peak has to stand clear of. A sector with + // no ring in it has a "peak" that is just the largest background fluctuation, and this is what + // rejects it - the alternative, an absolute intensity cut, would need a value per detector and beam. + float s = 0.0f; + int n = 0; + for (int i : {lo_bin, lo_bin + 1, hi_bin - 1, hi_bin}) { + const float r = value(i) - bkg_at(i); + s += r * r; + ++n; + } + const float noise = std::sqrt(s / static_cast(n)); + if (!(peak_height > min_peak_over_noise * noise)) + return NAN; + + const float half = 0.5f * peak_height; + double sum_wq = 0.0, sum_w = 0.0; + for (int i = peak; i >= lo_bin && value(i) - bkg_at(i) >= half; --i) { + const double w = value(i) - bkg_at(i); + sum_wq += w * q_of(i); + sum_w += w; + } + for (int i = peak + 1; i <= hi_bin && value(i) - bkg_at(i) >= half; ++i) { + const double w = value(i) - bkg_at(i); + sum_wq += w * q_of(i); + sum_w += w; + } + if (!(sum_w > 0.0)) + return NAN; + return static_cast(sum_wq / sum_w); +} + +} // namespace + +std::vector RingsFromAzimuthalProfile(const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + const UnitCell &calibrant, + float q_window_recipA, + float min_peak_over_noise) { + std::vector out; + + const int32_t q_bins = mapping.GetQBinCount(); + const int32_t azim_bins = mapping.GetAzimuthalBinCount(); + // One azimuthal bin is a plain radial profile: the ring is averaged over every direction at once, so + // nothing remains to say where its centre is. This needs the run to have been integrated with + // azimuthal bins (jfjoch_broker azim_int_settings.azimuthal_bins, rugnux --azim-phi-bins). + if (azim_bins < 4 || q_bins < 8 + || profile.size() != static_cast(q_bins) * static_cast(azim_bins)) + return out; + + const auto &settings = mapping.Settings(); + const float low_q = settings.GetLowQ_recipA(); + const float q_spacing = settings.GetQSpacing_recipA(); + const float high_q = low_q + static_cast(q_bins) * q_spacing; + const int window_bins = std::max(3, static_cast(std::lround(q_window_recipA / q_spacing))); + + for (const float q_ring : CalculateXtalRings(calibrant)) { + if (!(q_ring - q_window_recipA > low_q) || !(q_ring + q_window_recipA < high_q)) + continue; + const int centre_bin = static_cast((q_ring - low_q) / q_spacing); + const int lo_bin = std::max(0, centre_bin - window_bins); + const int hi_bin = std::min(q_bins - 1, centre_bin + window_bins); + if (hi_bin - lo_bin < 6) + continue; + + for (int phi_bin = 0; phi_bin < azim_bins; ++phi_bin) { + const float q_obs = SectorPeakQ(profile, q_bins, phi_bin, lo_bin, hi_bin, + low_q, q_spacing, min_peak_over_noise); + if (!std::isfinite(q_obs)) + continue; + + // The sector's CENTRE, not its lower edge: GetBin() floors phi into the sector, so a bin + // stands for [j, j+1) and taking its edge would rotate every ring point by half a sector - + // which is exactly the cos(phi) signal the beam centre is read from. + const float phi_rad = static_cast((static_cast(phi_bin) + 0.5) + * 2.0 * PI / static_cast(azim_bins)); + const auto [x, y] = geom.ResPhiToPxl(static_cast(2.0 * PI) / q_obs, phi_rad); + if (!std::isfinite(x) || !std::isfinite(y)) + continue; + out.push_back({x, y, q_ring}); + } + } + return out; +} diff --git a/image_analysis/geom_refinement/RingsFromProfile.h b/image_analysis/geom_refinement/RingsFromProfile.h new file mode 100644 index 00000000..9be6b30a --- /dev/null +++ b/image_analysis/geom_refinement/RingsFromProfile.h @@ -0,0 +1,37 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../../common/AzimuthalIntegrationMapping.h" +#include "../../common/DiffractionGeometry.h" +#include "../../common/UnitCell.h" +#include "RingOptimizer.h" + +// Turn an accumulated (q x azimuth) powder profile into ring points for RingOptimizer. +// +// The calibration this feeds already exists (AssignSpotsToRings + RingOptimizer); what it has always +// been given is a SPOT LIST from a single image. A powder ring is not a set of spots - it is a smooth +// arc - so a spot finder samples it wherever its threshold happens to bite, and one image carries only +// as much of the ring as that image's counting statistics allow. An azimuthally-binned profile summed +// over a run measures the same ring directly, at every azimuth, with the whole run's counts behind it. +// +// Where the ring falls is what carries the geometry. A powder ring is a conic centred on the beam, so +// if the beam centre is wrong its apparent radius oscillates once per turn (a cos(phi) term), and if +// the detector is tilted, twice (cos(2 phi)). Neither depends on the calibrant's d-spacings, which is +// why the beam centre is the one thing a powder pattern determines without assuming anything about the +// standard - the distance, by contrast, is only as good as the lattice constant it is measured against. +// +// profile is the mean intensity per bin (AzimuthalIntegrationProfile::GetResult()): q_bins x azimuthal +// bins, indexed bin = q_bin + phi_bin * q_bins. geom supplies the CURRENT geometry, used only to turn a +// measured (q, phi) back into the pixel it came from - RingOptimizer then refines that geometry so the +// q it predicts at that pixel matches the calibrant's. Rings outside the profile's q range, and sectors +// where no peak stands clear of the local background, are skipped rather than guessed at. +std::vector RingsFromAzimuthalProfile(const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + const UnitCell &calibrant, + float q_window_recipA = 0.06f, + float min_peak_over_noise = 3.0f); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 494b3ad5..7e04b3f1 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -59,6 +59,7 @@ ADD_EXECUTABLE(jfjoch_test JFJochReceiverPlotsTest.cpp GoniometerAxisTest.cpp DetGeomCalibTest.cpp + RingsFromProfileTest.cpp XtalOptimizerTest.cpp CrystalLatticeTest.cpp FPGAPTPTest.cpp diff --git a/tests/RingsFromProfileTest.cpp b/tests/RingsFromProfileTest.cpp new file mode 100644 index 00000000..433961c8 --- /dev/null +++ b/tests/RingsFromProfileTest.cpp @@ -0,0 +1,123 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include "../image_analysis/geom_refinement/RingsFromProfile.h" +#include "../image_analysis/geom_refinement/AssignSpotsToRings.h" +#include "../common/Definitions.h" +#include "../common/JFJochMath.h" + +namespace { + +constexpr UnitCell LAB6{LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90.0f, 90.0f, 90.0f}; + +// A (q x azimuth) powder profile as the azimuthal integration would build it: the rings sit where +// geom_true puts them, but every pixel is binned with geom_assumed - which is the whole point, since a +// wrong assumed geometry is what makes a ring's apparent q wander with azimuth. +std::vector SynthesiseProfile(const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom_assumed, + const DiffractionGeometry &geom_true) { + const auto &settings = mapping.Settings(); + const int32_t q_bins = mapping.GetQBinCount(); + const int32_t azim_bins = mapping.GetAzimuthalBinCount(); + std::vector profile(static_cast(q_bins) * azim_bins, 100.0f); // flat background + + for (const float q_ring : CalculateXtalRings(LAB6)) { + const float d = static_cast(2.0 * PI) / q_ring; + if (d <= geom_true.GetWavelength_A() / 2.0f) + continue; + for (int t = 0; t < 3600; ++t) { + const float phi_true = static_cast(2.0 * PI * t / 3600.0); + const auto [px, py] = geom_true.ResPhiToPxl(d, phi_true); + if (!std::isfinite(px) || !std::isfinite(py)) + continue; + const float q_obs = geom_assumed.PxlToQ(px, py); + float phi_deg = geom_assumed.Phi_rad(px, py) * 180.0f / static_cast(PI); + if (phi_deg < 0.0f) + phi_deg += 360.0f; + const uint16_t bin = settings.GetBin(q_obs, phi_deg); + if (bin == UINT16_MAX) + continue; + // Lay a narrow peak over the neighbouring q bins of this azimuthal row. + const int q_bin = bin % q_bins, phi_bin = bin / q_bins; + for (int k = -3; k <= 3; ++k) { + const int b = q_bin + k; + if (b < 0 || b >= q_bins) + continue; + profile[static_cast(phi_bin) * q_bins + b] += + 2000.0f * std::exp(-0.5f * static_cast(k * k) / (1.2f * 1.2f)); + } + } + } + return profile; +} + +} // namespace + +// The measurement this is for: a powder ring is a conic centred on the beam, so a wrong beam centre +// makes its apparent radius oscillate once per turn. Recovering the centre from that needs neither the +// calibrant's lattice constant nor the detector distance - only that the ring be round. +TEST_CASE("RingsFromProfile_RecoversBeamCenter", "[DetGeomCalib]") { + DiffractionExperiment x(DetJF4M()); + x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0); + auto azint = x.GetAzimuthalIntegrationSettings(); + azint.AzimuthalBinCount(32); + x.ImportAzimuthalIntegrationSettings(azint); + + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + const DiffractionGeometry geom_assumed = x.GetDiffractionGeometry(); + DiffractionGeometry geom_true = geom_assumed; + geom_true.BeamX_pxl(geom_assumed.GetBeamX_pxl() + 6.0f) + .BeamY_pxl(geom_assumed.GetBeamY_pxl() - 4.0f); + + const auto profile = SynthesiseProfile(mapping, geom_assumed, geom_true); + const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom_assumed, LAB6); + + // Several rings, sampled all the way round: without azimuthal coverage there is no centre to find. + REQUIRE(rings.size() > 64); + + RingOptimizer optimizer(geom_assumed); + const auto fitted = optimizer.Run(rings); + + CHECK(fitted.GetBeamX_pxl() == Catch::Approx(geom_true.GetBeamX_pxl()).margin(0.5)); + CHECK(fitted.GetBeamY_pxl() == Catch::Approx(geom_true.GetBeamY_pxl()).margin(0.5)); + // The starting point was wrong by 6 and 4 pixels, so a fit that did nothing would fail the above - + // but check explicitly that it moved toward the truth rather than merely landing near it. + CHECK(std::abs(fitted.GetBeamX_pxl() - geom_true.GetBeamX_pxl()) + < std::abs(geom_assumed.GetBeamX_pxl() - geom_true.GetBeamX_pxl())); +} + +// One azimuthal bin is a plain radial profile: the ring has been averaged over every direction, so +// nothing is left to say where its centre is. Refuse rather than return points that cannot constrain it. +TEST_CASE("RingsFromProfile_NeedsAzimuthalBins", "[DetGeomCalib]") { + DiffractionExperiment x(DetJF4M()); + x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0); + + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + REQUIRE(mapping.GetAzimuthalBinCount() == 1); + + const std::vector profile(static_cast(mapping.GetQBinCount()), 1000.0f); + CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6).empty()); +} + +// A profile with no rings in it must yield no ring points: the peak has to stand clear of the scatter +// of the background either side of it, or every azimuthal sector would contribute its largest noise +// excursion as though it were a measurement. +TEST_CASE("RingsFromProfile_FlatProfileGivesNothing", "[DetGeomCalib]") { + DiffractionExperiment x(DetJF4M()); + x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0); + auto azint = x.GetAzimuthalIntegrationSettings(); + azint.AzimuthalBinCount(32); + x.ImportAzimuthalIntegrationSettings(azint); + + PixelMask pixel_mask(x); + AzimuthalIntegrationMapping mapping(x, pixel_mask); + + const std::vector profile( + static_cast(mapping.GetQBinCount()) * mapping.GetAzimuthalBinCount(), 100.0f); + CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6).empty()); +}