ENH: Finalized cell interface.

src/Cell.cc

@@ -168,3 +168,26 @@ float vtkCell::GetLength2 ()
  
   return l;
 }
+
+void vtkCell::PrintSelf(ostream& os, vtkIndent indent)
+{
+  float *bounds;
+  int i;
+
+  vtkObject::PrintSelf(os,indent);
+
+  os << indent << "Number Of Points: " << this->PointIds.GetNumberOfIds() << "\n";
+  bounds = this->GetBounds();
+  os << indent << "Bounds: \n";
+  os << indent << "  Xmin,Xmax: (" << bounds[0] << ", " << bounds[1] << ")\n";
+  os << indent << "  Ymin,Ymax: (" << bounds[2] << ", " << bounds[3] << ")\n";
+  os << indent << "  Zmin,Zmax: (" << bounds[4] << ", " << bounds[5] << ")\n";
+
+  os << indent << "  Point ids are: ";
+  for (i=0; this->PointIds.GetNumberOfIds(); i++)
+    {
+    os << ", " << this->PointIds.GetId(i);
+    if ( i && !(i % 12) ) os << "\n\t";
+    }
+  os << indent << "\n";
+}

src/Hexa.cc

@@ -84,11 +84,7 @@ int vtkHexahedron::EvaluatePosition(float x[3], float closestPoint[3],
 //
 //  compute determinants and generate improvements
 //
-    if ( (d=math.Determinant3x3(rcol,scol,tcol)) == 0.0 )
-      {
-      dist2 = LARGE_FLOAT;
-      return 0;
-      }
+    if ( (d=math.Determinant3x3(rcol,scol,tcol)) == 0.0 ) return -1;
 
     pcoords[0] = params[0] - math.Determinant3x3 (fcol,scol,tcol) / d;
     pcoords[1] = params[1] - math.Determinant3x3 (rcol,fcol,tcol) / d;
@@ -116,33 +112,30 @@ int vtkHexahedron::EvaluatePosition(float x[3], float closestPoint[3],
 //  if not converged, set the parametric coordinates to arbitrary values
 //  outside of element
 //
-  if ( !converged )
+  if ( !converged ) return -1;
+
+  this->InterpolationFunctions(pcoords, weights);
+
+  if ( pcoords[0] >= -0.001 && pcoords[0] <= 1.001 &&
+  pcoords[1] >= -0.001 && pcoords[1] <= 1.001 &&
+  pcoords[2] >= -0.001 && pcoords[2] <= 1.001 )
     {
-    pcoords[0] = pcoords[1] =  pcoords[2] = 10.0;
-    dist2 = LARGE_FLOAT;
-    return 0;
+    closestPoint[0] = x[0]; closestPoint[1] = x[1]; closestPoint[2] = x[2];
+    dist2 = 0.0; //inside hexahedron
+    return 1;
     }
   else
     {
-    if ( pcoords[0] >= -0.001 && pcoords[0] <= 1.001 &&
-    pcoords[1] >= -0.001 && pcoords[1] <= 1.001 &&
-    pcoords[2] >= -0.001 && pcoords[2] <= 1.001 )
-      {
-      closestPoint[0] = x[0]; closestPoint[1] = x[1]; closestPoint[2] = x[2];
-      dist2 = 0.0; // inside hexahedron
-      return 1;
-      }
-    else
+    float pc[3], w[8];
+    for (i=0; i<3; i++) //only approximate, not really true for warped hexa
       {
-      for (i=0; i<3; i++)
-        {
-        if (pcoords[i] < 0.0) pcoords[i] = 0.0;
-        if (pcoords[i] > 1.0) pcoords[i] = 1.0;
-        }
-      this->EvaluateLocation(subId, pcoords, closestPoint, weights);
-      dist2 = math.Distance2BetweenPoints(closestPoint,x);
-      return 0;
+      if (pcoords[i] < 0.0) pc[i] = 0.0;
+      else if (pcoords[i] > 1.0) pc[i] = 1.0;
+      else pc[i] = pcoords[i];
       }
+    this->EvaluateLocation(subId, pc, closestPoint, w);
+    dist2 = math.Distance2BetweenPoints(closestPoint,x);
+    return 0;
     }
 }
 //
@@ -434,3 +427,75 @@ int vtkHexahedron::IntersectWithLine(float p1[3], float p2[3], float tol,
     }
   return intersection;
 }
