diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 76f972d8..6dfbb683 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -276,7 +276,7 @@ std::optional IndexAndRefine::Finalize() { return {}; } -std::optional IndexAndRefine::ScaleAllImages(const ScaleMergeOptions &opts) { +std::optional IndexAndRefine::ScaleAllImages(const ScaleMergeOptions &opts) { size_t nrefl = 0; for (const auto &i: reflections) nrefl += i.size(); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 2a5e5b86..4ef24b13 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -62,7 +62,7 @@ public: void ProcessImage(DataMessage &msg, const SpotFindingSettings &settings, const CompressedImage &image, BraggPrediction &prediction); - std::optional ScaleAllImages(const ScaleMergeOptions &opts = {}); + std::optional ScaleAllImages(const ScaleMergeOptions &opts = {}); /// Run scale-and-merge on accumulated reflections to refine per-image /// mosaicity (and optionally B-factors / scale factors). diff --git a/image_analysis/scale_merge/CMakeLists.txt b/image_analysis/scale_merge/CMakeLists.txt index 3f06312e..2bde4a67 100644 --- a/image_analysis/scale_merge/CMakeLists.txt +++ b/image_analysis/scale_merge/CMakeLists.txt @@ -6,5 +6,6 @@ ADD_LIBRARY(JFJochScaleMerge ScaleAll.cpp ScaleAll.h FrenchWilson.cpp FrenchWils ScaleOnTheFly.cpp ScaleOnTheFly.h HKLKey.cpp - HKLKey.h) + HKLKey.h + ScalingResult.h) TARGET_LINK_LIBRARIES(JFJochScaleMerge Ceres::ceres Eigen3::Eigen JFJochCommon) \ No newline at end of file diff --git a/image_analysis/scale_merge/ScaleAll.cpp b/image_analysis/scale_merge/ScaleAll.cpp index f83514dd..310392c1 100644 --- a/image_analysis/scale_merge/ScaleAll.cpp +++ b/image_analysis/scale_merge/ScaleAll.cpp @@ -410,7 +410,7 @@ namespace { } } -ScaleResult ScaleAll(std::vector > &observations, +ScalingResult ScaleAll(std::vector > &observations, const ScaleMergeOptions &opt) { if (opt.image_cluster <= 0) throw std::invalid_argument("image_cluster must be positive"); @@ -442,7 +442,7 @@ ScaleResult ScaleAll(std::vector > &observations, CalcCorrections(observations, opt, G, mosaicity, R_sq); - ScaleResult out{}; + ScalingResult out{}; out.image_scale_g.resize(observations.size(), NAN); out.mosaicity_deg.resize(observations.size(), NAN); diff --git a/image_analysis/scale_merge/ScaleAll.h b/image_analysis/scale_merge/ScaleAll.h index f06aaa26..45766914 100644 --- a/image_analysis/scale_merge/ScaleAll.h +++ b/image_analysis/scale_merge/ScaleAll.h @@ -4,13 +4,7 @@ #pragma once #include "Merge.h" +#include "ScalingResult.h" -struct ScaleResult { - std::vector image_scale_g; - std::vector mosaicity_deg; - std::vector image_bfactor_Ang2; -}; - -ScaleResult ScaleAll(std::vector>& observations, - const ScaleMergeOptions& opt = {}); \ No newline at end of file +ScalingResult ScaleAll(std::vector>& observations, const ScaleMergeOptions& opt = {}); \ No newline at end of file diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index 5bf75560..f225d6cc 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -3,6 +3,52 @@ #include "ScaleOnTheFly.h" +#include + +double SafeInv(double x, double fallback) { + if (!std::isfinite(x) || x == 0.0) + return fallback; + return 1.0 / x; +} + +struct IntensityRotationResidual { + IntensityRotationResidual(const Reflection &r, double Itrue, double sigma) + : Iobs(r.I), + Itrue(Itrue), + weight(SafeInv(sigma, 1.0)), + delta_phi_deg(r.delta_phi_deg), + lp(SafeInv(r.rlp, 1.0)), + c1(r.zeta / std::sqrt(2.0)), + b_resolution_coeff(SafeInv(-r.d * r.d / 4.0, 0.0)) { + } + + template + bool operator()(const T *const G, + const T *const B, + const T *const mosaicity, + const T *const wedge, + T *residual) const { + const T partiality = T(1.0); + if (mosaicity > 0) { + const T half_wedge = wedge[0] / T(2.0); + const T arg_plus = T(delta_phi_deg + half_wedge) * T(c1) / mosaicity[0]; + const T arg_minus = T(delta_phi_deg - half_wedge) * T(c1) / mosaicity[0]; + partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0); + } + const T B_term = ceres::exp(T(B[0]) * b_resolution_coeff); + residual[0] = (G[0] * partiality * T(lp) * Itrue - T(Iobs)) * T(weight); + return true; + } + + double Iobs; + double Itrue; + double weight; + double delta_phi_deg; + double lp; + double c1; + double b_resolution_coeff; +}; + ScaleOnTheFly::ScaleOnTheFly(const std::vector &ref, std::optional sg, bool merge_friedel) : sg(sg), merge_friedel(merge_friedel) { @@ -11,3 +57,7 @@ ScaleOnTheFly::ScaleOnTheFly(const std::vector &ref, reference_data[key] = r.I; } } + +std::optional ScaleOnTheFly::Scale(std::vector &reflections) { + +} diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h index 21b4a2ba..798698e0 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.h +++ b/image_analysis/scale_merge/ScaleOnTheFly.h @@ -5,6 +5,7 @@ #include "HKLKey.h" #include "Merge.h" +#include "ScalingResult.h" #include @@ -14,6 +15,6 @@ class ScaleOnTheFly { std::map reference_data; public: ScaleOnTheFly(const std::vector &ref, std::optional sg, bool merge_friedel = true); - void Scale(std::vector &reflections); + std::optional Scale(std::vector &reflections); }; diff --git a/image_analysis/scale_merge/ScalingResult.h b/image_analysis/scale_merge/ScalingResult.h new file mode 100644 index 00000000..e4701814 --- /dev/null +++ b/image_analysis/scale_merge/ScalingResult.h @@ -0,0 +1,12 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include + +struct ScalingResult { + std::vector image_scale_g; + std::vector mosaicity_deg; + std::vector image_bfactor_Ang2; +}; \ No newline at end of file