diff --git a/image_analysis/CMakeLists.txt b/image_analysis/CMakeLists.txt index 9572d0b0..913870c5 100644 --- a/image_analysis/CMakeLists.txt +++ b/image_analysis/CMakeLists.txt @@ -50,5 +50,6 @@ ADD_SUBDIRECTORY(lattice_search) ADD_SUBDIRECTORY(scale_merge) ADD_SUBDIRECTORY(image_preprocessing) ADD_SUBDIRECTORY(azint) +ADD_SUBDIRECTORY(pixel_refinement) TARGET_LINK_LIBRARIES(JFJochImageAnalysis JFJochAzIntEngine JFJochImagePreprocessing JFJochBraggPrediction JFJochBraggIntegration JFJochLatticeSearch JFJochIndexing JFJochSpotFinding JFJochCommon JFJochGeomRefinement JFJochScaleMerge gemmi) diff --git a/image_analysis/pixel_refinement/CMakeLists.txt b/image_analysis/pixel_refinement/CMakeLists.txt new file mode 100644 index 00000000..5d5e3c6c --- /dev/null +++ b/image_analysis/pixel_refinement/CMakeLists.txt @@ -0,0 +1,2 @@ +ADD_LIBRARY(JFJochPixelRefine PixelRefine.cpp PixelRefine.h) +TARGET_LINK_LIBRARIES(JFJochPixelRefine JFJochBraggPrediction JFJochCommon ceres_static) diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp new file mode 100644 index 00000000..1e547b08 --- /dev/null +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -0,0 +1,133 @@ +// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "PixelRefine.h" + +#include +#include + + + + +struct PixelRefineResidual { + ScalingPostRefResidual(int32_t h, int32_t k, int32_t l, + double x, double y, + double Itrue, + const DiffractionGeometry &geometry, + const CrystalLattice &lattice) + : + integration_center_x(r.predicted_x), + integration_center_y(r.predicted_y), + inv_lambda(SafeInv(geometry.GetWavelength_A(), 0.0)), + pixel_size(geometry.GetPixelSize_mm()), + det_dist_mm(geometry.GetDetectorDistance_mm()), + beam_x(geometry.GetBeamX_pxl()), + beam_y(geometry.GetBeamY_pxl()), + exp_h(r.h), + exp_k(r.k), + exp_l(r.l), + Astar(lattice.Astar()), + Bstar(lattice.Bstar()), + Cstar(lattice.Cstar()), + c1(std::cos(geometry.GetPoniRot1_rad())), + s1(std::sin(geometry.GetPoniRot1_rad())), + c2(std::cos(geometry.GetPoniRot2_rad())), + s2(std::sin(geometry.GetPoniRot2_rad())) { + } + + template + T CalcPartiality(const T *const R, + const T *const beam_corr, + const T *const p0) const { + // Detector coordinates in mm + const T det_x = (T(integration_center_x) - beam_x - beam_corr[0]) * T(pixel_size); + const T det_y = (T(integration_center_y) - beam_y - beam_corr[1]) * T(pixel_size); + const T det_z = T(det_dist_mm); + + // Apply Ry(rot1) first: rotate around Y + const T t1_x = T(c1) * det_x + T(s1) * det_z; + const T t1_y = det_y; + const T t1_z = T(-s1) * det_x + T(c1) * det_z; + + // Then apply Rx(-rot2): rotate around X + const T x = t1_x; + const T y = T(c2) * t1_y + T(s2) * t1_z; + const T z = - T(s2) * t1_y + T(c2) * t1_z; + + // convert to recip space + const T lab_norm = ceres::sqrt(x * x + y * y + z * z); + const T inv_norm = T(1) / lab_norm; + + T recip_obs[3]; + recip_obs[0] = x * inv_norm * inv_lambda; + recip_obs[1] = y * inv_norm * inv_lambda; + recip_obs[2] = (z * inv_norm - T(1.0)) * inv_lambda; + + const Eigen::Matrix e_obs_recip(recip_obs[0], recip_obs[1], recip_obs[2]); + + const T astar_unrot[3] = {T(Astar.x), T(Astar.y), T(Astar.z)}; + const T bstar_unrot[3] = {T(Bstar.x), T(Bstar.y), T(Bstar.z)}; + const T cstar_unrot[3] = {T(Cstar.x), T(Cstar.y), T(Cstar.z)}; + + T astar_rot[3], bstar_rot[3], cstar_rot[3]; + + ceres::AngleAxisRotatePoint(p0, astar_unrot, astar_rot); + ceres::AngleAxisRotatePoint(p0, bstar_unrot, bstar_rot); + ceres::AngleAxisRotatePoint(p0, cstar_unrot, cstar_rot); + + const Eigen::Matrix e_pred_recip(T(exp_h) * astar_rot[0] + T(exp_k) * bstar_rot[0] + T(exp_l) * cstar_rot[0], + T(exp_h) * astar_rot[1] + T(exp_k) * bstar_rot[1] + T(exp_l) * cstar_rot[1], + T(exp_h) * astar_rot[2] + T(exp_k) * bstar_rot[2] + T(exp_l) * cstar_rot[2] + ); + + // Ewald sphere centre is at -k_i = (0, 0, -inv_lambda) + // Radial direction: outward normal at g_hkl + const Eigen::Matrix S_pred( + e_pred_recip[0], + e_pred_recip[1], + e_pred_recip[2] + T(inv_lambda) // g_hkl + k_i + ); + const T S_pred_norm = S_pred.norm(); + if (S_pred_norm < T(1e-10)) + return T(0); + + const Eigen::Matrix n_radial = S_pred / S_pred_norm; + + const Eigen::Matrix delta_q = e_obs_recip - e_pred_recip; + const T eps_radial = delta_q.dot(n_radial); + const Eigen::Matrix dq_tang = delta_q - eps_radial * n_radial; + const T eps_tangential_sq = dq_tang.squaredNorm(); // guaranteed ≥ 0 + // ───────────────────────────────────────────────────────────── + + return ceres::exp(- eps_radial * eps_radial / (R[0] * R[0]) - eps_tangential_sq / (R[1] * R[1])); + } + + template + bool operator()(const T *const G, + const T *const B, + const T *const R, + const T *const beam_corr, + const T *const p0, + T *residual) const { + if (R[0] < T(1e-10) || R[1] < T(1e-10)) + return false; + + const T B_term = ceres::exp(B[0] * T(b_resolution_coeff)); + + const T partiality = CalcPartiality(R, beam_corr, p0); + residual[0] = (G[0] * partiality * B_term * T(lp) * T(Itrue) + - T(Iobs)) * T(weight); + return true; + } + + const double integration_center_x, integration_center_y; + const double inv_lambda; + const double pixel_size; + const double det_dist_mm; + const double beam_x, beam_y; + const double exp_h; + const double exp_k; + const double exp_l; + const Coord Astar, Bstar, Cstar; + const double c1,s1,c2,s2; +}; \ No newline at end of file diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h new file mode 100644 index 00000000..1610ea9c --- /dev/null +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -0,0 +1,36 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include "../bragg_prediction/BraggPrediction.h" +#include "../common/DiffractionExperiment.h" +#include "../common/AzimuthalIntegrationMapping.h" + +struct PixelRefineData { + DiffractionGeometry geom; + CrystalLattice latt; + gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic; + char centering = 'P'; + + double B_factor = 0.0; + double scale_factor = 1.0; + double R[2] = {0.001, 0.001}; + + bool refine_beam_center = false; + bool refine_unit_cell = false; +}; + +class PixelRefine { + BraggPrediction &prediction; + const AzimuthalIntegrationMapping &mapping; + const size_t xpixel, ypixel; + +public: + PixelRefine(const DiffractionExperiment &experiment, + const AzimuthalIntegrationMapping &mapping; + BraggPrediction &prediction); + + template + void Run(const T *image, PixelRefineData &data); +};