# HG changeset patch
# User Markus Mützel <markus.muetzel@gmx.de>
# Date 1527782355 -7200
#      Thu May 31 17:59:15 2018 +0200
# Node ID 12ccf01554120b0309dab6666523a740fe9943f9
# Parent  762eb2e33a7fb57fbf45d843c91a06d53cf5f9a1
Automatically calculate normals for patch objects (patch #8951).

* graphics.in.h, graphics.cc:
(patch properties): Add functions "calc_face_normals", "update_face_normals" and
"update_vertex_normals" to properties.
(base_properties): Add functions "has_child_of_type" and "get_children_of_type".
(patch): Override functions "initialize" and "reset_default_properties".
(axes properties): Add function "trigger_normals_calc".
(light): Override function "initialize".

* patch.m: Add BISTs.

diff -r 762eb2e33a7f -r 12ccf0155412 libinterp/corefcn/graphics.cc
--- a/libinterp/corefcn/graphics.cc	Mon Jul 18 18:33:48 2016 +0200
+++ b/libinterp/corefcn/graphics.cc	Thu May 31 17:59:15 2018 +0200
@@ -3397,6 +3397,69 @@
     p.delete_listener (val, mode);
 }
 
+bool
+base_properties::has_child_of_type (const caseless_str& chtype,
+                                    bool get_invisible = false,
+                                    bool traverse = true) const
+{
+  bool hcot = false;
+  Matrix ch = get_children ();
+  for (int i = 0; i < ch.numel (); i++)
+    {
+      graphics_handle hkid = gh_manager::lookup (ch(i));
+
+      if (hkid.ok ())
+        {
+          graphics_object go = gh_manager::get_object (ch(i));
+          if ( get_invisible || go.get_properties ().is_visible () )
+            {
+              if (go.isa (chtype))
+                {
+                  hcot = true;
+                  break;
+                }
+              else if (traverse & go.isa ("hggroup"))
+                {
+                  hcot = go.get_properties ().has_child_of_type (chtype,
+                                                                 get_invisible,
+                                                                 traverse);
+                  if (hcot)
+                    break;
+                }
+            }
+        }
+    }
+  return hcot;
+}
+
+void
+base_properties::get_children_of_type (const caseless_str& chtype,
+                                       bool get_invisible,
+                                       bool traverse,
+                                       std::list<graphics_object> &children_list) const
+{
+  Matrix ch = get_children ();
+  for (octave_idx_type i = 0; i < ch.numel (); i++)
+    {
+      graphics_handle hkid = gh_manager::lookup (ch(i));
+
+      if (hkid.ok ())
+        {
+          graphics_object go = gh_manager::get_object (hkid);
+          if ( get_invisible || go.get_properties ().is_visible () )
+            {
+              if (go.isa (chtype))
+                children_list.push_back (go);
+              else if (traverse && go.isa ("hggroup"))
+                go.get_properties ().get_children_of_type (chtype,
+                                                           get_invisible,
+                                                           traverse,
+                                                           children_list);
+            }
+        }
+    }
+}
+
 // ---------------------------------------------------------------------
 
 void
@@ -8627,6 +8690,38 @@
 }
 
 void
+axes::properties::trigger_normals_calc (void)
+{
+  // Find all patch (and surface) objects within axes
+  std::list<graphics_object> children_list;
+  std::list<graphics_object>::iterator children_list_iter;
+  get_children_of_type ("patch", false, true, children_list);
+  // FIXME: Un-comment when surface is ready:
+  // get_children_of_type ("surface", false, true, children_list);
+
+  // trigger normals calculation for these objects
+  for (children_list_iter = children_list.begin ();
+       children_list_iter != children_list.end (); children_list_iter++)
+    {
+      graphics_object kid = *children_list_iter;
+      if (kid.isa ("patch"))
+        {
+          patch::properties& patch_props =
+              dynamic_cast<patch::properties&> (kid.get_properties ());
+          patch_props.update_normals (false);
+        }
+      else
+        {
+          // FIXME: Un-comment when surface is ready:
+          // surface::properties& surface_props =
+          //     dynamic_cast<surface::properties&> (kid.get_properties ());
+          // surface_props.update_normals (false);
+        }
+    }
+
+}
+
+void
 axes::reset_default_properties (void)
 {
   // empty list of local defaults
@@ -8875,6 +8970,19 @@
 
 // ---------------------------------------------------------------------
 
+void
+light::initialize (const graphics_object& go)
+{
+  base_graphics_object::initialize (go);
+
+  // trigger normals calculation for the respective children of this axes object
+  axes::properties& parent_axes_prop =
+    dynamic_cast<axes::properties&> (go.get_ancestor ("axes").get_properties ());
+  parent_axes_prop.trigger_normals_calc ();
+}
+
+// ---------------------------------------------------------------------
+
 octave_value
 patch::properties::get_color_data (void) const
 {
@@ -9118,14 +9226,9 @@
   set_ydata (yd);
   set_zdata (zd);
   set_cdata (cd);
-}
-
-// ---------------------------------------------------------------------
-
-octave_value
-surface::properties::get_color_data (void) const
-{
-  return convert_cdata (*this, get_cdata (), cdatamapping_is ("scaled"), 3);
+
+  // Update normals
+  update_normals (true);
 }
 
 inline void
@@ -9139,6 +9242,267 @@
 }
 
 void
