PredictSpotsOnDetector: Add (but no correctness test so far)

image_analysis/CMakeLists.txt

@@ -4,7 +4,7 @@ ADD_LIBRARY(ImageAnalysis STATIC
         GPUImageAnalysis.h
         RadialIntegration.cpp RadialIntegration.h
         RadialIntegrationMapping.cpp RadialIntegrationMapping.h
-        StrongPixelSet.cpp StrongPixelSet.h GPUImageAnalysis.cpp RadialIntegrationProfile.cpp RadialIntegrationProfile.h)
+        StrongPixelSet.cpp StrongPixelSet.h GPUImageAnalysis.cpp RadialIntegrationProfile.cpp RadialIntegrationProfile.h PredictSpotsOnDetector.h)
 
 TARGET_LINK_LIBRARIES(ImageAnalysis CommonFunctions)
 

image_analysis/PredictSpotsOnDetector.h

@@ -0,0 +1,59 @@
+// Copyright (2019-2023) Paul Scherrer Institute
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#ifndef JUNGFRAUJOCH_PREDICTSPOTSONDETECTOR_H
+#define JUNGFRAUJOCH_PREDICTSPOTSONDETECTOR_H
+
+#include <vector>
+#include "../common/DiffractionGeometry.h"
+#include "../common/ROIFilter.h"
+#include "CrystalLattice.h"
+
+std::vector<std::pair<float, float>> PredictSpotsOnDetector(const DiffractionExperiment& experiment,
+                                                            const CrystalLattice& lattice,
+                                                            int32_t max_hkl = 30,
+                                                            float epsilon = 4e-4) {
+    CrystalLattice recip_l = lattice.ReciprocalLattice();
+    std::vector<std::pair<float, float>> ret;
+    for (int h = -max_hkl; h < max_hkl; h++) {
+        for (int k = -max_hkl; k < max_hkl; k++) {
+            for (int l = -max_hkl; l < max_hkl; l++) {
+                Coord recip = static_cast<float>(h) * recip_l.Vec0()
+                              + static_cast<float>(k) * recip_l.Vec1()
+                              + static_cast<float>(l) * recip_l.Vec2();
+
+                if (DistFromEwaldSphere(experiment, recip) < epsilon)
+                    ret.push_back(RecipToDector(experiment, recip));
+            }
+        }
+    }
+    return ret;
+}
+
+void PredictSpotsOnDetector(ROIFilter &filter,
+                            const DiffractionExperiment& experiment,
+                            const CrystalLattice& lattice,
+                            int32_t max_hkl = 30,
+                            float epsilon = 4e-4,
+                            uint16_t box_size = 7) {
+    CrystalLattice recip_l = lattice.ReciprocalLattice();
+    std::vector<std::pair<float, float>> ret;
+    for (int h = -max_hkl; h < max_hkl; h++) {
+        for (int k = -max_hkl; k < max_hkl; k++) {
+            for (int l = -max_hkl; l < max_hkl; l++) {
+                Coord recip = static_cast<float>(h) * recip_l.Vec0()
+                              + static_cast<float>(k) * recip_l.Vec1()
+                              + static_cast<float>(l) * recip_l.Vec2();
+
+                if (DistFromEwaldSphere(experiment, recip) < epsilon) {
+                    auto [x,y] = RecipToDector(experiment, recip);
+                    auto x0 = static_cast<int32_t>(std::lroundf(x - static_cast<float>(box_size)));
+                    auto y0 = static_cast<int32_t>(std::lroundf(y - static_cast<float>(box_size)));
+                    filter.SetRectangle(x0, y0, 2 * box_size + 1, 2 * box_size + 1);
+                }
+            }
+        }
+    }
+}
+
+#endif //JUNGFRAUJOCH_PREDICTSPOTSONDETECTOR_H

tools/DataAnalysisPerfTest.cpp

@@ -13,6 +13,7 @@
 
 #include "../image_analysis/RadialIntegration.h"
 #include "../common/Logger.h"
+#include "../image_analysis/PredictSpotsOnDetector.h"
 
 #define make_unit_cell(a1,a2,a3,a4,a5,a6) UnitCell{.a = a1, .b = a2, .c = a3, .alpha = a4, .beta = a5, .gamma = a6}
 
@@ -50,6 +51,46 @@ auto TestAll(const DiffractionExperiment &experiment, const JFJochProtoBuf::Data
     return strstream.str();
 }
 
+auto TestAllWithROI(const DiffractionExperiment &experiment, const JFJochProtoBuf::DataProcessingSettings &settings,
+             GPUImageAnalysis &spot_finder, int16_t* image, size_t nimages) {
+    IndexerWrapper indexer;
+    indexer.Setup(experiment.GetUnitCell());
+
+    spot_finder.SetInputBuffer(image);
+
+    std::vector<int16_t> roi_image(experiment.GetPixelsNum());
+
+    auto start_time = std::chrono::system_clock::now();
+    for (int i = 0; i < nimages; i++) {
+        ROIFilter filter(experiment.GetXPixelsNum(), experiment.GetYPixelsNum());
+
+        std::vector<DiffractionSpot> spots;
+        std::vector<float> result;
+        spot_finder.LoadDataToGPU();
+        spot_finder.RunSpotFinder(settings);
+        spot_finder.GetSpotFinderResults(experiment, settings, spots);
+        spot_finder.RunRadialIntegration();
+        std::vector<Coord> recip;
+        for (const auto& s: spots)
+            recip.emplace_back(s.ReciprocalCoord(experiment));
+        auto indexer_ret = indexer.Run(recip);
+        spot_finder.GetRadialIntegrationProfile(result);
+        if (!indexer_ret.empty()) {
+            PredictSpotsOnDetector(filter, experiment, indexer_ret[0].l);
+            filter.Apply<int16_t>(roi_image, INT16_MIN);
+        }
+    }
+
+    auto end_time = std::chrono::system_clock::now();
+    auto elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time);
+
+    std::ostringstream strstream;
+    logger.Info("{:20s} {:8.1f} ms/image", "Full+ROI",
+                elapsed.count() / (1000.0 * (double) nimages));
+
+    return strstream.str();
+}
+
 void TestIndexing() {
     constexpr const int nexec  = 5;
 
@@ -305,5 +346,6 @@ int main(int argc, char **argv) {
     if (GPUImageAnalysis::GPUPresent()) {
         GPUImageAnalysis local_peakfinder_gpu(x.GetXPixelsNum(), x.GetYPixelsNum(), one_byte_mask, mapping);
         TestAll(x, settings, local_peakfinder_gpu,image_conv.data(), nimages);
+        TestAllWithROI(x, settings, local_peakfinder_gpu,image_conv.data(), nimages);
     }
 }