leonarski_fandClaude Opus 5
ee0ad8d149
AnalyzeGridScan takes a finished ScanResult plus its GridScanSettings and returns the list of crystals the raster hit, sorted by score so element 0 is the one to collect. Pure function - no I/O, no FPGA, no JSON. Today the list holds nought or one entry; N is the point of the shape. The per-image protein score is scattered back onto the display grid through Rearrange, which already knows the snake order, the vertical flag and the step signs, thresholded, and labelled into blobs. Each blob is then measured: - Centre is a weighted centroid, pulled towards the cells that diffract best. The pull uses the RANK of the resolution inside the blob, never its value, so a salt grain reporting an absurd 0.8 A weighs exactly what a genuine best cell weighs and cannot drag the centre however extreme its number. A cell with no resolution gets the lowest weight rather than being dropped. The centroid of a banana- or L-shaped blob can land outside the blob, where no image exists, so image_number is snapped to the nearest cell that was actually collected. - Second moments are taken in micrometres, not in cells. A 20 x 16 um raster is ordinary and moments in cell units give a wrong angle - eight degrees wrong on the staircase in the tests. The angle is an axis, so it lives in [0,180) and wraps there. - The axis DIRECTION comes from the eigenvector but the LENGTH from the projected extent, because "how far do I scan" is an extent question and the constant taking a second moment to a length assumes a shape a blob of five cells does not have. Where the two disagree about which axis is longer - a moment dominated by clumps at the ends - the extents are swapped and the angle turned a quarter turn, so major_um >= minor_um with angle_deg along it is an invariant a consumer can draw a frame from. - score is the MEAN protein score over the blob, not the peak: the score saturates, so the peak is 1.0 for every real crystal and ranks nothing. res_A is the 25th percentile, not the minimum, the minimum being precisely where a salt spot or a hot pixel shows up; it is NaN when nothing in the blob measured a resolution. Sizes are measured and the beam is left in them. The beam is already in the file as incident_beam_size, so a consumer can deconvolve reproducibly and reversibly instead of inheriting ours; the result carries the beam size so it says what the extents contain. The header records that removing an anisotropic beam is a covariance-matrix subtraction followed by re-diagonalisation, not a per-axis quadrature removal, which is silently wrong whenever the crystal is not aligned with the grid - the needle case this design exists for. Labelling is a small dense flood fill in common/, 8-connected. StrongPixelSet::sparseccl is the wrong abstraction for a dense grid map: sparse union-find over raster-ordered strong pixels, hardcoded module dimensions, a 4000-pixel cap, spot-shape acceptance, and an FPGA header. 8-connected rather than 4 because where the step is coarser than the beam a needle at 45 degrees lands as corner-touching cells; under 4-connectivity that breaks into single cells and the minimum-size rule then discards the crystal entirely, which is the case oriented axes exist to catch. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
54 lines
2.0 KiB
C++
54 lines
2.0 KiB
C++
// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include "ConnectedComponents.h"
|
|
|
|
std::vector<int32_t> LabelConnectedComponents(const std::vector<uint8_t> &map,
|
|
int64_t nx, int64_t ny, int64_t min_size) {
|
|
std::vector<int32_t> label(map.size(), 0);
|
|
std::vector<int64_t> stack; // cells of the component being grown
|
|
std::vector<int64_t> members; // same, kept so the component can be erased if it is too small
|
|
int32_t next_label = 0;
|
|
|
|
for (int64_t seed = 0; seed < nx * ny; seed++) {
|
|
if (!map[seed] || label[seed] != 0)
|
|
continue;
|
|
|
|
next_label++;
|
|
stack = {seed};
|
|
members.clear();
|
|
label[seed] = next_label;
|
|
|
|
while (!stack.empty()) {
|
|
int64_t i = stack.back();
|
|
stack.pop_back();
|
|
members.push_back(i);
|
|
|
|
const int64_t x = i % nx;
|
|
const int64_t y = i / nx;
|
|
const bool left = x > 0, right = x < nx - 1, up = y > 0, down = y < ny - 1;
|
|
const int64_t neighbour[8] = {left ? i - 1 : -1,
|
|
right ? i + 1 : -1,
|
|
up ? i - nx : -1,
|
|
down ? i + nx : -1,
|
|
left && up ? i - nx - 1 : -1,
|
|
right && up ? i - nx + 1 : -1,
|
|
left && down ? i + nx - 1 : -1,
|
|
right && down ? i + nx + 1 : -1};
|
|
for (int64_t j : neighbour) {
|
|
if (j >= 0 && map[j] && label[j] == 0) {
|
|
label[j] = next_label;
|
|
stack.push_back(j);
|
|
}
|
|
}
|
|
}
|
|
|
|
if (static_cast<int64_t>(members.size()) < min_size) {
|
|
for (int64_t i : members)
|
|
label[i] = 0;
|
|
}
|
|
}
|
|
|
|
return label;
|
|
}
|