diff --git a/common/DetectorModules.cpp b/common/DetectorModules.cpp index 9e06edbb1..7b313ed65 100644 --- a/common/DetectorModules.cpp +++ b/common/DetectorModules.cpp @@ -112,6 +112,47 @@ std::pair DetectorModules::PlaneToImage(double px, double py) co return modules[m].PlaneToImage(px, py); } +std::pair DetectorModules::ImageToPlane(int64_t panel, double x, double y) const { + if (panel < 0) + return {x, y}; + return modules.at(panel).ImageToPlane(x, y); +} + +std::pair DetectorModules::PlaneToImage(int64_t panel, double px, double py) const { + if (panel < 0) + return {px, py}; + return modules.at(panel).PlaneToImage(px, py); +} + +std::optional DetectorModules::PlaneToPanel(double px, double py) const { + if (modules.empty()) + return PanelPoint{-1, px, py}; + for (size_t i = 0; i < modules.size(); i++) { + const auto [x, y] = modules[i].PlaneToImage(px, py); + if (modules[i].ContainsImagePixel(x, y)) + return PanelPoint{static_cast(i), x, y}; + } + return std::nullopt; +} + +PlaneBox DetectorModules::PlaneBounds(int64_t width, int64_t height) const { + if (modules.empty()) + return {0.0, 0.0, static_cast(width), static_cast(height)}; + PlaneBox box{INFINITY, INFINITY, -INFINITY, -INFINITY}; + for (const auto &m: modules) { + const double xa = static_cast(m.x0) - 0.5, xb = static_cast(m.x0 + m.width) - 0.5; + const double ya = static_cast(m.y0) - 0.5, yb = static_cast(m.y0 + m.height) - 0.5; + for (const auto &[x, y]: {std::pair{xa, ya}, std::pair{xb, ya}, std::pair{xa, yb}, std::pair{xb, yb}}) { + const auto [px, py] = m.ImageToPlane(x, y); + box.x0 = std::min(box.x0, px); + box.y0 = std::min(box.y0, py); + box.x1 = std::max(box.x1, px); + box.y1 = std::max(box.y1, py); + } + } + return box; +} + void DetectorModules::Apply(const std::vector &shifts) { for (auto &m: modules) { const auto [cx, cy] = m.PlaneCentre(); diff --git a/common/DetectorModules.h b/common/DetectorModules.h index fcadb5d09..d38649d4f 100644 --- a/common/DetectorModules.h +++ b/common/DetectorModules.h @@ -5,6 +5,7 @@ #include #include +#include #include #include @@ -39,7 +40,20 @@ struct DetectorModule { void Apply(const ModuleShift &shift); }; +// A point in the image rugnux reads, on panel `panel` (an index into DetectorModules; -1 when the table is +// empty and the whole image is one flat panel). +struct PanelPoint { + int64_t panel; + double x, y; +}; + +// An axis-aligned box on the detector plane, x0 <= x < x1, y0 <= y < y1. +struct PlaneBox { + double x0, y0, x1, y1; +}; + // The modules of one detector. Empty means the image is one flat detector: every mapping is the identity. +// A module is a panel: one rigid box of the image with one placement. Panel ids are indices into the table. class DetectorModules { std::vector modules; public: @@ -64,6 +78,39 @@ public: // Plane -> image; a point on no placed module is left as it is. [[nodiscard]] std::pair PlaneToImage(double px, double py) const; + // Image <-> plane on a known panel, no search. Panel -1 (the empty table) is the identity. + [[nodiscard]] std::pair ImageToPlane(int64_t panel, double x, double y) const; + [[nodiscard]] std::pair PlaneToImage(int64_t panel, double px, double py) const; + // The panel whose placed footprint holds a plane point, and where that point is in the image; nullopt in + // a gap or off every panel. Empty table: panel -1 and the point itself, wherever it is - the caller + // bounds it by the image, as without panels. + [[nodiscard]] std::optional PlaneToPanel(double px, double py) const; + // Box on the plane enclosing every panel's pixels (their footprints, half a pixel beyond the edge pixel + // centres): a quick reject before PlaneToPanel. Empty table: [0, width) x [0, height), the image + // rectangle the bounds tests without panels use. + [[nodiscard]] PlaneBox PlaneBounds(int64_t width, int64_t height) const; + + // Calls f(panel, index, px, py) for every image pixel on a panel - panel by panel, row by row - with + // index = y * width + x in the image and (px, py) its centre on the plane. Pixels on no panel (gaps) are + // not visited. Empty table: one panel -1 that is the whole image, plane == image. + template + void ForEachPixel(int64_t width, int64_t height, F &&f) const { + if (modules.empty()) { + for (int64_t y = 0; y < height; 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; + } + for (size_t p = 0; p < modules.size(); p++) { + const DetectorModule &m = modules[p]; + for (int64_t y = m.y0; y < m.y0 + m.height; y++) + for (int64_t x = m.x0; x < m.x0 + m.width; x++) { + const auto [px, py] = m.ImageToPlane(static_cast(x), static_cast(y)); + f(static_cast(p), static_cast(y * width + x), px, py); + } + } + } + // Each module moved by the shift whose centre is nearest its own plane centre (and within it). void Apply(const std::vector &shifts); // The modules as refinement starting points: their plane centres, no move. diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp index 8f96fe036..bdf16f9ca 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp +++ b/image_analysis/bragg_integration/BraggIntegrationEngineCPU.cpp @@ -78,7 +78,11 @@ BraggIntegrationEngineCPU::BraggIntegrationEngineCPU(const DiffractionExperiment : BraggIntegrationEngine(experiment), refl_mask(static_cast(xpixel), static_cast(ypixel), 0), owner(static_cast(xpixel), static_cast(ypixel), BRAGG_OWNER_NONE), - static_mask(mask.GetPackedMask()) {} + static_mask(mask.GetPackedMask()) { + if (mask.GetMask().size() != npixel) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Pixel mask size doesn't match the image"); +} void BraggIntegrationEngineCPU::MarkNeighbours(const std::vector &predicted, const std::vector &order) { diff --git a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu index 08215325b..c024db81c 100644 --- a/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu +++ b/image_analysis/bragg_integration/BraggIntegrationEngineGPU.cu @@ -834,6 +834,9 @@ BraggIntegrationEngineGPU::BraggIntegrationEngineGPU(const DiffractionExperiment d_global_n(1), d_invd2(2), d_counts(COUNT_SLOTS) { + if (mask.GetMask().size() != npixel) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Pixel mask size doesn't match the image"); threads = 128; // On the engine's stream (the member: the parameter was moved from), where the kernels counting // into it run - the NULL stream is not ordered before those. diff --git a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp index 9a344cd9a..85247cbd6 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp +++ b/image_analysis/image_preprocessing/ImagePreprocessorCPU.cpp @@ -9,7 +9,11 @@ ImagePreprocessorCPU::ImagePreprocessorCPU(const DiffractionExperiment &experiment, const PixelMask &mask) : ImagePreprocessor(experiment), - mask_bits(mask.GetPackedMask()) {} + mask_bits(mask.GetPackedMask()) { + if (mask.GetMask().size() != npixels) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Pixel mask size doesn't match the image"); +} ImageStatistics ImagePreprocessorCPU::Analyze(ImagePreprocessorBuffer &processed_image, const uint8_t *image_ptr, CompressedImageMode image_mode) { if (processed_image.size() != npixels) diff --git a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu index 273afbc23..c7ccb3793 100644 --- a/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu +++ b/image_analysis/image_preprocessing/ImagePreprocessorGPU.cu @@ -329,6 +329,9 @@ ImagePreprocessorGPU::ImagePreprocessorGPU(const DiffractionExperiment &experime copy_image_to_host(copy_image_to_host), gpu_stats(1), cpu_stats(1) { + if (mask.GetMask().size() != npixels) + throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, + "Pixel mask size doesn't match the image"); // Setup mask. The byte-per-pixel form and its checksum come from the PixelMask, which derives them // whenever the mask changes: they are the same for every worker, and deriving them walks every // pixel of the detector - 18 million of them on a 16 Mpx one, per engine, with an engine built per