diff --git a/rugnux/HotPixels.cpp b/rugnux/HotPixels.cpp index 88fdadc9e..c95765cfc 100644 --- a/rugnux/HotPixels.cpp +++ b/rugnux/HotPixels.cpp @@ -335,16 +335,13 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp } }); - // Masked: a persistent pixel standing alone - on its own or as one of a pair; a larger patch is - // a feature of the scattering, not a defect - whose mean is STRONG_RATIO times its ring's and - // whose excess is above the Poisson bound. A weaker one cannot make an outlier. - const auto persistent_neighbours = [&](size_t x, size_t y) { - int count = 0; - for (size_t yy = (y > 0 ? y - 1 : 0); yy <= std::min(height - 1, y + 1); yy++) - for (size_t xx = (x > 0 ? x - 1 : 0); xx <= std::min(width - 1, x + 1); xx++) - if ((xx != x || yy != y) && persistent[yy * width + xx]) count++; - return count; - }; + // Masked: a persistent pixel whose mean is STRONG_RATIO times its ring's and whose excess is above + // the Poisson bound - a weaker one cannot make an outlier. Alone or in a patch: a patch lit on frame + // after frame is stationary in the lab and so no reflection of a rotating crystal, whatever made it, + // and every reflection that crosses it is integrated on top of it. Measured: a ring of about twenty + // pixels at 1e5-3e5 counts on every frame (a dead centre in the middle) merged the reflections + // crossing it to 90 times the strongest real one; the space-group search, which merges without the + // equivalents an outlier test needs, read the crystal's four-fold off them as unrelated. ParallelChunks(static_cast(height), nthreads, [&](int y0, int y1) { for (size_t y = y0; y < static_cast(y1); y++) for (size_t x = 0; x < width; x++) { @@ -355,17 +352,9 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp continue; } if (!persistent[i]) continue; - const int neighbours = persistent_neighbours(x, y); - bool isolated = neighbours == 0; - if (neighbours == 1) - for (size_t yy = (y > 0 ? y - 1 : 0); yy <= std::min(height - 1, y + 1); yy++) - for (size_t xx = (x > 0 ? x - 1 : 0); xx <= std::min(width - 1, x + 1); xx++) - if ((xx != x || yy != y) && persistent[yy * width + xx]) - isolated = persistent_neighbours(xx, yy) == 1; const double mean_level = static_cast(sum_level(i)) / n_valid(i); const double excess = static_cast(sum_value[i] - sum_level(i)) / n_valid(i); - if (isolated - && static_cast(sum_value[i]) >= STRONG_RATIO * static_cast(sum_level(i)) + if (static_cast(sum_value[i]) >= STRONG_RATIO * static_cast(sum_level(i)) && excess >= LIT_NSIGMA * std::sqrt(std::max(mean_level, 0.0)) + LIT_OFFSET) ret.mask[i] = 1; } diff --git a/rugnux/HotPixels.h b/rugnux/HotPixels.h index b837ba32e..a27b7dd0d 100644 --- a/rugnux/HotPixels.h +++ b/rugnux/HotPixels.h @@ -27,9 +27,10 @@ // frames (measured, on the ring's pixels that are not candidates), and the number of lit frames // that fewer than 0.01 of all the detector's pixels would reach by chance is the binomial bound. // -// A persistent pixel is masked only if it stands alone (by itself or in a pair - a larger patch is a -// feature of the scattering, not a defect), reads on average at least ten times its ring, and its -// mean excess is above the Poisson bound: a weaker one cannot make an outlier. A pixel holding the +// A persistent pixel is masked if it reads on average at least ten times its ring and its mean excess +// is above the Poisson bound: a weaker one cannot make an outlier. Alone or as part of a patch - a +// patch lit through the whole sweep is stationary in the lab, so it is not a reflection of the +// rotating crystal, and a reflection crossing it reads the patch. A pixel holding the // detector's error value on most frames is masked with it: every frame already treats it as invalid, // but the static mask did not know. The detector's outermost row and column get no rule of their // own - a hot pixel there is caught like any other. diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 596ee0b0e..73b4ee4f4 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -7886,12 +7886,30 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const gemmi::Mat33 from_primitive = to_primitive.inverse(); // An operator counts as real on the same bound Stage A calls one present at - // SearchSpaceGroupOptions::min_operator_cc, applied by OperatorCorrelation over - // the same reflection population, so nothing new is calibrated here. + // the same reflection population - AND on the R contrast every promotion's added + // operators are held to (SearchSpaceGroupOptions::min_operator_r_contrast, + // against the same two ends the search above measured), so nothing new is + // calibrated here. The correlation alone cannot tell a pseudo-symmetric + // structure's near-operators from real ones: measured on a pseudo-cubic + // tetragonal crystal, two of the false cubic two-folds correlate at 0.31-0.33 and + // pass it, LePage closes them with the genuine ones into the whole cubic metric, + // and the tetragonal class between the two is never asked about. Their R reads + // 0.39 against 0.013 for the genuine operators and 0.49 for unrelated + // reflections - a contrast of 0.2 against the genuine operators' 1.0. + const double r_random = sg_search.random_pairing_r; + const double r_best = sg_search.global_best_operator_r; + // Read only where the search itself reads it: a best operator close to the + // random end leaves no scale to measure on (max_reference_over_random). + const bool contrast_known = std::isfinite(r_random) && std::isfinite(r_best) + && r_best <= sg_opts.max_reference_over_random * r_random; const auto keep_operator = [&](const gemmi::Mat33 &m_primitive) { const gemmi::Mat33 m_indexed = from_primitive.multiply(m_primitive).multiply(to_primitive); const auto score = OperatorCorrelation(sm.merged, m_indexed, sg_opts); - return score.has_value() && score->present; + if (!score.has_value() || !score->present) + return false; + return !contrast_known + || (r_random - score->r_stat) / (r_random - r_best) >= sg_opts.min_operator_r_contrast; }; const auto sub = LePageLattice(primitive, LATTICE_MAX_OBLIQUITY_DEG, keep_operator); if (sub && static_cast(sub->n_operators) > sg_search.point_group_order diff --git a/tests/HotPixelFinderTest.cpp b/tests/HotPixelFinderTest.cpp index a90adfb49..9d72f54b3 100644 --- a/tests/HotPixelFinderTest.cpp +++ b/tests/HotPixelFinderTest.cpp @@ -25,7 +25,8 @@ namespace { // high on every frame, one only 20 high - persistent, but too weak to make an outlier - one that holds // the error value on every frame, and one crossed by a genuine reflection. The last sits close to the rotation axis, where one reflection stays on a pixel for a // long stretch of rotation - here nine consecutive sampled frames, more than the chance bound allows -// but fewer than the one-reflection bound at that zeta. Only the first and the third are masked. +// but fewer than the one-reflection bound at that zeta. Only the first and the third are masked - and +// a 3x3 patch reading high on every frame, which is stationary in the lab and so no reflection either. TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") { DiffractionExperiment x(DetDECTRIS(W, H, "Test detector", "")); x.IncidentEnergy_keV(WVL_1A_IN_KEV).DetectorDistance_mm(10.0f); @@ -36,6 +37,7 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") { constexpr int HOT_X = 200, HOT_Y = 200, WARM_X = 190, WARM_Y = 60, ERR_X = 60, ERR_Y = 190; constexpr int BRAGG_X = 40, BRAGG_Y = 115; + constexpr int PATCH_X = 150, PATCH_Y = 225; std::mt19937 rng(1); // The powder ring is a smooth radial profile, as a real one is: 45 counts over the background at // 82 px, 4 px sigma. @@ -50,6 +52,9 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") { frame[I(HOT_X, HOT_Y)] += 60; frame[I(WARM_X, WARM_Y)] += 20; frame[I(ERR_X, ERR_Y)] = INT32_MIN; + for (int dy = -1; dy <= 1; dy++) + for (int dx = -1; dx <= 1; dx++) + frame[I(PATCH_X + dx, PATCH_Y + dy)] += 60; if (f >= 20 && f < 29) frame[I(BRAGG_X, BRAGG_Y)] += 500; finder.AddImage(frame.data(), scratch); @@ -61,7 +66,10 @@ TEST_CASE("HotPixelFinder_PersistentPixelNotBragg", "[HotPixelFinder]") { CHECK(result.mask[I(ERR_X, ERR_Y)] == 2); CHECK(result.mask[I(WARM_X, WARM_Y)] == 0); CHECK(result.mask[I(BRAGG_X, BRAGG_Y)] == 0); - CHECK(result.hot == 1); + for (int dy = -1; dy <= 1; dy++) + for (int dx = -1; dx <= 1; dx++) + CHECK(result.mask[I(PATCH_X + dx, PATCH_Y + dy)] == 1); + CHECK(result.hot == 10); CHECK(result.error == 1); } @@ -122,6 +130,8 @@ TEST_CASE("HotPixelFinder_DeviceMatchesHost", "[HotPixelFinder]") { REQUIRE(cudaMemcpy(device_image, image.data(), image.size() * sizeof(int32_t), cudaMemcpyHostToDevice) == cudaSuccess); device.AddDeviceImage(device_image, frame); + // The frame is queued, not processed: the next copy must not overwrite it before its kernels ran. + REQUIRE(cudaStreamSynchronize(*frame.stream) == cudaSuccess); } const auto expected = host.GetMask(OSC_DEG, SPACING_DEG, 4);