Use Zack's Eigen implementation of inverse distortion (#699)

* Pull in Zack's implementation of inverse distortion

* sqrt fix for mac install

* Use math defines

* replace constexpr with const to fix stuff

* Including cmath so use_math_defines does something

* Move up the math defines because eigen probably include cmath already.

* Accidentally included xsimd twice

Author: kevindlewis23

hexrd/transforms/cpp_sublibrary/src/inverse_distortion.cpp

@@ -1,39 +1,31 @@
-#define FORCE_IMPORT_ARRAY
-#define XTENSOR_USE_XSIMD
-
+#define _USE_MATH_DEFINES 
 #include <pybind11/pybind11.h>
-#include <xtensor/xmath.hpp>
-#include <xtensor/xarray.hpp>
-#include <xtensor-python/pyarray.hpp>
-#include <xtensor/xview.hpp>
+#include <pybind11/eigen.h>
+#include <Eigen/Core>
+#include <xsimd/xsimd.hpp>
+#include <cmath>
 
 namespace py = pybind11;
-constexpr double FOUR_THIRDS_PI = 4.1887902;
-constexpr double N_THREE_HALVES_SQRT_3 = -2.59807621;
-constexpr double TWO_OVER_SQRT_THREE = 1.154700538;
-
-xt::pyarray<double> ge_41rt_inverse_distortion(const xt::pyarray<double>& inputs, const double rhoMax, const xt::pyarray<double>& params) {
-    auto radii = xt::sqrt(xt::sum(xt::square(inputs), {1}));
-
-    if (xt::amax(radii)() < 1e-10) {
-        return xt::zeros<double>({inputs.shape()[0], inputs.shape()[1]});
-    }
-
-    auto inverted_radii = xt::pow(radii, -1);
-    xt::pyarray<double> cosines = xt::view(inputs, xt::all(), 0) * inverted_radii;
-    auto cosine_double_angles = 2 * xt::square(cosines) - 1;
-    auto cosine_quadruple_angles = 2 * xt::square(cosine_double_angles) - 1;
-    auto sqrt_p_is = rhoMax / xt::sqrt(-params[0] * cosine_double_angles - params[1] * cosine_quadruple_angles - params[2]);
-    auto solutions = TWO_OVER_SQRT_THREE * sqrt_p_is * xt::cos(xt::acos(N_THREE_HALVES_SQRT_3 * radii / sqrt_p_is) / 3 + FOUR_THIRDS_PI);
-    xt::pyarray<double> results = xt::view(solutions, xt::all(), xt::newaxis()) * inputs * xt::view(inverted_radii, xt::all(), xt::newaxis());
-
-    return results;
+const double FOUR_THIRDS_PI = M_PI * 4.0 / 3.0;
+const double N_THREE_HALVES_SQRT_3 = -3.0 / 2.0 * std::sqrt(3.0);
+const double TWO_OVER_SQRT_THREE = 2.0 / std::sqrt(3.0);
+
+Eigen::ArrayXXd ge_41rt_inverse_distortion(const Eigen::ArrayXXd& inputs, const double rhoMax, const Eigen::ArrayXd& params) {
+  Eigen::ArrayXd radii = inputs.matrix().rowwise().norm();
+  Eigen::ArrayXd inverted_radii = radii.cwiseInverse();
+  Eigen::ArrayXd cosines = inputs.col(0) * inverted_radii;
+  Eigen::ArrayXd cosine_double_angles = 2*cosines.square() - 1;
+  Eigen::ArrayXd cosine_quadruple_angles = 2*cosine_double_angles.square() - 1;
+  Eigen::ArrayXd sqrt_p_is = rhoMax / (-params[0]*cosine_double_angles - params[1]*cosine_quadruple_angles - params[2]).sqrt();
+  Eigen::ArrayXd solutions = TWO_OVER_SQRT_THREE*sqrt_p_is*(acos(N_THREE_HALVES_SQRT_3*radii/sqrt_p_is)/3 + FOUR_THIRDS_PI).cos();
+  Eigen::ArrayXXd results = solutions.rowwise().replicate(inputs.cols()).array() * inputs * inverted_radii.rowwise().replicate(inputs.cols()).array();
+
+  return results;
 }
 
 PYBIND11_MODULE(inverse_distortion, m)
 {
-    xt::import_numpy();
     m.doc() = "HEXRD inverse distribution plugin";
 
     m.def("ge_41rt_inverse_distortion", &ge_41rt_inverse_distortion, "Inverse distortion for ge_41rt");
-}
+}

setup.py

@@ -123,9 +123,7 @@ def get_cpp_extensions():
     # Define include directories
     include_dirs = [
         get_include_path('xsimd'),
-        get_include_path('xtensor'),
-        get_include_path('xtensor-python'),
-        get_include_path('xtl'),
+        get_include_path('eigen3'),
         get_pybind11_include_path(),
         numpy.get_include(),
     ]