Stills geometry: refine the free cell angles, -C as a prior and a hypothesis, place modules from the fitted geometry
Three defects of the stills global geometry/cell fit, each moving the detector distance: - The cell angles were held at the start cell (InitCellAngles, "constant in v1") while the lengths and the distance were free, so a wrong angle went into the lengths and the distance: a monoclinic cell given with beta 4 degrees off refined to lengths and a distance both ~4% long, and a cell given with beta held at 90 on a crystal at ~90.9 left the distance 1.3% short. The angles the crystal system leaves free (monoclinic beta, triclinic all three) are now a refined block. The cell given is a Gaussian prior (1% on a length, 1 degree on an angle) weighed against the spots' ~1 px scatter, replacing the unused length regularizer; the refined-vs-start deviation (scale, lengths, angles) and the distance move are logged. A cell move, not only a distance move, now triggers the re-indexing refit. - The module placement started from the file's geometry: on a run whose header distance was 3% short and centre 6 px off, the module scan absorbed the centre into every module and the placement pushed two second-ring modules 13 px outwards - a radius-dependent push the gauge does not hold, which trades with the distance/cell-scale curvature - and the distance settled 0.8% long. The placement now starts from a fit of the same frames with the modules held, each frame's lattice taken from it (GeometryRefinerResult::lattices). Same data, right cell: the distance went from 1.0% long to the value a reference model's R-free puts it at, the modules moved < 1.4 px. - A cell given with the wrong shape still indexes the frames at a distance of its own and the fit cannot leave it. With -C on stills the geometry refinement is now run both with the cell given and de novo, from the file's geometry, and each is indexed on the same 96-frame sample with its own refined cell and geometry. The cell given is kept where it puts >= 0.9x as many spots on the lattice, or where its fit puts the detector within 2% of the de novo fit (a crystal form of the sample, possibly a minority one the user asked for); otherwise the run is processed as without -C, with a warning. `rugnux -C <wrong cell> file` now gives the result of `rugnux file`. Costs one more geometry-refinement sample on -C stills runs (~+30-50 s on a 4000-frame run). Docs: RUGNUX_CLI (-C, --refine-geometry), CPU_DATA_ANALYSIS_INDEXING 7.7, RUGNUX_MODEL_AND_REFERENCE. Test: RefineGlobalGeometry_MonoclinicBeta. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK (cherry picked from commit fe442ff33aeec7569de1fdf90076dba6b39a8e8f)
This commit is contained in:
leonarski_f committed
1 parent
09707900f5
commit
0c74e6b666
10 files changed
+313
-96
No files matched your search
@@ -412,13 +412,28 @@ The extraction window around a ring is capped at half the gap to its neighbour,
|
||||
|
||||
### 7.7 Stills: a given cell and the detector modules
|
||||
|
||||
**A given cell is a starting point.** On stills a cell and a distance are degenerate - a spot's
|
||||
radius fixes only their ratio. A cell taken from a deposited structure is routinely a percent off the
|
||||
crystal at hand (another temperature, another batch), while a facility's stated distance is usually
|
||||
better than that, so the stills geometry refinement holds the distance and refines the cell. Only where
|
||||
the frames have rejected the file's distance (§1.4), or the file states the distance only nominally, is
|
||||
the given cell anchored and the distance refined from it. Frames are matched to the given cell in any
|
||||
setting of the same lattice, and are described on its axes.
|
||||
**A given cell is a prior, and a hypothesis.** On stills a cell and a distance are degenerate near the
|
||||
beam - a spot's radius fixes only their ratio - and are told apart only by the curvature of
|
||||
$r=D\tan2\theta$ against $\sin\theta$ at high angle. The stills geometry refinement therefore refines
|
||||
the beam centre, the distance and the cell together over the strong frames: the cell lengths and the
|
||||
angles its crystal system leaves free (the monoclinic $\beta$, all three triclinic angles), each held to
|
||||
the given cell by a Gaussian prior of 1 % on a length and 1° on an angle, weighed against the spots'
|
||||
scatter of about a pixel. On a bundle of strong frames the spots outweigh the prior wherever they
|
||||
determine the cell, so it decides only what they leave open. Holding the angles instead puts a wrong
|
||||
angle into the lengths and the distance: a monoclinic cell given with $\beta$ 4° off was taken up as a
|
||||
cell and a distance both 4 % long. The log states how far the refined cell is from the one it started
|
||||
from (scale, each length, each angle) and how far the distance moved. A distance given with
|
||||
`--detector-distance` is held, and the cell alone takes up its error. Frames are matched to the given
|
||||
cell in any setting of the same lattice, and are described on its axes.
|
||||
|
||||
With `-C` the cell given is also tested as a hypothesis against the cell the frames index to de novo.
|
||||
Both are fitted from the file's geometry - the cell given as above, the de novo one as without `-C` - and
|
||||
each, with its refined cell and geometry, indexes the same sample of frames. The cell given is kept where
|
||||
it puts at least 0.9 times as many of the sample's spots on the lattice, and also where it does not but
|
||||
its fit puts the detector within 2 % of where the de novo fit puts it: the detector is the same for every
|
||||
crystal of the run, so such a cell is a crystal form of the sample (possibly a minority one, which the
|
||||
user may be after). A cell that fits only with a detector of its own describes no crystal of the sample,
|
||||
and the run is then processed as without `-C`, with a warning.
|
||||
|
||||
**One crystal form.** A sample can hold two crystal forms, and a known-cell indexer fits both to the
|
||||
given cell's orientation search; their refined cells then fall into two volume populations, and
|
||||
@@ -457,7 +472,11 @@ fixed and with them placed, and scored on the other. The placement is used if it
|
||||
spots better as a whole, and then each module keeps its move only if the held-out spots on that module
|
||||
are predicted better; the others go back to where they were given and the placement is fitted again.
|
||||
On a detector of many modules the outer ones, reached only by a few strong high-resolution spots, are
|
||||
the ones this decides on. `--panel-refine` keeps every move without the tests, `--panel-refine=off` none. The modules move only relative to one another: a
|
||||
the ones this decides on. `--panel-refine` keeps every move without the tests, `--panel-refine=off` none. The placement starts from the beam, the distance and the cell fitted on the same frames with the modules held,
|
||||
each frame's lattice taken from that fit, not from the file's geometry: from a start off in distance and centre
|
||||
the search below moved every module by the centre's error and the placement then pushed two modules of the
|
||||
second ring outwards, where the distance's error had left their spots, and the distance settled 0.8 % long.
|
||||
The modules move only relative to one another: a
|
||||
gauge holds their spot-weighted mean shift, turn and outward push, and with the detector tilt refined
|
||||
also the two tilt-like gradients across the detector, at zero, so the beam centre, distance, tilt and
|
||||
crystal orientations carry what all modules share. A coarse search of each module's position comes
|
||||
|
||||
+2
-2
@@ -264,7 +264,7 @@ For very long axes the direction grid's angular resolution binds well below 500
|
||||
| `--force-still` | Treat a rotation (goniometer) dataset as independent stills instead of rotation |
|
||||
| `-X, --indexing-algorithm <txt>` | `Auto` (default) \| `FFT` \| `cuFFT` \| `FFTW` \| `FFBIDX` \| `None`, case-insensitive. `FFT` is the de-novo FFT search, on the GPU (cuFFT) where there is one and `--no-gpu` is not given, otherwise on the CPU (FFTW); `cuFFT` and `FFTW` name one implementation each; `cuFFT` needs a GPU and is refused without one. `FFBIDX` (known cell, `-C`) runs on the GPU where there is one and otherwise on the CPU. `Auto` takes `FFBIDX` on stills with a known cell, and `FFT` otherwise; rotation data always use the FFT search |
|
||||
| `-F` | Short for `-X FFT`, kept for old command lines |
|
||||
| `-C, --unit-cell <cell>` | Reference unit cell `"a,b,c,alpha,beta,gamma"` (required by `ffbidx`). On rotation data it also widens the FFT search to cover the cell given, at both ends: the longest axis looked for is raised to reach it, and `--fft-min-unit-cell` is lowered to admit it. A candidate that is an integer supercell of the cell given (a twin's coincidence lattice, a row taken at a multiple of its length) is stepped down to it |
|
||||
| `-C, --unit-cell <cell>` | Reference unit cell `"a,b,c,alpha,beta,gamma"` (required by `ffbidx`). On rotation data it also widens the FFT search to cover the cell given, at both ends: the longest axis looked for is raised to reach it, and `--fft-min-unit-cell` is lowered to admit it. A candidate that is an integer supercell of the cell given (a twin's coincidence lattice, a row taken at a multiple of its length) is stepped down to it. On stills the geometry refinement takes it as a prior (1 % on the lengths, 1° on the angles the crystal system leaves free) and tests it against the cell the frames index to de novo: a cell that fits the frames only with a detector distance of its own (more than 2 % from the de novo fit's) and puts fewer than 0.9 times as many sampled spots on the lattice is dropped with a warning, and the run is processed as without `-C` ([§7.7](CPU_DATA_ANALYSIS_INDEXING.md#77-stills-a-given-cell-and-the-detector-modules)) |
|
||||
| `--fft-min-unit-cell <num>` | Shortest cell axis the FFT search accepts, Å (default: 10). A candidate with a shorter axis is discarded — but on rotation data a row shorter than the floor, down to 3 Å, is recognised under the multiple of it the search reports past the floor, so a small-molecule cell is indexed without this flag. `-C` lowers the floor on its own to cover the cell given. On stills a crystal below the floor cannot be indexed unless this is lowered |
|
||||
| `--min-indexed-spots <num>` | Spots a frame must have on the lattice before it counts as indexed (default: 6, minimum 4). It sets the reported indexing rate and the count the rotation first pass ranks candidate lattices by; whether the run has a lattice at all is decided on the pooled spots instead ([§4.1](CPU_DATA_ANALYSIS_INDEXING.md#41-indexed-spot-decision-inlier-test)), and integration is not gated by it |
|
||||
| `-S, --space-group <num\|symbol>` | Space group number (`96`) or Hermann-Mauguin symbol (`P43212`) — for indexing and scaling. On serial stills whose group admits two indexings, it also lets the crystals be put on one from the data ([the indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity)) |
|
||||
@@ -273,7 +273,7 @@ For very long axes the direction grid's angular resolution binds well below 500
|
||||
| `--force-rotation-lattice <vec>` | Force rotation lattice (9 floats, Å), skipping the first pass |
|
||||
| `--rotation-no-postrefine` | Rotation: disable the default-on two-pass geometry post-refine (see [Rotation data](RUGNUX_TUTORIAL.md#rotation-data)) |
|
||||
| `--prepass-fraction <f>` | Rotation: the two-pass geometry pre-pass integrates the first fraction `f` of the sweep, but at least 120° of it (0.05-1, default **0.5**; `1` integrates the whole sweep twice). The second pass always integrates every frame. The pre-scan and the first-pass indexing sample the same part of the sweep |
|
||||
| `--refine-geometry[=N\|off]` | Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default 200) then re-indexes; default ON for stills, a stationary goniometer axis counting as stills, `=off` disables. The bundle anchors on the reference cell (`-C` / `-z`); without one, on the cell most of the sampled frames index to de novo, which is used for the geometry steps only (the run is still indexed de novo) |
|
||||
| `--refine-geometry[=N\|off]` | Stills: extra first pass that bundle-adjusts the shared beam/distance/cell from N strongly-indexed frames (default 200) then re-indexes; default ON for stills, a stationary goniometer axis counting as stills, `=off` disables. The cell's lengths and the angles its crystal system leaves free are refined with the distance; the bundle anchors them on the reference cell (`-C` / `-z`) by a prior of 1 % / 1°; without one, on the cell most of the sampled frames index to de novo, which is used for the geometry steps only (the run is still indexed de novo) |
|
||||
| `--index-ice-rings[=on\|off]` | Index on the spots flagged as sitting on an ice ring too, instead of setting them aside (default: **off**; no effect without `--detect-ice-rings`, which does the flagging) |
|
||||
|
||||
Indexer choice in brief: `ffbidx` (GPU or CPU) refines toward a **known cell** and is best for sparse
|
||||
|
||||
@@ -58,7 +58,7 @@ the rest walks through a run step by step.
|
||||
| `--reference-column <label>` | An MTZ column label | Which column of the reference is read | Through `-z` only |
|
||||
| `--free-flags <file>` | Merged data carrying an R-free set: an MTZ, or a structure-factor mmCIF - typically the deposited reflections of the `--model` | Its R-free set and nothing else: the model is validated and refined against it, on the model's axes, where the merge matches the file's intensities (`R_FREE_FLAGS= FILE`) | The free flags only |
|
||||
| `-S, --space-group` | A space group | Fixes it: no space-group search | **Yes**: the group is yours (a fitting model may still relabel its hand) |
|
||||
| `-C, --unit-cell` | A cell | Indexing looks for this cell (within tolerance); the FFT search range is widened to cover it; on stills the shared-geometry refinement anchors on it (without it, on the cell the frames index to) and, with `-S` and `--model`, the per-image use of the model | **Yes**: the cell, and the axis order the files are written in |
|
||||
| `-C, --unit-cell` | A cell | Indexing looks for this cell (within tolerance); the FFT search range is widened to cover it; on stills the shared-geometry refinement anchors on it as a prior (without it, on the cell the frames index to), and drops it, as without `-C`, where it fits the frames only with a detector distance of its own and, with `-S` and `--model`, the per-image use of the model | **Yes**: the cell, and the axis order the files are written in |
|
||||
| `--force-rotation-lattice` | Nine numbers: a lattice | Skips the first-pass indexing and uses this lattice | **Yes**: the lattice |
|
||||
| `--report-resolution <dmin>[,<dmax>]` | Another program's resolution range | A second statistics table at that range (`REFRES_*`) | No — report only |
|
||||
|
||||
|
||||
@@ -22,8 +22,8 @@ gemmi::CrystalSystem MapSystem(gemmi::CrystalSystem s) {
|
||||
}
|
||||
|
||||
// The cell angles (radians) that XtalResidual bakes into its B matrix for a given system. Only
|
||||
// monoclinic (beta) and triclinic (all three) read p2; the rest use fixed symmetry angles. We keep the
|
||||
// angles constant in v1, so these are set once from the input cell.
|
||||
// monoclinic (beta) and triclinic (all three) read p2, and only there are the angles refined (CellAnglesFree);
|
||||
// the rest use fixed symmetry angles.
|
||||
void InitCellAngles(gemmi::CrystalSystem system, const UnitCell &cell, double ang[3]) {
|
||||
const double d2r = PI / 180.0;
|
||||
switch (system) {
|
||||
@@ -43,6 +43,10 @@ void InitCellAngles(gemmi::CrystalSystem system, const UnitCell &cell, double an
|
||||
}
|
||||
}
|
||||
|
||||
bool CellAnglesFree(gemmi::CrystalSystem system) {
|
||||
return system == gemmi::CrystalSystem::Monoclinic || system == gemmi::CrystalSystem::Triclinic;
|
||||
}
|
||||
|
||||
// Effective (lengths, alpha,beta,gamma) that reproduce the lattice XtalResidual builds from (len, ang),
|
||||
// so a frame's lattice can be reconstructed with AngleAxisAndCellToLattice for HKL re-assignment.
|
||||
void EffectiveCell(gemmi::CrystalSystem system, const double len[3], const double ang[3],
|
||||
@@ -116,17 +120,21 @@ UnitCell BuildUnitCell(gemmi::CrystalSystem system, const double len[3], const d
|
||||
return uc;
|
||||
}
|
||||
|
||||
// Anchors the cell lengths to the input cell. Penalises each length's deviation from its input value
|
||||
// with a large weight, breaking the low-resolution distance<->cell-scale degeneracy (position ~
|
||||
// distance / cell_scale) so distance is determined by the data with the cell pinned.
|
||||
struct CellLengthRegularizer {
|
||||
double w0, w1, w2;
|
||||
double l0, l1, l2;
|
||||
// The prior on the cell (GeometryRefinerSettings::cell_length_sigma): each length and angle a Gaussian about
|
||||
// its input value. A spot's residual is a reciprocal-space distance, scattered by about a pixel's worth of it,
|
||||
// so the prior's residual is (value - input) / sigma in units of that pixel - the prior then weighs in the sum
|
||||
// as it would in a likelihood with the spots' scatter one pixel. Lengths and angles the crystal system does
|
||||
// not use (b of a tetragonal cell, alpha and gamma of a monoclinic one) are held at the input by it alone.
|
||||
struct CellPrior {
|
||||
double len0[3], ang0[3];
|
||||
double w_len[3], w_ang;
|
||||
|
||||
template <typename T>
|
||||
bool operator()(const T *const len, T *r) const {
|
||||
r[0] = T(w0) * (len[0] - T(l0));
|
||||
r[1] = T(w1) * (len[1] - T(l1));
|
||||
r[2] = T(w2) * (len[2] - T(l2));
|
||||
bool operator()(const T *const len, const T *const ang, T *r) const {
|
||||
for (int j = 0; j < 3; ++j) {
|
||||
r[j] = T(w_len[j]) * (len[j] - T(len0[j]));
|
||||
r[3 + j] = T(w_ang) * (ang[j] - T(ang0[j]));
|
||||
}
|
||||
return true;
|
||||
}
|
||||
};
|
||||
@@ -169,11 +177,11 @@ struct BeamBlock {
|
||||
|
||||
// The residual of one spot (XtalResidual) on the blocks this refinement can move, in this order: the beam
|
||||
// (read through BeamBlock) and the distance as one block, the tilt where it is refined, the frame's orientation,
|
||||
// the cell lengths, and the module's move (dx, dy, rot) where the spot's module is placed. The rotation axis, the
|
||||
// cell angles and a held tilt enter as plain numbers.
|
||||
// the cell lengths, the cell angles, and the module's move (dx, dy, rot) where the spot's module is placed. The
|
||||
// rotation axis and a held tilt enter as plain numbers.
|
||||
//
|
||||
// The residual is the difference of two halves that share no block: the observed reciprocal vector (beam,
|
||||
// distance, tilt, move) and the predicted one (orientation, cell). Ceres' AutoDiffCostFunction would carry the
|
||||
// distance, tilt, move) and the predicted one (orientation, cell lengths and angles). Ceres' AutoDiffCostFunction would carry the
|
||||
// derivatives of every block through both halves; here each half is differentiated on its own blocks only, which
|
||||
// gives the same values and the same derivatives with half as many of them per number.
|
||||
//
|
||||
@@ -183,10 +191,9 @@ struct BeamBlock {
|
||||
class SpotCost : public ceres::CostFunction {
|
||||
public:
|
||||
SpotCost(const XtalResidual &residual, const BeamBlock &beam_block, const double *detector_rot,
|
||||
const double *rotation_axis, const double *cell_ang, bool tilt, bool module, double ux = 0,
|
||||
double uy = 0)
|
||||
const double *rotation_axis, bool tilt, bool module, double ux = 0, double uy = 0)
|
||||
: residual(residual), beam_block(beam_block), detector_rot(detector_rot), rotation_axis(rotation_axis),
|
||||
cell_ang(cell_ang), tilt(tilt), module(module), ux(ux), uy(uy) {
|
||||
tilt(tilt), module(module), ux(ux), uy(uy) {
|
||||
set_num_residuals(3);
|
||||
std::vector<int32_t> &sizes = *mutable_parameter_block_sizes();
|
||||
sizes.push_back(3);
|
||||
@@ -194,6 +201,7 @@ public:
|
||||
sizes.push_back(2);
|
||||
sizes.push_back(3);
|
||||
sizes.push_back(3);
|
||||
sizes.push_back(3);
|
||||
if (module)
|
||||
sizes.push_back(3);
|
||||
}
|
||||
@@ -211,8 +219,9 @@ private:
|
||||
bool EvaluateWith(double const *const *p, double *out, double **jac) const {
|
||||
constexpr int N_OBS = 3 + (TILT ? 2 : 0) + (MODULE ? 3 : 0);
|
||||
using JO = ceres::Jet<double, N_OBS>;
|
||||
using JP = ceres::Jet<double, 6>;
|
||||
constexpr int TILT_BLOCK = 1, P0_BLOCK = TILT ? 2 : 1, P1_BLOCK = P0_BLOCK + 1, MOVE_BLOCK = P1_BLOCK + 1;
|
||||
using JP = ceres::Jet<double, 9>;
|
||||
constexpr int TILT_BLOCK = 1, P0_BLOCK = TILT ? 2 : 1, P1_BLOCK = P0_BLOCK + 1, P2_BLOCK = P1_BLOCK + 1,
|
||||
MOVE_BLOCK = P2_BLOCK + 1;
|
||||
|
||||
// The observed half, on the beam and distance, the tilt and the move.
|
||||
JO bd[3], t[2], m[3];
|
||||
@@ -230,13 +239,13 @@ private:
|
||||
else
|
||||
Observed(bd, detector_rot, MODULE ? m : nullptr, obs);
|
||||
|
||||
// The predicted half, on the orientation and the cell lengths.
|
||||
JP p0[3], p1[3];
|
||||
// The predicted half, on the orientation and the cell.
|
||||
JP p0[3], p1[3], p2[3];
|
||||
for (int j = 0; j < 3; ++j) {
|
||||
p0[j] = JP(p[P0_BLOCK][j], j);
|
||||
p1[j] = JP(p[P1_BLOCK][j], 3 + j);
|
||||
p2[j] = JP(p[P2_BLOCK][j], 6 + j);
|
||||
}
|
||||
const JP p2[3] = {JP(cell_ang[0]), JP(cell_ang[1]), JP(cell_ang[2])};
|
||||
JP recip_unrot[3];
|
||||
residual.PredictedRecipUnrot(p1, p2, recip_unrot);
|
||||
JP pred[3];
|
||||
@@ -261,6 +270,7 @@ private:
|
||||
fill(jac[MOVE_BLOCK], 3, [&](int r, int c) { return obs[r].v[N_OBS - 3 + c]; });
|
||||
fill(jac[P0_BLOCK], 3, [&](int r, int c) { return -pred[r].v[c]; });
|
||||
fill(jac[P1_BLOCK], 3, [&](int r, int c) { return -pred[r].v[3 + c]; });
|
||||
fill(jac[P2_BLOCK], 3, [&](int r, int c) { return -pred[r].v[6 + c]; });
|
||||
return true;
|
||||
}
|
||||
|
||||
@@ -283,7 +293,6 @@ private:
|
||||
const BeamBlock beam_block;
|
||||
const double *detector_rot; // the tilt where it is held
|
||||
const double *rotation_axis;
|
||||
const double *cell_ang;
|
||||
const bool tilt, module;
|
||||
const double ux, uy;
|
||||
};
|
||||
@@ -391,11 +400,13 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge
|
||||
y = p[1];
|
||||
};
|
||||
|
||||
// Shared cell blocks. Lengths refined (regularized); angles constant.
|
||||
// Shared cell blocks: the lengths, and the angles where the system leaves any free; both held to the
|
||||
// input cell by its prior (CellPrior).
|
||||
double cell_len[3] = {input_cell.a, input_cell.b, input_cell.c};
|
||||
const double cell_len0[3] = {input_cell.a, input_cell.b, input_cell.c};
|
||||
double cell_ang[3];
|
||||
InitCellAngles(system, input_cell, cell_ang);
|
||||
const double cell_ang0[3] = {cell_ang[0], cell_ang[1], cell_ang[2]};
|
||||
|
||||
double eff_len[3], eff_alpha, eff_beta, eff_gamma; // recomputed each round from cell_len/ang
|
||||
|
||||
@@ -526,21 +537,22 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge
|
||||
const auto &m = modules[s.module];
|
||||
const XtalResidual found(s.x, s.y, lambda, pixel, cos_rot3, sin_rot3, 0.0, h, k, l, system,
|
||||
orientation);
|
||||
auto *cost = new SpotCost(found, beam_block, detector_rot, rot_vec, cell_ang, refine_tilt, true,
|
||||
auto *cost = new SpotCost(found, beam_block, detector_rot, rot_vec, refine_tilt, true,
|
||||
s.x - m.x0, s.y - m.y0);
|
||||
if (refine_tilt)
|
||||
problem.AddResidualBlock(cost, &loss, beam_distance, detector_rot, orient[i].data(),
|
||||
cell_len, move[s.module].data());
|
||||
cell_len, cell_ang, move[s.module].data());
|
||||
else
|
||||
problem.AddResidualBlock(cost, &loss, beam_distance, orient[i].data(), cell_len,
|
||||
move[s.module].data());
|
||||
cell_ang, move[s.module].data());
|
||||
} else {
|
||||
auto *cost = new SpotCost(xr, beam_block, detector_rot, rot_vec, cell_ang, refine_tilt, false);
|
||||
auto *cost = new SpotCost(xr, beam_block, detector_rot, rot_vec, refine_tilt, false);
|
||||
if (refine_tilt)
|
||||
problem.AddResidualBlock(cost, &loss, beam_distance, detector_rot, orient[i].data(),
|
||||
cell_len);
|
||||
cell_len, cell_ang);
|
||||
else
|
||||
problem.AddResidualBlock(cost, &loss, beam_distance, orient[i].data(), cell_len);
|
||||
problem.AddResidualBlock(cost, &loss, beam_distance, orient[i].data(), cell_len,
|
||||
cell_ang);
|
||||
}
|
||||
++used_in_frame;
|
||||
++spots_used;
|
||||
@@ -557,21 +569,22 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge
|
||||
break; // keep the previous round's solution
|
||||
}
|
||||
|
||||
// Cell-length regularizer (weight ~ coeff * sqrt(N_obs) / L0^2 -> effectively pins the cell).
|
||||
const double sqrtN = std::sqrt(static_cast<double>(spots_used));
|
||||
if (settings.refine_beam_distance) {
|
||||
CellLengthRegularizer *reg = new CellLengthRegularizer{
|
||||
settings.cell_reg_coeff * sqrtN / (cell_len0[0] * cell_len0[0]),
|
||||
settings.cell_reg_coeff * sqrtN / (cell_len0[1] * cell_len0[1]),
|
||||
settings.cell_reg_coeff * sqrtN / (cell_len0[2] * cell_len0[2]),
|
||||
cell_len0[0], cell_len0[1], cell_len0[2]};
|
||||
problem.AddResidualBlock(
|
||||
new ceres::AutoDiffCostFunction<CellLengthRegularizer, 3, 3>(reg), nullptr, cell_len);
|
||||
auto *prior = new CellPrior;
|
||||
for (int j = 0; j < 3; ++j) {
|
||||
prior->len0[j] = cell_len0[j];
|
||||
prior->ang0[j] = cell_ang0[j];
|
||||
prior->w_len[j] = recip_per_px / (settings.cell_length_sigma * cell_len0[j]);
|
||||
}
|
||||
prior->w_ang = recip_per_px / (settings.cell_angle_sigma_deg * PI / 180.0);
|
||||
problem.AddResidualBlock(new ceres::AutoDiffCostFunction<CellPrior, 6, 3, 3>(prior), nullptr,
|
||||
cell_len, cell_ang);
|
||||
}
|
||||
|
||||
// Shared blocks go in the second elimination group (those the problem has; constant ones are
|
||||
// stripped by Ceres).
|
||||
for (double *b : {beam_distance, detector_rot, cell_len})
|
||||
for (double *b : {beam_distance, detector_rot, cell_len, cell_ang})
|
||||
if (problem.HasParameterBlock(b))
|
||||
ordering->AddElementToGroup(b, 1);
|
||||
|
||||
@@ -617,7 +630,7 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge
|
||||
}
|
||||
|
||||
// Beam + distance bounds around the nominal value; cell lengths loosely bounded (the
|
||||
// regularizer does the anchoring).
|
||||
// prior does the anchoring).
|
||||
if (!orientation_only) {
|
||||
problem.SetParameterLowerBound(beam_distance, 0, beam0[0] - settings.beam_range_px);
|
||||
problem.SetParameterUpperBound(beam_distance, 0, beam0[0] + settings.beam_range_px);
|
||||
@@ -635,8 +648,13 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge
|
||||
if (!settings.refine_beam_distance) {
|
||||
problem.SetParameterBlockConstant(beam_distance);
|
||||
problem.SetParameterBlockConstant(cell_len);
|
||||
} else if (!settings.refine_distance)
|
||||
problem.SetManifold(beam_distance, new ceres::SubsetManifold(3, {2}));
|
||||
problem.SetParameterBlockConstant(cell_ang);
|
||||
} else {
|
||||
if (!CellAnglesFree(system))
|
||||
problem.SetParameterBlockConstant(cell_ang);
|
||||
if (!settings.refine_distance)
|
||||
problem.SetManifold(beam_distance, new ceres::SubsetManifold(3, {2}));
|
||||
}
|
||||
}
|
||||
if (refine_tilt) {
|
||||
const double tilt_range = 2.0 * PI / 180.0;
|
||||
@@ -682,6 +700,7 @@ GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_ge
|
||||
for (size_t i = 0; i < frames.size(); ++i) {
|
||||
const CrystalLattice latt = AngleAxisAndCellToLattice(
|
||||
orient[i].data(), eff_len, eff_alpha, eff_beta, eff_gamma);
|
||||
result.lattices.push_back(latt);
|
||||
const Coord a = latt.Astar(), b = latt.Bstar(), c = latt.Cstar();
|
||||
const Coord v0 = latt.Vec0(), v1 = latt.Vec1(), v2 = latt.Vec2();
|
||||
for (const auto &s : frames[i].spots) {
|
||||
|
||||
@@ -15,11 +15,11 @@
|
||||
|
||||
// Offline global (multi-frame) geometry refinement for serial stills. Pass 1 indexes many frames
|
||||
// independently; this joint bundle adjustment then determines the ONE shared detector geometry (beam
|
||||
// centre, distance, cell scale) that per-image refinement cannot reach from a single sparse still, so
|
||||
// centre, distance, cell) that per-image refinement cannot reach from a single sparse still, so
|
||||
// re-indexing with it converts more frames. Mirrors the validated Python prototype: one Ceres problem
|
||||
// with shared beam/distance/cell blocks and a per-frame orientation block, robust loss, and a cell
|
||||
// regularizer that anchors the (known) cell to break the distance<->cell-scale degeneracy at low
|
||||
// resolution. The detector tilt (rot1/rot2) is refined only on request (refine_tilt), and then with the beam
|
||||
// with shared beam/distance/cell blocks and a per-frame orientation block, robust loss, and a prior
|
||||
// on the input cell (GeometryRefinerSettings::cell_length_sigma) that settles the distance<->cell-scale
|
||||
// degeneracy where the spots leave it open. The detector tilt (rot1/rot2) is refined only on request (refine_tilt), and then with the beam
|
||||
// written as the direct beam, which the tilt does not move.
|
||||
|
||||
struct GeomRefineSpot {
|
||||
@@ -39,9 +39,15 @@ struct GeomRefineFrame {
|
||||
|
||||
struct GeometryRefinerSettings {
|
||||
gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic;
|
||||
// Cell-scale regularizer strength (anchors the cell lengths to the input cell so the otherwise
|
||||
// degenerate distance is determined). Large => cell effectively fixed (the validated safe regime).
|
||||
double cell_reg_coeff = 30.0;
|
||||
// The input cell is a prior on the refined one, not a constraint: its lengths and the angles the crystal
|
||||
// system leaves free (beta of a monoclinic cell, all three of a triclinic one) are refined, each held to
|
||||
// the input value by a Gaussian of this width - relative for the lengths, in degrees for the angles -
|
||||
// weighed against the spots' scatter of about a pixel. A sigma of 1% / 1 degree covers what a cell from
|
||||
// another crystal, temperature or setting is usually off by. On a bundle of strong frames the spots
|
||||
// outweigh it by orders of magnitude wherever they determine the cell, so it decides only what the spots
|
||||
// leave undetermined: the cell scale near the beam, where it trades with the distance.
|
||||
double cell_length_sigma = 0.01;
|
||||
double cell_angle_sigma_deg = 1.0;
|
||||
// With the distance held, the cell is no longer degenerate with it and is refined from the spots.
|
||||
bool refine_distance = true;
|
||||
int rounds = 3; // hkl-reassignment / tolerance-tightening rounds
|
||||
@@ -90,10 +96,12 @@ struct GeometryRefinerResult {
|
||||
// Per frame, per spot (as frames[i].spots): where the refined geometry predicts the spot's reflection on
|
||||
// the detector, as the spot would be found (before its module's move); NaN where it predicts none.
|
||||
std::vector<std::vector<std::array<float, 2>>> predicted_px;
|
||||
// Per frame (as frames): its lattice at the end, the frame's orientation with the refined cell.
|
||||
std::vector<CrystalLattice> lattices;
|
||||
};
|
||||
|
||||
// nominal_geom supplies the starting beam/distance/tilt/wavelength/pixel size; input_cell is anchored by
|
||||
// the regularizer. Returns ok=false (and leaves the caller's geometry untouched) on any failure.
|
||||
// nominal_geom supplies the starting beam/distance/tilt/wavelength/pixel size; input_cell is where the cell
|
||||
// starts and the centre of its prior. Returns ok=false (and leaves the caller's geometry untouched) on any failure.
|
||||
GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_geom,
|
||||
const UnitCell &input_cell,
|
||||
const std::vector<GeomRefineFrame> &frames,
|
||||
|
||||
@@ -241,15 +241,16 @@ void RunModuleDiagnostics(const ModuleDiagnosticsInput &in) {
|
||||
}
|
||||
}
|
||||
|
||||
// The cell free (a weak anchor) instead of pinned.
|
||||
// The cell pinned (a prior a hundred times narrower) instead of refined.
|
||||
{
|
||||
GeometryRefinerSettings s2 = st;
|
||||
s2.cell_reg_coeff = 0.3;
|
||||
s2.cell_length_sigma = st.cell_length_sigma / 100.0;
|
||||
s2.cell_angle_sigma_deg = st.cell_angle_sigma_deg / 100.0;
|
||||
const Fit f = Place(in.nominal, in.cell, in.frames, in.start, s2, true);
|
||||
if (f.ok) {
|
||||
log.Info("ModuleDiag cell_free: beam ({:.3f}, {:.3f}) distance {:.4f} cell {:.3f} {:.3f} {:.3f}",
|
||||
log.Info("ModuleDiag cell_pinned: beam ({:.3f}, {:.3f}) distance {:.4f} cell {:.3f} {:.3f} {:.3f}",
|
||||
f.r.beam_x_px, f.r.beam_y_px, f.r.distance_mm, f.r.cell.a, f.r.cell.b, f.r.cell.c);
|
||||
LogShifts(log, "cell_free", in.centres, f.modules, Moved(f.r, st));
|
||||
LogShifts(log, "cell_pinned", in.centres, f.modules, Moved(f.r, st));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
+6
-2
@@ -975,8 +975,12 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b
|
||||
std::future<saturation_ceiling::Histogram> saturation;
|
||||
if (full && !experiment_.IsRotationIndexing() && !config_.live)
|
||||
saturation = std::async(std::launch::async, [&] { return SampleSaturation(start_image, images_to_process); });
|
||||
if (full && config_.refine_geometry.has_value() && !config_.detector_distance_unknown)
|
||||
RefineStillsGeometry(start_image, end_image, images_to_process, observer);
|
||||
if (full && config_.refine_geometry.has_value() && !config_.detector_distance_unknown) {
|
||||
if (stills_user_cell_ && !config_.detector_distance_given && !config_.live)
|
||||
RefineStillsGeometryOrDeNovo(start_image, end_image, images_to_process, observer);
|
||||
else
|
||||
RefineStillsGeometry(start_image, end_image, images_to_process, observer);
|
||||
}
|
||||
// The detector tilt's fits read only the frames just fitted and the geometry they left: started now, they
|
||||
// run beside the pre-scan and the geometry trial, and are used where neither moved the geometry.
|
||||
if (full)
|
||||
|
||||
@@ -1760,6 +1760,10 @@ class Rugnux {
|
||||
// cannot run. Stills only; rotation has its own two-pass.
|
||||
void RefineStillsGeometry(int start_image, int end_image, int images_to_process,
|
||||
RugnuxObserver *observer);
|
||||
// Stills with a cell given (-C): RefineStillsGeometry with the cell given, and again de novo, each from the
|
||||
// file's geometry, and the run keeps the hypothesis the frames bear out (see the definition).
|
||||
void RefineStillsGeometryOrDeNovo(int start_image, int end_image, int images_to_process,
|
||||
RugnuxObserver *observer);
|
||||
|
||||
// Stills: the beam centre and the detector distance as hypotheses, each indexed on the same sample
|
||||
// of frames and decided on what that sample shows (see DecideStillsGeometry).
|
||||
|
||||
+154
-27
@@ -198,6 +198,18 @@ size_t DropPowderRingSpots(std::vector<DataMessage> &frames, float beam_x, float
|
||||
return dropped;
|
||||
}
|
||||
|
||||
// How far a fit moved the geometry from (distance, cell): the largest of the distance's relative change, each
|
||||
// cell length's and each cell angle's in radians (a degree is 1.7% there, about what it moves a spot by at the
|
||||
// edge of a pattern as a length change of that size does).
|
||||
double StillsGeometryMove(double distance_mm, const UnitCell &cell, const GeometryRefinerResult &r) {
|
||||
constexpr double D2R = PI / 180.0;
|
||||
return std::max({std::abs(r.distance_mm / distance_mm - 1.0),
|
||||
std::abs(r.cell.a / cell.a - 1.0), std::abs(r.cell.b / cell.b - 1.0),
|
||||
std::abs(r.cell.c / cell.c - 1.0),
|
||||
D2R * std::abs(r.cell.alpha - cell.alpha), D2R * std::abs(r.cell.beta - cell.beta),
|
||||
D2R * std::abs(r.cell.gamma - cell.gamma)});
|
||||
}
|
||||
|
||||
} // namespace
|
||||
|
||||
void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_process,
|
||||
@@ -421,17 +433,19 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
GeometryRefinerSettings gr_settings;
|
||||
gr_settings.crystal_system = experiment_.GetCrystalSystem();
|
||||
gr_settings.num_threads = config_.nthreads;
|
||||
// The distance and the cell are refined together, both free. Near the beam they are degenerate - a
|
||||
// spot's radius fixes only their ratio - but not far from it: the radius goes as D tan(2 theta) while
|
||||
// the cell scales sin(theta), so the strong frames' high-angle spots fix both. Holding either one
|
||||
// puts its error on the other and from there on the high-angle predictions: a distance given 1.4%
|
||||
// short, with the cell absorbing it, mispredicted the spots past 2.5 A by enough to make the merged
|
||||
// intensities there negative. A deposited cell is often a percent off this crystal (another
|
||||
// temperature, or the cell the depositors started from), and a stated distance can be off as much,
|
||||
// so the cell (-C) and the file's distance (or DecideStillsGeometry's) are both only where the fit
|
||||
// starts. A distance the user gave (--detector-distance) is held: the cell alone absorbs any error.
|
||||
// The distance and the cell are refined together. Near the beam they are degenerate - a spot's radius
|
||||
// fixes only their ratio - but not far from it: the radius goes as D tan(2 theta) while the cell scales
|
||||
// sin(theta), so the strong frames' high-angle spots fix both. Holding either one puts its error on the
|
||||
// other and from there on the high-angle predictions: a distance given 1.4% short, with the cell
|
||||
// absorbing it, mispredicted the spots past 2.5 A by enough to make the merged intensities there
|
||||
// negative. The cell's shape is no different: with a monoclinic beta held 4 degrees off, the lengths
|
||||
// and the distance took it up (+4% in both) and the merge fell apart. A deposited cell is often a
|
||||
// percent off this crystal (another temperature, or the cell the depositors started from) and a stated
|
||||
// distance can be off as much, so the cell (-C) is a prior of a percent / a degree
|
||||
// (GeometryRefinerSettings::cell_length_sigma) that the spots overrule wherever they determine the
|
||||
// cell, and the file's distance (or DecideStillsGeometry's) is only where the fit starts. A distance
|
||||
// the user gave (--detector-distance) is held: the cell alone absorbs any error.
|
||||
gr_settings.refine_distance = !config_.detector_distance_given;
|
||||
gr_settings.cell_reg_coeff = 0.0;
|
||||
|
||||
const DiffractionGeometry nominal = experiment_.GetDiffractionGeometry();
|
||||
GeometryRefinerResult r = RefineGlobalGeometry(nominal, *cell, strong, gr_settings);
|
||||
@@ -442,13 +456,14 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
}
|
||||
|
||||
// The bundle's spots were indexed at the geometry the fit started from, and the high-angle ones only
|
||||
// where that geometry put them close enough to their reflections: a fit that has moved the distance far
|
||||
// carries the start in which spots it rests on (two starts a percent apart ended 0.2 % apart). So where it moved the distance by more than REFIT_ABOVE, the bundle's frames are indexed
|
||||
// again at the fitted geometry and cell and fitted again, until the distance stays put. Not where
|
||||
// modules are being placed: the placement below is fitted on spots indexed at the start, as before.
|
||||
// where that geometry put them close enough to their reflections: a fit that has moved the distance or the
|
||||
// cell far carries the start in which spots it rests on (two starts a percent apart ended 0.2 % apart). So
|
||||
// where it moved either by more than REFIT_ABOVE (StillsGeometryMove), the bundle's frames are indexed
|
||||
// again at the fitted geometry and cell and fitted again, until both stay put. Not where modules are
|
||||
// being placed: the placement below is fitted on spots indexed at the start, as before.
|
||||
constexpr double REFIT_ABOVE = 0.003, REFIT_CONVERGED = 0.0005;
|
||||
constexpr int MAX_REFITS = 3;
|
||||
double moved_by = std::abs(r.distance_mm / nominal.GetDetectorDistance_mm() - 1.0);
|
||||
double moved_by = StillsGeometryMove(nominal.GetDetectorDistance_mm(), *cell, r);
|
||||
std::vector<GeomRefineFrame> refitted; // the bundle as last indexed again, if it was
|
||||
for (int refit = 0; refit < MAX_REFITS && !placing_modules && moved_by > REFIT_CONVERGED
|
||||
&& (refit > 0 || moved_by > REFIT_ABOVE) && !cancelled_; refit++) {
|
||||
@@ -494,7 +509,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
const GeometryRefinerResult r2 = RefineGlobalGeometry(e.GetDiffractionGeometry(), r.cell, again, gr_settings);
|
||||
if (!r2.ok)
|
||||
break;
|
||||
moved_by = std::abs(r2.distance_mm / r.distance_mm - 1.0);
|
||||
moved_by = StillsGeometryMove(r.distance_mm, r.cell, r2);
|
||||
logger.Info("Geometry refinement: indexed again at the fitted geometry ({} frames), distance {:.4f} -> {:.4f} mm",
|
||||
again.size(), r.distance_mm, r2.distance_mm);
|
||||
r = r2;
|
||||
@@ -517,7 +532,26 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
modules_placed_ = true;
|
||||
if (!centres.empty()) {
|
||||
AssignModules(placing, given);
|
||||
AssignModules(strong, given);
|
||||
|
||||
// The modules are placed from where the frames put the beam, the distance and the cell with the modules
|
||||
// held, not from where the run started: a start off in centre or distance puts its error into the
|
||||
// modules, the scan below first. Started 3% short, with the centre 6 px off, the scan shifted every
|
||||
// module by the centre's error and the placement then moved the two modules of the second ring 13 px
|
||||
// outwards, where the distance's error had left their spots - a pattern of the cell's scale against the
|
||||
// distance that the gauge does not hold, so the distance settled 0.8% long. Every frame is fitted here
|
||||
// with the modules held, and each frame's lattice taken from that fit, so the scan reads the frames at
|
||||
// the geometry it starts from.
|
||||
const GeometryRefinerResult start_fit = RefineGlobalGeometry(nominal, *cell, placing, gr_settings);
|
||||
DiffractionGeometry start_geom = nominal;
|
||||
UnitCell start_cell = *cell;
|
||||
if (start_fit.ok) {
|
||||
start_geom.BeamX_pxl(start_fit.beam_x_px).BeamY_pxl(start_fit.beam_y_px)
|
||||
.DetectorDistance_mm(start_fit.distance_mm);
|
||||
start_cell = start_fit.cell;
|
||||
for (size_t i = 0; i < placing.size(); i++)
|
||||
placing[i].lattice = start_fit.lattices[i];
|
||||
}
|
||||
strong.assign(placing.begin(), placing.begin() + std::min(placing.size(), strong.size()));
|
||||
|
||||
// Where a module is off by more than the spots' spacing allows the refinement to pull, its spots
|
||||
// went unindexed: each module is first searched for over +-SCAN_RANGE_PX on the bundle's frames,
|
||||
@@ -525,7 +559,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
constexpr double SCAN_RANGE_PX = 25.0;
|
||||
TimingMark scan_mark(logger, "module scan");
|
||||
std::vector<ModuleShift> start = centres;
|
||||
const auto scan = ScanModuleShifts(nominal, strong, centres, SCAN_RANGE_PX, config_.nthreads);
|
||||
const auto scan = ScanModuleShifts(start_geom, strong, centres, SCAN_RANGE_PX, config_.nthreads);
|
||||
for (size_t m = 0; m < start.size(); m++)
|
||||
if (scan[m].on_lattice_at_best > scan[m].on_lattice_at_zero) {
|
||||
start[m].dx = scan[m].dx;
|
||||
@@ -560,7 +594,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
st.modules = start;
|
||||
st.refine_modules = true;
|
||||
auto placement = std::async(std::launch::async, [&] {
|
||||
return RefineGlobalGeometry(nominal, *cell, placing, st);
|
||||
return RefineGlobalGeometry(start_geom, start_cell, placing, st);
|
||||
});
|
||||
|
||||
// Held-out residual of the geometry fitted on half `fit`, scored on every spot of the other half.
|
||||
@@ -570,10 +604,10 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
sf.num_threads = std::max(1, config_.nthreads / 8);
|
||||
sf.modules = modules;
|
||||
sf.refine_modules = refine;
|
||||
const auto f = RefineGlobalGeometry(nominal, *cell, half[fit], sf);
|
||||
const auto f = RefineGlobalGeometry(start_geom, start_cell, half[fit], sf);
|
||||
if (!f.ok)
|
||||
return h;
|
||||
DiffractionGeometry g = nominal;
|
||||
DiffractionGeometry g = start_geom;
|
||||
g.BeamX_pxl(f.beam_x_px).BeamY_pxl(f.beam_y_px).DetectorDistance_mm(f.distance_mm);
|
||||
GeometryRefinerSettings ev = gr_settings;
|
||||
ev.num_threads = sf.num_threads;
|
||||
@@ -623,7 +657,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
GeometryRefinerSettings again = st;
|
||||
again.num_threads = config_.nthreads;
|
||||
again.modules_held = hold;
|
||||
rm = RefineGlobalGeometry(nominal, *cell, placing, again);
|
||||
rm = RefineGlobalGeometry(start_geom, start_cell, placing, again);
|
||||
apply = rm.ok;
|
||||
}
|
||||
}
|
||||
@@ -639,7 +673,7 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
// so every spot found and every prediction integrated from here on uses it. The image is
|
||||
// not touched.
|
||||
if (config_.module_diagnostics) {
|
||||
ModuleDiagnosticsInput diag{nominal, *cell, placing, centres, start, st, rm,
|
||||
ModuleDiagnosticsInput diag{start_geom, start_cell, placing, centres, start, st, rm,
|
||||
experiment_.GetDetectorSetup().GetSensorMaterial(),
|
||||
experiment_.GetDetectorSetup().GetSensorThickness_um()};
|
||||
RunModuleDiagnostics(diag);
|
||||
@@ -659,16 +693,110 @@ void Rugnux::RefineStillsGeometry(int start_image, int end_image, int images_to_
|
||||
logger.Info("Geometry refinement: beam ({:.2f}, {:.2f}) -> ({:.2f}, {:.2f}) px, distance {:.4f} -> {:.4f} mm",
|
||||
experiment_.GetBeamX_pxl(), experiment_.GetBeamY_pxl(), r.beam_x_px, r.beam_y_px,
|
||||
experiment_.GetDetectorDistance_mm(), r.distance_mm);
|
||||
logger.Info("Geometry refinement: cell a/b/c {:.3f}/{:.3f}/{:.3f} -> {:.3f}/{:.3f}/{:.3f} A "
|
||||
"(median residual {:.3f} px, {} frames, {} spots)",
|
||||
cell->a, cell->b, cell->c, r.cell.a, r.cell.b, r.cell.c,
|
||||
logger.Info("Geometry refinement: cell {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} {:.2f} -> {:.3f} {:.3f} {:.3f} {:.2f} {:.2f} "
|
||||
"{:.2f} (median residual {:.3f} px, {} frames, {} spots)",
|
||||
cell->a, cell->b, cell->c, cell->alpha, cell->beta, cell->gamma,
|
||||
r.cell.a, r.cell.b, r.cell.c, r.cell.alpha, r.cell.beta, r.cell.gamma,
|
||||
r.median_residual_px, r.frames_used, r.spots_used);
|
||||
// How far the refined cell is from the one the fit was anchored on: in scale (the cube root of the volume),
|
||||
// per length and per angle, beside the distance's move. Whether a cell given (-C) describes the frames at all
|
||||
// is decided apart (RefineStillsGeometryOrDeNovo).
|
||||
const double scale = std::cbrt(CrystalLattice(r.cell).CalcVolume() / CrystalLattice(*cell).CalcVolume()) - 1.0;
|
||||
logger.Info("Geometry refinement: refined cell against the {} cell: scale {:+.2f}%, lengths {:+.2f}/{:+.2f}/{:+.2f}%, "
|
||||
"angles {:+.2f}/{:+.2f}/{:+.2f} deg; distance {:+.2f}%", stills_user_cell_ ? "given" : "de novo",
|
||||
100 * scale, 100 * (r.cell.a / cell->a - 1.0), 100 * (r.cell.b / cell->b - 1.0),
|
||||
100 * (r.cell.c / cell->c - 1.0), r.cell.alpha - cell->alpha, r.cell.beta - cell->beta,
|
||||
r.cell.gamma - cell->gamma, 100 * (r.distance_mm / experiment_.GetDetectorDistance_mm() - 1.0));
|
||||
|
||||
experiment_.BeamX_pxl(r.beam_x_px).BeamY_pxl(r.beam_y_px).DetectorDistance_mm(r.distance_mm);
|
||||
experiment_.SetUnitCell(r.cell);
|
||||
stills_geometry_frames_ = r.frames_used;
|
||||
}
|
||||
|
||||
// The cell given (-C) is a hypothesis, and the cell the frames index to de novo another. A cell of the wrong shape
|
||||
// still indexes the frames, at a distance of its own: given 4% too long with its beta 4 degrees off, a cell
|
||||
// indexed the frames of a crystal at a distance 7% too long, refined towards its own shape there, and merged to
|
||||
// CC1/2 0.4 where the run without -C merged to 0.9. Each hypothesis is fitted from the file's geometry - the
|
||||
// cell given as RefineStillsGeometry fits it, and de novo as it does without -C - and each is then indexed with
|
||||
// its refined cell and geometry on the same sample of frames, as DecideStillsGeometry tries its centres and
|
||||
// distances. The cell given is kept where it puts as many of the sample's spots on the lattice, or nearly
|
||||
// (GIVEN_KEPT_ABOVE, as a file's centre is kept), and also where it does not but its fit puts the detector where
|
||||
// the de novo fit puts it: the detector is the same for every crystal of the run, so a cell that fits at the
|
||||
// same distance is a crystal form of the sample - a minority one, which the user may ask for - and only a cell
|
||||
// that needs a detector of its own describes no crystal of it. Otherwise the run goes on as it would without
|
||||
// -C, from the de novo hypothesis, and says so.
|
||||
void Rugnux::RefineStillsGeometryOrDeNovo(int start_image, int end_image, int images_to_process,
|
||||
RugnuxObserver *observer) {
|
||||
Logger logger("Rugnux");
|
||||
constexpr double GIVEN_KEPT_ABOVE = 0.9;
|
||||
// The fits' distances of one detector, from different cells and frames, agree within about a percent.
|
||||
constexpr double SAME_DETECTOR = 0.02;
|
||||
constexpr int SAMPLE_FRAMES = 96;
|
||||
|
||||
struct State {
|
||||
DiffractionExperiment experiment;
|
||||
std::vector<ModuleShift> module_shifts;
|
||||
bool modules_placed;
|
||||
std::optional<StillsTiltBundle> tilt_bundle;
|
||||
int geometry_frames;
|
||||
};
|
||||
const auto save = [&] {
|
||||
return State{experiment_, module_shifts_, modules_placed_, stills_tilt_bundle_, stills_geometry_frames_};
|
||||
};
|
||||
const auto restore = [&](const State &st) {
|
||||
experiment_ = st.experiment;
|
||||
module_shifts_ = st.module_shifts;
|
||||
modules_placed_ = st.modules_placed;
|
||||
stills_tilt_bundle_ = st.tilt_bundle;
|
||||
stills_geometry_frames_ = st.geometry_frames;
|
||||
};
|
||||
|
||||
const State start = save();
|
||||
const UnitCell given_cell = *stills_user_cell_;
|
||||
RefineStillsGeometry(start_image, end_image, images_to_process, observer);
|
||||
const State given = save();
|
||||
|
||||
restore(start);
|
||||
stills_user_cell_.reset();
|
||||
experiment_.SetUnitCell(std::nullopt);
|
||||
RefineStillsGeometry(start_image, end_image, images_to_process, observer);
|
||||
const State de_novo = save();
|
||||
stills_user_cell_ = given_cell;
|
||||
if (cancelled_ || !de_novo.experiment.GetUnitCell().has_value()) {
|
||||
restore(given);
|
||||
return;
|
||||
}
|
||||
|
||||
const std::vector<int> sample = select_equally_spaced_image_ordinals(
|
||||
images_to_process, std::min(images_to_process, SAMPLE_FRAMES));
|
||||
const std::vector<DataMessage> frames = FindStillsSampleSpots(sample, start_image);
|
||||
IndexingSettings pool_settings = experiment_.GetIndexingSettings();
|
||||
pool_settings.IndexingThreads(StillsIndexerThreads());
|
||||
IndexerThreadPool pool(pool_settings, IndexerConstruction::OnFirstUse, gpu_work_beside_);
|
||||
const auto trials = IndexStillsSample({{given.experiment, &frames}, {de_novo.experiment, &frames}}, pool);
|
||||
|
||||
const UnitCell g = given.experiment.GetUnitCell().value();
|
||||
const UnitCell d = de_novo.experiment.GetUnitCell().value();
|
||||
const double distance_g = given.experiment.GetDetectorDistance_mm();
|
||||
const double distance_d = de_novo.experiment.GetDetectorDistance_mm();
|
||||
logger.Info("Cell hypotheses: the cell given, refined to {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} at {:.2f} mm, puts "
|
||||
"{} sampled spots on the lattice; the cell indexed de novo, {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} {:.2f} at "
|
||||
"{:.2f} mm, {}", g.a, g.b, g.c, g.alpha, g.beta, g.gamma, distance_g, trials[0].on_lattice,
|
||||
d.a, d.b, d.c, d.alpha, d.beta, d.gamma, distance_d, trials[1].on_lattice);
|
||||
const bool same_detector = std::abs(distance_g / distance_d - 1.0) < SAME_DETECTOR;
|
||||
if (trials[0].on_lattice >= GIVEN_KEPT_ABOVE * trials[1].on_lattice || same_detector) {
|
||||
restore(given);
|
||||
return;
|
||||
}
|
||||
restore(de_novo);
|
||||
stills_user_cell_.reset();
|
||||
logger.Warning("The cell given (-C) does not describe these frames: with it the detector is at {:.1f} mm and {} "
|
||||
"sampled spots are on the lattice, with the cell they index to de novo ({:.2f} {:.2f} {:.2f} {:.2f} "
|
||||
"{:.2f} {:.2f}) it is at {:.1f} mm and {} are. The run is processed as without -C",
|
||||
distance_g, trials[0].on_lattice, d.a, d.b, d.c, d.alpha, d.beta, d.gamma, distance_d,
|
||||
trials[1].on_lattice);
|
||||
}
|
||||
|
||||
// The detector tilt of a stills run, as a hypothesis against the tilt the run has (the file's, normally none):
|
||||
// the strong frames of the geometry refinement are dealt alternately into two halves, each half is fitted with
|
||||
// the tilt held and with it refined - the beam, the distance and the cell refined in both, the modules held
|
||||
@@ -727,7 +855,6 @@ GeometryRefinerSettings Rugnux::StillsTiltSettings() const {
|
||||
st.crystal_system = experiment_.GetCrystalSystem();
|
||||
st.num_threads = config_.nthreads;
|
||||
st.refine_distance = !config_.detector_distance_given;
|
||||
st.cell_reg_coeff = 0.0;
|
||||
return st;
|
||||
}
|
||||
|
||||
|
||||
@@ -217,3 +217,38 @@ TEST_CASE("RefineGlobalGeometry_Tilt") {
|
||||
CHECK(tilted.capped_rms_px < held.capped_rms_px);
|
||||
}
|
||||
}
|
||||
|
||||
// A monoclinic crystal whose cell is given with the wrong beta and lengths a percent off: the angle is refined
|
||||
// with the lengths, so the distance is not left to absorb the shape - with beta held, the lengths and the distance
|
||||
// take it up instead.
|
||||
TEST_CASE("RefineGlobalGeometry_MonoclinicBeta") {
|
||||
DiffractionExperiment exp;
|
||||
exp.IncidentEnergy_keV(WVL_1A_IN_KEV).BeamX_pxl(1000).BeamY_pxl(1000).DetectorDistance_mm(150);
|
||||
const DiffractionGeometry truth = exp.GetDiffractionGeometry();
|
||||
|
||||
ModuleShift whole;
|
||||
whole.x0 = 1000;
|
||||
whole.y0 = 1000;
|
||||
const CrystalLattice latt(50, 60, 70, 90, 95, 90);
|
||||
std::vector<GeomRefineFrame> frames = FramesOnMovedModules(exp, latt, {whole, whole, whole, whole});
|
||||
REQUIRE(frames.size() > 30);
|
||||
for (auto &f: frames)
|
||||
for (auto &s: f.spots)
|
||||
s.module = -1;
|
||||
|
||||
DiffractionGeometry nominal = truth;
|
||||
nominal.DetectorDistance_mm(152);
|
||||
const UnitCell given{50.5, 59.5, 70.5, 90, 96, 90};
|
||||
GeometryRefinerSettings st;
|
||||
st.crystal_system = gemmi::CrystalSystem::Monoclinic;
|
||||
st.num_threads = 2;
|
||||
|
||||
const auto r = RefineGlobalGeometry(nominal, given, frames, st);
|
||||
REQUIRE(r.ok);
|
||||
CHECK(r.cell.beta == Catch::Approx(95).margin(0.01));
|
||||
CHECK(r.cell.a == Catch::Approx(50).margin(0.02));
|
||||
CHECK(r.cell.b == Catch::Approx(60).margin(0.02));
|
||||
CHECK(r.cell.c == Catch::Approx(70).margin(0.02));
|
||||
CHECK(r.distance_mm == Catch::Approx(150).margin(0.05));
|
||||
CHECK(r.median_residual_px < 0.1);
|
||||
}
|
||||
Reference in new issue
Block a user