# HG changeset patch
# User mmuetzel <markus.muetzel@gmx.de>
# Date 1458829942 -3600
#      Thu Mar 24 15:32:22 2016 +0100
# Node ID 408c51b83ae82299c9fddc096da8da7ed4743e62
# Parent  ffbd3e86e0be6a2ff8808d124290d8566313fbf1
Automatically calculate vertexnormals for patch objects.

* graphics.in.h: Add public method patch::properties::update_normals
for calculation of vertexnormals. Add and adapt update methods to
trigger re-calculation of vertexnormals when data has changed. Change
constraints for vertexnormals. Add public method
patch::reset_default_properties. Add protected method patch::initialize.
* graphics.cc: Implement method patch::properties::update_normals for
calculation of vertexnormals. Move inline function cross_product up to
be able to use it in update_normals. Use inherited methods
patch::reset_default_properties and patch::initialize to trigger
calculation of vertexnormals when patch was reset or created. Adapt test
for patch to add checks for vertexnormals.
* patch.m: Adapt test for patch to add checks for vertexnormals. The
test currently fails.

diff -r ffbd3e86e0be -r 408c51b83ae8 libinterp/corefcn/graphics.cc
--- a/libinterp/corefcn/graphics.cc	Tue Mar 22 23:14:20 2016 -0700
+++ b/libinterp/corefcn/graphics.cc	Thu Mar 24 15:32:22 2016 +0100
@@ -8257,6 +8257,8 @@
         }
     }
 
+  // Update vertex normals
+  update_normals ();
 
   unwind_protect frame;
   frame.protect_var (updating_patch_data);
@@ -8268,14 +8270,6 @@
   set_cdata (cd);
 }
 
-// ---------------------------------------------------------------------
-
-octave_value
-surface::properties::get_color_data (void) const
-{
-  return convert_cdata (*this, get_cdata (), cdatamapping_is ("scaled"), 3);
-}
-
 inline void
 cross_product (double x1, double y1, double z1,
                double x2, double y2, double z2,
@@ -8287,6 +8281,182 @@
 }
 
 void
