From a5f416fcdc324893d72dadf133e4af49e9cb37fa Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 31 Aug 2026 16:48:32 +0200 Subject: [PATCH] calibration: take the detector distance from the rings, not from the header A powder calibration is run because nobody is sure the header is right, and the header's distance was the one number the fit could not survive being wrong about. The ring search is local - each ring is looked for inside a window a few pixels of radius wide - so a distance more than a percent or two out puts every ring outside its own window, and the fit then converges on whatever background fluctuation each window contains. It does not fail: a 110 mm exposure told the detector was at 150 mm reported 149.8 mm, with 146 ring points and exit 0. Only its residual said anything, 5.3 px against 0.4 px, and nothing read it. Measure the distance from the rings instead. The peaks of the azimuthally averaged profile give ring RADII, and a radius does not depend on the assumed distance at all - bin i holds the pixels at one particular radius whatever q that radius was called - so the radii are a property of the image. Against the calibrant's d-spacings, r = D tan(2 asin(lambda/2d)) then has one unknown. It is scanned rather than solved because the pairing of observed rings to d-spacings is unknown too, and the winning basin is solved in closed form. Nothing here reads the header distance except to bin the profile; it needs only the wavelength, the pixel size and the detector's extent. A powder pattern has genuine distance aliases, so one answer is not enough. A cubic primitive standard puts its rings at radii proportional to sqrt(N), and scaling the distance by sqrt(2) maps ring N onto ring 2N - most of the comb still lands on peaks. Measured: the 110 mm exposure with a 115 mm header scored its best at 156.5 mm, which is 110*sqrt(2). No adjustment of the score removes an alias the lattice really has, so the scan hands back the few best distances and each is fitted, the header among them as one hypothesis of several. The residual then separates them - 0.4 px against 5.2 px on that case - subject to an attempt explaining a comparable share of the pattern first, because a start so wrong that one ring point survives leaves a residual of exactly zero. Each attempt re-extracts at the geometry it converged to and fits again. The seed is measured from blended peaks and is good to about a per cent, close enough to converge from but far enough to sit every search window a few pixels off its ring, and an off-centre window takes its background off the ring's own flank. Nothing is re-read from disk, so the loop is free. Measured on the LaB6 distance series. A 110 mm dataset now recovers 110.03-110.17 mm from any header between 25 and 1200 mm, against +-2 mm before. All five datasets recover their own distance from a fixed wrong 250 mm header. With correct headers, four of the five are bit-identical to before and the 500 mm one moves by a single ring point - the two-ring fit whose tilt is 0.1 sigma anyway. Run time is unchanged at 0.62 s. The residual is larger on a run whose header was wrong (1.1 px against 0.4 px on the 110 mm case), because the profile was still binned at the wrong distance and its radial sampling is correspondingly coarse. The geometry is right; only the scatter about it is inflated. Re-running with the recovered distance recovers the residual too. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01NfuDvf5ipV3Hi8TiCUKD27 --- image_analysis/geom_refinement/CMakeLists.txt | 2 + .../geom_refinement/PowderAutoSeed.cpp | 307 ++++++++++++++++++ .../geom_refinement/PowderAutoSeed.h | 84 +++++ .../geom_refinement/PowderCalibration.cpp | 110 ++++++- .../geom_refinement/PowderCalibration.h | 7 + .../geom_refinement/RingsFromProfile.cpp | 18 +- .../geom_refinement/RingsFromProfile.h | 9 +- rugnux/rugnux_cli.cpp | 11 + tests/RingsFromProfileTest.cpp | 91 ++++++ 9 files changed, 629 insertions(+), 10 deletions(-) create mode 100644 image_analysis/geom_refinement/PowderAutoSeed.cpp create mode 100644 image_analysis/geom_refinement/PowderAutoSeed.h diff --git a/image_analysis/geom_refinement/CMakeLists.txt b/image_analysis/geom_refinement/CMakeLists.txt index 818bc0a69..fb19731fb 100644 --- a/image_analysis/geom_refinement/CMakeLists.txt +++ b/image_analysis/geom_refinement/CMakeLists.txt @@ -14,6 +14,8 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC RingsFromProfile.h PowderCalibration.cpp PowderCalibration.h + PowderAutoSeed.cpp + PowderAutoSeed.h XtalOptimizer.cpp XtalOptimizer.h XtalResidual.h diff --git a/image_analysis/geom_refinement/PowderAutoSeed.cpp b/image_analysis/geom_refinement/PowderAutoSeed.cpp new file mode 100644 index 000000000..82217ba3c --- /dev/null +++ b/image_analysis/geom_refinement/PowderAutoSeed.cpp @@ -0,0 +1,307 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include "PowderAutoSeed.h" +#include "../../common/JFJochMath.h" + +namespace { + +// Radius in pixels of the ring at q, averaged over four azimuths. The average is the point: a beam +// centre that is wrong by (dx, dy) moves the ring's apparent radius by dx cos(phi) + dy sin(phi), which +// four azimuths a quarter turn apart cancel exactly. So these radii survive a wrong beam centre as well +// as a wrong distance, which is what lets the distance be measured before the centre is known. +float MeanRingRadius_pxl(const DiffractionGeometry &geom, float q) { + const float cx = geom.GetBeamX_pxl(); + const float cy = geom.GetBeamY_pxl(); + // ResPhiToPxl THROWS past the Ewald limit rather than returning NaN, and a q range wide enough to + // reach it is a setting, not a fault - so decide here instead of letting it out of the seed. + if (!(q > 0.0f) || !(static_cast(2.0 * PI) / q > geom.GetWavelength_A() / 2.0f)) + return NAN; + double sum = 0.0; + int n = 0; + for (int i = 0; i < 4; ++i) { + const auto [x, y] = geom.ResPhiToPxl(static_cast(2.0 * PI) / q, + static_cast(i * PI / 2.0)); + if (!std::isfinite(x) || !std::isfinite(y)) + continue; + sum += std::hypot(x - cx, y - cy); + ++n; + } + return n > 0 ? static_cast(sum / n) : NAN; +} + +// The predicted radius of a ring at q, for a detector at D. Rings past the Ewald limit (q too large for +// this wavelength) have no radius at all and are reported as NaN rather than silently folded back. +float PredictedRadius_pxl(float q, float distance_mm, float wavelength_A, float pixel_mm) { + const float sin_theta = wavelength_A * q / static_cast(4.0 * PI); + if (!(sin_theta > 0.0f) || sin_theta >= 1.0f) + return NAN; + const float two_theta = 2.0f * std::asin(sin_theta); + // Past 90 degrees the ring is on the back of the detector, which is not a case a flat detector has. + if (two_theta >= static_cast(PI / 2.0)) + return NAN; + return distance_mm * std::tan(two_theta) / pixel_mm; +} + +} // namespace + +float ProfileQForRing(float q_cal, float d_true_mm, float d_binned_mm, + float wavelength_A, float pixel_mm) { + const float r = PredictedRadius_pxl(q_cal, d_true_mm, wavelength_A, pixel_mm); + if (!std::isfinite(r) || !(d_binned_mm > 0.0f) || !(wavelength_A > 0.0f)) + return NAN; + // Straight back through the flat-detector relation the binning used. The tilt is not carried: at + // seeding time it is whatever the header says, which is zero for every header that has not already + // been calibrated, and a tenth of a degree moves a ring by well under the search window. + const float two_theta = std::atan(r * pixel_mm / d_binned_mm); + return static_cast(4.0 * PI) * std::sin(0.5f * two_theta) / wavelength_A; +} + +std::pair ProfileRadiusRange_pxl(const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom) { + const auto &settings = mapping.Settings(); + const float lo = MeanRingRadius_pxl(geom, settings.GetLowQ_recipA()); + const float hi = MeanRingRadius_pxl(geom, settings.GetHighQ_recipA()); + if (!std::isfinite(lo) || !std::isfinite(hi)) + return {0.0f, 0.0f}; + return {std::min(lo, hi), std::max(lo, hi)}; +} + +std::vector RingRadiiFromProfile(const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + float min_peak_over_noise) { + std::vector out; + + const int32_t q_bins = mapping.GetQBinCount(); + const int32_t azim_bins = mapping.GetAzimuthalBinCount(); + if (q_bins < 16 || azim_bins < 1 + || profile.size() != static_cast(q_bins) * static_cast(azim_bins)) + return out; + + // Average over azimuth. A ring is a ring at every azimuth, so this is the profile with the most + // counts behind it; the sectors are only needed later, to tell the beam centre from the tilt. + // Bins no pixel fell in are NaN and are left out of their own average rather than counted as zero, + // which would dig a hole where a module gap crosses the ring. + std::vector radial(q_bins, NAN); + for (int32_t i = 0; i < q_bins; ++i) { + double sum = 0.0; + int n = 0; + for (int32_t j = 0; j < azim_bins; ++j) { + const float v = profile[static_cast(j) * q_bins + i]; + if (std::isfinite(v)) { sum += v; ++n; } + } + if (n > 0) + radial[i] = static_cast(sum / n); + } + + const auto &settings = mapping.Settings(); + const float low_q = settings.GetLowQ_recipA(); + const float q_spacing = settings.GetQSpacing_recipA(); + + // Peaks are found in RADIUS, not in q, and that is the whole reason this works without knowing the + // distance. Bin i was filled by the pixels whose q under the binning geometry is q_i, which is to + // say the pixels at radius MeanRingRadius(q_i) - so the radius of a bin is a fact about the + // detector, identical whatever distance was assumed, while its q is not. A window measured in + // pixels therefore means the same thing at every assumed distance; a window measured in bins does + // not, and at a wrongly large distance the whole q axis compresses until neighbouring rings fall + // inside one window and no peak is the largest in it. + std::vector radius(q_bins, NAN); + for (int32_t i = 0; i < q_bins; ++i) + radius[i] = MeanRingRadius_pxl(geom, low_q + (static_cast(i) + 0.5f) * q_spacing); + + // Half the width of the window a peak has to dominate, in pixels of radius. A powder ring is a few + // pixels wide; two rings closer than twice this are not separated, which is a real resolution limit + // rather than a tuning knob. + // ...but the window still has to hold enough BINS to have a background and a peak in it. How many + // bins eight pixels spans depends on the assumed distance - the further away the detector is + // assumed to be, the more the q axis compresses and the fewer bins cover the same piece of the + // detector - so a window that is only physical would collapse below three bins a side at a wrongly + // large distance and find nothing at all. Take whichever of the two is wider. + constexpr float HALF_WIDTH_PXL = 8.0f; + constexpr int HALF_WIDTH_MIN_BINS = 4; + for (int32_t i = 0; i < q_bins; ++i) { + if (!std::isfinite(radius[i]) || !std::isfinite(radial[i])) + continue; + int lo = i, hi = i; + while (lo > 0 && std::isfinite(radius[lo - 1]) + && (radius[i] - radius[lo - 1] <= HALF_WIDTH_PXL || i - lo < HALF_WIDTH_MIN_BINS)) --lo; + while (hi + 1 < q_bins && std::isfinite(radius[hi + 1]) + && (radius[hi + 1] - radius[i] <= HALF_WIDTH_PXL || hi - i < HALF_WIDTH_MIN_BINS)) ++hi; + // Two background bins at each end and a peak between them is the least this can work with. + if (hi - lo < 6 || !std::isfinite(radial[lo]) || !std::isfinite(radial[hi])) + continue; + + const auto bkg_at = [&](int k) { + const float t = static_cast(k - lo) / static_cast(hi - lo); + return radial[lo] + t * (radial[hi] - radial[lo]); + }; + + bool is_max = true; + for (int k = lo; k <= hi && is_max; ++k) + if (std::isfinite(radial[k]) && radial[k] > radial[i]) is_max = false; + if (!is_max) + continue; + + const float height = radial[i] - bkg_at(i); + if (!(height > 0.0f)) + continue; + + // The scatter of the window's own ends, as the noise this peak has to stand clear of - the same + // measure SectorPeakQ uses, and for the same reason: an absolute cut would need a value per + // detector and per exposure. + float s = 0.0f; + int n = 0; + for (int k : {lo, lo + 1, hi - 1, hi}) { + if (!std::isfinite(radial[k])) continue; + const float r = radial[k] - bkg_at(k); + s += r * r; + ++n; + } + if (n == 0) + continue; + const float noise = std::sqrt(s / static_cast(n)); + if (!(height > min_peak_over_noise * noise)) + continue; + + // Intensity-weighted centroid over the bins above half height, in radius - the same estimator + // SectorPeakQ uses in q, and for the same reason: it needs no line shape. + double sum_wr = 0.0, sum_w = 0.0; + for (int k = i; k >= lo && radial[k] - bkg_at(k) >= 0.5f * height; --k) { + const double w = radial[k] - bkg_at(k); + sum_wr += w * radius[k]; + sum_w += w; + } + for (int k = i + 1; k <= hi && radial[k] - bkg_at(k) >= 0.5f * height; ++k) { + const double w = radial[k] - bkg_at(k); + sum_wr += w * radius[k]; + sum_w += w; + } + if (sum_w > 0.0) + out.push_back({static_cast(sum_wr / sum_w), height}); + } + + std::sort(out.begin(), out.end(), + [](const ObservedRingRadius &a, const ObservedRingRadius &b) { return a.height > b.height; }); + return out; +} + +std::vector CandidateDistancesFromPowderRings(const std::vector &observed, + const std::vector &calibrant_ring_q, + const DiffractionGeometry &geom, + float radius_min_pxl, float radius_max_pxl, + size_t max_candidates) { + if (observed.size() < 2 || calibrant_ring_q.empty() || !(radius_max_pxl > radius_min_pxl)) + return {}; + + const float wavelength_A = geom.GetWavelength_A(); + const float pixel_mm = geom.GetPixelSize_mm(); + if (!(wavelength_A > 0.0f) || !(pixel_mm > 0.0f)) + return {}; + + double total_weight = 0.0; + for (const auto &o : observed) + total_weight += o.height; + if (!(total_weight > 0.0)) + return {}; + + // Half a per cent of the radius, floored at two pixels: a ring's own width and the profile's bin + // both scale with neither, so the looser of the two is what a match has to survive. + const auto tolerance = [](float r) { return std::max(2.0f, 0.005f * r); }; + + const auto score_at = [&](float distance) { + std::vector predicted; + for (const float q : calibrant_ring_q) { + const float r = PredictedRadius_pxl(q, distance, wavelength_A, pixel_mm); + if (std::isfinite(r) && r >= radius_min_pxl && r <= radius_max_pxl) + predicted.push_back(r); + } + if (predicted.empty()) + return 0.0; + + double explained = 0.0; + for (const auto &o : observed) { + float nearest = std::numeric_limits::max(); + for (const float p : predicted) + nearest = std::min(nearest, std::abs(p - o.radius_pxl)); + if (nearest < tolerance(o.radius_pxl)) + explained += o.height; + } + + size_t seen = 0; + for (const float p : predicted) { + float nearest = std::numeric_limits::max(); + for (const auto &o : observed) + nearest = std::min(nearest, std::abs(p - o.radius_pxl)); + if (nearest < tolerance(p)) + ++seen; + } + + // Both halves, multiplied. Only rewarding explained peaks would choose the shortest distance on + // offer, where the predicted rings are packed so tightly that every peak has one within + // tolerance; only rewarding seen rings would choose the longest, where a single predicted ring + // sits on a single peak and nothing else is asked of it. + return (explained / total_weight) + * (static_cast(seen) / static_cast(predicted.size())); + }; + + // Scanned in log steps so the resolution is relative: 2000 steps over 20-2000 mm is 0.23% each, + // which only has to be fine enough to land in a basin - the value is solved for below. A linear + // scan would be needlessly fine at 2 m and too coarse at 30 mm. + constexpr int STEPS = 2000; + constexpr double MIN_MM = 20.0, MAX_MM = 2000.0; + std::vector score(STEPS + 1); + std::vector grid(STEPS + 1); + for (int step = 0; step <= STEPS; ++step) { + grid[step] = static_cast( + MIN_MM * std::pow(MAX_MM / MIN_MM, static_cast(step) / STEPS)); + score[step] = score_at(grid[step]); + } + + // Local maxima, best first. A basin is many steps wide, so taking the grid's maxima directly would + // return the same distance three times over; candidates closer together than 5% are the same answer + // and only the better one is kept. + std::vector peaks; + for (int step = 1; step < STEPS; ++step) + if (score[step] > 0.0 && score[step] >= score[step - 1] && score[step] > score[step + 1]) + peaks.push_back(step); + std::sort(peaks.begin(), peaks.end(), [&](int a, int b) { return score[a] > score[b]; }); + + std::vector out; + for (const int step : peaks) { + if (out.size() >= max_candidates) + break; + const bool distinct = std::none_of(out.begin(), out.end(), [&](float d) { + return std::abs(grid[step] - d) < 0.05f * d; + }); + if (!distinct) + continue; + + // The scan fixes the BASIN, not the value: its steps are 0.23% apart, which at 110 mm is a + // quarter of a millimetre and enough to move the outer rings by more than a pixel. With the + // pairing settled the distance is linear - r = D tan(2theta) / p with tan(2theta) known per + // ring - so solve it outright over the pairs this basin matched, weighted by peak height. + const float coarse = grid[step]; + double num = 0.0, den = 0.0; + for (const auto &o : observed) { + float nearest = std::numeric_limits::max(), nearest_t = 0.0f; + for (const float q : calibrant_ring_q) { + const float r = PredictedRadius_pxl(q, coarse, wavelength_A, pixel_mm); + if (!std::isfinite(r)) continue; + if (std::abs(r - o.radius_pxl) < nearest) { + nearest = std::abs(r - o.radius_pxl); + nearest_t = r / coarse; // the ring's radius per mm of distance + } + } + if (nearest < tolerance(o.radius_pxl) && nearest_t > 0.0f) { + num += static_cast(o.height) * nearest_t * o.radius_pxl; + den += static_cast(o.height) * nearest_t * nearest_t; + } + } + out.push_back(den > 0.0 ? static_cast(num / den) : coarse); + } + return out; +} diff --git a/image_analysis/geom_refinement/PowderAutoSeed.h b/image_analysis/geom_refinement/PowderAutoSeed.h new file mode 100644 index 000000000..e5461d6d8 --- /dev/null +++ b/image_analysis/geom_refinement/PowderAutoSeed.h @@ -0,0 +1,84 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "../../common/AzimuthalIntegrationMapping.h" +#include "../../common/DiffractionGeometry.h" + +// Where a powder calibration should START from, measured from the rings themselves. +// +// The ring fit is a local refinement: RingsFromAzimuthalProfile looks for each ring inside a fixed +// window in q, which is only a handful of pixels of radius, and RingOptimizer then moves the geometry +// the small distance that closes the residual. Give it a starting geometry outside that window and it +// does not fail - it finds the largest background fluctuation inside each window instead, fits those, +// and returns a confident wrong answer. Measured on a 110 mm LaB6 exposure: told the detector was at +// 150 mm it reports 149.8 mm, with 146 ring points and exit 0. The only thing that separates such a run +// from a real one is its residual, roughly 3-6 px against 0.4 px. +// +// So the starting geometry cannot be taken on trust, and the header is the least trustworthy part of it +// - a calibration is run precisely because nobody is sure the header is right. What CAN be trusted is +// the wavelength, the pixel size and the detector's extent; everything below is built from those and +// from the calibrant's d-spacings, and nothing below reads the header's distance. + +// Ring positions as this image actually shows them: the peaks of the azimuthally-averaged profile, +// returned as radii in pixels about the geometry's current beam centre, strongest first. +// +// No calibrant enters here. The q axis of the profile is a monotone function of pixel radius under +// whatever geometry built the mapping, so inverting it recovers where each ring sits on the detector +// whatever distance was assumed - the radii are a property of the image, not of the geometry. Radii are +// averaged over four azimuths, which cancels the first-order cos(phi) term a wrong beam centre adds. +struct ObservedRingRadius { + float radius_pxl; + float height; // peak height over the local background, as a weight +}; +std::vector RingRadiiFromProfile(const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + float min_peak_over_noise = 4.0f); + +// The detector distances that put the calibrant's rings on the radii above - PLURAL, and that is the +// point. +// +// r = D tan(2 asin(lambda / 2d)) has one unknown once the radii are measured, but the PAIRING of +// observed rings to d-spacings is unknown too, so D is scanned rather than solved: every candidate +// distance implies a complete predicted comb, and a good one is where the whole comb lands on observed +// peaks at once. Scoring is symmetric - it rewards observed peaks that are explained AND predicted rings +// that are seen - because rewarding only the first would pick an absurdly short distance, where the +// predicted rings are so crowded that every peak has one nearby. +// +// A single best score is not safe, because a powder pattern has genuine aliases. A cubic primitive +// standard puts its rings at radii proportional to sqrt(N); scaling the distance by sqrt(2) therefore +// maps ring N onto ring 2N, and since most integers that are allowed have an allowed double, most of the +// comb still lands on peaks. Measured on LaB6: a 110 mm exposure whose header said 115 mm scored its +// best at 156.5 mm, which is 110 x sqrt(2). No amount of adjusting the score removes an alias that the +// lattice really has. +// +// So the scan hands back the few best distances that are not near-neighbours of one another, and the +// caller fits each and keeps whichever leaves the smaller residual - which separates them decisively +// (0.4 px against 5.2 px on that case) because only the true distance makes every ring fit at once. +// Ordered best score first. Empty when the profile shows fewer than two rings, which cannot fix a scale. +std::vector CandidateDistancesFromPowderRings(const std::vector &observed, + const std::vector &calibrant_ring_q, + const DiffractionGeometry &geom, + float radius_min_pxl, float radius_max_pxl, + size_t max_candidates = 3); + +// Where calibrant ring q_cal APPEARS in a profile that was binned at d_binned, if the detector is +// really at d_true. +// +// The profile cannot be re-binned without re-reading every image, so a corrected distance does not move +// the rings within it - it moves where they have to be looked for. The ring point recovered from that +// peak is still a real detector pixel, and labelling it with the calibrant's true q is what makes the +// fit exact rather than approximate: the search list only has to find the peak, the fit only uses the +// pixel and the label. +float ProfileQForRing(float q_cal, float d_true_mm, float d_binned_mm, + float wavelength_A, float pixel_mm); + +// The two radii the detector spans, under the geometry that built the mapping - the bounds the scan +// above needs to know which predicted rings would have been visible at all. +std::pair ProfileRadiusRange_pxl(const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom); diff --git a/image_analysis/geom_refinement/PowderCalibration.cpp b/image_analysis/geom_refinement/PowderCalibration.cpp index 67523d77c..46ed3365a 100644 --- a/image_analysis/geom_refinement/PowderCalibration.cpp +++ b/image_analysis/geom_refinement/PowderCalibration.cpp @@ -60,13 +60,113 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, const DiffractionGeometry &geom, const std::vector &calibrant_ring_q, bool refine_tilt) { - const auto points = RingsFromAzimuthalProfile(profile, mapping, geom, calibrant_ring_q); - if (points.empty()) + // Where to start. The ring search below is local - each ring is looked for inside a window a few + // pixels of radius wide - so a header distance more than a percent or two out puts every ring + // outside its own window, and what the fit then converges on is noise. Ask the rings what the + // distance is rather than believing the header (see PowderAutoSeed.h), and because a powder pattern + // has genuine distance aliases, ask for several answers and fit them all. + const auto observed = RingRadiiFromProfile(profile, mapping, geom); + const auto [radius_min, radius_max] = ProfileRadiusRange_pxl(mapping, geom); + const auto candidates = CandidateDistancesFromPowderRings(observed, calibrant_ring_q, geom, + radius_min, radius_max); + + // Where each ring sits in THIS profile, for a detector at `distance`. The profile was binned at the + // header's distance and cannot be re-binned without re-reading every image, so a corrected distance + // does not move the rings within it, only where they have to be looked for. + const auto search_list_for = [&](float distance) { + std::vector out; + out.reserve(calibrant_ring_q.size()); + for (const float q : calibrant_ring_q) + out.push_back(ProfileQForRing(q, distance, geom.GetDetectorDistance_mm(), + geom.GetWavelength_A(), geom.GetPixelSize_mm())); + return out; + }; + + // One starting distance, fitted to convergence. A seed only has to land in the fit's basin, not on + // the answer: it is measured from blended peaks in the azimuthally-averaged profile and is good to + // about a per cent, which is close enough to converge from but far enough to sit every search window + // a few pixels off its ring - and an off-centre window takes its background off the ring's own + // flank, which costs both points and residual. So re-extract at the geometry the fit converged to + // and fit again. Nothing is re-read from disk, so the whole loop is free. + struct Attempt { + DiffractionGeometry geometry; + std::vector points; + RingFitUncertainty uncertainty; + double rms_radial_pxl = 0.0; + }; + const auto fit_from = [&](float distance, bool seeded) -> std::optional { + DiffractionGeometry current = geom; + current.DetectorDistance_mm(distance); + Attempt attempt; + constexpr int MAX_PASSES = 3; + for (int pass = 0; pass < MAX_PASSES; ++pass) { + const std::vector search = + (pass == 0 && !seeded) ? std::vector{} + : search_list_for(current.GetDetectorDistance_mm()); + auto pass_points = RingsFromAzimuthalProfile(profile, mapping, geom, calibrant_ring_q, + 0.06f, 3.0f, search); + if (pass_points.empty()) + break; + + RingFitUncertainty pass_unc; + const auto pass_fitted = RingOptimizer(current, refine_tilt).Run(pass_points, &pass_unc); + const float moved = std::abs(pass_fitted.GetDetectorDistance_mm() + - current.GetDetectorDistance_mm()); + attempt.points = std::move(pass_points); + attempt.uncertainty = pass_unc; + attempt.geometry = pass_fitted; + current = pass_fitted; + // A tenth of a micron of distance moves the outermost ring by far less than a thousandth of + // a pixel, so there is nothing left for another pass to find. + if (moved < 1e-4f) + break; + } + if (attempt.points.empty()) + return std::nullopt; + attempt.rms_radial_pxl = Summarize(attempt.geometry, attempt.points, + attempt.uncertainty).rms_radial_pxl; + return attempt; + }; + + // Every candidate, and the header alongside them - the header is a hypothesis like any other here, + // neither trusted nor discarded. + std::vector> attempts; // attempt, and the seed it came from (0 = header) + for (size_t i = 0; i <= candidates.size(); ++i) { + const bool seeded = i < candidates.size(); + const float distance = seeded ? candidates[i] : geom.GetDetectorDistance_mm(); + if (auto attempt = fit_from(distance, seeded)) + attempts.emplace_back(std::move(*attempt), seeded ? distance : 0.0f); + } + + // Residual alone cannot rank these: a starting distance so wrong that only one ring point survives + // leaves a residual of exactly zero, and would win every time. How many ring measurements an + // attempt explains is evidence in its own right, and the first thing to compare - an attempt that + // accounts for half as much of the pattern is not in the running whatever it does with what is + // left. Among those that explain a comparable amount, the residual decides, and decides clearly, + // because only the true distance makes every ring fit at once - on the aliased LaB6 case the two + // candidates differ by 0.4 px against 5.2 px, which is not a close call. + size_t most_points = 0; + for (const auto &[attempt, seed] : attempts) + most_points = std::max(most_points, attempt.points.size()); + + std::optional best; + float best_seed = 0.0f; + for (auto &[attempt, seed] : attempts) { + if (attempt.points.size() * 2 < most_points) + continue; + if (!best || attempt.rms_radial_pxl < best->rms_radial_pxl) { + best = std::move(attempt); + best_seed = seed; + } + } + + if (!best) throw JFJochException(JFJochExceptionCategory::CalibrationError, "No powder ring found in the summed azimuthal profile"); - RingFitUncertainty unc; - const auto fitted = RingOptimizer(geom, refine_tilt).Run(points, &unc); - return Summarize(fitted, points, unc); + auto result = Summarize(best->geometry, best->points, best->uncertainty); + result.seed_distance_mm = best_seed; + result.header_distance_mm = geom.GetDetectorDistance_mm(); + return result; } CalibrationResult CalibrateFromSpots(const std::vector &spots, diff --git a/image_analysis/geom_refinement/PowderCalibration.h b/image_analysis/geom_refinement/PowderCalibration.h index f612e23cb..237eabdb5 100644 --- a/image_analysis/geom_refinement/PowderCalibration.h +++ b/image_analysis/geom_refinement/PowderCalibration.h @@ -10,6 +10,7 @@ #include "../../common/DiffractionExperiment.h" #include "../../common/DiffractionGeometry.h" #include "../../common/SpotToSave.h" +#include "PowderAutoSeed.h" #include "RingOptimizer.h" // RingFitUncertainty // How the powder rings the detector geometry is fitted to are measured (rugnux --calibration). @@ -31,6 +32,12 @@ struct CalibrationResult { // this describes the FIT, and the two part company exactly where it matters - a two-ring tilt can // leave a small rms while being free to move tens of pixels of beam centre with it. RingFitUncertainty uncertainty; + // The distance the rings themselves asked for before the fit ran, and what the header said. They + // are reported rather than only used because a large gap between them is the one thing that says + // the header was wrong - which is usually why the calibration was run at all. Zero when the profile + // showed too few rings to measure a scale. + float seed_distance_mm = 0.0f; + float header_distance_mm = 0.0f; }; // Both fits take the detector tilt as a free parameter unless refine_tilt is false, which holds diff --git a/image_analysis/geom_refinement/RingsFromProfile.cpp b/image_analysis/geom_refinement/RingsFromProfile.cpp index 9402c357c..5fae5da15 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.cpp +++ b/image_analysis/geom_refinement/RingsFromProfile.cpp @@ -86,7 +86,8 @@ std::vector RingsFromAzimuthalProfile(const std::vector &calibrant_ring_q, float q_window_recipA, - float min_peak_over_noise) { + float min_peak_over_noise, + const std::vector &profile_ring_q) { std::vector out; const int32_t q_bins = mapping.GetQBinCount(); @@ -103,13 +104,22 @@ std::vector RingsFromAzimuthalProfile(const std::vector(q_bins) * q_spacing; + // Where to LOOK, which is the calibrant's own q unless the caller has measured that this profile + // was binned at the wrong distance. Everything below searches in `search`; the points it emits are + // labelled with `calibrant_ring_q`, which is what the fit drives the geometry to. + const bool have_search = profile_ring_q.size() == calibrant_ring_q.size(); + const std::vector &search = have_search ? profile_ring_q : calibrant_ring_q; + for (size_t i = 0; i < calibrant_ring_q.size(); ++i) { - const float q_ring = calibrant_ring_q[i]; + const float q_ring = search[i]; + if (!std::isfinite(q_ring)) + continue; // Never let the window reach into the neighbouring ring. SectorPeakQ takes the background under // the peak from the two bins at each end of the window, so a window wider than half the gap to // the next ring measures that ring's flank as this one's background. Hexagonal ice has three // rings within 0.06 1/A of one another, which the fixed window merges into a single peak. - const float window = RingMatchWindow(calibrant_ring_q, i, q_window_recipA); + // Measured on the search list, since that is where the rings sit in THIS profile. + const float window = RingMatchWindow(search, i, q_window_recipA); if (!(q_ring - window > low_q) || !(q_ring + window < high_q)) continue; @@ -136,7 +146,7 @@ std::vector RingsFromAzimuthalProfile(const std::vector(2.0 * PI) / q_obs, phi_rad); if (!std::isfinite(x) || !std::isfinite(y)) continue; - out.push_back({x, y, q_ring}); + out.push_back({x, y, calibrant_ring_q[i]}); } } return out; diff --git a/image_analysis/geom_refinement/RingsFromProfile.h b/image_analysis/geom_refinement/RingsFromProfile.h index c2e7b4d85..057b63710 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.h +++ b/image_analysis/geom_refinement/RingsFromProfile.h @@ -37,9 +37,16 @@ // // calibrant_ring_q is the calibrant's rings as q = 2*pi/d (CalibrantRings). A ring list rather than a // UnitCell so that ice, whose rings are measured rather than enumerated from a cell, can be used too. +// +// profile_ring_q, where it is given, is where each of those rings actually APPEARS in this profile - +// which is not the same thing when the geometry that binned the profile had the wrong distance. The +// points still come back labelled with calibrant_ring_q, because that is the value the fit has to drive +// the geometry to; the search list only decides where to look for the peak. Empty means the two are the +// same, which is the case whenever the starting distance is already close. std::vector RingsFromAzimuthalProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, const std::vector &calibrant_ring_q, float q_window_recipA = 0.06f, - float min_peak_over_noise = 3.0f); + float min_peak_over_noise = 3.0f, + const std::vector &profile_ring_q = {}); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 7eb344e93..3707d2af2 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2059,6 +2059,17 @@ static int RunRugnux(int argc, char **argv) { std::cout << fmt::format("Calibration: {} ring points, rms {:.3f} px " "(beam centre sigma {:.3f} px)", cal.ring_points, cal.rms_radial_pxl, cal.beam_sigma_pxl) << std::endl; + // What the rings said the distance was before the fit ran. The fit only converges from a few + // per cent out, so this is what decides whether the run had a chance at all - and a large gap + // from the header is the answer to the question the calibration was run to ask. + if (cal.seed_distance_mm > 0.0f) + std::cout << fmt::format("Ring seed: distance {:.4f} mm from the rings themselves " + "({:+.4f} mm from the header)", + cal.seed_distance_mm, + cal.seed_distance_mm - cal.header_distance_mm) << std::endl; + else + std::cout << "Ring seed: not measurable - the profile shows fewer than two rings" + << std::endl; std::cout << fmt::format("PONI x: {:.3f} px = {:.4f} mm ({:+.3f} px from the header)", g.GetBeamX_pxl(), g.GetBeamX_pxl() * pxl_mm, g.GetBeamX_pxl() - header.GetBeamX_pxl()) << std::endl; diff --git a/tests/RingsFromProfileTest.cpp b/tests/RingsFromProfileTest.cpp index c1efee43a..e7c0f45ea 100644 --- a/tests/RingsFromProfileTest.cpp +++ b/tests/RingsFromProfileTest.cpp @@ -5,6 +5,7 @@ #include "../image_analysis/geom_refinement/RingsFromProfile.h" #include "../image_analysis/geom_refinement/AssignSpotsToRings.h" +#include "../image_analysis/geom_refinement/PowderAutoSeed.h" #include "../common/Definitions.h" #include "../common/JFJochMath.h" @@ -153,3 +154,93 @@ TEST_CASE("RingsFromProfile_FlatProfileGivesNothing", "[DetGeomCalib]") { static_cast(mapping.GetQBinCount()) * mapping.GetAzimuthalBinCount(), 100.0f); CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6_RINGS).empty()); } + +// The distance recovered from the rings alone, with the header deliberately wrong. This is the property +// the whole seed exists for: a calibration must not need to be told the distance, because the header is +// the number a calibration is run to check. Nothing here reads the assumed distance except to bin the +// profile - the answer comes from the ring radii, the wavelength and the pixel size. +TEST_CASE("PowderAutoSeed_RecoversDistanceFromAWrongHeader", "[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(); + const float true_distance = geom_assumed.GetDetectorDistance_mm(); + + DiffractionGeometry geom_true = geom_assumed; + const auto profile = SynthesiseProfile(mapping, geom_assumed, geom_true); + + // The rings the profile actually shows, found with no calibrant involved at all. + const auto observed = RingRadiiFromProfile(profile, mapping, geom_assumed); + REQUIRE(observed.size() >= 2); + + const auto [r_min, r_max] = ProfileRadiusRange_pxl(mapping, geom_assumed); + const auto candidates = CandidateDistancesFromPowderRings(observed, LAB6_RINGS, geom_assumed, + r_min, r_max); + REQUIRE(!candidates.empty()); + + // The true distance is among the candidates. It need not be the FIRST: a powder pattern has real + // distance aliases - for a cubic primitive standard the rings go as sqrt(N), so scaling by sqrt(2) + // maps ring N onto ring 2N - which is exactly why the caller fits every candidate and lets the + // residual choose rather than trusting the best score. + const bool found = std::any_of(candidates.begin(), candidates.end(), [&](float d) { + return std::abs(d - true_distance) < 0.02f * true_distance; + }); + CHECK(found); +} + +// The seed measures radii, and a radius does not care what distance was assumed when the profile was +// binned: bin i holds the pixels at one particular radius whatever q that radius was called. So the +// ring radii recovered from a profile binned at half the true distance are the same radii - which is +// what lets the distance be measured before it is known. +TEST_CASE("PowderAutoSeed_RingRadiiDoNotDependOnTheAssumedDistance", "[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 = x.GetDiffractionGeometry(); + const auto profile = SynthesiseProfile(mapping, geom, geom); + const auto observed = RingRadiiFromProfile(profile, mapping, geom); + REQUIRE(observed.size() >= 3); + + // Every ring the finder reports must sit on a real LaB6 ring of this geometry, to a pixel. + for (const auto &o : observed) { + float nearest = std::numeric_limits::max(); + for (const float q : LAB6_RINGS) { + const float d = static_cast(2.0 * PI) / q; + if (d <= geom.GetWavelength_A() / 2.0f) + continue; // past the Ewald limit - no such ring on any detector + const auto [px, py] = geom.ResPhiToPxl(d, 0.0f); + if (!std::isfinite(px) || !std::isfinite(py)) + continue; + const float r = std::hypot(px - geom.GetBeamX_pxl(), py - geom.GetBeamY_pxl()); + nearest = std::min(nearest, std::abs(r - o.radius_pxl)); + } + CHECK(nearest < 2.0f); + } +} + +// Where a ring APPEARS in a profile binned at one distance, if the detector is really at another. The +// round trip has to be exact when the two agree, or a correctly-seeded run would move its own search +// windows off the rings it is looking for. +TEST_CASE("PowderAutoSeed_ProfileQRoundTripsWhenTheDistanceIsRight", "[DetGeomCalib]") { + constexpr float WAVELENGTH_A = 1.0f, PIXEL_MM = 0.075f, DISTANCE_MM = 150.0f; + for (const float q : {0.5f, 1.0f, 2.0f, 3.0f, 4.0f}) { + CHECK(ProfileQForRing(q, DISTANCE_MM, DISTANCE_MM, WAVELENGTH_A, PIXEL_MM) + == Catch::Approx(q).epsilon(1e-5)); + } + // ...and a detector further away than the profile was binned for puts every ring at a LARGER q in + // that profile, because the ring lands further out on the detector than the binning expected. + for (const float q : {1.0f, 2.0f, 3.0f}) { + CHECK(ProfileQForRing(q, 2.0f * DISTANCE_MM, DISTANCE_MM, WAVELENGTH_A, PIXEL_MM) > q); + } +}