From cdb60c29167fcc6d24829e66f34758f8af981e07 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 8 Sep 2026 07:15:30 +0200 Subject: [PATCH] grid scan: one cell is a crystal when that cell is decisive, and the list can be capped The minimum patch size was a flat floor, so a crystal that covers one or two grid points was thrown away however strong its diffraction. It is now a floor OR a piece of evidence: a patch is reported when it has min_blob_cells cells, or when its best cell reaches decisive_protein_score. One condition, and both halves are parameters. The bar is the patch PEAK, not the patch mean. A two-cell patch with one strong cell and one marginal one is the case this exists for, and the mean averages exactly that evidence away. Over the 67 labelled rasters the two populations do not overlap: no water raster reaches a peak of 0.15 and no ice raster reaches 0.50, while the weakest protein raster peaks at 0.67 - so 0.6, the middle of that gap, is the default. The peak is reported beside the mean, in the table and in the JSON, so an admission can be checked against the number that decided it. max_crystals caps the returned list after the sort, best first; 0, the default, is all of them, because a crystal that was found and then dropped is information the caller cannot get back. The four parameters now travel as one GridScanAnalysisParameters, which is the shape the GridScanAnalysisSettings class in common/ is to take: when it lands it replaces this struct in the signature and nothing else changes. On the corpus this reaches 17/17 protein loops (it was 16/17 - the miss was a crystal covering two grid points of a 4x4 raster) with water still 0/4 and ice still 0/3. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N --- docs/CHANGELOG.md | 1 + docs/RUGNUX_ADVANCED.md | 4 +- .../grid_scan_analysis/AnalyzeGridScan.cpp | 29 +++- .../grid_scan_analysis/AnalyzeGridScan.h | 34 ++++- .../grid_scan_analysis/GridScanResult.h | 6 +- rugnux/RasterReport.cpp | 41 +++--- rugnux/RasterReport.h | 4 +- rugnux/rugnux_cli.cpp | 38 ++++-- tests/AnalyzeGridScanTest.cpp | 7 +- tests/RasterReportTest.cpp | 129 ++++++++++++++++-- 10 files changed, 234 insertions(+), 59 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index cc5288c08..53935c0df 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` 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. * `rugnux` and `jfjoch_viewer` read PILATUS miniCBF sweeps natively, and open masters written at other facilities, including Eiger 1.x and third-party NXmx. diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index c2777fc98..303a968b8 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -224,7 +224,9 @@ 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-min-cells ` | Cells a connected patch must have before it is reported as a crystal (default `3`) | +| `--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) | Calibration (`--mode calibration`): diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp index 6fdbce02a..cd01b2fa3 100644 --- a/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp @@ -17,8 +17,7 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, const GridScanSettings &grid, float beam_size_x_um, float beam_size_y_um, - float protein_score_threshold, - int64_t min_blob_cells) { + const GridScanAnalysisParameters ¶ms) { const int64_t nx = grid.GetGridSizeX_step(); const int64_t ny = grid.GetGridSizeY_step(); const float step_x = fabsf(grid.GetGridStepX_um()); @@ -43,9 +42,12 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, std::vector above(nx * ny); for (int64_t i = 0; i < nx * ny; i++) - above[i] = (protein[i] > protein_score_threshold) ? 1 : 0; + above[i] = (protein[i] > params.protein_score_threshold) ? 1 : 0; - const std::vector label = LabelConnectedComponents(above, nx, ny, min_blob_cells); + // 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 int32_t n_label = label.empty() ? 0 : *std::max_element(label.begin(), label.end()); GridScanResult result; @@ -59,7 +61,7 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, cell.push_back(i); } if (cell.empty()) - continue; // label fell below min_blob_cells and was erased + continue; const auto n = static_cast(cell.size()); @@ -145,12 +147,20 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, 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; + float sum_protein = 0, sum_ice = 0, peak_protein = 0; for (int64_t i: cell) { sum_protein += protein[i]; sum_ice += std::max(ice[i], 0.0f); + peak_protein = std::max(peak_protein, protein[i]); } + // 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. + if (static_cast(cell.size()) < params.min_blob_cells + && peak_protein < params.decisive_protein_score) + continue; + GridScanCrystal crystal; crystal.nx = cx_um / step_x; crystal.ny = cy_um / step_y; @@ -181,6 +191,7 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, // real crystal and ranks nothing. The mean stays a detection confidence and compares. crystal.score = sum_protein / n; crystal.ice_score = sum_ice / n; + crystal.peak_score = peak_protein; // 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. @@ -196,5 +207,11 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, std::sort(result.crystals.begin(), result.crystals.end(), [](const GridScanCrystal &a, const GridScanCrystal &b) { return a.score > b.score; }); + // Kept to the best few only where a caller asked for that; the sort above is what makes the + // ones it keeps the right ones. + if (params.max_crystals > 0 + && result.crystals.size() > static_cast(params.max_crystals)) + result.crystals.resize(params.max_crystals); + return result; } diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.h b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h index d34d3a258..75b20f7aa 100644 --- a/image_analysis/grid_scan_analysis/AnalyzeGridScan.h +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h @@ -22,22 +22,42 @@ // 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. -// -// The last two arguments are the whole decision: which cells count and how many of them make a -// shape. They default to what the corpus was tuned on and are parameters rather than constants so a -// sweep can vary them without a rebuild (rugnux --raster-protein-threshold / --raster-min-cells). // 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_DEFAULT = 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. +// patch that is a shape rather than a coincidence - unless the diffraction in it is decisive, see +// below. constexpr int64_t MIN_BLOB_CELLS_DEFAULT = 3; +// 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 +// marginal one is exactly the case this exists for and the mean would average the evidence away. +// Measured over 67 labelled rasters the two populations do not overlap at all: no water raster +// reaches 0.15 and no ice raster reaches 0.50, while the weakest protein raster peaks at 0.67. 0.6 +// is the middle of that gap. +constexpr float DECISIVE_PROTEIN_SCORE_DEFAULT = 0.6f; + +// How many crystals are returned, best first; 0 is all of them, which is the default because a +// crystal that was found and then dropped is information the caller cannot get back. +constexpr int64_t MAX_CRYSTALS_DEFAULT = 0; + +// What counts as a crystal. Grouped rather than passed loose so that the GridScanAnalysisSettings +// class in common/ can adopt it whole when it lands: these are the members and the defaults it +// should carry, and AnalyzeGridScan then takes that class in place of this struct with no other +// change here or at any call site. +struct GridScanAnalysisParameters { + float protein_score_threshold = PROTEIN_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; +}; + GridScanResult AnalyzeGridScan(const ScanResult &scan, const GridScanSettings &grid, float beam_size_x_um, float beam_size_y_um, - float protein_score_threshold = PROTEIN_SCORE_THRESHOLD_DEFAULT, - int64_t min_blob_cells = MIN_BLOB_CELLS_DEFAULT); + const GridScanAnalysisParameters ¶ms = {}); diff --git a/image_analysis/grid_scan_analysis/GridScanResult.h b/image_analysis/grid_scan_analysis/GridScanResult.h index 770152659..1259f4846 100644 --- a/image_analysis/grid_scan_analysis/GridScanResult.h +++ b/image_analysis/grid_scan_analysis/GridScanResult.h @@ -23,7 +23,11 @@ struct GridScanCrystal { // 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 score = 0, ice_score = 0; // 0-1, the patch MEAN + // The best protein score in the patch. The mean above ranks the crystals; this says how strong + // the evidence in the patch ever gets, and it is what admits a patch too small to be a shape + // (AnalyzeGridScan). Reported so that admission can be checked against the number that decided it. + float peak_score = 0; float res_A = NAN; // robust best resolution in the blob, NaN if none was measured int64_t n_images = 0; }; diff --git a/rugnux/RasterReport.cpp b/rugnux/RasterReport.cpp index 09fdebc26..4a6de1c0b 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 = 1; + constexpr int RASTER_REPORT_VERSION = 2; const char *BANNER = " ******************************************************************************"; @@ -73,7 +73,7 @@ std::string RenderRasterReport(const std::string &input_file, const RasterSettings &settings, const RunProvenance &provenance) { std::ostringstream os; - const ScanSummary summary = Summarize(scan, settings.protein_score_threshold); + const ScanSummary summary = Summarize(scan, settings.analysis.protein_score_threshold); os << BANNER << "\n" << " RUGNUX RASTER REPORT\n" @@ -114,11 +114,15 @@ std::string RenderRasterReport(const std::string &input_file, 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. ICE_ABOVE_THRESHOLD counts the cells whose\n" - << " ICE score passes that same threshold - the one number the run already has, rather than a\n" - << " second one invented for ice - so it says how much of the loop is ice rather than sample.\n\n"; - Key(os, "PROTEIN_SCORE_THRESHOLD", fmt::format("{:.2f}", settings.protein_score_threshold)); - Key(os, "MIN_BLOB_CELLS", settings.min_blob_cells); + << " 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" + << " 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, "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); Key(os, "IMAGES_SCORED", summary.scored); Key(os, "IMAGES_ABOVE_THRESHOLD", summary.above); Key(os, "FRACTION_ABOVE_THRESHOLD", fmt::format("{:.4f}", Fraction(summary.above, summary.scored))); @@ -135,22 +139,24 @@ std::string RenderRasterReport(const std::string &input_file, << " extent along the crystal's own axes and ANGLE_DEG points along MAJOR_UM, counter-clockwise\n" << " from +x and folded into [0,180) because an axis has no sign. The sizes are MEASURED and\n" << " 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.\n\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"; Key(os, "CRYSTAL_COUNT", crystals.crystals.size()); os << "\n" - << " NX NY X_UM Y_UM IMAGE MAJOR_UM MINOR_UM ANGLE_DEG SCORE ICE_SCORE RES_A N_CELLS\n" - << " -------- ------ -------- -------- --------- -------- -------- --------- ------ --------- ------ -------\n"; + << " NX NY X_UM Y_UM IMAGE MAJOR_UM MINOR_UM ANGLE_DEG SCORE PEAK ICE_SCORE RES_A N_CELLS\n" + << " -------- ------ -------- -------- --------- -------- -------- --------- ------ ------ --------- ------ -------\n"; for (const auto &c : crystals.crystals) { - os << fmt::format(" {:8.2f} {:6.2f} {:8.1f} {:8.1f} {:9d} {:8.1f} {:8.1f} {:9.1f} {:6.3f} {:9.3f} {:>6} {:7d}\n", + os << fmt::format(" {:8.2f} {:6.2f} {:8.1f} {:8.1f} {:9d} {:8.1f} {:8.1f} {:9.1f} {:6.3f} {:6.3f} {:9.3f} {:>6} {:7d}\n", c.nx, c.ny, c.x_um, c.y_um, c.image_number, c.major_um, c.minor_um, - c.angle_deg, c.score, c.ice_score, + c.angle_deg, c.score, c.peak_score, c.ice_score, std::isfinite(c.res_A) ? fmt::format("{:.2f}", c.res_A) : std::string("-"), c.n_images); } if (crystals.crystals.empty()) os << " (none)\n"; else - os << " -------- ------ -------- -------- --------- -------- -------- --------- ------ --------- ------ -------\n"; + os << " -------- ------ -------- -------- --------- -------- -------- --------- ------ ------ --------- ------ -------\n"; os << "\n" << BANNER << "\n END OF REPORT\n" << BANNER << "\n"; return os.str(); @@ -162,7 +168,7 @@ std::string RenderRasterJson(const std::string &input_file, const GridScanResult &crystals, const RasterSettings &settings, const RunProvenance &provenance) { - const ScanSummary summary = Summarize(scan, settings.protein_score_threshold); + const ScanSummary summary = Summarize(scan, settings.analysis.protein_score_threshold); nlohmann::json out; out["raster_report_version"] = RASTER_REPORT_VERSION; @@ -189,8 +195,10 @@ std::string RenderRasterJson(const std::string &input_file, out["grid"] = g; nlohmann::json s; - s["protein_score_threshold"] = settings.protein_score_threshold; - s["min_blob_cells"] = settings.min_blob_cells; + s["protein_score_threshold"] = settings.analysis.protein_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; s["beam_size_x_um"] = settings.beam_size_x_um; s["beam_size_y_um"] = settings.beam_size_y_um; s["beam_size_source"] = settings.beam_size_source; @@ -220,6 +228,7 @@ std::string RenderRasterJson(const std::string &input_file, j["minor_um"] = c.minor_um; j["angle_deg"] = c.angle_deg; j["score"] = c.score; + j["peak_score"] = c.peak_score; j["ice_score"] = c.ice_score; // JSON has no NaN, so a resolution that was never measured is null rather than a number. if (std::isfinite(c.res_A)) diff --git a/rugnux/RasterReport.h b/rugnux/RasterReport.h index 5d039ab42..6aa9e8da3 100644 --- a/rugnux/RasterReport.h +++ b/rugnux/RasterReport.h @@ -8,6 +8,7 @@ #include "../common/GridScanSettings.h" #include "../common/Logger.h" #include "../common/ScanResult.h" +#include "../image_analysis/grid_scan_analysis/AnalyzeGridScan.h" #include "../image_analysis/grid_scan_analysis/GridScanResult.h" #include "ResultReport.h" // RunProvenance @@ -21,8 +22,7 @@ // runs apart. The beam size is here because it is the one setting a wrong value silently rotates the // answer with - see AnalyzeGridScan. struct RasterSettings { - float protein_score_threshold = 0; - int64_t min_blob_cells = 0; + GridScanAnalysisParameters analysis; float beam_size_x_um = 0, beam_size_y_um = 0; // Where the beam size came from: "FILE" (incident_beam_size), "COMMAND_LINE" (--beam-size) or // "NOT_STATED" (neither, and the sizes above are zero). diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 46d43828c..41b76a489 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -127,7 +127,9 @@ 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-min-cells Cells a connected patch must have before it is reported as a crystal (default: 3)" << 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; std::cout << std::endl; std::cout << " Calibration (--mode calibration)" << std::endl; @@ -342,7 +344,9 @@ enum { OPT_NO_REFINE_TILT, OPT_BEAM_SIZE, OPT_RASTER_PROTEIN_THRESHOLD, - OPT_RASTER_MIN_CELLS + OPT_RASTER_MIN_CELLS, + OPT_RASTER_DECISIVE_SCORE, + OPT_RASTER_MAX_CRYSTALS }; static option long_options[] = { @@ -370,6 +374,8 @@ 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-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}, {"no-merge", no_argument, nullptr, OPT_NO_MERGE}, {"no-p1-crosscheck", no_argument, nullptr, OPT_NO_P1_CROSSCHECK}, @@ -716,8 +722,8 @@ static int RunRugnux(int argc, char **argv) { // --mode raster. The beam size is empty until it is stated: the file's incident_beam_size is read // once the dataset is open, and --beam-size overrides it there. std::optional beam_size_x_um, beam_size_y_um; // --beam-size - float raster_protein_threshold = PROTEIN_SCORE_THRESHOLD_DEFAULT; // --raster-protein-threshold - int64_t raster_min_cells = MIN_BLOB_CELLS_DEFAULT; // --raster-min-cells + // What counts as a crystal (--raster-*). The defaults are the ones AnalyzeGridScan states. + GridScanAnalysisParameters raster_analysis; bool rotation_indexing = false; bool force_still = false; // --force-still: process a rotation dataset as stills (indexing + scaling) bool two_pass_rotation = true; @@ -1146,11 +1152,20 @@ static int RunRugnux(int argc, char **argv) { break; } case OPT_RASTER_PROTEIN_THRESHOLD: - raster_protein_threshold = parse_number_arg(optarg, "--raster-protein-threshold", - logger, 0.0f, 1.0f); + raster_analysis.protein_score_threshold = + parse_number_arg(optarg, "--raster-protein-threshold", logger, 0.0f, 1.0f); break; case OPT_RASTER_MIN_CELLS: - raster_min_cells = parse_number_arg(optarg, "--raster-min-cells", logger, 1); + raster_analysis.min_blob_cells = + parse_number_arg(optarg, "--raster-min-cells", logger, 1); + break; + case OPT_RASTER_DECISIVE_SCORE: + raster_analysis.decisive_protein_score = + parse_number_arg(optarg, "--raster-decisive-score", logger, 0.0f, 1.0f); + break; + case OPT_RASTER_MAX_CRYSTALS: + raster_analysis.max_crystals = + parse_number_arg(optarg, "--raster-max-crystals", logger, 0); break; case OPT_NO_MERGE: run_scaling = false; @@ -2145,8 +2160,7 @@ static int RunRugnux(int argc, char **argv) { // leaves the size at zero rather than at a guess: a made-up beam would be deconvolved as if // it had been measured, and an anisotropic beam given as square rotates the crystal axis. RasterSettings raster_settings; - raster_settings.protein_score_threshold = raster_protein_threshold; - raster_settings.min_blob_cells = raster_min_cells; + raster_settings.analysis = raster_analysis; if (beam_size_x_um) { raster_settings.beam_size_source = "COMMAND_LINE"; } else { @@ -2212,7 +2226,7 @@ static int RunRugnux(int argc, char **argv) { const GridScanResult crystals = AnalyzeGridScan(observer.scan, *grid, raster_settings.beam_size_x_um, raster_settings.beam_size_y_um, - raster_protein_threshold, raster_min_cells); + raster_analysis); provenance.wall_time_s = std::chrono::duration( std::chrono::steady_clock::now() - invocation_start).count(); @@ -2223,9 +2237,9 @@ static int RunRugnux(int argc, char **argv) { std::cout << fmt::format("Crystals found: {}", crystals.crystals.size()) << std::endl; for (const auto &c : crystals.crystals) std::cout << fmt::format(" cell ({:.1f},{:.1f}) image {} - {:.0f} x {:.0f} um at {:.0f} deg, " - "score {:.2f}, ice {:.2f}, {} cells", + "score {:.2f} (peak {:.2f}), ice {:.2f}, {} cells", c.nx, c.ny, c.image_number, c.major_um, c.minor_um, c.angle_deg, - c.score, c.ice_score, c.n_images) << std::endl; + c.score, c.peak_score, c.ice_score, c.n_images) << std::endl; std::cout << fmt::format("Total wall time: {:.2f} s", provenance.wall_time_s) << std::endl; if (result.cancelled) logger.Warning("Processing was cancelled after {} images", result.images_processed); diff --git a/tests/AnalyzeGridScanTest.cpp b/tests/AnalyzeGridScanTest.cpp index 279be156c..9fe5c0821 100644 --- a/tests/AnalyzeGridScanTest.cpp +++ b/tests/AnalyzeGridScanTest.cpp @@ -50,8 +50,11 @@ TEST_CASE("AnalyzeGridScan finds nothing in an empty scan", "[AnalyzeGridScan]") 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); + SECTION("a two-cell patch is below the minimum blob size and its score is not decisive") { + // Above the admission threshold, so the cells are in a patch at all, but under the score at + // which one cell would be evidence enough on its own. RasterReportTest covers the other + // side of that rule - a patch this small kept because its best cell IS decisive. + auto scan = MakeScan(6, 6, {{2, 2}, {3, 2}}, 0.55f); auto result = AnalyzeGridScan(scan, grid, 5.0f, 5.0f); CHECK(result.crystals.empty()); } diff --git a/tests/RasterReportTest.cpp b/tests/RasterReportTest.cpp index 86e69c5bf..6fbfa3817 100644 --- a/tests/RasterReportTest.cpp +++ b/tests/RasterReportTest.cpp @@ -1,6 +1,8 @@ // SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include + #include #include @@ -38,8 +40,6 @@ namespace { RasterSettings TestSettings() { RasterSettings settings; - settings.protein_score_threshold = PROTEIN_SCORE_THRESHOLD_DEFAULT; - settings.min_blob_cells = MIN_BLOB_CELLS_DEFAULT; settings.beam_size_x_um = 30.0f; settings.beam_size_y_um = 10.0f; settings.beam_size_source = "COMMAND_LINE"; @@ -52,16 +52,14 @@ TEST_CASE("RasterReport_Render", "[Diagnostics]") { const ScanResult scan = TestScan(); const RasterSettings settings = TestSettings(); const GridScanResult crystals = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, - settings.beam_size_y_um, - settings.protein_score_threshold, - settings.min_blob_cells); + settings.beam_size_y_um, settings.analysis); REQUIRE(crystals.crystals.size() == 1); const std::string report = RenderRasterReport("test_master.h5", grid, scan, crystals, settings); // 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= 1\n") != std::string::npos); + CHECK(report.find("RASTER_REPORT_VERSION= 2\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); @@ -71,6 +69,9 @@ 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("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); CHECK(report.find("CRYSTAL_COUNT= 1\n") != std::string::npos); CHECK(report.find("END OF REPORT") != std::string::npos); } @@ -80,9 +81,7 @@ TEST_CASE("RasterReport_Json", "[Diagnostics]") { const ScanResult scan = TestScan(); const RasterSettings settings = TestSettings(); const GridScanResult crystals = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, - settings.beam_size_y_um, - settings.protein_score_threshold, - settings.min_blob_cells); + settings.beam_size_y_um, settings.analysis); const nlohmann::json j = nlohmann::json::parse( RenderRasterJson("test_master.h5", grid, scan, crystals, settings)); @@ -101,6 +100,7 @@ TEST_CASE("RasterReport_Json", "[Diagnostics]") { CHECK(c["image_number"] == crystals.crystals[0].image_number); CHECK(c["res_A"].get() == Catch::Approx(2.0)); CHECK(c["score"].get() == Catch::Approx(0.9).margin(1e-5)); + CHECK(c["peak_score"].get() == Catch::Approx(0.9).margin(1e-5)); } TEST_CASE("RasterReport_NoCrystal", "[Diagnostics]") { @@ -113,9 +113,7 @@ TEST_CASE("RasterReport_NoCrystal", "[Diagnostics]") { const RasterSettings settings = TestSettings(); const GridScanResult crystals = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, - settings.beam_size_y_um, - settings.protein_score_threshold, - settings.min_blob_cells); + settings.beam_size_y_um, settings.analysis); CHECK(crystals.crystals.empty()); const std::string report = RenderRasterReport("test_master.h5", grid, scan, crystals, settings); @@ -127,3 +125,110 @@ TEST_CASE("RasterReport_NoCrystal", "[Diagnostics]") { CHECK(j["crystal_count"] == 0); CHECK(j["crystals"].empty()); } + +// The blob search returns every patch it finds, so a raster with two crystals in it has to come back +// with two rows in score order - and a caller that only wants the best few has to be able to say so. +namespace { + // A wider grid, so several crystals fit in it with empty cells between them. + constexpr int64_t WIDE_NX = 12; + constexpr int64_t WIDE_NY = 5; + + GridScanSettings WideGrid() { + GridScanSettings grid(WIDE_NX, 20.0f, 16.0f, /*snake=*/false, /*vertical=*/false); + grid.ImageNum(WIDE_NX * WIDE_NY); + return grid; + } + + // Each patch is a 2 x 2 block whose top-left cell is (x, y), all of its cells at score. + ScanResult WideScan(const std::vector> &patches) { + ScanResult scan; + for (int64_t i = 0; i < WIDE_NX * WIDE_NY; i++) { + ScanResultElem elem; + elem.number = i; + elem.protein_score = 0.1f; + scan.images.push_back(elem); + } + for (const auto &[x, y, score] : patches) { + for (int64_t dy = 0; dy < 2; dy++) { + for (int64_t dx = 0; dx < 2; dx++) + scan.images[(y + dy) * WIDE_NX + x + dx].protein_score = score; + } + } + return scan; + } +} + +TEST_CASE("RasterReport_TwoCrystals", "[Diagnostics]") { + const GridScanSettings grid = WideGrid(); + const ScanResult scan = WideScan({{1, 1, 0.9f}, {8, 2, 0.7f}}); + RasterSettings settings = TestSettings(); + + const GridScanResult crystals = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, + settings.beam_size_y_um, settings.analysis); + REQUIRE(crystals.crystals.size() == 2); + // Best first, and the two are far enough apart that neither took the other's cells. + CHECK(crystals.crystals[0].score == Catch::Approx(0.9).margin(1e-5)); + CHECK(crystals.crystals[1].score == Catch::Approx(0.7).margin(1e-5)); + CHECK(crystals.crystals[0].n_images == 4); + CHECK(crystals.crystals[1].n_images == 4); + CHECK(crystals.crystals[0].nx < crystals.crystals[1].nx); + + const std::string report = RenderRasterReport("test_master.h5", grid, scan, crystals, settings); + CHECK(report.find("CRYSTAL_COUNT= 2\n") != std::string::npos); + + const nlohmann::json j = nlohmann::json::parse( + RenderRasterJson("test_master.h5", grid, scan, crystals, settings)); + REQUIRE(j["crystals"].size() == 2); + CHECK(j["crystals"][0]["score"].get() > j["crystals"][1]["score"].get()); +} + +TEST_CASE("RasterReport_ThreeCrystalsAndTheCap", "[Diagnostics]") { + const GridScanSettings grid = WideGrid(); + const ScanResult scan = WideScan({{0, 0, 0.95f}, {5, 2, 0.75f}, {10, 0, 0.85f}}); + RasterSettings settings = TestSettings(); + + const GridScanResult all = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, + settings.beam_size_y_um, settings.analysis); + REQUIRE(all.crystals.size() == 3); + CHECK(all.crystals[0].score == Catch::Approx(0.95).margin(1e-5)); + CHECK(all.crystals[1].score == Catch::Approx(0.85).margin(1e-5)); + CHECK(all.crystals[2].score == Catch::Approx(0.75).margin(1e-5)); + + // The cap keeps the strongest, not the first found: the 0.85 patch is the last one in grid order. + settings.analysis.max_crystals = 2; + const GridScanResult capped = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, + settings.beam_size_y_um, settings.analysis); + REQUIRE(capped.crystals.size() == 2); + CHECK(capped.crystals[0].score == Catch::Approx(0.95).margin(1e-5)); + CHECK(capped.crystals[1].score == Catch::Approx(0.85).margin(1e-5)); + + const std::string report = RenderRasterReport("test_master.h5", grid, scan, capped, settings); + CHECK(report.find("MAX_CRYSTALS= 2\n") != std::string::npos); + CHECK(report.find("CRYSTAL_COUNT= 2\n") != std::string::npos); +} + +TEST_CASE("RasterReport_OneDecisiveCellIsACrystal", "[Diagnostics]") { + // A patch below MIN_BLOB_CELLS is kept when its best cell is decisive and dropped when it is not. + // Both patches here are one cell, so the only thing separating them is the evidence in them. + const GridScanSettings grid = WideGrid(); + ScanResult scan = WideScan({}); + scan.images[1 * WIDE_NX + 2].protein_score = 0.95f; // decisive + scan.images[3 * WIDE_NX + 9].protein_score = 0.55f; // above the admission threshold, not decisive + + RasterSettings settings = TestSettings(); + const GridScanResult crystals = AnalyzeGridScan(scan, grid, settings.beam_size_x_um, + settings.beam_size_y_um, settings.analysis); + REQUIRE(crystals.crystals.size() == 1); + CHECK(crystals.crystals[0].n_images == 1); + CHECK(crystals.crystals[0].peak_score == Catch::Approx(0.95).margin(1e-5)); + + // Raise the bar past the one cell that passed and nothing is left; a patch that meets the cell + // minimum is unaffected by the bar, which is what makes this an OR and not a second gate. + settings.analysis.decisive_protein_score = 0.99f; + CHECK(AnalyzeGridScan(scan, grid, settings.beam_size_x_um, settings.beam_size_y_um, + settings.analysis).crystals.empty()); + + const ScanResult big = WideScan({{1, 1, 0.55f}}); + CHECK(AnalyzeGridScan(big, grid, settings.beam_size_x_um, settings.beam_size_y_um, + settings.analysis).crystals.size() == 1); +}