diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp index 863843797..9a344cd9a 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp +++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp @@ -3,6 +3,10 @@ #include "ImagePreprocessorCPU.h" +#include +#include +#include + ImagePreprocessorCPU::ImagePreprocessorCPU(const DiffractionExperiment &experiment, const PixelMask &mask) : ImagePreprocessor(experiment), mask_bits(mask.GetPackedMask()) {} @@ -45,25 +49,50 @@ void ImagePreprocessorCPU::AnalyzeBlock(ImagePreprocessorBuffer &processed_image if (sat_pixel_val > saturation_limit) sat_pixel_val = static_cast(saturation_limit); - for (size_t i = 0; i < n; i++) { - if ((mask_bits[(first + i) / 32] >> ((first + i) % 32)) & 1U) { - out[i] = INT32_MIN; - ++ret.masked_pixel_count; - } else if (image[i] == err_pixel_val) { - // Error/invalid marker = the pixel type's extreme value (0xFFFFFFFF for EIGER uint32). - // Tested before saturation, since for unsigned types the marker also exceeds sat_pixel_val - // (which is clipped above to the HDF5 saturation_value). - out[i] = INT32_MIN; - ++ret.error_pixel_count; - } else if (image[i] >= sat_pixel_val) { - out[i] = INT32_MAX; - ++ret.saturated_pixel_count; + // Every pixel is classified without a branch, so that the loop over the pixels vectorises. The + // extremes are taken over keys that order as the values do (a signed type's sign bit flipped), and a + // pixel that does not count gives the key no maximum or minimum can be moved by; the counts are + // integers. So the statistics are those of the pixel-by-pixel tests, whatever the order. + using U = std::make_unsigned_t; + constexpr U flip = std::is_signed_v ? static_cast(U{1} << (8 * sizeof(T) - 1)) : U{0}; + uint64_t n_masked = 0, n_error = 0, n_saturated = 0; + U key_max = 0, key_min = std::numeric_limits::max(); + const auto pixel = [&](size_t i, bool masked) { + const T v = image[i]; + // Error/invalid marker = the pixel type's extreme value (0xFFFFFFFF for EIGER uint32). Tested + // before saturation, since for unsigned types the marker also exceeds sat_pixel_val (which is + // clipped above to the HDF5 saturation_value). + const bool error = !masked & (v == err_pixel_val); + const bool saturated = !masked & !error & (v >= sat_pixel_val); + const U counted = static_cast(-static_cast(!masked & !error & !saturated)); // all bits, or none + out[i] = masked | error ? INT32_MIN : saturated ? INT32_MAX : static_cast(v); + n_masked += masked; + n_error += error; + n_saturated += saturated; + const U key = static_cast(v) ^ flip; + key_max = std::max(key_max, key & counted); + key_min = std::min(key_min, key | static_cast(~counted)); + }; + // Runs of whole mask words with no pixel masked - nearly all of a detector - take a loop of their + // own, which the compiler vectorises; the rest is tested pixel by pixel. + for (size_t i = 0; i < n;) { + size_t end = i; + if ((first + i) % 32 == 0) + while (end + 32 <= n && mask_bits[(first + end) / 32] == 0) + end += 32; + if (end > i) { + for (; i < end; i++) + pixel(i, false); } else { - out[i] = static_cast(image[i]); - if (image[i] > ret.max_value) - ret.max_value = image[i]; - if (image[i] < ret.min_value) - ret.min_value = image[i]; + pixel(i, (mask_bits[(first + i) / 32] >> ((first + i) % 32)) & 1U); + i++; } } + ret.masked_pixel_count += n_masked; + ret.error_pixel_count += n_error; + ret.saturated_pixel_count += n_saturated; + if (n_masked + n_error + n_saturated < n) { // some pixel counted + ret.max_value = std::max(ret.max_value, static_cast(key_max ^ flip)); + ret.min_value = std::min(ret.min_value, static_cast(key_min ^ flip)); + } }