leonarski_fandClaude Opus 5
212fbf9bab
The stills first pass drew frames from a shared cursor and stopped when a shared counter reached its target, which got two things wrong at once. The cursor walked the equally spaced sample in ascending order, so stopping early read only its leading PREFIX - the beam centre, distance and cell were fitted to the beginning of the run, not across it, and the comment claiming otherwise was wrong. And where the stop landed depended on how the workers happened to interleave, so the set of frames varied run to run: on the same data at -N 32 and -N 8 the pass examined 483 and 457 frames and refined the detector distance to 168.0481 and 168.0530 mm. The sample is now cut into a fixed number of interleaved stripes, each stopping once it has contributed its share. Every stripe spans the whole run, so an early stop no longer biases the fit, and a stripe is processed identically whichever worker claims it - so what gets examined depends only on the data, not on timing and not on -N. The same three runs now give 451 frames examined and 168.0452 mm, identically. The bundle selection was order-dependent too: frames are collected in worker-completion order and sorted by spot count with a non-stable sort, so equally strong frames swapped places between runs. They carry their image ordinal now and it breaks the tie. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
62 lines
3.2 KiB
C++
62 lines
3.2 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#pragma once
|
|
|
|
#include <cstdint>
|
|
#include <vector>
|
|
|
|
#include "../../common/CrystalLattice.h"
|
|
#include "../../common/DiffractionGeometry.h"
|
|
#include "../../common/UnitCell.h"
|
|
#include "gemmi/symmetry.hpp"
|
|
|
|
// 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
|
|
// 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. Detector tilt (rot1/rot2) is NOT refined (gauge-coupled with the beam; zero indexing gain).
|
|
|
|
struct GeomRefineSpot {
|
|
float x = 0, y = 0; // observed (raw) spot position, pixels
|
|
int32_t h = 0, k = 0, l = 0; // integer Miller index assigned in pass 1
|
|
};
|
|
|
|
struct GeomRefineFrame {
|
|
CrystalLattice lattice; // per-frame indexed lattice (real space), for the orientation seed
|
|
std::vector<GeomRefineSpot> spots; // indexed spots on this frame
|
|
int32_t ordinal = 0; // image this came from; ties are broken on it, so the
|
|
// selected bundle does not depend on collection order
|
|
};
|
|
|
|
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;
|
|
int rounds = 3; // hkl-reassignment / tolerance-tightening rounds
|
|
double beam_range_px = 8.0; // beam bound around the nominal value
|
|
double distance_range_mm = 8.0; // distance bound around the nominal value
|
|
int min_spots_total = 100; // refuse to refine below this many indexed spots
|
|
int num_threads = 1;
|
|
};
|
|
|
|
struct GeometryRefinerResult {
|
|
bool ok = false;
|
|
double beam_x_px = 0, beam_y_px = 0; // refined beam centre (== direct beam, tilt not refined)
|
|
double distance_mm = 0;
|
|
UnitCell cell{}; // refined cell
|
|
int frames_used = 0;
|
|
int spots_used = 0;
|
|
double median_residual_px = 0; // final per-spot detector residual (fit quality)
|
|
};
|
|
|
|
// 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.
|
|
GeometryRefinerResult RefineGlobalGeometry(const DiffractionGeometry &nominal_geom,
|
|
const UnitCell &input_cell,
|
|
const std::vector<GeomRefineFrame> &frames,
|
|
const GeometryRefinerSettings &settings);
|