+
+int vtkHexahedron::Triangulate(int index, vtkFloatPoints &pts)
+{
+  pts.Reset();
+//
+// Create five tetrahedron. Triangulation varies depending upon index. This
+// is necessary to insure compatible voxel triangulations.
+//
+  if ( index % 2 )
+    {
+    pts.InsertNextPoint(this->Points.GetPoint(0));
+    pts.InsertNextPoint(this->Points.GetPoint(1));
+    pts.InsertNextPoint(this->Points.GetPoint(4));
+    pts.InsertNextPoint(this->Points.GetPoint(3));
+
+    pts.InsertNextPoint(this->Points.GetPoint(1));
+    pts.InsertNextPoint(this->Points.GetPoint(4));
+    pts.InsertNextPoint(this->Points.GetPoint(7));
+    pts.InsertNextPoint(this->Points.GetPoint(5));
+
+    pts.InsertNextPoint(this->Points.GetPoint(1));
+    pts.InsertNextPoint(this->Points.GetPoint(4));
+    pts.InsertNextPoint(this->Points.GetPoint(3));
+    pts.InsertNextPoint(this->Points.GetPoint(6));
+
+    pts.InsertNextPoint(this->Points.GetPoint(1));
+    pts.InsertNextPoint(this->Points.GetPoint(3));
+    pts.InsertNextPoint(this->Points.GetPoint(2));
+    pts.InsertNextPoint(this->Points.GetPoint(6));
+
+    pts.InsertNextPoint(this->Points.GetPoint(3));
+    pts.InsertNextPoint(this->Points.GetPoint(6));
+    pts.InsertNextPoint(this->Points.GetPoint(4));
+    pts.InsertNextPoint(this->Points.GetPoint(7));
+    }
+  else
+    {
+    pts.InsertNextPoint(this->Points.GetPoint(2));
+    pts.InsertNextPoint(this->Points.GetPoint(1));
+    pts.InsertNextPoint(this->Points.GetPoint(0));
+    pts.InsertNextPoint(this->Points.GetPoint(5));
+
+    pts.InsertNextPoint(this->Points.GetPoint(0));
+    pts.InsertNextPoint(this->Points.GetPoint(2));
+    pts.InsertNextPoint(this->Points.GetPoint(7));
+    pts.InsertNextPoint(this->Points.GetPoint(3));
+
+    pts.InsertNextPoint(this->Points.GetPoint(2));
+    pts.InsertNextPoint(this->Points.GetPoint(5));
+    pts.InsertNextPoint(this->Points.GetPoint(7));
+    pts.InsertNextPoint(this->Points.GetPoint(6));
+
+    pts.InsertNextPoint(this->Points.GetPoint(0));
+    pts.InsertNextPoint(this->Points.GetPoint(7));
+    pts.InsertNextPoint(this->Points.GetPoint(5));
+    pts.InsertNextPoint(this->Points.GetPoint(4));
+
+    pts.InsertNextPoint(this->Points.GetPoint(1));
+    pts.InsertNextPoint(this->Points.GetPoint(2));
+    pts.InsertNextPoint(this->Points.GetPoint(5));
+    pts.InsertNextPoint(this->Points.GetPoint(7));
+    }
+
+  return 1;
+}
+
+void vtkHexahedron::Derivatives(int subId, float pcoords[3], float *values, 
+                                int dim, float *derivs)
+{
+
+}
+

src/ImpMod.cc

@@ -148,9 +148,10 @@ void vtkImplicitModeller::Execute()
           x[0] = this->AspectRatio[0] * i + this->Origin[0];
           idx = jkFactor*k + this->SampleDimensions[0]*j + i;
           prevDistance2 = newScalars->GetScalar(idx);
-          cell->EvaluatePosition(x, closestPoint, subId, pcoords, 
-                                 distance2, weights);
-          if (distance2 < prevDistance2)
+
+          // union combination of distances
+          if ( cell->EvaluatePosition(x, closestPoint, subId, pcoords, 
+          distance2, weights) != -1 && distance2 < prevDistance2 )
             newScalars->SetScalar(idx,distance2);
           }
         }

