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