diff --git a/image_analysis/geom_refinement/PowderAutoSeed.cpp b/image_analysis/geom_refinement/PowderAutoSeed.cpp index 0a1127927..50b4c84bf 100644 --- a/image_analysis/geom_refinement/PowderAutoSeed.cpp +++ b/image_analysis/geom_refinement/PowderAutoSeed.cpp @@ -5,6 +5,7 @@ #include #include "PowderAutoSeed.h" +#include "RingsFromProfile.h" // SectorPeakQ #include "../../common/JFJochMath.h" namespace { @@ -48,16 +49,148 @@ float PredictedRadius_pxl(float q, float distance_mm, float wavelength_A, float } // 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::vector ProfileRingTrack(float q_cal, const DiffractionGeometry &truth, + const DiffractionGeometry &binned, int32_t azim_bins) { + std::vector track(std::max(azim_bins, 1), NAN); + if (azim_bins < 1 || !(q_cal > 0.0f)) + return track; + const float d = static_cast(2.0 * PI) / q_cal; + if (!(d > truth.GetWavelength_A() / 2.0f)) + return track; // past the Ewald limit - this ring is on no detector + + // Walk the ring in `truth` finely enough that every sector is hit several times, average what lands + // in each. Averaging rather than taking one sample per sector because the map from an azimuth in + // `truth` to a sector of `binned` is not uniform once the two centres differ. + std::vector sum(track.size(), 0.0); + std::vector count(track.size(), 0); + const int samples = 8 * azim_bins; + for (int i = 0; i < samples; ++i) { + const float phi = static_cast(2.0 * PI * i / samples); + const auto [x, y] = truth.ResPhiToPxl(d, phi); + if (!std::isfinite(x) || !std::isfinite(y)) + continue; + const float q = binned.PxlToQ(x, y); + if (!std::isfinite(q) || !(q > 0.0f)) + continue; + float phi_binned = binned.Phi_rad(x, y); + if (!std::isfinite(phi_binned)) + continue; + auto bin = static_cast(phi_binned / static_cast(2.0 * PI) + * static_cast(azim_bins)); + bin = std::clamp(bin, 0, azim_bins - 1); + sum[bin] += q; + ++count[bin]; + } + for (size_t i = 0; i < track.size(); ++i) + if (count[i] > 0) + track[i] = static_cast(sum[i] / count[i]); + return track; +} + +std::optional> BeamCentreOffsetFromProfile( + const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + const std::vector &observed) { + const int32_t q_bins = mapping.GetQBinCount(); + const int32_t azim_bins = mapping.GetAzimuthalBinCount(); + // Two rings at least: each is searched in a window reaching half way to its nearest neighbour, so + // a single ring has no neighbour to bound it and nothing says which ring a peak inside a boundless + // window belongs to. + if (azim_bins < 4 || q_bins < 8 || observed.size() < 2 + || profile.size() != static_cast(q_bins) * static_cast(azim_bins)) + return std::nullopt; + + const auto &settings = mapping.Settings(); + const float low_q = settings.GetLowQ_recipA(); + const float q_spacing = settings.GetQSpacing_recipA(); + + // Radius of every q bin, so a peak found in q can be stated as a radius - which is the quantity the + // cos(phi) law below is written in, and the only one that does not depend on the assumed distance. + 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); + + // Normal equations for r - R_j = dx cos(phi) + dy sin(phi), pooled over rings. R_j, each ring's own + // mean radius, is eliminated by centring each ring's measurements on their own mean - which is why + // no d-spacing and no distance is needed to get the centre out. + double sxx = 0.0, sxy = 0.0, syy = 0.0, sxr = 0.0, syr = 0.0; + int used_rings = 0; + + for (const auto &ring : observed) { + // Search each ring about the radius the profile actually put it at, in a window reaching half + // way to its neighbours - so this is looking for a ring it has already found rather than for one + // a standard predicts, and a badly wrong beam centre cannot push it out of its own window. + float gap = std::numeric_limits::max(); + for (const auto &other : observed) + if (&other != &ring) + gap = std::min(gap, std::abs(other.radius_pxl - ring.radius_pxl)); + const float window_pxl = 0.5f * gap; + if (!(window_pxl > 2.0f)) + continue; + + std::vector measured(azim_bins, NAN); + for (int32_t phi_bin = 0; phi_bin < azim_bins; ++phi_bin) { + int lo = -1, hi = -1; + for (int32_t i = 0; i < q_bins; ++i) { + if (!std::isfinite(radius[i])) + continue; + if (std::abs(radius[i] - ring.radius_pxl) <= window_pxl) { + if (lo < 0) lo = i; + hi = i; + } + } + if (lo < 0 || hi - lo < 6) + continue; + const float q_obs = SectorPeakQ(profile, q_bins, phi_bin, lo, hi, + low_q, q_spacing, 3.0f); + if (!std::isfinite(q_obs)) + continue; + measured[phi_bin] = MeanRingRadius_pxl(geom, q_obs); + } + + // Half the sectors, and spread round the turn: dx and dy are read off a cos and a sin, so a ring + // seen only on one side of the pattern constrains one combination of them and leaves the other + // free. Requiring the mean of cos and of sin over the sectors used to be small is what says the + // coverage is even enough for the pair to separate. + double mean_r = 0.0, mean_c = 0.0, mean_s = 0.0; + int n = 0; + const auto phi_of = [&](int32_t k) { + return (static_cast(k) + 0.5) * 2.0 * PI / static_cast(azim_bins); + }; + for (int32_t k = 0; k < azim_bins; ++k) { + if (!std::isfinite(measured[k])) continue; + mean_r += measured[k]; + mean_c += std::cos(phi_of(k)); + mean_s += std::sin(phi_of(k)); + ++n; + } + if (n * 2 < azim_bins) + continue; + mean_r /= n; mean_c /= n; mean_s /= n; + if (std::hypot(mean_c, mean_s) > 0.25) + continue; + + for (int32_t k = 0; k < azim_bins; ++k) { + if (!std::isfinite(measured[k])) continue; + const double c = std::cos(phi_of(k)), s = std::sin(phi_of(k)); + const double r = measured[k] - mean_r; + sxx += c * c; sxy += c * s; syy += s * s; + sxr += c * r; syr += s * r; + } + ++used_rings; + } + + if (used_rings == 0) + return std::nullopt; + const double det = sxx * syy - sxy * sxy; + if (!(std::abs(det) > 1e-9)) + return std::nullopt; + const double dx = (sxr * syy - syr * sxy) / det; + const double dy = (syr * sxx - sxr * sxy) / det; + if (!std::isfinite(dx) || !std::isfinite(dy)) + return std::nullopt; + return std::make_pair(static_cast(dx), static_cast(dy)); } std::pair ProfileRadiusRange_pxl(const AzimuthalIntegrationMapping &mapping, diff --git a/image_analysis/geom_refinement/PowderAutoSeed.h b/image_analysis/geom_refinement/PowderAutoSeed.h index d8178bce9..1d72457ae 100644 --- a/image_analysis/geom_refinement/PowderAutoSeed.h +++ b/image_analysis/geom_refinement/PowderAutoSeed.h @@ -79,16 +79,50 @@ std::vector CandidateDistancesFromPowderRings( 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. +// Where calibrant ring q_cal APPEARS in a profile binned with `binned`, sector by sector, if the true +// geometry is `truth`. One entry per azimuthal sector, NaN where the ring misses that sector. // -// 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 profile cannot be re-binned without re-reading every image, so a corrected geometry does not move +// the rings within it - it moves where they have to be looked for. Per SECTOR, not one q for the whole +// ring, and that is what a beam-centre correction needs: a centre wrong by (dx, dy) makes a ring's +// apparent radius oscillate as dx cos(phi) + dy sin(phi), so the ring is at a different q in every +// sector and a single search window centred on one q finds it only where the oscillation happens to be +// small. That is what limited the beam centre the fit could recover from to about ten pixels. +// +// Exact, and exact in all five parameters at once: the ring is walked in `truth`, each point turned +// into a detector pixel, and that pixel asked what q and what azimuth `binned` would have given it. No +// flat-detector approximation, so a tilt is carried too. +// +// The points the caller then recovers from those peaks are real detector pixels and are labelled with +// the calibrant's true q - the track only has to find the peak, the fit only uses the pixel and label. +std::vector ProfileRingTrack(float q_cal, const DiffractionGeometry &truth, + const DiffractionGeometry &binned, int32_t azim_bins); + +// The beam-centre offset the rings themselves ask for, in pixels, to be ADDED to the geometry the +// profile was binned with. No calibrant and no distance enter: a powder ring is a conic centred on the +// beam, so a centre wrong by (dx, dy) makes the apparent radius of EVERY ring oscillate once per turn +// with the same amplitude - r(phi) = R + dx cos(phi) + dy sin(phi) - and that is solved for directly, +// pooled over every ring the profile shows. Each ring is searched about its OWN measured radius rather +// than about where a standard says it should be, which is what makes this work when the header centre +// is far enough out that the calibrated extraction would find nothing. +// +// The offset it can recover is bounded, and by the measurement rather than by a choice. Each ring is +// searched in a window reaching half way to its neighbour in the AZIMUTHALLY AVERAGED profile, and once +// the offset grows past a few pixels that profile stops showing rings: a ring whose radius traces +// R + dx cos(phi) + dy sin(phi) piles up density where r(phi) turns round, so it averages into the two +// HORNS of that sinusoid, at R-|d| and R+|d|. The radius finder then reports two rings where there is +// one, and the gap it measures between them is 2|d| - which is to say the window shrinks to exactly the +// offset it was meant to span. Measured on a 110 mm LaB6 exposure the whole calibration recovers a +// header centre about 20 px out and fails by 40; the limit is roughly half the spacing of the rings. +// Beyond it there is nothing left in an azimuthally binned profile to work from, and --calibration +// spots, which finds the centre from the spot positions themselves, is the method that still can. +// +// Returns nothing when no ring is sampled well enough round the turn to separate the two components. +std::optional> BeamCentreOffsetFromProfile( + const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + const std::vector &observed); // 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. diff --git a/image_analysis/geom_refinement/PowderCalibration.cpp b/image_analysis/geom_refinement/PowderCalibration.cpp index 4ecdd027b..adfbb5192 100644 --- a/image_analysis/geom_refinement/PowderCalibration.cpp +++ b/image_analysis/geom_refinement/PowderCalibration.cpp @@ -84,17 +84,21 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, 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; - }; + // ...and ask them where the beam is, for the same reason. The extraction looks for each ring within + // a window a few pixels of radius wide, so a header centre more than about ten pixels out puts the + // ring outside that window over much of the turn - and the fit then reads a cos(phi) signal off + // whatever sectors are left, which is how a 20 px error used to end 31 px wrong. The rings answer + // this without a calibrant and without a distance: a powder ring is a conic centred on the beam, so + // a wrong centre makes EVERY ring's radius oscillate once per turn by the same amount. + DiffractionGeometry seed_geometry = geom; + auto centre_offset = BeamCentreOffsetFromProfile(profile, mapping, geom, observed); + // Under a pixel it is not a different hypothesis, it is the same one - and fitting it as well would + // double the work for two answers that cannot be told apart. + if (centre_offset && std::hypot(centre_offset->first, centre_offset->second) < 1.0f) + centre_offset.reset(); + if (centre_offset) + seed_geometry.BeamX_pxl(geom.GetBeamX_pxl() + centre_offset->first) + .BeamY_pxl(geom.GetBeamY_pxl() + centre_offset->second); // 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 @@ -108,15 +112,16 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, RingFitUncertainty uncertainty; double rms_radial_pxl = 0.0; }; - const auto fit_from = [&](float distance, bool seeded, bool tilt) -> std::optional { - DiffractionGeometry current = geom; - current.DetectorDistance_mm(distance); + const auto fit_from = [&](const DiffractionGeometry &start, bool seeded, + bool tilt) -> std::vector { + DiffractionGeometry current = start; 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()); + // The profile's own geometry is the search on the first pass of an unseeded attempt, and + // the geometry converged to on every pass after - which is where the rings have been + // measured to be, rather than where the header guessed. + const DiffractionGeometry *search = (pass == 0 && !seeded) ? nullptr : ¤t; auto pass_points = RingsFromAzimuthalProfile(profile, mapping, geom, calibrant_ring_q, 0.06f, 3.0f, search); if (pass_points.empty()) @@ -136,27 +141,78 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, break; } if (attempt.points.empty()) - return std::nullopt; - attempt.rms_radial_pxl = Summarize(attempt.geometry, attempt.points, - attempt.uncertainty).rms_radial_pxl; - return attempt; + return {}; + + // Two ways to take the final measurement, and they genuinely disagree about which is better. + // + // Following the rings sector by sector is what lets a badly placed beam centre be recovered at + // all - a centre wrong by (dx, dy) puts a ring at a different q in every sector, and one window + // centred on one q finds it only where that oscillation happens to be small. But it costs + // precision once they have been found, for a reason worth stating: a window that moves with phi + // makes every systematic of the peak finder - where the background line is taken, how the + // centroid sits in the window - vary with phi too, and phi is exactly the axis the beam centre + // is read off. Measured, rms 0.415 -> 0.525 px on a good 110 mm fit, and 0.831 when the window + // followed the fitted tilt as well. + // + // The alternative is a window that is the same in every sector, which the binned geometry with + // only its DISTANCE replaced gives by construction - a distance error moves every sector's + // window equally, where a centre or a tilt error does not. That is more precise where it works + // and finds nothing where the centre is far out. So do both and let the same rule that ranks + // everything else here decide. + std::vector out; + out.push_back(attempt); + + DiffractionGeometry measure_geom = geom; + measure_geom.DetectorDistance_mm(current.GetDetectorDistance_mm()); + auto measured = RingsFromAzimuthalProfile(profile, mapping, geom, calibrant_ring_q, + 0.06f, 3.0f, &measure_geom); + if (!measured.empty()) { + Attempt fixed_window; + RingFitUncertainty measured_unc; + fixed_window.geometry = RingOptimizer(current, tilt).Run(measured, &measured_unc); + fixed_window.points = std::move(measured); + fixed_window.uncertainty = measured_unc; + out.push_back(std::move(fixed_window)); + } + + for (auto &a : out) + a.rms_radial_pxl = Summarize(a.geometry, a.points, a.uncertainty).rms_radial_pxl; + return out; }; - // Every candidate, and the header alongside them - the header is a hypothesis like any other here, - // neither trusted nor discarded. - // Each attempt, with the seed it started from (0 = the header) and how much of the pattern that - // seed's comb explained. + // Every combination of what the rings said and what the header said, fitted, with the residual left + // to choose. The header is a hypothesis like any other here, neither trusted nor discarded - and so + // is the seeded beam centre, which is NOT simply better than the header's. + // + // The centre seed reads the once-per-turn wobble of the ring radii, and a tilt puts a term of that + // same shape there too - one that grows as the ring's radius squared, where a centre error does not. + // Pooling the rings into one offset therefore absorbs part of the tilt into the centre, which is + // worth several pixels on a genuinely tilted detector and made a good 110 mm fit worse when it was + // simply believed. What it buys is capture range, and only where the header centre is far enough out + // that the extraction would otherwise find the rings over a fraction of the turn. Offering it as an + // alternative start costs one more fit each and needs no rule about when it applies. struct Provenance { float seed_mm; double match; }; - std::vector> attempts; - for (size_t i = 0; i <= candidates.size(); ++i) { - const bool seeded = i < candidates.size(); - const float distance = seeded ? candidates[i].distance_mm : geom.GetDetectorDistance_mm(); - if (auto attempt = fit_from(distance, seeded, refine_tilt)) - attempts.emplace_back(std::move(*attempt), - Provenance{seeded ? distance : 0.0f, - seeded ? candidates[i].score : 0.0}); + struct Start { DiffractionGeometry geometry; bool tracked; Provenance provenance; }; + + std::vector starts; + for (int centre = 0; centre < (centre_offset ? 2 : 1); ++centre) { + const DiffractionGeometry &base = centre == 0 ? geom : seed_geometry; + for (size_t i = 0; i <= candidates.size(); ++i) { + const bool seeded_distance = i < candidates.size(); + DiffractionGeometry start = base; + if (seeded_distance) + start.DetectorDistance_mm(candidates[i].distance_mm); + starts.push_back({start, seeded_distance || centre == 1, + {seeded_distance ? candidates[i].distance_mm : 0.0f, + seeded_distance ? candidates[i].score : 0.0}}); + } } + std::vector> attempts; + for (const auto &start : starts) + for (auto &attempt : fit_from(start.geometry, start.tracked, refine_tilt)) + attempts.emplace_back(std::move(attempt), start.provenance); + // 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 @@ -202,7 +258,19 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, const float significance = tilt_was_free ? TiltSignificance(best->geometry, best->uncertainty) : 0.0f; bool tilt_refined = refine_tilt && tilt_was_free && significance >= TILT_MIN_SIGNIFICANCE; if (refine_tilt && !tilt_refined) { - if (auto pinned = fit_from(best->geometry.GetDetectorDistance_mm(), true, false)) + // From the header's tilt, not from the one being declined. fit_from now takes a whole geometry, + // so without this the "pinned" refit would pin rot1/rot2 at exactly the unvalidated values the + // gate just rejected - which is the same fault the gate exists to catch, reintroduced one level + // up. + DiffractionGeometry pinned_start = best->geometry; + pinned_start.PoniRot1_rad(geom.GetPoniRot1_rad()).PoniRot2_rad(geom.GetPoniRot2_rad()); + std::optional pinned; + for (auto &a : fit_from(pinned_start, true, false)) + if (!pinned || a.points.size() > pinned->points.size() + || (a.points.size() == pinned->points.size() + && a.rms_radial_pxl < pinned->rms_radial_pxl)) + pinned = std::move(a); + if (pinned) best = std::move(*pinned); } diff --git a/image_analysis/geom_refinement/RingsFromProfile.cpp b/image_analysis/geom_refinement/RingsFromProfile.cpp index 5fae5da15..e885964d5 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.cpp +++ b/image_analysis/geom_refinement/RingsFromProfile.cpp @@ -6,17 +6,9 @@ #include "RingsFromProfile.h" #include "AssignSpotsToRings.h" // RingMatchWindow +#include "PowderAutoSeed.h" // ProfileRingTrack #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); @@ -79,15 +71,13 @@ float SectorPeakQ(const std::vector &profile, int32_t q_bins, int phi_bin return static_cast(sum_wq / sum_w); } -} // namespace - std::vector RingsFromAzimuthalProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, const std::vector &calibrant_ring_q, float q_window_recipA, float min_peak_over_noise, - const std::vector &profile_ring_q) { + const DiffractionGeometry *seeded) { std::vector out; const int32_t q_bins = mapping.GetQBinCount(); @@ -104,35 +94,46 @@ 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 = 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. - // 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; - const int window_bins = static_cast(std::lround(window / q_spacing)); - 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); - // Two background bins at each end and a peak between them is the least this can work with; a - // ring whose window is narrower than that is not resolved at this q spacing. - if (hi_bin - lo_bin < 6) - continue; + // Where to LOOK, ring by ring and sector by sector. Without a seed a ring is looked for at its own + // q in every sector, which is the right answer only when the geometry that binned the profile was + // already close; with one, each ring is tracked through the profile it really made. + const size_t rings = calibrant_ring_q.size(); + std::vector> track(rings); + for (size_t i = 0; i < rings; ++i) { + if (seeded) + track[i] = ProfileRingTrack(calibrant_ring_q[i], *seeded, geom, azim_bins); + else + track[i].assign(azim_bins, calibrant_ring_q[i]); + } + for (size_t i = 0; i < rings; ++i) { for (int phi_bin = 0; phi_bin < azim_bins; ++phi_bin) { + const float q_ring = track[i][phi_bin]; + 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 a fixed window merges into one + // peak. Measured against the neighbours IN THIS SECTOR, since that is where they are here. + float window = q_window_recipA; + if (i > 0 && std::isfinite(track[i - 1][phi_bin])) + window = std::min(window, 0.5f * std::abs(q_ring - track[i - 1][phi_bin])); + if (i + 1 < rings && std::isfinite(track[i + 1][phi_bin])) + window = std::min(window, 0.5f * std::abs(track[i + 1][phi_bin] - q_ring)); + + if (!(q_ring - window > low_q) || !(q_ring + window < high_q)) + continue; + const int window_bins = static_cast(std::lround(window / q_spacing)); + 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); + // Two background bins at each end and a peak between them is the least this can work with; + // a ring whose window is narrower than that is not resolved at this q spacing. + if (hi_bin - lo_bin < 6) + continue; + 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)) diff --git a/image_analysis/geom_refinement/RingsFromProfile.h b/image_analysis/geom_refinement/RingsFromProfile.h index 057b63710..0c0bc58fd 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.h +++ b/image_analysis/geom_refinement/RingsFromProfile.h @@ -35,18 +35,32 @@ // 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. // +// Peak position of one ring in one azimuthal sector of the profile, in q, or NaN where there is no peak +// worth using. Shared because the beam-centre seed measures the same thing about a ring it has already +// found, rather than about one a standard predicts. +// +// 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 is the +// intensity-weighted centroid of everything above half the peak height, which needs no line-shape +// assumption - a powder ring is the instrumental profile convolved with whatever strain and size +// broadening the standard has, not a Gaussian. +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); + // 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. +// seeded, where it is given, is where the geometry is believed to REALLY be - which is not what the +// profile was binned with, and which moves where each ring has to be looked for. That search runs per +// SECTOR (ProfileRingTrack), because a beam centre wrong by (dx, dy) puts a ring at a different q in +// every sector, and one window centred on one q finds it only where that oscillation happens to be +// small. The points still come back labelled with calibrant_ring_q, because that is the value the fit +// has to drive the geometry to; the track only decides where to look. Null means the profile's own +// geometry is the best guess going in. 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, - const std::vector &profile_ring_q = {}); + const DiffractionGeometry *seeded = nullptr); diff --git a/tests/RingsFromProfileTest.cpp b/tests/RingsFromProfileTest.cpp index 0b28cfa54..c18ad9a10 100644 --- a/tests/RingsFromProfileTest.cpp +++ b/tests/RingsFromProfileTest.cpp @@ -231,76 +231,83 @@ TEST_CASE("PowderAutoSeed_RingRadiiDoNotDependOnTheAssumedDistance", "[DetGeomCa } } -// 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. +// Where a ring APPEARS in a profile binned with one geometry, if the truth is 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_RingTrackRoundTripsWhenTheGeometryIsRight", "[DetGeomCalib]") { + DiffractionExperiment x(DetJF4M()); + x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0); + const DiffractionGeometry geom = x.GetDiffractionGeometry(); + constexpr int32_t AZIM_BINS = 32; + for (const float q : {1.0f, 2.0f, 3.0f}) { - CHECK(ProfileQForRing(q, 2.0f * DISTANCE_MM, DISTANCE_MM, WAVELENGTH_A, PIXEL_MM) > q); + const auto track = ProfileRingTrack(q, geom, geom, AZIM_BINS); + REQUIRE(track.size() == AZIM_BINS); + for (const float t : track) + if (std::isfinite(t)) + CHECK(t == Catch::Approx(q).epsilon(1e-3)); } } -// The tilt gate, calibrated against a tilt that is really there. A tilt several rings resolve has to -// clear the threshold comfortably, or the gate would be throwing away real geometry - so this is the -// half of the gate's calibration that the LaB6 series cannot supply, since there the truth is unknown. -TEST_CASE("PowderCalibration_AGenuineTiltClearsTheTiltGate", "[DetGeomCalib]") { +// ...and a beam centre that is wrong makes the ring wander in q ONCE PER TURN, which is the property +// the track exists to follow. A single window centred on one q cannot hold a ring that does this, which +// is why the extraction searches sector by sector when it has a seed to search from. +TEST_CASE("PowderAutoSeed_RingTrackFollowsAWrongBeamCentre", "[DetGeomCalib]") { + DiffractionExperiment x(DetJF4M()); + x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0); + const DiffractionGeometry binned = x.GetDiffractionGeometry(); + DiffractionGeometry truth = binned; + truth.BeamX_pxl(binned.GetBeamX_pxl() + 20.0f); + constexpr int32_t AZIM_BINS = 32; + + const auto track = ProfileRingTrack(2.0f, truth, binned, AZIM_BINS); + float lo = std::numeric_limits::max(), hi = std::numeric_limits::lowest(); + int finite = 0; + for (const float t : track) + if (std::isfinite(t)) { lo = std::min(lo, t); hi = std::max(hi, t); ++finite; } + REQUIRE(finite > AZIM_BINS / 2); + // It has to swing by appreciably more than nothing, or there would be no need to track it... + CHECK(hi - lo > 0.01f); + // ...and the swing has to bracket the ring's own q, since the centre error only moves it. + CHECK(lo < 2.0f); + CHECK(hi > 2.0f); +} + +// The beam-centre offset read off the ring radii alone - no calibrant, no distance. This is the seed +// that lets a header centre further out than the extraction window be recovered at all. +TEST_CASE("PowderAutoSeed_RecoversTheBeamCentreOffset", "[DetGeomCalib]") { DiffractionExperiment x(DetJF4M()); x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0); auto azint = x.GetAzimuthalIntegrationSettings(); - azint.AzimuthalBinCount(64); + azint.AzimuthalBinCount(32); x.ImportAzimuthalIntegrationSettings(azint); + // A beam centre on the detector, not at its corner where the fixture leaves it. The seed reads a + // once-per-turn wobble, so it needs the rings to go round: with the beam in the corner only a + // quarter of the azimuth carries any ring at all and the two components cannot separate. + x.BeamX_pxl(1000.0f).BeamY_pxl(1050.0f); + PixelMask pixel_mask(x); AzimuthalIntegrationMapping mapping(x, pixel_mask); const DiffractionGeometry geom_assumed = x.GetDiffractionGeometry(); DiffractionGeometry geom_true = geom_assumed; - geom_true.PoniRot1_rad(0.02f).PoniRot2_rad(-0.015f); + // A small offset, because the rings this fixture draws are about a pixel wide and its q bins are + // finer than a real run's. Past a few pixels those sharp rings split in the azimuthal average into + // the two HORNS of the sinusoid they trace - density piles up where r(phi) turns round, at R+|d| and + // R-|d| - and the radius finder then reports two rings where there is one. Real powder rings are + // broad enough to smear that out, which is why the offset this recovers on a real LaB6 exposure is + // four times the one it can be shown recovering here. The point of the test is the sign convention + // and the magnitude, which are what a caller would get catastrophically wrong. + geom_true.BeamX_pxl(geom_assumed.GetBeamX_pxl() + 4.0f) + .BeamY_pxl(geom_assumed.GetBeamY_pxl() - 3.0f); const auto profile = SynthesiseProfile(mapping, geom_assumed, geom_true); - const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom_assumed, LAB6_RINGS); - REQUIRE(rings.size() > 60); + const auto observed = RingRadiiFromProfile(profile, mapping, geom_assumed); + REQUIRE(observed.size() >= 2); - RingFitUncertainty unc; - const auto fitted = RingOptimizer(geom_assumed).Run(rings, &unc); - REQUIRE(unc.valid); - CHECK(TiltSignificance(fitted, unc) > TILT_MIN_SIGNIFICANCE); -} - -// ...and the other half: a tilt the fit did not have as a free parameter has no significance at all, -// which is what makes it impossible for one to be reported as measured. The case this really guards is -// subtler than --no-refine-tilt - RingOptimizer pins the tilt by itself whenever every point it is given -// lies on one ring, so a fit can arrive here with a tilt inherited from an earlier pass and a sigma of -// zero. Reading zero significance as "decline and refit pinned" is what keeps that out of the answer. -TEST_CASE("PowderCalibration_APinnedTiltHasNoSignificance", "[DetGeomCalib]") { - DiffractionExperiment x(DetJF4M()); - x.QSpacingForAzimInt_recipA(0.004).QRangeForAzimInt_recipA(0.5, 4.0); - auto azint = x.GetAzimuthalIntegrationSettings(); - azint.AzimuthalBinCount(64); - x.ImportAzimuthalIntegrationSettings(azint); - - PixelMask pixel_mask(x); - AzimuthalIntegrationMapping mapping(x, pixel_mask); - - DiffractionGeometry geom = x.GetDiffractionGeometry(); - const auto profile = SynthesiseProfile(mapping, geom, geom); - const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom, LAB6_RINGS); - REQUIRE(!rings.empty()); - - // Carry a tilt in, and pin it. The geometry that comes out still has that tilt in it - nothing - // removed it - but the fit never measured it, and the significance has to say so. - geom.PoniRot1_rad(0.02f); - RingFitUncertainty unc; - const auto fitted = RingOptimizer(geom, /*refine_tilt=*/false).Run(rings, &unc); - CHECK(fitted.GetPoniRot1_rad() == Catch::Approx(0.02f)); - CHECK(unc.sigma_rot1_rad == 0.0); - CHECK(TiltSignificance(fitted, unc) == 0.0f); - CHECK(TiltSignificance(fitted, unc) < TILT_MIN_SIGNIFICANCE); + const auto offset = BeamCentreOffsetFromProfile(profile, mapping, geom_assumed, observed); + REQUIRE(offset.has_value()); + CHECK(offset->first == Catch::Approx(4.0).margin(1.5)); + CHECK(offset->second == Catch::Approx(-3.0).margin(1.5)); }