+patch::properties::update_normals (void)
+{
+  if (normalmode_is ("auto"))
+    {
+      Matrix v = get_vertices ().matrix_value ();
+      Matrix f = get_faces ().matrix_value ();
+
+      octave_idx_type num_v = v.rows (); // number of vertices
+      bool is_3d = (v.columns () == 3); // 2d or 3d patches
+      octave_idx_type num_f = f.rows (); // number of faces
+      octave_idx_type max_nc = f.columns (); // maximum number of polygon corners
+
+      // In which cases can we skip updating the normals?
+      if (max_nc < 3)
+        return;
+
+      // First step: Calculate the normals at the corners of all faces and
+      // assign to the respective vertices
+      std::vector<RowVector> vec_vn [num_v]; // list of normals for vertices
+      octave_idx_type i1, i2, i3;
+      octave_idx_type j1, j2, j3;
+      for (octave_idx_type i = 0; i < num_f; i++)
+        {
+          // check whether the face is a triangle
+          bool is_triangle = (max_nc == 3);
+          if (! is_triangle)
+            if (xisnan (f(i,3))) // 4th element is NaN
+              is_triangle = true;
+
+          if (is_triangle)
+            {
+              // fast way for triangles
+              i1 = f(i,0) - 1; i2 = f(i,1) - 1; i3 = f(i,2) - 1;
+
+              RowVector fn (3, 0.0);
+              double& nx = fn(0);
+              double& ny = fn(1);
+              double& nz = fn(2);
+              if (is_3d)
+                cross_product
+                  (v(i2,0) - v(i1,0), v(i2,1) - v(i1,1), v(i2,2) - v(i1,2),
+                   v(i3,0) - v(i1,0), v(i3,1) - v(i1,1), v(i3,2) - v(i1,2),
+                   nx, ny, nz);
+              else
+                {
+                  nz = (v(i2,0) - v(i1,0)) * (v(i3,1) - v(i1,1)) - 
+                       (v(i2,1) - v(i1,1)) * (v(i3,0) - v(i1,0));
+                  // 2-d vertices always point towards +z
+                  nz = (nz < 0) ? -nz : nz;
+                }
+
+              // assign normal to all edges of triangle
+              for (octave_idx_type j = 0; j < 3; j++)
+                vec_vn[static_cast<octave_idx_type> (f(i,j) - 1)].push_back (fn);
+            }
+          else
+            {
+              // more general for polygons with more than three corners
+
+              // find first element that is NaN to get number of corners
+              octave_idx_type nc = 3;
+              while (!xisnan (f(i,nc)) && nc < max_nc)
+                nc++;
+
+              // calculate normal separately for each corner
+              for (octave_idx_type j = 0; j < nc; j++)
+                {
+                  j1 = j; j2 = j-1; j3 = j+1;
+                  if (j==0)
+                    j2 = nc - 1;
+                  if (j == nc - 1)
+                    j3 = 0;
+
+                  i1 = f(i,j1) - 1; i2 = f(i,j2) - 1; i3 = f(i,j3) - 1;
+                  RowVector fn (3, 0.0);
+                  double& nx = fn(0);
+                  double& ny = fn(1);
+                  double& nz = fn(2);
+                  if (is_3d)
+                    cross_product
+                      (v(i2,0) - v(i1,0), v(i2,1) - v(i1,1), v(i2,2) - v(i1,2),
+                       v(i3,0) - v(i1,0), v(i3,1) - v(i1,1), v(i3,2) - v(i1,2),
+                       nx, ny, nz);
+                  else
+                    {
+                      nz = (v(i2,0) - v(i1,0)) * (v(i3,1) - v(i1,1)) - 
+                           (v(i2,1) - v(i1,1)) * (v(i3,0) - v(i1,0));
+                      // 2-d vertices always point towards +z
+                      nz = (nz < 0) ? -nz : nz;
+                    }
+
+                  // assign normal to this edge of the polygon
+                  vec_vn[static_cast<octave_idx_type> (i1)].push_back (fn);
+                }
+            }
+        }
+
+      // Second step: Calculate the normal for the vertices taking the weighted
+      // average of the normals determined from all adjacent faces (use area of
+      // faces as weights).
+      Matrix vn (num_v, 3, 0.0);
+      for (octave_idx_type i = 0; i < num_v; i++)
+        {
+          std::vector<RowVector>::iterator it = vec_vn[i].begin ();
+
+          // The normal of unused vertices is NaN.
+          RowVector vn0 (3, octave_NaN);
+
+          if (it != vec_vn[i].end ())
+            {
+              // FIXME: Currently, the first vector also determines the
+              // direction of the normal. How to determine the inner and outer
+              // faces of all parts of the patch and point the normals
+              // outwards? (Necessary for correct lighting with
+              // "backfacelighting" set to "lit" or "unlit".) Matlab does not
+              // seem to do it correctly either. So bother here?
+
+              vn0 = *it;
+              
+              for (++it; it != vec_vn[i].end (); ++it)
+                {
+                  RowVector vn1 = *it;
+                  // Use sign of dot product to point vectors in a similar
+                  // direction before taking the average.
+                  double dir = (vn0(0) * vn1(0) + vn0(1) * vn1(1) + vn0(2) * vn1(2) > 0)? 1: -1;
+                  for (octave_idx_type j = 0; j < 3; j++)
+                    vn0(j) += dir * vn1(j);
+                }
+              vn0 /= vec_vn[i].size ();
+              // VertexNormals are not normalized in Matlab.
+              // Should we bother here?
+            }
+          // save normal in matrix
+          for (octave_idx_type j = 0; j < 3; j++)
+            vn(i,j) = vn0(j);
+        }
+
+      vertexnormals = vn;
+    }
+}
+
+
+void
+patch::initialize (const graphics_object& go)
+{
+  base_graphics_object::initialize (go);
+
+  // calculate normals for default data
+  // This is done because the normals for the default data do not match
+  // get(0, "DefaultPatchVertexNormals").
+  xproperties.update_normals ();
+}
+
+
+void
+patch::reset_default_properties (void)
+{
+  // empty list of local defaults
+  default_properties = property_list ();
+  xreset_default_properties (get_handle (), xproperties.factory_defaults ());
+
+  // calculate normals for default data
+  // This is done because the normals for the default data do not match
+  // get(0, "DefaultPatchVertexNormals").
+  xproperties.update_normals ();
+}
+
+// ---------------------------------------------------------------------
+
+octave_value
+surface::properties::get_color_data (void) const
+{
+  return convert_cdata (*this, get_cdata (), cdatamapping_is ("scaled"), 3);
+}
+
+void
 surface::properties::update_normals (void)
 {
   if (normalmode_is ("auto"))
@@ -9651,6 +9821,9 @@
 %!   assert (get (hpa, "vertices"), get (0, "defaultpatchvertices"), tol);
 %!   assert (get (hpa, "facevertexcdata"),
 %!           get (0, "defaultpatchfacevertexcdata"), tol);
+%!   assert (get (0, "defaultpatchvertexnormals"), []);
+%!   assert (get (hpa, "vertexnormals"),
+%!           [0 0 1; 0 0 1; 0 0 1], tol);
 %! unwind_protect_cleanup
 %!   close (hf);
 %! end_unwind_protect
diff -r ffbd3e86e0be -r 408c51b83ae8 libinterp/corefcn/graphics.in.h
--- a/libinterp/corefcn/graphics.in.h	Tue Mar 22 23:14:20 2016 -0700
+++ b/libinterp/corefcn/graphics.in.h	Thu Mar 24 15:32:22 2016 +0100
@@ -4868,11 +4868,11 @@
       color_property markeredgecolor , color_property (radio_values ("none|{auto}|flat"), color_values (0, 0, 0))
       color_property markerfacecolor , color_property (radio_values ("{none}|auto|flat"), color_values (0, 0, 0))
       double_property markersize , 6
-      radio_property normalmode , "{auto}|manual"
+      radio_property normalmode u , "{auto}|manual"
       double_property specularcolorreflectance , 1.0
       double_property specularexponent , 10.0
       double_property specularstrength , 0.9
-      array_property vertexnormals , Matrix ()
+      array_property vertexnormals u , Matrix ()
       array_property vertices u , default_patch_vertices ()
       array_property xdata u , default_patch_xdata ()
       array_property ydata u , default_patch_ydata ()
@@ -4905,15 +4905,19 @@
       facevertexcdata.add_constraint (dim_vector (-1, 1));
       facevertexcdata.add_constraint (dim_vector (-1, 3));
       facevertexalphadata.add_constraint (dim_vector (-1, 1));
-      vertexnormals.add_constraint (dim_vector (-1, -1));
+      vertexnormals.add_constraint (dim_vector (-1, 3));
+      vertexnormals.add_constraint (dim_vector (0, 0));
     }
 
