From dcad4104cea8560a48575eb8093b3b8e00839b51 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 21 Oct 2026 10:26:30 +0200 Subject: [PATCH] Dense per-pixel maps through the panel table (panel layout stage 1) Every dense map that needs geometry now evaluates each image pixel at its place on the detector plane (DetectorModules), so it is right for a stacked panel layout and picks up refined module placements: - AzimuthalIntegrationMapping: bins, resolution, solid-angle/polarisation corrections built per pixel through the table; a pixel on no panel is left as a masked one. The geometry cache key now includes the panel table, so a mapping built after the stills module placement is rebuilt rather than served from the cache. ResolutionMaskBits follows (it is a lookup). - ROIMap / ExportROIMap: azimuthal ROIs on the plane (box and circle ROIs stay in image pixels). - ROIIntegration CPU+GPU: the weighted sums are taken per ROI and panel in image pixels and placed on the plane at the end (exact; the affine commutes with a weighted mean). - HotPixelFinder: ring/sector keys and the zeta of the persistence rule on the plane. - Background beam-centre walk CPU+GPU: pixels binned at their plane position; host and device share BackgroundPlanePoint so they stay bit-identical. New common/PanelRows.h: ForEachPixelInRows, ForEachPixel cut to a band of rows in index order, for the row-parallel builders. An empty table takes the old arithmetic in the old order: rotation sets and every single-panel input are unchanged. [PanelLayout] AzintMapping and ResolutionMask enabled (+ azimuthal ROI map and ROI centroid on the plane); new BeamCenterFromBackground_BinsOnThePlane and ROIIntegrationGPU_MatchesCPU_Panels. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_015CeEGgpKdfkCTJb3YNMJcu --- common/AzimuthalIntegrationMapping.cpp | 32 +++++-- common/AzimuthalIntegrationMapping.h | 8 +- common/CMakeLists.txt | 1 + common/DiffractionExperiment.cpp | 2 +- common/PanelRows.h | 39 ++++++++ common/ROIMap.cpp | 18 ++-- common/ROIMap.h | 6 +- .../geom_refinement/BackgroundBand.h | 20 +++- .../BeamCenterBackgroundGPU.cu | 33 ++++++- .../geom_refinement/BeamCenterBackgroundGPU.h | 3 +- .../BeamCenterFromBackground.cpp | 48 ++++++---- image_analysis/hot_pixels/HotPixels.cpp | 55 +++++------ image_analysis/hot_pixels/HotPixels.h | 1 + image_analysis/roi/ROIIntegration.cpp | 33 ++++++- image_analysis/roi/ROIIntegration.h | 16 +++- image_analysis/roi/ROIIntegrationCPU.h | 23 +++-- image_analysis/roi/ROIIntegrationGPU.cu | 91 ++++++++++++++----- image_analysis/roi/ROIIntegrationGPU.h | 11 ++- tests/BeamCenterFromBackgroundTest.cpp | 39 ++++++++ tests/PanelLayoutTest.cpp | 34 +++++-- tests/ROIIntegrationGPUTest.cpp | 46 +++++++--- tests/ROIMapTest.cpp | 2 +- 22 files changed, 427 insertions(+), 134 deletions(-) create mode 100644 common/PanelRows.h diff --git a/common/AzimuthalIntegrationMapping.cpp b/common/AzimuthalIntegrationMapping.cpp index f33ba2ca6..7027a5d24 100644 --- a/common/AzimuthalIntegrationMapping.cpp +++ b/common/AzimuthalIntegrationMapping.cpp @@ -14,6 +14,7 @@ #include "RawToConvertedGeometry.h" #include "TableChecksum.h" #include "ParallelFor.h" +#include "PanelRows.h" AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const DiffractionExperiment &experiment, const PixelMask &mask, @@ -36,7 +37,8 @@ AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const DiffractionExperi // transformed, raw when it is not. They have to match pixel_to_bin, which is sized per mode // below - the adaptive spot finders walk the image with these and index pixel_to_bin with it. width(experiment.GetXPixelsNum()), - height(experiment.GetYPixelsNum()) { + height(experiment.GetYPixelsNum()), + modules(experiment.GetModules()) { if (width <= 0) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, @@ -74,10 +76,13 @@ AzimuthalIntegrationMapping::AzimuthalIntegrationMapping(const DiffractionExperi void AzimuthalIntegrationMapping::SetupConvGeomRows(const DiffractionGeometry &geom, const std::vector *mask, size_t row0, size_t row_end) { - for (size_t row = row0; row < row_end && row < height; row++) { - for (size_t col = 0; col < width; col++) - SetupPixel(geom, mask, row * width + col, col, row); - } + // Each pixel at its place on the plane. A pixel on no panel keeps what a masked one has. + ForEachPixelInRows(modules, static_cast(width), static_cast(row0), + static_cast(std::min(row_end, height)), + [&](int64_t, size_t pxl, double px, double py) { + SetupPixel(geom, mask, static_cast(pxl), static_cast(px), + static_cast(py)); + }); } void AzimuthalIntegrationMapping::SetupConvGeom(const DiffractionGeometry &geom, const std::vector *mask) { @@ -127,6 +132,15 @@ void AzimuthalIntegrationMapping::SetupConvGeomCached(const DiffractionGeometry static_cast(settings.IsPolarizationCorrection()), static_cast(polarization_factor.has_value()), static_cast(width), static_cast(height)}); + // The panel table: the same geometry places the pixels elsewhere when a panel has moved. + for (const auto &m: modules.Get()) { + for (const int64_t i: {m.x0, m.y0, m.width, m.height}) + key.push_back(static_cast(i)); + for (const double d: {m.ox, m.oy, m.fx, m.fy, m.sx, m.sy}) { + const auto bits = std::bit_cast(d); + key.insert(key.end(), {static_cast(bits), static_cast(bits >> 32)}); + } + } std::lock_guard lock(cache.m); auto &entries = cache.entries; @@ -169,7 +183,7 @@ void AzimuthalIntegrationMapping::SetupRawGeom(const DiffractionExperiment &expe for (int m = 0; m < experiment.GetModulesNum(); m++) { for (int pxl = 0; pxl < RAW_MODULE_SIZE; pxl++) { auto [x,y] = RawToConvertedCoordinate(experiment, m, pxl); - SetupPixel(geom, &mask, m * RAW_MODULE_SIZE + pxl, x, y); + SetupPixel(geom, &mask, m * RAW_MODULE_SIZE + pxl, static_cast(x), static_cast(y)); } } } else { @@ -185,7 +199,7 @@ void AzimuthalIntegrationMapping::SetupRawGeom(const DiffractionExperiment &expe for (size_t m = module_begin; m < module_end; ++m) { for (int pxl = 0; pxl < RAW_MODULE_SIZE; ++pxl) { auto [x, y] = RawToConvertedCoordinate(experiment, m, pxl); - SetupPixel(geom, &mask, m * RAW_MODULE_SIZE + pxl, x, y); + SetupPixel(geom, &mask, m * RAW_MODULE_SIZE + pxl, static_cast(x), static_cast(y)); } } })); @@ -198,12 +212,10 @@ void AzimuthalIntegrationMapping::SetupRawGeom(const DiffractionExperiment &expe void AzimuthalIntegrationMapping::SetupPixel(const DiffractionGeometry &geom, const std::vector *mask, - uint32_t pxl, uint32_t col, uint32_t row) { + uint32_t pxl, float x, float y) { if (mask && (*mask)[pxl] != 0) return; - auto x = static_cast(col); - auto y = static_cast(row); float d = geom.PxlToRes(x, y); float phi_rad = geom.Phi_rad(x, y); diff --git a/common/AzimuthalIntegrationMapping.h b/common/AzimuthalIntegrationMapping.h index 23cd4f894..40d2a2a62 100644 --- a/common/AzimuthalIntegrationMapping.h +++ b/common/AzimuthalIntegrationMapping.h @@ -22,7 +22,7 @@ class AzimuthalIntegrationGeometryCache { friend class AzimuthalIntegrationMapping; static constexpr size_t MAX_ENTRIES = 6; struct Entry { - std::vector key; // bit patterns of everything SetupPixel reads + std::vector key; // bit patterns of everything SetupPixel reads, and the panel table std::vector pixel_to_bin; std::vector pixel_resolution; std::vector corrections; @@ -36,6 +36,9 @@ protected: const AzimuthalIntegrationSettings settings; const float wavelength; const size_t width, height; + // Where the image's pixels are on the plane the geometry describes; empty when they are where the + // image has them. + const DetectorModules modules; std::vector bin_to_q; std::vector bin_to_2theta; @@ -76,8 +79,9 @@ protected: void SetupConvGeomCached(const DiffractionGeometry &geom, const std::vector &mask, AzimuthalIntegrationGeometryCache &cache); + // Pixel pxl of the image, at (x, y) on the plane. void SetupPixel(const DiffractionGeometry &geom, const std::vector *mask, - uint32_t pxl, uint32_t col, uint32_t row); + uint32_t pxl, float x, float y); AzimuthalIntegrationMapping(const DiffractionExperiment& experiment, const PixelMask& mask, diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index dcd7f7ddd..63d53e882 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -132,6 +132,7 @@ ADD_LIBRARY(JFJochCommon STATIC DetectorGeometryModular.h DetectorModules.cpp DetectorModules.h + PanelRows.h ModuleShift.h DetectorGeometryFixed.cpp DetectorGeometryFixed.h diff --git a/common/DiffractionExperiment.cpp b/common/DiffractionExperiment.cpp index db1ce112b..3c0e059d7 100644 --- a/common/DiffractionExperiment.cpp +++ b/common/DiffractionExperiment.cpp @@ -1160,7 +1160,7 @@ const ROIMap &DiffractionExperiment::ROI() const { } std::vector DiffractionExperiment::ExportROIMap() const { - return roi_mask.GetROIMap(GetDiffractionGeometry(), GetXPixelsNumConv(), GetYPixelsNumConv()); + return roi_mask.GetROIMap(GetDiffractionGeometry(), GetModules(), GetXPixelsNumConv(), GetYPixelsNumConv()); } DiffractionExperiment &DiffractionExperiment::ImagesPerFile(const std::optional &input) { diff --git a/common/PanelRows.h b/common/PanelRows.h new file mode 100644 index 000000000..54b1649bd --- /dev/null +++ b/common/PanelRows.h @@ -0,0 +1,39 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include + +#include "DetectorModules.h" + +// DetectorModules::ForEachPixel cut to the image rows [row0, row1), so that a dense map can be built a band +// of rows per thread. Calls f(panel, index, px, py) for every pixel of those rows that is on a panel, in +// index order: row by row, and within a row panel by panel from left to right. Pixels on no panel are not +// visited. Empty table: every pixel of the rows, panel -1, plane == image. +template +void ForEachPixelInRows(const DetectorModules &modules, int64_t width, int64_t row0, int64_t row1, F &&f) { + if (modules.empty()) { + for (int64_t y = row0; y < row1; y++) + for (int64_t x = 0; x < width; x++) + f(int64_t(-1), static_cast(y * width + x), static_cast(x), static_cast(y)); + return; + } + const auto &m = modules.Get(); + std::vector on_row; + for (int64_t y = row0; y < row1; y++) { + on_row.clear(); + for (size_t p = 0; p < m.size(); p++) + if (y >= m[p].y0 && y < m[p].y0 + m[p].height) + on_row.push_back(static_cast(p)); + std::sort(on_row.begin(), on_row.end(), [&](int64_t a, int64_t b) { return m[a].x0 < m[b].x0; }); + for (const int64_t p: on_row) + for (int64_t x = m[p].x0; x < m[p].x0 + m[p].width; x++) { + const auto [px, py] = m[p].ImageToPlane(static_cast(x), static_cast(y)); + f(p, static_cast(y * width + x), px, py); + } + } +} diff --git a/common/ROIMap.cpp b/common/ROIMap.cpp index adb97ab5e..e578602ea 100644 --- a/common/ROIMap.cpp +++ b/common/ROIMap.cpp @@ -4,20 +4,22 @@ #include "JFJochMath.h" #include "ROIMap.h" #include "JFJochException.h" +#include "PanelRows.h" -std::vector ROIMap::GetROIMap(const DiffractionGeometry &geom, size_t xpixel, size_t ypixel) const { +std::vector ROIMap::GetROIMap(const DiffractionGeometry &geom, const DetectorModules &modules, + size_t xpixel, size_t ypixel) const { std::vector map(xpixel * ypixel, 0); std::vector resolution_map(xpixel * ypixel, 0); std::vector phi_map; if (!rois.azimuthal.empty()) { phi_map.resize(xpixel * ypixel, 0); - for (int y = 0; y < ypixel; y++) { - for (int x = 0; x < xpixel; x++) { - resolution_map[y * xpixel + x] = geom.PxlToRes(static_cast(x), static_cast(y)); - phi_map[y * xpixel + x] = geom.Phi_rad(static_cast(x), static_cast(y)) - * 180.0f / static_cast(PI); - } - } + // A pixel on no panel keeps resolution and angle 0. + ForEachPixelInRows(modules, static_cast(xpixel), 0, static_cast(ypixel), + [&](int64_t, size_t i, double px, double py) { + const auto x = static_cast(px), y = static_cast(py); + resolution_map[i] = geom.PxlToRes(x, y); + phi_map[i] = geom.Phi_rad(x, y) * 180.0f / static_cast(PI); + }); } uint16_t roi_id = 0; diff --git a/common/ROIMap.h b/common/ROIMap.h index 40be957dd..a3c459881 100644 --- a/common/ROIMap.h +++ b/common/ROIMap.h @@ -11,6 +11,7 @@ #include "DetectorSetup.h" #include "ROIDefinition.h" #include "DiffractionGeometry.h" +#include "DetectorModules.h" class ROIMap { ROIDefinition rois; @@ -18,7 +19,10 @@ class ROIMap { static constexpr const size_t roi_limit = 16; public: void SetROI(const ROIDefinition& input); - [[nodiscard]] std::vector GetROIMap(const DiffractionGeometry &geom, size_t width, size_t height) const; + // Box and circle ROIs are drawn in image pixels; azimuthal ROIs at each pixel's place on the plane + // (`modules`, empty when that is where the image has it). + [[nodiscard]] std::vector GetROIMap(const DiffractionGeometry &geom, const DetectorModules &modules, + size_t width, size_t height) const; [[nodiscard]] std::map GetROINameMap() const; [[nodiscard]] const ROIDefinition& GetROIDefinition() const; [[nodiscard]] bool empty() const; diff --git a/image_analysis/geom_refinement/BackgroundBand.h b/image_analysis/geom_refinement/BackgroundBand.h index 394ab3e84..a7124bc34 100644 --- a/image_analysis/geom_refinement/BackgroundBand.h +++ b/image_analysis/geom_refinement/BackgroundBand.h @@ -40,6 +40,20 @@ struct BackgroundBand { double tan_lo, tan_hi; // the band in tan(2theta), widened for the quick rejection }; +// A panel of the image and its place on the plane (DetectorModule), in a form the device can hold. +struct BackgroundPanel { + int x0, y0, width, height; + double ox, oy, fx, fy, sx, sy; +}; + +// Plane position of image pixel (x, y) of panel p - DetectorModule::ImageToPlane, written here so that the +// host and the device evaluate it the same way. +BACKGROUND_BAND_HD void BackgroundPlanePoint(const BackgroundPanel &p, int x, int y, float &px, float &py) { + const double u = x - p.x0, v = y - p.y0; + px = static_cast(p.ox + u * p.fx + v * p.sx); + py = static_cast(p.oy + u * p.fy + v * p.sy); +} + // atan2(y, x) from IEEE operations alone - add, multiply, divide, square root, all correctly rounded on // the host and on the device - so that, compiled without contraction (see CMakeLists.txt), the host // and the device return the same bits. The library atan2f does not: glibc's and CUDA's differ in the @@ -63,9 +77,9 @@ BACKGROUND_BAND_HD double BackgroundAtan2(double y, double x) { return y < 0 ? -a : a; } -// The cell of pixel (x, y) about (beam_x, beam_y), or -1 when it is outside the band; for a pixel in -// it, also the two components of the derivative of its 2theta with respect to the centre. -BACKGROUND_BAND_HD int BackgroundBandCell(const BackgroundBand &b, int x, int y, float beam_x, float beam_y, +// The cell of the pixel at (x, y) on the plane about (beam_x, beam_y), or -1 when it is outside the band; +// for a pixel in it, also the two components of the derivative of its 2theta with respect to the centre. +BACKGROUND_BAND_HD int BackgroundBandCell(const BackgroundBand &b, float x, float y, float beam_x, float beam_y, float &jac_x, float &jac_y) { const float *rot = b.rot; const float u = (x - beam_x) * b.pixel_size; diff --git a/image_analysis/geom_refinement/BeamCenterBackgroundGPU.cu b/image_analysis/geom_refinement/BeamCenterBackgroundGPU.cu index 6dfb87329..9901175ac 100644 --- a/image_analysis/geom_refinement/BeamCenterBackgroundGPU.cu +++ b/image_analysis/geom_refinement/BeamCenterBackgroundGPU.cu @@ -19,7 +19,9 @@ void check(cudaError_t err, const char *what) { // Every pixel's cell (BackgroundBand::CELLS where it is outside the band or unusable), its index, // and its derivatives. +// A pixel on no panel is outside the band. __global__ void bin_kernel(BackgroundBand band, int width, size_t npixels, float beam_x, float beam_y, + const BackgroundPanel *__restrict__ panels, int npanels, const char *__restrict__ usable, int32_t *__restrict__ key, int32_t *__restrict__ index, float *__restrict__ jac_x, float *__restrict__ jac_y) { @@ -29,10 +31,21 @@ __global__ void bin_kernel(BackgroundBand band, int width, size_t npixels, float int cell = BackgroundBand::CELLS; float jx = 0.0f, jy = 0.0f; if (usable[i]) { - const int c = BackgroundBandCell(band, static_cast(i % width), static_cast(i / width), - beam_x, beam_y, jx, jy); - if (c >= 0) - cell = c; + const int x = static_cast(i % width), y = static_cast(i / width); + float px = static_cast(x), py = static_cast(y); + bool on_panel = npanels == 0; + for (int p = 0; p < npanels && !on_panel; p++) { + const BackgroundPanel &m = panels[p]; + if (x >= m.x0 && x < m.x0 + m.width && y >= m.y0 && y < m.y0 + m.height) { + BackgroundPlanePoint(m, x, y, px, py); + on_panel = true; + } + } + if (on_panel) { + const int c = BackgroundBandCell(band, px, py, beam_x, beam_y, jx, jy); + if (c >= 0) + cell = c; + } } key[i] = cell; index[i] = static_cast(i); @@ -112,6 +125,8 @@ struct BeamCenterBackgroundGPU::Impl { CudaDevicePtr usable; CudaDevicePtr mean; CudaDevicePtr row_block; + CudaDevicePtr panels; + int npanels = 0; CudaDevicePtr key, sorted_key, index, sorted_index; CudaDevicePtr jac_x, jac_y; CudaDevicePtr offset; @@ -159,8 +174,15 @@ struct BeamCenterBackgroundGPU::Impl { }; BeamCenterBackgroundGPU::BeamCenterBackgroundGPU(int width, int height, const std::vector &block_row, + const std::vector &panels, const char *usable, const float *mean) : impl(std::make_unique(width, height)) { + impl->npanels = static_cast(panels.size()); + if (!panels.empty()) { + impl->panels = CudaDevicePtr(panels.size()); + check(cudaMemcpyAsync(impl->panels.get(), panels.data(), panels.size() * sizeof(BackgroundPanel), + cudaMemcpyHostToDevice, impl->stream), "upload"); + } std::vector row_block(height); for (size_t b = 0; b + 1 < block_row.size(); b++) for (int y = block_row[b]; y < block_row[b + 1]; y++) @@ -181,7 +203,8 @@ void BeamCenterBackgroundGPU::Bin(const BackgroundBand &band, float beam_x, floa std::vector &count) { Impl &d = *impl; const auto blocks = static_cast((d.npixels + THREADS - 1) / THREADS); - bin_kernel<<>>(band, d.width, d.npixels, beam_x, beam_y, d.usable, + bin_kernel<<>>(band, d.width, d.npixels, beam_x, beam_y, d.panels.get(), + d.npanels, d.usable, d.key, d.index, d.jac_x, d.jac_y); check(cudaGetLastError(), "bin"); // A radix sort is stable, so each cell's pixels come out in pixel order. diff --git a/image_analysis/geom_refinement/BeamCenterBackgroundGPU.h b/image_analysis/geom_refinement/BeamCenterBackgroundGPU.h index 72d47b266..23a0f2d9d 100644 --- a/image_analysis/geom_refinement/BeamCenterBackgroundGPU.h +++ b/image_analysis/geom_refinement/BeamCenterBackgroundGPU.h @@ -23,8 +23,9 @@ class BeamCenterBackgroundGPU { std::unique_ptr impl; public: // block_row: the first row of each of the host's row blocks, and one past the last row at the end. + // panels: where the image's panels are on the plane; empty when the image is the plane. BeamCenterBackgroundGPU(int width, int height, const std::vector &block_row, - const char *usable, const float *mean); + const std::vector &panels, const char *usable, const float *mean); ~BeamCenterBackgroundGPU(); // Bin the band about (beam_x, beam_y) and sum each cell: the values, their squares, the two diff --git a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp index 3f7512114..8fa39f795 100644 --- a/image_analysis/geom_refinement/BeamCenterFromBackground.cpp +++ b/image_analysis/geom_refinement/BeamCenterFromBackground.cpp @@ -11,6 +11,7 @@ #include "../../common/CompressedImage.h" #include "../../common/JFJochMath.h" #include "../../common/ParallelFor.h" +#include "../../common/PanelRows.h" #ifdef JFJOCH_USE_CUDA #include "../../common/CUDAWrapper.h" #include "BeamCenterBackgroundGPU.h" @@ -181,13 +182,20 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe band.tan_lo = tan_lo; band.tan_hi = tan_hi; + // The pixels are binned at their place on the plane. + const DetectorModules &modules = experiment.GetModules(); + std::vector panels; + for (const auto &m: modules.Get()) + panels.push_back({static_cast(m.x0), static_cast(m.y0), static_cast(m.width), + static_cast(m.height), m.ox, m.oy, m.fx, m.fy, m.sx, m.sy}); + #ifndef JFJOCH_USE_CUDA (void) allow_device; #else // With a GPU the two passes over the pixels run there, and only the cells come back. std::unique_ptr gpu; if (allow_device && get_gpu_count() > 0) - gpu = std::make_unique(W, H, block_row, usable.data(), mean.data()); + gpu = std::make_unique(W, H, block_row, panels, usable.data(), mean.data()); #endif float step_x = 0.0f, step_y = 0.0f, sigma_x = 0.0f, sigma_y = 0.0f; @@ -209,25 +217,25 @@ FindBeamCenterFromBackground(const DiffractionExperiment &experiment, const Pixe int32_t *cells = band_cell.data() + static_cast(block_row[b]) * W; float *values = band_value.data() + static_cast(block_row[b]) * W; size_t n_band = 0; - for (int y = block_row[b]; y < block_row[b + 1]; y++) { - for (int x = 0; x < W; x++) { - const size_t i = static_cast(y) * W + x; - if (!usable[i]) - continue; - float jac_x, jac_y; - const int cell = BackgroundBandCell(band, x, y, beam_x, beam_y, jac_x, jac_y); - if (cell < 0) - continue; - cells[n_band] = cell; - values[n_band] = mean[i]; - n_band++; - b_count[cell]++; - b_sum[cell] += mean[i]; - b_sum_sq[cell] += static_cast(mean[i]) * mean[i]; - b_jx[cell] += jac_x; - b_jy[cell] += jac_y; - } - } + ForEachPixelInRows(modules, W, block_row[b], block_row[b + 1], [&](int64_t p, size_t i, double px, double py) { + if (!usable[i]) + return; + auto x = static_cast(px), y = static_cast(py); + if (p >= 0) // as the device has it + BackgroundPlanePoint(panels[p], static_cast(i % W), static_cast(i / W), x, y); + float jac_x, jac_y; + const int cell = BackgroundBandCell(band, x, y, beam_x, beam_y, jac_x, jac_y); + if (cell < 0) + return; + cells[n_band] = cell; + values[n_band] = mean[i]; + n_band++; + b_count[cell]++; + b_sum[cell] += mean[i]; + b_sum_sq[cell] += static_cast(mean[i]) * mean[i]; + b_jx[cell] += jac_x; + b_jy[cell] += jac_y; + }); band_pixels[b] = n_band; }); fold(true); diff --git a/image_analysis/hot_pixels/HotPixels.cpp b/image_analysis/hot_pixels/HotPixels.cpp index 139e3caca..29937caec 100644 --- a/image_analysis/hot_pixels/HotPixels.cpp +++ b/image_analysis/hot_pixels/HotPixels.cpp @@ -12,6 +12,7 @@ #include "../../common/CUDAWrapper.h" #include "../../common/JFJochMath.h" #include "../../common/ParallelFor.h" +#include "../../common/PanelRows.h" #include @@ -75,6 +76,7 @@ int binomial_bound(int n, double q, double npix, double rate) { HotPixelFinder::HotPixelFinder(const DiffractionExperiment &experiment, const PixelMask &mask, size_t nthreads) : geometry(experiment.GetDiffractionGeometry()), + modules(experiment.GetModules()), axis(experiment.GetGoniometer() ? experiment.GetGoniometer()->GetAxis().Normalize() : Coord(0, 0, 0)), width(experiment.GetXPixelsNumConv()), height(experiment.GetYPixelsNumConv()), key(std::make_unique_for_overwrite(width * height)), @@ -87,19 +89,20 @@ HotPixelFinder::HotPixelFinder(const DiffractionExperiment &experiment, const Pi // Rings of equal 2theta, RING_WIDTH_PX wide where the beam meets the detector. const float ring_rad = RING_WIDTH_PX * geometry.GetPixelSize_mm() / geometry.GetDetectorDistance_mm(); ParallelChunks(static_cast(height), nthreads, [&](int y0, int y1) { - for (size_t y = y0; y < static_cast(y1); y++) - for (size_t x = 0; x < width; x++) { - const size_t i = y * width + x; - n_lit[i] = n_error[i] = n_error_ring_ok[i] = 0; - sum_value[i] = error_level_sum[i] = 0; - key[i] = -1; - if (m[i] != 0) continue; - const float fx = static_cast(x), fy = static_cast(y); - const int ring = static_cast(geometry.TwoTheta_rad(fx, fy) / ring_rad); - const int sector = std::min(SECTORS - 1, static_cast(geometry.Phi_rad(fx, fy) - / static_cast(2.0 * PI) * SECTORS)); - key[i] = ring * SECTORS + sector; - } + for (size_t i = y0 * width; i < y1 * width; i++) { + n_lit[i] = n_error[i] = n_error_ring_ok[i] = 0; + sum_value[i] = error_level_sum[i] = 0; + key[i] = -1; + } + // Each pixel keyed at its place on the plane; one on no panel keeps -1. + ForEachPixelInRows(modules, static_cast(width), y0, y1, [&](int64_t, size_t i, double px, double py) { + if (m[i] != 0) return; + const float fx = static_cast(px), fy = static_cast(py); + const int ring = static_cast(geometry.TwoTheta_rad(fx, fy) / ring_rad); + const int sector = std::min(SECTORS - 1, static_cast(geometry.Phi_rad(fx, fy) + / static_cast(2.0 * PI) * SECTORS)); + key[i] = ring * SECTORS + sector; + }); }); for (size_t i = 0; i < width * height; i++) if (key[i] >= 0) { @@ -377,20 +380,18 @@ HotPixelFinder::Result HotPixelFinder::GetMask(double oscillation_deg, double sp #endif std::vector persistent(width * height, 0); ParallelChunks(static_cast(height), nthreads, [&](int y0, int y1) { - for (size_t y = y0; y < static_cast(y1); y++) - for (size_t x = 0; x < width; x++) { - const size_t i = y * width + x; - if (key[i] < 0 || n_valid(i) < min_valid || n_lit[i] == 0 || !(spacing_deg > 0.0)) continue; - // One reflection's stay on this pixel, in sampled frames: 1 / |zeta| turns the rocking - // width into rotation, and zeta = |axis . (s1 x s0)| with s0 along z. - const Coord s1 = geometry.LabCoord(static_cast(x), static_cast(y)).Normalize(); - const double zeta = std::max(1e-4, static_cast(std::fabs(axis * (s1 % Coord(0, 0, 1))))); - const double stay = (oscillation_deg + ROCKING_WIDTH_DEG) / zeta; - const double k_bragg = 1.0 + std::ceil(stay / spacing_deg); - const double k = std::min(std::max(k_bragg, k_chance[key[i] / SECTORS]), - n_valid(i)); - persistent[i] = n_lit[i] >= k; - } + ForEachPixelInRows(modules, static_cast(width), y0, y1, [&](int64_t, size_t i, double px, double py) { + if (key[i] < 0 || n_valid(i) < min_valid || n_lit[i] == 0 || !(spacing_deg > 0.0)) return; + // One reflection's stay on this pixel, in sampled frames: 1 / |zeta| turns the rocking + // width into rotation, and zeta = |axis . (s1 x s0)| with s0 along z. + const Coord s1 = geometry.LabCoord(static_cast(px), static_cast(py)).Normalize(); + const double zeta = std::max(1e-4, static_cast(std::fabs(axis * (s1 % Coord(0, 0, 1))))); + const double stay = (oscillation_deg + ROCKING_WIDTH_DEG) / zeta; + const double k_bragg = 1.0 + std::ceil(stay / spacing_deg); + const double k = std::min(std::max(k_bragg, k_chance[key[i] / SECTORS]), + n_valid(i)); + persistent[i] = n_lit[i] >= k; + }); }); // Masked: a persistent pixel whose mean is STRONG_RATIO times its ring's and whose excess is above diff --git a/image_analysis/hot_pixels/HotPixels.h b/image_analysis/hot_pixels/HotPixels.h index 489d708f9..e364871e0 100644 --- a/image_analysis/hot_pixels/HotPixels.h +++ b/image_analysis/hot_pixels/HotPixels.h @@ -115,6 +115,7 @@ public: private: const DiffractionGeometry geometry; + const DetectorModules modules; // where the pixels are on the geometry's plane; empty = as imaged const Coord axis; // rotation axis; zero where there is none const size_t width, height; int nrings = 0; diff --git a/image_analysis/roi/ROIIntegration.cpp b/image_analysis/roi/ROIIntegration.cpp index 0cef8cca9..addaac25e 100644 --- a/image_analysis/roi/ROIIntegration.cpp +++ b/image_analysis/roi/ROIIntegration.cpp @@ -4,6 +4,9 @@ #include "ROIIntegration.h" #include "../../common/DiffractionExperiment.h" +#include +#include + ROIIntegration::ROIIntegration(const DiffractionExperiment &experiment) : roi_map(experiment.ExportROIMap()), width(experiment.GetXPixelsNumConv()), @@ -15,11 +18,39 @@ ROIIntegration::ROIIntegration(const DiffractionExperiment &experiment) roi_pixels(roi_count), roi_x_weighted(roi_count), roi_y_weighted(roi_count), - roi_max(roi_count) { + roi_max(roi_count), + modules(experiment.GetModules()), + panel_slots(std::max(modules.size(), 1)), + roi_panel_v(roi_count * panel_slots), + roi_panel_x(roi_count * panel_slots), + roi_panel_y(roi_count * panel_slots) { for (const auto &[name, id] : experiment.ROI().GetROINameMap()) roi_name[id] = name; } +void ROIIntegration::PlaceWeighted() { + for (uint16_t r = 0; r < roi_count; r++) { + if (modules.empty()) { + roi_x_weighted[r] = roi_panel_x[r]; + roi_y_weighted[r] = roi_panel_y[r]; + continue; + } + // Sum over panels of v * (o + (x - x0) f + (y - y0) s). + double x = 0, y = 0; + for (size_t p = 0; p < modules.size(); p++) { + const auto &m = modules[p]; + const size_t k = r * panel_slots + p; + const auto v = static_cast(roi_panel_v[k]); + const auto u = static_cast(roi_panel_x[k] - m.x0 * roi_panel_v[k]); + const auto w = static_cast(roi_panel_y[k] - m.y0 * roi_panel_v[k]); + x += m.ox * v + m.fx * u + m.sx * w; + y += m.oy * v + m.fy * u + m.sy * w; + } + roi_x_weighted[r] = std::llround(x); + roi_y_weighted[r] = std::llround(y); + } +} + void ROIIntegration::Export(std::map &out) const { for (uint16_t r = 0; r < roi_count; r++) out[roi_name[r]] = ROIMessage{ diff --git a/image_analysis/roi/ROIIntegration.h b/image_analysis/roi/ROIIntegration.h index 2d0ddd26e..8939f9521 100644 --- a/image_analysis/roi/ROIIntegration.h +++ b/image_analysis/roi/ROIIntegration.h @@ -8,6 +8,7 @@ #include #include +#include "../../common/DetectorModules.h" #include "../../common/JFJochMessages.h" #include "../image_preprocessing/ImagePreprocessorBuffer.h" @@ -17,7 +18,9 @@ class DiffractionExperiment; // good-pixel count, maximum and intensity-weighted centre of mass. Each pixel // carries a 16-bit mask selecting which of up to 16 ROIs it belongs to, so a // single pixel can contribute to none, one, or several ROIs at once. This is -// the software counterpart of the FPGA roi_calc kernel. +// the software counterpart of the FPGA roi_calc kernel. The centre of mass is on the detector plane: the +// weighted sums are taken per panel in image pixels and placed on the plane at the end - the panel's +// affine commutes with a weighted mean, so that is exact. class ROIIntegration { protected: std::vector roi_map; // per-pixel ROI bitmask, assembled (converted) geometry @@ -34,6 +37,17 @@ protected: std::vector roi_y_weighted; std::vector roi_max; + // Where the image's panels are on the plane (empty: the image is the plane), and per ROI and panel - + // index roi * panel_slots + panel, slot 0 for the empty table - the sums of value, value * x and + // value * y in image pixels. + DetectorModules modules; + size_t panel_slots; + std::vector roi_panel_v; + std::vector roi_panel_x; + std::vector roi_panel_y; + + // roi_x_weighted / roi_y_weighted on the plane, from the per-panel sums. + void PlaceWeighted(); void Export(std::map &out) const; public: explicit ROIIntegration(const DiffractionExperiment &experiment); diff --git a/image_analysis/roi/ROIIntegrationCPU.h b/image_analysis/roi/ROIIntegrationCPU.h index 0e5fb4eb6..b556a9874 100644 --- a/image_analysis/roi/ROIIntegrationCPU.h +++ b/image_analysis/roi/ROIIntegrationCPU.h @@ -3,6 +3,7 @@ #pragma once +#include #include #include @@ -25,23 +26,26 @@ public: roi_sum[r] = 0; roi_sum2[r] = 0; roi_pixels[r] = 0; - roi_x_weighted[r] = 0; - roi_y_weighted[r] = 0; roi_max[r] = INT64_MIN; } + std::fill(roi_panel_v.begin(), roi_panel_v.end(), 0); + std::fill(roi_panel_x.begin(), roi_panel_x.end(), 0); + std::fill(roi_panel_y.begin(), roi_panel_y.end(), 0); using pixel_t = std::remove_cv_t>; - for (size_t i = 0; i < npixel; i++) { + // Pixels on no panel are not visited. + modules.ForEachPixel(static_cast(width), static_cast(npixel / width), + [&](int64_t panel, size_t i, double, double) { const uint16_t mask = roi_map[i]; if (mask == 0) - continue; + return; const pixel_t v = image[i]; // masked/bad pixels (signed types only) are excluded entirely if constexpr (std::is_signed_v) { if (v == std::numeric_limits::min()) - continue; + return; } // saturated pixels still count towards the max, but not the sums const bool saturated = (v == std::numeric_limits::max()); @@ -49,6 +53,7 @@ public: const int64_t val = static_cast(v); const int64_t x = static_cast(i % width); const int64_t y = static_cast(i / width); + const size_t slot = panel < 0 ? 0 : static_cast(panel); for (uint16_t r = 0; r < roi_count; r++) { if (!(mask & (1u << r))) @@ -57,13 +62,15 @@ public: roi_sum[r] += val; roi_sum2[r] += static_cast(val * val); roi_pixels[r] += 1; - roi_x_weighted[r] += val * x; - roi_y_weighted[r] += val * y; + roi_panel_v[r * panel_slots + slot] += val; + roi_panel_x[r * panel_slots + slot] += val * x; + roi_panel_y[r * panel_slots + slot] += val * y; } if (val > roi_max[r]) roi_max[r] = val; } - } + }); + PlaceWeighted(); Export(out); } diff --git a/image_analysis/roi/ROIIntegrationGPU.cu b/image_analysis/roi/ROIIntegrationGPU.cu index d652c9297..4288d533c 100644 --- a/image_analysis/roi/ROIIntegrationGPU.cu +++ b/image_analysis/roi/ROIIntegrationGPU.cu @@ -13,6 +13,8 @@ inline void cuda_err(cudaError_t val) { // One pixel carries a 16-bit mask, so it can feed any subset of the ROIs. // Each block reduces into shared memory first to keep global atomics low. +// The weighted sums go per ROI and panel (roi * slots + panel; slots = 1 and panel 0 without a panel +// table); panel_box holds x0, y0, width, height of each panel, and a pixel on none is left out. __global__ void gpu_roi( const uint16_t *__restrict__ roi_map, @@ -20,28 +22,36 @@ void gpu_roi( size_t num_pixels, size_t width, int roi_count, + const int *__restrict__ panel_box, + int npanels, + int slots, unsigned long long *__restrict__ roi_sum, unsigned long long *__restrict__ roi_sum2, unsigned long long *__restrict__ roi_pixels, - unsigned long long *__restrict__ roi_x_weighted, - unsigned long long *__restrict__ roi_y_weighted, + unsigned long long *__restrict__ roi_panel_v, + unsigned long long *__restrict__ roi_panel_x, + unsigned long long *__restrict__ roi_panel_y, int *__restrict__ roi_max) { extern __shared__ unsigned long long shared[]; unsigned long long *s_sum = shared; unsigned long long *s_sum2 = &s_sum[roi_count]; unsigned long long *s_pixels = &s_sum2[roi_count]; - unsigned long long *s_xw = &s_pixels[roi_count]; - unsigned long long *s_yw = &s_xw[roi_count]; - int *s_max = (int *) &s_yw[roi_count]; + unsigned long long *s_pv = &s_pixels[roi_count]; + unsigned long long *s_xw = &s_pv[roi_count * slots]; + unsigned long long *s_yw = &s_xw[roi_count * slots]; + int *s_max = (int *) &s_yw[roi_count * slots]; for (int r = threadIdx.x; r < roi_count; r += blockDim.x) { s_sum[r] = 0; s_sum2[r] = 0; s_pixels[r] = 0; - s_xw[r] = 0; - s_yw[r] = 0; s_max[r] = INT_MIN; } + for (int k = threadIdx.x; k < roi_count * slots; k += blockDim.x) { + s_pv[k] = 0; + s_xw[k] = 0; + s_yw[k] = 0; + } __syncthreads(); for (size_t idx = blockIdx.x * blockDim.x + threadIdx.x; @@ -59,6 +69,16 @@ void gpu_roi( const long long val = v; const long long x = idx % width; const long long y = idx / width; + int panel = 0; + if (npanels > 0) { + panel = -1; + for (int p = 0; p < npanels && panel < 0; p++) + if (x >= panel_box[4 * p] && x < panel_box[4 * p] + panel_box[4 * p + 2] + && y >= panel_box[4 * p + 1] && y < panel_box[4 * p + 1] + panel_box[4 * p + 3]) + panel = p; + if (panel < 0) + continue; + } const unsigned long long val_u = (unsigned long long) val; const unsigned long long val2_u = (unsigned long long) (val * val); const unsigned long long vx_u = (unsigned long long) (val * x); @@ -71,8 +91,9 @@ void gpu_roi( atomicAdd(&s_sum[r], val_u); atomicAdd(&s_sum2[r], val2_u); atomicAdd(&s_pixels[r], 1ULL); - atomicAdd(&s_xw[r], vx_u); - atomicAdd(&s_yw[r], vy_u); + atomicAdd(&s_pv[r * slots + panel], val_u); + atomicAdd(&s_xw[r * slots + panel], vx_u); + atomicAdd(&s_yw[r * slots + panel], vy_u); } atomicMax(&s_max[r], v); } @@ -83,10 +104,13 @@ void gpu_roi( atomicAdd(&roi_sum[r], s_sum[r]); atomicAdd(&roi_sum2[r], s_sum2[r]); atomicAdd(&roi_pixels[r], s_pixels[r]); - atomicAdd(&roi_x_weighted[r], s_xw[r]); - atomicAdd(&roi_y_weighted[r], s_yw[r]); atomicMax(&roi_max[r], s_max[r]); } + for (int k = threadIdx.x; k < roi_count * slots; k += blockDim.x) { + atomicAdd(&roi_panel_v[k], s_pv[k]); + atomicAdd(&roi_panel_x[k], s_xw[k]); + atomicAdd(&roi_panel_y[k], s_yw[k]); + } } ROIIntegrationGPU::ROIIntegrationGPU(const DiffractionExperiment &experiment, std::shared_ptr stream) @@ -96,14 +120,16 @@ ROIIntegrationGPU::ROIIntegrationGPU(const DiffractionExperiment &experiment, st gpu_sum(roi_count), gpu_sum2(roi_count), gpu_pixels(roi_count), - gpu_x_weighted(roi_count), - gpu_y_weighted(roi_count), + gpu_panel_v(roi_count * panel_slots), + gpu_panel_x(roi_count * panel_slots), + gpu_panel_y(roi_count * panel_slots), gpu_max(roi_count), host_sum(roi_count), host_sum2(roi_count), host_pixels(roi_count), - host_x_weighted(roi_count), - host_y_weighted(roi_count), + host_panel_v(roi_count * panel_slots), + host_panel_x(roi_count * panel_slots), + host_panel_y(roi_count * panel_slots), host_max(roi_count), max_init(roi_count, INT_MIN) { @@ -115,7 +141,18 @@ ROIIntegrationGPU::ROIIntegrationGPU(const DiffractionExperiment &experiment, st cuda_err(cudaGetDeviceProperties(&prop, device)); threads = 128; blocks = 4 * prop.multiProcessorCount; - shared_needed = roi_count * (5 * sizeof(unsigned long long) + sizeof(int)); + shared_needed = roi_count * (3 * sizeof(unsigned long long) + sizeof(int)) + + roi_count * panel_slots * 3 * sizeof(unsigned long long); + + std::vector box; + for (const auto &m: modules.Get()) + box.insert(box.end(), {static_cast(m.x0), static_cast(m.y0), static_cast(m.width), + static_cast(m.height)}); + if (!box.empty()) { + gpu_panel_box = CudaDevicePtr(box.size()); + cuda_err(cudaMemcpyAsync(gpu_panel_box, box.data(), sizeof(int) * box.size(), + cudaMemcpyHostToDevice, *stream)); + } // On this engine's stream, like every other operation it issues: the streams are non-blocking, so a // NULL-stream copy is no longer ordered against the kernels that read the map. The one-time @@ -133,20 +170,24 @@ void ROIIntegrationGPU::Run(const ImagePreprocessorBuffer &image, std::map>>( gpu_roi_map, image.getGPUBuffer(), npixel, width, roi_count, - gpu_sum, gpu_sum2, gpu_pixels, gpu_x_weighted, gpu_y_weighted, gpu_max); + gpu_panel_box, static_cast(modules.size()), static_cast(panel_slots), + gpu_sum, gpu_sum2, gpu_pixels, gpu_panel_v, gpu_panel_x, gpu_panel_y, gpu_max); cuda_err(cudaGetLastError()); cudaMemcpyAsync(host_sum.data(), gpu_sum, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream); cudaMemcpyAsync(host_sum2.data(), gpu_sum2, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream); cudaMemcpyAsync(host_pixels.data(), gpu_pixels, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream); - cudaMemcpyAsync(host_x_weighted.data(), gpu_x_weighted, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream); - cudaMemcpyAsync(host_y_weighted.data(), gpu_y_weighted, sizeof(unsigned long long) * roi_count, cudaMemcpyDeviceToHost, *stream); + cudaMemcpyAsync(host_panel_v.data(), gpu_panel_v, sizeof(unsigned long long) * slots, cudaMemcpyDeviceToHost, *stream); + cudaMemcpyAsync(host_panel_x.data(), gpu_panel_x, sizeof(unsigned long long) * slots, cudaMemcpyDeviceToHost, *stream); + cudaMemcpyAsync(host_panel_y.data(), gpu_panel_y, sizeof(unsigned long long) * slots, cudaMemcpyDeviceToHost, *stream); cudaMemcpyAsync(host_max.data(), gpu_max, sizeof(int) * roi_count, cudaMemcpyDeviceToHost, *stream); cuda_err(cudaStreamSynchronize(*stream)); @@ -154,9 +195,13 @@ void ROIIntegrationGPU::Run(const ImagePreprocessorBuffer &image, std::map(host_sum[r]); roi_sum2[r] = host_sum2[r]; roi_pixels[r] = host_pixels[r]; - roi_x_weighted[r] = static_cast(host_x_weighted[r]); - roi_y_weighted[r] = static_cast(host_y_weighted[r]); roi_max[r] = (host_max[r] == INT_MIN) ? INT64_MIN : static_cast(host_max[r]); } + for (size_t k = 0; k < slots; k++) { + roi_panel_v[k] = static_cast(host_panel_v[k]); + roi_panel_x[k] = static_cast(host_panel_x[k]); + roi_panel_y[k] = static_cast(host_panel_y[k]); + } + PlaceWeighted(); Export(out); } diff --git a/image_analysis/roi/ROIIntegrationGPU.h b/image_analysis/roi/ROIIntegrationGPU.h index 1e0696ce8..b13bde79f 100644 --- a/image_analysis/roi/ROIIntegrationGPU.h +++ b/image_analysis/roi/ROIIntegrationGPU.h @@ -20,15 +20,18 @@ class ROIIntegrationGPU : public ROIIntegration { CudaDevicePtr gpu_sum; CudaDevicePtr gpu_sum2; CudaDevicePtr gpu_pixels; - CudaDevicePtr gpu_x_weighted; - CudaDevicePtr gpu_y_weighted; + CudaDevicePtr gpu_panel_v; // per ROI and panel, as roi_panel_v + CudaDevicePtr gpu_panel_x; + CudaDevicePtr gpu_panel_y; + CudaDevicePtr gpu_panel_box; // x0, y0, width, height per panel; unallocated without panels CudaDevicePtr gpu_max; std::vector host_sum; std::vector host_sum2; std::vector host_pixels; - std::vector host_x_weighted; - std::vector host_y_weighted; + std::vector host_panel_v; + std::vector host_panel_x; + std::vector host_panel_y; std::vector host_max; std::vector max_init; // INT_MIN seed copied into gpu_max each frame public: diff --git a/tests/BeamCenterFromBackgroundTest.cpp b/tests/BeamCenterFromBackgroundTest.cpp index d38db14e2..270004b77 100644 --- a/tests/BeamCenterFromBackgroundTest.cpp +++ b/tests/BeamCenterFromBackgroundTest.cpp @@ -308,6 +308,45 @@ TEST_CASE("FindBeamCenter_BlanksTheBeamStopOutOfTheCapture", "[BeamCenter]") { CHECK(masked_error <= unmasked_error); } +// The pixels are binned where the panel table puts them on the plane. Every module of the image moved by +// the same fractional step: the background is isotropic about the beam on the plane, and the walk, which +// works on the plane, finds the beam there - not where it would be on the image. +TEST_CASE("BeamCenterFromBackground_BinsOnThePlane", "[BeamCenter]") { + DiffractionExperiment x = TestExperiment(); + std::vector moved = ModulesAsImaged(x.GetDetectorSetup().GetGeometry()).Get(); + REQUIRE(!moved.empty()); + for (auto &m: moved) { + m.ox += 20.4; + m.oy -= 11.3; + } + x.Modules(DetectorModules(moved)); + PixelMask pixel_mask(x); + const auto &mask = pixel_mask.GetMask(x); + + const DiffractionGeometry geom_true = OffsetBy(x.GetDiffractionGeometry(), 3.0f, -2.5f); + std::vector projection(x.GetPixelsNum(), NAN); + std::mt19937 rng(20261021); + std::normal_distribution gauss(0.0f, 1.0f); + x.GetModules().ForEachPixel(x.GetXPixelsNumConv(), x.GetYPixelsNumConv(), + [&](int64_t, size_t i, double px, double py) { + if (mask[i] != 0) + return; + const float value = background(geom_true.TwoTheta_rad(static_cast(px), static_cast(py)), 60.0f); + projection[i] = value + gauss(rng) * std::sqrt(value / 60.0f); + }); + + const auto host = FindBeamCenterFromBackground(x, pixel_mask, projection, 0, {}, /*allow_device=*/false); + REQUIRE(host.has_value()); + CHECK(host->beam_x_pxl == Catch::Approx(geom_true.GetBeamX_pxl()).margin(0.5)); + CHECK(host->beam_y_pxl == Catch::Approx(geom_true.GetBeamY_pxl()).margin(0.5)); + + // The device places the pixels the same way, to the bit. + const auto device = FindBeamCenterFromBackground(x, pixel_mask, projection, 0, {}, /*allow_device=*/true); + REQUIRE(device.has_value()); + CHECK(device->beam_x_pxl == host->beam_x_pxl); + CHECK(device->beam_y_pxl == host->beam_y_pxl); +} + #ifdef JFJOCH_USE_CUDA #include "../common/CUDAWrapper.h" diff --git a/tests/PanelLayoutTest.cpp b/tests/PanelLayoutTest.cpp index 4c8f5a484..e4e1efbc6 100644 --- a/tests/PanelLayoutTest.cpp +++ b/tests/PanelLayoutTest.cpp @@ -11,7 +11,7 @@ // [PanelLayout] cases check what is already layout-correct. The [PanelLayoutTodo] cases are hidden: they // state the equivalence each later stage has to deliver, and fail until it does. A stage that converts a // consumer retags its case [PanelLayout]: -// stage 1 - dense maps: PanelLayout_AzintMapping, PanelLayout_ResolutionMask +// stage 1 - dense maps: PanelLayout_AzintMapping, PanelLayout_ResolutionMask (done) // 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, @@ -32,6 +32,7 @@ #include "../common/Reflection.h" #include "../image_analysis/bragg_integration/BraggIntegrationEngineCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" +#include "../image_analysis/roi/ROIIntegrationCPU.h" #include "../image_analysis/spot_finding/ImageSpotFinderCPU.h" #include "../image_analysis/spot_finding/SpotUtils.h" @@ -262,9 +263,10 @@ TEST_CASE("PanelLayout_IntegrationRoundTrip", "[PanelLayout]") { } } -// TODO stage 1: the azimuthal-integration map built per stored pixel through the panel table - same bin, -// resolution and correction as the assembled pixel it lands on. -TEST_CASE("PanelLayout_AzintMapping", "[.][PanelLayoutTodo]") { +// The azimuthal-integration map built per stored pixel through the panel table - same bin, resolution and +// correction as the assembled pixel it lands on. The same for the azimuthal ROI map, and the ROI centre of +// mass comes out on the plane. +TEST_CASE("PanelLayout_AzintMapping", "[PanelLayout]") { const Scene scene(INTERIOR_SPOTS); const AzimuthalIntegrationMapping a(scene.assembled, scene.assembled_mask); const AzimuthalIntegrationMapping s(scene.stacked, scene.stacked_mask); @@ -278,10 +280,30 @@ TEST_CASE("PanelLayout_AzintMapping", "[.][PanelLayoutTodo]") { CHECK(bin == 0); CHECK(resolution == 0); CHECK(correction == 0); + + DiffractionExperiment assembled = scene.assembled, stacked = scene.stacked; + for (auto *x: {&assembled, &stacked}) + x->ROI().SetROI(ROIDefinition{.azimuthal = {ROIAzimuthal("ring", 3.0f, 8.0f, 20.0f, 250.0f)}}); + const auto a_roi = assembled.ExportROIMap(), s_roi = stacked.ExportROIMap(); + size_t roi = 0, in_roi = 0; + scene.modules.ForEachPixel(PW, NPANELS * PH, [&](int64_t, size_t i, double px, double py) { + roi += s_roi[i] != a_roi[Scene::AssembledIndex(px, py)]; + in_roi += s_roi[i] != 0; + }); + CHECK(roi == 0); + CHECK(in_roi > 1000); + + std::map a_out, s_out; + ROIIntegrationCPU(assembled).Run(scene.assembled_image, a_out); + ROIIntegrationCPU(stacked).Run(scene.stacked_image, s_out); + CHECK(s_out.at("ring").sum == a_out.at("ring").sum); + CHECK(s_out.at("ring").pixels == a_out.at("ring").pixels); + CHECK(s_out.at("ring").x_weighted == a_out.at("ring").x_weighted); + CHECK(s_out.at("ring").y_weighted == a_out.at("ring").y_weighted); } -// TODO stage 1: the spot finders' resolution mask, per stored pixel through the panel table. -TEST_CASE("PanelLayout_ResolutionMask", "[.][PanelLayoutTodo]") { +// The spot finders' resolution mask, per stored pixel through the panel table. +TEST_CASE("PanelLayout_ResolutionMask", "[PanelLayout]") { const Scene scene(INTERIOR_SPOTS); const auto a = AzimuthalIntegrationMapping(scene.assembled, scene.assembled_mask).ResolutionMaskBits(3.0f, 8.0f); const auto s = AzimuthalIntegrationMapping(scene.stacked, scene.stacked_mask).ResolutionMaskBits(3.0f, 8.0f); diff --git a/tests/ROIIntegrationGPUTest.cpp b/tests/ROIIntegrationGPUTest.cpp index dc6d28412..9e134315b 100644 --- a/tests/ROIIntegrationGPUTest.cpp +++ b/tests/ROIIntegrationGPUTest.cpp @@ -46,22 +46,12 @@ void compare_results(const std::map &cpu, CHECK(g.pixels_masked == c.pixels_masked); } } -} // namespace -// The GPU kernel reduces with atomics and two's-complement unsigned accumulators, -// while the CPU path is a plain serial loop. On identical input every per-ROI -// statistic must be bit-for-bit identical, so we run both and compare. -TEST_CASE("ROIIntegrationGPU_MatchesCPU") { - if (get_gpu_count() == 0) { - WARN("No CUDA GPU present. Skipping ROIIntegrationGPU_MatchesCPU"); - return; - } - - const DiffractionExperiment experiment = make_roi_experiment(); +void run_and_compare(const DiffractionExperiment &experiment) { const auto roi_map = experiment.ExportROIMap(); const size_t npixel = roi_map.size(); const uint16_t roi_count = experiment.ROI().size(); - REQUIRE(roi_count == 4); + REQUIRE(roi_count > 0); // Deterministic image with both positive and negative values; negatives exercise // the signed weighted-sum path (val * x can be negative). @@ -110,5 +100,37 @@ TEST_CASE("ROIIntegrationGPU_MatchesCPU") { compare_results(out_cpu, out_gpu); } +} // namespace + +// The GPU kernel reduces with atomics and two's-complement unsigned accumulators, +// while the CPU path is a plain serial loop. On identical input every per-ROI +// statistic must be bit-for-bit identical, so we run both and compare. +TEST_CASE("ROIIntegrationGPU_MatchesCPU") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping ROIIntegrationGPU_MatchesCPU"); + return; + } + run_and_compare(make_roi_experiment()); +} + +// The same with a panel table: the centre of mass is summed per panel and placed on the plane. Two panels +// stacked in the image, the second turned by 180 degrees and moved by a fraction of a pixel; a strip of +// rows below them is on no panel. +TEST_CASE("ROIIntegrationGPU_MatchesCPU_Panels") { + if (get_gpu_count() == 0) { + WARN("No CUDA GPU present. Skipping ROIIntegrationGPU_MatchesCPU_Panels"); + return; + } + DetectorModule a = DetectorModule::AsImaged(0, 0, 256, 128); + DetectorModule b = DetectorModule::AsImaged(0, 128, 256, 128); + b.ox = 300.25; b.oy = 400.5; b.fx = -1; b.sy = -1; + DiffractionExperiment experiment(DetDECTRIS(256, 300, "panels", "")); + experiment.Modules(DetectorModules({a, b})); + experiment.ROI().SetROI(ROIDefinition{.boxes = { + ROIBox("roiA", 10, 210, 20, 220), + ROIBox("roiB", 0, 255, 0, 299), + }}); + run_and_compare(experiment); +} #endif diff --git a/tests/ROIMapTest.cpp b/tests/ROIMapTest.cpp index 1c4b52c89..d4a35d2e5 100644 --- a/tests/ROIMapTest.cpp +++ b/tests/ROIMapTest.cpp @@ -277,7 +277,7 @@ TEST_CASE("ROIMap") { // Check arrays are proper size REQUIRE(mask.GetROINameMap().size() == 2); - auto map = mask.GetROIMap(x.GetDiffractionGeometry(), + auto map = mask.GetROIMap(x.GetDiffractionGeometry(), x.GetModules(), x.GetXPixelsNum(), x.GetYPixelsNum());