IndexAndRefine: Simple optimization was not good enough - using a non-linear solver leads to considerably better predictions.

image_analysis/IndexAndRefine.cpp

@@ -93,49 +93,43 @@ void IndexAndRefine::RefineGeometryIfNeeded(DataMessage &msg, IndexAndRefine::In
     if (!outcome.lattice_candidate)
         return;
 
+    XtalOptimizerData data{
+        .geom = outcome.experiment.GetDiffractionGeometry(),
+        .latt = *outcome.lattice_candidate,
+        .crystal_system = outcome.symmetry.crystal_system,
+        .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
+        .refine_beam_center = true,
+        .refine_distance_mm = false,
+        .refine_detector_angles = false,
+        .max_time = 0.04 // 40 ms is max allowed time for the operation
+    };
+
     if (experiment.IsRotationIndexing()) {
-        SimpleRotXtalOptimizerData data{
-            .geom = outcome.experiment.GetDiffractionGeometry(),
-            .latt = *outcome.lattice_candidate,
-            .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
-        };
-        if (SimpleRotXtalOptimizer(data, msg.spots)) {
-            outcome.lattice_candidate = data.latt;
-        }
-    } else {
-        XtalOptimizerData data{
-            .geom = outcome.experiment.GetDiffractionGeometry(),
-            .latt = *outcome.lattice_candidate,
-            .crystal_system = outcome.symmetry.crystal_system,
-            .min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
-            .refine_beam_center = true,
-            .refine_distance_mm = false,
-            .refine_detector_angles = false,
-            .max_time = 0.04 // 40 ms is max allowed time for the operation
-        };
+        data.refine_beam_center = false;
+        data.refine_rotation_axis = false;
+        data.refine_unit_cell = false;
+    }
 
+    if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Trigonal)
+        data.crystal_system = gemmi::CrystalSystem::Hexagonal;
 
-        if (outcome.symmetry.crystal_system == gemmi::CrystalSystem::Trigonal)
-            data.crystal_system = gemmi::CrystalSystem::Hexagonal;
+    switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) {
+        case GeomRefinementAlgorithmEnum::None:
+            break;
+        case GeomRefinementAlgorithmEnum::BeamCenter:
+            if (XtalOptimizer(data, msg.spots)) {
+                outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
+                        .BeamY_pxl(data.geom.GetBeamY_pxl());
+                outcome.beam_center_updated = true;
+            }
+            break;
+    }
 
-        switch (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm()) {
-            case GeomRefinementAlgorithmEnum::None:
-                break;
-            case GeomRefinementAlgorithmEnum::BeamCenter:
-                if (XtalOptimizer(data, msg.spots)) {
-                    outcome.experiment.BeamX_pxl(data.geom.GetBeamX_pxl())
-                                          .BeamY_pxl(data.geom.GetBeamY_pxl());
-                    outcome.beam_center_updated = true;
-                }
-                break;
-        }
+    outcome.lattice_candidate = data.latt;
 
-        outcome.lattice_candidate = data.latt;
-
-        if (outcome.beam_center_updated) {
-            msg.beam_corr_x = data.beam_corr_x;
-            msg.beam_corr_y = data.beam_corr_y;
-        }
+    if (outcome.beam_center_updated) {
+        msg.beam_corr_x = data.beam_corr_x;
+        msg.beam_corr_y = data.beam_corr_y;
     }
 }
 

image_analysis/geom_refinement/CMakeLists.txt

@@ -6,8 +6,6 @@ ADD_LIBRARY(JFJochGeomRefinement STATIC
         AssignSpotsToRings.h
         XtalOptimizer.cpp
         XtalOptimizer.h
-        SimpleRotXtalOptimizer.cpp
-        SimpleRotXtalOptimizer.h
 )
 
 TARGET_LINK_LIBRARIES(JFJochGeomRefinement Ceres::ceres Eigen3::Eigen JFJochCommon)

image_analysis/geom_refinement/SimpleRotXtalOptimizer.cpp

