Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
434 lines (377 sloc) 14 KB
/*=========================================================================
Program: Visualization Toolkit
Module: vtkConvexPointSet.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 "vtkConvexPointSet.h"
#include "vtkCellArray.h"
#include "vtkDoubleArray.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkOrderedTriangulator.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkTetra.h"
#include "vtkTriangle.h"
vtkStandardNewMacro(vtkConvexPointSet);
//----------------------------------------------------------------------------
// Construct the hexahedron with eight points.
vtkConvexPointSet::vtkConvexPointSet()
{
this->Tetra = vtkTetra::New();
this->TetraIds = vtkIdList::New();
this->TetraPoints = vtkPoints::New();
this->TetraScalars = vtkDoubleArray::New();
this->TetraScalars->SetNumberOfTuples(4);
this->BoundaryTris = vtkCellArray::New();
this->BoundaryTris->Allocate(100);
this->Triangle = vtkTriangle::New();
this->Triangulator = vtkOrderedTriangulator::New();
this->Triangulator->PreSortedOff();
this->Triangulator->UseTemplatesOff();
this->ParametricCoords = NULL;
}
//----------------------------------------------------------------------------
vtkConvexPointSet::~vtkConvexPointSet()
{
this->Tetra->Delete();
this->TetraIds->Delete();
this->TetraPoints->Delete();
this->TetraScalars->Delete();
this->BoundaryTris->Delete();
this->Triangle->Delete();
if( this->ParametricCoords )
{
this->ParametricCoords->Delete();
}
}
//----------------------------------------------------------------------------
// Should be called by GetCell() prior to any other method invocation
void vtkConvexPointSet::Initialize()
{
// Initialize
vtkIdType numPts = this->GetNumberOfPoints();
if ( numPts < 1 ) return;
this->Triangulate(0, this->TetraIds,this->TetraPoints);
}
//----------------------------------------------------------------------------
int vtkConvexPointSet::GetNumberOfFaces()
{
this->BoundaryTris->Reset();
this->Triangulator->AddTriangles(this->BoundaryTris);
return this->BoundaryTris->GetNumberOfCells();
}
//----------------------------------------------------------------------------
vtkCell *vtkConvexPointSet::GetFace(int faceId)
{
int numCells = this->BoundaryTris->GetNumberOfCells();
if ( faceId < 0 || faceId >=numCells ) {return NULL;}
vtkIdType *cells = this->BoundaryTris->GetPointer();
// Each triangle has three points plus number of points
vtkIdType *cptr = cells + 4*faceId;
for (int i=0; i<3; i++)
{
this->Triangle->PointIds->SetId(i,this->PointIds->GetId(cptr[i+1]));
this->Triangle->Points->SetPoint(i,this->Points->GetPoint(cptr[i+1]));
}
return this->Triangle;
}
//----------------------------------------------------------------------------
int vtkConvexPointSet::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
vtkPoints *pts)
{
int numPts=this->GetNumberOfPoints();
int i;
double x[3];
vtkIdType ptId;
// Initialize
ptIds->Reset();
pts->Reset();
if ( numPts < 1 )
{
return 0;
}
// Initialize Delaunay insertion process.
// No more than numPts points can be inserted.
this->Triangulator->InitTriangulation(this->GetBounds(), numPts);
// Inject cell points into triangulation. Recall that the PreSortedOff()
// flag was set which means that the triangulator will order the points
// according to point id. We insert points with id == the index into the
// vtkConvexPointSet::PointIds and Points; but sort on the global point
// id.
for (i=0; i<numPts; i++)
{
ptId = this->PointIds->GetId(i);
this->Points->GetPoint(i, x);
this->Triangulator->InsertPoint(i, ptId, x, x, 0);
}//for all points
// triangulate the points
this->Triangulator->Triangulate();
// Add the triangulation to the mesh
this->Triangulator->AddTetras(0,ptIds,pts);
return 1;
}
//----------------------------------------------------------------------------
void vtkConvexPointSet::Contour(double value,
vtkDataArray *cellScalars,
vtkIncrementalPointLocator *locator,
vtkCellArray *verts, vtkCellArray *lines,
vtkCellArray *polys,
vtkPointData *inPd, vtkPointData *outPd,
vtkCellData *inCd, vtkIdType cellId,
vtkCellData *outCd)
{
// For each tetra, contour it
int i, j;
vtkIdType ptId, localId;
int numTets = this->TetraIds->GetNumberOfIds() / 4;
for (i=0; i<numTets; i++)
{
for (j=0; j<4; j++)
{
localId = this->TetraIds->GetId(4*i+j);
ptId = this->PointIds->GetId(localId);
this->Tetra->PointIds->SetId(j,ptId);
this->Tetra->Points->SetPoint(j,this->TetraPoints->GetPoint(4*i+j));
this->TetraScalars->SetValue(j,cellScalars->GetTuple1(localId));
}
this->Tetra->Contour(value,this->TetraScalars,locator,verts,lines,polys,
inPd,outPd,inCd,cellId,outCd);
}
}
//----------------------------------------------------------------------------
void vtkConvexPointSet::Clip(double value,
vtkDataArray *cellScalars,
vtkIncrementalPointLocator *locator, vtkCellArray *tets,
vtkPointData *inPD, vtkPointData *outPD,
vtkCellData *inCD, vtkIdType cellId,
vtkCellData *outCD, int insideOut)
{
// For each tetra, contour it
int i, j;
vtkIdType ptId, localId;
int numTets = this->TetraIds->GetNumberOfIds() / 4;
for (i=0; i<numTets; i++)
{
for (j=0; j<4; j++)
{
localId = this->TetraIds->GetId(4*i+j);
ptId = this->PointIds->GetId(localId);
this->Tetra->PointIds->SetId(j,ptId);
this->Tetra->Points->SetPoint(j,this->TetraPoints->GetPoint(4*i+j));
this->TetraScalars->SetValue(j,cellScalars->GetTuple1(localId));
}
this->Tetra->Clip(value,this->TetraScalars,locator,tets,inPD,outPD,inCD,
cellId, outCD, insideOut);
}
}
//----------------------------------------------------------------------------
int vtkConvexPointSet::CellBoundary(int subId, double pcoords[3],
vtkIdList *pts)
{
// This function was re-written to fix bug #9550.
// Thanks go to Bart Janssens.
vtkIdType pntIndx;
for ( int i = 0; i < 4; i ++ )
{
pntIndx = this->PointIds
->GetId( this->TetraIds->GetId( (subId << 2) + i ) );
this->Tetra->PointIds->SetId( i, pntIndx );
this->Tetra->Points
->SetPoint( i, this->TetraPoints->GetPoint( (subId << 2) + i ) );
}
// find the parametrically nearest triangle.
return this->Tetra->CellBoundary( subId, pcoords, pts );
}
//----------------------------------------------------------------------------
int vtkConvexPointSet::EvaluatePosition( double x[3],
double * vtkNotUsed(closestPoint),
int & subId, double pcoords[3],
double & minDist2, double * weights )
{
double pc[3], dist2;
int ignoreId, i, j, k, returnStatus=0, status;
double tempWeights[4];
double closest[3];
vtkIdType ptId;
int numPnts = this->GetNumberOfPoints();
int numTets = this->TetraIds->GetNumberOfIds() >> 2;
for ( minDist2 = VTK_DOUBLE_MAX, i = 0; i < numTets; i ++ )
{
for ( j = 0; j < 4; j ++ )
{
ptId = this->PointIds->GetId( this->TetraIds->GetId( (i << 2) + j ) );
this->Tetra->PointIds->SetId( j, ptId );
this->Tetra->Points->SetPoint
( j, this->TetraPoints->GetPoint( (i << 2) + j ) );
}
status = this->Tetra->EvaluatePosition( x, closest, ignoreId, pc, dist2,
tempWeights );
if ( status != -1 && dist2 < minDist2 )
{
// init (clear) all the weights since only the vertices of the closest
// tetrahedron are assigned with valid weights while the rest vertices
// (of those farther tetrahedra) are simply inited with zero weights
// (to make no any contribution). This fixes bug #9453
for ( k = 0; k < numPnts; k ++ )
{
weights[k] = 0.0;
}
returnStatus = status;
minDist2 = dist2;
subId = i;
pcoords[0] = pc[0];
pcoords[1] = pc[1];
pcoords[2] = pc[2];
// assign valid weights to the vertices of this closest tetrahedron only
// This fixes bug #9453.
weights[ this->TetraIds->GetId( (i << 2) ) ] = tempWeights[0];
weights[ this->TetraIds->GetId( (i << 2) + 1 ) ] = tempWeights[1];
weights[ this->TetraIds->GetId( (i << 2) + 2 ) ] = tempWeights[2];
weights[ this->TetraIds->GetId( (i << 2) + 3 ) ] = tempWeights[3];
}
}
return returnStatus;
}
//----------------------------------------------------------------------------
void vtkConvexPointSet::EvaluateLocation( int & subId, double pcoords[3],
double x[3], double * weights )
{
int i;
int numPnts;
double tmpWgts[4];
vtkIdType pntIndx;
for ( i = 0; i < 4; i ++ )
{
pntIndx = this->PointIds
->GetId( this->TetraIds->GetId( (subId << 2) + i ) );
this->Tetra->PointIds->SetId( i, pntIndx );
this->Tetra->Points
->SetPoint( i, this->TetraPoints->GetPoint( (subId << 2) + i ) );
}
// use tmpWgts to collect the valid weights of the tetra's four vertices
this->Tetra->EvaluateLocation( subId, pcoords, x, tmpWgts );
// init the actual array of weights (possibly greater than 4)
numPnts = this->GetNumberOfPoints();
for ( i = 0; i < numPnts; i ++ )
{
weights[i] = 0.0;
}
// update the target weights only
weights[ this->TetraIds->GetId( (subId << 2) ) ] = tmpWgts[0];
weights[ this->TetraIds->GetId( (subId << 2) + 1 ) ] = tmpWgts[1];
weights[ this->TetraIds->GetId( (subId << 2) + 2 ) ] = tmpWgts[2];
weights[ this->TetraIds->GetId( (subId << 2) + 3 ) ] = tmpWgts[3];
}
//----------------------------------------------------------------------------
int vtkConvexPointSet::IntersectWithLine(double p1[3], double p2[3], double tol,
double& minT, double x[3],
double pcoords[3], int& subId)
{
int subTest, i, j;
vtkIdType ptId;
double t, pc[3], xTemp[3];
int numTets = this->TetraIds->GetNumberOfIds() / 4;
int status = 0;
for (minT=VTK_DOUBLE_MAX, i=0; i<numTets; i++)
{
for (j=0; j<4; j++)
{
ptId = this->PointIds->GetId(this->TetraIds->GetId(4*i+j));
this->Tetra->PointIds->SetId(j,ptId);
this->Tetra->Points->SetPoint(j,this->TetraPoints->GetPoint(4*i+j));
}
if (this->Tetra->IntersectWithLine(p1,p2,tol,t,xTemp,pc,subTest) &&
t < minT )
{
status = 1;
subId = i;
minT = t;
x[0] = xTemp[0];
x[1] = xTemp[1];
x[2] = xTemp[2];
pcoords[0] = pc[0];
pcoords[1] = pc[1];
pcoords[2] = pc[2];
}
}
return status;
}
//----------------------------------------------------------------------------
void vtkConvexPointSet::Derivatives(int subId, double pcoords[3],
double *values, int dim, double *derivs)
{
vtkIdType ptId;
for (int j=0; j<4; j++)
{
ptId = this->PointIds->GetId(this->TetraIds->GetId(4*subId+j));
this->Tetra->PointIds->SetId(j,ptId);
this->Tetra->Points->SetPoint(j,this->TetraPoints->GetPoint(4*subId+j));
}
this->Tetra->Derivatives(subId, pcoords, values, dim, derivs);
}
//----------------------------------------------------------------------------
double *vtkConvexPointSet::GetParametricCoords()
{
int numPts = this->PointIds->GetNumberOfIds();
if ( ! this->ParametricCoords )
{
this->ParametricCoords = vtkDoubleArray::New();
}
this->ParametricCoords->SetNumberOfComponents(3);
this->ParametricCoords->SetNumberOfTuples(numPts);
double p[3], x[3], *bounds = this->GetBounds();
int i, j;
for (i=0; i < numPts; i++)
{
this->Points->GetPoint(i, x);
for (j=0; j<3; j++)
{
p[j] = (x[j] - bounds[2*j]) / (bounds[2*j+1] - bounds[2*j]);
}
this->ParametricCoords->SetTuple(i,p);
}
return this->ParametricCoords->GetPointer(0);
}
//----------------------------------------------------------------------------
void vtkConvexPointSet::InterpolateFunctions(double pcoords[3], double *sf)
{
(void)pcoords;
(void)sf;
}
//----------------------------------------------------------------------------
void vtkConvexPointSet::InterpolateDerivs(double pcoords[3], double *derivs)
{
(void)pcoords;
(void)derivs;
}
//----------------------------------------------------------------------------
void vtkConvexPointSet::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Tetra:\n";
this->Tetra->PrintSelf(os,indent.GetNextIndent());
os << indent << "TetraIds:\n";
this->TetraIds->PrintSelf(os,indent.GetNextIndent());
os << indent << "TetraPoints:\n";
this->TetraPoints->PrintSelf(os,indent.GetNextIndent());
os << indent << "TetraScalars:\n";
this->TetraScalars->PrintSelf(os,indent.GetNextIndent());
os << indent << "BoundaryTris:\n";
this->BoundaryTris->PrintSelf(os,indent.GetNextIndent());
os << indent << "Triangle:\n";
this->Triangle->PrintSelf(os,indent.GetNextIndent());
if ( this->ParametricCoords )
{
os << indent << "ParametricCoords " << this->ParametricCoords << "\n";
}
else
{
os << indent << "ParametricCoords: (null)\n";
}
}