You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1277 lines
39 KiB
1277 lines
39 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkDelaunay2D.cxx,v $
|
|
|
|
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 "vtkDelaunay2D.h"
|
|
|
|
#include "vtkAbstractTransform.h"
|
|
#include "vtkCellArray.h"
|
|
#include "vtkDoubleArray.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkMath.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkPlane.h"
|
|
#include "vtkPointData.h"
|
|
#include "vtkPolyData.h"
|
|
#include "vtkPolygon.h"
|
|
#include "vtkStreamingDemandDrivenPipeline.h"
|
|
#include "vtkTriangle.h"
|
|
#include "vtkTransform.h"
|
|
|
|
vtkCxxRevisionMacro(vtkDelaunay2D, "$Revision: 1.68.10.1 $");
|
|
vtkStandardNewMacro(vtkDelaunay2D);
|
|
vtkCxxSetObjectMacro(vtkDelaunay2D,Transform,vtkAbstractTransform);
|
|
|
|
// Construct object with Alpha = 0.0; Tolerance = 0.00001; Offset = 1.25;
|
|
// BoundingTriangulation turned off.
|
|
vtkDelaunay2D::vtkDelaunay2D()
|
|
{
|
|
this->Alpha = 0.0;
|
|
this->Tolerance = 0.00001;
|
|
this->BoundingTriangulation = 0;
|
|
this->Offset = 1.0;
|
|
this->Transform = NULL;
|
|
this->ProjectionPlaneMode = VTK_DELAUNAY_XY_PLANE;
|
|
|
|
// optional 2nd input
|
|
this->SetNumberOfInputPorts(2);
|
|
}
|
|
|
|
vtkDelaunay2D::~vtkDelaunay2D()
|
|
{
|
|
if (this->Transform)
|
|
{
|
|
this->Transform->UnRegister(this);
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Specify the input data or filter. Old style.
|
|
void vtkDelaunay2D::SetSource(vtkPolyData *input)
|
|
{
|
|
this->Superclass::SetInput(1, input);
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Specify the input data or filter. New style.
|
|
void vtkDelaunay2D::SetSourceConnection(vtkAlgorithmOutput *algOutput)
|
|
{
|
|
this->Superclass::SetInputConnection(1, algOutput);
|
|
}
|
|
|
|
vtkPolyData *vtkDelaunay2D::GetSource()
|
|
{
|
|
if (this->GetNumberOfInputConnections(1) < 1)
|
|
{
|
|
return NULL;
|
|
}
|
|
return vtkPolyData::SafeDownCast(
|
|
this->GetExecutive()->GetInputData(1, 0));
|
|
}
|
|
|
|
// Determine whether point x is inside of circumcircle of triangle
|
|
// defined by points (x1, x2, x3). Returns non-zero if inside circle.
|
|
// (Note that z-component is ignored.)
|
|
int vtkDelaunay2D::InCircle (double x[3], double x1[3], double x2[3],
|
|
double x3[3])
|
|
{
|
|
double radius2, center[2], dist2;
|
|
|
|
radius2 = vtkTriangle::Circumcircle(x1,x2,x3,center);
|
|
|
|
// check if inside/outside circumcircle
|
|
dist2 = (x[0]-center[0]) * (x[0]-center[0]) +
|
|
(x[1]-center[1]) * (x[1]-center[1]);
|
|
|
|
if ( dist2 < (0.999999999999*radius2) )
|
|
{
|
|
return 1;
|
|
}
|
|
else
|
|
{
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
#define VTK_DEL2D_TOLERANCE 1.0e-014
|
|
|
|
// Recursive method to locate triangle containing point. Starts with arbitrary
|
|
// triangle (tri) and "walks" towards it. Influenced by some of Guibas and
|
|
// Stolfi's work. Returns id of enclosing triangle, or -1 if no triangle
|
|
// found. Also, the array nei[3] is used to communicate info about points
|
|
// that lie on triangle edges: nei[0] is neighboring triangle id, and nei[1]
|
|
// and nei[2] are the vertices defining the edge.
|
|
vtkIdType vtkDelaunay2D::FindTriangle(double x[3], vtkIdType ptIds[3],
|
|
vtkIdType tri, double tol,
|
|
vtkIdType nei[3], vtkIdList *neighbors)
|
|
{
|
|
int i, j, ir, ic, inside, i2, i3;
|
|
vtkIdType *pts, npts, newNei;
|
|
double p[3][3], n[2], vp[2], vx[2], dp, minProj;
|
|
|
|
// get local triangle info
|
|
this->Mesh->GetCellPoints(tri,npts,pts);
|
|
for (i=0; i<3; i++)
|
|
{
|
|
ptIds[i] = pts[i];
|
|
this->GetPoint(ptIds[i], p[i]);
|
|
}
|
|
|
|
// Randomization (of find edge neighbora) avoids walking in
|
|
// circles in certain weird cases
|
|
srand(tri);
|
|
ir = rand() % 3;
|
|
// evaluate in/out of each edge
|
|
for (inside=1, minProj=0.0, ic=0; ic<3; ic++)
|
|
{
|
|
i = (ir+ic) % 3;
|
|
i2 = (i+1) % 3;
|
|
i3 = (i+2) % 3;
|
|
|
|
// create a 2D edge normal to define a "half-space"; evaluate points (i.e.,
|
|
// candiate point and other triangle vertex not on this edge).
|
|
n[0] = -(p[i2][1] - p[i][1]);
|
|
n[1] = p[i2][0] - p[i][0];
|
|
vtkMath::Normalize2D(n);
|
|
|
|
// compute local vectors
|
|
for (j=0; j<2; j++)
|
|
{
|
|
vp[j] = p[i3][j] - p[i][j];
|
|
vx[j] = x[j] - p[i][j];
|
|
}
|
|
|
|
//check for duplicate point
|
|
vtkMath::Normalize2D(vp);
|
|
if ( vtkMath::Normalize2D(vx) <= tol )
|
|
{
|
|
this->NumberOfDuplicatePoints++;
|
|
return -1;
|
|
}
|
|
|
|
// see if two points are in opposite half spaces
|
|
dp = vtkMath::Dot2D(n,vx) * (vtkMath::Dot2D(n,vp) < 0 ? -1.0 : 1.0);
|
|
if ( dp < VTK_DEL2D_TOLERANCE )
|
|
{
|
|
if ( dp < minProj ) //track edge most orthogonal to point direction
|
|
{
|
|
inside = 0;
|
|
nei[1] = ptIds[i];
|
|
nei[2] = ptIds[i2];
|
|
minProj = dp;
|
|
}
|
|
}//outside this edge
|
|
}//for each edge
|
|
|
|
if ( inside ) // all edges have tested positive
|
|
{
|
|
nei[0] = (-1);
|
|
return tri;
|
|
}
|
|
|
|
else if ( !inside && (fabs(minProj) < VTK_DEL2D_TOLERANCE) ) // on edge
|
|
{
|
|
this->Mesh->GetCellEdgeNeighbors(tri,nei[1],nei[2],neighbors);
|
|
nei[0] = neighbors->GetId(0);
|
|
return tri;
|
|
}
|
|
|
|
else //walk towards point
|
|
{
|
|
this->Mesh->GetCellEdgeNeighbors(tri,nei[1],nei[2],neighbors);
|
|
if ( (newNei=neighbors->GetId(0)) == nei[0] )
|
|
{
|
|
this->NumberOfDegeneracies++;
|
|
return -1;
|
|
}
|
|
else
|
|
{
|
|
nei[0] = tri;
|
|
return this->FindTriangle(x,ptIds,newNei,tol,nei,neighbors);
|
|
}
|
|
}
|
|
}
|
|
|
|
#undef VTK_DEL2D_TOLERANCE
|
|
|
|
// Recursive method checks whether edge is Delaunay, and if not, swaps edge.
|
|
// Continues until all edges are Delaunay. Points p1 and p2 form the edge in
|
|
// question; x is the coordinates of the inserted point; tri is the current
|
|
// triangle id.
|
|
void vtkDelaunay2D::CheckEdge(vtkIdType ptId, double x[3], vtkIdType p1,
|
|
vtkIdType p2, vtkIdType tri)
|
|
{
|
|
int i;
|
|
vtkIdType *pts, npts, numNei, nei, p3;
|
|
double x1[3], x2[3], x3[3];
|
|
vtkIdList *neighbors;
|
|
vtkIdType swapTri[3];
|
|
|
|
this->GetPoint(p1,x1);
|
|
this->GetPoint(p2,x2);
|
|
|
|
neighbors = vtkIdList::New();
|
|
neighbors->Allocate(2);
|
|
|
|
this->Mesh->GetCellEdgeNeighbors(tri,p1,p2,neighbors);
|
|
numNei = neighbors->GetNumberOfIds();
|
|
|
|
if ( numNei > 0 ) //i.e., not a boundary edge
|
|
{
|
|
// get neighbor info including opposite point
|
|
nei = neighbors->GetId(0);
|
|
this->Mesh->GetCellPoints(nei, npts, pts);
|
|
for (i=0; i<2; i++)
|
|
{
|
|
if ( pts[i] != p1 && pts[i] != p2 )
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
p3 = pts[i];
|
|
this->GetPoint(p3,x3);
|
|
|
|
// see whether point is in circumcircle
|
|
if ( this->InCircle (x3, x, x1, x2) )
|
|
{// swap diagonal
|
|
this->Mesh->RemoveReferenceToCell(p1,tri);
|
|
this->Mesh->RemoveReferenceToCell(p2,nei);
|
|
this->Mesh->ResizeCellList(ptId,1);
|
|
this->Mesh->AddReferenceToCell(ptId,nei);
|
|
this->Mesh->ResizeCellList(p3,1);
|
|
this->Mesh->AddReferenceToCell(p3,tri);
|
|
|
|
swapTri[0] = ptId; swapTri[1] = p3; swapTri[2] = p2;
|
|
this->Mesh->ReplaceCell(tri,3,swapTri);
|
|
|
|
swapTri[0] = ptId; swapTri[1] = p1; swapTri[2] = p3;
|
|
this->Mesh->ReplaceCell(nei,3,swapTri);
|
|
|
|
// two new edges become suspect
|
|
this->CheckEdge(ptId, x, p3, p2, tri);
|
|
this->CheckEdge(ptId, x, p1, p3, nei);
|
|
|
|
}//in circle
|
|
}//interior edge
|
|
|
|
neighbors->Delete();
|
|
}
|
|
|
|
// 2D Delaunay triangulation. Steps are as follows:
|
|
// 1. For each point
|
|
// 2. Find triangle point is in
|
|
// 3. Create 3 triangles from each edge of triangle that point is in
|
|
// 4. Recursively evaluate Delaunay criterion for each edge neighbor
|
|
// 5. If criterion not satisfied; swap diagonal
|
|
//
|
|
int vtkDelaunay2D::RequestData(
|
|
vtkInformation *vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
// get the info objects
|
|
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
|
|
vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
|
|
vtkInformation *outInfo = outputVector->GetInformationObject(0);
|
|
|
|
// get the input and ouptut
|
|
vtkPointSet *input = vtkPointSet::SafeDownCast(
|
|
inInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
vtkPolyData *source = 0;
|
|
if (sourceInfo)
|
|
{
|
|
source =
|
|
vtkPolyData::SafeDownCast(sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
}
|
|
vtkPolyData *output = vtkPolyData::SafeDownCast(
|
|
outInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
|
|
vtkIdType numPoints, i;
|
|
vtkIdType numTriangles = 0;
|
|
vtkIdType ptId, tri[4], nei[3];
|
|
vtkIdType p1 = 0;
|
|
vtkIdType p2 = 0;
|
|
vtkPoints *inPoints;
|
|
vtkPoints *points;
|
|
vtkPoints *tPoints = NULL;
|
|
vtkCellArray *triangles;
|
|
int ncells;
|
|
vtkIdType nodes[4][3], *neiPts;
|
|
vtkIdType *triPts = 0;
|
|
vtkIdType numNeiPts;
|
|
vtkIdType npts = 0;
|
|
vtkIdType pts[3];
|
|
vtkIdList *neighbors, *cells;
|
|
double center[3], radius, tol, x[3];
|
|
int *triUse = NULL;
|
|
double *bounds;
|
|
|
|
vtkDebugMacro(<<"Generating 2D Delaunay triangulation");
|
|
|
|
if (this->Transform && this->BoundingTriangulation)
|
|
{
|
|
vtkWarningMacro(<<"Bounding triangulation cannot be used when an input transform is specified. Output will not contain bounding triangulation.");
|
|
}
|
|
|
|
if (this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE && this->BoundingTriangulation)
|
|
{
|
|
vtkWarningMacro(<<"Bounding triangulation cannot be used when the best fitting plane option is on. Output will not contain bounding triangulation.");
|
|
}
|
|
|
|
// Initialize; check input
|
|
//
|
|
if ( (inPoints=input->GetPoints()) == NULL )
|
|
{
|
|
vtkErrorMacro("<<Cannot triangulate; no input points");
|
|
return 1;
|
|
}
|
|
|
|
if ( (numPoints=inPoints->GetNumberOfPoints()) <= 2 )
|
|
{
|
|
vtkErrorMacro("<<Cannot triangulate; need at least 3 input points");
|
|
return 1;
|
|
}
|
|
|
|
neighbors = vtkIdList::New(); neighbors->Allocate(2);
|
|
cells = vtkIdList::New(); cells->Allocate(64);
|
|
|
|
this->NumberOfDuplicatePoints = 0;
|
|
this->NumberOfDegeneracies = 0;
|
|
|
|
this->Mesh = vtkPolyData::New();
|
|
|
|
|
|
// If the user specified a transform, apply it to the input data.
|
|
//
|
|
// Only the input points are transformed. We do not bother
|
|
// transforming the source points (if specified). The reason is
|
|
// that only the topology of the Source is used during the constrain
|
|
// operation. The point ids in the Source topology are assumed to
|
|
// reference points in the input. So, when an input transform is
|
|
// used, only the input points are transformed. We do not bother
|
|
// with transforming the Source points since they are never
|
|
// referenced.
|
|
if (this->Transform)
|
|
{
|
|
tPoints = vtkPoints::New();
|
|
this->Transform->TransformPoints(inPoints, tPoints);
|
|
}
|
|
else
|
|
{
|
|
// If the user asked this filter to compute the best fitting plane,
|
|
// proceed to compute the plane and generate a transform that will
|
|
// map the input points into that plane.
|
|
if(this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE)
|
|
{
|
|
this->SetTransform( this->ComputeBestFittingPlane(input) );
|
|
tPoints = vtkPoints::New();
|
|
this->Transform->TransformPoints(inPoints, tPoints);
|
|
}
|
|
}
|
|
|
|
// Create initial bounding triangulation. Have to create bounding points.
|
|
// Initialize mesh structure.
|
|
//
|
|
points = vtkPoints::New();
|
|
// This will copy doubles to doubles if the input is double.
|
|
points->SetDataTypeToDouble();
|
|
points->SetNumberOfPoints(numPoints);
|
|
if (!this->Transform)
|
|
{
|
|
points->DeepCopy(inPoints);
|
|
}
|
|
else
|
|
{
|
|
points->DeepCopy(tPoints);
|
|
tPoints->Delete();
|
|
tPoints = NULL;
|
|
}
|
|
|
|
bounds = points->GetBounds();
|
|
center[0] = (bounds[0]+bounds[1])/2.0;
|
|
center[1] = (bounds[2]+bounds[3])/2.0;
|
|
center[2] = (bounds[4]+bounds[5])/2.0;
|
|
tol = input->GetLength();
|
|
radius = this->Offset * tol;
|
|
tol *= this->Tolerance;
|
|
|
|
for (ptId=0; ptId<8; ptId++)
|
|
{
|
|
x[0] = center[0]
|
|
+ radius*cos((double)(45.0*ptId)*vtkMath::DegreesToRadians());
|
|
x[1] = center[1]
|
|
+ radius*sin((double)(45.0*ptId)*vtkMath::DegreesToRadians());
|
|
x[2] = center[2];
|
|
points->InsertPoint(numPoints+ptId,x);
|
|
}
|
|
// We do this for speed accessing points
|
|
this->Points = ((vtkDoubleArray *)points->GetData())->GetPointer(0);
|
|
|
|
triangles = vtkCellArray::New();
|
|
triangles->Allocate(triangles->EstimateSize(2*numPoints,3));
|
|
|
|
//create bounding triangles (there are six)
|
|
pts[0] = numPoints; pts[1] = numPoints + 1; pts[2] = numPoints + 2;
|
|
triangles->InsertNextCell(3,pts);
|
|
pts[0] = numPoints + 2; pts[1] = numPoints + 3; pts[2] = numPoints + 4;
|
|
triangles->InsertNextCell(3,pts);
|
|
pts[0] = numPoints + 4; pts[1] = numPoints + 5; pts[2] = numPoints + 6;
|
|
triangles->InsertNextCell(3,pts);
|
|
pts[0] = numPoints + 6; pts[1] = numPoints + 7; pts[2] = numPoints + 0;
|
|
triangles->InsertNextCell(3,pts);
|
|
pts[0] = numPoints + 0; pts[1] = numPoints + 2; pts[2] = numPoints + 6;
|
|
triangles->InsertNextCell(3,pts);
|
|
pts[0] = numPoints + 2; pts[1] = numPoints + 4; pts[2] = numPoints + 6;
|
|
triangles->InsertNextCell(3,pts);
|
|
tri[0] = 0;
|
|
|
|
this->Mesh->SetPoints(points);
|
|
this->Mesh->SetPolys(triangles);
|
|
this->Mesh->BuildLinks(); //build cell structure
|
|
|
|
// For each point; find triangle containing point. Then evaluate three
|
|
// neighboring triangles for Delaunay criterion. Triangles that do not
|
|
// satisfy criterion have their edges swapped. This continues recursively
|
|
// until all triangles have been shown to be Delaunay.
|
|
//
|
|
for (ptId=0; ptId < numPoints; ptId++)
|
|
{
|
|
this->GetPoint(ptId,x);
|
|
nei[0] = (-1); //where we are coming from...nowhere initially
|
|
|
|
if ( (tri[0] = this->FindTriangle(x,pts,tri[0],tol,nei,neighbors)) >= 0 )
|
|
{
|
|
if ( nei[0] < 0 ) //in triangle
|
|
{
|
|
//delete this triangle; create three new triangles
|
|
//first triangle is replaced with one of the new ones
|
|
nodes[0][0] = ptId; nodes[0][1] = pts[0]; nodes[0][2] = pts[1];
|
|
this->Mesh->RemoveReferenceToCell(pts[2], tri[0]);
|
|
this->Mesh->ReplaceCell(tri[0], 3, nodes[0]);
|
|
this->Mesh->ResizeCellList(ptId,1);
|
|
this->Mesh->AddReferenceToCell(ptId,tri[0]);
|
|
|
|
//create two new triangles
|
|
nodes[1][0] = ptId; nodes[1][1] = pts[1]; nodes[1][2] = pts[2];
|
|
tri[1] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[1]);
|
|
|
|
nodes[2][0] = ptId; nodes[2][1] = pts[2]; nodes[2][2] = pts[0];
|
|
tri[2] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[2]);
|
|
|
|
// Check edge neighbors for Delaunay criterion. If not satisfied, flip
|
|
// edge diagonal. (This is done recursively.)
|
|
this->CheckEdge(ptId, x, pts[0], pts[1], tri[0]);
|
|
this->CheckEdge(ptId, x, pts[1], pts[2], tri[1]);
|
|
this->CheckEdge(ptId, x, pts[2], pts[0], tri[2]);
|
|
}
|
|
|
|
else // on triangle edge
|
|
{
|
|
//update cell list
|
|
this->Mesh->GetCellPoints(nei[0],numNeiPts,neiPts);
|
|
for (i=0; i<3; i++)
|
|
{
|
|
if ( neiPts[i] != nei[1] && neiPts[i] != nei[2] )
|
|
{
|
|
p1 = neiPts[i];
|
|
}
|
|
if ( pts[i] != nei[1] && pts[i] != nei[2] )
|
|
{
|
|
p2 = pts[i];
|
|
}
|
|
}
|
|
this->Mesh->ResizeCellList(p1,1);
|
|
this->Mesh->ResizeCellList(p2,1);
|
|
|
|
//replace two triangles
|
|
this->Mesh->RemoveReferenceToCell(nei[2],tri[0]);
|
|
this->Mesh->RemoveReferenceToCell(nei[2],nei[0]);
|
|
nodes[0][0] = ptId; nodes[0][1] = p2; nodes[0][2] = nei[1];
|
|
this->Mesh->ReplaceCell(tri[0], 3, nodes[0]);
|
|
nodes[1][0] = ptId; nodes[1][1] = p1; nodes[1][2] = nei[1];
|
|
this->Mesh->ReplaceCell(nei[0], 3, nodes[1]);
|
|
this->Mesh->ResizeCellList(ptId, 2);
|
|
this->Mesh->AddReferenceToCell(ptId,tri[0]);
|
|
this->Mesh->AddReferenceToCell(ptId,nei[0]);
|
|
|
|
tri[1] = nei[0];
|
|
|
|
//create two new triangles
|
|
nodes[2][0] = ptId; nodes[2][1] = p2; nodes[2][2] = nei[2];
|
|
tri[2] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[2]);
|
|
|
|
nodes[3][0] = ptId; nodes[3][1] = p1; nodes[3][2] = nei[2];
|
|
tri[3] = this->Mesh->InsertNextLinkedCell(VTK_TRIANGLE, 3, nodes[3]);
|
|
|
|
// Check edge neighbors for Delaunay criterion.
|
|
for ( i=0; i<4; i++ )
|
|
{
|
|
this->CheckEdge (ptId, x, nodes[i][1], nodes[i][2], tri[i]);
|
|
}
|
|
}
|
|
}//if triangle found
|
|
|
|
else
|
|
{
|
|
tri[0] = 0; //no triangle found
|
|
}
|
|
|
|
if ( ! (ptId % 1000) )
|
|
{
|
|
vtkDebugMacro(<<"point #" << ptId);
|
|
this->UpdateProgress ((double)ptId/numPoints);
|
|
if (this->GetAbortExecute())
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
|
|
}//for all points
|
|
|
|
vtkDebugMacro(<<"Triangulated " << numPoints <<" points, "
|
|
<< this->NumberOfDuplicatePoints
|
|
<< " of which were duplicates");
|
|
|
|
if ( this->NumberOfDegeneracies > 0 )
|
|
{
|
|
vtkDebugMacro(<< this->NumberOfDegeneracies
|
|
<< " degenerate triangles encountered, mesh quality suspect");
|
|
}
|
|
|
|
// Finish up by recovering the boundary, or deleting all triangles connected
|
|
// to the bounding triangulation points or not satisfying alpha criterion,
|
|
if ( !this->BoundingTriangulation || this->Alpha > 0.0 || source )
|
|
{
|
|
numTriangles = this->Mesh->GetNumberOfCells();
|
|
if ( source )
|
|
{
|
|
triUse = this->RecoverBoundary(source);
|
|
}
|
|
else
|
|
{
|
|
triUse = new int[numTriangles];
|
|
for (i=0; i<numTriangles; i++)
|
|
{
|
|
triUse[i] = 1;
|
|
}
|
|
}
|
|
}
|
|
|
|
// Delete triangles connected to boundary points (if not desired)
|
|
if ( ! this->BoundingTriangulation )
|
|
{
|
|
for (ptId=numPoints; ptId < (numPoints+8); ptId++)
|
|
{
|
|
this->Mesh->GetPointCells(ptId, cells);
|
|
ncells = cells->GetNumberOfIds();
|
|
for (i=0; i < ncells; i++)
|
|
{
|
|
triUse[cells->GetId(i)] = 0; //mark as deleted
|
|
}
|
|
}
|
|
}
|
|
|
|
// If non-zero alpha value, then figure out which parts of mesh are
|
|
// contained within alpha radius.
|
|
//
|
|
if ( this->Alpha > 0.0 )
|
|
{
|
|
double alpha2 = this->Alpha * this->Alpha;
|
|
double x1[3], x2[3], x3[3];
|
|
double xx1[3], xx2[3], xx3[3];
|
|
vtkIdType cellId, numNei, ap1, ap2, neighbor;
|
|
|
|
vtkCellArray *alphaVerts = vtkCellArray::New();
|
|
alphaVerts->Allocate(numPoints);
|
|
vtkCellArray *alphaLines = vtkCellArray::New();
|
|
alphaLines->Allocate(numPoints);
|
|
|
|
char *pointUse = new char[numPoints+8];
|
|
for (ptId=0; ptId < (numPoints+8); ptId++)
|
|
{
|
|
pointUse[ptId] = 0;
|
|
}
|
|
|
|
//traverse all triangles; evaluating Delaunay criterion
|
|
for (i=0; i < numTriangles; i++)
|
|
{
|
|
if ( triUse[i] == 1 )
|
|
{
|
|
this->Mesh->GetCellPoints(i, npts, triPts);
|
|
|
|
// if any point is one of the bounding points that was added
|
|
// at the beginning of the algorithm, then grab the points
|
|
// from the variable "points" (this list has the boundary
|
|
// points and the original points have been transformed by the
|
|
// input transform). if none of the points are bounding points,
|
|
// then grab the points from the variable "inPoints" so the alpha
|
|
// criterion is applied in the nontransformed space.
|
|
if (triPts[0]<numPoints && triPts[1]<numPoints && triPts[2]<numPoints)
|
|
{
|
|
inPoints->GetPoint(triPts[0],x1);
|
|
inPoints->GetPoint(triPts[1],x2);
|
|
inPoints->GetPoint(triPts[2],x3);
|
|
}
|
|
else
|
|
{
|
|
points->GetPoint(triPts[0],x1);
|
|
points->GetPoint(triPts[1],x2);
|
|
points->GetPoint(triPts[2],x3);
|
|
}
|
|
|
|
// evaluate the alpha criterion in 3D
|
|
vtkTriangle::ProjectTo2D(x1, x2, x3, xx1, xx2, xx3);
|
|
if ( vtkTriangle::Circumcircle(xx1,xx2,xx3,center) > alpha2 )
|
|
{
|
|
triUse[i] = 0;
|
|
}
|
|
else
|
|
{
|
|
for (int j=0; j<3; j++)
|
|
{
|
|
pointUse[triPts[j]] = 1;
|
|
}
|
|
}
|
|
}//if non-deleted triangle
|
|
}//for all triangles
|
|
|
|
//traverse all edges see whether we need to create some
|
|
for (cellId=0, triangles->InitTraversal();
|
|
triangles->GetNextCell(npts,triPts); cellId++)
|
|
{
|
|
if ( ! triUse[cellId] )
|
|
{
|
|
for (i=0; i < npts; i++)
|
|
{
|
|
ap1 = triPts[i];
|
|
ap2 = triPts[(i+1)%npts];
|
|
|
|
if (this->BoundingTriangulation || (ap1<numPoints && ap2<numPoints))
|
|
{
|
|
this->Mesh->GetCellEdgeNeighbors(cellId,ap1,ap2,neighbors);
|
|
numNei = neighbors->GetNumberOfIds();
|
|
|
|
if ( numNei < 1 || ((neighbor=neighbors->GetId(0)) > cellId
|
|
&& !triUse[neighbor]) )
|
|
{//see whether edge is shorter than Alpha
|
|
|
|
// same argument as above, if one is a boundary point, get
|
|
// it using this->GetPoint() which are transformed points. if
|
|
// neither of the points are boundary points, get the from
|
|
// inPoints (untransformed points) so alpha comparison done
|
|
// untransformed space
|
|
if (ap1 < numPoints && ap2 < numPoints)
|
|
{
|
|
inPoints->GetPoint(ap1,x1);
|
|
inPoints->GetPoint(ap2,x2);
|
|
}
|
|
else
|
|
{
|
|
this->GetPoint(ap1,x1);
|
|
this->GetPoint(ap2,x2);
|
|
}
|
|
if ( (vtkMath::Distance2BetweenPoints(x1,x2)*0.25) <= alpha2 )
|
|
{
|
|
pointUse[ap1] = 1; pointUse[ap2] = 1;
|
|
pts[0] = ap1;
|
|
pts[1] = ap2;
|
|
alphaLines->InsertNextCell(2,pts);
|
|
}//if passed test
|
|
}//test edge
|
|
}//if valid edge
|
|
}//for all edges of this triangle
|
|
}//if triangle not output
|
|
}//for all triangles
|
|
|
|
//traverse all points, create vertices if none used
|
|
for (ptId=0; ptId<(numPoints+8); ptId++)
|
|
{
|
|
if ( !pointUse[ptId]
|
|
&& (ptId < numPoints || this->BoundingTriangulation) )
|
|
{
|
|
pts[0] = ptId;
|
|
alphaVerts->InsertNextCell(1,pts);
|
|
}
|
|
}
|
|
|
|
// update output
|
|
delete [] pointUse;
|
|
output->SetVerts(alphaVerts);
|
|
alphaVerts->Delete();
|
|
output->SetLines(alphaLines);
|
|
alphaLines->Delete();
|
|
}
|
|
|
|
// Update output; free up supporting data structures.
|
|
//
|
|
if ( this->BoundingTriangulation && !this->Transform)
|
|
{
|
|
output->SetPoints(points);
|
|
}
|
|
else
|
|
{
|
|
output->SetPoints(inPoints);
|
|
output->GetPointData()->PassData(input->GetPointData());
|
|
}
|
|
|
|
if ( this->Alpha <= 0.0 && this->BoundingTriangulation && !source )
|
|
{
|
|
output->SetPolys(triangles);
|
|
}
|
|
else
|
|
{
|
|
vtkCellArray *alphaTriangles = vtkCellArray::New();
|
|
alphaTriangles->Allocate(numTriangles);
|
|
vtkIdType *alphaTriPts;
|
|
|
|
for (i=0; i<numTriangles; i++)
|
|
{
|
|
if ( triUse[i] )
|
|
{
|
|
this->Mesh->GetCellPoints(i,npts,alphaTriPts);
|
|
alphaTriangles->InsertNextCell(3,alphaTriPts);
|
|
}
|
|
}
|
|
output->SetPolys(alphaTriangles);
|
|
alphaTriangles->Delete();
|
|
delete [] triUse;
|
|
}
|
|
|
|
points->Delete();
|
|
triangles->Delete();
|
|
this->Mesh->Delete();
|
|
neighbors->Delete();
|
|
cells->Delete();
|
|
|
|
// If the best fitting option was ON, then the current transform
|
|
// is the one that was computed internally. We must now destroy it.
|
|
if( this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE )
|
|
{
|
|
if (this->Transform)
|
|
{
|
|
this->Transform->UnRegister(this);
|
|
this->Transform = NULL;
|
|
}
|
|
}
|
|
|
|
output->Squeeze();
|
|
|
|
return 1;
|
|
}
|
|
|
|
// Methods used to recover edges. Uses lines and polygons to determine boundary
|
|
// and inside/outside.
|
|
//
|
|
// Only the topology of the Source is used during the constrain operation.
|
|
// The point ids in the Source topology are assumed to reference points
|
|
// in the input. So, when an input transform is used, only the input points
|
|
// are transformed. We do not bother with transforming the Source points
|
|
// since they are never referenced.
|
|
int *vtkDelaunay2D::RecoverBoundary(vtkPolyData *source)
|
|
{
|
|
vtkCellArray *lines=source->GetLines();
|
|
vtkCellArray *polys=source->GetPolys();
|
|
vtkIdType *pts = 0;
|
|
vtkIdType npts = 0;
|
|
vtkIdType i, p1, p2;
|
|
int *triUse;
|
|
|
|
// Recover the edges of the mesh
|
|
for ( lines->InitTraversal(); lines->GetNextCell(npts,pts); )
|
|
{
|
|
for (i=0; i<(npts-1); i++)
|
|
{
|
|
p1 = pts[i];
|
|
p2 = pts[i+1];
|
|
if ( ! this->Mesh->IsEdge(p1,p2) )
|
|
{
|
|
this->RecoverEdge(p1, p2);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Recover the enclosed regions (polygons) of the mesh
|
|
for ( polys->InitTraversal(); polys->GetNextCell(npts,pts); )
|
|
{
|
|
for (i=0; i<npts; i++)
|
|
{
|
|
p1 = pts[i];
|
|
p2 = pts[(i+1)%npts];
|
|
if ( ! this->Mesh->IsEdge(p1,p2) )
|
|
{
|
|
this->RecoverEdge(p1, p2);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Generate inside/outside marks on mesh
|
|
int numTriangles = this->Mesh->GetNumberOfCells();
|
|
triUse = new int[numTriangles];
|
|
for (i=0; i<numTriangles; i++)
|
|
{
|
|
triUse[i] = 1;
|
|
}
|
|
|
|
// Use any polygons to mark inside and outside. (Note that if an edge was not
|
|
// recovered, we're going to have a problem.) The first polygon is assumed to
|
|
// define the outside of the polygon; additional polygons carve out inside
|
|
// holes.
|
|
this->FillPolygons(polys, triUse);
|
|
|
|
return triUse;
|
|
}
|
|
|
|
// Method attempts to recover an edge by retriangulating mesh around the edge.
|
|
// What we do is identify a "submesh" of triangles that includes the edge to recover.
|
|
// Then we split the submesh in two with the recovered edge, and triangulate each of
|
|
// the two halves. If any part of this fails, we leave things alone.
|
|
int vtkDelaunay2D::RecoverEdge(vtkIdType p1, vtkIdType p2)
|
|
{
|
|
vtkIdType cellId = 0;
|
|
int i, j, k;
|
|
double p1X[3], p2X[3], xyNormal[3], splitNormal[3], p21[3];
|
|
double x1[3], x2[3], sepNormal[3], v21[3];
|
|
int ncells, v1=0, v2=0, signX1=0, signX2, signP1, signP2;
|
|
vtkIdType *pts, *leftTris, *rightTris, npts, numRightTris, numLeftTris;
|
|
int success=0;
|
|
|
|
vtkIdList *cells=vtkIdList::New(); cells->Allocate(64);
|
|
vtkIdList *tris=vtkIdList::New(); tris->Allocate(64);
|
|
vtkPolygon *rightPoly=vtkPolygon::New();
|
|
vtkPolygon *leftPoly=vtkPolygon::New();
|
|
vtkIdList *leftChain=leftPoly->GetPointIds();
|
|
vtkIdList *rightChain=rightPoly->GetPointIds();
|
|
vtkPoints *leftChainX=leftPoly->GetPoints();
|
|
vtkPoints *rightChainX=rightPoly->GetPoints();
|
|
vtkIdList *neis=vtkIdList::New(); neis->Allocate(4);
|
|
vtkIdList *rightPtIds=vtkIdList::New(); rightPtIds->Allocate(64);
|
|
vtkIdList *leftPtIds=vtkIdList::New(); leftPtIds->Allocate(64);
|
|
vtkPoints *rightTriPts=vtkPoints::New(); rightTriPts->Allocate(64);
|
|
vtkPoints *leftTriPts=vtkPoints::New(); leftTriPts->Allocate(64);
|
|
|
|
// Compute a split plane along (p1,p2) and parallel to the z-axis.
|
|
//
|
|
this->GetPoint(p1,p1X); p1X[2] = 0.0; //split plane point
|
|
this->GetPoint(p2,p2X); p2X[2] = 0.0; //split plane point
|
|
xyNormal[0] = xyNormal[1] = 0.0; xyNormal[2] = 1.0;
|
|
for (i=0; i<3; i++ )
|
|
{
|
|
p21[i] = p2X[i] - p1X[i]; //working in x-y plane
|
|
}
|
|
|
|
vtkMath::Cross (p21,xyNormal,splitNormal);
|
|
if ( vtkMath::Normalize(splitNormal) == 0.0 )
|
|
{//Usually means coincident points
|
|
goto FAILURE;
|
|
}
|
|
|
|
// Identify a triangle connected to the point p1 containing a portion of the edge.
|
|
//
|
|
this->Mesh->GetPointCells(p1, cells);
|
|
ncells = cells->GetNumberOfIds();
|
|
for (i=0; i < ncells; i++)
|
|
{
|
|
cellId = cells->GetId(i);
|
|
this->Mesh->GetCellPoints(cellId, npts, pts);
|
|
for (j=0; j<3; j++)
|
|
{
|
|
if ( pts[j] == p1 )
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
v1 = pts[(j+1)%3];
|
|
v2 = pts[(j+2)%3];
|
|
this->GetPoint(v1,x1); x1[2] = 0.0;
|
|
this->GetPoint(v2,x2); x2[2] = 0.0;
|
|
signX1 = (vtkPlane::Evaluate(splitNormal, p1X, x1) > 0.0 ? 1 : -1);
|
|
signX2 = (vtkPlane::Evaluate(splitNormal, p1X, x2) > 0.0 ? 1 : -1);
|
|
if ( signX1 != signX2 ) //points of triangle on either side of edge
|
|
{
|
|
// determine if edge separates p1 from p2 - then we've found triangle
|
|
v21[0] = x2[0] - x1[0]; //working in x-y plane
|
|
v21[1] = x2[1] - x1[1];
|
|
v21[2] = 0.0;
|
|
|
|
vtkMath::Cross (v21,xyNormal,sepNormal);
|
|
if ( vtkMath::Normalize(sepNormal) == 0.0 )
|
|
{ //bad mesh
|
|
goto FAILURE;
|
|
}
|
|
|
|
signP1 = (vtkPlane::Evaluate(sepNormal, x1, p1X) > 0.0 ? 1 : -1);
|
|
signP2 = (vtkPlane::Evaluate(sepNormal, x1, p2X) > 0.0 ? 1 : -1);
|
|
if ( signP1 != signP2 ) //is a separation line
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
} //for all cells
|
|
|
|
if ( i >= ncells )
|
|
{//something is really screwed up
|
|
goto FAILURE;
|
|
}
|
|
|
|
// We found initial triangle; begin to track triangles containing edge. Also,
|
|
// the triangle defines the beginning of two "chains" which form a boundary of
|
|
// enclosing triangles around the edge. Create the two chains (from p1 to p2).
|
|
// (The chains are actually defining two polygons on either side of the edge.)
|
|
//
|
|
tris->InsertId(0, cellId);
|
|
rightChain->InsertId(0, p1); rightChainX->InsertPoint(0, p1X);
|
|
leftChain->InsertId(0, p1); leftChainX->InsertPoint(0, p1X);
|
|
if ( signX1 > 0 )
|
|
{
|
|
rightChain->InsertId(1, v1); rightChainX->InsertPoint(1, x1);
|
|
leftChain->InsertId(1, v2); leftChainX->InsertPoint(1, x2);
|
|
}
|
|
else
|
|
{
|
|
leftChain->InsertId(1, v1); leftChainX->InsertPoint(1, x1);
|
|
rightChain->InsertId(1, v2); rightChainX->InsertPoint(1, x2);
|
|
}
|
|
|
|
// Walk along triangles (edge neighbors) towards point p2.
|
|
while ( v1 != p2 )
|
|
{
|
|
this->Mesh->GetCellEdgeNeighbors(cellId, v1, v2, neis);
|
|
if ( neis->GetNumberOfIds() != 1 )
|
|
{//Mesh is folded or degenerate
|
|
goto FAILURE;
|
|
}
|
|
cellId = neis->GetId(0);
|
|
tris->InsertNextId(cellId);
|
|
this->Mesh->GetCellPoints(cellId, npts, pts);
|
|
for (j=0; j<3; j++)
|
|
{
|
|
if ( pts[j] != v1 && pts[j] != v2 )
|
|
{//found point opposite current edge (v1,v2)
|
|
if ( pts[j] == p2 )
|
|
{
|
|
v1 = p2; //this will cause the walk to stop
|
|
rightChain->InsertNextId(p2); rightChainX->InsertNextPoint(p2X);
|
|
leftChain->InsertNextId(p2); leftChainX->InsertNextPoint(p2X);
|
|
}
|
|
else
|
|
{//keep walking
|
|
this->GetPoint(pts[j], x1); x1[2] = 0.0;
|
|
if ( vtkPlane::Evaluate(splitNormal, p1X, x1) > 0.0 )
|
|
{
|
|
v1 = pts[j];
|
|
rightChain->InsertNextId(v1); rightChainX->InsertNextPoint(x1);
|
|
}
|
|
else
|
|
{
|
|
v2 = pts[j];
|
|
leftChain->InsertNextId(v2); leftChainX->InsertNextPoint(x1);
|
|
}
|
|
}
|
|
break;
|
|
}//else found opposite point
|
|
}//for all points in triangle
|
|
}//while walking
|
|
|
|
// Now that the to chains are formed, each chain forms a polygon (along with
|
|
// the edge (p1,p2)) that requires triangulation. If we can successfully
|
|
// triangulate the two polygons, we will delete the triangles contained within
|
|
// the chains and replace them with the new triangulation.
|
|
//
|
|
success = 1;
|
|
success &= (rightPoly->Triangulate(0, rightPtIds, rightTriPts));
|
|
numRightTris = rightPtIds->GetNumberOfIds() / 3;
|
|
|
|
success &= (leftPoly->Triangulate(0, leftPtIds, leftTriPts));
|
|
numLeftTris = leftPtIds->GetNumberOfIds() / 3;
|
|
|
|
if ( ! success )
|
|
{//polygons on either side of edge are poorly shaped
|
|
goto FAILURE;
|
|
}
|
|
|
|
// Okay, delete the old triangles and replace them with new ones. There should be
|
|
// the same number of new triangles as old ones.
|
|
leftTris = leftPtIds->GetPointer(0);
|
|
for ( j=i=0; i<numLeftTris; i++, j++, leftTris+=3)
|
|
{
|
|
cellId = tris->GetId(j);
|
|
this->Mesh->RemoveCellReference(cellId);
|
|
for (k=0; k<3; k++)
|
|
{//allocate new space for cell lists
|
|
this->Mesh->ResizeCellList(leftTris[k],1);
|
|
}
|
|
this->Mesh->ReplaceLinkedCell(cellId, 3, leftTris);
|
|
}
|
|
|
|
rightTris = rightPtIds->GetPointer(0);
|
|
for ( i=0; i<numRightTris; i++, j++, rightTris+=3)
|
|
{
|
|
cellId = tris->GetId(j);
|
|
this->Mesh->RemoveCellReference(cellId);
|
|
for (k=0; k<3; k++)
|
|
{//allocate new space for cell lists
|
|
this->Mesh->ResizeCellList(rightTris[k],1);
|
|
}
|
|
this->Mesh->ReplaceLinkedCell(cellId, 3, rightTris);
|
|
}
|
|
|
|
FAILURE:
|
|
tris->Delete(); cells->Delete();
|
|
leftPoly->Delete(); rightPoly->Delete(); neis->Delete();
|
|
rightPtIds->Delete(); leftPtIds->Delete();
|
|
rightTriPts->Delete(); leftTriPts->Delete();
|
|
return success;
|
|
}
|
|
|
|
void vtkDelaunay2D::FillPolygons(vtkCellArray *polys, int *triUse)
|
|
{
|
|
vtkIdType p1, p2, j, kk;
|
|
int i, k;
|
|
vtkIdType *pts = 0;
|
|
vtkIdType *triPts;
|
|
vtkIdType npts = 0;
|
|
vtkIdType numPts;
|
|
static double xyNormal[3]={0.0,0.0,1.0};
|
|
double negDir[3], x21[3], x1[3], x2[3], x[3];
|
|
vtkIdList *neis=vtkIdList::New();
|
|
vtkIdType cellId, numNeis;
|
|
vtkIdList *currentFront = vtkIdList::New(), *tmpFront;
|
|
vtkIdList *nextFront = vtkIdList::New();
|
|
vtkIdType numCellsInFront, neiId;
|
|
vtkIdType numTriangles=this->Mesh->GetNumberOfCells();
|
|
|
|
// Loop over edges of polygon, marking triangles on "outside" of polygon as outside.
|
|
// Then perform a fill.
|
|
for ( polys->InitTraversal(); polys->GetNextCell(npts,pts); )
|
|
{
|
|
currentFront->Reset();
|
|
for (i=0; i<npts; i++)
|
|
{
|
|
p1 = pts[i];
|
|
p2 = pts[(i+1)%npts];
|
|
if ( ! this->Mesh->IsEdge(p1,p2) )
|
|
{
|
|
vtkWarningMacro(<<"Edge not recovered, polygon fill suspect");
|
|
}
|
|
else //Mark the "outside" triangles
|
|
{
|
|
neis->Reset();
|
|
this->GetPoint(p1,x1);
|
|
this->GetPoint(p2,x2);
|
|
for (j=0; j<3; j++)
|
|
{
|
|
x21[j] = x2[j] - x1[j];
|
|
}
|
|
vtkMath::Cross (x21,xyNormal,negDir);
|
|
this->Mesh->GetCellEdgeNeighbors(-1, p1, p2, neis); //get both triangles
|
|
numNeis = neis->GetNumberOfIds();
|
|
for (j=0; j<numNeis; j++)
|
|
{//find the vertex not on the edge; evaluate it (and the cell) in/out
|
|
cellId = neis->GetId(j);
|
|
this->Mesh->GetCellPoints(cellId, numPts, triPts);
|
|
for (k=0; k<3; k++)
|
|
{
|
|
if ( triPts[k] != p1 && triPts[k] != p2 )
|
|
{
|
|
break;
|
|
}
|
|
}
|
|
this->GetPoint(triPts[k],x); x[2] = 0.0;
|
|
if ( vtkPlane::Evaluate(negDir, x1, x) > 0.0 )
|
|
{
|
|
triUse[cellId] = 0;
|
|
currentFront->InsertNextId(cellId);
|
|
}
|
|
else
|
|
{
|
|
triUse[cellId] = -1;
|
|
}
|
|
}
|
|
}//edge was recovered
|
|
}//for all edges in polygon
|
|
|
|
// Okay, now perform a fill operation (filling "outside" values).
|
|
//
|
|
while ( (numCellsInFront = currentFront->GetNumberOfIds()) > 0 )
|
|
{
|
|
for (j=0; j < numCellsInFront; j++)
|
|
{
|
|
cellId = currentFront->GetId(j);
|
|
|
|
this->Mesh->GetCellPoints(cellId, numPts, triPts);
|
|
for (k=0; k<3; k++)
|
|
{
|
|
p1 = triPts[k];
|
|
p2 = triPts[(k+1)%3];
|
|
|
|
this->Mesh->GetCellEdgeNeighbors(cellId, p1, p2, neis);
|
|
numNeis = neis->GetNumberOfIds();
|
|
for (kk=0; kk<numNeis; kk++)
|
|
{
|
|
neiId = neis->GetId(kk);
|
|
if ( triUse[neiId] == 1 ) //0 is what we're filling with
|
|
{
|
|
triUse[neiId] = 0;
|
|
nextFront->InsertNextId(neiId);
|
|
}
|
|
}//mark all neigbors
|
|
}//for all edges of cell
|
|
} //all cells in front
|
|
|
|
tmpFront = currentFront;
|
|
currentFront = nextFront;
|
|
nextFront = tmpFront;
|
|
nextFront->Reset();
|
|
} //while still advancing
|
|
|
|
}//for all polygons
|
|
|
|
//convert all unvisited to inside
|
|
for (i=0; i<numTriangles; i++)
|
|
{
|
|
if ( triUse[i] == -1 )
|
|
{
|
|
triUse[i] = 1;
|
|
}
|
|
}
|
|
|
|
currentFront->Delete();
|
|
nextFront->Delete();
|
|
neis->Delete();
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int vtkDelaunay2D::FillInputPortInformation(int port, vtkInformation* info)
|
|
{
|
|
if (port == 0)
|
|
{
|
|
info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPointSet");
|
|
}
|
|
else if (port == 1)
|
|
{
|
|
info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
|
|
info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
|
|
}
|
|
return 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
vtkAbstractTransform * vtkDelaunay2D::ComputeBestFittingPlane(
|
|
vtkPointSet *input)
|
|
{
|
|
vtkIdType numPts=input->GetNumberOfPoints();
|
|
double m[9], v[3], x[3];
|
|
vtkIdType ptId;
|
|
int i;
|
|
double *c1, *c2, *c3, det;
|
|
double normal[3];
|
|
double origin[3];
|
|
|
|
const double tolerance = 1.0e-03;
|
|
|
|
// This code was taken from the vtkTextureMapToPlane class
|
|
// and slightly modified.
|
|
//
|
|
for (i=0; i<3; i++)
|
|
{
|
|
normal[i] = 0.0;
|
|
}
|
|
|
|
// Compute least squares approximation.
|
|
// Compute 3x3 least squares matrix
|
|
v[0] = v[1] = v[2] = 0.0;
|
|
for (i=0; i<9; i++)
|
|
{
|
|
m[i] = 0.0;
|
|
}
|
|
|
|
for (ptId=0; ptId < numPts; ptId++)
|
|
{
|
|
input->GetPoint(ptId, x);
|
|
|
|
v[0] += x[0]*x[2];
|
|
v[1] += x[1]*x[2];
|
|
v[2] += x[2];
|
|
|
|
m[0] += x[0]*x[0];
|
|
m[1] += x[0]*x[1];
|
|
m[2] += x[0];
|
|
|
|
m[3] += x[0]*x[1];
|
|
m[4] += x[1]*x[1];
|
|
m[5] += x[1];
|
|
|
|
m[6] += x[0];
|
|
m[7] += x[1];
|
|
}
|
|
m[8] = numPts;
|
|
|
|
origin[0] = m[2] / numPts;
|
|
origin[1] = m[5] / numPts;
|
|
origin[2] = v[2] / numPts;
|
|
|
|
// Solve linear system using Kramers rule
|
|
//
|
|
c1 = m; c2 = m+3; c3 = m+6;
|
|
if ( (det = vtkMath::Determinant3x3 (c1,c2,c3)) > tolerance )
|
|
{
|
|
normal[0] = vtkMath::Determinant3x3 (v,c2,c3) / det;
|
|
normal[1] = vtkMath::Determinant3x3 (c1,v,c3) / det;
|
|
normal[2] = -1.0; // because of the formulation
|
|
}
|
|
|
|
vtkTransform * transform = vtkTransform::New();
|
|
|
|
// Set the new Z axis as the normal to the best fitting
|
|
// plane.
|
|
double zaxis[3];
|
|
zaxis[0] = 0;
|
|
zaxis[1] = 0;
|
|
zaxis[2] = 1;
|
|
|
|
double rotationAxis[3];
|
|
|
|
vtkMath::Normalize(normal);
|
|
vtkMath::Cross(normal,zaxis,rotationAxis);
|
|
vtkMath::Normalize(rotationAxis);
|
|
|
|
const double rotationAngle =
|
|
180.0*acos(vtkMath::Dot(zaxis,normal))/3.1415926;
|
|
|
|
transform->PreMultiply();
|
|
transform->Identity();
|
|
|
|
transform->RotateWXYZ(rotationAngle,
|
|
rotationAxis[0], rotationAxis[1], rotationAxis[2]);
|
|
|
|
// Set the center of mass as the origin of coordinates
|
|
transform->Translate( -origin[0], -origin[1], -origin[2] );
|
|
|
|
return transform;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkDelaunay2D::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
os << indent << "Alpha: " << this->Alpha << "\n";
|
|
os << indent << "ProjectionPlaneMode: "
|
|
<< ((this->ProjectionPlaneMode == VTK_BEST_FITTING_PLANE)? "Best Fitting Plane" : "XY Plane") << "\n";
|
|
os << indent << "Transform: " << (this->Transform ? "specified" : "none") << "\n";
|
|
os << indent << "Tolerance: " << this->Tolerance << "\n";
|
|
os << indent << "Offset: " << this->Offset << "\n";
|
|
os << indent << "Bounding Triangulation: "
|
|
<< (this->BoundingTriangulation ? "On\n" : "Off\n");
|
|
}
|
|
|