@@ -1,122 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
-// SPDX-License-Identifier: GPL-3.0-only
-
-#include "SimpleRotXtalOptimizer.h"
-
-#include <Eigen/Dense>
-
-static bool SimpleRotXtalOptimizerInternal(
-    SimpleRotXtalOptimizerData &data,
-    const std::vector<SpotToSave> &spots,
-    float tolerance,
-    Eigen::Matrix3d &r_total)
-{
-    Coord vec0 = data.latt.Vec0();
-    Coord vec1 = data.latt.Vec1();
-    Coord vec2 = data.latt.Vec2();
-
-    const double tol_sq = tolerance * tolerance;
-
-    std::vector<Eigen::Vector3d> obs;
-    std::vector<Eigen::Vector3d> exp;
-
-    // Collect reflections
-    for (const auto &pt : spots) {
-        if (!data.index_ice_rings && pt.ice_ring)
-            continue;
-
-        Coord recip = data.geom.DetectorToRecip(pt.x, pt.y);
-
-        double h_fp = recip * vec0;
-        double k_fp = recip * vec1;
-        double l_fp = recip * vec2;
-
-        double h = std::round(h_fp);
-        double k = std::round(k_fp);
-        double l = std::round(l_fp);
-
-        double d2 =
-            (h - h_fp)*(h - h_fp) +
-            (k - k_fp)*(k - k_fp) +
-            (l - l_fp)*(l - l_fp);
-
-        if (d2 > tol_sq)
-            continue;
-
-        obs.emplace_back(h_fp, k_fp, l_fp);
-        exp.emplace_back(h, k, l);
-    }
-
-    if (obs.size() < data.min_spots)
-        return false;
-
-    const int N = static_cast<int>(obs.size());
-
-    // Build linear system A * omega = b
-    Eigen::MatrixXd A(3*N, 3);
-    Eigen::VectorXd b(3*N);
-
-    for (int i = 0; i < N; i++) {
-        const auto &h = obs[i];
-        const auto &e = exp[i];
-
-        // Cross-product matrix of h_obs: [h]× such that [h]× · ω = h × ω
-        // But we need ω × h = -h × ω, so use -[h]×
-        A.row(3*i + 0) <<   0.0,  h.z(), -h.y();
-        A.row(3*i + 1) << -h.z(),  0.0,   h.x();
-        A.row(3*i + 2) <<  h.y(), -h.x(),  0.0;
-
-        b.segment<3>(3*i) = e - h;
-    }
-
-    // Solve least squares
-    Eigen::Vector3d omega = A.colPivHouseholderQr().solve(b);
-
-    // Optional sanity check
-    if (!omega.allFinite())
-        return false;
-
-    // Apply incremental rotation and accumulate it
-    const double angle = omega.norm();
-    Eigen::Matrix3d r_inc = Eigen::Matrix3d::Identity();
-    if (angle > 0.0)
-        r_inc = Eigen::AngleAxisd(angle, omega / angle).toRotationMatrix();
-
-    r_total = r_inc * r_total;
-
-    if (angle > 0.0) {
-        // Build the current lattice matrix (columns = a, b, c)
-        Eigen::Matrix3d L;
-        L.col(0) = Eigen::Vector3d(vec0.x, vec0.y, vec0.z);
-        L.col(1) = Eigen::Vector3d(vec1.x, vec1.y, vec1.z);
-        L.col(2) = Eigen::Vector3d(vec2.x, vec2.y, vec2.z);
-
-        // The HKL-space rotation R transforms: h_new = R * h_old
-        // This corresponds to: L_new = L_old * R^T  (real-space lattice)
-        Eigen::Matrix3d L_new = L * r_inc.transpose();
-
-        data.latt = CrystalLattice(
-            Coord(L_new(0,0), L_new(1,0), L_new(2,0)),
-            Coord(L_new(0,1), L_new(1,1), L_new(2,1)),
-            Coord(L_new(0,2), L_new(1,2), L_new(2,2))
-        );
-    }
-
-    return true;
-}
-
-bool SimpleRotXtalOptimizer(SimpleRotXtalOptimizerData &data, const std::vector<SpotToSave> &spots) {
-    Eigen::Matrix3d r_total = Eigen::Matrix3d::Identity();
-
-    if (!SimpleRotXtalOptimizerInternal(data, spots, 0.3, r_total))
-        return false;
-    SimpleRotXtalOptimizerInternal(data, spots, 0.2, r_total);
-    const bool ok = SimpleRotXtalOptimizerInternal(data, spots, 0.1, r_total);
-
-    Eigen::AngleAxisd aa(r_total);
-    data.rotation[0] = aa.axis().x();
-    data.rotation[1] = aa.axis().y();
-    data.rotation[2] = aa.axis().z();
-    data.angle = aa.angle() * 180.0 / M_PI;
-    return ok;
-}

image_analysis/geom_refinement/SimpleRotXtalOptimizer.h

@@ -1,24 +0,0 @@
-// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
-// SPDX-License-Identifier: GPL-3.0-only
-
-#ifndef JFJOCH_ROTOPTIMIZER_H
-#define JFJOCH_ROTOPTIMIZER_H
-
-#include <vector>
-#include "../common/SpotToSave.h"
-#include "../common/CrystalLattice.h"
-#include "../common/GoniometerAxis.h"
-
-struct SimpleRotXtalOptimizerData {
-    DiffractionGeometry geom;
-    CrystalLattice latt;
-    int64_t min_spots = 8;
-
-    bool index_ice_rings = true;
-    double rotation[3] = {0,0,0};
-    double angle = 0.0;
-};
-
-bool SimpleRotXtalOptimizer(SimpleRotXtalOptimizerData &data, const std::vector<SpotToSave> &spots);
-
-#endif //JFJOCH_ROTOPTIMIZER_H