Skip to content
Find file
Fetching contributors…
Cannot retrieve contributors at this time
2005 lines (1750 sloc) 55.4 KB
/*=========================================================================
Program: Visualization Toolkit
Module: vtkPolygon.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkPolygon.h"
#include "vtkCellArray.h"
#include "vtkDataSet.h"
#include "vtkDoubleArray.h"
#include "vtkLine.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPoints.h"
#include "vtkPriorityQueue.h"
#include "vtkQuad.h"
#include "vtkTriangle.h"
#include "vtkBox.h"
#include "vtkMergePoints.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkSmartPointer.h"
vtkStandardNewMacro(vtkPolygon);
//----------------------------------------------------------------------------
// Instantiate polygon.
vtkPolygon::vtkPolygon()
{
this->Tris = vtkIdList::New();
this->Tris->Allocate(VTK_CELL_SIZE);
this->Triangle = vtkTriangle::New();
this->Quad = vtkQuad::New();
this->TriScalars = vtkDoubleArray::New();
this->TriScalars->Allocate(3);
this->Line = vtkLine::New();
this->Tolerance = 0.0;
this->SuccessfulTriangulation = 0;
this->Normal[0] = this->Normal[1] = this->Normal[2] = 0.0;
this->UseMVCInterpolation = false;
}
//----------------------------------------------------------------------------
vtkPolygon::~vtkPolygon()
{
this->Tris->Delete();
this->Triangle->Delete();
this->Quad->Delete();
this->TriScalars->Delete();
this->Line->Delete();
}
//----------------------------------------------------------------------------
double vtkPolygon::ComputeArea()
{
double normal[3]; //not used, but required for the
//following ComputeArea call
return vtkPolygon::ComputeArea(this->GetPoints(),
this->GetNumberOfPoints(),
this->GetPointIds()->GetPointer(0),
normal);
}
#define VTK_POLYGON_FAILURE -1
#define VTK_POLYGON_OUTSIDE 0
#define VTK_POLYGON_INSIDE 1
#define VTK_POLYGON_INTERSECTION 2
#define VTK_POLYGON_ON_LINE 3
//----------------------------------------------------------------------------
//
// In many of the functions that follow, the Points and PointIds members
// of the Cell are assumed initialized. This is usually done indirectly
// through the GetCell(id) method in the DataSet objects.
//
// Compute the polygon normal from a points list, and a list of point ids
// that index into the points list. This version will handle non-convex
// polygons.
void vtkPolygon::ComputeNormal(vtkPoints *p, int numPts, vtkIdType *pts,
double *n)
{
int i;
double v0[3], v1[3], v2[3];
double ax, ay, az, bx, by, bz;
//
// Check for special triangle case. Saves extra work.
//
n[0] = n[1] = n[2] = 0.0;
if ( numPts == 2 || numPts == 1 )
{
return;
}
if ( numPts == 3 )
{
p->GetPoint(pts[0],v0);
p->GetPoint(pts[1],v1);
p->GetPoint(pts[2],v2);
vtkTriangle::ComputeNormal(v0, v1, v2, n);
return;
}
// Because polygon may be concave, need to accumulate cross products to
// determine true normal.
//
p->GetPoint(pts[0],v1); //set things up for loop
p->GetPoint(pts[1],v2);
for (i=0; i < numPts; i++)
{
v0[0] = v1[0]; v0[1] = v1[1]; v0[2] = v1[2];
v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2];
p->GetPoint(pts[(i+2)%numPts],v2);
// order is important!!! to maintain consistency with polygon vertex order
ax = v2[0] - v1[0]; ay = v2[1] - v1[1]; az = v2[2] - v1[2];
bx = v0[0] - v1[0]; by = v0[1] - v1[1]; bz = v0[2] - v1[2];
n[0] += (ay * bz - az * by);
n[1] += (az * bx - ax * bz);
n[2] += (ax * by - ay * bx);
}
vtkMath::Normalize(n);
}
//----------------------------------------------------------------------------
// Compute the polygon normal from a points list, and a list of point ids
// that index into the points list. This version will handle non-convex
// polygons.
void vtkPolygon::ComputeNormal(vtkIdTypeArray *ids, vtkPoints *p, double n[3])
{
vtkIdType i, numPts = ids->GetNumberOfTuples();
double v0[3], v1[3], v2[3];
double ax, ay, az, bx, by, bz;
//
// Check for special triangle case. Saves extra work.
//
n[0] = n[1] = n[2] = 0.0;
if ( numPts == 2 || numPts == 1 )
{
return;
}
if ( numPts == 3 )
{
p->GetPoint(ids->GetValue(0),v0);
p->GetPoint(ids->GetValue(1),v1);
p->GetPoint(ids->GetValue(2),v2);
vtkTriangle::ComputeNormal(v0, v1, v2, n);
return;
}
// Because polygon may be concave, need to accumulate cross products to
// determine true normal.
//
p->GetPoint(ids->GetValue(0),v1); //set things up for loop
p->GetPoint(ids->GetValue(1),v2);
for (i=0; i < numPts; i++)
{
v0[0] = v1[0]; v0[1] = v1[1]; v0[2] = v1[2];
v1[0] = v2[0]; v1[1] = v2[1]; v1[2] = v2[2];
p->GetPoint(ids->GetValue((i+2)%numPts),v2);
// order is important!!! to maintain consistency with polygon vertex order
ax = v2[0] - v1[0]; ay = v2[1] - v1[1]; az = v2[2] - v1[2];
bx = v0[0] - v1[0]; by = v0[1] - v1[1]; bz = v0[2] - v1[2];
n[0] += (ay * bz - az * by);
n[1] += (az * bx - ax * bz);
n[2] += (ax * by - ay * bx);
}
vtkMath::Normalize(n);
}
//----------------------------------------------------------------------------
// Compute the polygon normal from a list of doubleing points. This version
// will handle non-convex polygons.
void vtkPolygon::ComputeNormal(vtkPoints *p, double *n)
{
int i, numPts;
double v0[3], v1[3], v2[3];
double ax, ay, az, bx, by, bz;
// Polygon is assumed non-convex -> need to accumulate cross products to
// find correct normal.
//
numPts = p->GetNumberOfPoints();
p->GetPoint(0, v1); //set things up for loop
p->GetPoint(1, v2);
n[0] = n[1] = n[2] = 0.0;
for (i=0; i < numPts; i++)
{
memcpy(v0, v1, sizeof(v0));
memcpy(v1, v2, sizeof(v1));
p->GetPoint((i+2)%numPts, v2);
// order is important!!! to maintain consistency with polygon vertex order
ax = v2[0] - v1[0]; ay = v2[1] - v1[1]; az = v2[2] - v1[2];
bx = v0[0] - v1[0]; by = v0[1] - v1[1]; bz = v0[2] - v1[2];
n[0] += (ay * bz - az * by);
n[1] += (az * bx - ax * bz);
n[2] += (ax * by - ay * bx);
}//over all points
vtkMath::Normalize(n);
}
//----------------------------------------------------------------------------
// Compute the polygon normal from an array of points. This version assumes
// that the polygon is convex, and looks for the first valid normal.
void vtkPolygon::ComputeNormal (int numPts, double *pts, double n[3])
{
int i;
double *v1, *v2, *v3;
double length;
double ax, ay, az;
double bx, by, bz;
// Because some polygon vertices are colinear, need to make sure
// first non-zero normal is found.
//
v1 = pts;
v2 = pts + 3;
v3 = pts + 6;
for (i=0; i<numPts-2; i++)
{
ax = v2[0] - v1[0]; ay = v2[1] - v1[1]; az = v2[2] - v1[2];
bx = v3[0] - v1[0]; by = v3[1] - v1[1]; bz = v3[2] - v1[2];
n[0] = (ay * bz - az * by);
n[1] = (az * bx - ax * bz);
n[2] = (ax * by - ay * bx);
length = sqrt (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
if (length != 0.0)
{
n[0] /= length;
n[1] /= length;
n[2] /= length;
return;
}
else
{
v1 = v2;
v2 = v3;
v3 += 3;
}
} //over all points
}
//----------------------------------------------------------------------------
int vtkPolygon::EvaluatePosition(double x[3], double* closestPoint,
int& vtkNotUsed(subId), double pcoords[3],
double& minDist2, double *weights)
{
int i;
double p0[3], p10[3], l10, p20[3], l20, n[3], cp[3];
double ray[3];
this->ParameterizePolygon(p0, p10, l10, p20, l20, n);
this->InterpolateFunctions(x,weights);
vtkPlane::ProjectPoint(x,p0,n,cp);
for (i=0; i<3; i++)
{
ray[i] = cp[i] - p0[i];
}
pcoords[0] = vtkMath::Dot(ray,p10) / (l10*l10);
pcoords[1] = vtkMath::Dot(ray,p20) / (l20*l20);
if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
(this->PointInPolygon(cp, this->Points->GetNumberOfPoints(),
static_cast<vtkDoubleArray *>(
this->Points->GetData())->GetPointer(0),
this->GetBounds(),n) == VTK_POLYGON_INSIDE) )
{
if (closestPoint)
{
closestPoint[0] = cp[0];
closestPoint[1] = cp[1];
closestPoint[2] = cp[2];
minDist2 = vtkMath::Distance2BetweenPoints(x,closestPoint);
}
return 1;
}
// If here, point is outside of polygon, so need to find distance to boundary
//
else
{
double t, dist2;
int numPts;
double closest[3];
double pt1[3], pt2[3];
if (closestPoint)
{
numPts = this->Points->GetNumberOfPoints();
for (minDist2=VTK_DOUBLE_MAX,i=0; i<numPts; i++)
{
this->Points->GetPoint(i, pt1);
this->Points->GetPoint((i+1)%numPts, pt2);
dist2 = vtkLine::DistanceToLine(x, pt1, pt2, t, closest);
if ( dist2 < minDist2 )
{
closestPoint[0] = closest[0];
closestPoint[1] = closest[1];
closestPoint[2] = closest[2];
minDist2 = dist2;
}
}
}
return 0;
}
}
//----------------------------------------------------------------------------
void vtkPolygon::EvaluateLocation(int& vtkNotUsed(subId), double pcoords[3],
double x[3], double *weights)
{
int i;
double p0[3], p10[3], l10, p20[3], l20, n[3];
this->ParameterizePolygon(p0, p10, l10, p20, l20, n);
for (i=0; i<3; i++)
{
x[i] = p0[i] + pcoords[0]*p10[i] + pcoords[1]*p20[i];
}
this->InterpolateFunctions(x,weights);
}
//----------------------------------------------------------------------------
// Compute interpolation weights using 1/r**2 normalized sum or mean value
// coordinate.
void vtkPolygon::InterpolateFunctions(double x[3], double *weights)
{
// Compute interpolation weights using mean value coordinate.
if (this->UseMVCInterpolation)
{
this->InterpolateFunctionsUsingMVC(x, weights);
return;
}
// Compute interpolation weights using 1/r**2 normalized sum.
int i;
int numPts=this->Points->GetNumberOfPoints();
double sum, pt[3];
for (sum=0.0, i=0; i<numPts; i++)
{
this->Points->GetPoint(i, pt);
weights[i] = vtkMath::Distance2BetweenPoints(x,pt);
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 (i=0; i<numPts; i++)
{
weights[i] /= sum;
}
}
//----------------------------------------------------------------------------
// Compute interpolation weights using mean value coordinate.
void vtkPolygon::InterpolateFunctionsUsingMVC(double x[3], double *weights)
{
int numPts=this->Points->GetNumberOfPoints();
// Begin by initializing weights.
for (int i=0; i < numPts; i++)
{
weights[i] = static_cast<double>(0.0);
}
// create local array for storing point-to-vertex vectors and distances
double *dist = new double [numPts];
double *uVec = new double [3*numPts];
static const double eps = 0.00000001;
for (int i=0; i<numPts; i++)
{
double pt[3];
this->Points->GetPoint(i, pt);
// point-to-vertex vector
uVec[3*i] = pt[0] - x[0];
uVec[3*i+1] = pt[1] - x[1];
uVec[3*i+2] = pt[2] - x[2];
// distance
dist[i] = vtkMath::Norm(uVec+3*i);
// handle special case when the point is really close to a vertex
if (dist[i] < eps)
{
weights[i] = 1.0;
delete [] dist;
delete [] uVec;
return;
}
uVec[3*i] /= dist[i];
uVec[3*i+1] /= dist[i];
uVec[3*i+2] /= dist[i];
}
// Now loop over all vertices to compute weight
// w_i = ( tan(theta_i/2) + tan(theta_(i+1)/2) ) / dist_i
// To do consider the simplification of
// tan(alpha/2) = (1-cos(alpha))/sin(alpha)
// = (d0*d1 - cross(u0, u1))/(2*dot(u0,u1))
double *tanHalfTheta = new double [numPts];
for (int i = 0; i < numPts; i++)
{
int i1 = i+1;
if ( i1 == numPts )
{
i1 = 0;
}
double *u0 = uVec + 3*i;
double *u1 = uVec + 3*i1;
double l = sqrt(vtkMath::Distance2BetweenPoints(u0, u1));
double theta = 2.0*asin(l/2.0);
// special case where x lies on an edge
if (vtkMath::Pi() - theta < 0.001)
{
weights[i] = dist[i1] / (dist[i] + dist[i1]);
weights[i1] = 1 - weights[i];
delete [] dist;
delete [] uVec;
delete [] tanHalfTheta;
return;
}
tanHalfTheta[i] = tan(theta/2.0);
}
// Normal case
for (int i = 0; i < numPts; i++)
{
int i1 = i-1;
if ( i1 == -1 )
{
i1 = numPts-1;
}
weights[i] = (tanHalfTheta[i] + tanHalfTheta[i1]) / dist[i];
}
// clear memory
delete [] dist;
delete [] uVec;
delete [] tanHalfTheta;
// normalize weight
double sum = 0.0;
for (int i=0; i < numPts; i++)
{
sum += weights[i];
}
if (fabs(sum) < eps)
{
return;
}
for (int i=0; i < numPts; i++)
{
weights[i] /= sum;
}
return;
}
//----------------------------------------------------------------------------
void vtkPolygon::InterpolateDerivs(double pcoords[3], double *derivs)
{
(void)pcoords;
(void)derivs;
}
//----------------------------------------------------------------------------
// Create a local s-t coordinate system for a polygon. The point p0 is
// the origin of the local system, p10 is s-axis vector, and p20 is the
// t-axis vector. (These are expressed in the modelling coordinate system and
// are vectors of dimension [3].) The values l20 and l20 are the lengths of
// the vectors p10 and p20, and n is the polygon normal.
int vtkPolygon::ParameterizePolygon(double *p0, double *p10, double& l10,
double *p20,double &l20, double *n)
{
int i, j;
double s, t, p[3], p1[3], p2[3], sbounds[2], tbounds[2];
int numPts=this->Points->GetNumberOfPoints();
double x1[3], x2[3];
// This is a two pass process: first create a p' coordinate system
// that is then adjusted to insure that the polygon points are all in
// the range 0<=s,t<=1. The p' system is defined by the polygon normal,
// first vertex and the first edge.
//
this->ComputeNormal (this->Points,n);
this->Points->GetPoint(0, x1);
this->Points->GetPoint(1, x2);
for (i=0; i<3; i++)
{
p0[i] = x1[i];
p10[i] = x2[i] - x1[i];
}
vtkMath::Cross (n,p10,p20);
// Determine lengths of edges
//
if ( (l10=vtkMath::Dot(p10,p10)) == 0.0
|| (l20=vtkMath::Dot(p20,p20)) == 0.0 )
{
return 0;
}
// Now evalute all polygon points to determine min/max parametric
// coordinate values.
//
// first vertex has (s,t) = (0,0)
sbounds[0] = 0.0; sbounds[1] = 0.0;
tbounds[0] = 0.0; tbounds[1] = 0.0;
for(i=1; i<numPts; i++)
{
this->Points->GetPoint(i, x1);
for(j=0; j<3; j++)
{
p[j] = x1[j] - p0[j];
}
#ifdef BAD_WITH_NODEBUG
s = vtkMath::Dot(p,p10) / l10;
t = vtkMath::Dot(p,p20) / l20;
#endif
s = (p[0]*p10[0] + p[1]*p10[1] + p[2]*p10[2]) / l10;
t = (p[0]*p20[0] + p[1]*p20[1] + p[2]*p20[2]) / l20;
sbounds[0] = (s<sbounds[0]?s:sbounds[0]);
sbounds[1] = (s>sbounds[1]?s:sbounds[1]);
tbounds[0] = (t<tbounds[0]?t:tbounds[0]);
tbounds[1] = (t>tbounds[1]?t:tbounds[1]);
}
// Re-evaluate coordinate system
//
for (i=0; i<3; i++)
{
p1[i] = p0[i] + sbounds[1]*p10[i] + tbounds[0]*p20[i];
p2[i] = p0[i] + sbounds[0]*p10[i] + tbounds[1]*p20[i];
p0[i] = p0[i] + sbounds[0]*p10[i] + tbounds[0]*p20[i];
p10[i] = p1[i] - p0[i];
p20[i] = p2[i] - p0[i];
}
l10 = vtkMath::Norm(p10);
l20 = vtkMath::Norm(p20);
return 1;
}
#define VTK_POLYGON_CERTAIN 1
#define VTK_POLYGON_UNCERTAIN 0
#define VTK_POLYGON_RAY_TOL 1.e-03 //Tolerance for ray firing
#define VTK_POLYGON_MAX_ITER 10 //Maximum iterations for ray-firing
#define VTK_POLYGON_VOTE_THRESHOLD 2
#ifndef TRUE
#define FALSE 0
#define TRUE 1
#endif
//----------------------------------------------------------------------------
// Determine whether point is inside polygon. Function uses ray-casting
// to determine if point is inside polygon. Works for arbitrary polygon shape
// (e.g., non-convex). Returns 0 if point is not in polygon; 1 if it is.
// Can also return -1 to indicate degenerate polygon. Note: a point in
// bounding box check is NOT performed prior to in/out check. You may want
// to do this to improve performance.
int vtkPolygon::PointInPolygon (double x[3], int numPts, double *pts,
double bounds[6], double *n)
{
double *x1, *x2, xray[3], u, v;
double rayMag, mag=1, ray[3];
int testResult, rayOK, status, numInts, i;
int iterNumber;
int maxComp, comps[2];
int deltaVotes;
// do a quick bounds check
if ( x[0] < bounds[0] || x[0] > bounds[1] ||
x[1] < bounds[2] || x[1] > bounds[3] ||
x[2] < bounds[4] || x[2] > bounds[5])
{
return VTK_POLYGON_OUTSIDE;
}
//
// Define a ray to fire. The ray is a random ray normal to the
// normal of the face. The length of the ray is a function of the
// size of the face bounding box.
//
for (i=0; i<3; i++)
{
ray[i] = ( bounds[2*i+1] - bounds[2*i] )*1.1 +
fabs((bounds[2*i+1] + bounds[2*i])/2.0 - x[i]);
}
if ( (rayMag = vtkMath::Norm(ray)) == 0.0 )
{
return VTK_POLYGON_OUTSIDE;
}
// Get the maximum component of the normal.
//
if ( fabs(n[0]) > fabs(n[1]) )
{
if ( fabs(n[0]) > fabs(n[2]) )
{
maxComp = 0;
comps[0] = 1;
comps[1] = 2;
}
else
{
maxComp = 2;
comps[0] = 0;
comps[1] = 1;
}
}
else
{
if ( fabs(n[1]) > fabs(n[2]) )
{
maxComp = 1;
comps[0] = 0;
comps[1] = 2;
}
else
{
maxComp = 2;
comps[0] = 0;
comps[1] = 1;
}
}
// Check that max component is non-zero
//
if ( n[maxComp] == 0.0 )
{
return VTK_POLYGON_FAILURE;
}
// Enough information has been acquired to determine the random ray.
// Random rays are generated until one is satisfactory (i.e.,
// produces a ray of non-zero magnitude). Also, since more than one
// ray may need to be fired, the ray-firing occurs in a large loop.
//
// The variable iterNumber counts the number of iterations and is
// limited by the defined variable VTK_POLYGON_MAX_ITER.
//
// The variable deltaVotes keeps track of the number of votes for
// "in" versus "out" of the face. When delta_vote > 0, more votes
// have counted for "in" than "out". When delta_vote < 0, more votes
// have counted for "out" than "in". When the delta_vote exceeds or
// equals the defined variable VTK_POLYGON_VOTE_THRESHOLD, than the
// appropriate "in" or "out" status is returned.
//
for (deltaVotes = 0, iterNumber = 1;
(iterNumber < VTK_POLYGON_MAX_ITER)
&& (abs(deltaVotes) < VTK_POLYGON_VOTE_THRESHOLD);
iterNumber++)
{
//
// Generate ray
//
for (rayOK = FALSE; rayOK == FALSE; )
{
ray[comps[0]] = vtkMath::Random(-rayMag, rayMag);
ray[comps[1]] = vtkMath::Random(-rayMag, rayMag);
ray[maxComp] = -(n[comps[0]]*ray[comps[0]] +
n[comps[1]]*ray[comps[1]]) / n[maxComp];
if ( (mag = vtkMath::Norm(ray)) > rayMag*VTK_TOL )
{
rayOK = TRUE;
}
}
// The ray must be appropriately sized.
//
for (i=0; i<3; i++)
{
xray[i] = x[i] + (rayMag/mag)*ray[i];
}
// The ray may now be fired against all the edges
//
for (numInts=0, testResult=VTK_POLYGON_CERTAIN, i=0; i<numPts; i++)
{
x1 = pts + 3*i;
x2 = pts + 3*((i+1)%numPts);
// Fire the ray and compute the number of intersections. Be careful
// of degenerate cases (e.g., ray intersects at vertex).
//
if ((status=vtkLine::Intersection(x,xray,x1,x2,u,v)) == VTK_POLYGON_INTERSECTION)
{
if ( (VTK_POLYGON_RAY_TOL < v) && (v < 1.0-VTK_POLYGON_RAY_TOL) )
{
numInts++;
}
else
{
testResult = VTK_POLYGON_UNCERTAIN;
}
}
else if ( status == VTK_POLYGON_ON_LINE )
{
testResult = VTK_POLYGON_UNCERTAIN;
}
}
if ( testResult == VTK_POLYGON_CERTAIN )
{
if ( (numInts % 2) == 0)
{
--deltaVotes;
}
else
{
++deltaVotes;
}
}
} //try another ray
// If the number of intersections is odd, the point is in the polygon.
//
if ( deltaVotes < 0 )
{
return VTK_POLYGON_OUTSIDE;
}
else
{
return VTK_POLYGON_INSIDE;
}
}
#define VTK_POLYGON_TOLERANCE 1.0e-06
//----------------------------------------------------------------------------
// Triangulate polygon.
//
int vtkPolygon::Triangulate(vtkIdList *outTris)
{
int success;
double *bounds, d;
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]));
this->Tolerance = VTK_POLYGON_TOLERANCE * d;
this->SuccessfulTriangulation = 1;
this->Tris->Reset();
success = this->EarCutTriangulation();
if ( !success ) //degenerate triangle encountered
{
vtkDebugMacro(<< "Degenerate polygon encountered during triangulation");
}
outTris->DeepCopy(this->Tris);
return success;
}
//----------------------------------------------------------------------------
// Split into non-degenerate polygons prior to triangulation
//
int vtkPolygon::NonDegenerateTriangulate(vtkIdList *outTris)
{
double pt[3], bounds[6];
vtkIdType ptId, numPts;
// ComputeBounds does not give the correct bounds
// So we do it manually
bounds[0]= VTK_DOUBLE_MAX;
bounds[1]=-VTK_DOUBLE_MAX;
bounds[2]= VTK_DOUBLE_MAX;
bounds[3]=-VTK_DOUBLE_MAX;
bounds[4]= VTK_DOUBLE_MAX;
bounds[5]=-VTK_DOUBLE_MAX;
numPts = this->GetNumberOfPoints();
for(int i=0; i<numPts; i++)
{
this->Points->GetPoint(i,pt);
if(pt[0] < bounds[0]) { bounds[0] = pt[0]; }
if(pt[1] < bounds[2]) { bounds[2] = pt[1]; }
if(pt[2] < bounds[4]) { bounds[4] = pt[2]; }
if(pt[0] > bounds[1]) { bounds[1] = pt[0]; }
if(pt[1] > bounds[3]) { bounds[3] = pt[1]; }
if(pt[2] > bounds[5]) { bounds[5] = pt[2]; }
}
outTris->Reset();
outTris->Allocate(3*(2*numPts-4));
vtkPoints *newPts = vtkPoints::New();
newPts->Allocate(numPts);
vtkMergePoints *mergePoints = vtkMergePoints::New();
mergePoints->InitPointInsertion(newPts,bounds);
mergePoints->SetDivisions(10,10,10);
vtkIdTypeArray *matchingIds = vtkIdTypeArray::New();
matchingIds->SetNumberOfTuples(numPts);
int numDuplicatePts = 0;
for(int i=0; i<numPts; i++)
{
this->Points->GetPoint(i,pt);
if(mergePoints->InsertUniquePoint(pt,ptId))
{
matchingIds->SetValue(i,ptId+numDuplicatePts);
}
else
{
matchingIds->SetValue(i,ptId+numDuplicatePts);
numDuplicatePts++;
}
}
mergePoints->Delete();
newPts->Delete();
int numPtsRemoved = 0;
vtkIdType tri[3];
while(numPtsRemoved < numPts)
{
vtkIdType start = 0;
vtkIdType end = numPts-1;
for(;start<numPts;start++)
{
if(matchingIds->GetValue(start) >= 0)
{
break;
}
}
if(start >= end)
{
vtkErrorMacro("ERROR: start >= end");
break;
}
for(int i=start; i<numPts; i++)
{
if(matchingIds->GetValue(i) < 0)
{
continue;
}
if(matchingIds->GetValue(i) != i)
{
start = (matchingIds->GetValue(i) + 1) % numPts;
end = i;
while(matchingIds->GetValue(start) < 0)
{
start++;
}
break;
}
}
vtkPolygon *polygon = vtkPolygon::New();
polygon->Points->SetDataTypeToDouble();
int numPolygonPts = start<end ? end-start+1 : end-start+numPts+1;
for(int i=0; i<numPolygonPts; i++)
{
ptId = start+i;
if(matchingIds->GetValue(ptId) >= 0)
{
numPtsRemoved++;
matchingIds->SetValue(ptId,-1);
polygon->PointIds->InsertNextId(ptId);
polygon->Points->InsertNextPoint(this->Points->GetPoint(ptId));
}
}
vtkIdList *outTriangles = vtkIdList::New();
outTriangles->Allocate(3*(2*polygon->GetNumberOfPoints()-4));
polygon->Triangulate(outTriangles);
int outNumTris = outTriangles->GetNumberOfIds();
for(int i=0; i<outNumTris; i+=3)
{
tri[0] = outTriangles->GetId(i);
tri[1] = outTriangles->GetId(i+1);
tri[2] = outTriangles->GetId(i+2);
tri[0] = polygon->PointIds->GetId(tri[0]);
tri[1] = polygon->PointIds->GetId(tri[1]);
tri[2] = polygon->PointIds->GetId(tri[2]);
outTris->InsertNextId(tri[0]);
outTris->InsertNextId(tri[1]);
outTris->InsertNextId(tri[2]);
}
polygon->Delete();
outTriangles->Delete();
}
matchingIds->Delete();
return 1;
}
//----------------------------------------------------------------------------
// Special structures for building loops. This is a double-linked list.
typedef struct _vtkPolyVertex
{
int id;
double x[3];
double measure;
_vtkPolyVertex* next;
_vtkPolyVertex* previous;
} vtkLocalPolyVertex;
class vtkPolyVertexList { //structure to support triangulation
public:
vtkPolyVertexList(vtkIdList *ptIds, vtkPoints *pts, double tol2);
~vtkPolyVertexList();
int ComputeNormal();
double ComputeMeasure(vtkLocalPolyVertex *vtx);
void RemoveVertex(int i, vtkIdList *, vtkPriorityQueue *);
int CanRemoveVertex(int id, double tol);
int NumberOfVerts;
vtkLocalPolyVertex *Array;
vtkLocalPolyVertex *Head;
double Normal[3];
};
//----------------------------------------------------------------------------
// tolerance is squared
vtkPolyVertexList::vtkPolyVertexList(vtkIdList *ptIds, vtkPoints *pts,
double tol2)
{
int numVerts = ptIds->GetNumberOfIds();
this->NumberOfVerts = numVerts;
this->Array = new vtkLocalPolyVertex [numVerts];
int i;
// now load the data into the array
double x[3];
for (i=0; i<numVerts; i++)
{
this->Array[i].id = i;
pts->GetPoint(i,x);
this->Array[i].x[0] = x[0];
this->Array[i].x[1] = x[1];
this->Array[i].x[2] = x[2];
this->Array[i].next = this->Array + (i+1)%numVerts;
if ( i == 0 )
{
this->Array[i].previous = this->Array + numVerts - 1;
}
else
{
this->Array[i].previous = this->Array + i - 1;
}
}
// Make sure that there are no coincident vertices.
// Beware of multiple coincident vertices.
vtkLocalPolyVertex *vtx, *next;
this->Head = this->Array;
for (vtx=this->Head, i=0; i<numVerts; i++)
{
next = vtx->next;
if ( vtkMath::Distance2BetweenPoints(vtx->x,next->x) < tol2 )
{
next->next->previous = vtx;
vtx->next = next->next;
if ( next == this->Head )
{
this->Head = vtx;
}
this->NumberOfVerts--;
}
else //can move forward
{
vtx = next;
}
}
}
//----------------------------------------------------------------------------
vtkPolyVertexList::~vtkPolyVertexList()
{
delete [] this->Array;
}
//----------------------------------------------------------------------------
// Remove the vertex from the polygon (forming a triangle with
// its previous and next neighbors, and reinsert the neighbors
// into the priority queue.
void vtkPolyVertexList::RemoveVertex(int i, vtkIdList *tris,
vtkPriorityQueue *queue)
{
// Create triangle
tris->InsertNextId(this->Array[i].id);
tris->InsertNextId(this->Array[i].next->id);
tris->InsertNextId(this->Array[i].previous->id);
// remove vertex; special case if single triangle left
if ( --this->NumberOfVerts < 3 )
{
return;
}
if ( (this->Array + i) == this->Head )
{
this->Head = this->Array[i].next;
}
this->Array[i].previous->next = this->Array[i].next;
this->Array[i].next->previous = this->Array[i].previous;
// recompute measure, reinsert into queue
// note that id may have been previously deleted (with Pop()) if we
// are dealing with a concave polygon and vertex couldn't be split.
queue->DeleteId(this->Array[i].previous->id);
queue->DeleteId(this->Array[i].next->id);
if ( this->ComputeMeasure(this->Array[i].previous) > 0.0 )
{
queue->Insert(this->Array[i].previous->measure,
this->Array[i].previous->id);
}
if ( this->ComputeMeasure(this->Array[i].next) > 0.0 )
{
queue->Insert(this->Array[i].next->measure,
this->Array[i].next->id);
}
}
//----------------------------------------------------------------------------
int vtkPolyVertexList::ComputeNormal()
{
vtkLocalPolyVertex *vtx=this->Head;
double v1[3], v2[3], n[3], *anchor=vtx->x;
this->Normal[0] = this->Normal[1] = this->Normal[2] = 0.0;
for (vtx=vtx->next; vtx->next!=this->Head; vtx=vtx->next)
{
v1[0] = vtx->x[0] - anchor[0];
v1[1] = vtx->x[1] - anchor[1];
v1[2] = vtx->x[2] - anchor[2];
v2[0] = vtx->next->x[0] - anchor[0];
v2[1] = vtx->next->x[1] - anchor[1];
v2[2] = vtx->next->x[2] - anchor[2];
vtkMath::Cross(v1,v2,n);
this->Normal[0] += n[0];
this->Normal[1] += n[1];
this->Normal[2] += n[2];
}
if ( vtkMath::Normalize(this->Normal) == 0.0 )
{
return 0;
}
else
{
return 1;
}
}
//----------------------------------------------------------------------------
// The measure is the ratio of triangle perimeter^2 to area;
// the sign of the measure is determined by dotting the local
// vector with the normal (concave features return a negative
// measure).
double vtkPolyVertexList::ComputeMeasure(vtkLocalPolyVertex *vtx)
{
double v1[3], v2[3], v3[3], v4[3], area, perimeter;
for (int i=0; i<3; i++)
{
v1[i] = vtx->x[i] - vtx->previous->x[i];
v2[i] = vtx->next->x[i] - vtx->x[i];
v3[i] = vtx->previous->x[i] - vtx->next->x[i];
}
vtkMath::Cross(v1,v2,v4); //|v4| is twice the area
if ( (area=vtkMath::Dot(v4,this->Normal)) < 0.0 )
{
return (vtx->measure = -1.0); //concave or bad triangle
}
else if ( area == 0.0 )
{
return (vtx->measure = -VTK_DOUBLE_MAX); //concave or bad triangle
}
else
{
perimeter = vtkMath::Norm(v1) + vtkMath::Norm(v2) +
vtkMath::Norm(v3);
return (vtx->measure = perimeter*perimeter/area);
}
}
//----------------------------------------------------------------------------
// returns != 0 if vertex can be removed. Uses half-space
// comparison to determine whether ear-cut is valid, and may
// resort to line-plane intersections to resolve possible
// instersections with ear-cut.
int vtkPolyVertexList::CanRemoveVertex(int id, double tolerance)
{
int i, sign, currentSign;
double v[3], sN[3], *sPt, val, s, t;
vtkLocalPolyVertex *currentVtx, *previous, *next, *vtx;
// Check for simple case
if ( this->NumberOfVerts <= 3 )
{
return 1;
}
// Compute split plane, the point to be cut off
// is always on the positive side of the plane.
currentVtx = this->Array + id;
previous = currentVtx->previous;
next = currentVtx->next;
sPt = previous->x; //point on plane
for (i=0; i<3; i++)
{
v[i] = next->x[i] - previous->x[i]; //vector passing through point
}
vtkMath::Cross (v,this->Normal,sN);
if ( (vtkMath::Normalize(sN)) == 0.0 )
{
return 0; //bad split, indeterminant
}
// Traverse the other points to see if a) they are all on the
// other side of the plane; and if not b) whether they intersect
// the split line.
int oneNegative=0;
val = vtkPlane::Evaluate(sN,sPt,next->next->x);
currentSign = (val > tolerance ? 1 : (val < -tolerance ? -1 : 0));
oneNegative = (currentSign < 0 ? 1 : 0); //very important
// Intersections are only computed when the split half-space is crossed
for (vtx=next->next->next; vtx != previous; vtx = vtx->next)
{
val = vtkPlane::Evaluate(sN,sPt,vtx->x);
sign = (val > tolerance ? 1 : (val < -tolerance ? -1 : 0));
if ( sign != currentSign )
{
if ( !oneNegative )
{
oneNegative = (sign < 0 ? 1 : 0); //very important
}
if (vtkLine::Intersection(sPt,next->x,vtx->x,vtx->previous->x,s,t) != 0 )
{
return 0;
}
else
{
currentSign = sign;
}
}//if crossing occurs
}//for the rest of the loop
if ( !oneNegative )
{
return 0; //entire loop is on this side of plane
}
else
{
return 1;
}
}
//----------------------------------------------------------------------------
// Triangulation method based on ear-cutting. Triangles, or ears, are
// cut off from the polygon based on the angle of the vertex. Small
// angles (narrow triangles) are cut off first. This implementation uses
// a priority queue to cut off ears with smallest angles. Also, the
// algorithm works in 3D (the points don't have to be projected into
// 2D, and the ordering direction of the points is nor important as
// long as the polygon edges do not self intersect).
int vtkPolygon::EarCutTriangulation ()
{
vtkPolyVertexList poly(this->PointIds, this->Points,
this->Tolerance*this->Tolerance);
vtkLocalPolyVertex *vtx;
int i, id;
// First compute the polygon normal the correct way
//
if ( ! poly.ComputeNormal() )
{
return (this->SuccessfulTriangulation=0);
}
// Now compute the angles between edges incident to each
// vertex. Place the structure into a priority queue (those
// vertices with smallest angle are to be removed first).
//
vtkPriorityQueue *VertexQueue = vtkPriorityQueue::New();
VertexQueue->Allocate(poly.NumberOfVerts);
for (i=0, vtx=poly.Head; i < poly.NumberOfVerts; i++, vtx=vtx->next)
{
//concave (negative measure) vertices are not elgible for removal
if ( poly.ComputeMeasure(vtx) > 0.0 )
{
VertexQueue->Insert(vtx->measure, vtx->id);
}
}
// For each vertex in priority queue, and as long as there
// are three or more vertices, remove the vertex (if possible)
// and create a new triangle. If the number of vertices in the
// queue is equal to the number of vertices, then the polygon
// is convex and triangle removal can proceed without intersection
// checks.
//
int numInQueue;
while ( poly.NumberOfVerts > 2 &&
(numInQueue=VertexQueue->GetNumberOfItems()) > 0)
{
if ( numInQueue == poly.NumberOfVerts ) //convex, pop away
{
id = VertexQueue->Pop();
poly.RemoveVertex(id,this->Tris,VertexQueue);
}//convex
else
{
id = VertexQueue->Pop(); //removes it, even if can't be split
if ( poly.CanRemoveVertex(id,this->Tolerance) )
{
poly.RemoveVertex(id,this->Tris,VertexQueue);
}
}//concave
}//while
// Clean up
VertexQueue->Delete();
if ( poly.NumberOfVerts > 2 ) //couldn't triangulate
{
return (this->SuccessfulTriangulation=0);
}
return (this->SuccessfulTriangulation=1);
}
//----------------------------------------------------------------------------
int vtkPolygon::CellBoundary(int vtkNotUsed(subId), double pcoords[3],
vtkIdList *pts)
{
int i, numPts=this->PointIds->GetNumberOfIds();
double x[3], *weights;
int closestPoint=0, previousPoint, nextPoint;
double largestWeight=0.0;
double p0[3], p10[3], l10, p20[3], l20, n[3];
pts->Reset();
weights = new double[numPts];
// determine global coordinates given parametric coordinates
this->ParameterizePolygon(p0, p10, l10, p20, l20, n);
for (i=0; i<3; i++)
{
x[i] = p0[i] + pcoords[0]*p10[i] + pcoords[1]*p20[i];
}
//find edge with largest and next largest weight values. This will be
//the closest edge.
this->InterpolateFunctions(x,weights);
for ( i=0; i < numPts; i++ )
{
if ( weights[i] > largestWeight )
{
closestPoint = i;
largestWeight = weights[i];
}
}
pts->InsertId(0,this->PointIds->GetId(closestPoint));
previousPoint = closestPoint - 1;
nextPoint = closestPoint + 1;
if ( previousPoint < 0 )
{
previousPoint = numPts - 1;
}
if ( nextPoint >= numPts )
{
nextPoint = 0;
}
if ( weights[previousPoint] > weights[nextPoint] )
{
pts->InsertId(1,this->PointIds->GetId(previousPoint));
}
else
{
pts->InsertId(1,this->PointIds->GetId(nextPoint));
}
delete [] weights;
// determine whether point is inside of polygon
if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
(this->PointInPolygon(x, this->Points->GetNumberOfPoints(),
static_cast<vtkDoubleArray *>(
this->Points->GetData())->GetPointer(0),
this->GetBounds(),n) == VTK_POLYGON_INSIDE) )
{
return 1;
}
else
{
return 0;
}
}
//----------------------------------------------------------------------------
void vtkPolygon::Contour(double value, vtkDataArray *cellScalars,
vtkIncrementalPointLocator *locator,
vtkCellArray *verts, vtkCellArray *lines,
vtkCellArray *polys,
vtkPointData *inPd, vtkPointData *outPd,
vtkCellData *inCd, vtkIdType cellId,
vtkCellData *outCd)
{
int i, success;
double *bounds, d;
int p1, p2, p3;
this->TriScalars->SetNumberOfTuples(3);
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]));
this->Tolerance = VTK_POLYGON_TOLERANCE * d;
this->SuccessfulTriangulation = 1;
this->ComputeNormal(this->Points, this->Normal);
this->Tris->Reset();
success = this->EarCutTriangulation();
if ( !success ) // Just skip for now.
{
}
else // Contour triangle
{
for (i=0; i<this->Tris->GetNumberOfIds(); i += 3)
{
p1 = this->Tris->GetId(i);
p2 = this->Tris->GetId(i+1);
p3 = this->Tris->GetId(i+2);
this->Triangle->Points->SetPoint(0,this->Points->GetPoint(p1));
this->Triangle->Points->SetPoint(1,this->Points->GetPoint(p2));
this->Triangle->Points->SetPoint(2,this->Points->GetPoint(p3));
if ( outPd )
{
this->Triangle->PointIds->SetId(0,this->PointIds->GetId(p1));
this->Triangle->PointIds->SetId(1,this->PointIds->GetId(p2));
this->Triangle->PointIds->SetId(2,this->PointIds->GetId(p3));
}
this->TriScalars->SetTuple(0,cellScalars->GetTuple(p1));
this->TriScalars->SetTuple(1,cellScalars->GetTuple(p2));
this->TriScalars->SetTuple(2,cellScalars->GetTuple(p3));
this->Triangle->Contour(value, this->TriScalars, locator, verts,
lines, polys, inPd, outPd, inCd, cellId, outCd);
}
}
}
//----------------------------------------------------------------------------
vtkCell *vtkPolygon::GetEdge(int edgeId)
{
int numPts=this->Points->GetNumberOfPoints();
// load point id's
this->Line->PointIds->SetId(0,this->PointIds->GetId(edgeId));
this->Line->PointIds->SetId(1,this->PointIds->GetId((edgeId+1) % numPts));
// load coordinates
this->Line->Points->SetPoint(0,this->Points->GetPoint(edgeId));
this->Line->Points->SetPoint(1,this->Points->GetPoint((edgeId+1) % numPts));
return this->Line;
}
//----------------------------------------------------------------------------
//
// Intersect this plane with finite line defined by p1 & p2 with tolerance tol.
//
int vtkPolygon::IntersectWithLine(double p1[3], double p2[3], double tol,double& t,
double x[3], double pcoords[3], int& subId)
{
double pt1[3], n[3];
double tol2 = tol*tol;
double closestPoint[3];
double dist2;
int npts = this->GetNumberOfPoints();
double *weights;
subId = 0;
pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
// Define a plane to intersect with
//
this->Points->GetPoint(1, pt1);
this->ComputeNormal (this->Points,n);
// Intersect plane of the polygon with line
//
if ( ! vtkPlane::IntersectWithLine(p1,p2,n,pt1,t,x) )
{
return 0;
}
// Evaluate position
//
weights = new double[npts];
if ( this->EvaluatePosition(x, closestPoint, subId, pcoords, dist2, weights) >= 0)
{
if ( dist2 <= tol2 )
{
delete [] weights;
return 1;
}
}
delete [] weights;
return 0;
}
//----------------------------------------------------------------------------
int vtkPolygon::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
vtkPoints *pts)
{
int i, success;
double *bounds, d;
pts->Reset();
ptIds->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]));
this->Tolerance = VTK_POLYGON_TOLERANCE * d;
this->SuccessfulTriangulation = 1;
this->ComputeNormal(this->Points, this->Normal);
this->Tris->Reset();
success = this->EarCutTriangulation();
if ( !success ) // Indicate possible failure
{
vtkDebugMacro(<<"Possible triangulation failure");
}
for (i=0; i<this->Tris->GetNumberOfIds(); i++)
{
ptIds->InsertId(i,this->PointIds->GetId(this->Tris->GetId(i)));
pts->InsertPoint(i,this->Points->GetPoint(this->Tris->GetId(i)));
}
return this->SuccessfulTriangulation;
}
//----------------------------------------------------------------------------
// Samples at three points to compute derivatives in local r-s coordinate
// system and projects vectors into 3D model coordinate system.
// Note that the results are usually inaccurate because
// this method actually returns the derivative of the interpolation
// function which is obtained using 1/r**2 normalized sum.
#define VTK_SAMPLE_DISTANCE 0.01
void vtkPolygon::Derivatives(int vtkNotUsed(subId), double pcoords[3],
double *values, int dim, double *derivs)
{
int i, j, k, idx;
if ( this->Points->GetNumberOfPoints() == 4 )
{
for(i=0; i<4; i++)
{
this->Quad->Points->SetPoint(i,this->Points->GetPoint(i));
}
this->Quad->Derivatives(0, pcoords, values, dim, derivs);
return;
}
else if ( this->Points->GetNumberOfPoints() == 3 )
{
for(i=0; i<3; i++)
{
this->Triangle->Points->SetPoint(i,this->Points->GetPoint(i));
}
this->Triangle->Derivatives(0, pcoords, values, dim, derivs);
return;
}
double p0[3], p10[3], l10, p20[3], l20, n[3];
double x[3][3], l1, l2, v1[3], v2[3];
int numVerts=this->PointIds->GetNumberOfIds();
double *weights = new double[numVerts];
double *sample = new double[dim*3];
//setup parametric system and check for degeneracy
if ( this->ParameterizePolygon(p0, p10, l10, p20, l20, n) == 0 )
{
for ( j=0; j < dim; j++ )
{
for ( i=0; i < 3; i++ )
{
derivs[j*dim + i] = 0.0;
}
}
return;
}
//compute positions of three sample points
for (i=0; i<3; i++)
{
x[0][i] = p0[i] + pcoords[0]*p10[i] + pcoords[1]*p20[i];
x[1][i] = p0[i] + (pcoords[0]+VTK_SAMPLE_DISTANCE)*p10[i] +
pcoords[1]*p20[i];
x[2][i] = p0[i] + pcoords[0]*p10[i] +
(pcoords[1]+VTK_SAMPLE_DISTANCE)*p20[i];
}
//for each sample point, sample data values
for ( idx=0, k=0; k < 3; k++ ) //loop over three sample points
{
this->InterpolateFunctions(x[k],weights);
for ( j=0; j < dim; j++, idx++) //over number of derivates requested
{
sample[idx] = 0.0;
for ( i=0; i < numVerts; i++ )
{
sample[idx] += weights[i] * values[j + i*dim];
}
}
}
//compute differences along the two axes
for ( i=0; i < 3; i++ )
{
v1[i] = x[1][i] - x[0][i];
v2[i] = x[2][i] - x[0][i];
}
l1 = vtkMath::Normalize(v1);
l2 = vtkMath::Normalize(v2);
//compute derivatives along x-y-z axes
double ddx, ddy;
for ( j=0; j < dim; j++ )
{
ddx = (sample[dim+j] - sample[j]) / l1;
ddy = (sample[2*dim+j] - sample[j]) / l2;
//project onto global x-y-z axes
derivs[3*j] = ddx*v1[0] + ddy*v2[0];
derivs[3*j + 1] = ddx*v1[1] + ddy*v2[1];
derivs[3*j + 2] = ddx*v1[2] + ddy*v2[2];
}
delete [] weights;
delete [] sample;
}
//----------------------------------------------------------------------------
void vtkPolygon::Clip(double value, vtkDataArray *cellScalars,
vtkIncrementalPointLocator *locator, vtkCellArray *tris,
vtkPointData *inPD, vtkPointData *outPD,
vtkCellData *inCD, vtkIdType cellId, vtkCellData *outCD,
int insideOut)
{
int i, success;
double *bounds, d;
int p1, p2, p3;
this->TriScalars->SetNumberOfTuples(3);
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]));
this->Tolerance = VTK_POLYGON_TOLERANCE * d;
this->SuccessfulTriangulation = 1;
this->ComputeNormal(this->Points, this->Normal);
this->Tris->Reset();
success = this->EarCutTriangulation();
if ( success ) // clip triangles
{
for (i=0; i<this->Tris->GetNumberOfIds(); i += 3)
{
p1 = this->Tris->GetId(i);
p2 = this->Tris->GetId(i+1);
p3 = this->Tris->GetId(i+2);
this->Triangle->Points->SetPoint(0,this->Points->GetPoint(p1));
this->Triangle->Points->SetPoint(1,this->Points->GetPoint(p2));
this->Triangle->Points->SetPoint(2,this->Points->GetPoint(p3));
this->Triangle->PointIds->SetId(0,this->PointIds->GetId(p1));
this->Triangle->PointIds->SetId(1,this->PointIds->GetId(p2));
this->Triangle->PointIds->SetId(2,this->PointIds->GetId(p3));
this->TriScalars->SetTuple(0,cellScalars->GetTuple(p1));
this->TriScalars->SetTuple(1,cellScalars->GetTuple(p2));
this->TriScalars->SetTuple(2,cellScalars->GetTuple(p3));
this->Triangle->Clip(value, this->TriScalars, locator, tris,
inPD, outPD, inCD, cellId, outCD, insideOut);
}
}
}
//----------------------------------------------------------------------------
// Method intersects two polygons. You must supply the number of points and
// point coordinates (npts, *pts) and the bounding box (bounds) of the two
// polygons. Also supply a tolerance squared for controlling
// error. The method returns 1 if there is an intersection, and 0 if
// not. A single point of intersection x[3] is also returned if there
// is an intersection.
int vtkPolygon::IntersectPolygonWithPolygon(int npts, double *pts,double bounds[6],
int npts2, double *pts2,
double bounds2[6], double tol2,
double x[3])
{
double n[3], coords[3];
int i, j;
double *p1, *p2, ray[3];
double t;
// Intersect each edge of first polygon against second
//
vtkPolygon::ComputeNormal(npts2, pts2, n);
for (i=0; i<npts; i++)
{
p1 = pts + 3*i;
p2 = pts + 3*((i+1)%npts);
for (j=0; j<3; j++)
{
ray[j] = p2[j] - p1[j];
}
if ( ! vtkBox::IntersectBox(bounds2, p1, ray, coords, t) )
{
continue;
}
if ( (vtkPlane::IntersectWithLine(p1,p2,n,pts2,t,x)) == 1 )
{
if ( (npts2==3
&& vtkTriangle::PointInTriangle(x,pts2,pts2+3,pts2+6,tol2))
|| (npts2>3
&& vtkPolygon::PointInPolygon(x,npts2,pts2,bounds2,n)
==VTK_POLYGON_INSIDE))
{
return 1;
}
}
else
{
return 0;
}
}
// Intersect each edge of second polygon against first
//
vtkPolygon::ComputeNormal(npts, pts, n);
for (i=0; i<npts2; i++)
{
p1 = pts2 + 3*i;
p2 = pts2 + 3*((i+1)%npts2);
for (j=0; j<3; j++)
{
ray[j] = p2[j] - p1[j];
}
if ( ! vtkBox::IntersectBox(bounds, p1, ray, coords, t) )
{
continue;
}
if ( (vtkPlane::IntersectWithLine(p1,p2,n,pts,t,x)) == 1 )
{
if ( (npts==3 && vtkTriangle::PointInTriangle(x,pts,pts+3,pts+6,tol2))
|| (npts>3 && vtkPolygon::PointInPolygon(x,npts,pts,bounds,n)
==VTK_POLYGON_INSIDE))
{
return 1;
}
}
else
{
return 0;
}
}
return 0;
}
//----------------------------------------------------------------------------
// Compute the area of the polygon (oriented in 3D space). It uses an
// efficient approach where the area is computed in 2D and then projected into
// 3D space.
double vtkPolygon::ComputeArea(vtkPoints *p, vtkIdType numPts, vtkIdType *pts,
double n[3])
{
if ( numPts < 3 )
{
return 0.0;
}
else
{
double area=0.0;
double nx, ny, nz;
int coord, i;
vtkPolygon::ComputeNormal(p,numPts, pts, n);
// Select the projection direction
nx = (n[0] > 0.0 ? n[0] : -n[0]); // abs x-coord
ny = (n[1] > 0.0 ? n[1] : -n[1]); // abs y-coord
nz = (n[2] > 0.0 ? n[2] : -n[2]); // abs z-coord
coord = (nx > ny ? (nx > nz ? 0 : 2) : (ny > nz ? 1 : 2));
// compute area of the 2D projection
double x0[3], x1[3], x2[3], *v0, *v1, *v2;
v0 = x0;
v1 = x1;
v2 = x2;
for (i=0; i<numPts; i++)
{
p->GetPoint(pts[i],v0);
p->GetPoint(pts[(i+1)%numPts],v1);
p->GetPoint(pts[(i+2)%numPts],v2);
switch (coord)
{
case 0:
area += v1[1] * (v2[2] - v0[2]);
continue;
case 1:
area += v1[0] * (v2[2] - v0[2]);
continue;
case 2:
area += v1[0] * (v2[1] - v0[1]);
continue;
}
}
// scale to get area before projection
switch (coord)
{
case 0:
area /= (2.0*nx);
break;
case 1:
area /= (2.0*ny);
break;
case 2:
area /= (2.0*nz);
}
return fabs(area);
}//general polygon
}
//----------------------------------------------------------------------------
void vtkPolygon::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Tolerance: " << this->Tolerance << "\n";
os << indent << "SuccessfulTriangulation: " <<
this->SuccessfulTriangulation << "\n";
os << indent << "UseMVCInterpolation: " <<
this->UseMVCInterpolation << "\n";
os << indent << "Normal: (" << this->Normal[0] << ", "
<< this->Normal[1] << ", " << this->Normal[2] << ")\n";
os << indent << "Tris:\n";
this->Tris->PrintSelf(os,indent.GetNextIndent());
os << indent << "Triangle:\n";
this->Triangle->PrintSelf(os,indent.GetNextIndent());
os << indent << "Quad:\n";
this->Quad->PrintSelf(os,indent.GetNextIndent());
os << indent << "TriScalars:\n";
this->TriScalars->PrintSelf(os,indent.GetNextIndent());
os << indent << "Line:\n";
this->Line->PrintSelf(os,indent.GetNextIndent());
}
//----------------------------------------------------------------------------
// Compute the polygon centroid from a points list, and a list of point ids
// that index into the points list.
void vtkPolygon::ComputeCentroid(vtkIdTypeArray *ids,
vtkPoints *p,
double c[3])
{
vtkIdType i, numPts = ids->GetNumberOfTuples();
double p0[3];
double a = 1.0 / static_cast<double>(numPts);
c[0] = c[1] = c[2] = 0.0;
for (i=0; i < numPts; i++)
{
p->GetPoint(ids->GetValue(i),p0);
c[0] += p0[0];
c[1] += p0[1];
c[2] += p0[2];
}
c[0] *= a;
c[1] *= a;
c[2] *= a;
}
//----------------------------------------------------------------------------
double vtkPolygon::DistanceToPolygon(double x[3], int numPts, double *pts,
double bounds[6], double closest[3])
{
// First check to see if the point is inside the polygon
// do a quick bounds check
if ( x[0] >= bounds[0] && x[0] <= bounds[1] &&
x[1] >= bounds[2] && x[1] <= bounds[3] &&
x[2] >= bounds[4] && x[2] <= bounds[5])
{
double n[3];
vtkPolygon::ComputeNormal(numPts, pts, n);
if ( vtkPolygon::PointInPolygon(x,numPts,pts,bounds,n) )
{
closest[0] = x[0];
closest[1] = x[1];
closest[2] = x[2];
return 0.0;
}
}
// Not inside, compute the distance of the point to the edges.
double minDist2=VTK_LARGE_FLOAT;
double *p0, *p1, dist2, t, c[3];
for (int i=0; i<numPts; i++)
{
p0 = pts + 3*i;
p1 = pts + 3*((i+1)%numPts);
dist2 = vtkLine::DistanceToLine(x, p0, p1, t, c);
if ( dist2 < minDist2 )
{
minDist2 = dist2;
closest[0] = c[0];
closest[1] = c[1];
closest[2] = c[2];
}
}
return sqrt(minDist2);
}
//----------------------------------------------------------------------------
int vtkPolygon::IntersectConvex2DCells(vtkCell *cell1, vtkCell *cell2,
double tol, double p0[3], double p1[3])
{
// Intersect the six total edges of the two triangles against each other. Two points are
// all that are required.
double *x[2], pcoords[3], t, x0[3], x1[3];
x[0] = p0; x[1] = p1;
int subId, idx=0;
double t2 = tol*tol;
// Loop over edges of second polygon and intersect against first polygon
vtkIdType i, numPts = cell2->Points->GetNumberOfPoints();
for (i=0; i<numPts; i++)
{
cell2->Points->GetPoint(i,x0);
cell2->Points->GetPoint((i+1)%numPts,x1);
if ( cell1->IntersectWithLine(x0,x1,tol,t,x[idx],pcoords,subId) )
{
if ( idx == 0 )
{
idx++;
}
else if ( ((x[1][0]-x[0][0])*(x[1][0]-x[0][0]) + (x[1][1]-x[0][1])*(x[1][1]-x[0][1]) +
(x[1][2]-x[0][2])*(x[1][2]-x[0][2])) > t2 )
{
return 2;
}
}//if edge intersection
}//over all edges
// Loop over edges of first polygon and intersect against second polygon
numPts = cell1->Points->GetNumberOfPoints();
for (i=0; i<numPts; i++)
{
cell1->Points->GetPoint(i,x0);
cell1->Points->GetPoint((i+1)%numPts,x1);
if ( cell2->IntersectWithLine(x0,x1,tol,t,x[idx],pcoords,subId) )
{
if ( idx == 0 )
{
idx++;
}
else if ( ((x[1][0]-x[0][0])*(x[1][0]-x[0][0]) + (x[1][1]-x[0][1])*(x[1][1]-x[0][1]) +
(x[1][2]-x[0][2])*(x[1][2]-x[0][2])) > t2 )
{
return 2;
}
}//if edge intersection
}//over all edges
// Evaluate what we got
if ( idx == 1 )
{
return 1; //everything intersecting at single point
}
else
{
return 0;
}
}
//----------------------------------------------------------------------------
#ifndef VTK_LEGACY_REMOVE
void vtkPolygon::ComputeWeights(double x[3], double *weights)
{
VTK_LEGACY_REPLACED_BODY(vtkPolygon::ComputeWeights, "VTK 5.2",
vtkPolygon::InterpolateFunctions);
this->InterpolateFunctions(x, weights);
}
#endif
Jump to Line
Something went wrong with that request. Please try again.