diff --git a/rugnux/HotPixels.cpp b/rugnux/HotPixels.cpp index 65804f8ac..83e6450fd 100644 --- a/rugnux/HotPixels.cpp +++ b/rugnux/HotPixels.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include "../common/JFJochMath.h" #include "../common/ParallelFor.h" @@ -19,6 +20,41 @@ int32_t median_of(int32_t *v, size_t n) { return v[n / 2]; } +// A ring's background is a few counts, so its order statistics are read off a histogram of the values +// 0..RING_HIST_VALUES-1 - exactly, as long as the one asked for lands inside that range. +constexpr int32_t RING_HIST_VALUES = 1024; + +// The rank-th smallest value, if it is inside the histogram; `below` counts the values under 0. +std::optional hist_order_statistic(const std::vector &hist, size_t below, size_t rank) { + if (rank < below) + return std::nullopt; + size_t seen = below; + for (int32_t v = 0; v < RING_HIST_VALUES; v++) { + seen += hist[v]; + if (seen > rank) + return v; + } + return std::nullopt; +} + +// The rank-th smallest |value - m|, if every value that close to m is inside the histogram (none +// below 0 anywhere, and m + d still in range). +std::optional hist_deviation_order_statistic(const std::vector &hist, size_t below, + int32_t m, size_t rank) { + if (below > 0) + return std::nullopt; + size_t seen = hist[m]; + for (int32_t d = 0; seen <= rank; ) { + d++; + if (m + d >= RING_HIST_VALUES) + return std::nullopt; + seen += hist[m + d] + (m - d >= 0 ? hist[m - d] : 0); + if (seen > rank) + return d; + } + return 0; +} + // The smallest k with npix * P(Binomial(n, q) >= k) below the family-wise rate, n + 1 if none. int binomial_bound(int n, double q, double npix, double rate) { double tail = 0.0; @@ -89,21 +125,45 @@ void HotPixelFinder::AddImage(const int32_t *image, std::vector &scratc if (count[k] >= MIN_SECTOR_PIXELS) sector_level[k] = median_of(scratch.data() + key_begin[k], count[k]); - // Each ring's sectors packed together - forward into the space the ring owns, so nothing is - // overwritten before it is moved - then its median and its robust spread about that median. + // Each ring's median and its robust spread about that median: off the histogram of its values + // where both land inside it, otherwise with the ring's sectors packed together - forward into the + // space the ring owns, so nothing is overwritten before it is moved - and selected. std::vector ring_level(nrings, 0); std::vector ring_spread(nrings, 0.0f); std::vector ring_ok(nrings, 0); + std::vector hist(RING_HIST_VALUES); for (int r = 0; r < nrings; r++) { - int32_t *ring = scratch.data() + key_begin[r * SECTORS]; size_t n = 0; - for (int s = 0; s < SECTORS; s++) { - const size_t k = r * SECTORS + s; - std::memmove(ring + n, scratch.data() + key_begin[k], count[k] * sizeof(int32_t)); - n += count[k]; - } + for (int s = 0; s < SECTORS; s++) + n += count[r * SECTORS + s]; if (n < MIN_RING_PIXELS) continue; ring_ok[r] = 1; + + std::fill(hist.begin(), hist.end(), 0); + size_t below = 0; + for (int s = 0; s < SECTORS; s++) { + const size_t k = r * SECTORS + s; + for (size_t j = 0; j < count[k]; j++) { + const int32_t v = scratch[key_begin[k] + j]; + if (v < 0) below++; + else if (v < RING_HIST_VALUES) hist[v]++; + } + } + const auto median = hist_order_statistic(hist, below, n / 2); + const auto mad = median ? hist_deviation_order_statistic(hist, below, *median, n / 2) : std::nullopt; + if (median && mad) { + ring_level[r] = *median; + ring_spread[r] = 1.4826f * static_cast(*mad); + continue; + } + + int32_t *ring = scratch.data() + key_begin[r * SECTORS]; + size_t packed = 0; + for (int s = 0; s < SECTORS; s++) { + const size_t k = r * SECTORS + s; + std::memmove(ring + packed, scratch.data() + key_begin[k], count[k] * sizeof(int32_t)); + packed += count[k]; + } ring_level[r] = median_of(ring, n); for (size_t j = 0; j < n; j++) ring[j] = std::abs(ring[j] - ring_level[r]);