GivenCellPopulation: the population nearest the given cell among those holding a real share of the frames

The peak nearest the given cell could be a single stray frame when the given cell is a few percent off
the crystal (a cryo cell given for room-temperature lysozyme stills): the run then merged one frame.
Every population is now found from its own peak and valleys, and of those holding at least a tenth of
the frames the one centred nearest the given cell is kept.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

image_analysis/scale_merge/CellPopulation.cpp

@@ -14,6 +14,8 @@ std::pair<double, double> GivenCellPopulation(const std::vector<double> &relativ
     // A valley separates two populations only where it dips below this fraction of the kept peak; a
     // shallower dip is noise in one population's tail.
     constexpr double VALLEY_FRACTION = 0.5;
+    // A population holding less than this share of the frames is a stray group, not a crystal form.
+    constexpr double MIN_SHARE = 0.1;
 
     const size_t n = relative_volume.size();
     if (n < MIN_FRAMES)
@@ -47,18 +49,10 @@ std::pair<double, double> GivenCellPopulation(const std::vector<double> &relativ
         smooth[i] = s / (b - a + 1);
     }
 
-    // The peak nearest the given cell: from its bin, uphill.
-    int peak = std::clamp(static_cast<int>(std::floor(-lo / bin)), 0, nbins - 1);
-    while (true) {
-        if (peak + 1 < nbins && smooth[peak + 1] > smooth[peak]) peak++;
-        else if (peak > 0 && smooth[peak - 1] > smooth[peak]) peak--;
-        else break;
-    }
-
-    // Out from the peak on each side: the valley is the lowest point, below VALLEY_FRACTION of the peak,
+    // Out from a peak on each side: the valley is the lowest point, below VALLEY_FRACTION of the peak,
     // before the histogram rises again to twice its height by more than counting noise - the next
     // population, however small. None: the population runs out to the end of the data.
-    const auto valley = [&](int step) {
+    const auto valley = [&](int peak, int step) {
         int lowest = peak;
         for (int i = peak + step; i >= 0 && i < nbins; i += step) {
             if (smooth[i] < smooth[lowest])
@@ -71,7 +65,32 @@ std::pair<double, double> GivenCellPopulation(const std::vector<double> &relativ
         }
         return -1;
     };
-    const int left = valley(-1), right = valley(+1);
-    return {left < 0 ? -everything : lo + (left + 0.5) * bin,
-            right < 0 ? everything : lo + (right + 0.5) * bin};
+
+    // Every population: a local maximum of the histogram and the valleys about it. Of those that hold
+    // a real share of the frames, the one centred nearest the given cell is kept - not merely the peak
+    // nearest it, which can be a few stray frames when the given cell is a percent off the crystal.
+    double best_low = -everything, best_high = everything, best_distance = everything;
+    for (int i = 0; i < nbins; i++) {
+        if (smooth[i] <= 0 || (i > 0 && smooth[i - 1] > smooth[i]) || (i + 1 < nbins && smooth[i + 1] >= smooth[i]))
+            continue;
+        const int left = valley(i, -1), right = valley(i, +1);
+        const double low = left < 0 ? -everything : lo + (left + 0.5) * bin;
+        const double high = right < 0 ? everything : lo + (right + 0.5) * bin;
+        size_t count = 0;
+        double sum = 0;
+        for (double x : v)
+            if (x >= low && x <= high) {
+                count++;
+                sum += x;
+            }
+        if (count < MIN_SHARE * static_cast<double>(n))
+            continue;
+        const double distance = std::fabs(sum / static_cast<double>(count));
+        if (distance < best_distance) {
+            best_distance = distance;
+            best_low = low;
+            best_high = high;
+        }
+    }
+    return {best_low, best_high};
 }

tests/CellPopulationTest.cpp

@@ -43,6 +43,21 @@ TEST_CASE("CellPopulation_SecondFormLeftOut", "[CellPopulation][portable]") {
     CHECK(kept(r2, 0, 450) < 5);
 }
 
+TEST_CASE("CellPopulation_GivenCellOffTheCrystal", "[CellPopulation][portable]") {
+    // A given cell 3.5% smaller than the crystal's, and a single stray frame right at it: the crystal's
+    // population is kept, not the stray frame.
+    std::mt19937 rng(3);
+    std::normal_distribution<double> d(0.035, 0.003);
+    std::vector<double> r(2000);
+    for (auto &x : r) x = d(rng);
+    r.push_back(0.0);
+    const auto [low, high] = GivenCellPopulation(r);
+    size_t kept = 0;
+    for (double x : r)
+        kept += (x >= low && x <= high);
+    CHECK(kept >= 2000);
+}
+
 TEST_CASE("CellPopulation_TooFewFrames", "[CellPopulation][portable]") {
     const auto [low, high] = GivenCellPopulation({0.0, 0.1, 0.2});
     CHECK(std::isinf(low));