[changeset] surface normals (again)

[Kai Habel]

Hello all,

I have already send this patch in November [1], but got no reply from
the core developers. In the meantime I have renamed some of the
variables but made no major changes otherwise. So here is it again:

> Hello,
> 
> the attached changeset adds the calculation of normals for surface
> boundaries. For the calculation of interior normals four additional
> neighbors are being used.
> 
> Kai
> 
> 

I have changed the normalization to be more matlab compatible see this
remark also send in November [1]

> Just a quick comment on the change below
> 
> +              double d = - std::max(std::max(fabs(nx), fabs(ny)), fabs(nz));
> +              //double d = - sqrt (nx*nx + ny*ny + nz*nz);
> 
> It seems matlab calculates the normal vector in manner that the largest 
> element is set to 1.0. We have set the length to 1.0 previously. I think the 
> reason is to make the vector more human readable like
>> v=[1 -0.5 0.25]
> v =
> 
>   1.00000  -0.50000   0.25000
> 
> 
> versus
> 
>> v/norm(v)
> ans =
> 
>   0.87287  -0.43644   0.21822
> 
>>


Kai


[1] http://www.nabble.com/-changeset--surface-normals-to20407855.html

# HG changeset patch
# User Kai Habel
# Date 1231421722 -3600
# Node ID f53edfb28c12db747388c643b20bab1b334e7291
# Parent  a6b4d8fdbea1d2d725727f3db21e5a81013b5aa7
Calculate surface normals for boundaries, use more neighboring
points for interior normals

diff -r a6b4d8fdbea1 -r f53edfb28c12 src/ChangeLog
--- a/src/ChangeLog     Sun Dec 28 17:39:20 2008 +0100
+++ b/src/ChangeLog     Thu Jan 08 14:35:22 2009 +0100
@@ -1,3 +1,8 @@
+2009-01-08  Kai Habel <address@hidden>
+
+        * graphics.cc: Calculate normals for surface boundaries, use
+        more neighboring points to calculate interior normals
+
 2008-12-27  Jaroslav Hajek <address@hidden>
 
        * oct-obj.h, oct-obj.cc (octave_value_list::valid_scalar_indices): 
Remove.
diff -r a6b4d8fdbea1 -r f53edfb28c12 src/graphics.cc
--- a/src/graphics.cc   Sun Dec 28 17:39:20 2008 +0100
+++ b/src/graphics.cc   Thu Jan 08 14:35:22 2009 +0100
@@ -3594,62 +3594,69 @@
       Matrix y = get_ydata ().matrix_value ();
       Matrix z = get_zdata ().matrix_value ();
 
+
       int p = z.columns (), q = z.rows ();
-      int i1, i2, i3;
-      int j1, j2, j3;
+      int i1 = 0, i2 = 0, i3 = 0;
+      int j1 = 0, j2 = 0, j3 = 0;
 
       bool x_mat = (x.rows () == q);
       bool y_mat = (y.columns () == p);
 
       NDArray n (dim_vector (q, p, 3), 0.0);
 
-      i1 = i2 = i3 = 0;
-      j1 = j2 = j3 = 0;
-
-      // FIXME -- normal computation at boundaries.
-      for (int i = 1; i < (p-1); i++)
+      for (int i = 0; i < p; i++)
        {
          if (y_mat)
            {
-             i1 = i-1;
+             i1 = i - 1;
              i2 = i;
-             i3 = i+1;
-           }
-
-         for (int j = 1; j < (q-1); j++)
+             i3 = i + 1;
+           }
+
+         for (int j = 0; j < q; j++)
            {
              if (x_mat)
                {
-                 j1 = j-1;
+                 j1 = j - 1;
                  j2 = j;
-                 j3 = j+1;
-               }
-
-             double& nx = n(j,i,0);
-             double& ny = n(j,i,1);
-             double& nz = n(j,i,2);
-
-             cross_product (x(j3,i)-x(j2,i), y(j+1,i2)-y(j,i2), 
z(j+1,i)-z(j,i),
-                            x(j2,i+1)-x(j2,i), y(j,i3)-y(j,i2), 
z(j,i+1)-z(i,j),
-                            nx, ny, nz);
-             cross_product (x(j2,i-1)-x(j2,i), y(j,i1)-y(j,i2), 
z(j,i-1)-z(j,i),
-                            x(j3,i)-x(j2,i), y(j+1,i2)-y(j,i2), 
z(j+1,i)-z(i,j),
-                            nx, ny, nz);
-             cross_product (x(j1,i)-x(j2,i), y(j-1,i2)-y(j,i2), 
z(j-1,i)-z(j,i),
-                            x(j2,i-1)-x(j2,i), y(j,i1)-y(j,i2), 
z(j,i-1)-z(i,j),
-                            nx, ny, nz);
-             cross_product (x(j2,i+1)-x(j2,i), y(j,i3)-y(j,i2), 
z(j,i+1)-z(j,i),
-                            x(j1,i)-x(j2,i), y(j-1,i2)-y(j,i2), 
z(j-1,i)-z(i,j),
-                            nx, ny, nz);
-
-             double d = - sqrt (nx*nx + ny*ny + nz*nz);
+                 j3 = j + 1;
+               }
+
+             double& nx = n(j, i, 0);
+             double& ny = n(j, i, 1);
+             double& nz = n(j, i, 2);
+
+              if ((j > 0) && (i > 0))
+                  // upper left quadrangle
+                 cross_product (x(j1,i-1)-x(j2,i), y(j-1,i1)-y(j,i2), 
z(j-1,i-1)-z(j,i),
+                                x(j2,i-1)-x(j1,i), y(j,i1)-y(j-1,i2), 
z(j,i-1)-z(j-1,i),
+                                nx, ny, nz);
+
+              if ((j > 0) && (i < (p -1)))
+                  // upper right quadrangle
+                  cross_product (x(j1,i+1)-x(j2,i), y(j-1,i3)-y(j,i2), 
z(j-1,i+1)-z(j,i),
+                                x(j1,i)-x(j2,i+1), y(j-1,i2)-y(j,i3), 
z(j-1,i)-z(j,i+1),
+                                nx, ny, nz);
+
+              if ((j < (q - 1)) && (i > 0))
+                  // lower left quadrangle
+                  cross_product (x(j2,i-1)-x(j3,i), y(j,i1)-y(j+1,i2), 
z(j,i-1)-z(j+1,i),
+                                x(j3,i-1)-x(j2,i), y(j+1,i1)-y(j,i2), 
z(j+1,i-1)-z(j,i),
+                                nx, ny, nz);
+
+              if ((j < (q - 1)) && (i < (p -1)))
+                  // lower right quadrangle
+                 cross_product (x(j3,i)-x(j2,i+1), y(j+1,i2)-y(j,i3), 
z(j+1,i)-z(j,i+1),
+                                 x(j3,i+1)-x(j2,i), y(j+1,i3)-y(j,i2), 
z(j+1,i+1)-z(j,i),
+                                nx, ny, nz);
+
+              double d = - std::max(std::max(fabs(nx), fabs(ny)), fabs(nz));
 
              nx /= d;
              ny /= d;
              nz /= d;
            }
        }
-
       vertexnormals = n;
     }
 }