From ee0ad8d149f35059fc4e7e26b2c6d3757355ed17 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 7 Sep 2026 23:39:16 +0200 Subject: [PATCH] grid scan: a completed raster is reduced to the crystals in it, oriented and ranked 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 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- common/CMakeLists.txt | 2 + common/ConnectedComponents.cpp | 53 +++ common/ConnectedComponents.h | 19 + image_analysis/CMakeLists.txt | 3 +- .../grid_scan_analysis/AnalyzeGridScan.cpp | 206 +++++++++++ .../grid_scan_analysis/AnalyzeGridScan.h | 28 ++ .../grid_scan_analysis/CMakeLists.txt | 2 + .../grid_scan_analysis/GridScanResult.h | 39 +++ tests/AnalyzeGridScanTest.cpp | 330 ++++++++++++++++++ tests/CMakeLists.txt | 1 + 10 files changed, 682 insertions(+), 1 deletion(-) create mode 100644 common/ConnectedComponents.cpp create mode 100644 common/ConnectedComponents.h create mode 100644 image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp create mode 100644 image_analysis/grid_scan_analysis/AnalyzeGridScan.h create mode 100644 image_analysis/grid_scan_analysis/CMakeLists.txt create mode 100644 image_analysis/grid_scan_analysis/GridScanResult.h create mode 100644 tests/AnalyzeGridScanTest.cpp diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index 59aaf46c3..7cc3ac5c9 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -117,6 +117,8 @@ ADD_LIBRARY(JFJochCommon STATIC CrystalLattice.h ScanResult.cpp ScanResult.h + ConnectedComponents.cpp + ConnectedComponents.h ScanResultGenerator.cpp ScanResultGenerator.h BraggIntegrationSettings.cpp diff --git a/common/ConnectedComponents.cpp b/common/ConnectedComponents.cpp new file mode 100644 index 000000000..210c52f8f --- /dev/null +++ b/common/ConnectedComponents.cpp @@ -0,0 +1,53 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ConnectedComponents.h" + +std::vector LabelConnectedComponents(const std::vector &map, + int64_t nx, int64_t ny, int64_t min_size) { + std::vector label(map.size(), 0); + std::vector stack; // cells of the component being grown + std::vector 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(members.size()) < min_size) { + for (int64_t i : members) + label[i] = 0; + } + } + + return label; +} diff --git a/common/ConnectedComponents.h b/common/ConnectedComponents.h new file mode 100644 index 000000000..e7c5e06de --- /dev/null +++ b/common/ConnectedComponents.h @@ -0,0 +1,19 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +// 8-connected labelling of a dense boolean map of nx by ny cells, stored row-major. +// Returns one label per cell: 0 for a cell that is not set and for a component smaller +// than min_size, 1..n for the components that survive. Labels are not renumbered after +// the size cut, so the returned labels are not necessarily contiguous. +// +// Corner-touching cells are one component. On a grid scan that matters: where the step is +// coarser than the beam, a needle lying diagonally is sampled as a chain of cells that meet +// only at their corners, and 4-connectivity would break it into single cells that every +// minimum-size cut then throws away. +std::vector LabelConnectedComponents(const std::vector &map, + int64_t nx, int64_t ny, int64_t min_size); diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 2bf22bd3a..d5b28d225 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -87,5 +87,6 @@ ADD_SUBDIRECTORY(scale_merge) ADD_SUBDIRECTORY(image_preprocessing) ADD_SUBDIRECTORY(azint) ADD_SUBDIRECTORY(roi) +ADD_SUBDIRECTORY(grid_scan_analysis) -TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochROIIntegration JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi) +TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochROIIntegration JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge JFJochGridScanAnalysis gemmi) diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp new file mode 100644 index 000000000..5bce4d06b --- /dev/null +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp @@ -0,0 +1,206 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#include "../../common/ConnectedComponents.h" +#include "../../common/JFJochMath.h" +#include "AnalyzeGridScan.h" + +namespace { + // A cell counts as protein above this. The per-image score saturates, so this only has to + // separate "something diffracted here" from "nothing did". + constexpr float PROTEIN_SCORE_THRESHOLD = 0.5f; + + // Two cells can be the two ends of a single hit lying on a cell boundary; three is the + // smallest patch that is a shape rather than a coincidence. + constexpr int64_t MIN_BLOB_CELLS = 3; + + constexpr float NO_VALUE = -1.0f; +} + +GridScanResult AnalyzeGridScan(const ScanResult &scan, + const GridScanSettings &grid, + float beam_size_x_um, + float beam_size_y_um) { + const int64_t nx = grid.GetGridSizeX_step(); + const int64_t ny = grid.GetGridSizeY_step(); + const float step_x = fabsf(grid.GetGridStepX_um()); + const float step_y = fabsf(grid.GetGridStepY_um()); + + // Scatter the per-image quantities onto the display grid. Rearrange knows about snake order, + // the vertical flag and the step signs, so nothing here has to. + std::vector protein(nx * ny, NO_VALUE); + std::vector ice(nx * ny, NO_VALUE); + std::vector res(nx * ny, NO_VALUE); + std::vector image_no(nx * ny, -1); + + for (const auto &elem: scan.images) { + if (elem.number < 0 || elem.number >= grid.GetNElem()) + continue; + const int64_t i = grid.Rearrange(elem.number); + protein[i] = elem.protein_score.value_or(NO_VALUE); + ice[i] = elem.ice_score.value_or(NO_VALUE); + res[i] = elem.res.value_or(NO_VALUE); + image_no[i] = elem.number; + } + + std::vector above(nx * ny); + for (int64_t i = 0; i < nx * ny; i++) + above[i] = (protein[i] > PROTEIN_SCORE_THRESHOLD) ? 1 : 0; + + const std::vector label = LabelConnectedComponents(above, nx, ny, MIN_BLOB_CELLS); + const int32_t n_label = label.empty() ? 0 : *std::max_element(label.begin(), label.end()); + + GridScanResult result; + result.beam_size_x_um = beam_size_x_um; + result.beam_size_y_um = beam_size_y_um; + + for (int32_t l = 1; l <= n_label; l++) { + std::vector cell; + for (int64_t i = 0; i < nx * ny; i++) { + if (label[i] == l) + cell.push_back(i); + } + if (cell.empty()) + continue; // label fell below MIN_BLOB_CELLS and was erased + + const auto n = static_cast(cell.size()); + + // Cell centres in MICROMETRES, taken once. Every geometric quantity below is computed from + // these and never from cell units: step_x and step_y genuinely differ (20 x 16 um is an + // ordinary raster), so a second moment taken in cells gives the wrong axis angle. + // The i / nx and i % nx are the row and column of a row-major index - an index split, and + // truncation is the whole point of it; there is no precision to lose there. + std::vector px(cell.size()), py(cell.size()); + for (size_t k = 0; k < cell.size(); k++) { + px[k] = static_cast(cell[k] % nx) * step_x; + py[k] = static_cast(cell[k] / nx) * step_y; + } + + // Centre, pulled towards the cells that diffract best. The pull uses the RANK of the + // resolution within the blob, never its value: the best cell gets p = 1, the worst p = 0, + // and a cell with no resolution at all gets p = 0 rather than being dropped. A salt grain + // reporting an absurd 0.8 A is then weighted exactly like a genuine best cell, so no + // artefact can drag the centre however extreme its number is. + std::vector by_res; + for (int64_t i: cell) { + if (res[i] > 0) + by_res.push_back(i); + } + std::sort(by_res.begin(), by_res.end(), [&](int64_t a, int64_t b) { return res[a] < res[b]; }); + + std::vector weight(cell.size(), 1.0f); + for (size_t r = 0; r < by_res.size(); r++) { + const float p = (by_res.size() == 1) ? 1.0f + : 1.0f - static_cast(r) / static_cast(by_res.size() - 1); + // cell is built in ascending grid order, so it can be searched directly + weight[std::lower_bound(cell.begin(), cell.end(), by_res[r]) - cell.begin()] = 1.0f + 0.5f * p; + } + + float sum_w = 0, sum_wx = 0, sum_wy = 0; + for (size_t k = 0; k < cell.size(); k++) { + sum_w += weight[k]; + sum_wx += weight[k] * px[k]; + sum_wy += weight[k] * py[k]; + } + const float cx_um = sum_wx / sum_w; + const float cy_um = sum_wy / sum_w; + + // A weighted centroid of a banana- or L-shaped blob can land outside the blob, where no + // image was ever collected. The image number has to name a cell that exists, so snap. + size_t nearest = 0; + float nearest_d2 = INFINITY; + for (size_t k = 0; k < cell.size(); k++) { + const float d2 = (px[k] - cx_um) * (px[k] - cx_um) + (py[k] - cy_um) * (py[k] - cy_um); + if (d2 < nearest_d2) { + nearest_d2 = d2; + nearest = k; + } + } + + float sxx = 0, syy = 0, sxy = 0; + for (size_t k = 0; k < cell.size(); k++) { + const float dx = px[k] - cx_um; + const float dy = py[k] - cy_um; + sxx += dx * dx; + syy += dy * dy; + sxy += dx * dy; + } + const float angle = 0.5f * atan2f(2 * sxy, sxx - syy); + const float cos_a = cosf(angle); + const float sin_a = sinf(angle); + + // Direction from the eigenvector, LENGTH from the projected extent. "How far do I scan" + // is an extent question, and the constant taking a second moment to a length depends on + // an assumed shape that a blob of a few cells does not have. + float min_u = INFINITY, max_u = -INFINITY, min_v = INFINITY, max_v = -INFINITY; + for (size_t k = 0; k < cell.size(); k++) { + const float dx = px[k] - cx_um; + const float dy = py[k] - cy_um; + min_u = std::min(min_u, dx * cos_a + dy * sin_a); + max_u = std::max(max_u, dx * cos_a + dy * sin_a); + min_v = std::min(min_v, -dx * sin_a + dy * cos_a); + max_v = std::max(max_v, -dx * sin_a + dy * cos_a); + } + // The span runs between cell centres, so one cell has to be added back. A cell is a + // step_x by step_y rectangle, and its own width along a direction is that rectangle's + // support width - which is the plain step only when the axis lies along the grid. + const float cell_along_u = fabsf(step_x * cos_a) + fabsf(step_y * sin_a); + const float cell_along_v = fabsf(step_x * sin_a) + fabsf(step_y * cos_a); + + float sum_protein = 0, sum_ice = 0; + for (int64_t i: cell) { + sum_protein += protein[i]; + sum_ice += std::max(ice[i], 0.0f); + } + + GridScanCrystal crystal; + crystal.nx = cx_um / step_x; + crystal.ny = cy_um / step_y; + crystal.x_um = cx_um; + crystal.y_um = cy_um; + crystal.image_number = image_no[cell[nearest]]; + crystal.major_um = max_u - min_u + cell_along_u; + crystal.minor_um = max_v - min_v + cell_along_v; + crystal.angle_deg = angle * 180.0f / static_cast(PI); + + // The angle is the axis of the larger second MOMENT, the extents are MEASURED spans, and + // for a strongly non-convex blob the two can disagree about which axis is the longer. + // A consumer draws a major by minor frame rotated by angle_deg, so keep both facts by + // turning the frame a quarter turn rather than by dropping one of them. + if (crystal.major_um < crystal.minor_um) { + std::swap(crystal.major_um, crystal.minor_um); + crystal.angle_deg += 90.0f; + } + + // atan2 returns (-pi,pi], so the half-angle is in (-pi/2,pi/2] and the quarter turn above + // can carry it past 180; an axis has no sign, so fold it into [0,180). + if (crystal.angle_deg < 0) + crystal.angle_deg += 180.0f; + if (crystal.angle_deg >= 180.0f) + crystal.angle_deg -= 180.0f; + + // The MEAN protein score, not the peak: the score saturates, so the peak is 1.0 for every + // real crystal and ranks nothing. The mean stays a detection confidence and compares. + crystal.score = sum_protein / n; + crystal.ice_score = sum_ice / n; + + // The 25th percentile, not the minimum: the single best cell in a blob is precisely where + // a salt spot or a hot pixel shows up. + if (!by_res.empty()) { + const auto q = static_cast(0.25 * static_cast(by_res.size() - 1) + 0.5); + crystal.res_A = res[by_res[q]]; + } + crystal.n_images = static_cast(cell.size()); + + result.crystals.push_back(crystal); + } + + std::sort(result.crystals.begin(), result.crystals.end(), + [](const GridScanCrystal &a, const GridScanCrystal &b) { return a.score > b.score; }); + + return result; +} diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.h b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h new file mode 100644 index 000000000..6c52acf39 --- /dev/null +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h @@ -0,0 +1,28 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "../../common/GridScanSettings.h" +#include "../../common/ScanResult.h" +#include "GridScanResult.h" + +// Finds the crystals in a completed grid scan: the per-image protein score is scattered back onto +// the raster, the map is thresholded and labelled, and each blob is reported as one crystal. +// Pure - no I/O, no state, no detector. +// +// The sizes reported are MEASURED, and the beam is still in them: what the raster sees is the +// crystal convolved with the beam, and nothing here takes the beam back out. That is deliberate - +// the beam size is in the file as incident_beam_size, so a downstream consumer can do the +// deconvolution itself, reproducibly and reversibly, and is not stuck with ours. +// +// If it does: removing an anisotropic beam is a subtraction of the two 2x2 covariance matrices +// followed by re-diagonalising the difference, NOT a per-axis quadrature removal of beam_x from +// major_um and beam_y from minor_um. Per-axis is silently wrong the moment the crystal is not +// aligned with the grid axes - a needle at 45 deg has both beam widths mixed into both of its own +// axes - and a needle at an arbitrary angle is exactly the case this whole design exists for. +// beam_size_x_um/beam_size_y_um are only copied into the result, so it says what the sizes contain. +GridScanResult AnalyzeGridScan(const ScanResult &scan, + const GridScanSettings &grid, + float beam_size_x_um, + float beam_size_y_um); diff --git a/image_analysis/grid_scan_analysis/CMakeLists.txt b/image_analysis/grid_scan_analysis/CMakeLists.txt new file mode 100644 index 000000000..edc3ae5c5 --- /dev/null +++ b/image_analysis/grid_scan_analysis/CMakeLists.txt @@ -0,0 +1,2 @@ +ADD_LIBRARY(JFJochGridScanAnalysis STATIC AnalyzeGridScan.cpp AnalyzeGridScan.h GridScanResult.h) +TARGET_LINK_LIBRARIES(JFJochGridScanAnalysis JFJochCommon) diff --git a/image_analysis/grid_scan_analysis/GridScanResult.h b/image_analysis/grid_scan_analysis/GridScanResult.h new file mode 100644 index 000000000..770152659 --- /dev/null +++ b/image_analysis/grid_scan_analysis/GridScanResult.h @@ -0,0 +1,39 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +// One crystal found in a grid scan. Positions are in the display grid of GridScanSettings - +// column 0 is the lowest x, row 0 the lowest y, whatever direction the stage actually moved in - +// so they match the per-image positions the scan writes with GetXContainer_m/GetYContainer_m. +struct GridScanCrystal { + float nx = 0, ny = 0; // centre, grid coords, fractional, 0-based + float x_um = 0, y_um = 0; // centre, signed offset from centre of cell (0,0), along grid axes + int64_t image_number = -1; // nearest COLLECTED image, for the DAQ to address + // Extent along the crystal's own principal axes. major_um >= minor_um always, and + // angle_deg points along major_um - so a consumer can draw a major_um by minor_um frame + // rotated by angle_deg without checking which of the two is the longer. + float major_um = 0, minor_um = 0; + // Major axis from the +x grid axis, counter-clockwise. This is an AXIS, not a direction, so it + // lives in [0,180) and wraps there: 179 deg is adjacent to 0 deg, and code comparing two angles + // has to fold the difference into [0,90]. On a round blob the axis is arbitrary and the value is + // whatever the numerics produced - major_um/minor_um near 1 is what says so. + float angle_deg = 0; + float score = 0, ice_score = 0; // 0-1 + float res_A = NAN; // robust best resolution in the blob, NaN if none was measured + int64_t n_images = 0; +}; + +// Crystals found in one completed grid scan, sorted by score descending: element 0 is the one to +// collect. Empty when the raster hit nothing. +struct GridScanResult { + std::vector crystals; + + // The beam the sizes above were measured with, along the grid axes. The crystal extents still + // contain it (see AnalyzeGridScan), so this says what a consumer has to take back out. + float beam_size_x_um = 0, beam_size_y_um = 0; +}; diff --git a/tests/AnalyzeGridScanTest.cpp b/tests/AnalyzeGridScanTest.cpp new file mode 100644 index 000000000..279be156c --- /dev/null +++ b/tests/AnalyzeGridScanTest.cpp @@ -0,0 +1,330 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include + +#include "../common/GridScanSettings.h" +#include "../common/ScanResult.h" +#include "../image_analysis/grid_scan_analysis/AnalyzeGridScan.h" + +namespace { + // A plain raster with positive steps and no snake, so image number == grid index and a test + // can name cells by (column, row) directly. + GridScanSettings MakeGrid(int64_t nx, int64_t ny, float step_x_um, float step_y_um) { + GridScanSettings grid(nx, step_x_um, step_y_um, false, false); + grid.ImageNum(nx * ny); + return grid; + } + + // Every cell of the raster is collected; the ones named in hits get a protein score, the + // rest get nothing. + ScanResult MakeScan(int64_t nx, int64_t ny, + const std::vector> &hits, float score) { + ScanResult scan; + for (int64_t i = 0; i < nx * ny; i++) { + ScanResultElem elem; + elem.number = i; + scan.images.push_back(elem); + } + for (const auto &[x, y]: hits) + scan.images[y * nx + x].protein_score = score; + return scan; + } +} + +TEST_CASE("AnalyzeGridScan finds nothing in an empty scan", "[AnalyzeGridScan]") { + auto grid = MakeGrid(6, 6, 10.0f, 10.0f); + + SECTION("no image scored at all") { + auto result = AnalyzeGridScan(MakeScan(6, 6, {}, 0.0f), grid, 5.0f, 5.0f); + CHECK(result.crystals.empty()); + } + + SECTION("every image scored, none above threshold") { + auto scan = MakeScan(6, 6, {}, 0.0f); + for (auto &elem: scan.images) + elem.protein_score = 0.1f; + auto result = AnalyzeGridScan(scan, grid, 5.0f, 5.0f); + CHECK(result.crystals.empty()); + } + + SECTION("a two-cell patch is below the minimum blob size") { + auto scan = MakeScan(6, 6, {{2, 2}, {3, 2}}, 0.9f); + auto result = AnalyzeGridScan(scan, grid, 5.0f, 5.0f); + CHECK(result.crystals.empty()); + } +} + +TEST_CASE("AnalyzeGridScan measures a round blob", "[AnalyzeGridScan]") { +// A plus sign centred on cell (3,3): symmetric, so the centre is the centre cell exactly. + auto grid = MakeGrid(7, 7, 10.0f, 10.0f); + auto scan = MakeScan(7, 7, {{3, 2}, {2, 3}, {3, 3}, {4, 3}, {3, 4}}, 0.8f); + + auto result = AnalyzeGridScan(scan, grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + + CHECK(c.nx == Catch::Approx(3.0f)); + CHECK(c.ny == Catch::Approx(3.0f)); + CHECK(c.x_um == Catch::Approx(30.0f)); + CHECK(c.y_um == Catch::Approx(30.0f)); + CHECK(c.image_number == 3 * 7 + 3); + CHECK(c.n_images == 5); + +// Score is the mean over the blob, not the peak, and every cell here carries the same 0.8 + CHECK(c.score == Catch::Approx(0.8f)); + CHECK(c.ice_score == Catch::Approx(0.0f)); + CHECK(std::isnan(c.res_A)); // no resolution was measured anywhere + +// Round: the two axes come out equal, whatever arbitrary angle the numerics picked + CHECK(c.major_um == Catch::Approx(c.minor_um).margin(1e-3)); + CHECK(c.angle_deg >= 0.0f); + CHECK(c.angle_deg < 180.0f); + +// The beam is not removed - it is reported so a consumer can remove it itself + CHECK(result.beam_size_x_um == Catch::Approx(5.0f)); + CHECK(result.beam_size_y_um == Catch::Approx(5.0f)); +} + +TEST_CASE("AnalyzeGridScan keeps a corner-touching diagonal needle whole", "[AnalyzeGridScan]") { +// A needle lying at 45 degrees across a raster whose step is coarser than the beam is sampled as +// cells that meet only at their corners. 4-connectivity would cut this into four blobs of one +// cell each and the minimum-size rule would then discard the crystal entirely; 8-connectivity +// keeps it as the one oriented needle it is. + const std::vector> diagonal = {{1, 1}, {2, 2}, {3, 3}, {4, 4}}; + + SECTION("square 20 x 20 um steps") { + auto grid = MakeGrid(7, 7, 20.0f, 20.0f); + auto result = AnalyzeGridScan(MakeScan(7, 7, diagonal, 0.9f), grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + + CHECK(c.n_images == 4); + CHECK(c.angle_deg == Catch::Approx(45.0f).margin(0.01)); + CHECK(c.nx == Catch::Approx(2.5f).margin(1e-4)); + CHECK(c.ny == Catch::Approx(2.5f).margin(1e-4)); + CHECK(c.major_um == Catch::Approx(113.137f).margin(0.01)); + CHECK(c.minor_um == Catch::Approx(28.284f).margin(0.01)); + } + + SECTION("anisotropic 20 x 16 um steps") { + auto grid = MakeGrid(7, 7, 20.0f, 16.0f); + auto result = AnalyzeGridScan(MakeScan(7, 7, diagonal, 0.9f), grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + + CHECK(c.n_images == 4); +// atan(16/20) - the same cells, read on a raster with a shorter y step + CHECK(c.angle_deg == Catch::Approx(38.6598f).margin(0.01)); + CHECK(c.x_um == Catch::Approx(50.0f).margin(1e-3)); + CHECK(c.y_um == Catch::Approx(40.0f).margin(1e-3)); + CHECK(c.major_um == Catch::Approx(102.450f).margin(0.01)); + CHECK(c.minor_um == Catch::Approx(24.988f).margin(0.01)); + } +} + +TEST_CASE("AnalyzeGridScan takes the needle angle from micrometres, not cells", "[AnalyzeGridScan]") { +// The SAME five cells - a staircase running up and to the right - read on two rasters that +// differ only in their y step. A purely diagonal line would not do: it is not 4-connected, and +// a real needle is sampled as a staircase anyway, the beam being wider than one step. + const std::vector> staircase = {{1, 1}, {2, 1}, {2, 2}, {3, 2}, {3, 3}}; + + SECTION("square 20 x 20 um steps") { + auto grid = MakeGrid(7, 7, 20.0f, 20.0f); + auto result = AnalyzeGridScan(MakeScan(7, 7, staircase, 0.9f), grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + +// Symmetric about the diagonal, so the axis is at exactly 45 degrees + CHECK(c.angle_deg == Catch::Approx(45.0f).margin(0.01)); + CHECK(c.nx == Catch::Approx(2.2f).margin(1e-4)); + CHECK(c.ny == Catch::Approx(1.8f).margin(1e-4)); + CHECK(c.x_um == Catch::Approx(44.0f).margin(1e-3)); + CHECK(c.y_um == Catch::Approx(36.0f).margin(1e-3)); + CHECK(c.major_um == Catch::Approx(84.853f).margin(0.01)); + CHECK(c.minor_um == Catch::Approx(42.426f).margin(0.01)); + } + + SECTION("anisotropic 20 x 16 um steps") { + auto grid = MakeGrid(7, 7, 20.0f, 16.0f); + auto result = AnalyzeGridScan(MakeScan(7, 7, staircase, 0.9f), grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + +// The cells are identical, so second moments taken in CELL units would still say 45 degrees. +// In micrometres the shorter y step tilts the axis down to 37.01 - eight degrees away, which +// is what makes this assertion worth making. + CHECK(c.angle_deg == Catch::Approx(37.010f).margin(0.01)); + CHECK(c.angle_deg < 44.0f); + +// The centre sits a fifth of a cell off a grid node in both directions + CHECK(c.nx == Catch::Approx(2.2f).margin(1e-4)); + CHECK(c.ny == Catch::Approx(1.8f).margin(1e-4)); + CHECK(c.x_um == Catch::Approx(44.0f).margin(1e-3)); + CHECK(c.y_um == Catch::Approx(28.8f).margin(1e-3)); + + CHECK(c.major_um == Catch::Approx(76.806f).margin(0.01)); + CHECK(c.minor_um == Catch::Approx(38.329f).margin(0.01)); + CHECK(c.major_um > c.minor_um); + } +} + +TEST_CASE("AnalyzeGridScan snaps a centre that falls outside the blob", "[AnalyzeGridScan]") { +// An L: a vertical arm at column 0 (rows 0..5) and a horizontal arm along row 5 (columns 1..3). +// The centroid of that is (0.667, 3.333), which rounds to cell (1,3) - a cell that is NOT in +// the blob and was never collected as part of this crystal. + auto grid = MakeGrid(6, 6, 10.0f, 10.0f); + auto scan = MakeScan(6, 6, {{0, 0}, {0, 1}, {0, 2}, {0, 3}, {0, 4}, {0, 5}, + {1, 5}, {2, 5}, {3, 5}}, 0.7f); + + auto result = AnalyzeGridScan(scan, grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + + CHECK(c.nx == Catch::Approx(6.0f / 9.0f).margin(1e-4)); + CHECK(c.ny == Catch::Approx(30.0f / 9.0f).margin(1e-4)); + CHECK(c.n_images == 9); + +// The reported image is the nearest cell that is actually in the blob, (0,3), not (1,3) + CHECK(c.image_number == 3 * 6 + 0); +} + +TEST_CASE("AnalyzeGridScan centres on resolution rank, not resolution value", "[AnalyzeGridScan]") { +// A horizontal bar of five cells. The best resolution sits at one end, so it pulls the centre +// off the geometric middle - but by exactly the same amount whether that end reports a +// plausible 2.0 A or an absurd 0.8 A, because only the rank is used. + auto grid = MakeGrid(7, 3, 10.0f, 10.0f); + const std::vector> bar = {{1, 1}, {2, 1}, {3, 1}, {4, 1}, {5, 1}}; + + auto with_best_at_left = [&](float best_res) { + auto scan = MakeScan(7, 3, bar, 0.6f); + const float res[5] = {best_res, 3.0f, 3.5f, 4.0f, 4.5f}; + for (int i = 0; i < 5; i++) + scan.images[1 * 7 + 1 + i].res = res[i]; + return AnalyzeGridScan(scan, grid, 5.0f, 5.0f).crystals.at(0); + }; + + const auto plausible = with_best_at_left(2.0f); + const auto artefact = with_best_at_left(0.8f); + +// Pulled towards the good end, so left of the geometric centre at 3.0 + CHECK(plausible.nx < 3.0f); +// ...and the absurd number moves it not one bit + CHECK(artefact.nx == Catch::Approx(plausible.nx)); + CHECK(artefact.ny == Catch::Approx(plausible.ny)); + +// res_A is the 25th percentile of the five, i.e. the second best (3.0), never the best + CHECK(plausible.res_A == Catch::Approx(3.0f)); + CHECK(artefact.res_A == Catch::Approx(3.0f)); +} + +TEST_CASE("AnalyzeGridScan separates two crystals and sorts them by score", "[AnalyzeGridScan]") { +// Two blobs with a clear gap between them; the weaker one is collected first. Labelling is +// 8-connected, so the gap has to be more than one cell in EVERY direction, diagonals included - +// three empty columns here. + auto grid = MakeGrid(9, 5, 10.0f, 10.0f); + auto scan = MakeScan(9, 5, {}, 0.0f); + + for (auto [x, y]: std::vector>{{1, 1}, {1, 2}, {2, 1}, {2, 2}}) + scan.images[y * 9 + x].protein_score = 0.6f; + for (auto [x, y]: std::vector>{{6, 2}, {7, 2}, {6, 3}}) + scan.images[y * 9 + x].protein_score = 0.95f; + + scan.images[2 * 9 + 6].ice_score = 0.4f; + + auto result = AnalyzeGridScan(scan, grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 2); + +// Element 0 is the one to collect + CHECK(result.crystals[0].score == Catch::Approx(0.95f)); + CHECK(result.crystals[0].n_images == 3); + CHECK(result.crystals[0].nx > 4.0f); + CHECK(result.crystals[0].ice_score == Catch::Approx(0.4f / 3.0f)); + + CHECK(result.crystals[1].score == Catch::Approx(0.6f)); + CHECK(result.crystals[1].n_images == 4); + CHECK(result.crystals[1].nx < 4.0f); + CHECK(result.crystals[1].ice_score == Catch::Approx(0.0f)); +} + +TEST_CASE("AnalyzeGridScan reports major_um along angle_deg", "[AnalyzeGridScan]") { +// A hook: two cells down, a step across, one more down. On a 20 x 10 um raster the axis of the +// larger second MOMENT runs at 58.3 degrees, but the blob's measured extent along that axis +// (46.5 um) is SHORTER than across it (49.8 um) - the moment is dominated by two clumps at the +// ends, the extent is not. A consumer draws a major by minor frame rotated by angle_deg, so the +// two facts are kept by turning the frame a quarter turn, not by discarding one of them. + auto grid = MakeGrid(5, 5, 20.0f, 10.0f); + auto scan = MakeScan(5, 5, {{0, 0}, {0, 1}, {1, 1}, {1, 2}, {0, 3}}, 0.8f); + + auto result = AnalyzeGridScan(scan, grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + + CHECK(c.n_images == 5); + CHECK(c.major_um == Catch::Approx(49.798f).margin(0.01)); + CHECK(c.minor_um == Catch::Approx(46.549f).margin(0.01)); +// 58.28 + 90, folded back into [0,180) + CHECK(c.angle_deg == Catch::Approx(148.283f).margin(0.01)); +} + +TEST_CASE("AnalyzeGridScan keeps major_um >= minor_um for every blob", "[AnalyzeGridScan]") { +// The invariant a consumer relies on, over every shape the other cases exercise. + const std::vector> shapes[] = { + {{3, 2}, {2, 3}, {3, 3}, {4, 3}, {3, 4}}, // round + {{1, 1}, {2, 2}, {3, 3}, {4, 4}}, // diagonal needle + {{1, 1}, {2, 1}, {2, 2}, {3, 2}, {3, 3}}, // staircase needle + {{0, 0}, {0, 1}, {1, 1}, {1, 2}, {0, 3}}, // hook + {{0, 0}, {0, 1}, {0, 2}, {0, 3}, {0, 4}, {0, 5}, + {1, 5}, {2, 5}, {3, 5}}, // L + {{1, 1}, {2, 1}, {3, 1}, {4, 1}, {5, 1}}, // bar + }; + + for (const auto &shape: shapes) { + for (auto [sx, sy]: std::vector>{{10.0f, 10.0f}, + {20.0f, 16.0f}, + {20.0f, 10.0f}, + {10.0f, 20.0f}}) { + auto grid = MakeGrid(7, 7, sx, sy); + auto result = AnalyzeGridScan(MakeScan(7, 7, shape, 0.8f), grid, 5.0f, 5.0f); + REQUIRE(result.crystals.size() == 1); + const auto &c = result.crystals[0]; + + CHECK(c.major_um >= c.minor_um); + CHECK(c.angle_deg >= 0.0f); + CHECK(c.angle_deg < 180.0f); + } + } +} + +TEST_CASE("AnalyzeGridScan follows the snake and the step signs", "[AnalyzeGridScan]") { +// The same crystal, described once on a plain raster and once on a snake raster with a +// negative slow step. Rearrange puts both on the same display grid, so both have to report +// the same position - only the image numbers differ. + const int64_t nx = 6, ny = 6; + GridScanSettings plain(nx, 10.0f, 10.0f, false, false); + plain.ImageNum(nx * ny); + GridScanSettings snake(nx, 10.0f, -10.0f, true, false); + snake.ImageNum(nx * ny); + + const std::vector> blob = {{2, 1}, {3, 1}, {2, 2}, {3, 2}}; + + auto scan_plain = MakeScan(nx, ny, blob, 0.8f); + ScanResult scan_snake; + for (int64_t i = 0; i < nx * ny; i++) { + ScanResultElem elem; + elem.number = i; + elem.protein_score = scan_plain.images[snake.Rearrange(i)].protein_score; + scan_snake.images.push_back(elem); + } + + auto a = AnalyzeGridScan(scan_plain, plain, 5.0f, 5.0f).crystals.at(0); + auto b = AnalyzeGridScan(scan_snake, snake, 5.0f, 5.0f).crystals.at(0); + + CHECK(a.nx == Catch::Approx(b.nx)); + CHECK(a.ny == Catch::Approx(b.ny)); + CHECK(a.x_um == Catch::Approx(b.x_um)); + CHECK(a.y_um == Catch::Approx(b.y_um)); + CHECK(snake.Rearrange(b.image_number) == plain.Rearrange(a.image_number)); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index a4af137b6..df33f23af 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -60,6 +60,7 @@ ADD_EXECUTABLE(jfjoch_test ImageMetadataTest.cpp JFJochReceiverLiteTest.cpp GridScanSettingsTest.cpp + AnalyzeGridScanTest.cpp JFJochReceiverPlotsTest.cpp GoniometerAxisTest.cpp DetGeomCalibTest.cpp