diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 9c44be2b2..a330571d7 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,6 +3,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 --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 24daa9691..1a484bc73 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -216,6 +216,15 @@ Mode — `--mode ` (default `mx`): | `azint` | Only azimuthal integration (no spot finding/indexing); writes `_process.h5` | | `scale` | Only re-scale/merge the already-integrated reflections in the input `_process.h5` (no re-integration) | | `calibration` | Determine the detector geometry from powder rings; writes `.poni` and `.json` | +| `raster` | Score every image of a stored grid scan and report the crystals in it; writes `_raster_report.txt` and `_raster.json` | + +Raster (`--mode raster`): + +| Option | Description | +| --- | --- | +| `--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`) | Calibration (`--mode calibration`): diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp index 5bce4d06b..6fdbce02a 100644 --- a/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.cpp @@ -10,21 +10,15 @@ #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) { + float beam_size_y_um, + float protein_score_threshold, + int64_t min_blob_cells) { const int64_t nx = grid.GetGridSizeX_step(); const int64_t ny = grid.GetGridSizeY_step(); const float step_x = fabsf(grid.GetGridStepX_um()); @@ -49,9 +43,9 @@ 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] > protein_score_threshold) ? 1 : 0; - const std::vector label = LabelConnectedComponents(above, nx, ny, MIN_BLOB_CELLS); + 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; @@ -65,7 +59,7 @@ GridScanResult AnalyzeGridScan(const ScanResult &scan, cell.push_back(i); } if (cell.empty()) - continue; // label fell below MIN_BLOB_CELLS and was erased + continue; // label fell below min_blob_cells and was erased const auto n = static_cast(cell.size()); diff --git a/image_analysis/grid_scan_analysis/AnalyzeGridScan.h b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h index 6c52acf39..d34d3a258 100644 --- a/image_analysis/grid_scan_analysis/AnalyzeGridScan.h +++ b/image_analysis/grid_scan_analysis/AnalyzeGridScan.h @@ -22,7 +22,22 @@ // 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. +constexpr int64_t MIN_BLOB_CELLS_DEFAULT = 3; + GridScanResult AnalyzeGridScan(const ScanResult &scan, const GridScanSettings &grid, float beam_size_x_um, - float beam_size_y_um); + float beam_size_y_um, + float protein_score_threshold = PROTEIN_SCORE_THRESHOLD_DEFAULT, + int64_t min_blob_cells = MIN_BLOB_CELLS_DEFAULT); diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index a8d0e7c24..fd2015581 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -10,6 +10,8 @@ ADD_LIBRARY(Rugnux STATIC RugnuxCommandLine.h ModelValidation.cpp ModelValidation.h + RasterReport.cpp + RasterReport.h RigidBodyRefine.cpp RigidBodyRefine.h SigmaA.cpp @@ -22,12 +24,13 @@ ADD_LIBRARY(Rugnux STATIC WriteModel.h ) -TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochWriter gemmi) +TARGET_LINK_LIBRARIES(Rugnux JFJochReader JFJochImageAnalysis JFJochGridScanAnalysis JFJochWriter gemmi) # rugnux is the single offline analysis CLI; --mode picks what it does: mx (the full pipeline, the -# default), azint, scale (re-scale/merge stored reflections) or calibration (powder-ring geometry). +# default), azint, scale (re-scale/merge stored reflections), calibration (powder-ring geometry) or +# raster (the crystals in a stored grid scan). ADD_EXECUTABLE(rugnux rugnux_cli.cpp) -TARGET_LINK_LIBRARIES(rugnux Rugnux JFJochReader JFJochImageAnalysis JFJochWriter) +TARGET_LINK_LIBRARIES(rugnux Rugnux JFJochReader JFJochImageAnalysis JFJochGridScanAnalysis JFJochWriter) INSTALL(TARGETS rugnux RUNTIME COMPONENT rugnux) # libcufft_static.a carries a relocatable-device-code object (separate_callback.o), so an executable diff --git a/rugnux/RasterReport.cpp b/rugnux/RasterReport.cpp new file mode 100644 index 000000000..09fdebc26 --- /dev/null +++ b/rugnux/RasterReport.cpp @@ -0,0 +1,266 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include +#include + +#include +#include + +#include "../common/GitInfo.h" +#include "../common/time_utc.h" + +#include "RasterReport.h" + +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; + + const char *BANNER = " ******************************************************************************"; + + void Section(std::ostream &os, const std::string &title) { + os << "\n" << BANNER << "\n " << title << "\n" << BANNER << "\n\n"; + } + + template void Key(std::ostream &os, const char *key, const T &value) { + os << key << "= " << value << "\n"; + } + + // What the per-image scores add up to over the whole raster. Both renderers say the same thing, + // so both read it from here. + struct ScanSummary { + int64_t scored = 0; // images that came back with a protein score at all + int64_t above = 0; // ... of which above the threshold the blob search used + int64_t ice_above = 0; // ... whose ICE score is above that same threshold + float mean_protein = 0; + float mean_ice = 0; + }; + + ScanSummary Summarize(const ScanResult &scan, float threshold) { + ScanSummary s; + float sum_protein = 0, sum_ice = 0; + for (const auto &e : scan.images) { + if (!e.protein_score) + continue; + s.scored++; + sum_protein += *e.protein_score; + if (*e.protein_score > threshold) + s.above++; + if (e.ice_score) { + sum_ice += *e.ice_score; + if (*e.ice_score > threshold) + s.ice_above++; + } + } + if (s.scored > 0) { + s.mean_protein = sum_protein / static_cast(s.scored); + s.mean_ice = sum_ice / static_cast(s.scored); + } + return s; + } + + float Fraction(int64_t part, int64_t whole) { + return whole > 0 ? static_cast(part) / static_cast(whole) : 0.0f; + } +} + +std::string RenderRasterReport(const std::string &input_file, + const GridScanSettings &grid, + const ScanResult &scan, + const GridScanResult &crystals, + const RasterSettings &settings, + const RunProvenance &provenance) { + std::ostringstream os; + const ScanSummary summary = Summarize(scan, settings.protein_score_threshold); + + os << BANNER << "\n" + << " RUGNUX RASTER REPORT\n" + << BANNER << "\n\n" + << " What this grid scan found: where the crystals are, how large they are and how well they\n" + << " diffract. The `KEY= value` lines and the table below are a stable interface - a script\n" + << " greps them, and RASTER_REPORT_VERSION says when that interface last changed. The same\n" + << " content, typed, is in the _raster.json beside this file.\n\n"; + + Key(os, "RASTER_REPORT_VERSION", RASTER_REPORT_VERSION); + Key(os, "RUGNUX_VERSION", jfjoch_version()); + if (!jfjoch_git_sha1().empty()) + Key(os, "RUGNUX_GIT", jfjoch_git_sha1().substr(0, 6) + " " + jfjoch_git_date()); + Key(os, "DATE", time_UTC(std::chrono::system_clock::now())); + Key(os, "INPUT_FILE", input_file); + if (!provenance.command_line.empty()) + Key(os, "COMMAND_LINE", provenance.command_line); + if (provenance.wall_time_s > 0.0) + Key(os, "WALL_TIME", fmt::format("{:.2f}", provenance.wall_time_s)); + + // ------------------------------------------------------------------------------ 1. GRID + Section(os, "1. GRID"); + os << " The raster as it was collected. Positions everywhere below are in the DISPLAY grid -\n" + << " column 0 the lowest x, row 0 the lowest y - whichever way the stage moved, so a snake\n" + << " scan and a raster scan of the same area name the same cell the same way.\n\n"; + Key(os, "GRID_SIZE", fmt::format("{} {}", grid.GetGridSizeX_step(), grid.GetGridSizeY_step())); + Key(os, "GRID_STEP_UM", fmt::format("{:.2f} {:.2f}", std::fabs(grid.GetGridStepX_um()), + std::fabs(grid.GetGridStepY_um()))); + Key(os, "GRID_EXTENT_UM", fmt::format("{:.1f} {:.1f}", grid.GetGridSizeX_um(), grid.GetGridSizeY_um())); + Key(os, "GRID_SNAKE", grid.IsSnakeScan() ? "TRUE" : "FALSE"); + Key(os, "GRID_VERTICAL", grid.IsVerticalScan() ? "TRUE" : "FALSE"); + Key(os, "GRID_N_ELEM", grid.GetNElem()); + Key(os, "IMAGES_PROCESSED", scan.images.size()); + Key(os, "BEAM_SIZE_UM", fmt::format("{:.2f} {:.2f}", settings.beam_size_x_um, settings.beam_size_y_um)); + Key(os, "BEAM_SIZE_SOURCE", settings.beam_size_source); + + // -------------------------------------------------------------------- 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. 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); + 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))); + Key(os, "ICE_ABOVE_THRESHOLD", summary.ice_above); + Key(os, "ICE_FRACTION", fmt::format("{:.4f}", Fraction(summary.ice_above, summary.scored))); + Key(os, "MEAN_PROTEIN_SCORE", fmt::format("{:.4f}", summary.mean_protein)); + Key(os, "MEAN_ICE_SCORE", fmt::format("{:.4f}", summary.mean_ice)); + + // ------------------------------------------------------------------------- 3. CRYSTALS + Section(os, "3. CRYSTALS"); + os << " One row per crystal, best first. NX/NY and X_UM/Y_UM are the centre, the micrometres\n" + << " measured from the centre of cell (0,0) along the grid axes; IMAGE is the nearest image\n" + << " that was actually collected, so a DAQ can drive back to it. MAJOR_UM x MINOR_UM is the\n" + << " 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"; + 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"; + 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", + 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, + 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" << BANNER << "\n END OF REPORT\n" << BANNER << "\n"; + return os.str(); +} + +std::string RenderRasterJson(const std::string &input_file, + const GridScanSettings &grid, + const ScanResult &scan, + const GridScanResult &crystals, + const RasterSettings &settings, + const RunProvenance &provenance) { + const ScanSummary summary = Summarize(scan, settings.protein_score_threshold); + + nlohmann::json out; + out["raster_report_version"] = RASTER_REPORT_VERSION; + out["rugnux_version"] = jfjoch_version(); + if (!jfjoch_git_sha1().empty()) + out["rugnux_git"] = jfjoch_git_sha1().substr(0, 6); + out["date"] = time_UTC(std::chrono::system_clock::now()); + out["input_file"] = input_file; + if (!provenance.command_line.empty()) + out["command_line"] = provenance.command_line; + if (provenance.wall_time_s > 0.0) + out["wall_time_s"] = provenance.wall_time_s; + + nlohmann::json g; + g["size_x"] = grid.GetGridSizeX_step(); + g["size_y"] = grid.GetGridSizeY_step(); + g["step_x_um"] = std::fabs(grid.GetGridStepX_um()); + g["step_y_um"] = std::fabs(grid.GetGridStepY_um()); + g["extent_x_um"] = grid.GetGridSizeX_um(); + g["extent_y_um"] = grid.GetGridSizeY_um(); + g["snake"] = grid.IsSnakeScan(); + g["vertical"] = grid.IsVerticalScan(); + g["n_elem"] = grid.GetNElem(); + out["grid"] = g; + + nlohmann::json s; + s["protein_score_threshold"] = settings.protein_score_threshold; + s["min_blob_cells"] = settings.min_blob_cells; + 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; + out["settings"] = s; + + nlohmann::json sc; + sc["images_processed"] = scan.images.size(); + sc["images_scored"] = summary.scored; + sc["images_above_threshold"] = summary.above; + sc["fraction_above_threshold"] = Fraction(summary.above, summary.scored); + sc["ice_above_threshold"] = summary.ice_above; + sc["ice_fraction"] = Fraction(summary.ice_above, summary.scored); + sc["mean_protein_score"] = summary.mean_protein; + sc["mean_ice_score"] = summary.mean_ice; + out["scores"] = sc; + + out["crystal_count"] = crystals.crystals.size(); + nlohmann::json list = nlohmann::json::array(); + for (const auto &c : crystals.crystals) { + nlohmann::json j; + j["nx"] = c.nx; + j["ny"] = c.ny; + j["x_um"] = c.x_um; + j["y_um"] = c.y_um; + j["image_number"] = c.image_number; + j["major_um"] = c.major_um; + j["minor_um"] = c.minor_um; + j["angle_deg"] = c.angle_deg; + j["score"] = c.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)) + j["res_A"] = c.res_A; + else + j["res_A"] = nullptr; + j["n_images"] = c.n_images; + list.push_back(j); + } + out["crystals"] = list; + + return out.dump(1) + "\n"; +} + +void WriteRasterReport(const std::string &output_prefix, + const std::string &input_file, + const GridScanSettings &grid, + const ScanResult &scan, + const GridScanResult &crystals, + const RasterSettings &settings, + Logger &logger, + const RunProvenance &provenance) { + if (output_prefix.empty()) + return; // "analyse the raster, persist nothing" + + // As for the results report: a run that found the crystals must not fail because a side file + // could not be written. Report the failure and carry on. + const std::string txt_path = output_prefix + "_raster_report.txt"; + const std::string json_path = output_prefix + "_raster.json"; + try { + std::ofstream txt(txt_path); + txt.exceptions(std::ios::failbit | std::ios::badbit); + txt << RenderRasterReport(input_file, grid, scan, crystals, settings, provenance); + } catch (const std::exception &e) { + logger.Warning("Could not write the raster report {}: {}", txt_path, e.what()); + } + try { + std::ofstream json(json_path); + json.exceptions(std::ios::failbit | std::ios::badbit); + json << RenderRasterJson(input_file, grid, scan, crystals, settings, provenance); + } catch (const std::exception &e) { + logger.Warning("Could not write the raster JSON {}: {}", json_path, e.what()); + } +} diff --git a/rugnux/RasterReport.h b/rugnux/RasterReport.h new file mode 100644 index 000000000..5d039ab42 --- /dev/null +++ b/rugnux/RasterReport.h @@ -0,0 +1,58 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +#include "../common/GridScanSettings.h" +#include "../common/Logger.h" +#include "../common/ScanResult.h" +#include "../image_analysis/grid_scan_analysis/GridScanResult.h" +#include "ResultReport.h" // RunProvenance + +// _raster_report.txt and _raster.json - what a grid scan found, written by +// rugnux --mode raster. Same idiom as the results report next door: prose a person reads top to +// bottom, `KEY= value` assignment lines and one fixed-width table a script greps without parsing +// prose, and a RASTER_REPORT_VERSION that is bumped if any of that moves. The JSON carries the same +// content typed, for a battery that aggregates many rasters at once. + +// What the run analysed with, so the report says what produced its numbers and a sweep can tell two +// 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; + 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). + std::string beam_size_source; +}; + +// Renders the human-readable report. Exposed for testing; the CLI writes it with WriteRasterReport. +std::string RenderRasterReport(const std::string &input_file, + const GridScanSettings &grid, + const ScanResult &scan, + const GridScanResult &crystals, + const RasterSettings &settings, + const RunProvenance &provenance = {}); + +// Renders the same content as JSON. Exposed for testing. +std::string RenderRasterJson(const std::string &input_file, + const GridScanSettings &grid, + const ScanResult &scan, + const GridScanResult &crystals, + const RasterSettings &settings, + const RunProvenance &provenance = {}); + +// Writes _raster_report.txt and _raster.json. Does nothing when the +// prefix is empty. Never throws: a run that found the crystals must not be lost because a side file +// could not be written, so a failure is logged and swallowed. +void WriteRasterReport(const std::string &output_prefix, + const std::string &input_file, + const GridScanSettings &grid, + const ScanResult &scan, + const GridScanResult &crystals, + const RasterSettings &settings, + Logger &logger, + const RunProvenance &provenance = {}); diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index d672d9a80..6c5458a1a 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -41,10 +42,12 @@ #include "../image_analysis/scale_merge/AnisotropyAnalysis.h" #include "../image_analysis/scale_merge/TwinningAnalysis.h" #include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "../image_analysis/grid_scan_analysis/AnalyzeGridScan.h" #include "Rugnux.h" #include "RugnuxDefaults.h" #include "ModelValidation.h" #include "WriteModel.h" +#include "RasterReport.h" #include "ResultReport.h" // Spots kept per image (the strongest ones) and handed to indexing. Offline reprocessing is not @@ -67,7 +70,7 @@ constexpr double RELATIVE_B_DEFAULT_DEG = 10.0; // What the run does (--mode). Each of these was once a flag of its own; one option keeps them // mutually exclusive by construction. -enum class RugnuxMode { MX, Azint, Scale, Calibration }; +enum class RugnuxMode { MX, Azint, Scale, Calibration, Raster }; // Azimuthal sectors used by --mode calibration --calibration rings when --azim-phi-bins was not given. // The default of 1 is a plain radial profile, in which a ring has been averaged over every direction so @@ -77,6 +80,26 @@ enum class RugnuxMode { MX, Azint, Scale, Calibration }; // buys nothing and runs into the 65534-bin limit on q x azimuth, so 32 is the default. constexpr int CALIBRATION_AZIM_BINS_DEFAULT = 32; +// --mode raster: what a grid scan needs off the per-image loop. AnalyzeGridScan reads elem.number as +// the image's place in the grid, so it is the file's OWN image number that goes there - the loop's +// ordinal counts only the images -s/-e/--stride selected and means nothing to a grid. +class RasterObserver : public RugnuxObserver { + std::mutex m; +public: + ScanResult scan; + + void OnImageProcessed(const DataMessage &msg) override { + std::lock_guard lock(m); + ScanResultElem elem; + elem.number = msg.original_number.value_or(msg.number); + elem.protein_score = msg.protein_score; + elem.ice_score = msg.ice_score; + elem.res = msg.resolution_estimate; + elem.spot_count = msg.spot_count; + scan.images.push_back(elem); + } +}; + void print_usage() { std::cout << "Usage rugnux {} " << std::endl; std::cout << "Options:" << std::endl; @@ -96,6 +119,13 @@ void print_usage() { std::cout << " azint Only azimuthal integration (no spot finding/indexing); writes _process.h5" << std::endl; std::cout << " scale Only re-scale/merge the already-integrated reflections in (no re-integration)" << std::endl; std::cout << " calibration Determine the detector geometry from powder rings; writes .poni for pyFAI and .json, whose dataset_settings member is a jfjoch_broker dataset_settings body as it stands" << std::endl; + std::cout << " raster Score every image of a stored grid scan and report the crystals in it; writes _raster_report.txt and _raster.json" << std::endl; + std::cout << std::endl; + + 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 << std::endl; std::cout << " Calibration (--mode calibration)" << std::endl; @@ -301,7 +331,10 @@ enum { OPT_DETECTOR_QUARTER_TURNS, OPT_FFT_MIN_UNIT_CELL, OPT_POLARIZATION, - OPT_NO_REFINE_TILT + OPT_NO_REFINE_TILT, + OPT_BEAM_SIZE, + OPT_RASTER_PROTEIN_THRESHOLD, + OPT_RASTER_MIN_CELLS }; static option long_options[] = { @@ -326,6 +359,9 @@ static option long_options[] = { {"mode", required_argument, nullptr, OPT_MODE}, {"calibrant", required_argument, nullptr, OPT_CALIBRANT}, {"calibration", required_argument, nullptr, OPT_CALIBRATION}, + {"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}, {"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}, @@ -666,6 +702,11 @@ static int RunRugnux(int argc, char **argv) { std::string calibrant = "lab6"; // --calibrant CalibrationMethod calibration_method = CalibrationMethod::Rings; // --calibration bool calibration_refine_tilt = true; // --no-refine-tilt + // --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 bool rotation_indexing = false; bool force_still = false; // --force-still: process a rotation dataset as stills (indexing + scaling) bool two_pass_rotation = true; @@ -1056,8 +1097,9 @@ static int RunRugnux(int argc, char **argv) { else if (m == "azint") mode = RugnuxMode::Azint; else if (m == "scale") mode = RugnuxMode::Scale; else if (m == "calibration") mode = RugnuxMode::Calibration; + else if (m == "raster") mode = RugnuxMode::Raster; else { - logger.Error("Invalid --mode: {} (expected mx|azint|scale|calibration)", m); + logger.Error("Invalid --mode: {} (expected mx|azint|scale|calibration|raster)", m); return 1; } break; @@ -1080,6 +1122,23 @@ static int RunRugnux(int argc, char **argv) { case OPT_NO_REFINE_TILT: calibration_refine_tilt = false; break; + case OPT_BEAM_SIZE: { + // "x" is a square beam, "x,y" a rectangular one. Nothing else: a beam has two sizes. + const std::string arg = optarg ? optarg : ""; + const size_t comma = arg.find(','); + beam_size_x_um = parse_number_arg(arg.substr(0, comma).c_str(), "--beam-size", logger, 0.0f); + beam_size_y_um = comma == std::string::npos + ? beam_size_x_um + : parse_number_arg(arg.substr(comma + 1).c_str(), "--beam-size", logger, 0.0f); + break; + } + case OPT_RASTER_PROTEIN_THRESHOLD: + raster_protein_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); + break; case OPT_NO_MERGE: run_scaling = false; break; @@ -1997,6 +2056,111 @@ static int RunRugnux(int argc, char **argv) { return 0; } + // --mode raster: score every image of a stored grid scan and report the crystals in it. Spot + // finding runs, because the per-cell resolution comes from it; indexing does not, because a + // raster cell is one short still and the question a raster answers is where the crystal is, + // not what its lattice is. + if (mode == RugnuxMode::Raster) { + ApplyRugnuxExperimentDefaults(experiment); + if (polarization_factor) experiment.PolarizationFactor(polarization_factor.value()); + + const auto grid = experiment.GetGridScan(); + if (!grid) { + logger.Error("--mode raster needs a grid scan and {} does not describe one " + "(no /entry/sample/grid_scan)", input_file); + return 1; + } + + // The beam is not taken out of the reported crystal sizes, but it is reported alongside them + // so a consumer can do it. --beam-size wins over the file, and a file that states nothing + // 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; + if (beam_size_x_um) { + raster_settings.beam_size_source = "COMMAND_LINE"; + } else { + beam_size_x_um = experiment.GetBeamSizeX_um(); + beam_size_y_um = experiment.GetBeamSizeY_um(); + raster_settings.beam_size_source = beam_size_x_um ? "FILE" : "NOT_STATED"; + } + raster_settings.beam_size_x_um = beam_size_x_um.value_or(0.0f); + raster_settings.beam_size_y_um = beam_size_y_um.value_or(0.0f); + + ProcessConfig config; + config.mode = ProcessMode::FullAnalysis; + config.start_image = start_image; + config.end_image = end_image; + config.stride = image_stride; + config.nthreads = nthreads; + config.nthreads_auto = nthreads_auto; + config.output_prefix = output_prefix; + config.detect_beam_stop = detect_beam_stop; + config.estimate_beam_center = estimate_beam_center; + config.beam_center_check = beam_center_check; + config.fit_spindle = fit_spindle; + config.run_scaling = false; + config.write_merged = false; + // The two raster files below are the output of this mode; the per-image _process.h5 is + // written only when it is asked for. + config.write_process_h5 = write_process_h5_flag; + + config.spot_finding.enable = true; + config.spot_finding.indexing = false; + config.spot_finding.signal_to_noise_threshold = sigma_spot_finding; + config.spot_finding.photon_count_threshold = photon_count_threshold_spot_finding; + config.spot_finding.min_pix_per_spot = min_pix_per_spot.value_or(2); + config.spot_finding.adaptive_threshold = adaptive_spots.value_or(true); + config.spot_finding.high_resolution_limit = d_min_spot_finding; + config.spot_finding.false_pixels_per_frame = false_pixels_per_frame; + config.spot_finding.low_resolution_limit = d_max_spot_finding; + experiment.MaxSpotCount(max_spot_count_override.value_or(RUGNUX_MAX_SPOT_COUNT)); + + logger.Info("Grid scan: {} x {} cells of {:.1f} x {:.1f} um, {} scan, beam {:.1f} x {:.1f} um ({})", + grid->GetGridSizeX_step(), grid->GetGridSizeY_step(), + std::fabs(grid->GetGridStepX_um()), std::fabs(grid->GetGridStepY_um()), + grid->IsSnakeScan() ? "snake" : "raster", + raster_settings.beam_size_x_um, raster_settings.beam_size_y_um, + raster_settings.beam_size_source); + + RasterObserver observer; + Rugnux process(reader, experiment, *dataset->pixel_mask, config); + g_active_process = &process; + std::signal(SIGINT, handle_sigint); + + ProcessResult result; + try { + result = process.Run(&observer); + } catch (const std::exception &e) { + logger.Error("Processing failed: {}", e.what()); + exit(EXIT_FAILURE); + } + g_active_process = nullptr; + + 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); + + provenance.wall_time_s = std::chrono::duration( + std::chrono::steady_clock::now() - invocation_start).count(); + WriteRasterReport(output_prefix, input_file, *grid, observer.scan, crystals, + raster_settings, logger, provenance); + + std::cout << fmt::format("Images processed: {}", result.images_processed) << std::endl; + 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", + 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; + 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); + return 0; + } + // --mode calibration: fit the detector geometry to the calibrant's powder rings, measured over // every selected image, and write it as a PONI file. if (mode == RugnuxMode::Calibration) { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 2f37be496..24dc297cc 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -53,6 +53,7 @@ ADD_EXECUTABLE(jfjoch_test ShadowFinderTest.cpp ParallelForTest.cpp RugnuxTest.cpp + RasterReportTest.cpp ResultReportTest.cpp RugnuxLargeTest.cpp TestData.h diff --git a/tests/RasterReportTest.cpp b/tests/RasterReportTest.cpp new file mode 100644 index 000000000..86e69c5bf --- /dev/null +++ b/tests/RasterReportTest.cpp @@ -0,0 +1,129 @@ +// 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" +#include "../rugnux/RasterReport.h" + +// A 5 x 4 raster with one 2 x 2 patch of diffracting cells in it - the smallest scan that has a +// crystal to report and cells around it that are not one. +namespace { + constexpr int64_t GRID_NX = 5; + constexpr int64_t GRID_NY = 4; + + GridScanSettings TestGrid() { + GridScanSettings grid(GRID_NX, 20.0f, 16.0f, /*snake=*/false, /*vertical=*/false); + grid.ImageNum(GRID_NX * GRID_NY); + return grid; + } + + ScanResult TestScan() { + ScanResult scan; + for (int64_t i = 0; i < GRID_NX * GRID_NY; i++) { + const bool crystal = (i == 6 || i == 7 || i == 11 || i == 12); + ScanResultElem elem; + elem.number = i; + elem.protein_score = crystal ? 0.9f : 0.1f; + elem.ice_score = 0.05f; + if (crystal) + elem.res = 2.0f; + scan.images.push_back(elem); + } + return scan; + } + + 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"; + return settings; + } +} + +TEST_CASE("RasterReport_Render", "[Diagnostics]") { + const GridScanSettings grid = TestGrid(); + 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); + 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("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); + CHECK(report.find("GRID_N_ELEM= 20\n") != std::string::npos); + CHECK(report.find("BEAM_SIZE_UM= 30.00 10.00\n") != std::string::npos); + CHECK(report.find("BEAM_SIZE_SOURCE= COMMAND_LINE\n") != std::string::npos); + 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("CRYSTAL_COUNT= 1\n") != std::string::npos); + CHECK(report.find("END OF REPORT") != std::string::npos); +} + +TEST_CASE("RasterReport_Json", "[Diagnostics]") { + const GridScanSettings grid = TestGrid(); + 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); + + const nlohmann::json j = nlohmann::json::parse( + RenderRasterJson("test_master.h5", grid, scan, crystals, settings)); + + CHECK(j["grid"]["size_x"] == GRID_NX); + CHECK(j["grid"]["size_y"] == GRID_NY); + CHECK(j["scores"]["images_above_threshold"] == 4); + CHECK(j["settings"]["beam_size_x_um"].get() == Catch::Approx(30.0)); + CHECK(j["settings"]["beam_size_y_um"].get() == Catch::Approx(10.0)); + REQUIRE(j["crystal_count"] == 1); + REQUIRE(j["crystals"].size() == 1); + + // The crystal the JSON carries is the one the report's table carries - same numbers, typed. + const auto &c = j["crystals"][0]; + CHECK(c["n_images"] == 4); + 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)); +} + +TEST_CASE("RasterReport_NoCrystal", "[Diagnostics]") { + // A raster over an empty loop: the report still has to say so, with a table that is empty rather + // than absent - a battery that greps CRYSTAL_COUNT must find it on every run. + const GridScanSettings grid = TestGrid(); + ScanResult scan = TestScan(); + for (auto &elem : scan.images) + elem.protein_score = 0.1f; + + 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); + CHECK(crystals.crystals.empty()); + + const std::string report = RenderRasterReport("test_master.h5", grid, scan, crystals, settings); + CHECK(report.find("CRYSTAL_COUNT= 0\n") != std::string::npos); + CHECK(report.find("IMAGES_ABOVE_THRESHOLD= 0\n") != std::string::npos); + + const nlohmann::json j = nlohmann::json::parse( + RenderRasterJson("test_master.h5", grid, scan, crystals, settings)); + CHECK(j["crystal_count"] == 0); + CHECK(j["crystals"].empty()); +}