diff --git a/image_analysis/geom_refinement/PowderCalibration.cpp b/image_analysis/geom_refinement/PowderCalibration.cpp index adfbb5192..a9f514ac5 100644 --- a/image_analysis/geom_refinement/PowderCalibration.cpp +++ b/image_analysis/geom_refinement/PowderCalibration.cpp @@ -10,7 +10,7 @@ #include "PowderCalibration.h" #include "../../common/GitInfo.h" #include "../../common/JFJochMath.h" -#include "AssignSpotsToRings.h" +#include "AssignSpotsToRings.h" // FindCircleCenter #include "RingOptimizer.h" #include "RingsFromProfile.h" @@ -73,7 +73,8 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, const std::vector &calibrant_ring_q, - bool refine_tilt) { + bool refine_tilt, + const std::vector &spots) { // 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 @@ -90,15 +91,47 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, // 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); + // Three places a beam centre can come from, and they fail in different regimes - which is the whole + // reason to carry all of them. The header is right on a well-configured instrument and is what a + // calibration is run to check. The rings' own wobble is exact while the error stays under about half + // a ring spacing, and stops meaning anything beyond that, where each ring reaches the azimuthally + // averaged profile as two horns rather than one peak. The circle through the spots reads nothing + // from the header whatsoever, so it holds where both of the others have given up - measured, it + // returns the same answer from a header 400 px out in the centre and eight times out in distance. + // None of them is trusted; each is fitted and the residual chooses. + std::vector centres = {geom}; + const auto add_centre = [&](float x, float y) { + // Within a pixel of one already on the list it is not another hypothesis, it is the same one, + // and fitting it again would only double the work for two answers nothing can tell apart. + for (const auto &c : centres) + if (std::hypot(x - c.GetBeamX_pxl(), y - c.GetBeamY_pxl()) < 1.0f) + return; + DiffractionGeometry candidate = geom; + centres.push_back(candidate.BeamX_pxl(x).BeamY_pxl(y)); + }; + + if (const auto offset = BeamCentreOffsetFromProfile(profile, mapping, geom, observed)) + add_centre(geom.GetBeamX_pxl() + offset->first, geom.GetBeamY_pxl() + offset->second); + + // What the spots make of it, as a COMPLETE geometry rather than only a centre. Taking the circle + // centre alone is not enough: the distance candidates are read off the azimuthally averaged profile, + // which is averaged about the header's centre, and once that is badly wrong the profile no longer + // shows rings at all - a ring smeared over hundreds of pixels cannot be un-smeared by reading its + // bins differently. The spots have no such problem, because a spot's position is a fact about the + // image; GuessGeometry votes for the circle centre, clusters the radii into rings and takes the + // distance from the innermost one, which is exactly the path that makes --calibration spots immune + // to the header. Let it produce one whole starting geometry of its own. + std::optional from_spots; + if (spots.size() >= 3) { + try { + DiffractionGeometry guess = geom; + GuessGeometry(guess, spots, calibrant_ring_q, refine_tilt); + from_spots = guess; + } catch (const std::exception &) { + // No circle, or no ring clusters - the spots simply have nothing to say here. The profile's + // own hypotheses stand on their own, so there is nothing to report and nothing to stop. + } + } // 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 @@ -195,14 +228,17 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, 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; + // The spots' geometry is a start in its own right, not one to cross with the profile's distances - + // its centre and its distance are measured together and belong together. + if (from_spots) + starts.push_back({*from_spots, true, {from_spots->GetDetectorDistance_mm(), 0.0}}); + for (size_t centre = 0; centre < centres.size(); ++centre) { for (size_t i = 0; i <= candidates.size(); ++i) { const bool seeded_distance = i < candidates.size(); - DiffractionGeometry start = base; + DiffractionGeometry start = centres[centre]; if (seeded_distance) start.DetectorDistance_mm(candidates[i].distance_mm); - starts.push_back({start, seeded_distance || centre == 1, + starts.push_back({start, seeded_distance || centre > 0, {seeded_distance ? candidates[i].distance_mm : 0.0f, seeded_distance ? candidates[i].score : 0.0}}); } @@ -275,6 +311,14 @@ CalibrationResult CalibrateFromProfile(const std::vector &profile, } auto result = Summarize(best->geometry, best->points, best->uncertainty); + if (from_spots) { + result.spots_available = true; + result.spots_beam_x_pxl = from_spots->GetBeamX_pxl(); + result.spots_beam_y_pxl = from_spots->GetBeamY_pxl(); + result.spots_disagreement_pxl = + std::hypot(from_spots->GetBeamX_pxl() - best->geometry.GetBeamX_pxl(), + from_spots->GetBeamY_pxl() - best->geometry.GetBeamY_pxl()); + } result.tilt_refined = tilt_refined; result.tilt_significance = significance; result.seed_distance_mm = best_provenance.seed_mm; diff --git a/image_analysis/geom_refinement/PowderCalibration.h b/image_analysis/geom_refinement/PowderCalibration.h index c4d4cf392..16769e5d9 100644 --- a/image_analysis/geom_refinement/PowderCalibration.h +++ b/image_analysis/geom_refinement/PowderCalibration.h @@ -45,6 +45,17 @@ struct CalibrationResult { // the tilt was never a free parameter. bool tilt_refined = false; float tilt_significance = 0.0f; + // How far the beam centre the SPOTS vote for is from the one the rings were fitted to, in pixels, + // and where the spots put it. An independent cross-check that costs nothing, because the spots have + // already been found: the circle through them reads nothing from the header, so it holds in exactly + // the regime the profile does not - a profile binned about a badly wrong centre shows each ring + // smeared across its sectors, and no seeding recovers that. A large disagreement is therefore not a + // close call to arbitrate but a statement that this profile could not have been fitted, whatever + // the residual says. Zero when no spots were available. + float spots_beam_x_pxl = 0.0f; + float spots_beam_y_pxl = 0.0f; + float spots_disagreement_pxl = 0.0f; + bool spots_available = false; }; // How many of its own sigmas the fitted tilt is away from zero - the number the gate above reads. @@ -69,11 +80,20 @@ constexpr float TILT_MIN_SIGNIFICANCE = 3.0f; // given a geometry measured with the tilt pinned than one measured tilted and then flattened. // Fit the geometry to the rings found in a run-summed azimuthal profile (CalibrationMethod::Rings). +// +// spots, where they are given, add one more beam-centre hypothesis: the centre of the circle through +// them, voted for by every triple (FindCircleCenter). It reads nothing from the header at all - only +// where the spots landed - so it is the one estimate that survives a header the profile itself cannot +// correct from. The rings' own estimate stops working once the beam centre is out by more than about +// half a ring spacing, because past that the azimuthally averaged profile shows each ring as the two +// horns of a sinusoid rather than as a ring; the circle through the spots has no such limit. Measured +// on a 110 mm LaB6 exposure it returns the same geometry from a header 400 px and 8x in distance wrong. CalibrationResult CalibrateFromProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, const std::vector &calibrant_ring_q, - bool refine_tilt = true); + bool refine_tilt = true, + const std::vector &spots = {}); // Fit the geometry to a pooled spot list (CalibrationMethod::Spots): the beam centre from scratch off // the Hough circle centre, then the same ring fit. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 90a6e41bb..98be31ba6 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1611,12 +1611,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } const bool full = (config_.mode == ProcessMode::FullAnalysis); - // Powder calibration by spots finds spots and nothing else, so it runs the same per-image engine as - // the full analysis (spot finding lives inside MXAnalysisWithoutFPGA) with indexing switched off in - // its settings; by rings it only needs the azimuthal profile, which is the azint worker exactly. + // Powder calibration runs the same per-image engine as the full analysis (spot finding lives inside + // MXAnalysisWithoutFPGA) with indexing switched off in its settings. By rings it needs only the + // azimuthal profile and the azint worker would do - but the spots are worth their half second even + // then, because the circle through them places the beam centre without reference to the header at + // all, which is the one hypothesis that survives a header the profile cannot correct. const bool calibration = (config_.mode == ProcessMode::Calibration); const bool calibration_spots = calibration && config_.calibration_method == CalibrationMethod::Spots; - const bool per_image_analysis = full || calibration_spots; + const bool per_image_analysis = full || calibration; const bool write_files = write_output && !config_.output_prefix.empty(); // Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization, @@ -2812,7 +2814,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b msg.run_number = experiment_.GetRunNumber(); msg.run_name = experiment_.GetRunName(); - if (calibration_spots) { + if (calibration) { std::lock_guard lock(calibration_spot_m); for (size_t i = 0; i < msg.spots.size() && i < calibration_spots_per_image; ++i) calibration_spot_list.push_back(msg.spots[i]); @@ -2926,7 +2928,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.calibration = CalibrateFromProfile(summed->GetResult(), mapping, experiment_.GetDiffractionGeometry(), config_.calibrant_ring_q, - config_.calibration_refine_tilt); + config_.calibration_refine_tilt, + calibration_spot_list); } } diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index f9d271808..711dc1393 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2092,6 +2092,18 @@ static int RunRugnux(int argc, char **argv) { const auto [beam_x, beam_y] = g.GetDirectBeam_pxl(); std::cout << fmt::format("Direct beam: {:.3f}, {:.3f} px", beam_x, beam_y) << std::endl; + // Where the spots put the beam, independently of everything above. The circle through them + // reads nothing from the header, so it holds in exactly the regime the summed profile does not: + // a profile binned about a badly wrong centre shows each ring smeared across its sectors, and no + // amount of seeding recovers that - the fit then converges, prettily, on the wrong geometry. + // Printed as a fact rather than gated on, because the reader can see at a glance whether two + // methods that share no assumption agree. + if (cal.spots_available) + std::cout << fmt::format("Spot cross-check: the circle through the spots puts the beam at " + "{:.1f}, {:.1f} px - {:.1f} px from the fitted centre", + cal.spots_beam_x_pxl, cal.spots_beam_y_pxl, + cal.spots_disagreement_pxl) << std::endl; + // What the fit says about itself. The scatter line above is about the MEASUREMENTS; this is // about the PARAMETERS, and the two disagree exactly where the calibration is worth doubting: // a fit with few rings can sit tightly on the points it has while leaving the tilt free to