[changeset] surface normals

[Kai Habel]

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

# HG changeset patch
# User Kai Habel
# Date 1226242704 -3600
# Node ID 36c59eb57b3764c762d273499c2730c47093ab04
# Parent  5fe0f4dfdbecf82cac0fa6112ab2f8e99b13167c
Calculate surface normals for boundaries, use more neighboring
points for interior normals

diff -r 5fe0f4dfdbec -r 36c59eb57b37 src/ChangeLog
--- a/src/ChangeLog     Sun Nov 09 12:59:42 2008 +0100
+++ b/src/ChangeLog     Sun Nov 09 15:58:24 2008 +0100
@@ -1,3 +1,8 @@
+2008-11-09  Kai Habel <address@hidden>
+
+       * graphics.cc: Calculate normals for surface boundaries, use
+        more neighboring points to calculate interior normals
+
 2008-10-31  Jaroslav Hajek <address@hidden>
 
        * xnorm.cc: New source.
diff -r 5fe0f4dfdbec -r 36c59eb57b37 src/graphics.cc
--- a/src/graphics.cc   Sun Nov 09 12:59:42 2008 +0100
+++ b/src/graphics.cc   Sun Nov 09 15:58:24 2008 +0100
@@ -3594,55 +3594,66 @@
       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;
-
-      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++)
+      int nc = z.columns (), nr = z.rows ();
+      int c1, c2, c3;
+      int r1, r2, r3;
+
+      bool x_mat = (x.rows () == nr);
+      bool y_mat = (y.columns () == nc);
+
+      NDArray n (dim_vector (nr, nc, 3), 0.0);
+
+      c1 = c2 = c3 = 0;
+      r1 = r2 = r3 = 0;
+
+      for (int c = 0; c < nc; c++)
        {
          if (y_mat)
            {
-             i1 = i-1;
-             i2 = i;
-             i3 = i+1;
-           }
-
-         for (int j = 1; j < (q-1); j++)
+             c1 = c-1;
+             c2 = c;
+             c3 = c+1;
+           }
+
+         for (int r = 0; r < nr; r++)
            {
              if (x_mat)
                {
-                 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);
+                 r1 = r-1;
+                 r2 = r;
+                 r3 = r+1;
+               }
+
+             double& nx = n(r,c,0);
+             double& ny = n(r,c,1);
+             double& nz = n(r,c,2);
+
+              if ((r > 0) && (c > 0))
+                  // upper left quadrangle
+                 cross_product (x(r1,c-1)-x(r2,c), y(r-1,c1)-y(r,c2), 
z(r-1,c-1)-z(r,c),
+                                x(r2,c-1)-x(r1,c), y(r,c1)-y(r-1,c2), 
z(r,c-1)-z(r-1,c),
+                                nx, ny, nz);
+
+              if ((r > 0) && (c < (nc -1)))
+                  // upper right quadrangle
+                  cross_product (x(r1,c+1)-x(r2,c), y(r-1,c3)-y(r,c2), 
z(r-1,c+1)-z(r,c),
+                                x(r1,c)-x(r2,c+1), y(r-1,c2)-y(r,c3), 
z(r-1,c)-z(r,c+1),
+                                nx, ny, nz);
+
+              if ((r < (nr - 1)) && (c > 0))
+                  // lower left quadrangle
+
+                  cross_product (x(r2,c-1)-x(r3,c), y(r,c1)-y(r+1,c2), 
z(r,c-1)-z(r+1,c),
+                                x(r3,c-1)-x(r2,c), y(r+1,c1)-y(r,c2), 
z(r+1,c-1)-z(r,c),
+                                nx, ny, nz);
+
+              if ((r < (nr - 1)) && (c < (nc -1)))
+                  // lower right quadrangle
+                 cross_product (x(r3,c)-x(r2,c+1), y(r+1,c2)-y(r,c3), 
z(r+1,c)-z(r,c+1),
+                                 x(r3,c+1)-x(r2,c), y(r+1,c3)-y(r,c2), 
z(r+1,c+1)-z(r,c),
+                                nx, ny, nz);
+              double d = - std::max(std::max(fabs(nx), fabs(ny)), fabs(nz));
+              //double d = - sqrt (nx*nx + ny*ny + nz*nz);
 
              nx /= d;
              ny /= d;