From d1ed975cf82f8ad31b35a46b0c2dbfc0b2d277e4 Mon Sep 17 00:00:00 2001
From: Stephan Schulz <stephan.schulz-x2q@ruhr-uni-bochum.de>
Date: Fri, 10 Sep 2021 13:54:14 +0200
Subject: [PATCH] Fix point distance calculation

---
 CHANGELOG.rst                   |  7 ++++++
 include/ALL_Point.hpp           | 41 +++++++++++++++------------------
 tests/unit/point_arithmetic.cpp | 15 ++++++++++++
 3 files changed, 41 insertions(+), 22 deletions(-)

diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 8972a24..fd07d4c 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -11,6 +11,13 @@ Version 1.0.0
 Version 0.9
 -----------
 
+Version 0.9.3
+*************
+
+- Bug: ALL::Point.dist_plane(..) had a type bug. This manifestes with a cuda
+  compiler even without explicit usage of the function. Fixed and added
+  corresponding test.
+
 Version 0.9.2
 *************
 
diff --git a/include/ALL_Point.hpp b/include/ALL_Point.hpp
index 760cbdb..d6db689 100644
--- a/include/ALL_Point.hpp
+++ b/include/ALL_Point.hpp
@@ -220,31 +220,17 @@ public:
       throw PointDimensionMissmatchException(__FILE__, __func__, __LINE__);
 
     // vectors spanning the plane from vertex 'a'
-    T vb[3];
-    T vc[3];
-
-    for (int i = 0; i < 3; ++i) {
-      vb[i] = B[i] - A[i];
-      vc[i] = C[i] - A[i];
-    }
+    Point<T> vb = B - A;
+    Point<T> vc = C - A;
 
     // normal vector of plane
-    T n[3];
-
-    n[0] = vb[1] * vc[2] - vb[2] * vc[1];
-    n[1] = vb[2] * vc[0] - vb[0] * vc[2];
-    n[2] = vb[0] * vc[1] - vb[1] * vc[0];
-
-    // length of normal vector
-    T dn = n.norm();
+    Point<T> n = vb.cross(vc);
+    n = n/n.norm();
 
-    // distance from origin
-    T d = n * A;
-
-    // compute distance of test vertex to plane
-    return std::abs((n[0] * coordinates.at(0) + n[1] * coordinates.at(1) +
-                     n[2] * coordinates.at(2) - d) /
-                    dn);
+    // return r.n - A.n = (r-A).n as distance from plane to r
+    return std::abs( (coordinates.at(0) - A[0]) * n[0] +
+		     (coordinates.at(1) - A[1]) * n[1] +
+		     (coordinates.at(2) - A[2]) * n[2] );
   }
 
   /// operator for the addition of two point objects
@@ -301,6 +287,17 @@ public:
     return result;
   }
 
+  /// operator to scale the local point object by a provided factor
+  /// @param rhs scaling factor
+  /// @return scaled point object
+  Point<T> operator/(const T &rhs) const {
+    Point<T> result(dimension);
+    for (int d = 0; d < dimension; ++d) {
+      result[d] = coordinates.at(d) / rhs;
+    }
+    return result;
+  }
+
   /// operator to compute the cross product between two point objects
   /// @param rhs point object to compute the cross product with
   /// @return cross product between the two point objects
diff --git a/tests/unit/point_arithmetic.cpp b/tests/unit/point_arithmetic.cpp
index ba528c5..2df27ba 100644
--- a/tests/unit/point_arithmetic.cpp
+++ b/tests/unit/point_arithmetic.cpp
@@ -172,5 +172,20 @@ BOOST_AUTO_TEST_CASE(point_scale_long_double) {
     BOOST_CHECK_CLOSE(C2[2],4.4l,1e-14);
 }    
 
+BOOST_AUTO_TEST_CASE(point_dist_plane) {
+    std::vector<double> A_data = {1.,0.,0.};
+    std::vector<double> B_data = {0.,1.,0.};
+    std::vector<double> C_data = {0.,0.,1.};
+    std::vector<double> p_data = {0.5,0.5,0.5};
+    ALL::Point<double> A(A_data);
+    ALL::Point<double> B(B_data);
+    ALL::Point<double> C(C_data);
+    ALL::Point<double> p(p_data);
+
+    double dist = p.dist_plane(A,B,C);
+
+    BOOST_CHECK_CLOSE(dist, 0.28867513459481292, 1e-14);
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 
-- 
GitLab