Fix bug in half_space_box (#2269) [ci skip]

* test

* fix 1384 and add test

include/igl/copyleft/cgal/half_space_box.cpp

@@ -13,30 +13,25 @@ IGL_INLINE void igl::copyleft::cgal::half_space_box(
   typedef CGAL::Point_3<CGAL::Epeck> Point;
   typedef CGAL::Vector_3<CGAL::Epeck> Vector;
   typedef CGAL::Epeck::FT EScalar;
-  Eigen::Matrix<typename DerivedV::Scalar,1,3> avg(0,0,0);
-  for(int v = 0;v<V.rows();v++) for(int c = 0;c<V.cols();c++) avg(c) += V(v,c);
-  avg /= V.rows();
-
-  Point o3(avg(0),avg(1),avg(2));
+  Point o3(V(0,0),V(0,1),V(0,2));
   Point o2 = P.projection(o3);
-  Vector u;
-  EScalar max_sqrd = -1;
-  for(int v = 0;v<V.rows();v++)
+    //https://math.stackexchange.com/a/4112622/35376
+  const auto copysign = [](const EScalar & a, const EScalar & b)->EScalar
   {
-    Vector v2 = P.projection(Point(V(v,0),V(v,1),V(v,2))) - o2;
-    const EScalar sqrd = v2.squared_length();
-    if(max_sqrd<0 || sqrd < max_sqrd)
-    {
-      u = v2;
-      max_sqrd = sqrd;
-    }
-  }
-  // L1 bbd 
+    return (b>=0 ? abs(a) : -abs(a));
+  };
+  Vector n = P.orthogonal_vector();
+  Vector u(
+    copysign(n.z(),n.x()),
+    copysign(n.z(),n.y()),
+   -copysign(n.x(),n.z()) - copysign(n.y(),n.z()));
+  
+  // L1 of bounding box diagonal. ‖x‖₂ ≤ ‖x‖₁
   const EScalar bbd = 
     (EScalar(V.col(0).maxCoeff())- EScalar(V.col(0).minCoeff())) + 
     (EScalar(V.col(1).maxCoeff())- EScalar(V.col(1).minCoeff())) + 
     (EScalar(V.col(2).maxCoeff())- EScalar(V.col(2).minCoeff()));
-  Vector n = P.orthogonal_vector();
+  const EScalar bbd_sqr = bbd*bbd;
   // now we have a center o2 and a vector u to the farthest point on the plane 
   std::vector<Point> vBV;vBV.reserve(8);
   Vector v = CGAL::cross_product(u,n);
@@ -49,9 +44,9 @@ IGL_INLINE void igl::copyleft::cgal::half_space_box(
       x = 2.*x;
     }
   };
-  longer_than(bbd*bbd,u);
-  longer_than(bbd*bbd,v);
-  longer_than(bbd*bbd,n);
+  longer_than(bbd_sqr,u);
+  longer_than(bbd_sqr,v);
+  longer_than(bbd_sqr,n);
   vBV.emplace_back( o2 + u + v);
   vBV.emplace_back( o2 - u + v);
   vBV.emplace_back( o2 - u - v);
@@ -109,8 +104,11 @@ IGL_INLINE void igl::copyleft::cgal::half_space_box(
   Plane P(equ(0),equ(1),equ(2),equ(3));
   return half_space_box(P,V,BV,BF);
 }
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::copyleft::cgal::half_space_box<Eigen::Matrix<double, 1, 4, 1, 1, 4>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 4, 1, 1, 4>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::Matrix<CGAL::Epeck::FT, 8, 3, 0, 8, 3>&, Eigen::Matrix<int, 12, 3, 0, 12, 3>&);
 template void igl::copyleft::cgal::half_space_box<Eigen::Matrix<double, -1, 3, 1, -1, 3> >(CGAL::Plane_3<CGAL::Epeck> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::Matrix<CGAL::Epeck::FT, 8, 3, 0, 8, 3>&, Eigen::Matrix<int, 12, 3, 0, 12, 3>&);
 // generated by autoexplicit.sh
 template void igl::copyleft::cgal::half_space_box<Eigen::Matrix<float, -1, 3, 1, -1, 3> >(CGAL::Plane_3<CGAL::Epeck> const&, Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::Matrix<CGAL::Epeck::FT, 8, 3, 0, 8, 3>&, Eigen::Matrix<int, 12, 3, 0, 12, 3>&);

tests/include/igl/copyleft/cgal/half_space_box.cpp

@@ -0,0 +1,20 @@
+#include <igl/copyleft/cgal/half_space_box.h>
+#include <test_common.h>
+
+TEST_CASE("half_space_box: simple", "[igl]")
+{
+  Eigen::MatrixXd V(2,3);
+  V<<0,0,0,
+    1,1,1;
+  Eigen::RowVector4d equ;
+  Eigen::Matrix<CGAL::Epeck::FT,8,3> BV;
+  Eigen::Matrix<int,12,3> BF;
+
+  equ << 1,-1,0,0;
+  igl::copyleft::cgal::half_space_box(equ,V,BV,BF);
+  REQUIRE((BV.colwise().maxCoeff() - BV.colwise().minCoeff()).squaredNorm() > 0);
+
+  equ << 1,1,1,0;
+  igl::copyleft::cgal::half_space_box(equ,V,BV,BF);
+  REQUIRE((BV.colwise().maxCoeff() - BV.colwise().minCoeff()).squaredNorm() > 0);
+}