diff --git a/image_analysis/scale_merge/CellPopulation.cpp b/image_analysis/scale_merge/CellPopulation.cpp index ff8f7b145..601708b4e 100644 --- a/image_analysis/scale_merge/CellPopulation.cpp +++ b/image_analysis/scale_merge/CellPopulation.cpp @@ -14,6 +14,8 @@ std::pair GivenCellPopulation(const std::vector &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 GivenCellPopulation(const std::vector &relativ smooth[i] = s / (b - a + 1); } - // The peak nearest the given cell: from its bin, uphill. - int peak = std::clamp(static_cast(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 GivenCellPopulation(const std::vector &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(n)) + continue; + const double distance = std::fabs(sum / static_cast(count)); + if (distance < best_distance) { + best_distance = distance; + best_low = low; + best_high = high; + } + } + return {best_low, best_high}; } diff --git a/tests/CellPopulationTest.cpp b/tests/CellPopulationTest.cpp index e7fc0f248..31335a5e1 100644 --- a/tests/CellPopulationTest.cpp +++ b/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 d(0.035, 0.003); + std::vector 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));