diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 6b927ae9..0ebfd5a4 100644 --- a/image_analysis/CMakeLists.txt +++ b/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) diff --git a/image_analysis/PredictSpotsOnDetector.h b/image_analysis/PredictSpotsOnDetector.h new file mode 100644 index 00000000..37c79640 --- /dev/null +++ b/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 +#include "../common/DiffractionGeometry.h" +#include "../common/ROIFilter.h" +#include "CrystalLattice.h" + +std::vector> PredictSpotsOnDetector(const DiffractionExperiment& experiment, + const CrystalLattice& lattice, + int32_t max_hkl = 30, + float epsilon = 4e-4) { + CrystalLattice recip_l = lattice.ReciprocalLattice(); + std::vector> 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(h) * recip_l.Vec0() + + static_cast(k) * recip_l.Vec1() + + static_cast(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> 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(h) * recip_l.Vec0() + + static_cast(k) * recip_l.Vec1() + + static_cast(l) * recip_l.Vec2(); + + if (DistFromEwaldSphere(experiment, recip) < epsilon) { + auto [x,y] = RecipToDector(experiment, recip); + auto x0 = static_cast(std::lroundf(x - static_cast(box_size))); + auto y0 = static_cast(std::lroundf(y - static_cast(box_size))); + filter.SetRectangle(x0, y0, 2 * box_size + 1, 2 * box_size + 1); + } + } + } + } +} + +#endif //JUNGFRAUJOCH_PREDICTSPOTSONDETECTOR_H diff --git a/tools/DataAnalysisPerfTest.cpp b/tools/DataAnalysisPerfTest.cpp index 6c976355..d3c2ac3e 100644 --- a/tools/DataAnalysisPerfTest.cpp +++ b/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 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 spots; + std::vector result; + spot_finder.LoadDataToGPU(); + spot_finder.RunSpotFinder(settings); + spot_finder.GetSpotFinderResults(experiment, settings, spots); + spot_finder.RunRadialIntegration(); + std::vector 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(roi_image, INT16_MIN); + } + } + + auto end_time = std::chrono::system_clock::now(); + auto elapsed = std::chrono::duration_cast(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); } }