+patch::properties::calc_face_normals (Matrix& fn)
+{
+  Matrix v = get_vertices ().matrix_value ();
+  Matrix f = get_faces ().matrix_value ();
+
+  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)
+    {
+      fn = Matrix ();
+      return;
+    }
+
+  // Calculate normals for all faces
+  std::list<std::list<octave_idx_type>>::const_iterator
+    cp_it = coplanar_last_idx.begin ();
+  octave_idx_type i1, i2, i3;
+  octave_idx_type j1, j2;
+  for (octave_idx_type i = 0; i < num_f; i++)
+    {
+      bool is_coplanar = true;
+      if ((*cp_it).size () > 1)
+      {
+        is_coplanar = false;
+      }
+      cp_it++;
+
+      // get number of corners
+      octave_idx_type nc = 3;
+      if (max_nc > 3)
+        {
+          while (! octave::math::isnan (f(i,nc)) && nc < max_nc)
+            nc++;
+        }
+
+      RowVector fnc (3, 0.0);
+      double& nx = fnc(0);
+      double& ny = fnc(1);
+      double& nz = fnc(2);
+
+      if (is_coplanar)
+        {
+          // fast way for coplanar polygons
+          i1 = f(i,0) - 1; i2 = f(i,1) - 1; i3 = f(i,nc-1) - 1;
+
+          if (is_3d)
+            cross_product
+              (v(i3,0) - v(i1,0), v(i3,1) - v(i1,1), v(i3,2) - v(i1,2),
+               v(i2,0) - v(i1,0), v(i2,1) - v(i1,1), v(i2,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;
+            }
+        }
+      else
+        {
+          // more general for non-planar polygons
+
+          // calculate face normal with Newill method
+          // https://courses.cit.cornell.edu/cs417-land/SECTIONS/normals.html
+
+          j1 = nc - 1; j2 = 0;
+          i1 = f(i,j1) - 1; i2 = f(i,j2) - 1;
+
+          nx = (v(i2,1) - v(i1,1)) * (v(i1,2) + v(i2,2));
+          ny = (v(i2,2) - v(i1,2)) * (v(i1,0) + v(i2,0));
+          nz = (v(i2,0) - v(i1,0)) * (v(i1,1) + v(i2,1));
+
+          for (octave_idx_type j = 1; j < nc; j++)
+            {
+              j1 = j-1; j2 = j;
+              i1 = f(i,j1) - 1; i2 = f(i,j2) - 1;
+
+              nx += (v(i2,1) - v(i1,1)) * (v(i1,2) + v(i2,2));
+              ny += (v(i2,2) - v(i1,2)) * (v(i1,0) + v(i2,0));
+              nz += (v(i2,0) - v(i1,0)) * (v(i1,1) + v(i2,1));
+            }
+        }
+
+      // normalize normal vector
+      double n_len = sqrt (nx*nx+ny*ny+nz*nz);
+
+      // assign normal to current face
+      if ( n_len < std::numeric_limits<double>::epsilon () )
+        for (octave_idx_type j = 0; j < 3; j++)
+          fn(i,j) = 0.0;
+      else
+        for (octave_idx_type j = 0; j < 3; j++)
+          fn(i,j) = fnc(j) / n_len;
+    }
+
+}
+
+void
+patch::properties::update_face_normals (bool reset)
+{
+  if (updating_patch_data)
+    return;
+
+  if ( facenormalsmode_is ("auto") && 
+       (facelighting_is ("flat") || edgelighting_is ("flat")) )
+    {
+      graphics_object go = gh_manager::get_object (get___myhandle__ ());
+
+      if (go.get_ancestor ("axes").get_properties ().has_child_of_type ("light", false, true))
+        {
+          Matrix f = get_faces ().matrix_value ();
+
+          octave_idx_type num_f = f.rows (); // number of faces
+          Matrix fn (num_f, 3, 0.0);
+
+          calc_face_normals (fn);
+          facenormals = fn;
+        }
+    }
+    else if ( facenormalsmode_is ("auto") && reset )
+      facenormals = Matrix ();
+
+}
+
+void
+patch::properties::update_vertex_normals (bool reset)
+{
+  if (updating_patch_data)
+    return;
+
+  if ( vertexnormalsmode_is ("auto") &&
+       (facelighting_is ("gouraud") || facelighting_is ("phong") ||
+        edgelighting_is ("gouraud") || edgelighting_is ("phong")) )
+    {
+      graphics_object go = gh_manager::get_object (get___myhandle__ ());
+
+      if (go.get_ancestor ("axes").get_properties ().has_child_of_type ("light", false, true))
+        {
+          Matrix v = get_vertices ().matrix_value ();
+          Matrix f = get_faces ().matrix_value ();
+
+          octave_idx_type num_v = v.rows (); // number of vertices
+          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 for all faces
+          Matrix fn = get_facenormals ().matrix_value ();
+          if ( fn.isempty () )
+            {
+              // calculate facenormals here
+              fn = Matrix (num_f, 3, 0.0);
+              calc_face_normals (fn);
+            }
+
+          // Second step: assign normals to the respective vertices
+          std::vector<RowVector> vec_vn [num_v]; // list of normals for vertices
+          for (octave_idx_type i = 0; i < num_f; i++)
+            {
+              // get number of corners
+              octave_idx_type nc = 3;
+              if (max_nc > 3)
+                {
+                  while (! octave::math::isnan (f(i,nc)) && nc < max_nc)
+                    nc++;
+                }
+
+              for (octave_idx_type j = 0; j < nc; j++)
+                vec_vn[static_cast<octave_idx_type> (f(i,j) - 1)].push_back (fn.row (i));
+            }
+
+          // Third step: Calculate the normal for the vertices taking the average
+          // of the normals determined from all adjacent faces
+          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);
+                    }
+
+                  // normalize normal vector
+                  double n_len = sqrt (vn0(0)*vn0(0)+vn0(1)*vn0(1)+vn0(2)*vn0(2));
+
+                  // save normal in matrix
+                  for (octave_idx_type j = 0; j < 3; j++)
+                    vn(i,j) = vn0(j)/n_len;
+                }
+            }
+
+          vertexnormals = vn;
+        }
+    }
+    else if ( vertexnormalsmode_is ("auto") && reset )
+      vertexnormals = Matrix ();
+
+}
+
+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") in Matlab.
+  xproperties.update_normals (true);
+}
+
+
+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") in Matlab.
+  xproperties.update_normals (true);
+}
+
+// ---------------------------------------------------------------------
+
+octave_value
+surface::properties::get_color_data (void) const
+{
+  return convert_cdata (*this, get_cdata (), cdatamapping_is ("scaled"), 3);
+}
+
+void
 surface::properties::update_vertex_normals (void)
 {
   if (vertexnormalsmode_is ("auto"))
diff -r 762eb2e33a7f -r 12ccf0155412 libinterp/corefcn/graphics.in.h
--- a/libinterp/corefcn/graphics.in.h	Mon Jul 18 18:33:48 2016 +0200
+++ b/libinterp/corefcn/graphics.in.h	Thu May 31 17:59:15 2018 +0200
@@ -2287,6 +2287,10 @@
     return children.get_hidden ();
   }
 
