Files
Jungfraujoch/image_analysis/grid_scan_analysis/AnalyzeGridScan.h
T
leonarski_fandClaude Opus 5 273591af11 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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
2026-09-08 07:34:36 +02:00

72 lines
4.2 KiB
C++

// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
// SPDX-License-Identifier: GPL-3.0-only
#pragma once
#include "../../common/GridScanSettings.h"
#include "../../common/ScanResult.h"
#include "GridScanResult.h"
// Finds the crystals in a completed grid scan: the per-image protein score is scattered back onto
// the raster, the map is thresholded and labelled, and each blob is reported as one crystal.
// Pure - no I/O, no state, no detector.
//
// The sizes reported are MEASURED, and the beam is still in them: what the raster sees is the
// crystal convolved with the beam, and nothing here takes the beam back out. That is deliberate -
// the beam size is in the file as incident_beam_size, so a downstream consumer can do the
// deconvolution itself, reproducibly and reversibly, and is not stuck with ours.
//
// If it does: removing an anisotropic beam is a subtraction of the two 2x2 covariance matrices
// followed by re-diagonalising the difference, NOT a per-axis quadrature removal of beam_x from
// major_um and beam_y from minor_um. Per-axis is silently wrong the moment the crystal is not
// aligned with the grid axes - a needle at 45 deg has both beam widths mixed into both of its own
// axes - and a needle at an arbitrary angle is exactly the case this whole design exists for.
// beam_size_x_um/beam_size_y_um are only copied into the result, so it says what the sizes contain.
// 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 - unless the diffraction in it is decisive, see
// 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
// 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;
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;
};
GridScanResult AnalyzeGridScan(const ScanResult &scan,
const GridScanSettings &grid,
float beam_size_x_um,
float beam_size_y_um,
const GridScanAnalysisParameters &params = {});