From 273591af11a2fdaa49e1b7ca9d0eb6271b8ceeb4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 8 Sep 2026 07:34:36 +0200 Subject: [PATCH] grid scan: a crystal is grown out of its seeds, so one weak grid point no longer splits it A crystal in the corpus was reported as two because the grid point joining its halves scored 0.498 - two thousandths under the threshold. That is an artefact of putting a hard edge through a continuous quantity, not a gap in a crystal. The patch search is now hysteresis. A patch is labelled at a lower grow threshold (0.35) and kept only if it holds a cell above the seed threshold (0.50), which is the same thing as growing out of the seeds in one pass of the labeller rather than two. A cell between the two levels joins a crystal that already exists but can never start one, so no lowering of the grow threshold can turn a background into a crystal - and the negatives are safe by construction, since water peaks at 0.146 and ice at 0.490, both under the seed level. Both halves of the small-patch rule read the GROWN patch - the count includes the cells growth added and the peak is the patch's best cell wherever it lies - which is stated at the test, because a reader will otherwise wonder whether a grown cell can rescue a patch no seed would have admitted. It cannot: the patch is discarded before that test unless it holds a seed. The reported extents are those of the grown patch, so they reach the 0.35 contour. Over the corpus that is major_um x1.04 and n_images x1.07 at the median (x1.11 and x1.14 at the mean, +7% cells overall), and the cells it adds are spread fairly evenly over 0.35-0.50 rather than piled at the bottom - a crystal edge sampled at one grid step, not bleed into a neighbour that never diffracted. The patch mean falls with them, and can now sit below the seed threshold; the peak beside it is what the admission was decided on. RASTER_REPORT_VERSION is 3, because N_CELLS and the extents are a different measurement under the same names. The corpus table is unchanged: 17/17 protein, 0/4 water, 0/3 ice, 10/10 heldout. The split heals into one crystal, and one raster's five patches become three. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- docs/CHANGELOG.md | 1 + docs/RUGNUX_ADVANCED.md | 1 + .../grid_scan_analysis/AnalyzeGridScan.cpp | 25 ++++++++++-- .../grid_scan_analysis/AnalyzeGridScan.h | 8 ++++ rugnux/RasterReport.cpp | 16 +++++--- rugnux/rugnux_cli.cpp | 9 ++++- tests/RasterReportTest.cpp | 40 ++++++++++++++++++- 7 files changed, 90 insertions(+), 10 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 53935c0df..26b7ebd27 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -28,6 +28,7 @@ ### 1.0.0-rc.166 * `rugnux --mode raster` scores every image of a stored grid scan and reports the crystals in it - grid position, nearest image, measured size and orientation, score, ice score and resolution - as `_raster_report.txt` and `_raster.json`; `--beam-size` states the beam the sizes still contain. +* `rugnux --mode raster` grows a crystal out to `--raster-grow-score` once one of its grid points passes `--raster-protein-threshold`, so a crystal is no longer split in two by a single point that fell just under the higher bar; the reported extents and grid-point count are those of the grown patch. * `rugnux --mode raster` reports a patch smaller than `--raster-min-cells` when its best cell reaches `--raster-decisive-score`, so a crystal covering a single grid point is found where that point is clearly protein, and `--raster-max-crystals` keeps the report to the best few. * `rugnux --model` treats the model as a hypothesis: it decides the enantiomorph and the indexing only where its R-work beats that of the same model in random orientations, and a model the data reject is still scored, placed and mapped, but leaves the reflection files byte for byte what a run with no model writes. * `rugnux --model` places the model against the data as a rigid body before scoring it, writes sigma_A-weighted 2mFo-DFc and mFo-DFc maps in place of the unweighted 2Fo-Fc and Fo-Fc, and writes the model as it was placed - `_model.cif`, and `_model.pdb` where the PDB format can express the cell - in the cell and space group of the reflection files beside it. diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 303a968b8..679d252af 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -224,6 +224,7 @@ Raster (`--mode raster`): | --- | --- | | `--beam-size [,]` | Beam size at the sample in um, x then y; one value is a square beam. Default: the file's `incident_beam_size`. The reported crystal sizes are measured and still contain the beam, so this says what a consumer has to take back out — and since that removal is a subtraction of two covariance matrices, a beam given as square when it is not rotates the reported crystal axis | | `--raster-protein-threshold ` | Protein score a cell must exceed to count as diffracting (default `0.50`) | +| `--raster-grow-score ` | Once a patch has a cell above `--raster-protein-threshold` it is grown out to this score (default `0.35`), so a crystal is not broken in two by one cell that fell just under the higher bar. Growth cannot start on its own — a patch that never reaches the higher bar is not a patch — so lowering this cannot invent a crystal | | `--raster-min-cells ` | Cells a connected patch must have before it is reported as a crystal (default `3`). A smaller patch is still reported when its best cell reaches `--raster-decisive-score` | | `--raster-decisive-score ` | Protein score at which one cell is evidence enough on its own, so a patch below `--raster-min-cells` is still reported (default `0.60`). Measured over 67 labelled rasters the populations do not overlap: no water raster reaches 0.15 and no ice raster 0.50, while the weakest protein raster peaks at 0.67 | | `--raster-max-crystals ` | Report at most this many crystals, the best-scoring ones (default `0`, all of them) | diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp index cd01b2fa3..fc140701d 100644 --- a/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp @@ -40,14 +40,19 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, image_no[i] = elem.number; } - std::vector above(nx * ny); + // Hysteresis. The patches are labelled at the GROW level and then only those holding at least + // one cell at the SEED level are kept, which is the same thing as growing out of the seeds and + // is one pass of the labeller rather than two. A cell between the two levels therefore joins a + // crystal that already exists but can never start one, so however low the grow level is set, a + // background that never reaches the seed level produces nothing. + std::vector grown(nx * ny); for (int64_t i = 0; i < nx * ny; i++) - above[i] = (protein[i] > params.protein_score_threshold) ? 1 : 0; + grown[i] = (protein[i] > params.grow_score_threshold) ? 1 : 0; // Labelled with no size cut of its own: whether a patch is big enough is no longer a plain // floor - a small patch survives on the strength of its diffraction - and that test needs the // patch's scores, which the labeller does not have. - const std::vector label = LabelConnectedComponents(above, nx, ny, 1); + const std::vector label = LabelConnectedComponents(grown, nx, ny, 1); const int32_t n_label = label.empty() ? 0 : *std::max_element(label.begin(), label.end()); GridScanResult result; @@ -63,6 +68,14 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, if (cell.empty()) continue; + // Hysteresis, the second half: a patch that never reaches the seed level is not a patch, + // whatever its shape, so nothing about it is worth measuring. + bool has_seed = false; + for (int64_t i: cell) + has_seed = has_seed || protein[i] > params.protein_score_threshold; + if (!has_seed) + continue; + const auto n = static_cast(cell.size()); // Cell centres in MICROMETRES, taken once. Every geometric quantity below is computed from @@ -157,6 +170,12 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, // A patch is a crystal when it is big enough to be a shape rather than a coincidence, OR // when it is smaller than that but the diffraction in it is decisive on its own. One cell // is enough where that cell is clearly protein; a weak patch still has to be a shape. + // + // Both halves are read over the GROWN patch: the count includes the cells hysteresis added + // and the peak is the patch's best cell wherever it lies. That cannot let growth rescue a + // patch no seed would have admitted, because the patch has already been discarded above + // unless it holds a seed cell - and a seed cell is by definition the strongest kind there + // is, so the peak of a grown patch is the peak of its seeds. if (static_cast(cell.size()) < params.min_blob_cells && peak_protein < params.decisive_protein_score) continue; diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.h b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h index 75b20f7aa..c8f15041e 100644 --- a/image_analysis/grid_scan_analysis/AnalyzeGridScan.h +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h @@ -32,6 +32,13 @@ constexpr float PROTEIN_SCORE_THRESHOLD_DEFAULT = 0.5f; // below. constexpr int64_t MIN_BLOB_CELLS_DEFAULT = 3; +// A patch is GROWN out to this score once it has started, so a crystal is not broken in two by a +// single cell that fell just under the seed threshold - a real split in the corpus had its bridging +// cell at 0.498, two thousandths under it, which is a threshold artefact and not a gap in a crystal. +// Growth can never start on its own: a patch reaching only this level and never the seed level is +// not a patch at all, so no lowering of this can turn a weak background into a crystal. +constexpr float GROW_SCORE_THRESHOLD_DEFAULT = 0.35f; + // A patch smaller than the minimum is kept anyway when its BEST cell scores at least this, so one // cell is enough where that cell is clearly protein while a weak patch still has to be a shape. // The bar is the peak and not the patch mean, because a two-cell patch with one strong cell and one @@ -51,6 +58,7 @@ constexpr int64_t MAX_CRYSTALS_DEFAULT = 0; // change here or at any call site. struct GridScanAnalysisParameters { float protein_score_threshold = PROTEIN_SCORE_THRESHOLD_DEFAULT; + float grow_score_threshold = GROW_SCORE_THRESHOLD_DEFAULT; int64_t min_blob_cells = MIN_BLOB_CELLS_DEFAULT; float decisive_protein_score = DECISIVE_PROTEIN_SCORE_DEFAULT; int64_t max_crystals = MAX_CRYSTALS_DEFAULT; diff --git a/rugnux/RasterReport.cpp b/rugnux/RasterReport.cpp index 4a6de1c0b..69f97f6d5 100644 --- a/rugnux/RasterReport.cpp +++ b/rugnux/RasterReport.cpp @@ -16,7 +16,7 @@ namespace { // The version of this file format. Bumped when a key is renamed or removed or a table column // moves - a consumer can gate on it. - constexpr int RASTER_REPORT_VERSION = 2; + constexpr int RASTER_REPORT_VERSION = 3; const char *BANNER = " ******************************************************************************"; @@ -113,13 +113,17 @@ std::string RenderRasterReport(const std::string &input_file, // -------------------------------------------------------------------- 2. PER-IMAGE SCORES Section(os, "2. PER-IMAGE SCORES"); os << " Every image carries a protein score and an ice score, each between 0 and 1. A cell counts\n" - << " as diffracting above PROTEIN_SCORE_THRESHOLD, and the crystals below are the connected\n" - << " patches of at least MIN_BLOB_CELLS such cells - or smaller patches whose best cell reaches\n" - << " DECISIVE_PROTEIN_SCORE, so one cell is enough where that cell is clearly protein.\n" + << " as diffracting above PROTEIN_SCORE_THRESHOLD, and a patch holding one grows out to\n" + << " GROW_SCORE_THRESHOLD, so a crystal is not broken in two by a single cell that fell just\n" + << " under the higher bar; a patch that never reaches the higher bar is not a patch at all.\n" + << " The crystals below are the grown patches of at least MIN_BLOB_CELLS cells - or smaller\n" + << " patches whose best cell reaches DECISIVE_PROTEIN_SCORE, so one cell is enough where that\n" + << " cell is clearly protein.\n" << " ICE_ABOVE_THRESHOLD counts the cells whose ICE score passes the protein threshold - the one\n" << " number the run already has, rather than a second one invented for ice - so it says how much\n" << " of the loop is ice rather than sample.\n\n"; Key(os, "PROTEIN_SCORE_THRESHOLD", fmt::format("{:.2f}", settings.analysis.protein_score_threshold)); + Key(os, "GROW_SCORE_THRESHOLD", fmt::format("{:.2f}", settings.analysis.grow_score_threshold)); Key(os, "MIN_BLOB_CELLS", settings.analysis.min_blob_cells); Key(os, "DECISIVE_PROTEIN_SCORE", fmt::format("{:.2f}", settings.analysis.decisive_protein_score)); Key(os, "MAX_CRYSTALS", settings.analysis.max_crystals); @@ -141,7 +145,8 @@ std::string RenderRasterReport(const std::string &input_file, << " still contain the beam (BEAM_SIZE_UM above says which one) - nothing here deconvolves it.\n" << " RES_A is the robust best resolution inside the patch, blank where none was measured. SCORE\n" << " is the patch mean and ranks the rows; PEAK is its best cell, which is what admits a patch\n" - << " with fewer than MIN_BLOB_CELLS cells in it.\n\n"; + << " with fewer than MIN_BLOB_CELLS cells in it. N_CELLS and the extents describe the GROWN\n" + << " patch, so they reach the GROW_SCORE_THRESHOLD contour and not the seed one.\n\n"; Key(os, "CRYSTAL_COUNT", crystals.crystals.size()); os << "\n" << " NX NY X_UM Y_UM IMAGE MAJOR_UM MINOR_UM ANGLE_DEG SCORE PEAK ICE_SCORE RES_A N_CELLS\n" @@ -196,6 +201,7 @@ std::string RenderRasterJson(const std::string &input_file, nlohmann::json s; s["protein_score_threshold"] = settings.analysis.protein_score_threshold; + s["grow_score_threshold"] = settings.analysis.grow_score_threshold; s["min_blob_cells"] = settings.analysis.min_blob_cells; s["decisive_protein_score"] = settings.analysis.decisive_protein_score; s["max_crystals"] = settings.analysis.max_crystals; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 41b76a489..57ed77010 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -127,6 +127,7 @@ void print_usage() { std::cout << " Raster (--mode raster)" << std::endl; std::cout << " --beam-size [,] Beam size at the sample in um, x then y; one value is a square beam. Default: the file's incident_beam_size. The reported crystal sizes are MEASURED and still contain the beam, so this is what a consumer takes back out - and because that removal is a covariance subtraction, a beam given as square when it is not rotates the reported crystal axis" << std::endl; std::cout << " --raster-protein-threshold Protein score a cell must exceed to count as diffracting (default: 0.50)" << std::endl; + std::cout << " --raster-grow-score Once a patch has a cell above --raster-protein-threshold it is grown out to this score (default: 0.35), so a crystal is not broken in two by one cell that fell just under the higher bar. A patch that never reaches the higher bar is not a patch at all, so lowering this cannot invent one" << std::endl; std::cout << " --raster-min-cells Cells a connected patch must have before it is reported as a crystal (default: 3). A smaller patch is still reported when its best cell reaches --raster-decisive-score" << std::endl; std::cout << " --raster-decisive-score Protein score at which one cell is evidence enough on its own, so a patch below --raster-min-cells is still reported (default: 0.60). Over 67 labelled rasters no water raster reaches 0.15 and no ice raster 0.50, while the weakest protein raster peaks at 0.67" << std::endl; std::cout << " --raster-max-crystals Report at most this many crystals, the best-scoring ones (default: 0, all of them)" << std::endl; @@ -346,7 +347,8 @@ enum { OPT_RASTER_PROTEIN_THRESHOLD, OPT_RASTER_MIN_CELLS, OPT_RASTER_DECISIVE_SCORE, - OPT_RASTER_MAX_CRYSTALS + OPT_RASTER_MAX_CRYSTALS, + OPT_RASTER_GROW_SCORE }; static option long_options[] = { @@ -374,6 +376,7 @@ static option long_options[] = { {"beam-size", required_argument, nullptr, OPT_BEAM_SIZE}, {"raster-protein-threshold", required_argument, nullptr, OPT_RASTER_PROTEIN_THRESHOLD}, {"raster-min-cells", required_argument, nullptr, OPT_RASTER_MIN_CELLS}, + {"raster-grow-score", required_argument, nullptr, OPT_RASTER_GROW_SCORE}, {"raster-decisive-score", required_argument, nullptr, OPT_RASTER_DECISIVE_SCORE}, {"raster-max-crystals", required_argument, nullptr, OPT_RASTER_MAX_CRYSTALS}, {"no-refine-tilt", no_argument, nullptr, OPT_NO_REFINE_TILT}, @@ -1159,6 +1162,10 @@ static int RunRugnux(int argc, char **argv) { raster_analysis.min_blob_cells = parse_number_arg(optarg, "--raster-min-cells", logger, 1); break; + case OPT_RASTER_GROW_SCORE: + raster_analysis.grow_score_threshold = + parse_number_arg(optarg, "--raster-grow-score", logger, 0.0f, 1.0f); + break; case OPT_RASTER_DECISIVE_SCORE: raster_analysis.decisive_protein_score = parse_number_arg(optarg, "--raster-decisive-score", logger, 0.0f, 1.0f); diff --git a/tests/RasterReportTest.cpp b/tests/RasterReportTest.cpp index 6fbfa3817..32accbf9a 100644 --- a/tests/RasterReportTest.cpp +++ b/tests/RasterReportTest.cpp @@ -59,7 +59,7 @@ TEST_CASE("RasterReport_Render", "[Diagnostics]") { // The keys are the interface a battery greps; pin the spellings and the values that say what // the run saw. - CHECK(report.find("RASTER_REPORT_VERSION= 2\n") != std::string::npos); + CHECK(report.find("RASTER_REPORT_VERSION= 3\n") != std::string::npos); CHECK(report.find("INPUT_FILE= test_master.h5\n") != std::string::npos); CHECK(report.find("GRID_SIZE= 5 4\n") != std::string::npos); CHECK(report.find("GRID_STEP_UM= 20.00 16.00\n") != std::string::npos); @@ -69,6 +69,7 @@ TEST_CASE("RasterReport_Render", "[Diagnostics]") { CHECK(report.find("IMAGES_SCORED= 20\n") != std::string::npos); CHECK(report.find("IMAGES_ABOVE_THRESHOLD= 4\n") != std::string::npos); CHECK(report.find("ICE_ABOVE_THRESHOLD= 0\n") != std::string::npos); + CHECK(report.find("GROW_SCORE_THRESHOLD= 0.35\n") != std::string::npos); CHECK(report.find("MIN_BLOB_CELLS= 3\n") != std::string::npos); CHECK(report.find("DECISIVE_PROTEIN_SCORE= 0.60\n") != std::string::npos); CHECK(report.find("MAX_CRYSTALS= 0\n") != std::string::npos); @@ -232,3 +233,40 @@ TEST_CASE("RasterReport_OneDecisiveCellIsACrystal", "[Diagnostics]") { CHECK(AnalyzeGridScan(big, grid, settings.beam_size_x_um, settings.beam_size_y_um, settings.analysis).crystals.size() == 1); } + +TEST_CASE("RasterReport_HysteresisHealsASplit", "[Diagnostics]") { + // The case this exists for, taken from the corpus: one crystal whose middle grid point fell two + // thousandths under the seed threshold, which broke it into two crystals that were never two. + const GridScanSettings grid = WideGrid(); + ScanResult scan = WideScan({}); + for (int64_t x = 2; x <= 8; x++) + scan.images[2 * WIDE_NX + x].protein_score = 0.85f; + scan.images[2 * WIDE_NX + 5].protein_score = 0.498f; // above the grow level, under the seed + + RasterSettings settings = TestSettings(); + const GridScanResult healed = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, + settings.beam_size_y_um, settings.analysis); + REQUIRE(healed.crystals.size() == 1); + CHECK(healed.crystals[0].n_images == 7); // the bridging cell is part of the crystal + + // Growing no further than the seed level is the behaviour hysteresis replaced, and it splits. + settings.analysis.grow_score_threshold = settings.analysis.protein_score_threshold; + const GridScanResult split = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, + settings.beam_size_y_um, settings.analysis); + CHECK(split.crystals.size() == 2); +} + +TEST_CASE("RasterReport_GrowthCannotStartOnItsOwn", "[Diagnostics]") { + // A patch that reaches the grow level and never the seed level is not a crystal however large it + // is or however low the grow level is set - which is what makes the level safe to lower. + const GridScanSettings grid = WideGrid(); + const ScanResult scan = WideScan({{1, 1, 0.45f}, {8, 2, 0.40f}}); + + RasterSettings settings = TestSettings(); + CHECK(AnalyzeGridScan(scan, grid, settings.beam_size_x_um, settings.beam_size_y_um, + settings.analysis).crystals.empty()); + + settings.analysis.grow_score_threshold = 0.05f; + CHECK(AnalyzeGridScan(scan, grid, settings.beam_size_x_um, settings.beam_size_y_um, + settings.analysis).crystals.empty()); +}