Add c++ version of xyz_to_rlp

src/dials/algorithms/indexing/boost_python/indexing_ext.cc

@@ -11,12 +11,16 @@
 #include <boost/python.hpp>
 #include <boost/python/def.hpp>
 #include <dials/algorithms/indexing/index.h>
+#include <dials/algorithms/indexing/fastindex.h>
 
 namespace dials { namespace algorithms { namespace boost_python {
 
   using namespace boost::python;
 
   void export_fft3d();
+  void export_xyz_to_rlp() {
+    def("xyz_to_rlp", &xyz_to_rlp);
+  }
 
   void export_assign_indices() {
     typedef AssignIndices w_t;
@@ -56,6 +60,7 @@ namespace dials { namespace algorithms { namespace boost_python {
   BOOST_PYTHON_MODULE(dials_algorithms_indexing_ext) {
     export_fft3d();
     export_assign_indices();
+    export_xyz_to_rlp();
     export_assign_indices_local();
   }
 

src/dials/algorithms/indexing/fastindex.h

@@ -0,0 +1,94 @@
+#include <scitbx/array_family/flex_types.h>
+#include <scitbx/vec3.h>
+#include <scitbx/mat3.h>
+#include <scitbx/math/utils.h>
+#include <cctbx/miller.h>
+#include <scitbx/constants.h>
+#include <dials/array_family/scitbx_shared_and_versa.h>
+
+class SimpleBeam {
+public:
+  double wavelength;
+  scitbx::vec3<double> s0;
+
+  SimpleBeam(double wavelength) {
+    this->wavelength = wavelength;
+    s0 = {0.0, 0.0, -1.0 / wavelength};
+  }
+};
+
+class SimpleDetector {
+public:
+  scitbx::mat3<double> d_matrix;
+  double pixel_size;  // in mm
+
+  SimpleDetector(scitbx::mat3<double> d_matrix, double pixel_size) {
+    this->d_matrix = d_matrix;
+    this->pixel_size = pixel_size;
+  }
+};
+
+class SimpleScan {
+public:
+  int image_range_start;
+  double osc_start;
+  double osc_width;
+
+  SimpleScan(int image_range_start, double osc_start, double osc_width) {
+    this->image_range_start = image_range_start;
+    this->osc_start = osc_start;
+    this->osc_width = osc_width;
+  }
+};
+
+class SimpleGonio {
+public:
+  scitbx::mat3<double> sample_rotation;
+  scitbx::vec3<double> rotation_axis;
+  scitbx::mat3<double> setting_rotation;
+  scitbx::mat3<double> sample_rotation_inverse;
+  scitbx::mat3<double> setting_rotation_inverse;
+
+  SimpleGonio(scitbx::mat3<double> sample_rotation,
+              scitbx::vec3<double> rotation_axis,
+              scitbx::mat3<double> setting_rotation) {
+    this->sample_rotation = sample_rotation;
+    this->rotation_axis = rotation_axis;
+    this->setting_rotation = setting_rotation;
+    sample_rotation_inverse = sample_rotation.inverse();
+    setting_rotation_inverse = setting_rotation.inverse();
+  }
+};
+
+scitbx::af::shared<scitbx::vec3<double>> xyz_to_rlp(
+  scitbx::af::shared<scitbx::vec3<double>> xyzobs_px,
+  scitbx::mat3<double> sample_rotation,
+  scitbx::mat3<double> detector_d_matrix) {
+  SimpleBeam beam(1.23985);
+  SimpleDetector detector(detector_d_matrix, 0.172);
+  SimpleScan scan(1, 0.0, 0.5);
+  SimpleGonio gonio(
+    sample_rotation, {1.0, 0.0, 0.0}, {1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0});
+
+  float DEG2RAD = scitbx::constants::pi / 180.0;
+
+  scitbx::af::shared<scitbx::vec3<double>> rlp(xyzobs_px.size());
+  for (int i = 0; i < xyzobs_px.size(); ++i) {
+    scitbx::vec3<double> xyzobs_i = xyzobs_px[i];
+    double x_mm = xyzobs_i[0] * detector.pixel_size;
+    double y_mm = xyzobs_i[1] * detector.pixel_size;
+    double rot_angle =
+      (((xyzobs_i[2] + 1 - scan.image_range_start) * scan.osc_width) + scan.osc_start)
+      * DEG2RAD;
+    scitbx::vec3<double> m = {x_mm, y_mm, 1.0};
+    scitbx::vec3<double> s1_i = detector.d_matrix * m;
+    double norm =
+      std::pow(std::pow(s1_i[0], 2) + std::pow(s1_i[1], 2) + std::pow(s1_i[2], 2), 0.5);
+    scitbx::vec3<double> s1_this = (s1_i / norm) * (1.0 / beam.wavelength);
+    scitbx::vec3<double> S = gonio.setting_rotation_inverse * (s1_this - beam.s0);
+    scitbx::vec3<double> rlp_this =
+      S.rotate_around_origin(gonio.rotation_axis, -1.0 * rot_angle);
+    rlp[i] = gonio.sample_rotation_inverse * rlp_this;
+  }
+  return rlp;
+}

src/dials/command_line/autoindex.py

@@ -15,6 +15,7 @@
     group_vectors,
     is_approximate_integer_multiple,
 )
+from dials_algorithms_indexing_ext import xyz_to_rlp as xyz_to_rlp_cpp
 
 
 def _find_peaks(
@@ -169,6 +170,7 @@ def sites_to_vecs(sites, volumes, fft_cell, min_cell=3, max_cell=92.3):
 
 r = flex.reflection_table.from_file("strong_1_60.refl")
 xyzobs_px = r["xyzobs.px.value"]
+print(xyzobs_px)
 st = time.time()
 from dxtbx.serialize import load
 
@@ -190,7 +192,9 @@ def xyz_to_rlp(xyzobs_px, expt):
     setting_rotation = matrix.sqr(expt.goniometer.get_setting_rotation())
     rotation_axis = expt.goniometer.get_rotation_axis_datum()
     sample_rotation = matrix.sqr(expt.goniometer.get_fixed_rotation())
-
+    print(setting_rotation)
+    print(rotation_axis)
+    print(type(sample_rotation))
     # centroid_px_to_mm_panel
     x, y, z = xyzobs_px.parts()
     # fixme this is simple px to mm, but really should be ParallaxCorrectedPxMmStrategy.to_millimeter
@@ -236,8 +240,41 @@ def xyz_to_rlp(xyzobs_px, expt):
     return rlp
 
 
+import time
+
+st = time.time()
 rlp = xyz_to_rlp(xyzobs_px, expt)
-print(len(rlp))
+end = time.time()
+print(end - st)
+
+i_panel = 0
+f = expt.detector[i_panel].get_fast_axis()
+s = expt.detector[i_panel].get_slow_axis()
+n = expt.detector[i_panel].get_normal()
+origin = expt.detector[i_panel].get_origin()
+d_ = matrix.sqr(
+    (
+        f[0],
+        s[0],
+        n[0] + origin[0],
+        f[1],
+        s[1],
+        n[1] + origin[1],
+        f[2],
+        s[2],
+        n[2] + origin[2],
+    )
+)
+st = time.time()
+rlp2 = xyz_to_rlp_cpp(xyzobs_px, matrix.sqr(expt.goniometer.get_fixed_rotation()), d_)
+end = time.time()
+print(end - st)
+print(rlp[100])
+print(rlp2[100])
+for r1, r2 in zip(rlp, rlp2):
+    for i in range(0, 3):
+        assert abs(r1[i] - r2[i]) < 1e-8, f"{r1[i]}, {r2[i]}"
+assert 0
 
 # then find_basis_vectors - fft3d
 candidate_basis_vectors, used_in_indexing = do_fft3d(rlp, d_min=1.8)