src/ImpTC.cc

@@ -100,16 +100,16 @@ void vtkImplicitTextureCoords::Execute()
     scale[i] = 1.0;
     if ( max[i] > 0.0 && min[i] < 0.0 ) //have positive & negative numbers
       {
-      if ( max[i] > (-min[i]) ) scale[i] = 0.5 / max[i]; //scale into 0.5->1
-      else scale[i] = -0.5 / min[i]; //scale into 0->0.5
+      if ( max[i] > (-min[i]) ) scale[i] = 0.499 / max[i]; //scale into 0.5->1
+      else scale[i] = -0.499 / min[i]; //scale into 0->0.5
       }
     else if ( max[i] > 0.0 ) //have positive numbers only
       {
-      scale[i] = 0.5 / max[i]; //scale into 0.5->1.0
+      scale[i] = 0.499 / max[i]; //scale into 0.5->1.0
       }
     else if ( min[i] < 0.0 ) //have negative numbers only
       {
-      scale[i] = -0.5 / min[i]; //scale into 0.0->0.5
+      scale[i] = -0.499 / min[i]; //scale into 0.0->0.5
       }
     }
 

src/Line.cc

@@ -17,7 +17,7 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
 #include "vtkMath.hh"
 #include "CellArr.hh"
 
-static vtkMath math;
+static vtkMath math; //avoid lots of construction
 
 // Description:
 // Deep copy of cell.
@@ -35,60 +35,22 @@ int vtkLine::EvaluatePosition(float x[3], float closestPoint[3],
                              int& subId, float pcoords[3],
                              float& dist2, float weights[MAX_CELL_SIZE])
 {
-  float *a1, *a2, a21[3], denom, num;
-  int i, return_status;
-  float *closest;
+  float *a1, *a2;
+  int i;
 
   subId = 0;
   pcoords[1] = pcoords[2] = 0.0;
 
   a1 = this->Points.GetPoint(0);
   a2 = this->Points.GetPoint(1);
-//
-//   Determine appropriate vectors
-// 
-  for (i=0; i<3; i++) a21[i] = a2[i] - a1[i];
-//
-//   Get parametric location
-//
-  num = a21[0]*(x[0]-a1[0]) + a21[1]*(x[1]-a1[1]) + a21[2]*(x[2]-a1[2]);
-  denom = math.Dot(a21,a21);
 
-  if ( (denom = math.Dot(a21,a21)) < fabs(TOL*num) )
-    {
-    dist2 = LARGE_FLOAT;
-    }
-  else 
-    {
-    pcoords[0] = num / denom;
-    }
-//
-// If parametric coordinate is within 0<=p<=1, then the point is closest to
-// the line.  Otherwise, it's closest to a point at the end of the line.
-//
-  if ( pcoords[0] < 0.0 )
-    {
-    closest = a1;
-    return_status = 0;
-    }
-  else if ( pcoords[0] > 1.0 )
-    {
-    closest = a2;
-    return_status = 0;
-    }
-  else
-    {
-    closest = a21;
-    for (i=0; i<3; i++) a21[i] = a1[i] + pcoords[0]*a21[i];
-    return_status = 1;
-    }
+  dist2 = this->DistanceToLine(x,a1,a2,pcoords[0],closestPoint);
 
-  dist2 = math.Distance2BetweenPoints(closest,x);
-  closestPoint[0] = closest[0]; closestPoint[1] = closest[1]; closestPoint[2] = closest[2]; 
   weights[0] = pcoords[0];
   weights[1] = 1.0 - pcoords[0];
 
-  return return_status;
+  if ( pcoords[0] < 0.0 || pcoords[0] > 1.0 ) return 0;
+  else return 1;
 }
 
 void vtkLine::EvaluateLocation(int& subId, float pcoords[3], float x[3],
@@ -234,9 +196,58 @@ void vtkLine::Contour(float value, vtkFloatScalars *cellScalars,
     }
 }
 
