Add pinned polymer diffraction peak inference - #267
Draft
NJMarchese wants to merge 21 commits into
Draft
Conversation
Wire the existing quantem calibration tools into preprocess() (previously a
stub): image-center finding, ellipticity fitting, descan (CoM plane fit), and
detector-rotation estimation. Follows the class's lazy design -- parameters are
measured and cached (image_centers, ellipse_params, descan_origin,
detector_rotation_deg, sampling_inv_A), then consumed downstream by the polar
transforms rather than re-warping the raw 4D data.
- Ellipticity via fit_probe_ellipse on the mean DP -> metadata["ellipticity"]
- Descan + rotation via CenterOfMassOriginModel (calculate_origin,
fit_origin_background, estimate_detector_rotation)
- image_centers via find_central_beams_4d ("descent"/"grid"/"peaks") or the
CoM/descan field, selectable with center_source
- Register the new calibration attributes in __init__
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Render a boustrophedon cursor walk over the scan to an animated GIF: each frame pairs the real-space intensity map (cursor crosshair) with that position's diffraction pattern and detected Bragg-peak overlay, reusing the save_peak_figures rendering primitives. PIL-based writer.
…spacing axis - estimate_peak_windows(mode='intensity'|'count', log_scale=False): detect peak windows on the peak-count radial profile as well as intensity, and optionally on log1p(profile) so small peaks are not dominated by large ones. peak_info now carries 'mode', 'log_scale', and 'profile'. - peak_radial_intensity_plot / peak_radial_count_plot: add log_scale (log y-axis with a positive fill baseline) and show_d_spacing (top axis in real-space d-spacing (Å) = 1/q) when plotting. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- estimate_peak_windows: printed per-peak summary now reports d-spacing (Å) = 1/q alongside the q values (center and window). - plot_peak_count_map: panel titles now include the d-spacing range (Å) under the q range (1/Å). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Reorder BraggPeaksPolymer.preprocess() so calibration runs descan/rotation -> centering -> ellipticity (was ellipse first). Fitting the ellipse on the raw mean DP smeared the diffraction ring by the descan drift and biased the fit toward the central beam. The ellipse is now fit on a centered mean DP built by the new _centered_dp_mean() helper: when a fitted CoM model is available it reuses shift_origin_to()'s sub-pixel grid-sampler to align every pattern to the detector center on-device, otherwise it falls back to translating the plain mean DP by the average center offset. Centering (find_central_beams_4d) now runs without ellipse_params since the ellipse is not yet known. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Replace fit_probe_ellipse (which Otsu-thresholds and fits the brightest blob = the central beam, so it measured probe shape and latched onto a smeared/off-center beam) with a diffraction-RING fit using Karen Ehrhardt's angular-uniformity criterion. New _fit_ellipse_from_ring() holds the found center fixed and searches (b/a, theta) to minimise the azimuthal variance of a ring annulus in the polar transform (quantem.diffraction.polar_transform). It samples out at the ring radius and never touches the central beam, so a doubled / off-center / drift-smeared central beam no longer biases the ellipse. Details: - ring band auto-detected from the circular radial profile (skip central beam to first trough, take strongest ring beyond), or set explicitly via new ellipse_radial_min / ellipse_radial_max preprocess params. - coarse (b/a, theta) grid + local refine; output canonicalised to a/b >= 1. - centered mean DP now cached on self.dp_mean_centered (was ephemeral); preprocess caches self.ellipse_ring_band; show=True plots the DP plus before/after polar so the ring flattening is visible. - ellipse_threshold kept but unused (back-compat). Validated against polar_transform on synthetic elliptical rings (a/b and theta recovered to <0.5%). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
_centered_dp_mean used integer torch.round rolls. The per-pattern rounding residual (up to 0.5 px) does NOT average out when the origin spread is narrow -- every pattern rounds to the same integer, so the fractional part becomes a constant bias and the mean beam lands off the detector center by a fraction of a pixel (visible as the central beam sitting up/right of the geometric-center crosshair, and it also throws off the fixed-center ring ellipse fit). Replace integer rolls with bilinear splatting: each pattern's continuous shift (center - fitted origin) is split into floor + fraction, and its intensity is distributed over the 4 integer-shift corners with bilinear weights, grouped by floor shift for efficiency. Same low memory (one (Qy,Qx) accumulator + one batch), now centers to <0.01 px. Also fix the no-CoM fallback to average image_centers over valid (non-zero) positions only -- it is 0 outside the scan mask, which otherwise biases the fallback center toward the origin. Validated on synthetic beams: bilinear hits the detector center exactly at every origin spread; integer rolls drift up to ~0.4 px at narrow spread. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The centered mean DP still showed the central beam a few px off the detector center. Root cause: _centered_dp_mean shifted each pattern by the descan CoM origin (origin_fitted). The CoM is the centroid of the WHOLE pattern, so ring/background asymmetry (e.g. a bright corner) pulls it several px off the actual beam -- centering by it puts the CoM at the center and leaves the beam off. Center instead by image_centers, the angular-uniformity BEAM center (found by ring symmetry, immune to intensity asymmetry, and the center everything downstream already uses). Only ROI (in-mask, non-zero) patterns are averaged. Still sub-pixel bilinear, same low memory. Validated on a synthetic pattern with an asymmetric background pulling the CoM to (130,132): centering by CoM leaves the beam at (114.5,114.3) (~7 px off); centering by image_centers lands it at (107.9,107.9) on the (107.5,107.5) detector center. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The descan residual (origin_com_measured - descan_origin) showed a uniform ~1-2 px offset inside the ROI. Cause: CenterOfMassOriginModel fits fit_origin_background over the WHOLE scan, so out-of-ROI patterns (vacuum/substrate, meaningless CoM) drag the global plane and leave a constant offset inside the mask; hot/dead pixels contaminate it too. CenterOfMassOriginModel has no mask support, so preprocess now fits the plane itself over ROI patterns via new _fit_origin_roi (ordinary LS z = a + b*row + c*col per component, with sigma-clip passes to reject outlier CoM), then assigns com_model.origin_fitted and zeroes the residual outside the ROI (measured := fitted there) so estimate_detector_rotation's curl isn't contaminated by junk patterns either. Falls back to the full scan when scan_mask is None. Other centering steps already respect the ROI: find_central_beams_4d takes scan_mask, and the ellipse fit runs on the image_centers-centered (ROI-only) mean DP. Validated: with out-of-ROI junk + hot-pixel outliers, the global fit leaves a -13 px uniform ROI residual; the ROI fit brings it to -0.005 px and recovers the true plane exactly. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
save_cartesian_peaks / save_polar_peaks / save_polar_data / save_peak_intensities did np.save(vector). A canon Vector is array-like, so numpy flattened it into an (Ry,Rx) object array of cells, losing the Vector's fields/units. On load, .item() only unwraps a size-1 array, so you got back a size-(Ry*Rx) object array instead of a Vector, and polar_transform_peaks raised "can only convert an array of size 1 to a Python scalar". Add _save_object() which wraps the object in a 0-d object array so np.save pickles the whole Vector; the existing size-1/.item() unwrap in the load_* methods then restores it. Route all four save_* through it. Note: files saved by the old code are already flattened and can't be reconstructed (metadata lost) -- re-run find_peaks_model and re-save. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
calculate_orientation_correlation sized its angular-mode batches from torch.cuda.mem_get_info(), which reports memory free on the *device* and is blind to torch.cuda.set_per_process_memory_fraction. Under a cap it budgeted headroom the process was not allowed to allocate and then died in fft2 with the card still mostly free: a 23.74 GiB cap against 60.65 GiB device-free raised OutOfMemoryError on a 100x100 scan at orientation_upsample=8. Headroom is now min(device_free, cap - memory_allocated). It is measured against *allocated*, not *reserved*: cached blocks no tensor is using are reusable, and after a failed run they can occupy most of the cap, which would otherwise report a zero budget and refuse to run. Also reserve room for theta_spectrum, which is allocated after the estimate and was previously unaccounted, and retry once after empty_cache() before giving up, so a failed run no longer poisons the next attempt. The MemoryError now reports the cap, live usage, spectrum reservation and device-free rather than a bare "0.00 GiB". Tests: tests/diffraction/test_orientation_correlation.py, 10 passed. Numerics unchanged, max relative difference 8.96e-08. Both failure modes reproduced on GPU and confirmed fixed: a cap saturated with stale cache now completes at 2.74 GiB under a 7.60 GiB cap, and a 1.14 GiB cap degrades to 182 small batches instead of failing. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Replaces the central-beam/probe-blob ellipse fit with a ridge fit on the diffuse ring. It measures the ring radius around every azimuth, jointly refines the residual center and the ellipse, and accepts the correction only when held-out angular sectors improve over a plain circle. Because it never touches the central beam, descan-drift smearing and doubling of the beam no longer bias the ellipse. Adds _fit_ellipse_from_ridge with polar_at / extract_ridge helpers and ellipse and circle residual models, wires it into preprocess() via ellipse_fit_method, and records the outcome in ellipse_fit_diagnostics (method, accepted, selected) alongside ellipse_ring_band. Tests: tests/diffraction/test_ellipse_ring_fit.py plus tests/diffraction/test_origin_finding.py, 21 passed. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The correlation-equals-one boundary saturates: on real data (pg3T2 07-03-2024 scan 61) ring pair (2,2) climbs from 33.8 degrees at zero separation to 71.6 by mid-lobe, then flattens to 78.9. Fitting an unweighted straight line over the whole lobe let that flat tail dominate, which displaced the intercept off the measured boundary and biased the slope low. The fitted intercept landed up to 6.99 degrees away, visibly inside the blue region rather than on the gray baseline. Weight the fit exponentially towards short distances so the reported slope is the near-origin tangent, which is the physically meaningful quantity. Intercept error against the measured boundary at zero separation, across all six ring pairs: before +0.66 +3.14 -6.99 +4.70 -3.66 +6.89 (max 6.99) after -1.11 +1.40 -1.69 -1.09 +1.31 -1.39 (max 1.69) The 1/e decay defaults to SLOPE_WEIGHT_FRACTION (0.10) times the fitted distance span, chosen by sweeping the fraction against this dataset. slope_weight_scale=numpy.inf restores the previous unweighted fit exactly. R-squared is now weighted with the same weights and is not comparable to the old value. Adds slope_fit_intercept_degrees, slope_fit_effective_point_count (Kish) and slope_weight_scale to the metrics so the fit can be checked numerically. The fit line is drawn only out to three decay lengths, past which the boundary has saturated away from the tangent. Note this changes reported slopes, which were biased low: pair (2,2) 0.231 -> 0.371 deg/px (+61%), (1,1) +49%, (1,2) +46%, (0,2) -20%. Previously exported slope values need regenerating. Tests: tests/diffraction/test_orientation_correlation.py, 10 passed. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The polymer ice flagger only tested "bright and in the q band, on one
six-fold lattice". Three additions, each off by default:
Sharpness gate. Ice reflections are small and sharp; polymer signal is a
larger dot or a broad diffuse region. Peaks carry no width, so measure it
from the polar volume: radial/annular FWHM above a local baseline, at each
candidate's (r, theta). max_width_r_invA / max_width_theta_deg gate the q
band before the lattice search, so broad peaks neither get flagged nor drag
the phi estimate. sharpness_mode picks intersection ("both", compact dots)
or union ("either", which also keeps thin streaks). sharpness_mask is public
so a tuning preview applies the same rule instead of reimplementing it, and
collect_peak_widths / measure_ice_peak_widths report the widths to tune from.
Multiple crystallites. A pattern can hold several crystallites at unrelated
orientations; one pass only ever saw the strongest. The matcher now peels:
claim the best-supported lattice, remove its peaks, look again, up to
max_crystallites. min_phi_separation_deg stops a single lattice being
re-found as a near-duplicate. IceFlaggerDebug.phi_deg lists what was found.
Folded theta. process_polar(two_fold_symmetry=True) folds theta to [0,180),
collapsing every Friedel pair onto one angle -- so a pair scored one bin and
was rejected by min_matches=2, and arms 3-5 were unreachable, making
min_matches>3 silently impossible. The matcher is now period-aware
(theta_period_deg, read off the BraggPeaksPolymer), and detect_ice raises on
an unsatisfiable min_matches rather than matching nothing. min_peaks_per_arm
recovers the distinction folding destroys: 2 demands a Friedel pair on an
arm and rejects a lone peak.
Also: plot_q_intensity_density shades the candidate region as a rectangle
instead of drawing bare window edges.
Tests: 77 pass in tests/diffraction (28 new in test_polymer_ice.py).
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Folding theta to [0,180) for two-fold symmetry discarded which half of the circle a peak came from, so a Friedel pair became indistinguishable from two peaks that merely sit close together. Peaks are small, so keep the full angle: polar_transform_peaks now emits a "theta_unfolded" column alongside the folded "theta". Nothing else reads it by position -- lookups go through fields.index() -- so the extra column is transparent. IceFlaggerParams.require_friedel_pair uses it: an arm counts only if it holds two peaks whose unfolded angles differ by 180 +/- dtheta_deg. That is the test min_peaks_per_arm can only approximate, since a count of two is also satisfied by two neighbours within dtheta_deg. The field is optional -- detect_ice raises a clear "re-run polar_transform_peaks" error only when the strict test is requested on a vector that predates it. Tests: 80 pass in tests/diffraction. test_bragg_peak_polar_transform_inverts_ ellipse_mapping updated for the new column, and now pins the field order. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Follow-up to the sharpness gate, from working it against a real polymer scan. Annular width is the discriminator; drop sharpness_mode. Ice is annularly sharp both as compact dots and as radial streaks (narrow in theta, extended in q), while polymer at the same q is an annularly broad arc. The radial axis does not separate ice from polymer -- dots and arcs measure alike there -- it only separates streaks from dots, so a radial ceiling costs the streaks and buys nothing. That makes the "both"/"either" combining rule pointless: with one ceiling the modes are identical, and "either" silently defeated a tight annular ceiling for any radially sharp peak. Ceilings are now simply ANDed. Width measurement is noise-robust. The half-maximum walk stopped at the first sample below half, so one downward fluctuation on a weak diffuse arc ended it early: a true 40 degree arc measured 10-30 degrees at realistic SNR, i.e. as narrow as ice. The profile is now Gaussian-smoothed (sigma 1 bin) and the crossing must persist for 3 samples. The kernel is removed in quadrature -- exact for a Gaussian convolved with a Gaussian, which is why it is a Gaussian and not a boxcar; a boxcar leaves a residual bias on sharp peaks. Median over 8 seeds, true 40 degrees at SNR 2: 10.5 -> 43.6, while a true 5 degree peak stays at 6.0. Also widened the annular window to 90 degrees and dropped the baseline quantile to 0.05, since a window the feature fills puts the baseline partway up it and saturates the measurement near 39 degrees. Parameters validate on construction. Every confusion this cost came back as "nothing was flagged" rather than an error. Module docstring now states the four criteria and the folding constraints. Tests: 89 pass in tests/diffraction. The synthetic polar fixture sampled q at 0.025 1/A per bin, too coarse for the default radial window; now 0.005. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…ce flagger Choosing IceFlaggerParams against a real dataset needs plots, not just the algorithm: a selection box over the sharpness ceilings, ice-band orientation histograms for all/kept/removed peaks, a widget view of either side of the split, and a per-peak FWHM probe that shows the cuts each width is measured from. These started as notebook cells, which meant every scan directory carried its own drifting copy; in the package any notebook can import them. Free functions taking (bp, params) explicitly rather than methods on BraggPeaksPolymer, which is already 79 methods. The algorithms stay in polymer_ice; this is only the interactive layer. Nothing mutates bp -- IcePeakView delegates to it so the widget can show a peak subset without the caller's object being altered. The quantem.widget import is deferred into the one function that needs it, so the module imports without that package. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
BraggPeaksPolymer, peak detection, angular-origin/polar transforms, and calibrated flowline/figure helpersdev/polymer overlap reconciliationThe default paper record remains intentionally unavailable for network download until its immutable public Zenodo record exists. This PR should remain draft until that URL is pinned.
Verification
pytest -p no:cacheprovider -q tests/diffraction tests/datastructures/test_dataset4dstem.py(32 passed)(1, 2, 256, 256)inference passedScope exclusions
Grain clustering, generator/training launchers, research notebooks, generated results, and the private Zenodo archive are not included.