From 83d7a89f22585a2da956b411bcd1817c80908f49 Mon Sep 17 00:00:00 2001
From: Giles Bathgate <giles.bathgate@gmail.com>
Date: Fri, 29 Jan 2021 17:32:05 +0000
Subject: [PATCH] Calculate normals for facets more efficiently

---
 Nef_3/include/CGAL/normal_vector_newell_3.h | 100 +++++++++++++++-----
 1 file changed, 78 insertions(+), 22 deletions(-)

diff --git a/Nef_3/include/CGAL/normal_vector_newell_3.h b/Nef_3/include/CGAL/normal_vector_newell_3.h
index 282f8e5b5cb..5e34f350d33 100644
--- a/Nef_3/include/CGAL/normal_vector_newell_3.h
+++ b/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,74 @@ 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 = phw * phw * qhw * 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& phy = p.hy();
+  const FT& qhy = q.hy();
+  const FT& phx = p.hx();
+  const FT& qhx = q.hx();
+  const FT& phz = p.hz();
+  const FT& qhz = q.hz();
+
+  n = Vector(
+    n.hx()
+        + ( phy - qhy )
+        * ( phz + qhz ),
+    n.hy()
+        + ( phz - qhz )
+        * ( phx + qhx ),
+    n.hz()
+        + ( phx - qhx )
+        * ( phy + qhy )
+    );
+}
+
+template <class IC>
+bool is_triangle_3( const IC& first )
+{
+    IC last(first);
+    ++last; ++last; ++last;
+    return first==last;
+}
+
 }
 
 template <class IC, class Vector>
@@ -70,6 +114,12 @@ 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,*(++first),*(++first));
+      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,12 @@ 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,*(++first)),get(vpm,*(++first)));
+      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 +158,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);
 }