# HG changeset patch
# User Eddy Xiao <bewantbe@gmail.com>
# Date 1551796113 -28800
#      Tue Mar 05 22:28:33 2019 +0800
# Branch stable
# Node ID c682e9e89b8623fc3b92bf2e136538412dd8eda0
# Parent  5d27954485ede62941b40d1852af15077ff210c7
Speed up patch() by improving coplanar point set search.(Bug #55751)

* graphics.cc (patch::properties::update_data): Collect coplanar point set by increamentally appending points and perform coplanar test by eig of rank one updated covariance matrix.

diff -r 5d27954485ed -r c682e9e89b86 libinterp/corefcn/graphics.cc
--- a/libinterp/corefcn/graphics.cc	Mon Feb 25 14:30:46 2019 +0000
+++ b/libinterp/corefcn/graphics.cc	Tue Mar 05 22:28:33 2019 +0800
@@ -9438,7 +9438,7 @@
             }
         }
     }
-
+  
   // check coplanarity for 3D-faces with more than 3 corners
   int fcmax = idx.rows ();
   if (fcmax > 3 && vert.columns () > 2 && 
@@ -9455,36 +9455,64 @@
 
               std::list<octave_idx_type> coplanar_ends;
 
+              double tol = 100 * std::numeric_limits<double>::epsilon ();
+
+              Matrix plane_pivot = Matrix (1, 3, 0.0);
+              for (octave_idx_type i = 0; i < 3; i++)
+                plane_pivot(0, i) = vert(idx(0, jj) - 1, i);
+
+              Matrix fc = Matrix (3, 3, 0.0);
+              Matrix fa = Matrix (1, 3, 0.0);
+              Matrix coor_cov = Matrix (3, 3, 0.0);
+
               octave_idx_type i_start = 1;
               octave_idx_type i_end = 2;
-              while (i_end < nc - 1)
+
+              while (i_start < nc - 1)
                 {
-                  // look for coplanar subsets
-                  for (i_end = nc-1; i_end > i_start+1; i_end--)
+                  i_end = i_start + 2;
+                  if (i_end > nc - 1)
                     {
-                      Matrix fc = Matrix (i_end - i_start + 1, 3, 0.0);
-                      for (octave_idx_type j = 0; j <= i_end-i_start; j++)
-                        for (octave_idx_type i = 0; i < 3; i++)
-                          fc(j,i) = vert(idx(j + i_start,jj)-1,i)
-                                    - vert(idx(0,jj)-1,i);
-
-                      // FIXME: Using  svd's to check for co-planarity is slow
-                      // for faces with many vertices. Is there a better way to
-                      // check this?
-
-                      // calculate rank of matrix
-                      octave::math::svd<Matrix> result
-                        (fc,
-                         octave::math::svd<Matrix>::Type::sigma_only,
-                         octave::math::svd<Matrix>::Driver::GESDD);
-                      DiagMatrix sigma = result.singular_values ();
-                      double tol = nc * sigma(0,0)
-                                   * std::numeric_limits<double>::epsilon ();
-                      if (sigma(2,2) < tol)
+                      coplanar_ends.push_back (nc - 1);
+                      break;
+                    }
+                  
+                  // Algorithm: start from 3 points, keep adding 
+                  // points until the point set is no more in
+                  // a plane. Record the coplanar point set then
+                  // advance i_start.
+
+                  // prepare 1+3 points for coplanar test
+                  // the first point is implicitly included
+                  for (octave_idx_type j = 0; j < 3; j++)
+                    for (octave_idx_type i = 0; i < 3; i++)
+                      fc(j, i) = vert(idx(j + i_start, jj) - 1, i)
+                                 - plane_pivot(i);
+                  
+                  // covariance matrix between coordinates of vertex.
+                  coor_cov = fc.transpose () * fc;
+                  
+                  while (true)
+                    {
+                      // coplanar test
+                      EIG eig (coor_cov, false, false, true);
+                      ColumnVector ev = real (eig.eigenvalues ());
+                      if (ev.min () > tol * ev.max ())
                         break;
+                      
+                      i_end ++;
+                      if (i_end > nc - 1)
+                        break;
+                      
+                      // add a point to plane
+                      for (octave_idx_type i = 0; i < 3; i++)
+                        fa(0, i) = vert(idx(i_end, jj) - 1, i)
+                                   - plane_pivot(i);
+                      coor_cov += fa.transpose () * fa;
                     }
-                  coplanar_ends.push_back (i_end);
-                  i_start = i_end;
+                  
+                  i_start = i_end - 1;
+                  coplanar_ends.push_back (i_start);
                 }
               coplanar_last_idx.push_back (coplanar_ends);
             }