+// Description:
+// Compute distance to finite line.
+float vtkLine::DistanceToLine(float x[3], float p1[3], float p2[3], 
+                              float &t, float closestPoint[3])
+{
+  int i;
+  float p21[3], denom, num;
+  float *closest;
+//
+//   Determine appropriate vectors
+// 
+  for (i=0; i<3; i++) p21[i] = p2[i] - p1[i];
+//
+//   Get parametric location
+//
+  num = p21[0]*(x[0]-p1[0]) + p21[1]*(x[1]-p1[1]) + p21[2]*(x[2]-p1[2]);
+  denom = math.Dot(p21,p21);
+
+  if ( (denom = math.Dot(p21,p21)) < fabs(TOL*num) ) //numerically bad!
+    {
+    closest = p1; //arbitrary, point is (numerically) far away
+    }
+//
+// If parametric coordinate is within 0<=p<=1, then the point is closest to
+// the line.  Otherwise, it's closest to a point at the end of the line.
+//
+  else if ( (t=num/denom) < 0.0 )
+    {
+    closest = p1;
+    }
+  else if ( t > 1.0 )
+    {
+    closest = p2;
+    }
+  else
+    {
+    closest = p21;
+    for (i=0; i<3; i++) p21[i] = p1[i] + t*p21[i];
+    }
+
+  closestPoint[0] = closest[0]; 
+  closestPoint[1] = closest[1]; 
+  closestPoint[2] = closest[2]; 
+
+  return math.Distance2BetweenPoints(closest,x);
+}
+
 //
-//  Determine the distance of the current vertex to the edge defined by
-//  the vertices provided.  Returns distance squared.
+// Description:
+// Determine the distance of the current vertex to the edge defined by
+// the vertices provided.  Returns distance squared. Note: line is assumed
+// infinite in extent.
 //
 float vtkLine::DistanceToLine (float x[3], float p1[3], float p2[3])
 {
@@ -294,3 +305,19 @@ int vtkLine::IntersectWithLine(float p1[3], float p2[3], float tol, float& t,
       return 0;
     }
 }
+
+int vtkLine::Triangulate(int index, vtkFloatPoints &pts)
+{
+  pts.Reset();
+  pts.InsertPoint(0,this->Points.GetPoint(0));
+  pts.InsertPoint(1,this->Points.GetPoint(1));
+
+  return 1;
+}
+
+void vtkLine::Derivatives(int subId, float pcoords[3], float *values, 
+                          int dim, float *derivs)
+{
+
+}
+

src/Pixel.cc

