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<uint32_t> *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<int64_t>(width), static_cast<int64_t>(row0),
+                       static_cast<int64_t>(std::min(row_end, height)),
+                       [&](int64_t, size_t pxl, double px, double py) {
+                           SetupPixel(geom, mask, static_cast<uint32_t>(pxl), static_cast<float>(px),
+                                      static_cast<float>(py));
+                       });
 }
 
 void AzimuthalIntegrationMapping::SetupConvGeom(const DiffractionGeometry &geom, const std::vector<uint32_t> *mask) {
@@ -127,6 +132,15 @@ void AzimuthalIntegrationMapping::SetupConvGeomCached(const DiffractionGeometry
                            static_cast<uint32_t>(settings.IsPolarizationCorrection()),
                            static_cast<uint32_t>(polarization_factor.has_value()),
                            static_cast<uint32_t>(width), static_cast<uint32_t>(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<uint32_t>(i));
+        for (const double d: {m.ox, m.oy, m.fx, m.fy, m.sx, m.sy}) {
+            const auto bits = std::bit_cast<uint64_t>(d);
+            key.insert(key.end(), {static_cast<uint32_t>(bits), static_cast<uint32_t>(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<float>(x), static_cast<float>(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<float>(x), static_cast<float>(y));
                     }
                 }
             }));
@@ -198,12 +212,10 @@ void AzimuthalIntegrationMapping::SetupRawGeom(const DiffractionExperiment &expe
 
 void AzimuthalIntegrationMapping::SetupPixel(const DiffractionGeometry &geom,
                                       const std::vector<uint32_t> *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<float>(col);
-    auto y = static_cast<float>(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<uint32_t> key;   // bit patterns of everything SetupPixel reads
+        std::vector<uint32_t> key;   // bit patterns of everything SetupPixel reads, and the panel table
         std::vector<uint16_t> pixel_to_bin;
         std::vector<float> pixel_resolution;
         std::vector<float> 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<float> bin_to_q;
     std::vector<float> bin_to_2theta;
@@ -76,8 +79,9 @@ protected:
     void SetupConvGeomCached(const DiffractionGeometry &geom, const std::vector<uint32_t> &mask,
                              AzimuthalIntegrationGeometryCache &cache);
 
+    // Pixel pxl of the image, at (x, y) on the plane.
     void SetupPixel(const DiffractionGeometry &geom, const std::vector<uint32_t> *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<uint16_t> DiffractionExperiment::ExportROIMap() const {
-    return roi_mask.GetROIMap(GetDiffractionGeometry(), GetXPixelsNumConv(), GetYPixelsNumConv());
+    return roi_mask.GetROIMap(GetDiffractionGeometry(), GetModules(), GetXPixelsNumConv(), GetYPixelsNumConv());
 }
 
 DiffractionExperiment &DiffractionExperiment::ImagesPerFile(const std::optional<int64_t> &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 <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
+#include <vector>
+
+#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<class F>
+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<size_t>(y * width + x), static_cast<double>(x), static_cast<double>(y));
+        return;
+    }
+    const auto &m = modules.Get();
+    std::vector<int64_t> 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<int64_t>(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<double>(x), static_cast<double>(y));
+                f(p, static_cast<size_t>(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<uint16_t> ROIMap::GetROIMap(const DiffractionGeometry &geom, size_t xpixel, size_t ypixel) const {
+std::vector<uint16_t> ROIMap::GetROIMap(const DiffractionGeometry &geom, const DetectorModules &modules,
+                                        size_t xpixel, size_t ypixel) const {
     std::vector<uint16_t> map(xpixel * ypixel, 0);
     std::vector<float> resolution_map(xpixel * ypixel, 0);
     std::vector<float> 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<float>(x), static_cast<float>(y));
-                phi_map[y * xpixel + x] = geom.Phi_rad(static_cast<float>(x), static_cast<float>(y))
-                                          * 180.0f / static_cast<float>(PI);
-            }
-        }
+        // A pixel on no panel keeps resolution and angle 0.
+        ForEachPixelInRows(modules, static_cast<int64_t>(xpixel), 0, static_cast<int64_t>(ypixel),
+                           [&](int64_t, size_t i, double px, double py) {
+                               const auto x = static_cast<float>(px), y = static_cast<float>(py);
+                               resolution_map[i] = geom.PxlToRes(x, y);
+                               phi_map[i] = geom.Phi_rad(x, y) * 180.0f / static_cast<float>(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<uint16_t> 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<uint16_t> GetROIMap(const DiffractionGeometry &geom, const DetectorModules &modules,
+                                                  size_t width, size_t height) const;
     [[nodiscard]] std::map<std::string, uint16_t> 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<float>(p.ox + u * p.fx + v * p.sx);
+    py = static_cast<float>(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<int>(i % width), static_cast<int>(i / width),
-                                         beam_x, beam_y, jx, jy);
-        if (c >= 0)
-            cell = c;
+        const int x = static_cast<int>(i % width), y = static_cast<int>(i / width);
+        float px = static_cast<float>(x), py = static_cast<float>(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<int32_t>(i);
@@ -112,6 +125,8 @@ struct BeamCenterBackgroundGPU::Impl {
     CudaDevicePtr<char> usable;
     CudaDevicePtr<float> mean;
     CudaDevicePtr<int32_t> row_block;
+    CudaDevicePtr<BackgroundPanel> panels;
+    int npanels = 0;
     CudaDevicePtr<int32_t> key, sorted_key, index, sorted_index;
     CudaDevicePtr<float> jac_x, jac_y;
     CudaDevicePtr<int32_t> offset;
@@ -159,8 +174,15 @@ struct BeamCenterBackgroundGPU::Impl {
 };
 
 BeamCenterBackgroundGPU::BeamCenterBackgroundGPU(int width, int height, const std::vector<int> &block_row,
+                                                 const std::vector<BackgroundPanel> &panels,
                                                  const char *usable, const float *mean)
     : impl(std::make_unique<Impl>(width, height)) {
+    impl->npanels = static_cast<int>(panels.size());
+    if (!panels.empty()) {
+        impl->panels = CudaDevicePtr<BackgroundPanel>(panels.size());
+        check(cudaMemcpyAsync(impl->panels.get(), panels.data(), panels.size() * sizeof(BackgroundPanel),
+                              cudaMemcpyHostToDevice, impl->stream), "upload");
+    }
     std::vector<int32_t> 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<int32_t> &count) {
     Impl &d = *impl;
     const auto blocks = static_cast<unsigned>((d.npixels + THREADS - 1) / THREADS);
-    bin_kernel<<<blocks, THREADS, 0, d.stream>>>(band, d.width, d.npixels, beam_x, beam_y, d.usable,
+    bin_kernel<<<blocks, THREADS, 0, d.stream>>>(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> 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<int> &block_row,
-                            const char *usable, const float *mean);
+                            const std::vector<BackgroundPanel> &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<BackgroundPanel> panels;
+    for (const auto &m: modules.Get())
+        panels.push_back({static_cast<int>(m.x0), static_cast<int>(m.y0), static_cast<int>(m.width),
+                          static_cast<int>(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<BeamCenterBackgroundGPU> gpu;
     if (allow_device && get_gpu_count() > 0)
-        gpu = std::make_unique<BeamCenterBackgroundGPU>(W, H, block_row, usable.data(), mean.data());
+        gpu = std::make_unique<BeamCenterBackgroundGPU>(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<size_t>(block_row[b]) * W;
             float *values = band_value.data() + static_cast<size_t>(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<size_t>(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<double>(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<float>(px), y = static_cast<float>(py);
+                if (p >= 0)   // as the device has it
+                    BackgroundPlanePoint(panels[p], static_cast<int>(i % W), static_cast<int>(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<double>(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 <future>
 
@@ -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<int32_t[]>(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<int>(height), nthreads, [&](int y0, int y1) {
-        for (size_t y = y0; y < static_cast<size_t>(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<float>(x), fy = static_cast<float>(y);
-                const int ring = static_cast<int>(geometry.TwoTheta_rad(fx, fy) / ring_rad);
-                const int sector = std::min(SECTORS - 1, static_cast<int>(geometry.Phi_rad(fx, fy)
-                                                                          / static_cast<float>(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<int64_t>(width), y0, y1, [&](int64_t, size_t i, double px, double py) {
+            if (m[i] != 0) return;
+            const float fx = static_cast<float>(px), fy = static_cast<float>(py);
+            const int ring = static_cast<int>(geometry.TwoTheta_rad(fx, fy) / ring_rad);
+            const int sector = std::min(SECTORS - 1, static_cast<int>(geometry.Phi_rad(fx, fy)
+                                                                      / static_cast<float>(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<uint8_t> persistent(width * height, 0);
     ParallelChunks(static_cast<int>(height), nthreads, [&](int y0, int y1) {
-        for (size_t y = y0; y < static_cast<size_t>(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<float>(x), static_cast<float>(y)).Normalize();
-                const double zeta = std::max(1e-4, static_cast<double>(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<double>(std::max<double>(k_bragg, k_chance[key[i] / SECTORS]),
-                                                  n_valid(i));
-                persistent[i] = n_lit[i] >= k;
-            }
+        ForEachPixelInRows(modules, static_cast<int64_t>(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<float>(px), static_cast<float>(py)).Normalize();
+            const double zeta = std::max(1e-4, static_cast<double>(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<double>(std::max<double>(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 <algorithm>
+#include <cmath>
+
 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<size_t>(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<double>(roi_panel_v[k]);
+            const auto u = static_cast<double>(roi_panel_x[k] - m.x0 * roi_panel_v[k]);
+            const auto w = static_cast<double>(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<std::string, ROIMessage> &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 <string>
 #include <vector>
 
+#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<uint16_t> roi_map;   // per-pixel ROI bitmask, assembled (converted) geometry
@@ -34,6 +37,17 @@ protected:
     std::vector<int64_t>  roi_y_weighted;
     std::vector<int64_t>  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<int64_t> roi_panel_v;
+    std::vector<int64_t> roi_panel_x;
+    std::vector<int64_t> roi_panel_y;
+
+    // roi_x_weighted / roi_y_weighted on the plane, from the per-panel sums.
+    void PlaceWeighted();
     void Export(std::map<std::string, ROIMessage> &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 <algorithm>
 #include <limits>
 #include <type_traits>
 
@@ -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<std::remove_reference_t<decltype(image[0])>>;
 
-        for (size_t i = 0; i < npixel; i++) {
+        // Pixels on no panel are not visited.
+        modules.ForEachPixel(static_cast<int64_t>(width), static_cast<int64_t>(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<pixel_t>) {
                 if (v == std::numeric_limits<pixel_t>::min())
-                    continue;
+                    return;
             }
             // saturated pixels still count towards the max, but not the sums
             const bool saturated = (v == std::numeric_limits<pixel_t>::max());
@@ -49,6 +53,7 @@ public:
             const int64_t val = static_cast<int64_t>(v);
             const int64_t x = static_cast<int64_t>(i % width);
             const int64_t y = static_cast<int64_t>(i / width);
+            const size_t slot = panel < 0 ? 0 : static_cast<size_t>(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<uint64_t>(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<CudaStream> 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<int> box;
+    for (const auto &m: modules.Get())
+        box.insert(box.end(), {static_cast<int>(m.x0), static_cast<int>(m.y0), static_cast<int>(m.width),
+                               static_cast<int>(m.height)});
+    if (!box.empty()) {
+        gpu_panel_box = CudaDevicePtr<int>(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<std::
     cuda_err(cudaMemsetAsync(gpu_sum, 0, sizeof(unsigned long long) * roi_count, *stream));
     cuda_err(cudaMemsetAsync(gpu_sum2, 0, sizeof(unsigned long long) * roi_count, *stream));
     cuda_err(cudaMemsetAsync(gpu_pixels, 0, sizeof(unsigned long long) * roi_count, *stream));
-    cuda_err(cudaMemsetAsync(gpu_x_weighted, 0, sizeof(unsigned long long) * roi_count, *stream));
-    cuda_err(cudaMemsetAsync(gpu_y_weighted, 0, sizeof(unsigned long long) * roi_count, *stream));
+    const size_t slots = roi_count * panel_slots;
+    cuda_err(cudaMemsetAsync(gpu_panel_v, 0, sizeof(unsigned long long) * slots, *stream));
+    cuda_err(cudaMemsetAsync(gpu_panel_x, 0, sizeof(unsigned long long) * slots, *stream));
+    cuda_err(cudaMemsetAsync(gpu_panel_y, 0, sizeof(unsigned long long) * slots, *stream));
     cuda_err(cudaMemcpyAsync(gpu_max, max_init.data(), sizeof(int) * roi_count, cudaMemcpyHostToDevice, *stream));
 
     gpu_roi<<<blocks, threads, shared_needed, *stream>>>(
         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<int>(modules.size()), static_cast<int>(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<std::
         roi_sum[r] = static_cast<int64_t>(host_sum[r]);
         roi_sum2[r] = host_sum2[r];
         roi_pixels[r] = host_pixels[r];
-        roi_x_weighted[r] = static_cast<int64_t>(host_x_weighted[r]);
-        roi_y_weighted[r] = static_cast<int64_t>(host_y_weighted[r]);
         roi_max[r] = (host_max[r] == INT_MIN) ? INT64_MIN : static_cast<int64_t>(host_max[r]);
     }
+    for (size_t k = 0; k < slots; k++) {
+        roi_panel_v[k] = static_cast<int64_t>(host_panel_v[k]);
+        roi_panel_x[k] = static_cast<int64_t>(host_panel_x[k]);
+        roi_panel_y[k] = static_cast<int64_t>(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<unsigned long long> gpu_sum;
     CudaDevicePtr<unsigned long long> gpu_sum2;
     CudaDevicePtr<unsigned long long> gpu_pixels;
-    CudaDevicePtr<unsigned long long> gpu_x_weighted;
-    CudaDevicePtr<unsigned long long> gpu_y_weighted;
+    CudaDevicePtr<unsigned long long> gpu_panel_v;   // per ROI and panel, as roi_panel_v
+    CudaDevicePtr<unsigned long long> gpu_panel_x;
+    CudaDevicePtr<unsigned long long> gpu_panel_y;
+    CudaDevicePtr<int> gpu_panel_box;   // x0, y0, width, height per panel; unallocated without panels
     CudaDevicePtr<int> gpu_max;
 
     std::vector<unsigned long long> host_sum;
     std::vector<unsigned long long> host_sum2;
     std::vector<unsigned long long> host_pixels;
-    std::vector<unsigned long long> host_x_weighted;
-    std::vector<unsigned long long> host_y_weighted;
+    std::vector<unsigned long long> host_panel_v;
+    std::vector<unsigned long long> host_panel_x;
+    std::vector<unsigned long long> host_panel_y;
     std::vector<int> host_max;
     std::vector<int> 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<DetectorModule> 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<float> projection(x.GetPixelsNum(), NAN);
+    std::mt19937 rng(20261021);
+    std::normal_distribution<float> 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<float>(px), static_cast<float>(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<std::string, ROIMessage> 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<std::string, ROIMessage> &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());