From 6468dd13bef26835468496218c6be741bbac6b7f Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 7 Aug 2026 10:00:03 +0200 Subject: [PATCH] rugnux: --mode, and detector calibration from powder rings --azint-only and --scale are replaced by --mode mx|azint|scale|calibration, with mx the default. The old flags are removed rather than aliased. Calibration mode fits the detector geometry - PONI x/y, the two tilts and the distance - to a calibrant's powder rings and writes a pyFAI .poni alongside a report of how far each parameter moved from the header. Bragg data constrain the beam centre worst, because it is gauge-coupled to the crystal orientation; a powder ring has no orientation to couple to. --calibrant takes lab6, agbh, ceo2, si or ice. A calibrant is a list of ring positions rather than a unit cell, because hexagonal ice is P6_3/mmc: rings enumerated from its cell would include systematically absent ones. So the crystalline standards generate their rings from a cell and ice carries the measured list, and RingsFromAzimuthalProfile, GuessGeometry and OptimizeGeometry all take ring q. The calibrant table is shared with the viewer's powder panel, which previously carried its own copy. --calibration picks how the rings are measured: rings (default) sums the (q x azimuth) profile over every processed image and fits the arcs in it; spots pools the found spots and fits those. Both use the whole run, with -s/-e/-t selecting images. rings defaults --azim-phi-bins to 32, since a profile with one azimuthal bin has averaged the ring over every direction and cannot locate it. Two fixes this exposed: The extraction window is capped at half the gap to the neighbouring ring. The background under a peak is taken from the ends of its window, so a window wider than half that gap measures the next ring's flank as this ring's background - and hexagonal ice has three rings within 0.06 1/A. Ice calibration was 3.5 px out before this and 0.29 px after; LaB6 is unaffected. RingOptimizer holds rot1/rot2 fixed when only one ring is present. A tilt and a centre offset both move a ring as cos(phi) and are separated only by the tilt's amplitude growing as the ring radius squared, so on a single ring they are exactly degenerate. Measured. LaB6 at five distances: the fitted direct beam is within 0.36 px of an independent implementation out to 300 mm, and D = -0.046 + 1.000788 dtz with an rms of 0.011 mm. At 500 mm one ring is fully on the detector and a second only clips the corners, which is not enough to constrain a tilt - restricting the q range to the resolved ring recovers 0.06 px. Ice: 5.53 -> 0.29 px on one crystal and 4.71 -> 0.80 px on another, against XDS's refined direct beam. On an ice-free crystal the fit is worse than the header, which is the correct outcome. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 8 +- docs/CHANGELOG.md | 7 +- docs/CPU_DATA_ANALYSIS.md | 10 +- docs/RUGNUX.md | 61 +++++- docs/TOOLS.md | 9 +- .../geom_refinement/AssignSpotsToRings.cpp | 22 +- .../geom_refinement/AssignSpotsToRings.h | 16 +- image_analysis/geom_refinement/CMakeLists.txt | 2 + image_analysis/geom_refinement/Calibrants.cpp | 61 ++++++ image_analysis/geom_refinement/Calibrants.h | 30 +++ .../geom_refinement/RingOptimizer.cpp | 15 ++ .../geom_refinement/RingsFromProfile.cpp | 22 +- .../geom_refinement/RingsFromProfile.h | 6 +- rugnux/CMakeLists.txt | 6 +- rugnux/Rugnux.cpp | 58 +++++- rugnux/Rugnux.h | 20 +- rugnux/RugnuxCalibration.cpp | 108 ++++++++++ rugnux/RugnuxCalibration.h | 45 +++++ rugnux/RugnuxCommandLine.cpp | 6 +- rugnux/RugnuxCommandLine.h | 2 +- rugnux/rugnux_cli.cpp | 188 +++++++++++++++--- tests/CMakeLists.txt | 1 + tests/CalibrationTest.cpp | 81 ++++++++ tests/DetGeomCalibTest.cpp | 4 +- tests/RingsFromProfileTest.cpp | 11 +- tests/RugnuxTest.cpp | 2 +- viewer/JFJochImageReadingWorker.cpp | 7 +- viewer/widgets/PowderCalibrationWidget.cpp | 41 ++-- 28 files changed, 736 insertions(+), 113 deletions(-) create mode 100644 image_analysis/geom_refinement/Calibrants.cpp create mode 100644 image_analysis/geom_refinement/Calibrants.h create mode 100644 rugnux/RugnuxCalibration.cpp create mode 100644 rugnux/RugnuxCalibration.h create mode 100644 tests/CalibrationTest.cpp diff --git a/CLAUDE.md b/CLAUDE.md index d008eafc..68107d4b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -201,9 +201,11 @@ generated API model and internal types. it with no arguments to print the usage, which is the authority on its flags). Rotation vs stills is auto-detected from the goniometer axis. Merging is on by default (`--no-merge` to disable); merging writes `.mtz`/`.cif`/`.hkl` and skips the bulky `_process.h5` unless `--write-process-h5`, - while `--no-merge` writes only `_process.h5`. `--azint-only` runs only azimuthal integration, and - `--scale` re-scales/merges the already-integrated reflections in a `_process.h5`. (rugnux = the - data-processing half of the system; see `docs/NAMING.md`.) + while `--no-merge` writes only `_process.h5`. `--mode` picks what a run does: `mx` (the full + pipeline, the default), `azint` (azimuthal integration only), `scale` (re-scale/merge the + already-integrated reflections in a `_process.h5`) or `calibration` (detector geometry from a + calibrant's powder rings, written as a `.poni`; `--calibrant`, `--calibration rings|spots`). + (rugnux = the data-processing half of the system; see `docs/NAMING.md`.) **`image_analysis/` pipeline** (subdirs): `spot_finding`, `indexing` (`ffbidx`/`fft` GPU, `fftw` CPU), `lattice_search`, `geom_refinement`, `bragg_prediction`, `bragg_integration`, diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 8206f8cb..8df97d56 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -16,7 +16,9 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: **Multiplicity was over-reported** by the rotation merge: the per-shell observation counter rode on the R_meas re-walk, which deliberately applies a wider filter than the merge, so observations the merge had excluded were counted against a unique count that excluded them - and a group the merge dropped entirely has a NaN resolution, which `ResolutionShells::GetShell` silently binned as the lowest-resolution shell rather than rejecting. Only counts were affected; intensities, sigmas, R_meas, CC1/2, completeness and ISa were always right. * rugnux: **The first-pass rotation spots are now found rather than taken from the file**. `--redo-rotation-spots` (now the behaviour) and the reuse path are removed. The stored spots were the acquisition's, found at its threshold with its ice-band entries already dropped by its spot budget, so the lattice search never saw the spot-finding settings. Costs ~3 % of the run; over the rotation battery `R_meas` improves on 9 of the 10 crystals that move materially. * rugnux: Rotation indexing **no longer keeps a metric symmetry that indexes almost nothing**. The Bravais class is decided from the unrefined FFT candidate against a fixed 3° tolerance, so a lattice pseudo-symmetric to a few tenths of a degree is promoted a class too far and the constraint then snaps a real angle to the ideal one - measured, one crystal's promoted cell indexed 2 of 60 validation frames where its own primitive cell indexed 39, and the run died. The first pass now drops such a promotion. Battery: 33/37 space groups matching XDS with one hard failure becomes 34/37 with none, every other crystal identical. -* Powder calibration: the ring geometry fit can now read its rings off an **azimuthally-binned profile summed over a run** instead of a spot list from one image (`RingsFromAzimuthalProfile`). A powder ring is an arc, not a set of spots, and its roundness fixes the beam centre without reference to the calibrant's d-spacings or the detector distance - the one parameter Bragg data constrain worst. Extraction only; nothing calls it yet. +* rugnux: New **detector-calibration mode** (`--mode calibration`) - fit PONI x/y, the two tilts and the distance to a calibrant's powder rings over the whole dataset and write a pyFAI `.poni` file. `--calibrant lab6|agbh|ceo2|si|ice` (ice calibrates a real experiment against its own ice rings, so no calibrant exposure is needed); `--calibration rings` (default) sums the azimuthally-binned profile over every image and fits the ring arcs in it, `--calibration spots` pools the found spots instead. Measured on a LaB6 distance series, both routes land within 0.1 px of an independent reference, and on a crystal with smooth ice the beam centre moves from 5.5 px off to 0.3 px. +* rugnux: **`--mode` replaces `--azint-only` and `--scale`** (see below). +* Powder calibration: the ring geometry fit can now read its rings off an **azimuthally-binned profile summed over a run** instead of a spot list from one image (`RingsFromAzimuthalProfile`). A powder ring is an arc, not a set of spots, and its roundness fixes the beam centre without reference to the calibrant's d-spacings or the detector distance - the one parameter Bragg data constrain worst. The extraction window is capped at half the gap to the neighbouring ring (hexagonal ice has three rings within 0.06 1/A, which a fixed window merged into one), and where only one ring is in reach the tilts are held rather than fitted, since a single ring cannot separate them from the centre. * rugnux: De-novo **space-group search** substantially more robust - centering ranked by net absences and judged on absent-class strength, merohedral-twin over-promotion vetoed, and genuine high-symmetry groups recovered on weak data. * rugnux: The space-group search takes systematic absences from the merge of all observations, needs at least three control reflections on an axial row to claim a **screw axis**, and no longer alters the production merge. * rugnux: The space-group search no longer **starves on a low-ISa merge**: its fixed `I/sigma >= 3` cut could select nothing at all, since a merged sigma is floored so that no reflection reads above `ISa = 1/b`, leaving every operator correlation undefined and the point group at 1. The cut is now capped at the merge's own quantile, and is unchanged on healthy merges. @@ -44,6 +46,9 @@ This is an UNSTABLE release. It includes many experimental features, as well as * Build: `PROJECT_VERSION` is derived from the `VERSION` file rather than hardcoded, and the version is read from it once instead of in two places; and the Bragg integrator no longer uses `M_PI`, which is not standard C++ and left `image_analysis` - which the Windows viewer builds - failing to compile under MSVC. * Docs: Document device-side bitshuffle+LZ4 image decoding in CPU_DATA_ANALYSIS.md; add a Release contents page (what each artefact ships, its CPU baseline, and how the CUDA builds behave without a GPU); correct the RHEL 9 repository URLs and package names; bring the rugnux (`--model`), viewer (processing jobs) and writer pages up to date. +**Breaking change to the rugnux command line:** +* `--azint-only` and `--scale` are **removed**, replaced by `--mode azint` and `--mode scale`; the full pipeline is `--mode mx` and remains the default. A script passing the old flags now fails with the list of valid modes rather than silently running the wrong one. + **Breaking changes to OpenAPI** - regenerate the client (`jfjoch-client` 1.0.0-rc.161, `frontend/src/client`) or read the affected fields as optional: * `image_scale_b` is removed from the `plot_type` enum, so a client requesting that plot now gets an error rather than a curve. * `azim_int_settings.high_q_recipA` and `spot_finding_settings.high_resolution_limit` are no longer `required`. Both mean "as far as the detector reaches" when unset and are omitted from the response instead of carrying a placeholder value, which raises in a client generated from an rc.160-or-earlier spec. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 5365a6ce..befd54f2 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -482,13 +482,19 @@ The refined pass is written as the canonical `_*` output; the pass-1 (he Everything above fits the geometry to *Bragg* data, where the beam centre is the weakest parameter: it is gauge-coupled to the crystal orientation, which is why §7.5 restrains it toward the header value and commits only a sub-1 % move. A **powder ring has no orientation to be coupled to**. Where it falls on the detector depends on the geometry and on nothing else, which makes a calibrant — LaB₆, silver behenate, CeO₂, silicon — or even ice an independent constraint on exactly the quantity Bragg data cannot pin. -The ring positions are generated from the calibrant's cell, matched to the observed rings, and the geometry is refined (Ceres, five parameters: beam centre, distance, and the two detector tilts) so that the $|s|$ predicted at each observed ring point matches the ring it belongs to. +The ring positions are matched to the observed rings and the geometry is refined (Ceres, five parameters: beam centre, distance, and the two detector tilts) so that the $|s|$ predicted at each observed ring point matches the ring it belongs to. + +**Calibrants.** LaB₆, silver behenate, CeO₂ and silicon are held as unit cells and their rings enumerated from them. Ice is held as the eleven **measured** hexagonal-ice ring positions of §3.3 instead, because hexagonal ice is $P6_3/mmc$ and enumerating $hkl$ from its cell would emit rings that are systematically absent. A calibrant is therefore a list of ring $q$ values throughout, not a cell. **What a ring can and cannot determine.** A ring is a conic centred on the beam, so a wrong centre makes its apparent radius oscillate once per turn, $r(\phi)=R+\delta_x\cos\phi+\delta_y\sin\phi$, with the **same amplitude on every ring**. A detector tilt $\beta$ produces a $\cos\phi$ term too — not the $\cos2\phi$ one might expect — but one that grows as the ring's radius *squared*, $r(\phi)=R+(R^2/F)(\beta_x\cos\phi+\beta_y\sin\phi)$; the true $\cos2\phi$ term is $O(R^3\beta^2/F^2)$, hundredths of a pixel. The two are therefore separated by how the amplitude scales with radius, which needs **at least two rings** — on a single ring they are exactly degenerate. None of this uses the calibrant's $d$-spacings, so the centre is determined without assuming anything about the standard. The **distance** is different: it follows from $r=F\tan2\theta$ with $\sin\theta=\lambda/2d$, so a fractional error in the lattice constant passes straight into it, and the $\lambda$–$F$ pair is separated only by the curvature of $\tan(2\arcsin(\lambda/2d))$ across the rings — $\partial\ln r/\partial\ln F=1$ at every ring against $\partial\ln r/\partial\ln\lambda=4\tan\theta/\sin4\theta$, which runs from about 1.05 at low angle to 1.43 at high. That lever collapses as the detector moves back and the rings crowd into small $2\theta$, so distance is a short-distance measurement and the wavelength is better calibrated by other means. -**Reading the rings.** Historically the ring points came from the spot list of a single image, which samples an arc wherever the spot finder's threshold happens to bite. They can instead be read from the **azimuthally-binned profile (§2) summed over a run**: for each ring and each azimuthal sector, the radial peak is fitted against a locally interpolated background and the measured $(q,\phi)$ mapped back through the current geometry to the pixel it came from. Sixteen to thirty-two sectors are enough — beyond that the limit is the ring's own texture, not counting statistics — and the accumulated profile is the same size however many images went into it. +**Reading the rings.** The ring points come from one of two measurements, both accumulated over **every processed image** rather than one. The default reads the **azimuthally-binned profile (§2) summed over the run**: for each ring and each azimuthal sector, the radial peak is fitted against a locally interpolated background and the measured $(q,\phi)$ mapped back through the current geometry to the pixel it came from. The alternative pools the **spot lists**, which samples each arc wherever the spot finder's threshold happens to bite. The accumulated profile is the same size however many images went into it; the pooled spot list is capped, each image contributing an equal share. + +A plain radial profile — one azimuthal sector — has averaged the ring over every direction and carries no centre at all, so the profile route requires at least four sectors and uses 32 by default. Sixteen to thirty-two are enough; beyond that the limit is the ring's own texture, not counting statistics. + +The extraction window around a ring is capped at half the gap to its neighbour, because the background under a peak is taken from the ends of that window: hexagonal ice has three rings within 0.06 Å⁻¹ of one another, which a fixed window merges into a single peak. Where only one ring is in reach the two tilts are held at their input values rather than fitted, since on a single ring they are degenerate with the centre (above) and the fit would otherwise trade the centre away for them. --- diff --git a/docs/RUGNUX.md b/docs/RUGNUX.md index 27672701..955c5168 100644 --- a/docs/RUGNUX.md +++ b/docs/RUGNUX.md @@ -82,14 +82,56 @@ merged in P41212 against a P43212 mo into the model's hand), and — when no reference MTZ has already fixed it — a merohedral indexing ambiguity, by keeping the candidate reindexing with the lowest R-free. -## Re-scaling and re-merging (`rugnux --scale`) +## Re-scaling and re-merging (`rugnux --mode scale`) -The `--scale` mode re-scales and merges the *already-integrated* reflections stored in a +The `scale` mode re-scales and merges the *already-integrated* reflections stored in a `_process.h5` file, without re-running spot finding or integration. Use it to re-merge quickly with a different space group, resolution limit, anomalous setting or reference MTZ. It reuses the same `-o/-N/-s/-e/-S/-A/-B/-z/--scaling-*` options as the full run, and (unlike the full pipeline) does not run a space-group search, so pass `-S` for the correct symmetry. +## Detector calibration from powder rings (`rugnux --mode calibration`) + +The `calibration` mode determines the detector geometry — PONI *x*/*y*, the two tilts +rot1/rot2 and the distance — from the powder rings of a calibrant, and writes it as a +pyFAI **`.poni`** file alongside a printed report of how far each parameter moved from +the header. Bragg data pin the beam centre worst (it is gauge-coupled to the crystal orientation); +a powder ring has no orientation to be coupled to, so this is the measurement that fixes it. + +``` +rugnux --mode calibration --calibrant lab6 -N 8 -o det LaB6_master.h5 +``` + +`--calibrant` takes `lab6`, `agbh` (silver behenate), `ceo2`, `si` or `ice`, case-insensitively. +**`ice`** calibrates a real experiment against its own ice rings — no calibrant exposure needed — +and is the reason a calibrant is a list of ring positions rather than a unit cell: hexagonal ice +is P63/mmc, so rings enumerated from its cell would include systematically absent ones. + +`--calibration` picks how the rings are measured, and **both use the whole dataset** — `-s`/`-e`/`-t` +select which images: + +- **`rings`** (default) sums the (*q* × azimuth) azimuthal profile over every processed image into + one map and fits the ring arcs in it. A powder ring is an arc, not a set of spots, and the summed + profile measures it at every azimuth with all the run's counts behind it. It needs the profile to + be binned in azimuth, so this mode defaults `--azim-phi-bins` to 32. +- **`spots`** pools the found spots of every processed image and fits those. It determines the + centre from scratch (a Hough circle vote, which quantises it to a whole pixel) and then refines. + +Both routes read the ring position out of a binned profile or a spot centroid, so the radial +sampling matters: at a long detector distance the default 0.01 Å⁻¹ *q* bin is several pixels wide +and quantises the `rings` route accordingly — pass a finer `--azim-q-spacing` there (the total +*q* × azimuth bin count must stay under 65534). + +The report prints the fitted geometry, the scatter of the ring points about the fitted rings and the +standard error that implies on the centre. That error is *formal*: it measures the scatter of the +points, not whether the rings themselves are trustworthy, so it stays small when a fit goes wrong +for a structural reason — one visible ring, or ice that is textured rather than smooth. + +Both the PONI (the point of normal incidence, which is what a `.poni` file stores) and the direct +beam (where the beam lands, which is what most other programs call the beam centre) are printed. +They differ by *distance* × tan(rot) once the detector is tilted, which on a 0.3° tilt at 300 mm is +several pixels — enough to look like a disagreement with another program when there is none. + ## Quick start ### Rotation data @@ -192,12 +234,21 @@ General: | `-t, --stride ` | Process every *n*-th image (default: 1) | | `-v, --verbose` | Verbose output | -Modes (default: full analysis — spot finding, indexing, integration and merging): +Mode — `--mode ` (default `mx`): + +| Value | Description | +| --- | --- | +| `mx` | Full analysis — spot finding, indexing, integration and merging | +| `azint` | Only azimuthal integration (no spot finding/indexing); writes `_process.h5` | +| `scale` | Only re-scale/merge the already-integrated reflections in the input `_process.h5` (no re-integration) | +| `calibration` | Determine the detector geometry from powder rings; writes `.poni` | + +Calibration (`--mode calibration`): | Option | Description | | --- | --- | -| `--azint-only` | Only run azimuthal integration (no spot finding/indexing); writes `_process.h5` | -| `--scale` | Only re-scale/merge the already-integrated reflections in the input `_process.h5` (no re-integration) | +| `--calibrant ` | Powder standard: `lab6` \| `agbh` \| `ceo2` \| `si` \| `ice` (default `lab6`, case-insensitive) | +| `--calibration ` | How the rings are measured: `rings` \| `spots` (default `rings`; see above). `rings` defaults `--azim-phi-bins` to 32 | Spot finding: diff --git a/docs/TOOLS.md b/docs/TOOLS.md index d9f74dcf..5e9fa0ea 100644 --- a/docs/TOOLS.md +++ b/docs/TOOLS.md @@ -9,10 +9,11 @@ number of command-line tools. Each prints its own usage when run with `-h` or wi ### rugnux Offline CLI tool that runs the full crystallographic analysis pipeline (spot finding, indexing, integration, scaling/merging) on a stored HDF5 dataset, producing a `_process.h5` file and, when -merging, reflection files. Merging is on by default (`--no-merge` disables it). Two extra modes -narrow the work: `--azint-only` runs only azimuthal integration (no spot finding/indexing), and -`--scale` re-scales/merges the already-integrated reflections in a `_process.h5` without -re-integrating. See [rugnux](RUGNUX.md). +merging, reflection files. Merging is on by default (`--no-merge` disables it). `--mode` picks what a +run does: `mx` (the above, the default), `azint` (only azimuthal integration, no spot +finding/indexing), `scale` (re-scale/merge the already-integrated reflections in a `_process.h5` +without re-integrating) or `calibration` (detector geometry from a calibrant's powder rings, written +as a `.poni` file). See [rugnux](RUGNUX.md). ### jfjoch_extract_hkl Extracts reflections (HKL list) from a Jungfraujoch master file; can sum the same HKL across diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.cpp b/image_analysis/geom_refinement/AssignSpotsToRings.cpp index 01691ddd..07dcdbf1 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.cpp +++ b/image_analysis/geom_refinement/AssignSpotsToRings.cpp @@ -194,7 +194,7 @@ float GuessDetectorDistance(const DiffractionGeometry& geom, float ring_radius_p return det_dist_mm; } -std::vector GuessInitialGeometry(DiffractionGeometry &geom, const std::vector &v, float calibrant_a_A) { +std::vector GuessInitialGeometry(DiffractionGeometry &geom, const std::vector &v, float largest_ring_d_A) { // Reset rotations. The model assumes these are very small in any case! geom.PoniRot1_rad(0.0).PoniRot2_rad(0.0).PoniRot3_rad(0.0); @@ -214,15 +214,14 @@ std::vector GuessInitialGeometry(DiffractionGeometry &geom, const auto cluster_annot = AnalyzeClusters(radii, clusters); - float det_distance = GuessDetectorDistance(geom, cluster_annot[0].R_obs, calibrant_a_A); + float det_distance = GuessDetectorDistance(geom, cluster_annot[0].R_obs, largest_ring_d_A); geom.DetectorDistance_mm(det_distance); return cluster_annot; } -void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant_uc) { - std::vector ring_q = CalculateXtalRings(calibrant_uc); +void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q) { if (ring_q.empty()) - throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Unit cell problem"); + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "No calibrant rings given"); auto cluster_annot = GuessInitialGeometry(geom, v, 2 * PI / ring_q[0]); std::vector optimizer_input; @@ -258,21 +257,24 @@ void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, geom = optimizer.Run(optimizer_input); } -void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant) { - std::vector q_ring = CalculateXtalRings(calibrant); - +std::vector AssignSpotsToRings(const DiffractionGeometry &geom, + const std::vector &v, + const std::vector &ring_q) { std::vector optimizer_input; for (const auto& s: v) { float q_obs = 2 * PI / geom.PxlToRes(s.x, s.y); - for (const auto &q : q_ring) { + for (const auto &q : ring_q) { if (std::fabs(q - q_obs) < 0.1) { optimizer_input.push_back({s.x, s.y, q}); break; } } } + return optimizer_input; +} +void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q) { RingOptimizer optimizer(geom); - geom = optimizer.Run(optimizer_input); + geom = optimizer.Run(AssignSpotsToRings(geom, v, ring_q)); } diff --git a/image_analysis/geom_refinement/AssignSpotsToRings.h b/image_analysis/geom_refinement/AssignSpotsToRings.h index f4abff2e..e2f170da 100644 --- a/image_analysis/geom_refinement/AssignSpotsToRings.h +++ b/image_analysis/geom_refinement/AssignSpotsToRings.h @@ -32,8 +32,18 @@ std::vector CalculateXtalRings(const UnitCell &cell, int hkl_max = 6); std::vector CalculateCubicXtalRings( float a, int hkl_max = 4); float GuessDetectorDistance(const DiffractionGeometry& geom, float ring_radius_pxl, float d_A); -std::vector GuessInitialGeometry(DiffractionGeometry &geom, const std::vector &v, float calibrant_a_A); -void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant); -void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const UnitCell &calibrant); +// The calibrant enters the three functions below as its ring q values (2*pi/d, 1/A, ascending - what +// CalculateXtalRings and CalibrantRings return), not as a UnitCell, so a standard whose rings are +// measured rather than enumerated from a cell - hexagonal ice - can be calibrated against as well. +// GuessInitialGeometry needs only the outermost d-spacing, to seed the distance from the innermost ring. +std::vector GuessInitialGeometry(DiffractionGeometry &geom, const std::vector &v, float largest_ring_d_A); +void GuessGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q); +void OptimizeGeometry(DiffractionGeometry &geom, const std::vector &v, const std::vector &ring_q); + +// Each spot paired with the calibrant ring nearest its observed q, as the points RingOptimizer fits. +// Spots more than 0.1 1/A from every ring are dropped rather than forced onto the closest one. +std::vector AssignSpotsToRings(const DiffractionGeometry &geom, + const std::vector &v, + const std::vector &ring_q); diff --git a/image_analysis/geom_refinement/CMakeLists.txt b/image_analysis/geom_refinement/CMakeLists.txt index 164a96d4..c364dad5 100644 --- a/image_analysis/geom_refinement/CMakeLists.txt +++ b/image_analysis/geom_refinement/CMakeLists.txt @@ -4,6 +4,8 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC RingOptimizer.h AssignSpotsToRings.cpp AssignSpotsToRings.h + Calibrants.cpp + Calibrants.h RingsFromProfile.cpp RingsFromProfile.h XtalOptimizer.cpp diff --git a/image_analysis/geom_refinement/Calibrants.cpp b/image_analysis/geom_refinement/Calibrants.cpp new file mode 100644 index 00000000..e3140200 --- /dev/null +++ b/image_analysis/geom_refinement/Calibrants.cpp @@ -0,0 +1,61 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include "Calibrants.h" +#include "AssignSpotsToRings.h" // CalculateXtalRings +#include "../../common/Definitions.h" +#include "../../common/JFJochMath.h" + +namespace { + std::string lower_case(const std::string &s) { + std::string out = s; + std::transform(out.begin(), out.end(), out.begin(), + [](unsigned char c) { return static_cast(std::tolower(c)); }); + return out; + } +} + +const std::vector &Calibrants() { + static const std::vector table = { + {"LaB6", UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90, 90, 90)}, + // Silver behenate. T. C. Huang, H. Toraya, T. N. Blanton, Y. Wu, J. Appl. Cryst. 26 (1993), 180-184. + {"AgBh", UnitCell(5.1769, 4.7218, 58.380, 89.440, 89.634, 75.854)}, + // CeO2, cubic fluorite, NIST SRM 674b. + {"CeO2", UnitCell(5.4115, 5.4115, 5.4115, 90, 90, 90)}, + // Silicon, cubic, NIST SRM 640. + {"Si", UnitCell(5.43102, 5.43102, 5.43102, 90, 90, 90)}, + // Hexagonal ice: measured ring positions, not a cell (see the header). + {"ice", std::nullopt} + }; + return table; +} + +std::vector CalibrantRings(const std::string &name) { + const std::string key = lower_case(name); + for (const auto &c : Calibrants()) { + if (lower_case(c.name) != key) + continue; + if (c.cell) + return CalculateXtalRings(*c.cell); + + std::vector q; + for (const float d : ICE_RING_RES_A) + q.push_back(static_cast(2.0 * PI) / d); + std::sort(q.begin(), q.end()); + return q; + } + return {}; +} + +std::string CalibrantNameList() { + std::string out; + for (const auto &c : Calibrants()) { + if (!out.empty()) + out += ", "; + out += lower_case(c.name); + } + return out; +} diff --git a/image_analysis/geom_refinement/Calibrants.h b/image_analysis/geom_refinement/Calibrants.h new file mode 100644 index 00000000..faf04add --- /dev/null +++ b/image_analysis/geom_refinement/Calibrants.h @@ -0,0 +1,30 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "../../common/UnitCell.h" + +// The powder standards the detector geometry can be calibrated against, shared by +// rugnux --mode=calibration and the viewer's powder-calibration panel. +struct Calibrant { + std::string name; // as offered on the command line and in the GUI (matched case-insensitively) + std::optional cell; // unset for ice - see CalibrantRings() +}; + +const std::vector &Calibrants(); + +// The rings a calibrant produces, as q = 2*pi/d in 1/A - the unit CalculateXtalRings returns and +// RingOptimizerInput::q_expected expects - sorted ascending. Empty if the name is not a calibrant. +// +// A ring LIST rather than a UnitCell because ice cannot be given as one: hexagonal ice is P6_3/mmc, so +// enumerating hkl from its cell would emit rings that are systematically absent, while the eleven +// entries of ICE_RING_RES_A are measured ring positions. +std::vector CalibrantRings(const std::string &name); + +// "lab6, agbh, ceo2, si, ice" - for the error message when a name is not one of them. +std::string CalibrantNameList(); diff --git a/image_analysis/geom_refinement/RingOptimizer.cpp b/image_analysis/geom_refinement/RingOptimizer.cpp index 1a120180..1e845d24 100644 --- a/image_analysis/geom_refinement/RingOptimizer.cpp +++ b/image_analysis/geom_refinement/RingOptimizer.cpp @@ -1,6 +1,8 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include + #include "../../common/JFJochMath.h" #include "RingOptimizer.h" #include "ceres/ceres.h" @@ -82,6 +84,19 @@ DiffractionGeometry RingOptimizer::Run(const std::vector &in ); } + // A single ring cannot tell the beam centre from the detector tilt: both displace its radius as + // cos(phi), and what separates them is only how that amplitude scales with the ring's radius, which + // takes two rings. Hold the tilt where it was given, so the one thing a single ring does fix - where + // its centre lies - comes out rather than being traded away against an unconstrained tilt. + const bool single_ring = !input.empty() + && std::all_of(input.begin(), input.end(), [&](const RingOptimizerInput &p) { + return p.q_expected == input.front().q_expected; + }); + if (single_ring) { + problem.SetParameterBlockConstant(&rot1); + problem.SetParameterBlockConstant(&rot2); + } + // Configure solver ceres::Solver::Options options; options.linear_solver_type = ceres::DENSE_QR; diff --git a/image_analysis/geom_refinement/RingsFromProfile.cpp b/image_analysis/geom_refinement/RingsFromProfile.cpp index c9372845..b33f26d1 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.cpp +++ b/image_analysis/geom_refinement/RingsFromProfile.cpp @@ -5,7 +5,6 @@ #include #include "RingsFromProfile.h" -#include "AssignSpotsToRings.h" // CalculateXtalRings #include "../../common/JFJochMath.h" namespace { @@ -75,7 +74,7 @@ float SectorPeakQ(const std::vector &profile, int32_t q_bins, int phi_bin std::vector RingsFromAzimuthalProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, - const UnitCell &calibrant, + const std::vector &calibrant_ring_q, float q_window_recipA, float min_peak_over_noise) { std::vector out; @@ -93,14 +92,27 @@ std::vector RingsFromAzimuthalProfile(const std::vector(q_bins) * q_spacing; - const int window_bins = std::max(3, static_cast(std::lround(q_window_recipA / q_spacing))); - for (const float q_ring : CalculateXtalRings(calibrant)) { - if (!(q_ring - q_window_recipA > low_q) || !(q_ring + q_window_recipA < high_q)) + for (size_t i = 0; i < calibrant_ring_q.size(); ++i) { + const float q_ring = calibrant_ring_q[i]; + // 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. + float window = q_window_recipA; + if (i > 0) + window = std::min(window, 0.5f * (q_ring - calibrant_ring_q[i - 1])); + if (i + 1 < calibrant_ring_q.size()) + window = std::min(window, 0.5f * (calibrant_ring_q[i + 1] - 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; diff --git a/image_analysis/geom_refinement/RingsFromProfile.h b/image_analysis/geom_refinement/RingsFromProfile.h index 5fc4e863..c2e7b4d8 100644 --- a/image_analysis/geom_refinement/RingsFromProfile.h +++ b/image_analysis/geom_refinement/RingsFromProfile.h @@ -7,7 +7,6 @@ #include "../../common/AzimuthalIntegrationMapping.h" #include "../../common/DiffractionGeometry.h" -#include "../../common/UnitCell.h" #include "RingOptimizer.h" // Turn an accumulated (q x azimuth) powder profile into ring points for RingOptimizer. @@ -35,9 +34,12 @@ // measured (q, phi) back into the pixel it came from - RingOptimizer then refines that geometry so the // 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. +// +// 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. std::vector RingsFromAzimuthalProfile(const std::vector &profile, const AzimuthalIntegrationMapping &mapping, const DiffractionGeometry &geom, - const UnitCell &calibrant, + const std::vector &calibrant_ring_q, float q_window_recipA = 0.06f, float min_peak_over_noise = 3.0f); diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index a618af95..d9328a29 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -6,6 +6,8 @@ ADD_LIBRARY(Rugnux STATIC RugnuxDefaults.cpp RugnuxDefaults.h Rugnux.h + RugnuxCalibration.cpp + RugnuxCalibration.h RugnuxCommandLine.cpp RugnuxCommandLine.h ModelValidation.cpp @@ -14,8 +16,8 @@ ADD_LIBRARY(Rugnux STATIC TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochWriter gemmi) -# rugnux is the single offline analysis CLI: full pipeline by default, plus --azint-only -# (azimuthal integration) and --scale (re-scale/merge stored reflections) modes. +# rugnux is the single offline analysis CLI; --mode picks what it does: mx (the full pipeline, the +# default), azint, scale (re-scale/merge stored reflections) or calibration (powder-ring geometry). ADD_EXECUTABLE(rugnux rugnux_cli.cpp) TARGET_LINK_LIBRARIES(rugnux Rugnux JFJochReader JFJochImageAnalysis JFJochWriter) INSTALL(TARGETS rugnux RUNTIME COMPONENT viewer) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index d5c90599..de7116a0 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -455,6 +455,12 @@ 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. + 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 write_files = write_output && !config_.output_prefix.empty(); // Output/runtime invariants. Algorithm settings (indexing, scaling, integration, polarization, @@ -526,16 +532,17 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b writer = std::make_unique(start_message, /*check_overwrite_at_start=*/true, /*trusted_path=*/true); + const char *mode_name = full ? "full analysis" + : calibration ? "powder calibration" : "azimuthal integration"; logger.Info("Processing {} images (range {}-{}, stride {}) using {} threads [{}]", - images_to_process, start_image, end_image, config_.stride, config_.nthreads, - full ? "full analysis" : "azimuthal integration"); + images_to_process, start_image, end_image, config_.stride, config_.nthreads, mode_name); if (observer) - observer->OnPhase(full ? "Full analysis" : "Azimuthal integration"); + observer->OnPhase(full ? "Full analysis" : calibration ? "Powder calibration" : "Azimuthal integration"); // Full-analysis shared engines. std::unique_ptr indexer_pool; std::unique_ptr indexer; - if (full) { + if (per_image_analysis) { auto indexing_settings = experiment_.GetIndexingSettings(); // Two-pass rotation fires the FFT indexer only twice - the two first-pass schemes, run in // parallel (the main per-image loop then uses the forced lattice, no FFT). Size the pool to @@ -845,7 +852,7 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // give the same remedies. constexpr const char *next_steps = " Nothing was integrated. Check the beam centre (--beam-x / --beam-y), raise " - "--max-spots if the frames are well populated with spots, and use --azint-only to " + "--max-spots if the frames are well populated with spots, and use --mode azint to " "write the per-image analysis for inspection."; if (!best.result.has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, @@ -928,6 +935,15 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b }; std::atomic total_uncompressed_bytes = 0; + // Calibration by spots: the pooled spot list of the whole run. Ring clustering is O(n^2) in it and + // the Hough circle centre O(n^3) in its first few hundred, so each image contributes a fair share of + // a fixed budget instead of everything it found - a powder ring is over-determined either way. + constexpr size_t CALIBRATION_MAX_SPOTS = 20000; + const size_t calibration_spots_per_image = + std::max(1, CALIBRATION_MAX_SPOTS / static_cast(images_to_process)); + std::mutex calibration_spot_m; + std::vector calibration_spot_list; + auto azint_worker = [&]() { std::vector decompression_buffer; ImagePreprocessorCPU preprocessor(experiment_, pixel_mask_); @@ -1025,6 +1041,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b msg.run_number = experiment_.GetRunNumber(); msg.run_name = experiment_.GetRunName(); + if (calibration_spots) { + 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]); + } + plots.Add(msg, profile); if (writer) writer->Write(msg); note_written(ordinal); @@ -1037,8 +1059,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (observer) observer->OnPhase("Processing images"); - std::function worker = full ? std::function(full_worker) - : std::function(azint_worker); + std::function worker = per_image_analysis ? std::function(full_worker) + : std::function(azint_worker); std::vector > futures; futures.reserve(config_.nthreads); const auto image_loop_start = std::chrono::steady_clock::now(); @@ -1064,6 +1086,28 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b result.mean_processing_time = plots.GetMeanProcessingTime(); result.indexing_rate = plots.GetIndexingRate(); + // Powder calibration: the whole run has now been accumulated - the summed (q x azimuth) profile in + // plots, or the pooled spot list - so fit the detector geometry to it. The caller reports it and + // writes the PONI file; nothing here is stored in the _process.h5. + if (calibration && !cancelled_) { + if (observer) + observer->OnPhase("Powder calibration"); + if (calibration_spots) { + logger.Info("Powder calibration from {} pooled spots", calibration_spot_list.size()); + result.calibration = CalibrateFromSpots(calibration_spot_list, + experiment_.GetDiffractionGeometry(), + config_.calibrant_ring_q); + } else { + const auto *summed = plots.GetAzIntProfileObject(); + if (!summed) + throw JFJochException(JFJochExceptionCategory::CalibrationError, + "No azimuthal profile was accumulated"); + result.calibration = CalibrateFromProfile(summed->GetResult(), mapping, + experiment_.GetDiffractionGeometry(), + config_.calibrant_ring_q); + } + } + // End message (also written to the file). EndMessage end_msg; // The EXTENT of the per-image datasets (see note_written), not the success count: those two diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 7da26936..91aed5dc 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -20,18 +20,20 @@ #include "../image_analysis/scale_merge/TwinningAnalysis.h" // TwinningAnalysisResult #include "../image_analysis/scale_merge/SearchSpaceGroup.h" // SearchSpaceGroupResult #include "../image_analysis/rotation_indexer/RotationIndexer.h" // RotationIndexerResult +#include "RugnuxCalibration.h" // CalibrationMethod, CalibrationResult class JFJochHDF5Reader; -// Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by rugnux (including its -// --azint-only mode) and (later) the viewer. The full processing workflow lives here, not in the CLIs: +// Offline reprocessing of a stored Jungfraujoch HDF5 dataset, shared by rugnux (all of its --mode +// values that read images) and (later) the viewer. The full workflow lives here, not in the CLIs: // setup, an optional two-pass rotation-indexing pre-pass, a parallel per-image loop (std::thread), // an optional scaling/merging post-pass, and the _process.h5 output. The detector geometry and all // algorithm settings are configured on the DiffractionExperiment by the caller; ProcessConfig only // carries run control. Cancellable from any thread (e.g. SIGINT or a GUI button) via Cancel(). enum class ProcessMode { - AzimuthalIntegration, // preprocess + azimuthal integration only (rugnux --azint-only) - FullAnalysis // spot finding + indexing + refinement + integration (rugnux) + AzimuthalIntegration, // preprocess + azimuthal integration only (rugnux --mode=azint) + FullAnalysis, // spot finding + indexing + refinement + integration (rugnux --mode=mx) + Calibration // detector geometry from powder rings (rugnux --mode=calibration) }; struct ProcessConfig { @@ -54,7 +56,12 @@ struct ProcessConfig { // only wants the ISa/R-meas numbers). Only relevant to FullAnalysis with scaling. bool write_merged = true; - SpotFindingSettings spot_finding; // FullAnalysis spot finding + SpotFindingSettings spot_finding; // FullAnalysis spot finding (and Calibration by spots) + + // Calibration: how the powder rings are measured, and the calibrant's rings as q = 2*pi/d (1/A). + // The library takes the ring list, not a calibrant name - the CLI resolves it (CalibrantRings). + CalibrationMethod calibration_method = CalibrationMethod::Rings; + std::vector calibrant_ring_q; // Rotation indexing (FullAnalysis) bool rotation_indexing = false; @@ -136,6 +143,9 @@ struct ProcessResult { // search when the user did not fix one, otherwise the fixed group. std::optional space_group_number; + // Calibration mode: the geometry fitted to the powder rings of the whole run. + std::optional calibration; + // Structured result of the de-novo space-group search (point group + ranked candidate scores), // populated only when a search ran (no space group was fixed). The library no longer renders it to // text - the CLI formats it for stdout, the viewer draws it as a table. diff --git a/rugnux/RugnuxCalibration.cpp b/rugnux/RugnuxCalibration.cpp new file mode 100644 index 00000000..a14aa517 --- /dev/null +++ b/rugnux/RugnuxCalibration.cpp @@ -0,0 +1,108 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include + +#include + +#include "RugnuxCalibration.h" +#include "../common/GitInfo.h" +#include "../image_analysis/geom_refinement/AssignSpotsToRings.h" +#include "../image_analysis/geom_refinement/RingOptimizer.h" +#include "../image_analysis/geom_refinement/RingsFromProfile.h" + +namespace { + +// How well the ring points sit on the fitted rings, in a unit a user can judge: the radial distance in +// pixels between where a point is and where the fitted geometry puts its ring. The fit's own residual +// is in q, so it is divided by the local dq/dr - measured by stepping one pixel outward along the radius +// rather than assumed, since dq/dr varies with two-theta and with the tilt. +// +// The beam centre enters a ring's radius as r(phi) = R + dx cos(phi) + dy sin(phi), so fitting it to n +// points of scatter s leaves the textbook var = 2 s^2 / n on each of dx and dy. That is the number that +// separates a beam centre that was measured from one that was merely reported. +CalibrationResult Summarize(const DiffractionGeometry &fitted, + const std::vector &points) { + CalibrationResult result; + result.geometry = fitted; + + const float cx = fitted.GetBeamX_pxl(); + const float cy = fitted.GetBeamY_pxl(); + double sum_sq = 0.0; + for (const auto &p : points) { + const float r = std::hypot(p.x - cx, p.y - cy); + if (!(r > 1.0f)) + continue; + const float q = fitted.PxlToQ(p.x, p.y); + const float dq_dr = fitted.PxlToQ(p.x + (p.x - cx) / r, p.y + (p.y - cy) / r) - q; + if (!(std::abs(dq_dr) > 0.0f)) + continue; + const double dr = (q - p.q_expected) / dq_dr; + sum_sq += dr * dr; + ++result.ring_points; + } + if (result.ring_points > 0) { + result.rms_radial_pxl = std::sqrt(sum_sq / static_cast(result.ring_points)); + result.beam_sigma_pxl = result.rms_radial_pxl + * std::sqrt(2.0 / static_cast(result.ring_points)); + } + return result; +} + +} // namespace + +CalibrationResult CalibrateFromProfile(const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + const std::vector &calibrant_ring_q) { + const auto points = RingsFromAzimuthalProfile(profile, mapping, geom, calibrant_ring_q); + if (points.empty()) + throw JFJochException(JFJochExceptionCategory::CalibrationError, + "No powder ring found in the summed azimuthal profile"); + return Summarize(RingOptimizer(geom).Run(points), points); +} + +CalibrationResult CalibrateFromSpots(const std::vector &spots, + const DiffractionGeometry &geom, + const std::vector &calibrant_ring_q) { + DiffractionGeometry fitted = geom; + // From scratch (Hough circle centre + ring clustering), then refined: the guess pins the centre to a + // whole pixel and only sees the spots its clustering kept, so the refine re-matches every spot at + // that geometry. + GuessGeometry(fitted, spots, calibrant_ring_q); + OptimizeGeometry(fitted, spots, calibrant_ring_q); + return Summarize(fitted, AssignSpotsToRings(fitted, spots, calibrant_ring_q)); +} + +void WritePoniFile(const std::string &path, const DiffractionExperiment &experiment, + const DiffractionGeometry &geom) { + std::ofstream f(path); + if (!f) + throw JFJochException(JFJochExceptionCategory::FileWriteError, "Cannot write " + path); + + const double pixel_m = geom.GetPixelSize_mm() * 1e-3; + // pyFAI's axis convention is the trap: Poni1 (and pixel1) is the SLOW axis - rows, our y - and + // Poni2 the FAST axis - columns, our x - both in metres from the detector origin. A transposed PONI + // file is silently wrong, so the mapping is spelled out here rather than left to the reader. + // + // DiffractionGeometry's beam_x/beam_y IS the PONI: LabCoord rotates the vector measured FROM that + // pixel, i.e. it is the point of normal incidence, so it maps straight across with no correction. + // GetDirectBeam_pxl() is a different quantity - where the direct beam lands - and parts from the + // PONI as soon as rot1/rot2 are non-zero. + f << fmt::format("# Calibration done by Jungfraujoch rugnux {}\n", jfjoch_version()); + f << "poni_version: 2\n"; + f << "Detector: Detector\n"; + f << fmt::format("Detector_config: {{\"pixel1\": {:g}, \"pixel2\": {:g}, \"max_shape\": [{}, {}]}}\n", + pixel_m, pixel_m, experiment.GetYPixelsNumConv(), experiment.GetXPixelsNumConv()); + f << fmt::format("Distance: {:.9g}\n", geom.GetDetectorDistance_mm() * 1e-3); + f << fmt::format("Poni1: {:.9g}\n", geom.GetBeamY_pxl() * pixel_m); + f << fmt::format("Poni2: {:.9g}\n", geom.GetBeamX_pxl() * pixel_m); + f << fmt::format("Rot1: {:.9g}\n", geom.GetPoniRot1_rad()); + f << fmt::format("Rot2: {:.9g}\n", geom.GetPoniRot2_rad()); + f << fmt::format("Rot3: {:.9g}\n", geom.GetPoniRot3_rad()); + f << fmt::format("Wavelength: {:.9g}\n", geom.GetWavelength_A() * 1e-10); + f.flush(); + if (!f) + throw JFJochException(JFJochExceptionCategory::FileWriteError, "Error writing " + path); +} diff --git a/rugnux/RugnuxCalibration.h b/rugnux/RugnuxCalibration.h new file mode 100644 index 00000000..3f98d664 --- /dev/null +++ b/rugnux/RugnuxCalibration.h @@ -0,0 +1,45 @@ +// 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/DiffractionExperiment.h" +#include "../common/DiffractionGeometry.h" +#include "../common/SpotToSave.h" + +// How the powder rings the detector geometry is fitted to are measured (rugnux --calibration). +enum class CalibrationMethod { + Rings, // the run-summed (q x azimuth) azimuthal profile: the ring measured at every azimuth + Spots // the pooled per-image spot lists: the ring sampled wherever the spot finder bit +}; + +struct CalibrationResult { + DiffractionGeometry geometry; // the fitted geometry + size_t ring_points = 0; // ring measurements the fit used + // Scatter of those measurements about the fitted rings, and the standard error it implies on the + // beam centre. Both in pixels - a calibration that has gone wrong (textured ice, one visible ring) + // says so here, and that is the only warning a user gets. + double rms_radial_pxl = 0.0; + double beam_sigma_pxl = 0.0; +}; + +// Fit the geometry to the rings found in a run-summed azimuthal profile (CalibrationMethod::Rings). +CalibrationResult CalibrateFromProfile(const std::vector &profile, + const AzimuthalIntegrationMapping &mapping, + const DiffractionGeometry &geom, + const std::vector &calibrant_ring_q); + +// 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. +CalibrationResult CalibrateFromSpots(const std::vector &spots, + const DiffractionGeometry &geom, + const std::vector &calibrant_ring_q); + +// Write the geometry as a pyFAI PONI file, the interchange format every azimuthal-integration tool +// reads. Throws if the file cannot be written. +void WritePoniFile(const std::string &path, const DiffractionExperiment &experiment, + const DiffractionGeometry &geom); diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index 16dec5ad..fbeabf84 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -55,8 +55,10 @@ std::string RugnuxCommandLine(const ProcessConfig &config, std::vector args; const bool azint = (config.mode == ProcessMode::AzimuthalIntegration); args.emplace_back("rugnux"); - if (azint) - args.emplace_back("--azint-only"); + if (azint) { + args.emplace_back("--mode"); + args.emplace_back("azint"); + } auto add = [&](const std::string &flag, const std::string &val) { args.push_back(flag); diff --git a/rugnux/RugnuxCommandLine.h b/rugnux/RugnuxCommandLine.h index b8362f70..203204d7 100644 --- a/rugnux/RugnuxCommandLine.h +++ b/rugnux/RugnuxCommandLine.h @@ -9,7 +9,7 @@ class DiffractionExperiment; -// Reconstruct an equivalent rugnux command line (including --azint-only) for a configured run, so a +// Reconstruct an equivalent rugnux command line (including --mode azint) for a configured run, so a // job set up in the GUI can be handed off to a cluster. Covers the settings that matter for the // run, not every obscure flag; geometry is taken from the input file, so geometry overrides are // not emitted. diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ad878983..61ddfc1b 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -22,6 +22,8 @@ #include "../common/DiffractionExperiment.h" #include "../common/PixelMask.h" #include "../common/print_license.h" +#include "../common/JFJochMath.h" +#include "../image_analysis/geom_refinement/Calibrants.h" #include "../image_analysis/LoadFCalcFromMtz.h" #include "../image_analysis/UpdateReflectionResolution.h" #include "../image_analysis/WriteReflections.h" @@ -51,6 +53,18 @@ constexpr double SMOOTH_G_DEFAULT_DEG = 5.0; // Default rot3d per-batch relative-B batch width (bare --relative-b), in degrees of rotation. constexpr double RELATIVE_B_DEFAULT_DEG = 10.0; +// What the run does (--mode). Each of these was once a flag of its own; one option keeps them +// mutually exclusive by construction. +enum class RugnuxMode { MX, Azint, Scale, Calibration }; + +// Azimuthal sectors used by --mode calibration --calibration rings when --azim-phi-bins was not given. +// The default of 1 is a plain radial profile, in which a ring has been averaged over every direction so +// nothing is left to say where its centre is (RingsFromAzimuthalProfile refuses below 4). Sweeping the +// count on a LaB6 exposure and on a crystal's ice rings, 32 and above agree to about 0.1 px and the +// residual is flat, while 8 sectors is visibly coarser (0.3 px away, twice the residual). Going higher +// buys nothing and runs into the 65534-bin limit on q x azimuth, so 32 is the default. +constexpr int CALIBRATION_AZIM_BINS_DEFAULT = 32; + void print_usage() { std::cout << "Usage rugnux {} " << std::endl; std::cout << "Options:" << std::endl; @@ -62,9 +76,17 @@ void print_usage() { std::cout << " -v, --verbose Verbose output" << std::endl; std::cout << std::endl; - std::cout << " Modes (default: full analysis - spot finding, indexing, integration and merging)" << std::endl; - std::cout << " --azint-only Only run azimuthal integration (no spot finding/indexing); writes _process.h5" << std::endl; - std::cout << " --scale Only re-scale/merge the already-integrated reflections in (no re-integration)" << std::endl; + std::cout << " Mode" << std::endl; + std::cout << " --mode What this run does (default: mx)" << std::endl; + std::cout << " mx Full analysis - spot finding, indexing, integration and merging" << std::endl; + std::cout << " azint Only azimuthal integration (no spot finding/indexing); writes _process.h5" << std::endl; + std::cout << " scale Only re-scale/merge the already-integrated reflections in (no re-integration)" << std::endl; + std::cout << " calibration Determine the detector geometry from powder rings; writes .poni" << std::endl; + std::cout << std::endl; + + std::cout << " Calibration (--mode calibration)" << std::endl; + std::cout << " --calibrant Powder standard: lab6|agbh|ceo2|si|ice (default: lab6, case-insensitive)" << std::endl; + std::cout << " --calibration How the rings are measured: rings|spots (default: rings). rings sums the (q x azimuth) azimuthal profile over every processed image and fits the ring arcs in it; spots pools the found spots and fits those. -s/-e/--stride select the images; rings defaults --azim-phi-bins to 32" << std::endl; std::cout << std::endl; std::cout << " Spot finding" << std::endl; @@ -207,8 +229,9 @@ enum { OPT_AZIM_MAX_Q, OPT_AZIM_PHI_BINS, OPT_AZIM_SIGMA_CLIP, - OPT_AZINT_ONLY, - OPT_SCALE, + OPT_MODE, + OPT_CALIBRANT, + OPT_CALIBRATION, OPT_NO_MERGE, OPT_POLARIZATION_CORRECTION, OPT_SOLID_ANGLE_CORRECTION, @@ -236,8 +259,9 @@ static option long_options[] = { {"dump-observations", required_argument, nullptr, OPT_DUMP_OBSERVATIONS}, {"space-group", required_argument, nullptr, 'S'}, {"anomalous", no_argument, nullptr, 'A'}, - {"azint-only", no_argument, nullptr, OPT_AZINT_ONLY}, - {"scale", no_argument, nullptr, OPT_SCALE}, + {"mode", required_argument, nullptr, OPT_MODE}, + {"calibrant", required_argument, nullptr, OPT_CALIBRANT}, + {"calibration", required_argument, nullptr, OPT_CALIBRATION}, {"no-merge", no_argument, nullptr, OPT_NO_MERGE}, {"scale-fulls", no_argument, nullptr, OPT_SCALE_FULLS}, {"no-scale-fulls", no_argument, nullptr, OPT_NO_SCALE_FULLS}, @@ -518,8 +542,9 @@ static int RunRugnux(int argc, char **argv) { int image_stride = 1; bool verbose = false; - bool azint_only = false; // --azint-only: azimuthal integration only (no spot finding/indexing) - bool scale_only = false; // --scale: re-scale/merge stored reflections only (no re-integration) + RugnuxMode mode = RugnuxMode::MX; // --mode + std::string calibrant = "lab6"; // --calibrant + CalibrationMethod calibration_method = CalibrationMethod::Rings; // --calibration bool rotation_indexing = false; bool force_still = false; // --force-still: process a rotation dataset as stills (indexing + scaling) bool two_pass_rotation = true; @@ -819,11 +844,32 @@ static int RunRugnux(int argc, char **argv) { case OPT_MAX_SPOTS: max_spot_count_override = parse_number_arg(optarg, "--max-spots", logger, 1); break; - case OPT_AZINT_ONLY: - azint_only = true; + case OPT_MODE: { + const std::string m = optarg ? optarg : ""; + if (m == "mx") mode = RugnuxMode::MX; + else if (m == "azint") mode = RugnuxMode::Azint; + else if (m == "scale") mode = RugnuxMode::Scale; + else if (m == "calibration") mode = RugnuxMode::Calibration; + else { + logger.Error("Invalid --mode: {} (expected mx|azint|scale|calibration)", m); + return 1; + } break; - case OPT_SCALE: - scale_only = true; + } + case OPT_CALIBRANT: + calibrant = optarg ? optarg : ""; + if (CalibrantRings(calibrant).empty()) { + logger.Error("Unknown --calibrant: {} (expected {})", calibrant, CalibrantNameList()); + return 1; + } + break; + case OPT_CALIBRATION: + if (strcmp(optarg, "rings") == 0) calibration_method = CalibrationMethod::Rings; + else if (strcmp(optarg, "spots") == 0) calibration_method = CalibrationMethod::Spots; + else { + logger.Error("Invalid --calibration: {} (expected rings|spots)", optarg); + return 1; + } break; case OPT_NO_MERGE: run_scaling = false; @@ -1052,11 +1098,6 @@ static int RunRugnux(int argc, char **argv) { nthreads = hw > 0 ? static_cast(hw) : 1; } - if (azint_only && scale_only) { - logger.Error("--azint-only and --scale are mutually exclusive"); - exit(EXIT_FAILURE); - } - // Validate space group number early const gemmi::SpaceGroup *space_group = nullptr; if (space_group_number.has_value()) { @@ -1122,7 +1163,7 @@ static int RunRugnux(int argc, char **argv) { // only takes effect on the full-analysis path, which otherwise determines its cell later by // indexing (unknown here, so nothing can be checked yet). const std::optional data_cell = - scale_only ? dataset->experiment.GetUnitCell() : fixed_reference_unit_cell; + mode == RugnuxMode::Scale ? dataset->experiment.GetUnitCell() : fixed_reference_unit_cell; const auto warning = ReferenceConsistencyWarning( reference, data_cell, space_group_number.has_value() ? std::optional(static_cast(*space_group_number)) @@ -1148,9 +1189,9 @@ static int RunRugnux(int argc, char **argv) { } } - // --scale: re-scale and merge the already-integrated reflections stored in the input file, + // --mode scale: re-scale and merge the already-integrated reflections stored in the input file, // without re-running spot finding or integration (folded in from the former rugnux_scale tool). - if (scale_only) { + if (mode == RugnuxMode::Scale) { const auto total_images = static_cast(reader.GetNumberOfImages()); const int last_image = (end_image < 0 || end_image >= total_images) ? total_images - 1 : end_image; auto reflections = reader.ReadReflections(start_image, last_image); @@ -1415,8 +1456,14 @@ static int RunRugnux(int argc, char **argv) { if (rot2_rad) experiment.PoniRot2_rad(rot2_rad.value()); // --polarization is applied after configure_offline_output below, which sets the rugnux default. + // Calibrating from the run-summed profile needs the profile to be binned in azimuth - see + // CALIBRATION_AZIM_BINS_DEFAULT. An explicit --azim-phi-bins always wins. + if (mode == RugnuxMode::Calibration && calibration_method == CalibrationMethod::Rings + && !azimuthal_bins.has_value()) + azimuthal_bins = CALIBRATION_AZIM_BINS_DEFAULT; + // Azimuthal integration (default q-spacing 0.01 1/A, from AzimuthalIntegrationSettings): the profile - // resolves the narrow ice rings for the ice-ring score. Shared by --azint-only and full analysis. + // resolves the narrow ice rings for the ice-ring score. Shared by every mode that reads images. // -q / --azim-* / correction flags override; defaults come from the input file. { AzimuthalIntegrationSettings azint_settings = experiment.GetAzimuthalIntegrationSettings(); @@ -1440,9 +1487,9 @@ static int RunRugnux(int argc, char **argv) { azint_settings.GetAzimuthalBinCount()); } - // --azint-only: azimuthal integration only (no spot finding / indexing / scaling). Rugnux reads + // --mode azint: azimuthal integration only (no spot finding / indexing / scaling). Rugnux reads // the geometry and azimuthal-integration settings configured above off the experiment. - if (azint_only) { + if (mode == RugnuxMode::Azint) { ProcessConfig config; config.mode = ProcessMode::AzimuthalIntegration; config.start_image = start_image; @@ -1472,6 +1519,99 @@ static int RunRugnux(int argc, char **argv) { return 0; } + // --mode calibration: fit the detector geometry to the calibrant's powder rings, measured over + // every selected image, and write it as a PONI file. + if (mode == RugnuxMode::Calibration) { + if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value()); + + ProcessConfig config; + config.mode = ProcessMode::Calibration; + config.calibration_method = calibration_method; + config.calibrant_ring_q = CalibrantRings(calibrant); + config.start_image = start_image; + config.end_image = end_image; + config.stride = image_stride; + config.nthreads = nthreads; + config.output_prefix = output_prefix; + config.write_process_h5 = false; // the .poni below is the output of this mode + + // Spot finding for --calibration spots. Indexing is off: a calibration wants the spot positions + // and nothing else, and the calibrant is a powder with no lattice to index. + config.spot_finding.enable = true; + config.spot_finding.indexing = false; + config.spot_finding.signal_to_noise_threshold = sigma_spot_finding; + config.spot_finding.photon_count_threshold = photon_count_threshold_spot_finding; + config.spot_finding.min_pix_per_spot = min_pix_per_spot.value_or(2); + config.spot_finding.adaptive_threshold = adaptive_spots.value_or(true); + config.spot_finding.high_resolution_limit = d_min_spot_finding; + config.spot_finding.false_pixels_per_frame = false_pixels_per_frame; + if (d_max_spot_finding > 0.0f) + config.spot_finding.low_resolution_limit = d_max_spot_finding; + experiment.MaxSpotCount(max_spot_count_override.value_or(RUGNUX_MAX_SPOT_COUNT)); + + logger.Info("Powder calibration: calibrant {} ({} rings), method {}, {} azimuthal bins", + calibrant, config.calibrant_ring_q.size(), + calibration_method == CalibrationMethod::Rings ? "rings" : "spots", + experiment.GetAzimuthalIntegrationSettings().GetAzimuthalBinCount()); + + Rugnux process(reader, experiment, *dataset->pixel_mask, config); + g_active_process = &process; + std::signal(SIGINT, handle_sigint); + + ProcessResult result; + try { + result = process.Run(); + } catch (const std::exception &e) { + logger.Error("Calibration failed: {}", e.what()); + return 1; + } + g_active_process = nullptr; + + if (!result.calibration.has_value()) { + logger.Error("Calibration produced no geometry"); + return 1; + } + const auto &cal = *result.calibration; + const auto &g = cal.geometry; + const DiffractionGeometry header = experiment.GetDiffractionGeometry(); + const float pxl_mm = g.GetPixelSize_mm(); + constexpr double RAD_TO_DEG = 180.0 / PI; + + 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; + 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; + std::cout << fmt::format("PONI y: {:.3f} px = {:.4f} mm ({:+.3f} px from the header)", + g.GetBeamY_pxl(), g.GetBeamY_pxl() * pxl_mm, + g.GetBeamY_pxl() - header.GetBeamY_pxl()) << std::endl; + std::cout << fmt::format("Rot1: {:+.6f} rad = {:+.4f} deg ({:+.6f} rad from the header)", + g.GetPoniRot1_rad(), g.GetPoniRot1_rad() * RAD_TO_DEG, + g.GetPoniRot1_rad() - header.GetPoniRot1_rad()) << std::endl; + std::cout << fmt::format("Rot2: {:+.6f} rad = {:+.4f} deg ({:+.6f} rad from the header)", + g.GetPoniRot2_rad(), g.GetPoniRot2_rad() * RAD_TO_DEG, + g.GetPoniRot2_rad() - header.GetPoniRot2_rad()) << std::endl; + std::cout << fmt::format("Distance: {:.4f} mm ({:+.4f} mm from the header)", + g.GetDetectorDistance_mm(), + g.GetDetectorDistance_mm() - header.GetDetectorDistance_mm()) << std::endl; + // The PONI is the point of normal incidence, which is what the .poni file wants; the direct beam + // is where the beam actually lands, which is the number every other program reports as the beam + // centre. They part company by distance*tan(rot) as soon as the detector is tilted, so print both. + const auto [beam_x, beam_y] = g.GetDirectBeam_pxl(); + std::cout << fmt::format("Direct beam: {:.3f}, {:.3f} px", beam_x, beam_y) << std::endl; + + const std::string poni_path = output_prefix + ".poni"; + try { + WritePoniFile(poni_path, experiment, g); + } catch (const std::exception &e) { + logger.Error("{}", e.what()); + return 1; + } + logger.Info("Wrote {}", poni_path); + return 0; + } + configure_offline_output(experiment, output_prefix); // configure_offline_output applies the rugnux analysis defaults, one of which is the polarization // factor, so an explicit --polarization has to land after it or it is silently overwritten. diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 7e04b3f1..9c47f38e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -60,6 +60,7 @@ ADD_EXECUTABLE(jfjoch_test GoniometerAxisTest.cpp DetGeomCalibTest.cpp RingsFromProfileTest.cpp + CalibrationTest.cpp XtalOptimizerTest.cpp CrystalLatticeTest.cpp FPGAPTPTest.cpp diff --git a/tests/CalibrationTest.cpp b/tests/CalibrationTest.cpp new file mode 100644 index 00000000..068c44bf --- /dev/null +++ b/tests/CalibrationTest.cpp @@ -0,0 +1,81 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include +#include + +#include "../common/Definitions.h" +#include "../common/JFJochMath.h" +#include "../image_analysis/geom_refinement/Calibrants.h" +#include "../rugnux/RugnuxCalibration.h" + +TEST_CASE("Calibrants_LookupIsCaseInsensitive", "[DetGeomCalib]") { + CHECK(CalibrantRings("LaB6") == CalibrantRings("lab6")); + CHECK(CalibrantRings("AgBh") == CalibrantRings("agbh")); + CHECK(CalibrantRings("nonsense").empty()); +} + +// The innermost ring of a cubic standard is its (100), so the first q is 2*pi/a. This is what fixes +// the distance in GuessInitialGeometry, so a wrong table would put every calibration off by that scale. +TEST_CASE("Calibrants_CubicStandardsHaveTheirLatticeConstant", "[DetGeomCalib]") { + const std::map a_A = {{"lab6", LAB6_CELL_A}, {"ceo2", 5.4115}, {"si", 5.43102}}; + for (const auto &[name, a] : a_A) { + const auto q = CalibrantRings(name); + REQUIRE(!q.empty()); + CHECK(q.front() == Catch::Approx(2.0 * PI / a).epsilon(1e-5)); + } +} + +// Ice is the reason the calibrant abstraction is a ring list and not a UnitCell: its entries are +// measured ring positions, and enumerating hkl from the hexagonal cell would add rings that are +// systematically absent in P6_3/mmc. +TEST_CASE("Calibrants_IceIsTheMeasuredRingList", "[DetGeomCalib]") { + const auto q = CalibrantRings("ice"); + REQUIRE(q.size() == ICE_RING_RES_A.size()); + CHECK(std::is_sorted(q.begin(), q.end())); + CHECK(q.front() == Catch::Approx(2.0 * PI / ICE_RING_RES_A[0]).epsilon(1e-5)); // 3.895 A, the widest +} + +// pyFAI's Poni1 is the SLOW axis (rows, our y) and Poni2 the FAST axis (columns, our x), both in +// metres. Transposing them produces a file that is silently wrong, so pin the mapping with a geometry +// whose two axes differ. +TEST_CASE("Calibration_PoniFileAxisConvention", "[DetGeomCalib]") { + DiffractionExperiment x(DetJF4M()); + x.BeamX_pxl(1000.0f).BeamY_pxl(1275.0f).DetectorDistance_mm(150.0f); + + DiffractionGeometry geom = x.GetDiffractionGeometry(); + geom.PoniRot1_rad(0.01f).PoniRot2_rad(-0.02f); + + const std::string path = "poni_test.poni"; + WritePoniFile(path, x, geom); + + std::map keys; + std::ifstream f(path); + std::string line; + while (std::getline(f, line)) { + const auto colon = line.find(':'); + if (line.empty() || line[0] == '#' || colon == std::string::npos) + continue; + keys[line.substr(0, colon)] = line.substr(colon + 2); + } + f.close(); + std::remove(path.c_str()); + + const double pixel_m = geom.GetPixelSize_mm() * 1e-3; + CHECK(keys["poni_version"] == "2"); + CHECK(std::stod(keys["Poni1"]) == Catch::Approx(1275.0 * pixel_m)); // slow axis = y + CHECK(std::stod(keys["Poni2"]) == Catch::Approx(1000.0 * pixel_m)); // fast axis = x + CHECK(std::stod(keys["Distance"]) == Catch::Approx(0.150)); + CHECK(std::stod(keys["Rot1"]) == Catch::Approx(0.01)); + CHECK(std::stod(keys["Rot2"]) == Catch::Approx(-0.02)); + CHECK(std::stod(keys["Rot3"]) == Catch::Approx(0.0)); + CHECK(std::stod(keys["Wavelength"]) == Catch::Approx(geom.GetWavelength_A() * 1e-10)); + // max_shape is [rows, cols] - the same slow-then-fast order as Poni1/Poni2. + const std::string shape = "[" + std::to_string(x.GetYPixelsNumConv()) + ", " + + std::to_string(x.GetXPixelsNumConv()) + "]"; + CHECK(keys["Detector_config"].find(shape) != std::string::npos); +} diff --git a/tests/DetGeomCalibTest.cpp b/tests/DetGeomCalibTest.cpp index 7328bf87..b1ec526e 100644 --- a/tests/DetGeomCalibTest.cpp +++ b/tests/DetGeomCalibTest.cpp @@ -190,7 +190,7 @@ TEST_CASE("DetGeomCalib_GuessGeometry") { DiffractionGeometry geom_out; geom_out.Wavelength_A(1.0).DetectorDistance_mm(200.0); - GuessGeometry(geom_out, spots, UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90)); + GuessGeometry(geom_out, spots, CalculateXtalRings(UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90))); CHECK(fabsf(geom_out.GetBeamX_pxl() - geom.GetBeamX_pxl()) < 0.001f); CHECK(fabsf(geom_out.GetBeamY_pxl() - geom.GetBeamY_pxl()) < 0.001f); @@ -258,7 +258,7 @@ TEST_CASE("DetGeomCalib_OptimizeGeometry") { geom_i.Wavelength_A(1.0).BeamX_pxl(995.0).BeamY_pxl(1277.0) .DetectorDistance_mm(98).PoniRot1_rad(0.0975).PoniRot2_rad(0.055); - OptimizeGeometry(geom_i, spots, UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90)); + OptimizeGeometry(geom_i, spots, CalculateXtalRings(UnitCell(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90,90,90))); CHECK(geom_i.GetBeamX_pxl() == Catch::Approx(geom.GetBeamX_pxl())); CHECK(geom_i.GetBeamY_pxl() == Catch::Approx(geom.GetBeamY_pxl())); diff --git a/tests/RingsFromProfileTest.cpp b/tests/RingsFromProfileTest.cpp index 8bb02d59..c1efee43 100644 --- a/tests/RingsFromProfileTest.cpp +++ b/tests/RingsFromProfileTest.cpp @@ -11,6 +11,7 @@ namespace { constexpr UnitCell LAB6{LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90.0f, 90.0f, 90.0f}; +const std::vector LAB6_RINGS = CalculateXtalRings(LAB6); // A (q x azimuth) powder profile as the azimuthal integration would build it: the rings sit where // geom_true puts them, but every pixel is binned with geom_assumed - which is the whole point, since a @@ -23,7 +24,7 @@ std::vector SynthesiseProfile(const AzimuthalIntegrationMapping &mapping, const int32_t azim_bins = mapping.GetAzimuthalBinCount(); std::vector profile(static_cast(q_bins) * azim_bins, 100.0f); // flat background - for (const float q_ring : CalculateXtalRings(LAB6)) { + for (const float q_ring : LAB6_RINGS) { const float d = static_cast(2.0 * PI) / q_ring; if (d <= geom_true.GetWavelength_A() / 2.0f) continue; @@ -74,7 +75,7 @@ TEST_CASE("RingsFromProfile_RecoversBeamCenter", "[DetGeomCalib]") { .BeamY_pxl(geom_assumed.GetBeamY_pxl() - 4.0f); const auto profile = SynthesiseProfile(mapping, geom_assumed, geom_true); - const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom_assumed, LAB6); + const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom_assumed, LAB6_RINGS); // Several rings, sampled all the way round: without azimuthal coverage there is no centre to find. REQUIRE(rings.size() > 64); @@ -111,7 +112,7 @@ TEST_CASE("RingsFromProfile_RecoversTilt", "[DetGeomCalib]") { geom_true.PoniRot1_rad(0.02f).PoniRot2_rad(-0.015f); const auto profile = SynthesiseProfile(mapping, geom_assumed, geom_true); - const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom_assumed, LAB6); + const auto rings = RingsFromAzimuthalProfile(profile, mapping, geom_assumed, LAB6_RINGS); REQUIRE(rings.size() > 60); RingOptimizer optimizer(geom_assumed); @@ -132,7 +133,7 @@ TEST_CASE("RingsFromProfile_NeedsAzimuthalBins", "[DetGeomCalib]") { REQUIRE(mapping.GetAzimuthalBinCount() == 1); const std::vector profile(static_cast(mapping.GetQBinCount()), 1000.0f); - CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6).empty()); + CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6_RINGS).empty()); } // A profile with no rings in it must yield no ring points: the peak has to stand clear of the scatter @@ -150,5 +151,5 @@ TEST_CASE("RingsFromProfile_FlatProfileGivesNothing", "[DetGeomCalib]") { const std::vector profile( static_cast(mapping.GetQBinCount()) * mapping.GetAzimuthalBinCount(), 100.0f); - CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6).empty()); + CHECK(RingsFromAzimuthalProfile(profile, mapping, x.GetDiffractionGeometry(), LAB6_RINGS).empty()); } diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 5783b7da..dd801a37 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -178,7 +178,7 @@ TEST_CASE("RugnuxCommandLine_AzInt", "[process]") { const std::string cmd = RugnuxCommandLine(config, x, "in.h5"); CHECK(cmd.rfind("rugnux", 0) == 0); - CHECK(cmd.find("--azint-only") != std::string::npos); + CHECK(cmd.find("--mode azint") != std::string::npos); CHECK(cmd.find("--azim-phi-bins 4") != std::string::npos); CHECK(cmd.find("--azim-min-q") != std::string::npos); CHECK(cmd.find("in.h5") != std::string::npos); diff --git a/viewer/JFJochImageReadingWorker.cpp b/viewer/JFJochImageReadingWorker.cpp index 94a5a5da..adcf4f34 100644 --- a/viewer/JFJochImageReadingWorker.cpp +++ b/viewer/JFJochImageReadingWorker.cpp @@ -539,11 +539,12 @@ void JFJochImageReadingWorker::FindCenter(const UnitCell& calibrant, bool guess) logger.Info("Finding center"); DiffractionGeometry geom = current_image_ptr->Dataset().experiment.GetDiffractionGeometry(); + const std::vector ring_Q = CalculateXtalRings(calibrant); try { if (guess) - GuessGeometry(geom, current_image_ptr->ImageData().spots, calibrant); + GuessGeometry(geom, current_image_ptr->ImageData().spots, ring_Q); else - OptimizeGeometry(geom, current_image_ptr->ImageData().spots, calibrant); + OptimizeGeometry(geom, current_image_ptr->ImageData().spots, ring_Q); } catch (const JFJochException &e) { logger.ErrorException(e); return; @@ -560,8 +561,6 @@ void JFJochImageReadingWorker::FindCenter(const UnitCell& calibrant, bool guess) .PoniRot3_rad(geom.GetPoniRot3_rad()); UpdateDataset_i(new_experiment); - std::vector ring_Q = CalculateXtalRings(calibrant); - QVector rings; for (int i = 0; i < 15 && i < ring_Q.size(); i++) { rings.push_back(2 * PI / ring_Q[i]); diff --git a/viewer/widgets/PowderCalibrationWidget.cpp b/viewer/widgets/PowderCalibrationWidget.cpp index 8e24bf94..77c73c79 100644 --- a/viewer/widgets/PowderCalibrationWidget.cpp +++ b/viewer/widgets/PowderCalibrationWidget.cpp @@ -8,6 +8,7 @@ #include "../../common/JFJochMath.h" #include "../image_analysis/geom_refinement/AssignSpotsToRings.h" +#include "../image_analysis/geom_refinement/Calibrants.h" PowderCalibrationWidget::PowderCalibrationWidget(QWidget *parent) : QWidget(parent) { auto layout = new QVBoxLayout(this); @@ -45,37 +46,27 @@ PowderCalibrationWidget::PowderCalibrationWidget(QWidget *parent) : QWidget(pare } UnitCell PowderCalibrationWidget::GetCalibrant() const { - UnitCell uc(LAB6_CELL_A, LAB6_CELL_A, LAB6_CELL_A, 90, 90, 90); - switch (calibrantCombo->currentIndex()) { - case 1: - // T. C. Huang , H. Toraya, T. N. Blanton, Y. Wu, J. Appl. Cryst. 26 (1993), 180-184. - uc = UnitCell(5.1769, 4.7218, 58.380, 89.440, 89.634, 75.854); - break; - case 2: - // CeO2, cubic fluorite (a = 5.4115 Å), NIST SRM 674b. - uc = UnitCell(5.4115, 5.4115, 5.4115, 90, 90, 90); - break; - case 3: - // Silicon, cubic (a = 5.43102 Å), NIST SRM 640. - uc = UnitCell(5.43102, 5.43102, 5.43102, 90, 90, 90); - break; - case 4: - if (sample_cell) - uc = sample_cell.value(); - break; - default: - break; + // The combo lists the shared calibrant table (in its order), then the current sample's own cell. + // Ice is skipped: it is a ring list rather than a cell, and this panel works in cells. + int idx = calibrantCombo->currentIndex(); + for (const auto &c : Calibrants()) { + if (!c.cell) + continue; + if (idx == 0) + return *c.cell; + idx--; } - return uc; + if (sample_cell) + return sample_cell.value(); + return *Calibrants().front().cell; } void PowderCalibrationWidget::updateCalibrantList() { calibrantCombo->clear(); - calibrantCombo->addItem("LaB6"); - calibrantCombo->addItem("Silver Behenate"); - calibrantCombo->addItem("CeO2"); - calibrantCombo->addItem("Silicon"); + for (const auto &c : Calibrants()) + if (c.cell) + calibrantCombo->addItem(QString::fromStdString(c.name)); if (sample_cell) calibrantCombo->addItem(QString("Current sample (%1 %2 %3 %4 %5 %6)") .arg(QString::number(sample_cell->a, 'f', 1))