@@ -70,21 +70,24 @@ int vtkPixel::EvaluatePosition(float x[3], float closestPoint[3],
   pcoords[0] = math.Dot(p21,p) / (l21*l21);
   pcoords[1] = math.Dot(p31,p) / (l31*l31);
 
+  this->InterpolationFunctions(pcoords, weights);
+
   if ( pcoords[0] >= 0.0 && pcoords[1] <= 1.0 &&
   pcoords[1] >= 0.0 && pcoords[1] <= 1.0 )
     {
     dist2 = math.Distance2BetweenPoints(closestPoint,x); //projection distance
-    this->InterpolationFunctions(pcoords, weights);
     return 1;
     }
   else
     {
+    float pc[3], w[4];
     for (i=0; i<2; i++)
       {
-      if (pcoords[i] < 0.0) pcoords[i] = 0.0;
-      if (pcoords[i] > 1.0) pcoords[i] = 1.0;
+      if (pcoords[i] < 0.0) pc[i] = 0.0;
+      else if (pcoords[i] > 1.0) pc[i] = 1.0;
+      else pc[i] = pcoords[i];
       }
-    this->EvaluateLocation(subId, pcoords, closestPoint, weights);
+    this->EvaluateLocation(subId, pc, closestPoint, w);
     dist2 = math.Distance2BetweenPoints(closestPoint,x);
     return 0;
     }
@@ -291,3 +294,24 @@ int vtkPixel::IntersectWithLine(float p1[3], float p2[3], float tol, float& t,
 
   return 0;
 }
+
+int vtkPixel::Triangulate(int index, vtkFloatPoints &pts)
+{
+  pts.Reset();
+  pts.InsertPoint(0,this->Points.GetPoint(0));
+  pts.InsertPoint(1,this->Points.GetPoint(1));
+  pts.InsertPoint(2,this->Points.GetPoint(2));
+
+  pts.InsertPoint(3,this->Points.GetPoint(1));
+  pts.InsertPoint(4,this->Points.GetPoint(3));
+  pts.InsertPoint(5,this->Points.GetPoint(2));
+
+  return 1;
+}
+
+void vtkPixel::Derivatives(int subId, float pcoords[3], float *values, 
+                            int dim, float *derivs)
+{
+
+}
+

src/PointSet.cc

@@ -108,8 +108,8 @@ int vtkPointSet::FindCell(float x[3], vtkCell *cell, float tol2, int& subId,
       {
       cellId = cellIds.GetId(i);
       cell = this->GetCell(cellId);
-      cell->EvaluatePosition(x,closestPoint,sId,pc,dist2,w);
-      if ( dist2 <= tol2 && dist2 < minDist2 )
+      if ( cell->EvaluatePosition(x,closestPoint,sId,pc,dist2,w) != -1 &&
+      dist2 <= tol2 && dist2 < minDist2 )
         {
         minDist2 = dist2;
         closestCell = cellId;

src/PolyLine.cc

@@ -356,7 +356,7 @@ int vtkPolyLine::EvaluatePosition(float x[3], float closestPoint[3],
     line.Points.SetPoint(0,this->Points.GetPoint(i));
     line.Points.SetPoint(1,this->Points.GetPoint(i+1));
     status = line.EvaluatePosition(x,closest,ignoreId,pc,dist2,lineWeights);
-    if ( dist2 < minDist2 )
+    if ( status != -1 && dist2 < minDist2 )
       {
       return_status = status;
       closestPoint[0] = closest[0]; closestPoint[1] = closest[1]; closestPoint[2] = closest[2]; 
@@ -446,3 +446,22 @@ int vtkPolyLine::IntersectWithLine(float p1[3], float p2[3],float tol,float& t,
 
   return 0;
 }
+
+int vtkPolyLine::Triangulate(int index, vtkFloatPoints &pts)
+{
+  pts.Reset();
+  for (int subId=0; subId<this->Points.GetNumberOfPoints()-1; subId++)
+    {
+    pts.InsertPoint(subId,this->Points.GetPoint(subId));
+    pts.InsertPoint(subId+1,this->Points.GetPoint(subId+1));
+    }
+
+  return 1;
+}
+
+void vtkPolyLine::Derivatives(int subId, float pcoords[3], float *values, 
+                              int dim, float *derivs)
+{
+
+}
+

src/PolyVert.cc

@@ -126,3 +126,30 @@ int vtkPolyVertex::IntersectWithLine(float p1[3], float p2[3],
 
   return 0;
 }
+
+int vtkPolyVertex::Triangulate(int index, vtkFloatPoints &pts)
+{
+  int subId;
+
+  pts.Reset();
+  for (subId=0; subId<this->Points.GetNumberOfPoints(); subId++)
+    {
+    pts.InsertPoint(subId,this->Points.GetPoint(subId));
+    }
+  return 1;
+}
+
+void vtkPolyVertex::Derivatives(int subId, float pcoords[3], float *values, 
+                            int dim, float *derivs)
+{
+  int i, idx;
+
+  for (i=0; i<dim; i++)
+    {
+    idx = i*dim;
+    derivs[idx] = 0.0;
+    derivs[idx+1] = 0.0;
+    derivs[idx+2] = 0.0;
+    }
+}
+

src/Polygon.cc

@@ -181,25 +181,27 @@ int vtkPolygon::EvaluatePosition(float x[3], float closestPoint[3],
 //
 // If here, point is outside of polygon, so need to find distance to boundary
 //
-  float pc[3], dist2;
-  int ignoreId, numPts;
-  float closest[3];
-  float dummyWeights[MAX_CELL_SIZE];
-
-  numPts = this->Points.GetNumberOfPoints();
-  for (minDist2=LARGE_FLOAT,i=0; i<numPts - 1; i++)
+  else
     {
-    line.Points.SetPoint(0,this->Points.GetPoint(i));
-    line.Points.SetPoint(1,this->Points.GetPoint(i+1));
-    line.EvaluatePosition(x, closest, ignoreId, pc, dist2, dummyWeights);
-    if ( dist2 < minDist2 )
+    float t, dist2;
+    int numPts;
+    float closest[3];
+
+    numPts = this->Points.GetNumberOfPoints();
+    for (minDist2=LARGE_FLOAT,i=0; i<numPts - 1; i++)
       {
-      closestPoint[0] = closest[0]; closestPoint[1] = closest[1]; closestPoint[2] = closest[2];
-      minDist2 = dist2;
+      dist2 = line.DistanceToLine(x,this->Points.GetPoint(i),
+                                  this->Points.GetPoint(i+1),t,closest);
+      if ( dist2 < minDist2 )
+        {
+        closestPoint[0] = closest[0]; 
+        closestPoint[1] = closest[1]; 
+        closestPoint[2] = closest[2];
+        minDist2 = dist2;
+        }
       }
+    return 0;
     }
-
-  return 0;
 }
 
 void vtkPolygon::EvaluateLocation(int& subId, float pcoords[3], float x[3],
@@ -729,7 +731,8 @@ vtkCell *vtkPolygon::GetEdge(int edgeId)
 }
 
 //
-// Compute interpolation weights using 1/(1-r**2) normalized sum.
+// Description:
+// Compute interpolation weights using 1/r**2 normalized sum.
 //
 void vtkPolygon::ComputeWeights(float x[3], float weights[MAX_CELL_SIZE])
 {
@@ -737,20 +740,28 @@ void vtkPolygon::ComputeWeights(float x[3], float weights[MAX_CELL_SIZE])
   int numPts=this->Points.GetNumberOfPoints();
   float maxDist2, sum, *pt;
 
-  for (maxDist2=0.0, i=0; i<numPts; i++)
+  for (sum=0.0, maxDist2=0.0, i=0; i<numPts; i++)
     {
     pt = this->Points.GetPoint(i);
     weights[i] = math.Distance2BetweenPoints(x,pt);
-    if ( weights[i] > maxDist2 ) maxDist2 = weights[i];
+    if ( weights[i] == 0.0 ) //exact hit
+      {
+      for (int j=0; j<numPts; j++) weights[j] = 0.0;
+      weights[i] = 1.0;
+      return;
+      }
+    else
+      {
+      weights[i] = 1.0 / (weights[i]*weights[i]);
+      sum += weights[i];
+      }
     }
 
-  for (sum=0.0, i=0; i<numPts; i++) sum += 1.0 / (1.0 - weights[i]/maxDist2);
-
   for (i=0; i<numPts; i++) weights[i] /= sum;
 }
 
-// 
-// Intersect plane; see whether point is inside.
+//
+// Intersect this plane with finite line defined by p1 & p2 with tolerance tol.
 //
 int vtkPolygon::IntersectWithLine(float p1[3], float p2[3], float tol,float& t,
                                  float x[3], float pcoords[3], int& subId)
@@ -783,3 +794,47 @@ int vtkPolygon::IntersectWithLine(float p1[3], float p2[3], float tol,float& t,
   return 0;
 
 }
+
+int vtkPolygon::Triangulate(int index, vtkFloatPoints &pts)
+{
+  int i, success;
+  float *bounds, d;
+  int verts[MAX_CELL_SIZE];
+  int numVerts=this->PointIds.GetNumberOfIds();
+  static vtkIdList Tris((MAX_CELL_SIZE-2)*3);
+
+  pts.Reset();
+
+  bounds = this->GetBounds();
+  d = sqrt((bounds[1]-bounds[0])*(bounds[1]-bounds[0]) +
+           (bounds[3]-bounds[2])*(bounds[3]-bounds[2]) +
+           (bounds[5]-bounds[4])*(bounds[5]-bounds[4]));
+  Tolerance = TOLERANCE * d;
+  SuccessfulTriangulation = 1;
+  this->ComputeNormal(&this->Points, Normal);
+
+  for (i=0; i<numVerts; i++) verts[i] = i;
+  Tris.Reset();