diff --git a/common/AssembleToPlane.cpp b/common/AssembleToPlane.cpp new file mode 100644 index 000000000..49afa4d89 --- /dev/null +++ b/common/AssembleToPlane.cpp @@ -0,0 +1,16 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "AssembleToPlane.h" + +#include + +PlaneRaster PlaneRasterOf(const DetectorModules &modules, int64_t width, int64_t height) { + if (modules.empty()) + return {0, 0, width, height}; + // Raster pixel centres are whole plane coordinates; the first and last that lie on the footprints. + const PlaneBox box = modules.PlaneBounds(width, height); + const auto x0 = std::min(0, static_cast(std::ceil(box.x0))); + const auto y0 = std::min(0, static_cast(std::ceil(box.y0))); + return {x0, y0, static_cast(std::ceil(box.x1)) - x0, static_cast(std::ceil(box.y1)) - y0}; +} diff --git a/common/AssembleToPlane.h b/common/AssembleToPlane.h new file mode 100644 index 000000000..42d6a416f --- /dev/null +++ b/common/AssembleToPlane.h @@ -0,0 +1,68 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +#include "DetectorModules.h" + +// The detector plane as a picture: raster pixel (i, j) is the plane point (x0 + i, y0 + j). For the +// algorithms that need the detector as it really is - a beam-stop shadow and the gaps it crosses, the point +// symmetry about the beam, a picture of the detector - and not as the panels are stored. +struct PlaneRaster { + int64_t x0 = 0, y0 = 0, width = 0, height = 0; +}; + +// The raster holding every panel. It starts at the plane origin, so plane and raster coordinates agree, +// unless a panel lies left of or above it. Empty table: the image itself. +PlaneRaster PlaneRasterOf(const DetectorModules &modules, int64_t width, int64_t height); + +// The image (width x height) as the plane raster: each raster pixel takes the image pixel nearest its centre +// on the panel it falls on, and `gap` where it falls on none. Empty table, with the image's own raster: a copy. +template +std::vector AssembleToPlane(const DetectorModules &modules, const std::vector &image, + int64_t width, int64_t height, const PlaneRaster &raster, T gap) { + std::vector out(static_cast(raster.width * raster.height), gap); + const auto put = [&](const DetectorModule &m, int64_t i0, int64_t i1, int64_t j0, int64_t j1) { + for (int64_t j = std::max(j0, 0); j < std::min(j1, raster.height); j++) + for (int64_t i = std::max(i0, 0); i < std::min(i1, raster.width); i++) { + const auto [x, y] = m.PlaneToImage(static_cast(raster.x0 + i), + static_cast(raster.y0 + j)); + if (m.ContainsImagePixel(x, y)) + out[j * raster.width + i] = image[static_cast(std::floor(y + 0.5)) * width + + static_cast(std::floor(x + 0.5))]; + } + }; + if (modules.empty()) { + put(DetectorModule::AsImaged(0, 0, width, height), -raster.x0, width - raster.x0, -raster.y0, + height - raster.y0); + return out; + } + // Panel by panel, over the raster pixels of its footprint's box. Panels do not overlap. + for (const auto &m: modules.Get()) { + const PlaneBox box = DetectorModules({m}).PlaneBounds(width, height); + put(m, static_cast(std::floor(box.x0)) - raster.x0, static_cast(std::ceil(box.x1)) - raster.x0, + static_cast(std::floor(box.y0)) - raster.y0, static_cast(std::ceil(box.y1)) - raster.y0); + } + return out; +} + +// The reverse, for a result computed on the raster: each image pixel on a panel takes the raster pixel +// nearest its centre on the plane, `outside` where that is off the raster; a pixel on no panel (a gap) +// is `outside` too. Empty table, with the image's own raster: a copy. +template +std::vector GatherFromPlane(const DetectorModules &modules, const std::vector &plane, + const PlaneRaster &raster, int64_t width, int64_t height, T outside) { + std::vector out(static_cast(width * height), outside); + modules.ForEachPixel(width, height, [&](int64_t, size_t index, double px, double py) { + const auto i = static_cast(std::floor(px + 0.5)) - raster.x0; + const auto j = static_cast(std::floor(py + 0.5)) - raster.y0; + if (i >= 0 && i < raster.width && j >= 0 && j < raster.height) + out[index] = plane[j * raster.width + i]; + }); + return out; +} diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index dcd7f7ddd..fd541c47b 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -132,6 +132,10 @@ ADD_LIBRARY(JFJochCommon STATIC DetectorGeometryModular.h DetectorModules.cpp DetectorModules.h + PanelBoxes.cpp + PanelBoxes.h + AssembleToPlane.cpp + AssembleToPlane.h ModuleShift.h DetectorGeometryFixed.cpp DetectorGeometryFixed.h diff --git a/common/DiffractionSpot.cpp b/common/DiffractionSpot.cpp index ba58bbe2a..14573df44 100644 --- a/common/DiffractionSpot.cpp +++ b/common/DiffractionSpot.cpp @@ -24,8 +24,9 @@ DiffractionSpot& DiffractionSpot::operator+=(const DiffractionSpot &other) { } DiffractionSpot::DiffractionSpot(int64_t x_sum, int64_t y_sum, int64_t in_pixel_count, - int64_t in_photons, int64_t in_max_photons) - : x(x_sum), y(y_sum), pixel_count(in_pixel_count), photons(in_photons), max_photons(in_max_photons) {} + int64_t in_photons, int64_t in_max_photons, int64_t in_panel) + : x(x_sum), y(y_sum), pixel_count(in_pixel_count), photons(in_photons), max_photons(in_max_photons), + panel(in_panel) {} int64_t DiffractionSpot::Count() const { return photons; @@ -44,6 +45,14 @@ Coord DiffractionSpot::RawCoord() const { static_cast(static_cast(y) / static_cast(photons)), 0}; } +int64_t DiffractionSpot::Panel() const { + return panel; +} + +void DiffractionSpot::Panel(int64_t in_panel) { + panel = in_panel; +} + int64_t DiffractionSpot::PixelCount() const { return pixel_count; } diff --git a/common/DiffractionSpot.h b/common/DiffractionSpot.h index 482f988c9..2c26533bf 100644 --- a/common/DiffractionSpot.h +++ b/common/DiffractionSpot.h @@ -21,17 +21,21 @@ class DiffractionSpot { int64_t pixel_count = 0; int64_t photons = 0; // total photon count int64_t max_photons = INT64_MIN; // maximum number of counts per pixel in the spot + int64_t panel = -1; // the panel (DetectorModules index) the spot was found on; -1 = no panel table public: DiffractionSpot() = default; DiffractionSpot(uint32_t col, uint32_t line, int64_t photons); // From already-summed quantities, i.e. the members as AddPixel leaves them, NOT a centroid. // Used by the GPU spot extractor, which computes those sums on the device. - DiffractionSpot(int64_t x_sum, int64_t y_sum, int64_t pixel_count, int64_t photons, int64_t max_photons); + DiffractionSpot(int64_t x_sum, int64_t y_sum, int64_t pixel_count, int64_t photons, int64_t max_photons, + int64_t panel = -1); DiffractionSpot& operator+=(const DiffractionSpot& spot); int64_t PixelCount() const; int64_t Count() const; int64_t MaxCount() const; Coord RawCoord() const; + int64_t Panel() const; + void Panel(int64_t panel); void ConvertToImageCoordinates(const DiffractionExperiment& experiment, uint16_t module_number); std::optional Export(const DiffractionGeometry &geometry, int64_t image_num = 0) const; void AddPixel(uint32_t col, uint32_t line, int64_t photons); diff --git a/common/PanelBoxes.cpp b/common/PanelBoxes.cpp new file mode 100644 index 000000000..ad7ca7d32 --- /dev/null +++ b/common/PanelBoxes.cpp @@ -0,0 +1,17 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "PanelBoxes.h" + +std::vector PanelBoxes(const DetectorModules &modules, int64_t width, int64_t height) { + if (modules.empty()) + return {PanelBox{0, 0, static_cast(width), static_cast(height), -1}}; + std::vector out; + for (size_t p = 0; p < modules.size(); p++) { + const DetectorModule &m = modules[p]; + out.push_back(PanelBox{static_cast(m.x0), static_cast(m.y0), + static_cast(m.x0 + m.width), static_cast(m.y0 + m.height), + static_cast(p)}); + } + return out; +} diff --git a/common/PanelBoxes.h b/common/PanelBoxes.h new file mode 100644 index 000000000..132769530 --- /dev/null +++ b/common/PanelBoxes.h @@ -0,0 +1,23 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "DetectorModules.h" + +// One panel's box in the image rugnux reads, x0 <= x < x1, y0 <= y < y1, and its id in DetectorModules +// (-1 for the one panel of an empty table). What the neighbourhood algorithms on the image - the spot +// finders' windows, the connected-component search, the squares read about a spot - stay inside, the way +// the FPGA spot finder stays inside a JUNGFRAU module: stacked panels lie next to each other in the image +// and far apart on the detector. +struct PanelBox { + int32_t x0, y0, x1, y1; + int32_t panel; +}; + +// One box per panel. An empty table gives one box, the whole image, so whatever is clamped to it is +// clamped to the image exactly as without panels. +std::vector PanelBoxes(const DetectorModules &modules, int64_t width, int64_t height); diff --git a/image_analysis/MXAnalysisWithoutFPGA.cpp b/image_analysis/MXAnalysisWithoutFPGA.cpp index 0e50fbc00..70c4e0dd1 100644 --- a/image_analysis/MXAnalysisWithoutFPGA.cpp +++ b/image_analysis/MXAnalysisWithoutFPGA.cpp @@ -106,6 +106,11 @@ MXAnalysisWithoutFPGA::MXAnalysisWithoutFPGA(const DiffractionExperiment &in_exp } adaptiveSpotFinder = std::move(cpu); } + // The panels the image is made of: the spot finders' windows and components, and the squares read + // about a spot, stay inside one (one, the whole image, without a module table). + adaptiveSpotFinder->SetPanels(experiment.GetModules()); + preprocessor_buffer->SetPanels(PanelBoxes(experiment.GetModules(), experiment.GetXPixelsNum(), + experiment.GetYPixelsNum())); } void MXAnalysisWithoutFPGA::Analyze(DataMessage &output, @@ -387,6 +392,7 @@ ImageSpotFinder &MXAnalysisWithoutFPGA::FixedThresholdFinder() { #endif spotFinder = std::make_unique(experiment.GetXPixelsNum(), experiment.GetYPixelsNum()); + spotFinder->SetPanels(experiment.GetModules()); // It missed every mask update that happened before it existed, so it takes the current one // now. Without this it would find spots outside the resolution limits. if (mask_resolution) diff --git a/image_analysis/beam_stop/ShadowFinder.cpp b/image_analysis/beam_stop/ShadowFinder.cpp index e07d57aff..3fa23b110 100644 --- a/image_analysis/beam_stop/ShadowFinder.cpp +++ b/image_analysis/beam_stop/ShadowFinder.cpp @@ -15,6 +15,7 @@ #include +#include "../../common/AssembleToPlane.h" #include "../../common/CUDAWrapper.h" #include "../../common/ParallelFor.h" #include "../../common/JFJochException.h" @@ -407,7 +408,8 @@ ShadowFinder::ShadowFinder(const DiffractionExperiment &experiment, const PixelM beam_y(experiment.GetBeamY_pxl()), geometry(experiment.GetDiffractionGeometry()), polarization(experiment.GetPolarizationFactor()), - pixel_mask(mask.GetMask(experiment)) { + pixel_mask(mask.GetMask(experiment)), + modules(experiment.GetModules()) { if (pixel_mask.size() != static_cast(width) * height) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "ShadowFinder: pixel mask does not match the detector"); @@ -650,8 +652,9 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { if (nthreads == 0) nthreads = AvailableCpus(); #ifdef JFJOCH_USE_CUDA - // Where every frame went to the device the mask is made there too, from the projection as it lies. - if (ShadowAccumulatorGPU *gpu = Gpu(); gpu && gpu->GetFrameCount() > 0 && host.frames == 0) { + // Where every frame went to the device the mask is made there too, from the projection as it lies - + // when it lies as the detector does (no panel table). + if (ShadowAccumulatorGPU *gpu = Gpu(); gpu && gpu->GetFrameCount() > 0 && host.frames == 0 && modules.empty()) { const float diag = std::hypot(static_cast(width), static_cast(height)); if (!std::isfinite(beam_x) || !std::isfinite(beam_y) || std::fabs(beam_x - width * 0.5f) > 4.0f * diag || std::fabs(beam_y - height * 0.5f) > 4.0f * diag) @@ -673,12 +676,32 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { } #endif const Projection &p = Reduced(); + if (modules.empty()) + return MaskOn(p, pixel_mask, width, height, beam_x, beam_y, nthreads); + + // Panels: the shadow and the gaps it crosses are on the detector plane, not in the image the panels are + // stored in, so the mask is made on the plane raster (the projection and the mask assembled to it, gaps + // counted as nothing and masked) and each stored pixel takes the answer of the raster pixel it lies on. + const PlaneRaster raster = PlaneRasterOf(modules, width, height); + Projection plane; + plane.max_value = AssembleToPlane(modules, p.max_value, width, height, raster, int64_t{0}); + plane.sum_value = AssembleToPlane(modules, p.sum_value, width, height, raster, int64_t{0}); + plane.valid_count = AssembleToPlane(modules, p.valid_count, width, height, raster, uint32_t{0}); + plane.frames = p.frames; + const auto plane_mask = AssembleToPlane(modules, pixel_mask, width, height, raster, + uint32_t{1u << PixelMask::ModuleGapPixelBit}); + const auto mask = MaskOn(plane, plane_mask, static_cast(raster.width), static_cast(raster.height), + beam_x - static_cast(raster.x0), beam_y - static_cast(raster.y0), nthreads); + return GatherFromPlane(modules, mask, raster, width, height, uint32_t{0}); +} + +std::vector ShadowFinder::MaskOn(const Projection &p, const std::vector &known_mask, int W, int H, + float centre_x, float centre_y, size_t nthreads) const { const auto &max_value = p.max_value; const auto &sum_value = p.sum_value; const auto &valid_count = p.valid_count; const uint32_t frames = p.frames; - const int W = width, H = height; const int n_pixels = W * H; std::vector mask(n_pixels, 0); @@ -690,8 +713,8 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { // makes max_radius, and every per-ring table sized by it, arbitrarily large. No centre, no // shadow. const float diag = std::hypot(static_cast(W), static_cast(H)); - if (!std::isfinite(beam_x) || !std::isfinite(beam_y) - || std::fabs(beam_x - W * 0.5f) > 4.0f * diag || std::fabs(beam_y - H * 0.5f) > 4.0f * diag) + if (!std::isfinite(centre_x) || !std::isfinite(centre_y) + || std::fabs(centre_x - W * 0.5f) > 4.0f * diag || std::fabs(centre_y - H * 0.5f) > 4.0f * diag) return mask; // mean projection, usable pixels and the ring each pixel is on (ring_position). The mean is @@ -702,7 +725,7 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { // in-plane rotation - about the centre the caller named. // Following Kahn, Fourme, Gadet, Janin, Dumas & Andre (1982) J. Appl. Cryst. 15, 330-337 DiffractionGeometry pol_geometry = geometry; - pol_geometry.BeamX_pxl(beam_x).BeamY_pxl(beam_y); + pol_geometry.BeamX_pxl(centre_x).BeamY_pxl(centre_y); const std::vector det_matrix = geometry.GetDetectorMatrix().arr(); const float distance_px = geometry.GetDetectorDistance_mm() / geometry.GetPixelSize_mm(); @@ -716,7 +739,7 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { for (int y = ylo; y < yhi; y++) for (int x = 0; x < W; x++) { const int i = y * W + x; - const RingPosition ring = ring_position(det_matrix.data(), distance_px, x - beam_x, y - beam_y); + const RingPosition ring = ring_position(det_matrix.data(), distance_px, x - centre_x, y - centre_y); pol[i] = ring.obliquity * (polarization.has_value() ? pol_geometry.CalcAzIntPolarizationCorr(static_cast(x), @@ -729,7 +752,7 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { // the pooled means, which the running box sums then turn into NaN for a whole row. mean[i] = 0.0f; valid[i] = 0; - if (valid_count[i] > 0 && pixel_mask[i] == 0 && pol[i] > 0.0f) { + if (valid_count[i] > 0 && known_mask[i] == 0 && pol[i] > 0.0f) { mean[i] = static_cast(static_cast(sum_value[i]) / valid_count[i] / pol[i]); valid[i] = 1; @@ -970,8 +993,8 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { if (!valid[i] || region[i]) continue; const RingPosition ring = ring_position(det_matrix.data(), distance_px, - static_cast(i % W) - beam_x, - static_cast(i / W) - beam_y); + static_cast(i % W) - centre_x, + static_cast(i / W) - centre_y); const double phi = std::atan2(ring.ly, ring.lx) + std::numbers::pi; const int sector = std::min(HARMONIC_SECTORS - 1, static_cast(phi / (2.0 * std::numbers::pi) * HARMONIC_SECTORS)); values[static_cast(radius[i] / HARMONIC_BAND_PX) * HARMONIC_SECTORS + sector].push_back(ratio[i]); @@ -1002,8 +1025,8 @@ std::vector ShadowFinder::GetMask(size_t nthreads) const { continue; // cos 2phi and sin 2phi from the offset to the beam. const RingPosition ring = ring_position(det_matrix.data(), distance_px, - static_cast(i % W) - beam_x, - static_cast(i / W) - beam_y); + static_cast(i % W) - centre_x, + static_cast(i / W) - centre_y); const float dx = ring.lx, dy = ring.ly; const float r2 = std::max(dx * dx + dy * dy, 1e-6f); const int band = radius[i] / HARMONIC_BAND_PX; diff --git a/image_analysis/beam_stop/ShadowFinder.h b/image_analysis/beam_stop/ShadowFinder.h index 53efb99b3..1c844f5cd 100644 --- a/image_analysis/beam_stop/ShadowFinder.h +++ b/image_analysis/beam_stop/ShadowFinder.h @@ -63,6 +63,8 @@ class ShadowFinder { const std::optional polarization; // unset leaves the background as it was measured std::vector pixel_mask; // pixels already masked carry no background to test + // The panels the image is made of. With any, the mask is made on the detector plane (GetMask). + const DetectorModules modules; // Per-pixel projection over the frames added so far (converted geometry). The sums and counts are // integers and the maximum is a maximum, so the result does not depend on the order the frames @@ -120,6 +122,11 @@ class ShadowFinder { // took frames. Called with `m` held. [[nodiscard]] const Projection &Reduced() const; + // The mask of a W x H projection, its known-masked pixels and the centre the rings are drawn about - + // the image's own, or the plane raster's where the image is made of panels. + [[nodiscard]] std::vector MaskOn(const Projection &p, const std::vector &known_mask, + int W, int H, float centre_x, float centre_y, size_t nthreads) const; + public: static constexpr uint32_t SHADOW = 1; static constexpr uint32_t TRANSMITTING = 2; diff --git a/image_analysis/geom_refinement/BeamCenterFFT.cpp b/image_analysis/geom_refinement/BeamCenterFFT.cpp index 53c4468ac..36e7685ef 100644 --- a/image_analysis/geom_refinement/BeamCenterFFT.cpp +++ b/image_analysis/geom_refinement/BeamCenterFFT.cpp @@ -12,6 +12,7 @@ #include "BeamCenterFFTCPU.h" #include "BeamCenterFFTEngine.h" +#include "../../common/AssembleToPlane.h" #include "../../common/CUDAWrapper.h" #include #ifdef JFJOCH_USE_CUDA @@ -241,6 +242,20 @@ BeamCenterFFTResult FindBeamCenterFFT(const DiffractionExperiment &experiment, const std::vector &mean, const BeamCenterFFTSettings &settings) { // The pre-scan projection is in converted geometry (ShadowFinder's frame). - return BeamCenterFFTScore(experiment.GetXPixelsNumConv(), experiment.GetYPixelsNumConv(), - mean, settings); + const DetectorModules &modules = experiment.GetModules(); + const int64_t width = experiment.GetXPixelsNumConv(), height = experiment.GetYPixelsNumConv(); + if (modules.empty()) + return BeamCenterFFTScore(width, height, mean, settings); + + // Panels: the point symmetry is about the beam on the detector plane, so the projection is scored as + // the plane raster (gaps hold nothing), and the centres come back in plane coordinates. + const PlaneRaster raster = PlaneRasterOf(modules, width, height); + BeamCenterFFTResult result = BeamCenterFFTScore(raster.width, raster.height, + AssembleToPlane(modules, mean, width, height, raster, NAN), + settings); + for (auto &c: result.point) { + c.beam_x_pxl += static_cast(raster.x0); + c.beam_y_pxl += static_cast(raster.y0); + } + return result; } diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp index 01dd73c43..688606726 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp +++ b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.cpp @@ -23,17 +23,25 @@ const int32_t *ImagePreprocessorBuffer::getGPUBuffer() const { return nullptr; } +PanelBox ImagePreprocessorBuffer::TileBox(int32_t x, int32_t y, int width) const { + for (const auto &p: panels) + if (x >= p.x0 && x < p.x1 && y >= p.y0 && y < p.y1) + return p; + return PanelBox{0, 0, width, static_cast(npixels / width), -1}; +} + const int32_t *ImagePreprocessorBuffer::GatherTiles(const std::vector &x, const std::vector &y, int width, int half) { const int side = 2 * half + 1; - const int height = static_cast(npixels / width); tiles.resize(x.size() * side * side); int32_t *out = tiles.data(); - for (size_t i = 0; i < x.size(); i++) + for (size_t i = 0; i < x.size(); i++) { + const PanelBox b = TileBox(x[i], y[i], width); for (int py = y[i] - half; py <= y[i] + half; py++) for (int px = x[i] - half; px <= x[i] + half; px++) - *out++ = (px < 0 || py < 0 || px >= width || py >= height) + *out++ = (px < b.x0 || py < b.y0 || px >= b.x1 || py >= b.y1) ? INT32_MIN : buffer[static_cast(py) * width + px]; + } return tiles.data(); } diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h index 40fb715be..0b8372661 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h +++ b/image_analysis/image_preprocessing/ImagePreprocessorBuffer.h @@ -8,8 +8,15 @@ #include #include +#include "../../common/PanelBoxes.h" + class ImagePreprocessorBuffer { protected: + // The panels the image is made of (PanelBoxes); empty means one, the whole image. + std::vector panels; + // The box the square about pixel (x, y) is clipped to: its panel's, or the image's. + PanelBox TileBox(int32_t x, int32_t y, int width) const; + // The image size, kept apart from `buffer` so size() still answers when the host copy was not // allocated at all (see the host_mirror constructor). const size_t npixels; @@ -46,9 +53,13 @@ public: virtual void Gather(const std::vector &npixel, std::vector &values) const; // The (2*half+1)^2 pixels about each (x[i], y[i]) of the image, `width` pixels wide: square after - // square, each row by row, with a pixel off the image read as INT32_MIN. Valid until the next call. + // square, each row by row, with a pixel off the image - or off the panel (x[i], y[i]) is on - read as + // INT32_MIN. Valid until the next call. virtual const int32_t *GatherTiles(const std::vector &x, const std::vector &y, int width, int half); + // The panels of the image, for GatherTiles. Default: none, the whole image is one. + void SetPanels(const std::vector &in_panels) { panels = in_panels; } + // The whole image on the host - on the GPU copied there first, where it was not already. virtual const int32_t *HostImage(); private: diff --git a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu index 6c5981ada..f2f11ecec 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.cu @@ -19,7 +19,8 @@ __global__ void gather_kernel(const int32_t *__restrict__ image, values[i] = image[npixel[i]]; } -__global__ void gather_tiles_kernel(const int32_t *__restrict__ image, int width, int height, +// centre holds six values per tile: its centre x, y and the box x0, y0, x1, y1 it is clipped to. +__global__ void gather_tiles_kernel(const int32_t *__restrict__ image, int width, const int32_t *__restrict__ centre, int half, int32_t *__restrict__ tiles, size_t count) { const int side = 2 * half + 1; @@ -27,9 +28,10 @@ __global__ void gather_tiles_kernel(const int32_t *__restrict__ image, int width for (size_t i = blockIdx.x * blockDim.x + threadIdx.x; i < count; i += static_cast(blockDim.x) * gridDim.x) { const size_t tile = i / area; const int p = static_cast(i - tile * area); - const int px = centre[2 * tile] - half + p % side; - const int py = centre[2 * tile + 1] - half + p / side; - tiles[i] = (px < 0 || py < 0 || px >= width || py >= height) + const int32_t *c = centre + 6 * tile; + const int px = c[0] - half + p % side; + const int py = c[1] - half + p / side; + tiles[i] = (px < c[2] || py < c[3] || px >= c[4] || py >= c[5]) ? INT32_MIN : image[static_cast(py) * width + px]; } } @@ -79,18 +81,24 @@ const int32_t *ImagePreprocessorBufferGPU::GatherTiles(const std::vector(tile_capacity); host_tiles.resize(tile_capacity); } - if (2 * x.size() > tile_centre.size()) { - tile_centre.resize(std::max(2 * x.size(), 2 * tile_centre.size())); + if (6 * x.size() > tile_centre.size()) { + tile_centre.resize(std::max(6 * x.size(), 2 * tile_centre.size())); gpu_tile_centre = CudaDevicePtr(tile_centre.size()); } for (size_t i = 0; i < x.size(); i++) { - tile_centre[2 * i] = x[i]; - tile_centre[2 * i + 1] = y[i]; + const PanelBox b = TileBox(x[i], y[i], width); + int32_t *c = tile_centre.data() + 6 * i; + c[0] = x[i]; + c[1] = y[i]; + c[2] = b.x0; + c[3] = b.y0; + c[4] = b.x1; + c[5] = b.y1; } - cuda_err(cudaMemcpyAsync(gpu_tile_centre.get(), tile_centre.data(), 2 * x.size() * sizeof(int32_t), + cuda_err(cudaMemcpyAsync(gpu_tile_centre.get(), tile_centre.data(), 6 * x.size() * sizeof(int32_t), cudaMemcpyHostToDevice, gather_stream)); const int blocks = static_cast(std::min((count + 255) / 256, 4096)); - gather_tiles_kernel<<>>(gpu_image.get(), width, static_cast(npixels / width), + gather_tiles_kernel<<>>(gpu_image.get(), width, gpu_tile_centre.get(), half, gpu_tiles.get(), count); cuda_err(cudaGetLastError()); cuda_err(cudaMemcpyAsync(host_tiles.data(), gpu_tiles.get(), count * sizeof(int32_t), diff --git a/image_analysis/spot_finding/ImageSpotFinder.cpp b/image_analysis/spot_finding/ImageSpotFinder.cpp index 7eab95ccf..82efb22ed 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.cpp +++ b/image_analysis/spot_finding/ImageSpotFinder.cpp @@ -11,7 +11,8 @@ ImageSpotFinder::ImageSpotFinder(int32_t width, int32_t height, bool host_bit_bu : width(width), height(height), output_buffer(host_bit_buffer ? width * height / 32 + 1 : 0), - res_mask_bits(OutputSize(), 0) { + res_mask_bits(OutputSize(), 0), + panels(PanelBoxes({}, width, height)) { // Exclude the padding bits of the last word up front, so neither the host scan nor the GPU // compaction needs a separate "is this bit still inside the image?" test. const size_t npixel = static_cast(width) * height; @@ -47,6 +48,20 @@ void ImageSpotFinder::SetResolutionMaskBits(const std::vector &packed_ const size_t npixel = static_cast(width) * height; if (npixel % 32 != 0) res_mask_bits.back() |= ~((1u << (npixel % 32)) - 1u); + for (size_t i = 0; i < off_panel_bits.size(); i++) + res_mask_bits[i] |= off_panel_bits[i]; +} + +void ImageSpotFinder::SetPanels(const DetectorModules &modules) { + panels = PanelBoxes(modules, width, height); + off_panel_bits.clear(); + if (!modules.empty()) { + off_panel_bits.assign(OutputSize(), UINT32_MAX); + modules.ForEachPixel(width, height, [&](int64_t, size_t i, double, double) { + off_panel_bits[i / 32] &= ~(1u << (i % 32)); + }); + } + SetResolutionMaskBits(res_mask_bits); } const std::vector &ImageSpotFinder::GetRingBackground() const { @@ -79,10 +94,22 @@ void ImageSpotFinder::ExtractComponentsHost(const ImagePreprocessorBuffer &image image.Gather(strong_pixel, strong_pixel_value); - StrongPixelSet pixel_set; - for (size_t i = 0; i < strong_pixel.size(); i++) - pixel_set.AddStrongPixel(strong_pixel[i] % width, strong_pixel[i] / width, strong_pixel_value[i]); - pixel_set.FindComponentsImage(settings, components); + // Panel by panel, as the FPGA path labels module by module: a component never joins pixels of two + // panels, however close the image stores them. The strong pixels are in raster order, and so are + // those of one box. Every strong pixel is on a panel (off_panel_bits). + for (const auto &p: panels) { + StrongPixelSet pixel_set; + for (size_t i = 0; i < strong_pixel.size(); i++) { + const auto col = static_cast(strong_pixel[i] % width); + const auto line = static_cast(strong_pixel[i] / width); + if (col >= p.x0 && col < p.x1 && line >= p.y0 && line < p.y1) + pixel_set.AddStrongPixel(col, line, strong_pixel_value[i]); + } + const size_t first = components.size(); + pixel_set.FindComponentsImage(settings, components); + for (size_t i = first; i < components.size(); i++) + components[i].Panel(p.panel); + } } const std::vector &ImageSpotFinder::ExtractComponents(const ImagePreprocessorBuffer &image, diff --git a/image_analysis/spot_finding/ImageSpotFinder.h b/image_analysis/spot_finding/ImageSpotFinder.h index 14ec681bf..801988a68 100644 --- a/image_analysis/spot_finding/ImageSpotFinder.h +++ b/image_analysis/spot_finding/ImageSpotFinder.h @@ -7,6 +7,7 @@ #include #include "../../common/DiffractionSpot.h" +#include "../../common/PanelBoxes.h" #include "../image_preprocessing/ImagePreprocessorBuffer.h" @@ -26,6 +27,13 @@ protected: // The connected components of the last extraction. A member so ExtractComponents can hand out a // reference and reuse the allocation from image to image. std::vector components; + // The panels the image is made of (PanelBoxes): the local windows and the connected components stay + // inside a panel and every spot carries its panel - the FPGA spot finder's model, per JUNGFRAU module. + // Default, and for an empty module table: one panel, the whole image. + std::vector panels; + // Pixels on no panel, packed like res_mask_bits; empty when every pixel is on one. Folded into the + // resolution mask, so such a pixel never reaches a component. + std::vector off_panel_bits; // host_bit_buffer = false leaves output_buffer empty: the GPU finders extract on the device and // never read the bit buffer on the host, so allocating and pinning 2.26 MB per engine (at 18 MP) @@ -65,6 +73,10 @@ public: // resolution map rather than by each of them (AzimuthalIntegrationMapping::ResolutionMaskBits). virtual void SetResolutionMaskBits(const std::vector &packed_mask); + // The panels of the image, from the experiment's module table. Set once, before the first image; + // the resolution mask is set again after it (SetResolutionMaskBits folds the panels in). + virtual void SetPanels(const DetectorModules &modules); + // Every connected component of the last Detect() with at most max-pix pixels. min-pix is NOT // applied here on purpose - it is the only spot setting that changes between the passes of the // per-image min-pix search, so ONE extraction serves all three of them. diff --git a/image_analysis/spot_finding/ImageSpotFinderCPU.cpp b/image_analysis/spot_finding/ImageSpotFinderCPU.cpp index 4062f5b10..e97d3eafb 100644 --- a/image_analysis/spot_finding/ImageSpotFinderCPU.cpp +++ b/image_analysis/spot_finding/ImageSpotFinderCPU.cpp @@ -81,11 +81,12 @@ void ImageSpotFinderCPU::DetectAt(const ImagePreprocessorBuffer &image, const Sp for (const auto &c : candidate_windows) { const int32_t line = c.pxl / width; const int32_t col = c.pxl % width; + const PanelBox &b = *c.panel; int64_t sum = c.sum, sum2 = c.sum2, valid = c.valid; // Take out of the window what the second pass does not count: the first pass's strong pixels. - for (int32_t y = std::max(line - NBX, 0); y <= std::min(line + NBX, height - 1); y++) { - const int32_t first = y * width + std::max(col - NBX, 0); - const int32_t last = y * width + std::min(col + NBX, width - 1); + for (int32_t y = std::max(line - NBX, b.y0); y <= std::min(line + NBX, b.y1 - 1); y++) { + const int32_t first = y * width + std::max(col - NBX, b.x0); + const int32_t last = y * width + std::min(col + NBX, b.x1 - 1); for (int32_t w = first / 32; w <= last / 32; w++) { uint32_t bits = first_pass_buffer[w]; if (w == first / 32) @@ -119,34 +120,6 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, for (int i = 0; i < OutputSize(); i++) out_buffer[i] = 0; - // The first pass's bits are read only inside the window of a candidate (and at the candidate), so - // with candidates it tests only the pixels of the row/column blocks such a window reaches; the bits - // it leaves unset there are never read. Those blocks are marked from each row's candidates as the - // row enters the vertical sums below, which is before any row within NBX of it is tested - so the - // candidates themselves can be filled in a row at a time (fill_row), while the row is in cache. - const int32_t nblocks = (width + 31) / 32; - const auto new_row = [&](int32_t r) { - if (!candidates) - return; - if (fill_row) - fill_row(r); - const int32_t first = r * width; - const int32_t last = first + width - 1; - for (int32_t w = first / 32; w <= last / 32; w++) { - uint32_t bits = candidates[w]; - if (w == first / 32) - bits &= UINT32_MAX << (first % 32); - if (w == last / 32 && last % 32 != 31) - bits &= (1U << (last % 32 + 1)) - 1; - for (; bits; bits &= bits - 1) { - const int32_t col = w * 32 + std::countr_zero(bits) - first; - for (int32_t y = std::max(r - NBX, 0); y <= std::min(r + NBX, height - 1); y++) - for (int32_t b = std::max(col - NBX, 0) / 32; b <= std::min(col + NBX, width - 1) / 32; b++) - first_pass_needed[static_cast(y) * nblocks + b] = 1; - } - } - }; - // A pixel found strong by the previous pass reads as INT32_MAX, which the accumulation below // already skips and the acceptance test below already takes as strong - the same substitution // the GPU kernel makes when it reads prev_out. @@ -156,50 +129,98 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, return image[pxl]; }; + // Panel by panel (one panel, the whole image, without a module table): every window, sum and + // run below is clipped to the panel's box, so nothing on one panel is measured against another. + // The bits of a row are gathered 32 at a time and or-ed in at a word end and at the panel's edge. std::bitset<32> out = 0; + const auto put = [&](int32_t pxl, bool strong, bool row_end) { + if (strong) + out.set(pxl % 32); + if (pxl % 32 == 31 || row_end) { + out_buffer[pxl / 32] |= out.to_ulong(); + out.reset(); + } + }; if (settings.signal_to_noise_threshold <= 0.0) { if (settings.photon_count_threshold > 0) { - for (int pxl = 0; pxl < height * width; pxl++) { - int32_t bit = pxl % 32; - int32_t pxl_val = value_at(pxl); - if (pxl_val == INT32_MAX || (pxl_val > settings.photon_count_threshold && pxl_val != INT32_MIN)) - out.set(bit); + for (const auto &b : panels) + for (int32_t line = b.y0; line < b.y1; line++) + for (int32_t col = b.x0; col < b.x1; col++) { + const int32_t pxl = line * width + col; + const int32_t pxl_val = value_at(pxl); + put(pxl, pxl_val == INT32_MAX || (pxl_val > settings.photon_count_threshold && pxl_val != INT32_MIN), + col == b.x1 - 1); + } + } + return; + } - if (bit == 31) { - out_buffer[pxl / 32] = out.to_ulong(); - out.reset(); + float strong2 = settings.signal_to_noise_threshold * settings.signal_to_noise_threshold; + + // Sum and sum of squares of (2*NBY+1) vertical elements + // These are updated after each line is finished + // 64-bit integer guarantees calculations are made without rounding errors + std::vector sum_vert(width, 0); + std::vector sum2_vert(width, 0); + std::vector valid_vert(width, 0); + 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(); + const std::vector no_row(width, INT32_MIN); + const int32_t nblocks = (width + 31) / 32; + // The rows fill_row has been asked for; a row shared by panels side by side is filled once. + std::vector row_filled(fill_row ? height : 0, 0); + + for (const auto &b : panels) { + const int32_t bw = b.x1 - b.x0; + + // The first pass's bits are read only inside the window of a candidate (and at the candidate), + // so with candidates it tests only the pixels of the row/column blocks such a window reaches; + // the bits it leaves unset there are never read. Those blocks are marked from each row's + // candidates as the row enters the vertical sums below, which is before any row within NBX of + // it is tested - so the candidates themselves can be filled in a row at a time (fill_row), + // while the row is in cache. + const auto new_row = [&](int32_t r) { + if (!candidates) + return; + if (fill_row && !row_filled[r]) { + fill_row(r); + row_filled[r] = 1; + } + const int32_t first = r * width + b.x0; + const int32_t last = r * width + b.x1 - 1; + for (int32_t word = first / 32; word <= last / 32; word++) { + uint32_t bits = candidates[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 col = word * 32 + std::countr_zero(bits) - r * width; + for (int32_t y = std::max(r - NBX, b.y0); y <= std::min(r + NBX, b.y1 - 1); y++) + for (int32_t k = std::max(col - NBX, b.x0) / 32; k <= std::min(col + NBX, b.x1 - 1) / 32; k++) + first_pass_needed[static_cast(y) * nblocks + k] = 1; } } - } - } else { - float strong2 = settings.signal_to_noise_threshold * settings.signal_to_noise_threshold; - - // Sum and sum of squares of (2*NBY+1) vertical elements - // These are updated after each line is finished - // 64-bit integer guarantees calculations are made without rounding errors - std::vector sum_vert(width, 0); - std::vector sum2_vert(width, 0); - std::vector valid_vert(width, 0); + }; // Add the row entering the window (line_in) and take out the one leaving it (line_out); -1 is // no row. A pixel strong in the previous pass is not counted either; SlideVertical adds it with // the rest of its row and it is 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(); - const std::vector no_row(width, INT32_MIN); auto slide_vert = [&](int line_in, int line_out) { - SlideVertical(line_in >= 0 ? img + static_cast(line_in) * w : no_row.data(), - line_out >= 0 ? img + static_cast(line_out) * w : no_row.data(), w, sv, sv2, vv); + SlideVertical(line_in >= 0 ? img + static_cast(line_in) * w + b.x0 : no_row.data(), + line_out >= 0 ? img + static_cast(line_out) * w + b.x0 : no_row.data(), bw, + sv + b.x0, sv2 + b.x0, vv + b.x0); 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; + const int32_t row = line * w; + const int32_t first = row + b.x0, last = row + b.x1 - 1; for (int32_t word = first / 32; word <= last / 32; word++) { uint32_t bits = prev_strong[word]; if (word == first / 32) @@ -211,9 +232,9 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, 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); + sv[pxl - row] -= sign * v; + sv2[pxl - row] -= sign * v * v; + vv[pxl - row] -= static_cast(sign); } } }; @@ -223,15 +244,19 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, uncount_strong(line_out, -1); }; - for (int line = 0; line < NBX; line++) { + std::fill(sum_vert.begin() + b.x0, sum_vert.begin() + b.x1, 0); + std::fill(sum2_vert.begin() + b.x0, sum2_vert.begin() + b.x1, 0); + std::fill(valid_vert.begin() + b.x0, valid_vert.begin() + b.x1, 0); + + for (int line = b.y0; line < std::min(b.y0 + NBX, b.y1); line++) { new_row(line); slide_vert(line, -1); } - for (int line = 0; line < height; line++) { - if (line < height - NBX) + for (int line = b.y0; line < b.y1; line++) { + if (line + NBX < b.y1) new_row(line + NBX); - slide_vert(line < height - NBX ? line + NBX : -1, line >= NBX + 1 ? line - (NBX + 1) : -1); + slide_vert(line + NBX < b.y1 ? line + NBX : -1, line - (NBX + 1) >= b.y0 ? line - (NBX + 1) : -1); if (candidates) { // Only the runs of 32-column blocks whose result is read (first_pass_needed; every @@ -239,25 +264,25 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, // columns it covers there - integers, so the same sums the running window holds at // that column - and slides it as below; the bits outside the runs stay 0, as below. const uint8_t *needed = first_pass_needed.data() + static_cast(line) * nblocks; - for (int32_t b = 0; b < nblocks;) { - if (!needed[b]) { b++; continue; } - const int32_t c0 = b * 32; - while (b < nblocks && needed[b]) b++; - const int32_t c1 = std::min(b * 32, width); + for (int32_t k = b.x0 / 32; k <= (b.x1 - 1) / 32;) { + if (!needed[k]) { k++; continue; } + const int32_t c0 = std::max(k * 32, b.x0); + while (k <= (b.x1 - 1) / 32 && needed[k]) k++; + const int32_t c1 = std::min(k * 32, b.x1); int64_t sum = 0, sum2 = 0, valid = 0; - for (int32_t c = std::max(c0 - NBX, 0); c <= std::min(c0 + NBX, width - 1); c++) { + for (int32_t c = std::max(c0 - NBX, b.x0); c <= std::min(c0 + NBX, b.x1 - 1); c++) { sum += sum_vert[c]; sum2 += sum2_vert[c]; valid += valid_vert[c]; } for (int32_t col = c0; col < c1; col++) { if (col > c0) { - if (col < width - NBX) { + if (col + NBX < b.x1) { sum += sum_vert[col + NBX]; sum2 += sum2_vert[col + NBX]; valid += valid_vert[col + NBX]; } - if (col >= NBX + 1) { + if (col - NBX - 1 >= b.x0) { sum -= sum_vert[col - NBX - 1]; sum2 -= sum2_vert[col - NBX - 1]; valid -= valid_vert[col - NBX - 1]; @@ -265,7 +290,7 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, } const int32_t pxl = line * width + col; if (candidates[pxl / 32] & (1U << (pxl % 32))) - candidate_windows.push_back({pxl, sum, sum2, valid}); + candidate_windows.push_back({pxl, sum, sum2, valid, &b}); if (StrongInWindow(value_at(pxl), sum, sum2, valid, settings, strong2)) out_buffer[pxl / 32] |= 1U << (pxl % 32); } @@ -277,38 +302,27 @@ void ImageSpotFinderCPU::DetectPass(const ImagePreprocessorBuffer &image, int64_t sum2 = 0; int64_t valid = 0; - for (int col = 0; col < NBX; col++) { + for (int col = b.x0; col < std::min(b.x0 + NBX, b.x1); col++) { sum += sum_vert[col]; sum2 += sum2_vert[col]; valid += valid_vert[col]; } - for (int col = 0; col < width; col++) { - if (col < width - NBX) { + for (int col = b.x0; col < b.x1; col++) { + if (col + NBX < b.x1) { sum += sum_vert[col + NBX]; sum2 += sum2_vert[col + NBX]; valid += valid_vert[col + NBX]; } - if (col >= NBX + 1) { + if (col - NBX - 1 >= b.x0) { sum -= sum_vert[col - NBX - 1]; sum2 -= sum2_vert[col - NBX - 1]; valid -= valid_vert[col - NBX - 1]; } const int32_t pxl = line * width + col; - const int32_t bit = pxl % 32; - - if (StrongInWindow(value_at(pxl), sum, sum2, valid, settings, strong2)) - out.set(bit); - - if (bit == 31) { - out_buffer[pxl / 32] = out.to_ulong(); - out.reset() ; - } + put(pxl, StrongInWindow(value_at(pxl), sum, sum2, valid, settings, strong2), col == b.x1 - 1); } } } - - if (height * width % 32 != 0) - out_buffer[OutputSize() - 1] |= out.to_ulong(); // zeroed above; the candidates path sets its bits directly } diff --git a/image_analysis/spot_finding/ImageSpotFinderCPU.h b/image_analysis/spot_finding/ImageSpotFinderCPU.h index 5b754b2da..37a51430a 100644 --- a/image_analysis/spot_finding/ImageSpotFinderCPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderCPU.h @@ -37,6 +37,7 @@ class ImageSpotFinderCPU : public ImageSpotFinder { struct CandidateWindow { int32_t pxl; int64_t sum, sum2, valid; + const PanelBox *panel; // the window is clipped to it }; std::vector candidate_windows; // Where the first pass's own result is read by DetectAt - within NBX of a candidate - per row and diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.cu b/image_analysis/spot_finding/ImageSpotFinderGPU.cu index 4d8ece36c..6d836f4c0 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.cu +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.cu @@ -4,6 +4,8 @@ // GPU Spot finding developed by Hans-Christian Stadler (PSI) // Copyright (2019-2023) Paul Scherrer Institute +#include + #include "ImageSpotFinderGPU.h" #include "../../common/JFJochException.h" @@ -135,36 +137,48 @@ __device__ __forceinline__ T warp_window_sum(T a, T b) { // read 31 rows earlier and sits in cache) and the row entering it is added. The horizontal window // sums come from warp prefix scans, so there is no shared memory and no block-wide synchronisation. // All sums are integer (modular 64-bit), so the result is exact and independent of the order. +// +// Panel by panel (one panel, the whole image, without a module table): the warps are numbered panel after +// panel, each panel tiled the same way, and every window is clipped to the warp's panel - the FPGA's spot +// finder works module by module the same way. Warp-uniform, so the shuffles and ballots stay collective. __global__ void analyze_pixel(const int32_t *in, const uint32_t *prev_out, uint32_t *out, - const uint32_t *needed, const spot_parameters params, int32_t rowsPerWave) + const uint32_t *needed, const spot_parameters params, int32_t rowsPerWave, + const PanelBox *panels, int32_t npanels) { constexpr int32_t NBX = ImageSpotFinder::NBX; const int32_t lane = threadIdx.x & (WARP_SIZE - 1); - const int32_t warp = (blockIdx.x * blockDim.x + threadIdx.x) / WARP_SIZE; - const int32_t ngroups = (params.width + WARP_SIZE - 1) / WARP_SIZE; - const int32_t nWaves = (params.height + rowsPerWave - 1) / rowsPerWave; + int32_t warp = (blockIdx.x * blockDim.x + threadIdx.x) / WARP_SIZE; + int32_t p = 0; + for (; p < npanels; p++) { + const int32_t n = (panels[p].x1 - panels[p].x0 + WARP_SIZE - 1) / WARP_SIZE + * ((panels[p].y1 - panels[p].y0 + rowsPerWave - 1) / rowsPerWave); + if (warp < n) + break; + warp -= n; + } + if (p == npanels) + return; + const PanelBox b = panels[p]; + const int32_t ngroups = (b.x1 - b.x0 + WARP_SIZE - 1) / WARP_SIZE; const int32_t group = warp % ngroups; const int32_t wave = warp / ngroups; - const int32_t rmin = wave * rowsPerWave; - const int32_t rmax = min(rmin + rowsPerWave, params.height); - // Warp-uniform, so the shuffles and ballots below stay collective. - if (wave >= nWaves) - return; + const int32_t rmin = b.y0 + wave * rowsPerWave; + const int32_t rmax = min(rmin + rowsPerWave, b.y1); - const int32_t out_col = group * WARP_SIZE + lane; // the column this lane writes - const int32_t col_a = group * WARP_SIZE - NBX - 1 + lane; // the two input columns this lane sums + const int32_t out_col = b.x0 + group * WARP_SIZE + lane; // the column this lane writes + const int32_t col_a = b.x0 + group * WARP_SIZE - NBX - 1 + lane; // the two input columns this lane sums const int32_t col_b = col_a + WARP_SIZE; - const bool use_a = col_a >= 0 && col_a < params.width; - const bool use_b = lane < WARP_SIZE - 1 && col_b >= 0 && col_b < params.width; + const bool use_a = col_a >= b.x0 && col_a < b.x1; + const bool use_b = lane < WARP_SIZE - 1 && col_b >= b.x0 && col_b < b.x1; if (needed != nullptr) { // Is any needed pixel within NBX of the tile, rows [rmin, rmax) x columns [c0, c0 + 32)? // Each lane checks rows of the dilated band, the columns [c0 - NBX, c0 + 32 + NBX) as the // (at most three) words they fall in, masked to exactly that range. - const int64_t c_lo = max(group * WARP_SIZE - NBX, 0); - const int64_t c_hi = min(group * WARP_SIZE + WARP_SIZE + NBX, params.width); // past the last + const int64_t c_lo = max(b.x0 + group * WARP_SIZE - NBX, b.x0); + const int64_t c_hi = min(b.x0 + group * WARP_SIZE + WARP_SIZE + NBX, b.x1); // past the last bool any = false; - for (int32_t r = max(rmin - NBX, 0) + lane; r < min(rmax + NBX, params.height) && !any; r += WARP_SIZE) { + for (int32_t r = max(rmin - NBX, b.y0) + lane; r < min(rmax + NBX, b.y1) && !any; r += WARP_SIZE) { const int64_t first = static_cast(r) * params.width + c_lo; const int64_t last = static_cast(r) * params.width + c_hi - 1; for (int64_t w = first / 32; w <= last / 32; w++) { @@ -203,7 +217,7 @@ __global__ void analyze_pixel(const int32_t *in, const uint32_t *prev_out, uint3 } }; - for (int32_t r = max(rmin - NBX, 0); r < min(rmin + NBX, params.height - 1) + 1; r++) { + for (int32_t r = max(rmin - NBX, b.y0); r < min(rmin + NBX, b.y1 - 1) + 1; r++) { slide(read_a(r), 1, sum_a, sum2_a, count_a); slide(read_b(r), 1, sum_b, sum2_b, count_b); } @@ -211,19 +225,19 @@ __global__ void analyze_pixel(const int32_t *in, const uint32_t *prev_out, uint3 for (int32_t row = rmin; row < rmax; row++) { // The rows leaving and entering the window are read before the window sums are formed, so the // loads are in flight while the scans run rather than after them. - const bool leave = row - NBX >= 0; - const bool enter = row + NBX + 1 < params.height; + const bool leave = row - NBX >= b.y0; + const bool enter = row + NBX + 1 < b.y1; const int32_t leave_a = leave ? read_a(row - NBX) : INT32_MIN; const int32_t leave_b = leave ? read_b(row - NBX) : INT32_MIN; const int32_t enter_a = enter ? read_a(row + NBX + 1) : INT32_MIN; const int32_t enter_b = enter ? read_b(row + NBX + 1) : INT32_MIN; - const int32_t value = (out_col < params.width) ? value_at(row, out_col) : INT32_MIN; + const int32_t value = (out_col < b.x1) ? value_at(row, out_col) : INT32_MIN; const auto sum = static_cast(warp_window_sum(sum_a, sum_b)); const auto sum2 = static_cast(warp_window_sum(sum2_a, sum2_b)); const auto count = static_cast(warp_window_sum(count_a, count_b)); uint8_t val = 0; - if (out_col < params.width) + if (out_col < b.x1) val = pixel_result(params, value, sum, sum2, count); write_result(params, out, row * params.width + out_col, val); @@ -245,8 +259,10 @@ __global__ void analyze_pixel(const int32_t *in, const uint32_t *prev_out, uint3 // One warp per candidate: lane dx (0..30) sums column c + dx - NBX over the window's rows, and the // warp reduces. The candidate words are scanned 32 at a time, one per lane; a word is owned by one // warp, which writes it whole, so out needs no atomics (it must be zeroed beforehand). +// The window is clipped to the candidate's panel, as in analyze_pixel; a candidate on no panel is not strong. __global__ void analyze_candidates(const int32_t *in, const uint32_t *prev_out, const uint32_t *cand, - uint32_t *out, const spot_parameters params, size_t nwords) + uint32_t *out, const spot_parameters params, size_t nwords, + const PanelBox *panels, int32_t npanels) { constexpr unsigned ALL_THREADS = UINT32_MAX; constexpr int NBX = ImageSpotFinder::NBX; @@ -277,10 +293,17 @@ __global__ void analyze_candidates(const int32_t *in, const uint32_t *prev_out, break; const int32_t row = pixel / params.width; const int32_t col = pixel % params.width + lane - NBX; + int32_t p = 0; + while (p < npanels && !(pixel % params.width >= panels[p].x0 && pixel % params.width < panels[p].x1 + && row >= panels[p].y0 && row < panels[p].y1)) + p++; + if (p == npanels) + continue; + const PanelBox box = panels[p]; int64_t sum = 0, sum2 = 0; int32_t count = 0; - if (lane < 2 * NBX + 1 && col >= 0 && col < params.width) { + if (lane < 2 * NBX + 1 && col >= box.x0 && col < box.x1) { // Every read of the column first, all of them unconditional (rows past an edge are // read at the edge and dropped below), so the 2 * (2 * NBX + 1) loads are issued // together. Read through value_at, each pixel's load waited for its previous-pass @@ -290,7 +313,7 @@ __global__ void analyze_candidates(const int32_t *in, const uint32_t *prev_out, uint32_t prev[2 * NBX + 1]; #pragma unroll for (int k = 0; k < 2 * NBX + 1; k++) { - const int32_t r = min(max(row + k - NBX, 0), params.height - 1); + const int32_t r = min(max(row + k - NBX, box.y0), box.y1 - 1); const int32_t pix = r * params.width + col; vals[k] = in[pix]; prev[k] = prev_out[pix / 32] >> (pix % 32); @@ -300,7 +323,7 @@ __global__ void analyze_candidates(const int32_t *in, const uint32_t *prev_out, const int32_t r = row + k - NBX; // A previous-pass pixel reads as INT32_MAX, so like the sentinels it is not counted. const int32_t val = vals[k]; - if (r >= 0 && r < params.height && (prev[k] & 1u) == 0 + if (r >= box.y0 && r < box.y1 && (prev[k] & 1u) == 0 && val != INT32_MAX && val != INT32_MIN) { sum += val; sum2 += static_cast(val) * val; @@ -329,6 +352,7 @@ ImageSpotFinderGPU::ImageSpotFinderGPU(int32_t in_width, int32_t in_height, extractor(in_width, in_height, std::move(in_stream)) { gpu_out_0 = CudaDevicePtr(OutputSize()); gpu_out_1 = CudaDevicePtr(OutputSize()); + UploadPanels(); int device; cuda_err(cudaGetDevice(&device)); int sm_count; @@ -336,6 +360,19 @@ ImageSpotFinderGPU::ImageSpotFinderGPU(int32_t in_width, int32_t in_height, candidateBlocks = 8 * sm_count; } +void ImageSpotFinderGPU::UploadPanels() { + gpu_panels = CudaDevicePtr(panels.size()); + cuda_err(cudaMemcpyAsync(gpu_panels, panels.data(), panels.size() * sizeof(PanelBox), + cudaMemcpyHostToDevice, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); +} + +void ImageSpotFinderGPU::SetPanels(const DetectorModules &modules) { + ImageSpotFinder::SetPanels(modules); + extractor.SetPanels(panels); + UploadPanels(); +} + void ImageSpotFinderGPU::SetResolutionMaskBits(const std::vector &packed_mask) { ImageSpotFinder::SetResolutionMaskBits(packed_mask); extractor.SetResolutionMask(res_mask_bits); @@ -344,6 +381,11 @@ void ImageSpotFinderGPU::SetResolutionMaskBits(const std::vector &pack const std::vector &ImageSpotFinderGPU::ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) { extractor.Extract(gpu_out_1, image.getGPUBuffer(), settings, components); + // The device labels the whole image at once, so its components come in the order of their first + // pixel; the host labels panel by panel. Put them panel by panel too (a no-op with one panel), and + // the two lists are the same list. + std::stable_sort(components.begin(), components.end(), + [](const DiffractionSpot &a, const DiffractionSpot &b) { return a.Panel() < b.Panel(); }); return components; } @@ -383,22 +425,25 @@ void ImageSpotFinderGPU::RunDetect(const ImagePreprocessorBuffer &image, const S // every tile is worked, and long tiles amortise the 2 * NBX rows each one spends filling its window. const int32_t rowsPerWave = (gpu_candidates == nullptr) ? (spot_params.height + numberOfWaves - 1) / numberOfWaves : sparseRowsPerWave; - const int32_t nWaves = (spot_params.height + rowsPerWave - 1) / rowsPerWave; - const int32_t ngroups = (spot_params.width + 31) / 32; + int32_t nWarps = 0; + for (const auto &b : panels) + nWarps += (b.x1 - b.x0 + 31) / 32 * ((b.y1 - b.y0 + rowsPerWave - 1) / rowsPerWave); const int32_t warpsPerBlock = numberOfCudaThreads / 32; - const int32_t nBlocks = (ngroups * nWaves + warpsPerBlock - 1) / warpsPerBlock; + const int32_t nBlocks = (nWarps + warpsPerBlock - 1) / warpsPerBlock; + const auto npanels = static_cast(panels.size()); cuda_err(cudaMemsetAsync(gpu_out_0, 0, OutputSize() * sizeof(uint32_t), *stream)); cuda_err(cudaMemsetAsync(gpu_out_1, 0, OutputSize() * sizeof(uint32_t), *stream)); // The first pass has no previous one (the classic engine used to hand it the zeroed gpu_out_1). analyze_pixel<<>> - (image.getGPUBuffer(), nullptr, gpu_out_0, gpu_candidates, spot_params, rowsPerWave); + (image.getGPUBuffer(), nullptr, gpu_out_0, gpu_candidates, spot_params, rowsPerWave, gpu_panels, npanels); cuda_err(cudaGetLastError()); if (gpu_candidates == nullptr) analyze_pixel<<>> - (image.getGPUBuffer(), gpu_out_0, gpu_out_1, nullptr, spot_params, rowsPerWave); + (image.getGPUBuffer(), gpu_out_0, gpu_out_1, nullptr, spot_params, rowsPerWave, gpu_panels, npanels); else analyze_candidates<<>> - (image.getGPUBuffer(), gpu_out_0, gpu_candidates, gpu_out_1, spot_params, OutputSize()); + (image.getGPUBuffer(), gpu_out_0, gpu_candidates, gpu_out_1, spot_params, OutputSize(), + gpu_panels, npanels); cuda_err(cudaGetLastError()); } diff --git a/image_analysis/spot_finding/ImageSpotFinderGPU.h b/image_analysis/spot_finding/ImageSpotFinderGPU.h index f363fb5ca..fd7be7e7b 100644 --- a/image_analysis/spot_finding/ImageSpotFinderGPU.h +++ b/image_analysis/spot_finding/ImageSpotFinderGPU.h @@ -21,6 +21,7 @@ protected: CudaDevicePtr gpu_out_0; CudaDevicePtr gpu_out_1; // holds the strong-pixel bits after Detect() SpotExtractorGPU extractor; + CudaDevicePtr gpu_panels; // ImageSpotFinder::panels on the device private: const int numberOfCudaThreads = 128; // #threads per block of analyze_pixel (one warp per 32 columns) @@ -29,6 +30,8 @@ private: const int windowSizeLimit = 32; // limit on the window size (2nby+1, 2nbx+1): a warp holds 2 x 32 columns int candidateBlocks = 0; // grid of analyze_candidates (256 threads per block) + void UploadPanels(); + void RunDetect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings, const uint32_t *gpu_candidates); protected: @@ -44,6 +47,7 @@ public: void Detect(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; void SetResolutionMaskBits(const std::vector &packed_mask) override; + void SetPanels(const DetectorModules &modules) override; [[nodiscard]] uint32_t StrongPixelCount() const override { return extractor.StrongPixelCount(); } const std::vector &ExtractComponents(const ImagePreprocessorBuffer &image, const SpotFindingSettings &settings) override; diff --git a/image_analysis/spot_finding/SpotExtractorGPU.cu b/image_analysis/spot_finding/SpotExtractorGPU.cu index b91e97205..6b2e40637 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.cu +++ b/image_analysis/spot_finding/SpotExtractorGPU.cu @@ -154,23 +154,38 @@ __global__ void init_parent(uint32_t *__restrict__ parent, const uint32_t *__res parent[i] = i; } +// The panel (index into panels) holding pixel `flat`, or npanels for none - the resolution mask keeps such +// a pixel out of the list, so that does not happen. +__device__ __forceinline__ int panel_of(const PanelBox *panels, int npanels, uint32_t flat, int width) { + const int col = static_cast(flat % width), line = static_cast(flat / width); + int p = 0; + while (p < npanels && !(col >= panels[p].x0 && col < panels[p].x1 && line >= panels[p].y0 && line < panels[p].y1)) + p++; + return p; +} + // The list is sorted by flat index, so a pixel's 8-neighbours that come EARLIER in it are exactly // four: (line, col-1), (line-1, col-1), (line-1, col) and (line-1, col+1). Each is one binary -// search away, which is what makes the sparse formulation cheap. +// search away, which is what makes the sparse formulation cheap. Only those on the pixel's own panel +// count: a component never joins two panels (the host labels panel by panel). __global__ void union_neighbours(const uint32_t *__restrict__ index, uint32_t *__restrict__ parent, - const uint32_t *__restrict__ nstrong, uint32_t capacity, int width) { + const uint32_t *__restrict__ nstrong, uint32_t capacity, int width, + const PanelBox *__restrict__ panels, int npanels) { const int n = static_cast(min(*nstrong, capacity)); if (static_cast(n) >= capacity) return; for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) { const uint32_t flat = index[i]; - const uint32_t col = flat % width; - const uint32_t line = flat / width; + const int p = panel_of(panels, npanels, flat, width); + if (p == npanels) continue; + const int col = static_cast(flat % width); + const int line = static_cast(flat / width); + const PanelBox b = panels[p]; uint32_t candidate[4]; int ncandidate = 0; - if (col > 0) candidate[ncandidate++] = flat - 1; - if (line > 0 && col > 0) candidate[ncandidate++] = flat - width - 1; - if (line > 0) candidate[ncandidate++] = flat - width; - if (line > 0 && col + 1 < static_cast(width)) candidate[ncandidate++] = flat - width + 1; + if (col > b.x0) candidate[ncandidate++] = flat - 1; + if (line > b.y0 && col > b.x0) candidate[ncandidate++] = flat - width - 1; + if (line > b.y0) candidate[ncandidate++] = flat - width; + if (line > b.y0 && col + 1 < b.x1) candidate[ncandidate++] = flat - width + 1; for (int k = 0; k < ncandidate; k++) { const int j = lower_bound_device(index, i, candidate[k]); if (j < i && index[j] == candidate[k]) merge_roots(parent, static_cast(i), static_cast(j)); @@ -234,6 +249,7 @@ __global__ void label_components(const uint32_t *__restrict__ root, uint32_t *__ s.max_photons = LLONG_MIN; s.pixel_count = 0; s.bbox_side = 0; + s.panel = -1; } } @@ -251,7 +267,8 @@ __global__ void init_boxes(int4 *__restrict__ box, const uint32_t *__restrict__ __global__ void sum_components(const uint32_t *__restrict__ index, const int32_t *__restrict__ value, const uint32_t *__restrict__ root, const uint32_t *__restrict__ label, SpotExtractorGPUSpot *__restrict__ scratch, int4 *__restrict__ box, - const uint32_t *__restrict__ nstrong, uint32_t capacity, int width) { + const uint32_t *__restrict__ nstrong, uint32_t capacity, int width, + const PanelBox *__restrict__ panels, int npanels) { const int n = static_cast(min(*nstrong, capacity)); if (static_cast(n) >= capacity) return; for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += blockDim.x * gridDim.x) { @@ -268,6 +285,10 @@ __global__ void sum_components(const uint32_t *__restrict__ index, const int32_t atomicMax(&box[l].y, col); atomicMin(&box[l].z, line); atomicMax(&box[l].w, line); + // The same panel for every pixel of the component, so any of them may write it. + const int p = panel_of(panels, npanels, index[i], width); + if (p < npanels) + s.panel = panels[p].panel; } } @@ -351,6 +372,15 @@ SpotExtractorGPU::SpotExtractorGPU(int32_t in_width, int32_t in_height, std::sha cuda_err(cudaMemcpyAsync(gpu_res_mask, mask.data(), nwords * sizeof(uint32_t), cudaMemcpyHostToDevice, *stream)); cuda_err(cudaStreamSynchronize(*stream)); + SetPanels(PanelBoxes({}, in_width, in_height)); +} + +void SpotExtractorGPU::SetPanels(const std::vector &panels) { + npanels = static_cast(panels.size()); + gpu_panels = CudaDevicePtr(panels.size()); + cuda_err(cudaMemcpyAsync(gpu_panels, panels.data(), panels.size() * sizeof(PanelBox), + cudaMemcpyHostToDevice, *stream)); + cuda_err(cudaStreamSynchronize(*stream)); } void SpotExtractorGPU::SetResolutionMask(const std::vector &packed_mask) { @@ -378,7 +408,8 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im // never has to synchronise in the middle of the frame. init_parent<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_nstrong, max_strong); cuda_err(cudaGetLastError()); - union_neighbours<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_parent, gpu_nstrong, max_strong, width); + union_neighbours<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_parent, gpu_nstrong, max_strong, width, + gpu_panels, npanels); cuda_err(cudaGetLastError()); resolve_roots<<<512, THREADS, 0, *stream>>>(gpu_parent, gpu_root, gpu_nstrong, max_strong); cuda_err(cudaGetLastError()); @@ -388,7 +419,7 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im init_boxes<<<512, THREADS, 0, *stream>>>(gpu_box, gpu_nlabel); cuda_err(cudaGetLastError()); sum_components<<<512, THREADS, 0, *stream>>>(gpu_index, gpu_value, gpu_root, gpu_label, gpu_spot, gpu_box, - gpu_nstrong, max_strong, width); + gpu_nstrong, max_strong, width, gpu_panels, npanels); cuda_err(cudaGetLastError()); filter_components<<<1, FINISH_THREADS, 0, *stream>>>(gpu_spot, gpu_box, gpu_spot_out, gpu_nspot, gpu_nlabel, max_pix, static_cast(SPOT_SHAPE_FREE_PIXELS), @@ -415,5 +446,5 @@ void SpotExtractorGPU::Extract(const uint32_t *gpu_strong, const int32_t *gpu_im spots.clear(); spots.reserve(nspot); for (uint32_t i = 0; i < nspot; i++) - spots.emplace_back(s[i].x, s[i].y, s[i].pixel_count, s[i].photons, s[i].max_photons); + spots.emplace_back(s[i].x, s[i].y, s[i].pixel_count, s[i].photons, s[i].max_photons, s[i].panel); } diff --git a/image_analysis/spot_finding/SpotExtractorGPU.h b/image_analysis/spot_finding/SpotExtractorGPU.h index 8639a46f9..7aebc2130 100644 --- a/image_analysis/spot_finding/SpotExtractorGPU.h +++ b/image_analysis/spot_finding/SpotExtractorGPU.h @@ -32,6 +32,7 @@ #include #include "../../common/DiffractionSpot.h" +#include "../../common/PanelBoxes.h" #include "../indexing/CUDAMemHelpers.h" #include "SpotFindingSettings.h" @@ -46,6 +47,7 @@ struct SpotExtractorGPUSpot { // Longer side of the component's bounding box, for SpotShapeAccepted. It fits in what used to be // padding, so carrying it costs nothing. int32_t bbox_side; + int32_t panel; // PanelBox::panel of the component's pixels }; class SpotExtractorGPU { @@ -77,6 +79,8 @@ class SpotExtractorGPU { CudaDevicePtr gpu_spot; CudaDevicePtr gpu_spot_out; CudaDevicePtr gpu_nspot; + CudaDevicePtr gpu_panels; // the panels components stay inside + int npanels = 0; CudaHostPtr host_nstrong; CudaHostPtr host_nspot; @@ -89,6 +93,8 @@ public: ~SpotExtractorGPU() { if (stream) cudaStreamSynchronize(*stream); } void SetResolutionMask(const std::vector &packed_mask); + // The panels of the image (PanelBoxes); default one, the whole image. + void SetPanels(const std::vector &panels); // gpu_strong is the finder's device bit buffer, gpu_image the preprocessed image it was built // from. Fills spots with every component of at most max-pix pixels, in the same order the host diff --git a/image_analysis/spot_finding/SpotUtils.cpp b/image_analysis/spot_finding/SpotUtils.cpp index 1ada9ee7f..2931a355c 100644 --- a/image_analysis/spot_finding/SpotUtils.cpp +++ b/image_analysis/spot_finding/SpotUtils.cpp @@ -372,12 +372,12 @@ void SpotAnalyze(const DiffractionExperiment &experiment, std::vector spots_out; - // A spot is found on the image; from here on it is where its module puts it on the detector plane. + // A spot is found on the image; from here on it is where its panel puts it on the detector plane. const auto &modules = experiment.GetModules(); for (const auto &spot: spots) { if (auto s = spot.Export(geom, output.number); s.has_value()) { if (!modules.empty()) { - const auto [px, py] = modules.ImageToPlane(s->x, s->y); + const auto [px, py] = modules.ImageToPlane(spot.Panel(), s->x, s->y); s->x = static_cast(px); s->y = static_cast(py); s->d_A = geom.PxlToRes(s->x, s->y); diff --git a/rugnux/DiagnosticOutput.cpp b/rugnux/DiagnosticOutput.cpp index 7083bccaa..9d6940405 100644 --- a/rugnux/DiagnosticOutput.cpp +++ b/rugnux/DiagnosticOutput.cpp @@ -8,6 +8,7 @@ #include +#include "../common/AssembleToPlane.h" #include "../common/ColorScale.h" #include "../common/CompressedImage.h" #include "../common/JFJochException.h" @@ -95,10 +96,16 @@ std::vector DetectorDiagnosticRGB(const PixelMask &pixel_mask, std::string RenderDetectorDiagnostic(const DiffractionExperiment &experiment, const PixelMask &pixel_mask, const std::vector &mean_projection) { + // The picture is of the detector as it is: with panels, the stored layout put together on the plane + // raster, the space between panels in the gap colour. Without, the image itself. + const DetectorModules &modules = experiment.GetModules(); + const int64_t width = experiment.GetXPixelsNumConv(), height = experiment.GetYPixelsNumConv(); + const PlaneRaster raster = PlaneRasterOf(modules, width, height); + const rgb gap_color = ColorScale().Apply(ColorScaleSpecial::Gap); // Named, because CompressedImage does not own the pixels it is handed. - const std::vector image = DetectorDiagnosticRGB(pixel_mask, mean_projection); - const CompressedImage rgb_image(image, experiment.GetXPixelsNumConv(), - experiment.GetYPixelsNumConv()); + const std::vector image = AssembleToPlane(modules, DetectorDiagnosticRGB(pixel_mask, mean_projection), + width, height, raster, gap_color); + const CompressedImage rgb_image(image, raster.width, raster.height); return WriteJPEGToMem(rgb_image, 85); } diff --git a/rugnux/RugnuxPreScan.cpp b/rugnux/RugnuxPreScan.cpp index 0ae59a893..9e26bd3fd 100644 --- a/rugnux/RugnuxPreScan.cpp +++ b/rugnux/RugnuxPreScan.cpp @@ -319,6 +319,12 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru std::unique_ptr preprocessed; std::unique_ptr spot_finder; }; + // The spot finder's windows and components, and the squares read about a spot, stay inside a panel. + const auto set_panels = [&](PreScanWorker &w) { + w.spot_finder->SetPanels(prescan_x.GetModules()); + w.preprocessed->SetPanels(PanelBoxes(prescan_x.GetModules(), prescan_x.GetXPixelsNumConv(), + prescan_x.GetYPixelsNumConv())); + }; const auto make_worker = [&] { PreScanWorker w; #ifdef JFJOCH_USE_CUDA @@ -338,6 +344,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru else w.spot_finder = std::make_unique(prescan_x.GetXPixelsNumConv(), prescan_x.GetYPixelsNumConv(), stream); + set_panels(w); return w; } #endif @@ -349,6 +356,7 @@ void Rugnux::PreScan(int start_image, int images_to_process, int frame_count, Ru else w.spot_finder = std::make_unique(prescan_x.GetXPixelsNumConv(), prescan_x.GetYPixelsNumConv()); + set_panels(w); } return w; }; diff --git a/tests/PanelLayoutTest.cpp b/tests/PanelLayoutTest.cpp index 4c8f5a484..963493a3b 100644 --- a/tests/PanelLayoutTest.cpp +++ b/tests/PanelLayoutTest.cpp @@ -13,9 +13,10 @@ // consumer retags its case [PanelLayout]: // stage 1 - dense maps: PanelLayout_AzintMapping, PanelLayout_ResolutionMask // stage 2 - sparse/Bragg: PanelLayout_BraggIntegration -// stage 3 - per-panel spot finding: PanelLayout_SpotsNeverCrossPanels; the raster algorithms -// (ShadowFinder, BeamCenterFFT, BeamCenterFromBackground, DiagnosticOutput) add their case here, -// comparing the result on the stacked layout with the result on the assembled image. +// stage 3 - per-panel spot finding: PanelLayout_SpotsNeverCrossPanels (+ PanelLayout_SpotFinderParity on a +// GPU); the raster algorithms on the plane raster: PanelLayout_AssembleToPlane, +// PanelLayout_ShadowFinder, PanelLayout_BeamCenterFFT, PanelLayout_DetectorDiagnostic - each +// comparing the stacked layout with the assembled image cut to the raster. #include @@ -23,6 +24,7 @@ #include #include +#include "../common/AssembleToPlane.h" #include "../common/AzimuthalIntegrationMapping.h" #include "../common/BraggIntegrationSettings.h" #include "../common/DetectorModules.h" @@ -30,10 +32,20 @@ #include "../common/DiffractionExperiment.h" #include "../common/PixelMask.h" #include "../common/Reflection.h" +#include "../common/CUDAWrapper.h" +#include "../image_analysis/beam_stop/ShadowFinder.h" #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h" +#include "../image_analysis/geom_refinement/BeamCenterFFT.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderCPU.h" #include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" #include "../image_analysis/spot_finding/SpotUtils.h" +#include "../rugnux/DiagnosticOutput.h" +#ifdef JFJOCH_USE_CUDA +#include "../image_analysis/image_preprocessing/ImagePreprocessorBufferGPU.h" +#include "../image_analysis/spot_finding/AdaptiveSpotFinderGPU.h" +#include "../image_analysis/spot_finding/ImageSpotFinderGPU.h" +#endif namespace { @@ -123,6 +135,7 @@ std::vector FindSpots(const DiffractionExperiment &x, const ImagePre settings.min_pix_per_spot = 1; settings.high_res_gap_Q_recipA = std::nullopt; ImageSpotFinderCPU finder(static_cast(x.GetXPixelsNum()), static_cast(x.GetYPixelsNum())); + finder.SetPanels(x.GetModules()); finder.SetResolutionMask(std::vector(x.GetPixelsNum(), false)); DataMessage msg{}; msg.number = 0; @@ -307,17 +320,238 @@ TEST_CASE("PanelLayout_BraggIntegration", "[.][PanelLayoutTodo]") { } } -// TODO stage 3: spot finding per panel. Two half spots on the facing edges of panels 0 and 1 as stored - -// adjacent rows of the stacked array, 139 px apart on the plane - stay two spots, and the local-background -// (SNR) finder agrees with the assembled image away from panel edges. -TEST_CASE("PanelLayout_SpotsNeverCrossPanels", "[.][PanelLayoutTodo]") { +// Spot finding per panel. Two half spots on the facing edges of panels 0 and 1 as stored - adjacent rows of +// the stacked array, 139 px apart on the plane - stay two spots, and the local-background (SNR) finder +// agrees with the assembled image: its windows stop at the panel edge, as they stop at the masked gap on +// the assembled image (16 px there, wider than the window's 15), instead of reaching into the other spot. +TEST_CASE("PanelLayout_SpotsNeverCrossPanels", "[PanelLayout]") { // Panel 0 stored row 63 and panel 1 stored row 0, both at stored column 30. const std::vector edge = {{40.0, 73.0}, {179.0, 73.0}}; const Scene scene(edge); - CheckSameSpots(FindSpots(scene.assembled, scene.assembled_image), - FindSpots(scene.stacked, scene.stacked_image)); + const auto a = FindSpots(scene.assembled, scene.assembled_image); + CHECK(a.size() == 2); + CheckSameSpots(a, FindSpots(scene.stacked, scene.stacked_image)); + const auto a_snr = FindSpots(scene.assembled, scene.assembled_image, 3.0f); + CHECK(a_snr.size() == 2); + CheckSameSpots(a_snr, FindSpots(scene.stacked, scene.stacked_image, 3.0f)); const Scene interior(INTERIOR_SPOTS); CheckSameSpots(FindSpots(interior.assembled, interior.assembled_image, 3.0f), FindSpots(interior.stacked, interior.stacked_image, 3.0f)); } + +#ifdef JFJOCH_USE_CUDA +// The GPU finders (local box, adaptive, and the device connected components) work panel by panel as the CPU +// ones do: on the stacked edge scene the same components, in the same order, on the same panels. +TEST_CASE("PanelLayout_SpotFinderParity", "[PanelLayout]") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping PanelLayout_SpotFinderParity"); + return; + } + const Scene scene({{40.0, 73.0}, {179.0, 73.0}, {50.3, 40.7}, {40.4, 140.8}}); + const DiffractionExperiment &x = scene.stacked; + const auto w = static_cast(x.GetXPixelsNum()), h = static_cast(x.GetYPixelsNum()); + ImagePreprocessorBufferGPU buffer(x.GetPixelsNum()); + for (size_t i = 0; i < buffer.size(); i++) + buffer[i] = scene.stacked_image[i]; + REQUIRE(cudaMemcpy(buffer.getGPUBuffer(), buffer.getBuffer().data(), x.GetPixelsNum() * sizeof(int32_t), + cudaMemcpyHostToDevice) == cudaSuccess); + REQUIRE(cudaDeviceSynchronize() == cudaSuccess); + + const auto same = [&](const std::vector &cpu, const std::vector &gpu) { + REQUIRE(cpu.size() == gpu.size()); + for (size_t i = 0; i < cpu.size(); i++) { + CHECK(cpu[i].RawCoord().x == gpu[i].RawCoord().x); + CHECK(cpu[i].RawCoord().y == gpu[i].RawCoord().y); + CHECK(cpu[i].Count() == gpu[i].Count()); + CHECK(cpu[i].Panel() == gpu[i].Panel()); + const auto c = cpu[i].RawCoord(); + CHECK(cpu[i].Panel() == x.GetModules().ModuleOfImagePixel(c.x, c.y)); + } + }; + const std::vector no_mask(x.GetPixelsNum(), false); + auto stream = std::make_shared(); + + SpotFindingSettings settings; + settings.signal_to_noise_threshold = 3.0f; + settings.photon_count_threshold = 20; + settings.min_pix_per_spot = 1; + ImageSpotFinderCPU cpu(w, h); + ImageSpotFinderGPU gpu(w, h, stream); + for (ImageSpotFinder *f: {static_cast(&cpu), static_cast(&gpu)}) { + f->SetPanels(x.GetModules()); + f->SetResolutionMask(no_mask); + } + const auto cpu_spots = cpu.Run(buffer, settings); + CHECK(cpu_spots.size() == 4); + same(cpu_spots, gpu.Run(buffer, settings)); + + const AzimuthalIntegrationMapping mapping(x, scene.stacked_mask); + SpotFindingSettings adaptive = settings; + adaptive.adaptive_threshold = true; + adaptive.false_pixels_per_frame = 10.0f; + AdaptiveSpotFinderCPU adaptive_cpu(mapping); + AdaptiveSpotFinderGPU adaptive_gpu(mapping, stream); + for (ImageSpotFinder *f: {static_cast(&adaptive_cpu), static_cast(&adaptive_gpu)}) { + f->SetPanels(x.GetModules()); + f->SetResolutionMask(no_mask); + } + const auto adaptive_spots = adaptive_cpu.Run(buffer, adaptive); + CHECK(adaptive_spots.size() == 4); + same(adaptive_spots, adaptive_gpu.Run(buffer, adaptive)); +} +#endif + +namespace { + +// The raster the stacked panels are assembled to starts at the plane origin and ends at the last panel +// pixel, so it is the assembled image cut to its size: the same picture, in which the raster algorithms +// must find the same answer for both layouts. +template +std::vector CutToRaster(const std::vector &assembled, const PlaneRaster &r) { + std::vector out; + for (int64_t j = 0; j < r.height; j++) + for (int64_t i = 0; i < r.width; i++) + out.push_back(assembled[(r.y0 + j) * PLANE_W + r.x0 + i]); + return out; +} + +PlaneRaster StackedRaster() { + return PlaneRasterOf(StackedPanels(), PW, NPANELS * PH); +} + +// A plane picture: background b, nothing inside the disk of radius r about the beam and along an arm running +// right of it, and a ring of radius 60 about (cx, cy). +std::vector Frame(const DetectorModules &modules, int64_t width, int64_t height, bool stacked, + double r, double cx, double cy) { + const auto value = [&](double px, double py) { + const double dx = px - 110.0, dy = py - 82.0; + if (dx * dx + dy * dy <= r * r || (dx >= 0 && std::fabs(dy) <= 4.0)) + return int32_t{0}; + const double ring = std::hypot(px - cx, py - cy) - 60.0; + return static_cast(4 + std::lround(20.0 * std::exp(-ring * ring / 8.0))); + }; + std::vector out(static_cast(width * height), INT32_MIN); + if (stacked) + modules.ForEachPixel(width, height, [&](int64_t, size_t i, double px, double py) { out[i] = value(px, py); }); + else + StackedPanels().ForEachPixel(PW, NPANELS * PH, [&](int64_t, size_t, double px, double py) { + out[std::lround(py) * width + std::lround(px)] = value(px, py); + }); + return out; +} + +// The assembled detector at the raster's size, gaps masked. +struct RasterScene { + PlaneRaster raster = StackedRaster(); + DiffractionExperiment assembled = Experiment(raster.width, raster.height, {}); + PixelMask assembled_mask; + DiffractionExperiment stacked = Experiment(PW, NPANELS * PH, StackedPanels()); + PixelMask stacked_mask{PW, NPANELS * PH}; + + RasterScene() { + std::vector gaps(raster.width * raster.height, 1u << PixelMask::ModuleGapPixelBit); + StackedPanels().ForEachPixel(PW, NPANELS * PH, [&](int64_t, size_t, double px, double py) { + gaps[std::lround(py) * raster.width + std::lround(px)] = 0; + }); + assembled_mask = PixelMask(gaps); + } +}; + +} // namespace + +// The panels put together on the plane raster are the assembled image, gaps and all; an empty table is a copy. +TEST_CASE("PanelLayout_AssembleToPlane", "[PanelLayout]") { + const Scene scene(INTERIOR_SPOTS); + const PlaneRaster r = StackedRaster(); + CHECK(r.x0 == 0); + CHECK(r.y0 == 0); + CHECK(r.width == 210); + CHECK(r.height == 186); + const std::vector stacked(scene.stacked_image.data(), scene.stacked_image.data() + scene.stacked_image.size()); + const std::vector assembled(scene.assembled_image.data(), + scene.assembled_image.data() + scene.assembled_image.size()); + const auto plane = AssembleToPlane(scene.modules, stacked, PW, NPANELS * PH, r, INT32_MIN); + CHECK(plane == CutToRaster(assembled, r)); + CHECK(GatherFromPlane(scene.modules, plane, r, PW, NPANELS * PH, 0) == stacked); + + const DetectorModules none; + const PlaneRaster image = PlaneRasterOf(none, PLANE_W, PLANE_H); + CHECK(image.width == PLANE_W); + CHECK(image.height == PLANE_H); + CHECK(AssembleToPlane(none, assembled, PLANE_W, PLANE_H, image, 0) == assembled); + CHECK(GatherFromPlane(none, assembled, image, PLANE_W, PLANE_H, 0) == assembled); +} + +// The beam stop - a disk on the beam in the gap cross, reaching onto all four panels, and an arm across the +// gap between panels 0 and 1 - is found on the plane: each stored pixel gets the mask of the assembled pixel +// it is. +TEST_CASE("PanelLayout_ShadowFinder", "[PanelLayout]") { + const RasterScene rs; + ShadowFinder assembled(rs.assembled, rs.assembled_mask); + ShadowFinder stacked(rs.stacked, rs.stacked_mask); + std::vector> frames; + std::vector buffer; + for (int f = 0; f < 12; f++) { + frames.push_back(Frame({}, rs.raster.width, rs.raster.height, false, 25.0, 110.0, 82.0)); + DataMessage a{}; + a.image = CompressedImage(frames.back(), rs.raster.width, rs.raster.height); + assembled.AddImage(a, buffer); + frames.push_back(Frame(StackedPanels(), PW, NPANELS * PH, true, 25.0, 110.0, 82.0)); + DataMessage s{}; + s.image = CompressedImage(frames.back(), PW, NPANELS * PH); + stacked.AddImage(s, buffer); + } + const auto a = assembled.GetMask(); + const auto s = stacked.GetMask(); + REQUIRE(s.size() == static_cast(PW * NPANELS * PH)); + size_t shadowed = 0, differ = 0; + StackedPanels().ForEachPixel(PW, NPANELS * PH, [&](int64_t, size_t i, double px, double py) { + shadowed += s[i] != 0; + differ += s[i] != a[std::lround(py) * rs.raster.width + std::lround(px)]; + }); + CHECK(shadowed > 500); + CHECK(differ == 0); +} + +// The point-symmetry capture scores the plane raster and reports plane coordinates: the same shortlist from +// both layouts, led by the centre of the ring. +TEST_CASE("PanelLayout_BeamCenterFFT", "[PanelLayout]") { + const RasterScene rs; + const auto to_float = [](const std::vector &v) { + std::vector out; + for (int32_t i: v) + out.push_back(i == INT32_MIN ? NAN : static_cast(i)); + return out; + }; + const auto a = FindBeamCenterFFT(rs.assembled, to_float(Frame({}, rs.raster.width, rs.raster.height, false, + 0.0, 104.0, 95.0))); + const auto s = FindBeamCenterFFT(rs.stacked, to_float(Frame(StackedPanels(), PW, NPANELS * PH, true, + 0.0, 104.0, 95.0))); + REQUIRE(!a.point.empty()); + REQUIRE(s.point.size() == a.point.size()); + for (size_t i = 0; i < a.point.size(); i++) { + CHECK(s.point[i].beam_x_pxl == a.point[i].beam_x_pxl); + CHECK(s.point[i].beam_y_pxl == a.point[i].beam_y_pxl); + CHECK(s.point[i].score == a.point[i].score); + } + CHECK(s.margin_point == a.margin_point); +} + +// The diagnostic picture shows the detector as it is, gaps in the gap colour: the same JPEG from both layouts. +TEST_CASE("PanelLayout_DetectorDiagnostic", "[PanelLayout]") { + const RasterScene rs; + const auto to_float = [](const std::vector &v) { + std::vector out; + for (int32_t i: v) + out.push_back(i == INT32_MIN ? NAN : static_cast(i)); + return out; + }; + const auto a = RenderDetectorDiagnostic(rs.assembled, rs.assembled_mask, + to_float(Frame({}, rs.raster.width, rs.raster.height, false, + 25.0, 104.0, 95.0))); + const auto s = RenderDetectorDiagnostic(rs.stacked, rs.stacked_mask, + to_float(Frame(StackedPanels(), PW, NPANELS * PH, true, + 25.0, 104.0, 95.0))); + CHECK(s == a); +}