+  bool has_child_of_type (const caseless_str& type, bool get_invisible, bool traverse) const;
+
+  void get_children_of_type (const caseless_str& type, bool get_invisible, bool traverse, std::list<graphics_object> &children_list) const;
+
   void set_modified (const octave_value& val) { set___modified__ (val); }
 
   void set___modified__ (const octave_value& val) { __modified__ = val; }
@@ -4175,6 +4179,8 @@
       update_axes_layout ();
     }
 
+    void trigger_normals_calc (void);
+
   };
 
 private:
@@ -4804,6 +4811,9 @@
       retval = base_properties::has_readonly_property (pname);
     return retval;
   }
+
+protected:
+  void initialize (const graphics_object& go);
 };
 
 // ---------------------------------------------------------------------
@@ -4850,13 +4860,13 @@
       string_property displayname , ""
       double_radio_property edgealpha , double_radio_property (1.0, radio_values ("flat|interp"))
       color_property edgecolor , color_property (color_values (0, 0, 0), radio_values ("none|flat|interp"))
-      radio_property edgelighting , "{none}|flat|gouraud|phong"
+      radio_property edgelighting u , "{none}|flat|gouraud|phong"
       radio_property erasemode h , "{normal}|none|xor|background"
       double_radio_property facealpha , double_radio_property (1.0, radio_values ("flat|interp"))
       color_property facecolor , color_property (color_values (0, 0, 0), radio_values ("none|flat|interp"))
