diff --git a/common/PixelMask.cpp b/common/PixelMask.cpp index a5338f4ee..8e455195e 100644 --- a/common/PixelMask.cpp +++ b/common/PixelMask.cpp @@ -59,6 +59,9 @@ void PixelMask::UpdateBinaryMask() { for (size_t i = 0; i < mask.size(); i++) binary_mask[i] = (mask[i] != 0); binary_mask_checksum = TableChecksum(binary_mask.data(), binary_mask.size()); + packed_mask.assign((mask.size() + 31) / 32, 0); + for (size_t i = 0; i < mask.size(); i++) + packed_mask[i / 32] |= static_cast(binary_mask[i]) << (i % 32); } void PixelMask::CalcEdgePixels_i(const DiffractionExperiment &experiment) { @@ -127,6 +130,10 @@ const std::vector &PixelMask::GetBinaryMask() const { return binary_mask; } +const std::vector &PixelMask::GetPackedMask() const { + return packed_mask; +} + uint64_t PixelMask::GetBinaryMaskChecksum() const { return binary_mask_checksum; } diff --git a/common/PixelMask.h b/common/PixelMask.h index 5408b0789..ddcf50be2 100644 --- a/common/PixelMask.h +++ b/common/PixelMask.h @@ -26,6 +26,9 @@ class PixelMask { // derived here once instead of in every one of those engines. std::vector binary_mask; uint64_t binary_mask_checksum = 0; + // The same flag packed 32 pixels to a word (pixel i is bit i % 32 of word i / 32) - the form the CPU + // image preprocessor reads, derived here for the same reason. + std::vector packed_mask; uint32_t LoadMask(const std::vector& mask, uint8_t bit); // Everything that follows from the mask, recomputed wherever the mask changes. @@ -68,6 +71,7 @@ public: [[nodiscard]] const std::vector &GetMask(const DiffractionExperiment& experiment) const; [[nodiscard]] const std::vector &GetMask() const; [[nodiscard]] const std::vector &GetBinaryMask() const; + [[nodiscard]] const std::vector &GetPackedMask() const; [[nodiscard]] uint64_t GetBinaryMaskChecksum() const; [[nodiscard]] std::vector GetUserMask(const DiffractionExperiment& experiment) const; diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp index 8803a73db..bb3a98931 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp +++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp @@ -5,10 +5,7 @@ ImagePreprocessorCPU::ImagePreprocessorCPU(const DiffractionExperiment &experiment, const PixelMask &mask) : ImagePreprocessor(experiment), - mask_1bit(npixels, false) { - for (int i = 0; i < npixels; i++) - mask_1bit[i] = (mask.GetMask().at(i) != 0); -} + mask_bits(mask.GetPackedMask()) {} ImageStatistics ImagePreprocessorCPU::Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *image_ptr, CompressedImageMode image_mode) { switch (image_mode) { @@ -43,7 +40,7 @@ ImageStatistics ImagePreprocessorCPU::Analyze(ImagePreprocessorBuffer &processed sat_pixel_val = static_cast(saturation_limit); for (int i = 0; i < npixels; i++) { - if (mask_1bit[i] != 0) { + if ((mask_bits[i / 32] >> (i % 32)) & 1U) { processed_image[i] = INT32_MIN; ++ret.masked_pixel_count; } else if (image[i] == err_pixel_val) { diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.h b/image_analysis/image_preprocessing/ImagePreprocessorCPU.h index aa326d77b..2b07bd54a 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorCPU.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.h @@ -7,7 +7,9 @@ class ImagePreprocessorCPU : public ImagePreprocessor { - std::vector mask_1bit; + // The pixel mask one bit per pixel, 32 to a word - PixelMask derives it once, so each engine only + // copies it. + std::vector mask_bits; template ImageStatistics Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *input, T err_value, T sat_value); public: diff --git a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp index 93b7b9597..ef66199a3 100644 --- a/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/AdaptiveSpotFinderCPU.cpp @@ -55,26 +55,58 @@ void AdaptiveSpotFinderCPU::AccumulateRings(const ImagePreprocessorBuffer &image std::fill(azint_sum2.begin(), azint_sum2.end(), 0.0f); std::fill(azint_count.begin(), azint_count.end(), 0); } + // Consecutive pixels mostly share a ring, so the ring's sums are held in locals while they do + // and written back when the ring changes: the same additions in the same order, without a + // store and a reload of the same address on every pixel. + uint32_t *hist = ring_hist.data(); + size_t cur = nbins; // the ring held in the locals below; nbins = none + int64_t sum = 0, cnt = 0; + uint64_t sum2 = 0; + float az_sum = 0.0f, az_sum2 = 0.0f; + uint32_t az_cnt = 0; + const auto store = [&] { + if (cur == nbins) + return; + ring_sum[cur] = sum; + ring_sum2[cur] = sum2; + ring_cnt[cur] = cnt; + if (fuse_azint) { + azint_sum[cur] = az_sum; + azint_sum2[cur] = az_sum2; + azint_count[cur] = az_cnt; + } + }; for (size_t pxl = 0; pxl < npix; ++pxl) { const int32_t v = image[pxl]; if (v == INT32_MIN || v == INT32_MAX) continue; // bad / saturated const uint16_t b = pixel_to_bin[pxl]; if (b >= nbins) continue; // masked / out of range (UINT16_MAX) + if (b != cur) { + store(); + cur = b; + sum = ring_sum[b]; + sum2 = ring_sum2[b]; + cnt = ring_cnt[b]; + az_sum = azint_sum[b]; + az_sum2 = azint_sum2[b]; + az_cnt = azint_count[b]; + } if (fuse_azint) { const float val = static_cast(v) * corrections[pxl]; const float val_sq = val * val; - azint_sum[b] += val; - azint_sum2[b] += val_sq; - ++azint_count[b]; + az_sum += val; + az_sum2 += val_sq; + ++az_cnt; } - ring_sum[b] += v; - ring_sum2[b] += static_cast(static_cast(v) * v); - ring_cnt[b] += 1; + sum += v; + sum2 += static_cast(static_cast(v) * v); + cnt += 1; if (v >= 0 && v < HIST_VALUES) - ring_hist[b * HIST_VALUES + v] += 1; + hist[b * HIST_VALUES + v] += 1; else ring_overflow.emplace_back(b, v); } + store(); } else { // A sigma-clip pass over the plain pass's values: each distinct value of a ring meets the same // test the pixels holding it would, and its pixels are added as a count. diff --git a/image_analysis/spot_finding/ImageSpotFinderCPU.cpp b/image_analysis/spot_finding/ImageSpotFinderCPU.cpp index e2a41bc15..b69993fa6 100644 --- a/image_analysis/spot_finding/ImageSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/ImageSpotFinderCPU.cpp @@ -146,42 +146,61 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, std::vector sum2_vert(width, 0); std::vector valid_vert(width, 0); - for (int line = 0; line < NBX; line++) { - for (int col = 0; col < width; col++) { - auto pxl = line * width + col; - - int64_t tmp = value_at(pxl); - if (tmp != INT32_MAX && tmp != INT32_MIN) { - sum_vert[col] += tmp; - sum2_vert[col] += tmp * tmp; - valid_vert[col] += 1; - } + // Add the row entering the window (line_in) and take out the one leaving it (line_out); -1 is + // no row. A bad or saturated pixel adds 0. Written without branches, so the loop vectorises. + // A pixel strong in the previous pass is not counted either; it is added here with the rest + // of its row and taken out again below - there are few, and the sums are integers, so the + // totals are the same as never adding it. + const int32_t *img = image.data(); + const int32_t w = width; + int64_t *sv = sum_vert.data(); + int64_t *sv2 = sum2_vert.data(); + uint16_t *vv = valid_vert.data(); + auto slide_vert = [&](int line_in, int line_out) { + for (int32_t col = 0; col < w; col++) { + const int64_t in = line_in >= 0 ? img[line_in * w + col] : INT32_MIN; + const int64_t out = line_out >= 0 ? img[line_out * w + col] : INT32_MIN; + const bool in_ok = in != INT32_MAX && in != INT32_MIN; + const bool out_ok = out != INT32_MAX && out != INT32_MIN; + const int64_t in_v = in_ok ? in : 0; + const int64_t out_v = out_ok ? out : 0; + sv[col] += in_v - out_v; + sv2[col] += in_v * in_v - out_v * out_v; + vv[col] += static_cast(in_ok - out_ok); } - } + if (!prev_strong) + return; + // sign = +1 takes out what line_in added, -1 puts back what line_out took out. + const auto uncount_strong = [&](int line, int64_t sign) { + const int32_t first = line * w, last = first + w - 1; + for (int32_t word = first / 32; word <= last / 32; word++) { + uint32_t bits = prev_strong[word]; + if (word == first / 32) + bits &= UINT32_MAX << (first % 32); + if (word == last / 32 && last % 32 != 31) + bits &= (1U << (last % 32 + 1)) - 1; + for (; bits; bits &= bits - 1) { + const int32_t pxl = word * 32 + std::countr_zero(bits); + const int64_t v = img[pxl]; + if (v == INT32_MAX || v == INT32_MIN) + continue; + sv[pxl - first] -= sign * v; + sv2[pxl - first] -= sign * v * v; + vv[pxl - first] -= static_cast(sign); + } + } + }; + if (line_in >= 0) + uncount_strong(line_in, 1); + if (line_out >= 0) + uncount_strong(line_out, -1); + }; + + for (int line = 0; line < NBX; line++) + slide_vert(line, -1); for (int line = 0; line < height; line++) { - for (int col = 0; col < width; col++) { - if (line < height - NBX) { - auto pxl = (line + NBX) * width + col; - - int64_t tmp = value_at(pxl); - if (tmp != INT32_MAX && tmp != INT32_MIN) { - sum_vert[col] += tmp; - sum2_vert[col] += tmp * tmp; - valid_vert[col] += 1; - } - } - - if (line >= NBX + 1) { - auto pxl = (line - (NBX + 1)) * width + col; - int64_t tmp = value_at(pxl); - if (tmp != INT32_MAX && tmp != INT32_MIN) { - sum_vert[col] -= tmp; - sum2_vert[col] -= tmp * tmp; - valid_vert[col] -= 1; - } - } - } + slide_vert(line < height - NBX ? line + NBX : -1, line >= NBX + 1 ? line - (NBX + 1) : -1); if (candidates) { // Only the runs of 32-column blocks whose result is read (first_pass_needed; every