diff --git a/image_analysis/SpotAnalyze.cpp b/image_analysis/SpotAnalyze.cpp index 1e1b9dfa..d3fd817e 100644 --- a/image_analysis/SpotAnalyze.cpp +++ b/image_analysis/SpotAnalyze.cpp @@ -41,13 +41,22 @@ void SpotAnalyze(const DiffractionExperiment &experiment, output.spots = spots_out; if ((indexer != nullptr) && spot_finding_settings.indexing) { - auto latt_f = indexer->Run(experiment, output); - auto latt = latt_f.get(); + std::vector recip; + recip.reserve(output.spots.size()); + bool index_ice_rings = experiment.GetIndexingSettings().GetIndexIceRings(); + for (const auto &i: output.spots) { + if (index_ice_rings || !i.ice_ring) + recip.push_back(i.ReciprocalCoord(geom)); + } - if (!latt) + auto latt_f = indexer->Run(experiment, recip); + auto indexer_result = latt_f.get(); + output.indexing_time_s = indexer_result.indexing_time_s; + + if (indexer_result.lattice.empty()) output.indexing_result = false; else { - auto uc = latt->GetUnitCell(); + auto latt = indexer_result.lattice[0]; bool beam_center_updated = false; @@ -71,7 +80,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment, }; } } else { - auto sym_result = LatticeSearch(latt.value()); + auto sym_result = LatticeSearch(latt); symmetry = LatticeMessage{ .centering = sym_result.centering, .niggli_class = sym_result.niggli_class, @@ -85,7 +94,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment, DiffractionExperiment experiment_copy(experiment); XtalOptimizerData data{ .geom = experiment_copy.GetDiffractionGeometry(), - .latt = latt.value(), + .latt = latt, .crystal_system = symmetry.crystal_system, .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(), }; @@ -105,7 +114,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment, break; } - if (AnalyzeIndexing(output, experiment_copy, latt.value())) { + if (AnalyzeIndexing(output, experiment_copy, latt)) { float ewald_dist_cutoff = 0.001f; if (output.profile_radius) @@ -120,7 +129,7 @@ void SpotAnalyze(const DiffractionExperiment &experiment, } if (spot_finding_settings.quick_integration) { - auto res = BraggIntegrate2D(experiment_copy, image, latt.value(), + auto res = BraggIntegrate2D(experiment_copy, image, latt, prediction, ewald_dist_cutoff, output.number, symmetry.centering); constexpr size_t kMaxReflections = 10000; diff --git a/image_analysis/indexing/FFBIDXIndexer.cpp b/image_analysis/indexing/FFBIDXIndexer.cpp index d09d58c0..d3a6ed97 100644 --- a/image_analysis/indexing/FFBIDXIndexer.cpp +++ b/image_analysis/indexing/FFBIDXIndexer.cpp @@ -18,7 +18,7 @@ void FFBIDXIndexer::SetupUnitCell(const std::optional &cell) { l1.Vec2().x, l1.Vec2().y, l1.Vec2().z; } -std::vector FFBIDXIndexer::Run(const std::vector &coord, size_t nspots) { +std::vector FFBIDXIndexer::RunInternal(const std::vector &coord, size_t nspots) { std::vector ret; if (nspots > coord.size()) diff --git a/image_analysis/indexing/FFBIDXIndexer.h b/image_analysis/indexing/FFBIDXIndexer.h index 1c48a98e..cd0308a4 100644 --- a/image_analysis/indexing/FFBIDXIndexer.h +++ b/image_analysis/indexing/FFBIDXIndexer.h @@ -38,7 +38,7 @@ public: FFBIDXIndexer(const FFBIDXIndexer &i) = delete; const FFBIDXIndexer& operator=(const FFBIDXIndexer &i) = delete; - std::vector Run(const std::vector &coord, size_t nspots) override; + std::vector RunInternal(const std::vector &coord, size_t nspots) override; }; #endif //JUNGFRAUJOCH_INDEXERWRAPPER_H diff --git a/image_analysis/indexing/FFTIndexer.cpp b/image_analysis/indexing/FFTIndexer.cpp index eae67466..0465c3f5 100644 --- a/image_analysis/indexing/FFTIndexer.cpp +++ b/image_analysis/indexing/FFTIndexer.cpp @@ -159,7 +159,7 @@ std::vector FFTIndexer::FilterFFTResults() const { } -std::vector FFTIndexer::Run(const std::vector &coord, size_t nspots) { +std::vector FFTIndexer::RunInternal(const std::vector &coord, size_t nspots) { if (nspots > coord.size()) nspots = coord.size(); diff --git a/image_analysis/indexing/FFTIndexer.h b/image_analysis/indexing/FFTIndexer.h index 0f38ea73..ad907080 100644 --- a/image_analysis/indexing/FFTIndexer.h +++ b/image_analysis/indexing/FFTIndexer.h @@ -43,7 +43,7 @@ public: explicit FFTIndexer(const IndexingSettings& settings); ~FFTIndexer() override = default; - std::vector Run(const std::vector &coord, size_t nspots) override; + std::vector RunInternal(const std::vector &coord, size_t nspots) override; }; #endif //JFJOCH_FFTINDEXER_H \ No newline at end of file diff --git a/image_analysis/indexing/Indexer.cpp b/image_analysis/indexing/Indexer.cpp index c1a97390..a24f7697 100644 --- a/image_analysis/indexing/Indexer.cpp +++ b/image_analysis/indexing/Indexer.cpp @@ -12,25 +12,13 @@ void Indexer::Setup(const DiffractionExperiment& experiment) { SetupUnitCell(experiment.GetUnitCell()); } -std::optional Indexer::Run(DataMessage &message) { +IndexerResult Indexer::Run(const std::vector &coord) { + IndexerResult ret; auto start = std::chrono::high_resolution_clock::now(); - - std::vector recip; - recip.reserve(message.spots.size()); - - for (const auto &i: message.spots) { - if (index_ice_rings || !i.ice_ring) - recip.push_back(i.ReciprocalCoord(geom)); - } - - auto ret = Run(recip, recip.size()); + ret.lattice = RunInternal(coord, coord.size()); auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration duration = end - start; - message.indexing_time_s = duration.count(); - - if (ret.empty()) - return {}; - - return ret[0]; + ret.indexing_time_s = duration.count(); + return ret; } diff --git a/image_analysis/indexing/Indexer.h b/image_analysis/indexing/Indexer.h index e3207b3c..361a0223 100644 --- a/image_analysis/indexing/Indexer.h +++ b/image_analysis/indexing/Indexer.h @@ -12,6 +12,10 @@ #include "../../common/JFJochMessages.h" #include "../../common/SpotToSave.h" +struct IndexerResult { + std::vector lattice; + float indexing_time_s; +}; class Indexer { protected: @@ -24,11 +28,11 @@ protected: std::optional reference_unit_cell; virtual void SetupUnitCell(const std::optional& cell) = 0; + virtual std::vector RunInternal(const std::vector &coord, size_t nspots) = 0; public: virtual ~Indexer() = default; void Setup(const DiffractionExperiment& experiment); - virtual std::vector Run(const std::vector &coord, size_t nspots) = 0; - std::optional Run(DataMessage &message); + IndexerResult Run(const std::vector &coord); }; #endif //JFJOCH_INDEXER_H diff --git a/image_analysis/indexing/IndexerThreadPool.cpp b/image_analysis/indexing/IndexerThreadPool.cpp index dacf1bc2..ee6ff3ff 100644 --- a/image_analysis/indexing/IndexerThreadPool.cpp +++ b/image_analysis/indexing/IndexerThreadPool.cpp @@ -47,11 +47,11 @@ IndexerThreadPool::~IndexerThreadPool() { { } } -std::future > IndexerThreadPool::Run(const DiffractionExperiment &experiment, - DataMessage &message) { +std::future IndexerThreadPool::Run(const DiffractionExperiment &experiment, + const std::vector& recip) { // Create a promise/future pair - auto promise = std::make_shared > >(); - std::future > result = promise->get_future(); { + auto promise = std::make_shared >(); + std::future result = promise->get_future(); { std::unique_lock lock(m); // Don't allow enqueueing after stopping the pool @@ -60,7 +60,7 @@ std::future > IndexerThreadPool::Run(const Diffrac } // Create a task package with the data message and coordinates - taskQueue.emplace(TaskPackage{promise, &experiment, &message}); + taskQueue.emplace(TaskPackage{promise, &experiment, &recip}); } cond.notify_one(); @@ -128,7 +128,7 @@ void IndexerThreadPool::Worker(int32_t threadIndex, const NUMAHWPolicy &numa_pol } } try { - std::optional result; + IndexerResult result; auto algorithm = task.experiment->GetIndexingAlgorithm(); Indexer *indexer = nullptr; @@ -143,7 +143,7 @@ void IndexerThreadPool::Worker(int32_t threadIndex, const NUMAHWPolicy &numa_pol if (indexer) { indexer->Setup(*task.experiment); - result = indexer->Run(*task.message); + result = indexer->Run(*task.recip); } // Set the result via the promise diff --git a/image_analysis/indexing/IndexerThreadPool.h b/image_analysis/indexing/IndexerThreadPool.h index 5eae63ef..c2055aac 100644 --- a/image_analysis/indexing/IndexerThreadPool.h +++ b/image_analysis/indexing/IndexerThreadPool.h @@ -26,9 +26,9 @@ class IndexerThreadPool { std::atomic failed_start = false; struct TaskPackage { - std::shared_ptr>> promise; - const DiffractionExperiment* experiment; - DataMessage* message; + std::shared_ptr> promise; + const DiffractionExperiment *experiment; + const std::vector *recip; }; std::vector workers; @@ -44,7 +44,7 @@ public: IndexerThreadPool(const IndexingSettings& settings, const NUMAHWPolicy &numa_policy = NUMAHWPolicy()); ~IndexerThreadPool(); - std::future> Run(const DiffractionExperiment& experiment, DataMessage& message); + std::future Run(const DiffractionExperiment& experiment, const std::vector& recip); }; diff --git a/tests/FFTIndexerTest.cpp b/tests/FFTIndexerTest.cpp index 5ba4d450..b4cd67e8 100644 --- a/tests/FFTIndexerTest.cpp +++ b/tests/FFTIndexerTest.cpp @@ -36,11 +36,11 @@ TEST_CASE("FFTIndexer") { logger.Info("Spots {}", vec.size()); auto start = std::chrono::high_resolution_clock::now(); - auto result = indexer->Run(vec, vec.size()); + auto result = indexer->Run(vec); auto end = std::chrono::high_resolution_clock::now(); - REQUIRE(result.size() == 1); - auto uc_out = result[0].GetUnitCell(); + REQUIRE(result.lattice.size() == 1); + auto uc_out = result.lattice[0].GetUnitCell(); CHECK(uc_out.a == Catch::Approx(uc.a));; CHECK(uc_out.b == Catch::Approx(uc.b));; CHECK(uc_out.c == Catch::Approx(uc.c));; diff --git a/tests/IndexingUnitTest.cpp b/tests/IndexingUnitTest.cpp index c8c21804..cd8cc992 100644 --- a/tests/IndexingUnitTest.cpp +++ b/tests/IndexingUnitTest.cpp @@ -54,8 +54,8 @@ TEST_CASE("FastFeedbackIndexer","[Indexing]") { experiment.SetUnitCell(c); indexer->Setup(experiment); - auto ret = indexer->Run(recip, recip.size()); - REQUIRE(!ret.empty()); + auto ret = indexer->Run(recip); + REQUIRE(!ret.lattice.empty()); //auto uc = ret[0].GetUnitCell(); //REQUIRE(c.a == Catch::Approx(uc.a)); @@ -64,9 +64,9 @@ TEST_CASE("FastFeedbackIndexer","[Indexing]") { double err[3] = {0.0, 0.0, 0.0}; for (const auto &iter: recip) { - err[0] += round_err(ret[0].Vec0() * iter); - err[1] += round_err(ret[0].Vec1() * iter); - err[2] += round_err(ret[0].Vec2() * iter); + err[0] += round_err(ret.lattice[0].Vec0() * iter); + err[1] += round_err(ret.lattice[0].Vec1() * iter); + err[2] += round_err(ret.lattice[0].Vec2() * iter); } REQUIRE (err[0] < 0.001 * recip.size()); REQUIRE (err[1] < 0.001 * recip.size()); diff --git a/tools/jfjoch_indexing_test.cpp b/tools/jfjoch_indexing_test.cpp index 5d14e978..b975dc6a 100644 --- a/tools/jfjoch_indexing_test.cpp +++ b/tools/jfjoch_indexing_test.cpp @@ -8,6 +8,7 @@ int main(int argc, char** argv) { Logger logger("jfjoch_indexing_test"); + /* if (argc != 2) { throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Usage: ./jfjoch_indexing_test "); } @@ -21,6 +22,7 @@ int main(int argc, char** argv) { experiment.IndexingAlgorithm(IndexingAlgorithmEnum::Auto); // experiment.SetUnitCell(UnitCell(37,78,78,90,90,90)); auto indexer = CreateIndexer(experiment); + auto geom = experiment.GetDiffractionGeometry(); for (int i = 0; i < reader.GetNumberOfImages(); i++) { auto img = reader.LoadImage(i); @@ -29,7 +31,12 @@ int main(int argc, char** argv) { continue; DataMessage msg = img->ImageData(); - + std::vector recip; + recip.reserve(msg.spots.size()); + for (const auto &s: msg.spots) { + if (!s.ice_ring) + recip.emplace_back(s.ReciprocalCoord(geom)); + } auto output = indexer->Run(msg); logger.Info("Result {} {}", msg.number, msg.indexing_result.value_or(0)); @@ -42,5 +49,6 @@ int main(int argc, char** argv) { logger.Info("Lattice {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f} {:8.02f}", a[0],a[1],a[2],b[0],b[1],b[2],c[0],c[1], c[2]); } - } + } */ + } \ No newline at end of file