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))