-      radio_property facelighting , "none|{flat}|gouraud|phong"
+      radio_property facelighting u , "none|{flat}|gouraud|phong"
       array_property facenormals m , Matrix ()
-      radio_property facenormalsmode , "{auto}|manual"
+      radio_property facenormalsmode u , "{auto}|manual"
       array_property faces u , default_patch_faces ()
       array_property facevertexalphadata , Matrix ()
       array_property facevertexcdata u , Matrix ()
@@ -4872,7 +4882,7 @@
       double_property specularexponent , 10.0
       double_property specularstrength , 0.9
       array_property vertexnormals m , Matrix ()
-      radio_property vertexnormalsmode , "{auto}|manual"
+      radio_property vertexnormalsmode u , "{auto}|manual"
       array_property vertices u , default_patch_vertices ()
       array_property xdata u , default_patch_xdata ()
       array_property ydata u , default_patch_ydata ()
@@ -4923,12 +4933,20 @@
       specularstrength.add_constraint ("max", 1.0, true);
     }
 
+  public:
+    void update_normals (bool reset)
+    {
+      update_face_normals (reset);
+      update_vertex_normals (reset);
+    }
+
+
   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 ();}
 
@@ -4947,7 +4965,10 @@
           set_faces (Matrix ());
         }
       else
-        update_fvc ();
+        {
+          update_fvc ();
+          update_normals (true);
+        }
 
       set_xlim (xdata.get_limits ());
     }
@@ -4962,7 +4983,10 @@
           set_faces (Matrix ());
         }
       else
-        update_fvc ();
+        {
+          update_fvc ();
+          update_normals (true);
+        }
 
       set_ylim (ydata.get_limits ());
     }
@@ -4970,12 +4994,14 @@
     void update_zdata (void)
     {
       update_fvc ();
+      update_normals (true);
       set_zlim (zdata.get_limits ());
     }
 
     void update_cdata (void)
     {
       update_fvc ();
+      update_normals (false);
 
       if (cdatamapping_is ("scaled"))
         set_clim (cdata.get_limits ());
@@ -4984,10 +5010,35 @@
     }
 
     void update_data (void);
+
+    void calc_face_normals (Matrix& normals);
+    void update_face_normals (bool reset);
+    void update_vertex_normals (bool reset);
+
+    void update_edgelighting (void)
+    {
+      update_normals (false);
+    }
+
+    void update_facelighting (void)
+    {
+      update_normals (false);
+    }
+
+    void update_facenormalsmode (void)
+    {
+      update_face_normals (false);
+    }
+
+    void update_vertexnormalsmode (void)
+    {
+      update_vertex_normals (false);
+    }
   };
 
 private:
   properties xproperties;
+  property_list default_properties;
 
 public:
   patch (const graphics_handle& mh, const graphics_handle& p)
@@ -5009,6 +5060,12 @@
       retval = base_properties::has_readonly_property (pname);
     return retval;
   }
+
+  void reset_default_properties (void);
+
+protected:
+  void initialize (const graphics_object& go);
+
 };
 
 // ---------------------------------------------------------------------
diff -r 762eb2e33a7f -r 12ccf0155412 scripts/plot/draw/patch.m
--- a/scripts/plot/draw/patch.m	Mon Jul 18 18:33:48 2016 +0200
+++ b/scripts/plot/draw/patch.m	Thu May 31 17:59:15 2018 +0200
@@ -282,12 +282,22 @@
 %!   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, "facenormalsmode"), "auto");
+%!   assert (get (h, "facenormals"), get (0, "defaultpatchfacenormals"));
+%!   assert (get (h, "vertexnormalsmode"), "auto");
+%!   assert (get (h, "vertexnormals"), get (0, "defaultpatchvertexnormals"));
 %!   assert (get (h, "type"), "patch");
 %!   assert (get (h, "facecolor"), [0 0 0]);
 %!   assert (get (h, "linestyle"), get (0, "defaultpatchlinestyle"));
 %!   assert (get (h, "linewidth"), get (0, "defaultpatchlinewidth"), eps);
 %!   assert (get (h, "marker"), get (0, "defaultpatchmarker"));
 %!   assert (get (h, "markersize"), get (0, "defaultpatchmarkersize"));
+%!   hl = light;
+%!   assert (get (h, "facelighting"), "flat");
+%!   assert (get (h, "facenormals"), [0 0 1], eps);
+%!   assert (get (h, "vertexnormals"), []);
+%!   set (h, "facelighting", "gouraud")
+%!   assert (get (h, "vertexnormals"), [0 0 1; 0 0 1; 0 0 1], eps);
 %! unwind_protect_cleanup
 %!   close (hf);
 %! end_unwind_protect
