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); + } +}