Merge pull request #5412 from GilesBathgate/performance-calculate-normals

Calculate normals for facets more efficiently

Nef_3/include/CGAL/normal_vector_newell_3.h

@@ -23,6 +23,8 @@
 #undef CGAL_NEF_DEBUG
 #define CGAL_NEF_DEBUG 79
 #include <CGAL/Nef_2/debug.h>
+#include <CGAL/representation_tags.h>
+#include <CGAL/Kernel_traits.h>
 
 namespace CGAL {
 
@@ -30,32 +32,72 @@ namespace internal_nef
 {
 template <class Handle, class Vector>
 CGAL_MEDIUM_INLINE
-void newell_single_step_3( const Handle& p, const Handle& q, Vector& n )
+void newell_single_step_3( const Handle& p, const Handle& q, Vector& n, const Homogeneous_tag&)
 {
+  typedef typename Kernel_traits<Vector>::Kernel::RT RT;
+  const RT& phw = p.hw();
+  const RT& qhw = q.hw();
+  const RT& nhw = n.hw();
+  const RT& sq = square(phw) * square(qhw);
+  const RT& phyqhw = p.hy() * qhw;
+  const RT& qhyphw = q.hy() * phw;
+  const RT& phxqhw = p.hx() * qhw;
+  const RT& qhxphw = q.hx() * phw;
+  const RT& phzqhw = p.hz() * qhw;
+  const RT& qhzphw = q.hz() * phw;
+
   n = Vector(
     n.hx()
-        * p.hw() * p.hw()
-        * q.hw() * q.hw()
-      +   n.hw()
-        * ( p.hy() * q.hw() - q.hy() * p.hw())
-        * ( p.hz() * q.hw() + q.hz() * p.hw()),
+        * sq
+      + nhw
+        * ( phyqhw - qhyphw )
+        * ( phzqhw + qhzphw ),
     n.hy()
-        * p.hw() * p.hw()
-        * q.hw() * q.hw()
-      +   n.hw()
-        * ( p.hz() * q.hw() - q.hz() * p.hw())
-        * ( p.hx() * q.hw() + q.hx() * p.hw()),
+        * sq
+      + nhw
+        * ( phzqhw - qhzphw )
+        * ( phxqhw + qhxphw ),
     n.hz()
-        * p.hw() * p.hw()
-        * q.hw() * q.hw()
-      +   n.hw()
-        * ( p.hx() * q.hw() - q.hx() * p.hw())
-        * ( p.hy() * q.hw() + q.hy() * p.hw()),
+        * sq
+      + nhw
+        * ( phxqhw - qhxphw )
+        * ( phyqhw + qhyphw ),
     n.hw()
-        * p.hw() * p.hw()
-        * q.hw() * q.hw()
+        * sq
     );
 }
+
+template <class Handle, class Vector>
+CGAL_MEDIUM_INLINE
+void newell_single_step_3( const Handle& p, const Handle& q, Vector& n, const Cartesian_tag&)
+{
+  typedef typename Kernel_traits<Vector>::Kernel::FT FT;
+  const FT& py = p.y();
+  const FT& qy = q.y();
+  const FT& px = p.x();
+  const FT& qx = q.x();
+  const FT& pz = p.z();
+  const FT& qz = q.z();
+
+  n = Vector(
+    n.x()
+        + ( py - qy )
+        * ( pz + qz ),
+    n.y()
+        + ( pz - qz )
+        * ( px + qx ),
+    n.z()
+        + ( px - qx )
+        * ( py + qy )
+    );
+}
+
+template <class IC>
+bool is_triangle_3( const IC& first )
+{
+    return std::next(first,3) == first;
+}
+
 }
 
 template <class IC, class Vector>
@@ -70,6 +112,14 @@ void normal_vector_newell_3( IC first, IC last, Vector& n )
     // three.
 {
     CGAL_assertion( !CGAL::is_empty_range( first, last));
+
+    if(internal_nef::is_triangle_3(first)) {
+      n = orthogonal_vector(*first,*(std::next(first)),*(std::next(first,2)));
+      return;
+    }
+
+
+    typedef typename Kernel_traits<Vector>::Kernel R;
     // Compute facet normals via the Newell-method as described in
     // Filippo Tampieri: Newell's Method for Computing the Plane
     // Equation of a Polygon. Graphics Gems III, David Kirk,
@@ -80,11 +130,11 @@ void normal_vector_newell_3( IC first, IC last, Vector& n )
     IC prev = first;
     ++first;
     while( first != last) {
-        internal_nef::newell_single_step_3( *prev, *first, n);
+        internal_nef::newell_single_step_3( *prev, *first, n, typename R::Kernel_tag());
         prev = first;
         ++first;
     }
-    internal_nef::newell_single_step_3( *prev, *start_point, n);
+    internal_nef::newell_single_step_3( *prev, *start_point, n, typename R::Kernel_tag());
     CGAL_NEF_TRACEN("newell normal vector "<<n);
 }
 
@@ -92,6 +142,15 @@ template <class IC, class Vector, class VertexPointMap>
 void normal_vector_newell_3( IC first, IC last, VertexPointMap vpm, Vector& n )
 {
     CGAL_assertion( !CGAL::is_empty_range( first, last));
+
+    if(internal_nef::is_triangle_3(first)) {
+
+      n = orthogonal_vector(get(vpm,*first),get(vpm,*(std::next(first))),get(vpm,*(std::next(first,2))));
+      return;
+    }
+
+
+    typedef typename Kernel_traits<Vector>::Kernel R;
     // Compute facet normals via the Newell-method as described in
     // Filippo Tampieri: Newell's Method for Computing the Plane
     // Equation of a Polygon. Graphics Gems III, David Kirk,
@@ -102,11 +161,11 @@ void normal_vector_newell_3( IC first, IC last, VertexPointMap vpm, Vector& n )
     IC prev = first;
     ++first;
     while( first != last) {
-        internal_nef::newell_single_step_3( get(vpm,*prev), get(vpm,*first), n);
+        internal_nef::newell_single_step_3( get(vpm,*prev), get(vpm,*first), n, typename R::Kernel_tag());
         prev = first;
         ++first;
     }
-    internal_nef::newell_single_step_3( get(vpm,*prev), get(vpm,*start_point), n);
+    internal_nef::newell_single_step_3( get(vpm,*prev), get(vpm,*start_point), n, typename R::Kernel_tag());
     CGAL_NEF_TRACEN("newell normal vector "<<n);
 }