+  public:
+    void update_normals (void);
+
   private:
     std::string bad_data_msg;
 
     void update_faces (void) { update_data ();}
 
-    void update_vertices (void)  {  update_data ();}
+    void update_vertices (void) { update_data ();}
 
     void update_facevertexcdata (void) { update_data ();}
 
@@ -4932,7 +4936,10 @@
           set_faces (Matrix ());
         }
       else
-        update_fvc ();
+        {
+          update_fvc ();
+          update_normals ();
+        }
 
       set_xlim (xdata.get_limits ());
     }
@@ -4947,7 +4954,10 @@
           set_faces (Matrix ());
         }
       else
-        update_fvc ();
+        {
+          update_fvc ();
+          update_normals ();
+        }
 
       set_ylim (ydata.get_limits ());
     }
@@ -4955,12 +4965,14 @@
     void update_zdata (void)
     {
       update_fvc ();
+      update_normals ();
       set_zlim (zdata.get_limits ());
     }
 
     void update_cdata (void)
     {
       update_fvc ();
+      update_normals ();
 
       if (cdatamapping_is ("scaled"))
         set_clim (cdata.get_limits ());
@@ -4968,12 +4980,19 @@
         clim = cdata.get_limits ();
     }
 
-
     void update_data (void);
+
+    void update_normalmode (void)
+    { update_normals (); }
+
+    void update_vertexnormals (void)
+    { set_normalmode ("manual"); }
+
   };
 
 private:
   properties xproperties;
+  property_list default_properties;
 
 public:
   patch (const graphics_handle& mh, const graphics_handle& p)
@@ -4995,6 +5014,12 @@
       retval = base_properties::has_readonly_property (pname);
     return retval;
   }
+
+  void reset_default_properties (void);
+
+protected:
+  void initialize (const graphics_object& go);
+
 };
 
 // ---------------------------------------------------------------------
diff -r ffbd3e86e0be -r 408c51b83ae8 scripts/plot/draw/patch.m
--- a/scripts/plot/draw/patch.m	Tue Mar 22 23:14:20 2016 -0700
+++ b/scripts/plot/draw/patch.m	Thu Mar 24 15:32:22 2016 +0100
@@ -258,6 +258,9 @@
 %!   assert (isempty (get (h, "cdata")));
 %!   assert (get (h, "faces"), [1, 2, 3], eps);
 %!   assert (get (h, "vertices"), [0 1; 1 1; 0 0], eps);
+%!   assert (get (h, "normalmode"), "auto");
+%!   assert (get (0, "defaultpatchvertexnormals"), []);
+%!   assert (get (h, "vertexnormals"), [0 0 1; 0 0 1; 0 0 1], eps);
 %!   assert (get (h, "type"), "patch");
 %!   assert (get (h, "facecolor"), [0 0 0]);
 %!   assert (get (h, "linestyle"), get (0, "defaultpatchlinestyle"));
