leonarski_fandClaude Opus 5
28a98e1a14
A rotation sweep never records a double cone of half-angle asin(lambda/2d) about the spindle. That loss is normally repaired by the point group; it is not repaired when an operator's axis lies inside the cone, because the cone then maps onto itself. A still cannot know the point group, but it can see where the crystal's short lattice rows are, and a symmetry axis is always one of them - measured over 107 solved cells its length is 1.4x the shortest row at the median. The number is the fraction of the blind cone that a 2-fold about the nearest short row carries back into the blind cone: two equal caps of angular radius theta_max whose centres are 2*beta apart, which is the circle-lens area in beta/theta_max and within 0.035 of the spherical value even for a 55 deg cone. 0 means one sweep reaches everything symmetry could give; 1 means the row is on the spindle and the whole cone is lost coherently. No goniometer geometry enters, so the number describes the problem and leaves the remedy to the beamline. Free: it rides on the FilterFFTResults shortlist the indexer already builds, needing only the spindle, the wavelength and the frame's own resolution. FilterFFTResults gains an optional out-parameter for each row's peak prominence, which the length window is taken over. Measured on 384 stills of 22 solved crystals against the frames' own symmetry axes, with the same procedure re-run along five decoy directions per frame as the null: 0.93 of severe orientations reported severe, 0.010 of harmless ones reported severe, AUC 0.948, and the value tracks the true severity to 0.10 at the 90th percentile. Below 60 spots the misses triple, so that is where it stops answering. It is blind to a symmetry axis much longer than the crystal's shortest row: measured on synthetic stills the search grid finds a 150 A row every time, a 200 A row half the time and a 300 A row once in sixteen. Two cheap consistency tests refuse to answer for part of that regime and cost nothing on real frames, but they do not cover it - the reach is a documented property of the number, not something every frame can detect. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
88 lines
3.8 KiB
C++
88 lines
3.8 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <algorithm>
|
|
#include <cmath>
|
|
|
|
#include "SpindleBlindFraction.h"
|
|
|
|
namespace {
|
|
// Only rows short enough to be a plausible symmetry axis count. The cut is relative to the
|
|
// crystal's own shortest row, not an absolute length, so it works the same for a 40 A cell and
|
|
// a 200 A one: measured over 107 solved cells, 2.5 x the shortest row covers 88% of all
|
|
// symmetry axes and 98% of crystals' shortest one.
|
|
constexpr float MAX_ROW_LENGTH_RATIO = 2.5f;
|
|
// The reference length is taken over the strong rows only. A long-cell still makes the pass
|
|
// invent short spurious rows - measured at 0.4-0.6 of the true row's peak - and taking the
|
|
// reference over every row would let one of those shrink the window until the real aligned row
|
|
// fell outside it, turning a severe orientation into a reported zero.
|
|
constexpr float MIN_MAGNITUDE_RATIO = 0.5f;
|
|
// Where the search grid can no longer resolve the crystal's rows, the pass stops returning them
|
|
// and starts returning short spurious ones instead, and the shortlist becomes internally
|
|
// inconsistent: its strong rows are many times longer than its shortest entry. Measured over 384
|
|
// frames of 22 solved crystals the ratio never exceeded 2.34 and was 1.00 at the median; on a
|
|
// synthetic still whose cell is past the grid's reach it runs 4-13. Past this the score would be
|
|
// reporting a cone it cannot see into, so it reports nothing instead.
|
|
constexpr float MAX_REFERENCE_LENGTH_RATIO = 3.0f;
|
|
}
|
|
|
|
float BlindConeSelfOverlap(float x) {
|
|
if (x >= 1.0f)
|
|
return 0.0f;
|
|
if (x <= 0.0f)
|
|
return 1.0f;
|
|
return static_cast<float>(2.0 / M_PI) *
|
|
(std::acos(x) - x * std::sqrt(1.0f - x * x));
|
|
}
|
|
|
|
std::optional<SpindleSeverity> SpindleBlindFraction(const std::vector<Coord> &rows,
|
|
const std::vector<float> &magnitudes,
|
|
const Coord &spindle,
|
|
float theta_max_deg) {
|
|
if (rows.empty() || rows.size() != magnitudes.size() || theta_max_deg <= 0)
|
|
return {};
|
|
|
|
const float axis_length = spindle.Length();
|
|
if (axis_length < 1e-6f)
|
|
return {};
|
|
const Coord axis = spindle.Normalize();
|
|
|
|
float max_magnitude = 0;
|
|
for (const auto &m : magnitudes)
|
|
max_magnitude = std::max(max_magnitude, m);
|
|
|
|
float reference_length = 0;
|
|
for (size_t i = 0; i < rows.size(); i++)
|
|
if (magnitudes[i] >= MIN_MAGNITUDE_RATIO * max_magnitude) {
|
|
const float l = rows[i].Length();
|
|
if (reference_length == 0 || l < reference_length)
|
|
reference_length = l;
|
|
}
|
|
if (reference_length == 0)
|
|
return {};
|
|
|
|
float shortest_length = reference_length;
|
|
for (const auto &r : rows)
|
|
shortest_length = std::min(shortest_length, r.Length());
|
|
if (reference_length > MAX_REFERENCE_LENGTH_RATIO * shortest_length)
|
|
return {};
|
|
|
|
SpindleSeverity ret;
|
|
for (size_t i = 0; i < rows.size(); i++) {
|
|
const float length = rows[i].Length();
|
|
if (length <= 0 || length > MAX_ROW_LENGTH_RATIO * reference_length)
|
|
continue;
|
|
if (magnitudes[i] < MIN_MAGNITUDE_RATIO * max_magnitude)
|
|
continue;
|
|
const float cos_beta = std::min(1.0f, std::fabs(rows[i] * axis) / length);
|
|
const float beta_deg = static_cast<float>(std::acos(cos_beta) * 180.0 / M_PI);
|
|
const float score = BlindConeSelfOverlap(beta_deg / theta_max_deg);
|
|
if (score > ret.score || ret.row_length_A == 0) {
|
|
ret.score = score;
|
|
ret.row_length_A = length;
|
|
ret.miss_angle_deg = beta_deg;
|
|
}
|
|
}
|
|
return ret;
|
|
}
|