v1.0.0-rc.40

image_analysis/QuickIntegrate.cpp

@@ -0,0 +1,120 @@
+// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "QuickIntegrate.h"
+#include <iostream>
+
+QuickIntegrate::QuickIntegrate(const DiffractionExperiment &in_experiment)
+        : experiment(in_experiment) {}
+
+template<class T>
+std::vector<Reflection> QuickIntegrate::IntegrateInternal(const CompressedImage &image, const CrystalLattice &latt,
+                                                          int64_t special_value, int64_t saturation) {
+    auto start = std::chrono::high_resolution_clock::now();
+    std::vector<uint8_t> buffer;
+    auto ptr = reinterpret_cast<const T*>(image.GetUncompressedPtr(buffer));
+
+    std::vector<Reflection> ret;
+
+    Coord Astar = latt.Astar();
+    Coord Bstar = latt.Bstar();
+    Coord Cstar = latt.Cstar();
+
+    DiffractionGeometry geom = experiment.GetDiffractionGeometry();
+    for (int h = -max_value; h < max_value; h++) {
+        for (int k = -max_value; k < max_value; k++) {
+            for (int l = -max_value; l < max_value; l++) {
+                Coord recip = Astar * h + Bstar * k + Cstar * l;
+                float d = 1.0f / recip.Length();
+                float dist_ewald_sphere = geom.DistFromEwaldSphere(recip);
+                if (d > d_min && dist_ewald_sphere < ewald_sphere_dist_cutoff) {
+                    auto [x,y] = geom.RecipToDector(recip);
+                    if ((x >= 0) && (x < image.GetWidth()) && (y >= 0) && (y < image.GetHeight())) {
+                        Reflection r{.h = h, .k = k, .l = l, .center_x = x, .center_y = y, .d = d};
+                        bool observed = IntegrateInternal(r, ptr, image.GetWidth(), image.GetHeight(),
+                                          x, y, special_value, saturation);
+                        if (observed)
+                            ret.push_back(r);
+                    }
+                }
+            }
+        }
+    }
+    return ret;
+}
+
+std::vector<Reflection> QuickIntegrate::Integrate(const CompressedImage& image, const CrystalLattice &latt) {
+    if (image.GetCompressedSize() == 0)
+        return {};
+
+    switch (image.GetMode()) {
+        case CompressedImageMode::Int8:
+            return IntegrateInternal<int8_t>(image, latt, INT8_MIN, INT8_MAX);
+        case CompressedImageMode::Int16:
+            return IntegrateInternal<int16_t>(image, latt, INT16_MIN, INT16_MAX);
+        case CompressedImageMode::Int32:
+            return IntegrateInternal<int32_t>(image, latt, INT32_MIN, INT32_MAX);
+        case CompressedImageMode::Uint8:
+            return IntegrateInternal<uint8_t>(image, latt, UINT8_MAX, UINT8_MAX);
+        case CompressedImageMode::Uint16:
+            return IntegrateInternal<uint16_t>(image, latt, UINT16_MAX, UINT16_MAX);
+        case CompressedImageMode::Uint32:
+            return IntegrateInternal<uint16_t>(image, latt, UINT32_MAX, UINT32_MAX);
+        default:
+            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                                  "Image mode not supported");
+    }
+
+}
+
+template<class T>
+bool QuickIntegrate::IntegrateInternal(Reflection &r, const T *image, size_t xpixel, size_t ypixel,
+                                       float center_x, float center_y,
+                                       int64_t special_value, int64_t saturation) {
+    r.I = 0;
+    r.bkg = 0;
+
+    int64_t x0 = std::floor(center_x - r_3 - 1.0);
+    int64_t x1 = std::ceil(center_x + r_3 + 1.0);
+    int64_t y0 = std::floor(center_y - r_3 - 1.0);
+    int64_t y1 = std::ceil(center_y + r_3 + 1.0);
+    if (x0 < 0)
+        x0 = 0;
+    if (y0 < 0)
+        y0 = 0;
+    if (x1 >= xpixel)
+        x1 = xpixel - 1;
+    if (y1 >= ypixel)
+        y1 = ypixel - 1;
+
+    int64_t I_sum = 0;
+    int64_t I_npixel = 0;
+    int64_t bkg_sum = 0;
+    int64_t bkg_npixel = 0;
+
+    for (int64_t y = y0; y <= y1; y++) {
+        for (int64_t x = x0; x <= x1; x++) {
+            float dist_sq = (x - center_x) * (x - center_x) + (y - center_y) * (y - center_y);
+            if (image[y * xpixel + x] == special_value || image[y * xpixel + x] == saturation) {
+                continue;
+            } else if (dist_sq < r_1_sq) {
+                I_sum += image[y * xpixel + x];
+                I_npixel++;
+            } else if (dist_sq >= r_2_sq && dist_sq < r_3_sq) {
+                bkg_sum += image[y * xpixel + x];
+                bkg_npixel++;
+            }
+        }
+    }
+
+    if ((I_npixel > 2) && (bkg_npixel > 5)) {
+        r.bkg = static_cast<float>(bkg_sum) / static_cast<float>(bkg_npixel);
+        r.I = static_cast<float>(I_sum) - static_cast<float>(I_npixel) * r.bkg;
+        if (I_sum > 0)
+            r.sigma = std::sqrt(static_cast<float>(I_sum));
+        else
+            r.sigma = 0;
+        return true;
+    }
+    return false;
+}