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.
1096 lines
30 KiB
1096 lines
30 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkTriangle.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 "vtkTriangle.h"
|
|
|
|
#include "vtkCellArray.h"
|
|
#include "vtkCellData.h"
|
|
#include "vtkLine.h"
|
|
#include "vtkMath.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkPlane.h"
|
|
#include "vtkPointData.h"
|
|
#include "vtkPointLocator.h"
|
|
#include "vtkPoints.h"
|
|
#include "vtkPolygon.h"
|
|
#include "vtkQuadric.h"
|
|
|
|
vtkCxxRevisionMacro(vtkTriangle, "$Revision: 1.4.10.1 $");
|
|
vtkStandardNewMacro(vtkTriangle);
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Construct the triangle with three points.
|
|
vtkTriangle::vtkTriangle()
|
|
{
|
|
this->Points->SetNumberOfPoints(3);
|
|
this->PointIds->SetNumberOfIds(3);
|
|
for (int i = 0; i < 3; i++)
|
|
{
|
|
this->Points->SetPoint(i, 0.0, 0.0, 0.0);
|
|
this->PointIds->SetId(i,0);
|
|
}
|
|
this->Line = vtkLine::New();
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
vtkTriangle::~vtkTriangle()
|
|
{
|
|
this->Line->Delete();
|
|
}
|
|
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Create a new cell and copy this triangle's information into the cell.
|
|
// Returns a poiner to the new cell created.
|
|
int vtkTriangle::EvaluatePosition(double x[3], double* closestPoint,
|
|
int& subId, double pcoords[3],
|
|
double& dist2, double *weights)
|
|
{
|
|
int i, j;
|
|
double pt1[3], pt2[3], pt3[3], n[3], fabsn;
|
|
double rhs[2], c1[2], c2[2];
|
|
double det;
|
|
double maxComponent;
|
|
int idx=0, indices[2];
|
|
double dist2Point, dist2Line1, dist2Line2;
|
|
double *closest, closestPoint1[3], closestPoint2[3], cp[3];
|
|
|
|
subId = 0;
|
|
|
|
// Get normal for triangle, only the normal direction is needed, i.e. the
|
|
// normal need not be normalized (unit length)
|
|
//
|
|
this->Points->GetPoint(1, pt1);
|
|
this->Points->GetPoint(2, pt2);
|
|
this->Points->GetPoint(0, pt3);
|
|
|
|
vtkTriangle::ComputeNormalDirection(pt1, pt2, pt3, n);
|
|
|
|
// Project point to plane
|
|
//
|
|
vtkPlane::GeneralizedProjectPoint(x,pt1,n,cp);
|
|
|
|
// Construct matrices. Since we have over determined system, need to find
|
|
// which 2 out of 3 equations to use to develop equations. (Any 2 should
|
|
// work since we've projected point to plane.)
|
|
//
|
|
for (maxComponent=0.0, i=0; i<3; i++)
|
|
{
|
|
// trying to avoid an expensive call to fabs()
|
|
if (n[i] < 0)
|
|
{
|
|
fabsn = -n[i];
|
|
}
|
|
else
|
|
{
|
|
fabsn = n[i];
|
|
}
|
|
if (fabsn > maxComponent)
|
|
{
|
|
maxComponent = fabsn;
|
|
idx = i;
|
|
}
|
|
}
|
|
for (j=0, i=0; i<3; i++)
|
|
{
|
|
if ( i != idx )
|
|
{
|
|
indices[j++] = i;
|
|
}
|
|
}
|
|
|
|
for (i=0; i<2; i++)
|
|
{
|
|
rhs[i] = cp[indices[i]] - pt3[indices[i]];
|
|
c1[i] = pt1[indices[i]] - pt3[indices[i]];
|
|
c2[i] = pt2[indices[i]] - pt3[indices[i]];
|
|
}
|
|
|
|
if ( (det = vtkMath::Determinant2x2(c1,c2)) == 0.0 )
|
|
{
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
|
|
return -1;
|
|
}
|
|
|
|
pcoords[0] = vtkMath::Determinant2x2(rhs,c2) / det;
|
|
pcoords[1] = vtkMath::Determinant2x2(c1,rhs) / det;
|
|
pcoords[2] = 1.0 - (pcoords[0] + pcoords[1]);
|
|
|
|
// Okay, now find closest point to element
|
|
//
|
|
weights[0] = pcoords[2];
|
|
weights[1] = pcoords[0];
|
|
weights[2] = pcoords[1];
|
|
|
|
if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
|
|
pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
|
|
pcoords[2] >= 0.0 && pcoords[2] <= 1.0 )
|
|
{
|
|
//projection distance
|
|
if (closestPoint)
|
|
{
|
|
dist2 = vtkMath::Distance2BetweenPoints(cp,x);
|
|
closestPoint[0] = cp[0];
|
|
closestPoint[1] = cp[1];
|
|
closestPoint[2] = cp[2];
|
|
}
|
|
return 1;
|
|
}
|
|
else
|
|
{
|
|
double t;
|
|
if (closestPoint)
|
|
{
|
|
if ( pcoords[0] < 0.0 && pcoords[1] < 0.0 )
|
|
{
|
|
dist2Point = vtkMath::Distance2BetweenPoints(x,pt3);
|
|
dist2Line1 = vtkLine::DistanceToLine(x,pt1,pt3,t,closestPoint1);
|
|
dist2Line2 = vtkLine::DistanceToLine(x,pt3,pt2,t,closestPoint2);
|
|
if (dist2Point < dist2Line1)
|
|
{
|
|
dist2 = dist2Point;
|
|
closest = pt3;
|
|
}
|
|
else
|
|
{
|
|
dist2 = dist2Line1;
|
|
closest = closestPoint1;
|
|
}
|
|
if (dist2Line2 < dist2)
|
|
{
|
|
dist2 = dist2Line2;
|
|
closest = closestPoint2;
|
|
}
|
|
for (i=0; i<3; i++)
|
|
{
|
|
closestPoint[i] = closest[i];
|
|
}
|
|
}
|
|
else if ( pcoords[1] < 0.0 && pcoords[2] < 0.0 )
|
|
{
|
|
dist2Point = vtkMath::Distance2BetweenPoints(x,pt1);
|
|
dist2Line1 = vtkLine::DistanceToLine(x,pt1,pt3,t,closestPoint1);
|
|
dist2Line2 = vtkLine::DistanceToLine(x,pt1,pt2,t,closestPoint2);
|
|
if (dist2Point < dist2Line1)
|
|
{
|
|
dist2 = dist2Point;
|
|
closest = pt1;
|
|
}
|
|
else
|
|
{
|
|
dist2 = dist2Line1;
|
|
closest = closestPoint1;
|
|
}
|
|
if (dist2Line2 < dist2)
|
|
{
|
|
dist2 = dist2Line2;
|
|
closest = closestPoint2;
|
|
}
|
|
for (i=0; i<3; i++)
|
|
{
|
|
closestPoint[i] = closest[i];
|
|
}
|
|
}
|
|
else if ( pcoords[0] < 0.0 && pcoords[2] < 0.0 )
|
|
{
|
|
dist2Point = vtkMath::Distance2BetweenPoints(x,pt2);
|
|
dist2Line1 = vtkLine::DistanceToLine(x,pt2,pt3,t,closestPoint1);
|
|
dist2Line2 = vtkLine::DistanceToLine(x,pt1,pt2,t,closestPoint2);
|
|
if (dist2Point < dist2Line1)
|
|
{
|
|
dist2 = dist2Point;
|
|
closest = pt2;
|
|
}
|
|
else
|
|
{
|
|
dist2 = dist2Line1;
|
|
closest = closestPoint1;
|
|
}
|
|
if (dist2Line2 < dist2)
|
|
{
|
|
dist2 = dist2Line2;
|
|
closest = closestPoint2;
|
|
}
|
|
for (i=0; i<3; i++)
|
|
{
|
|
closestPoint[i] = closest[i];
|
|
}
|
|
}
|
|
else if ( pcoords[0] < 0.0 )
|
|
{
|
|
dist2 = vtkLine::DistanceToLine(x,pt2,pt3,t,closestPoint);
|
|
}
|
|
else if ( pcoords[1] < 0.0 )
|
|
{
|
|
dist2 = vtkLine::DistanceToLine(x,pt1,pt3,t,closestPoint);
|
|
}
|
|
else if ( pcoords[2] < 0.0 )
|
|
{
|
|
dist2 = vtkLine::DistanceToLine(x,pt1,pt2,t,closestPoint);
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkTriangle::EvaluateLocation(int& vtkNotUsed(subId), double pcoords[3],
|
|
double x[3], double *weights)
|
|
{
|
|
double u3;
|
|
double pt0[3], pt1[3], pt2[3];
|
|
int i;
|
|
|
|
this->Points->GetPoint(0, pt0);
|
|
this->Points->GetPoint(1, pt1);
|
|
this->Points->GetPoint(2, pt2);
|
|
|
|
u3 = 1.0 - pcoords[0] - pcoords[1];
|
|
|
|
for (i=0; i<3; i++)
|
|
{
|
|
x[i] = pt0[i]*u3 + pt1[i]*pcoords[0] + pt2[i]*pcoords[1];
|
|
}
|
|
|
|
weights[0] = u3;
|
|
weights[1] = pcoords[0];
|
|
weights[2] = pcoords[1];
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Compute iso-parametric interpolation functions
|
|
//
|
|
void vtkTriangle::InterpolationFunctions(double pcoords[3], double sf[3])
|
|
{
|
|
sf[0] = 1. - pcoords[0] - pcoords[1];
|
|
sf[1] = pcoords[0];
|
|
sf[2] = pcoords[1];
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkTriangle::InterpolationDerivs(double *, double derivs[6])
|
|
{
|
|
//r-derivatives
|
|
derivs[0] = -1;
|
|
derivs[1] = 1;
|
|
derivs[2] = 0;
|
|
|
|
//s-derivatives
|
|
derivs[3] = -1;
|
|
derivs[4] = 0;
|
|
derivs[5] = 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int vtkTriangle::CellBoundary(int vtkNotUsed(subId), double pcoords[3],
|
|
vtkIdList *pts)
|
|
{
|
|
double t1=pcoords[0]-pcoords[1];
|
|
double t2=0.5*(1.0-pcoords[0])-pcoords[1];
|
|
double t3=2.0*pcoords[0]+pcoords[1]-1.0;
|
|
|
|
pts->SetNumberOfIds(2);
|
|
|
|
// compare against three lines in parametric space that divide element
|
|
// into three pieces
|
|
if ( t1 >= 0.0 && t2 >= 0.0 )
|
|
{
|
|
pts->SetId(0,this->PointIds->GetId(0));
|
|
pts->SetId(1,this->PointIds->GetId(1));
|
|
}
|
|
|
|
else if ( t2 < 0.0 && t3 >= 0.0 )
|
|
{
|
|
pts->SetId(0,this->PointIds->GetId(1));
|
|
pts->SetId(1,this->PointIds->GetId(2));
|
|
}
|
|
|
|
else //( t1 < 0.0 && t3 < 0.0 )
|
|
{
|
|
pts->SetId(0,this->PointIds->GetId(2));
|
|
pts->SetId(1,this->PointIds->GetId(0));
|
|
}
|
|
|
|
if ( pcoords[0] < 0.0 || pcoords[1] < 0.0 ||
|
|
pcoords[0] > 1.0 || pcoords[1] > 1.0 ||
|
|
(1.0 - pcoords[0] - pcoords[1]) < 0.0 )
|
|
{
|
|
return 0;
|
|
}
|
|
else
|
|
{
|
|
return 1;
|
|
}
|
|
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
//
|
|
// Marching triangles
|
|
//
|
|
typedef int EDGE_LIST;
|
|
typedef struct {
|
|
EDGE_LIST edges[3];
|
|
} LINE_CASES;
|
|
|
|
static LINE_CASES lineCases[] = {
|
|
{{-1, -1, -1}},
|
|
{{0, 2, -1}},
|
|
{{1, 0, -1}},
|
|
{{1, 2, -1}},
|
|
{{2, 1, -1}},
|
|
{{0, 1, -1}},
|
|
{{2, 0, -1}},
|
|
{{-1, -1, -1}}
|
|
};
|
|
|
|
static int edges[3][2] = { {0,1}, {1,2}, {2,0} };
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkTriangle::Contour(double value, vtkDataArray *cellScalars,
|
|
vtkPointLocator *locator,
|
|
vtkCellArray *verts,
|
|
vtkCellArray *lines,
|
|
vtkCellArray *vtkNotUsed(polys),
|
|
vtkPointData *inPd, vtkPointData *outPd,
|
|
vtkCellData *inCd, vtkIdType cellId,
|
|
vtkCellData *outCd)
|
|
{
|
|
static int CASE_MASK[3] = {1,2,4};
|
|
LINE_CASES *lineCase;
|
|
EDGE_LIST *edge;
|
|
int i, j, index, *vert;
|
|
vtkIdType pts[2];
|
|
int e1, e2, newCellId;
|
|
double t, x1[3], x2[3], x[3], deltaScalar;
|
|
vtkIdType offset = verts->GetNumberOfCells();
|
|
|
|
// Build the case table
|
|
for ( i=0, index = 0; i < 3; i++)
|
|
{
|
|
if (cellScalars->GetComponent(i,0) >= value)
|
|
{
|
|
index |= CASE_MASK[i];
|
|
}
|
|
}
|
|
|
|
lineCase = lineCases + index;
|
|
edge = lineCase->edges;
|
|
|
|
for ( ; edge[0] > -1; edge += 2 )
|
|
{
|
|
for (i=0; i<2; i++) // insert line
|
|
{
|
|
vert = edges[edge[i]];
|
|
// calculate a preferred interpolation direction
|
|
deltaScalar = (cellScalars->GetComponent(vert[1],0)
|
|
- cellScalars->GetComponent(vert[0],0));
|
|
if (deltaScalar > 0)
|
|
{
|
|
e1 = vert[0]; e2 = vert[1];
|
|
}
|
|
else
|
|
{
|
|
e1 = vert[1]; e2 = vert[0];
|
|
deltaScalar = -deltaScalar;
|
|
}
|
|
|
|
// linear interpolation
|
|
if (deltaScalar == 0.0)
|
|
{
|
|
t = 0.0;
|
|
}
|
|
else
|
|
{
|
|
t = (value - cellScalars->GetComponent(e1,0)) / deltaScalar;
|
|
}
|
|
|
|
this->Points->GetPoint(e1, x1);
|
|
this->Points->GetPoint(e2, x2);
|
|
|
|
for (j=0; j<3; j++)
|
|
{
|
|
x[j] = x1[j] + t * (x2[j] - x1[j]);
|
|
}
|
|
if ( locator->InsertUniquePoint(x, pts[i]) )
|
|
{
|
|
if ( outPd )
|
|
{
|
|
vtkIdType p1 = this->PointIds->GetId(e1);
|
|
vtkIdType p2 = this->PointIds->GetId(e2);
|
|
outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
|
|
}
|
|
}
|
|
}
|
|
// check for degenerate line
|
|
if ( pts[0] != pts[1] )
|
|
{
|
|
newCellId = offset + lines->InsertNextCell(2,pts);
|
|
outCd->CopyData(inCd,cellId,newCellId);
|
|
}
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Get the edge specified by edgeId (range 0 to 2) and return that edge's
|
|
// coordinates.
|
|
vtkCell *vtkTriangle::GetEdge(int edgeId)
|
|
{
|
|
int edgeIdPlus1 = (edgeId > 1 ? 0 : (edgeId+1) );
|
|
|
|
// load point id's
|
|
this->Line->PointIds->SetId(0,this->PointIds->GetId(edgeId));
|
|
this->Line->PointIds->SetId(1,this->PointIds->GetId(edgeIdPlus1));
|
|
|
|
// load coordinates
|
|
this->Line->Points->SetPoint(0,this->Points->GetPoint(edgeId));
|
|
this->Line->Points->SetPoint(1,this->Points->GetPoint(edgeIdPlus1));
|
|
|
|
return this->Line;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Plane intersection plus in/out test on triangle. The in/out test is
|
|
// performed using tol as the tolerance.
|
|
int vtkTriangle::IntersectWithLine(double p1[3], double p2[3], double tol,
|
|
double& t, double x[3], double pcoords[3],
|
|
int& subId)
|
|
{
|
|
double pt1[3], pt2[3], pt3[3], n[3];
|
|
double tol2 = tol*tol;
|
|
double closestPoint[3];
|
|
double dist2, weights[3];
|
|
|
|
subId = 0;
|
|
|
|
// Get normal for triangle
|
|
//
|
|
this->Points->GetPoint(1, pt1);
|
|
this->Points->GetPoint(2, pt2);
|
|
this->Points->GetPoint(0, pt3);
|
|
|
|
vtkTriangle::ComputeNormal (pt1, pt2, pt3, n);
|
|
|
|
// Intersect plane of triangle with line
|
|
//
|
|
if ( ! vtkPlane::IntersectWithLine(p1,p2,n,pt1,t,x) )
|
|
{
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
|
|
return 0;
|
|
}
|
|
|
|
// Evaluate position
|
|
//
|
|
int inside;
|
|
if ( (inside = this->EvaluatePosition(x, closestPoint, subId, pcoords,
|
|
dist2, weights)) >= 0)
|
|
{
|
|
if ( dist2 <= tol2 )
|
|
{
|
|
pcoords[2] = 0.0;
|
|
return 1;
|
|
}
|
|
return inside;
|
|
}
|
|
|
|
// so the easy test failed. The line is not intersecting the triangle.
|
|
// Let's now do the 3d case check to see how close the line comes.
|
|
// basically we just need to test against the three lines of the triangle
|
|
this->Line->PointIds->InsertId(0,0);
|
|
this->Line->PointIds->InsertId(1,1);
|
|
|
|
if (pcoords[2] < 0.0)
|
|
{
|
|
this->Line->Points->InsertPoint(0,pt1);
|
|
this->Line->Points->InsertPoint(1,pt2);
|
|
if (this->Line->IntersectWithLine(p1,p2,tol,t,x,pcoords,subId))
|
|
{
|
|
pcoords[2] = 0.0;
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
if (pcoords[0] < 0.0)
|
|
{
|
|
this->Line->Points->InsertPoint(0,pt2);
|
|
this->Line->Points->InsertPoint(1,pt3);
|
|
if (this->Line->IntersectWithLine(p1,p2,tol,t,x,pcoords,subId))
|
|
{
|
|
pcoords[2] = 0.0;
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
if (pcoords[1] < 0.0)
|
|
{
|
|
this->Line->Points->InsertPoint(0,pt3);
|
|
this->Line->Points->InsertPoint(1,pt1);
|
|
if (this->Line->IntersectWithLine(p1,p2,tol,t,x,pcoords,subId))
|
|
{
|
|
pcoords[2] = 0.0;
|
|
return 1;
|
|
}
|
|
}
|
|
|
|
pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
|
|
return 0;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int vtkTriangle::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
|
|
vtkPoints *pts)
|
|
{
|
|
pts->Reset();
|
|
ptIds->Reset();
|
|
|
|
for ( int i=0; i < 3; i++ )
|
|
{
|
|
ptIds->InsertId(i,this->PointIds->GetId(i));
|
|
pts->InsertPoint(i,this->Points->GetPoint(i));
|
|
}
|
|
|
|
return 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Used a staged computation: first compute derivatives in local x'-y'
|
|
// coordinate system; then convert into x-y-z modelling system.
|
|
void vtkTriangle::Derivatives(int vtkNotUsed(subId), double vtkNotUsed(pcoords)[3],
|
|
double *values, int dim, double *derivs)
|
|
{
|
|
double v0[2], v1[2], v2[2], v[3], v10[3], v20[3], lenX;
|
|
double x0[3], x1[3], x2[3], n[3];
|
|
double *J[2], J0[2], J1[2];
|
|
double *JI[2], JI0[2], JI1[2];
|
|
double functionDerivs[6], sum[2], dBydx, dBydy;
|
|
int i, j;
|
|
|
|
// Project points of triangle into 2D system
|
|
this->Points->GetPoint(0, x0);
|
|
this->Points->GetPoint(1, x1);
|
|
this->Points->GetPoint(2, x2);
|
|
vtkTriangle::ComputeNormal (x0, x1, x2, n);
|
|
|
|
for (i=0; i < 3; i++)
|
|
{
|
|
v10[i] = x1[i] - x0[i];
|
|
v[i] = x2[i] - x0[i];
|
|
}
|
|
|
|
vtkMath::Cross(n,v10,v20); //creates local y' axis
|
|
|
|
if ( (lenX=vtkMath::Normalize(v10)) <= 0.0
|
|
|| vtkMath::Normalize(v20) <= 0.0 ) //degenerate
|
|
{
|
|
for ( j=0; j < dim; j++ )
|
|
{
|
|
for ( i=0; i < 3; i++ )
|
|
{
|
|
derivs[j*dim + i] = 0.0;
|
|
}
|
|
}
|
|
return;
|
|
}
|
|
|
|
v0[0] = v0[1] = 0.0; //convert points to 2D (i.e., local system)
|
|
v1[0] = lenX; v1[1] = 0.0;
|
|
v2[0] = vtkMath::Dot(v,v10);
|
|
v2[1] = vtkMath::Dot(v,v20);
|
|
|
|
// Compute interpolation function derivatives
|
|
vtkTriangle::InterpolationDerivs(NULL,functionDerivs);
|
|
|
|
// Compute Jacobian: Jacobian is constant for a triangle.
|
|
J[0] = J0; J[1] = J1;
|
|
JI[0] = JI0; JI[1] = JI1;
|
|
|
|
J[0][0] = v1[0] - v0[0];
|
|
J[1][0] = v2[0] - v0[0];
|
|
J[0][1] = v1[1] - v0[1];
|
|
J[1][1] = v2[1] - v0[1];
|
|
|
|
// Compute inverse Jacobian
|
|
vtkMath::InvertMatrix(J,JI,2);
|
|
|
|
// Loop over "dim" derivative values. For each set of values, compute
|
|
// derivatives in local system and then transform into modelling system.
|
|
// First compute derivatives in local x'-y' coordinate system
|
|
for ( j=0; j < dim; j++ )
|
|
{
|
|
sum[0] = sum[1] = 0.0;
|
|
for ( i=0; i < 3; i++) //loop over interp. function derivatives
|
|
{
|
|
sum[0] += functionDerivs[i] * values[dim*i + j];
|
|
sum[1] += functionDerivs[3 + i] * values[dim*i + j];
|
|
}
|
|
dBydx = sum[0]*JI[0][0] + sum[1]*JI[0][1];
|
|
dBydy = sum[0]*JI[1][0] + sum[1]*JI[1][1];
|
|
|
|
// Transform into global system (dot product with global axes)
|
|
derivs[3*j] = dBydx * v10[0] + dBydy * v20[0];
|
|
derivs[3*j + 1] = dBydx * v10[1] + dBydy * v20[1];
|
|
derivs[3*j + 2] = dBydx * v10[2] + dBydy * v20[2];
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Compute the triangle normal from a points list, and a list of point ids
|
|
// that index into the points list.
|
|
void vtkTriangle::ComputeNormal(vtkPoints *p, int vtkNotUsed(numPts),
|
|
vtkIdType *pts, double n[3])
|
|
{
|
|
double v1[3], v2[3], v3[3];
|
|
|
|
p->GetPoint(pts[0],v1);
|
|
p->GetPoint(pts[1],v2);
|
|
p->GetPoint(pts[2],v3);
|
|
|
|
vtkTriangle::ComputeNormal(v1,v2,v3,n);
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Compute the circumcenter (center[3]) and radius squared (method
|
|
// return value) of a triangle defined by the three points x1, x2, and
|
|
// x3. (Note that the coordinates are 2D. 3D points can be used but
|
|
// the z-component will be ignored.)
|
|
double vtkTriangle::Circumcircle(double x1[2], double x2[2], double x3[2],
|
|
double center[2])
|
|
{
|
|
double n12[2], n13[2], x12[2], x13[2];
|
|
double *A[2], rhs[2], sum, diff;
|
|
int i;
|
|
|
|
// calculate normals and intersection points of bisecting planes.
|
|
//
|
|
for (i=0; i<2; i++)
|
|
{
|
|
n12[i] = x2[i] - x1[i];
|
|
n13[i] = x3[i] - x1[i];
|
|
x12[i] = (x2[i] + x1[i])/2.0;
|
|
x13[i] = (x3[i] + x1[i])/2.0;
|
|
}
|
|
|
|
// Compute solutions to the intersection of two bisecting lines
|
|
// (2-eqns. in 2-unknowns).
|
|
//
|
|
// form system matrices
|
|
//
|
|
A[0] = n12;
|
|
A[1] = n13;
|
|
|
|
rhs[0] = vtkMath::Dot2D(n12,x12);
|
|
rhs[1] = vtkMath::Dot2D(n13,x13);
|
|
|
|
// Solve system of equations
|
|
//
|
|
if ( vtkMath::SolveLinearSystem(A,rhs,2) == 0 )
|
|
{
|
|
center[0] = center[1] = 0.0;
|
|
return VTK_DOUBLE_MAX;
|
|
}
|
|
else
|
|
{
|
|
center[0] = rhs[0]; center[1] = rhs[1];
|
|
}
|
|
|
|
//determine average value of radius squared
|
|
for (sum=0, i=0; i<2; i++)
|
|
{
|
|
diff = x1[i] - center[i];
|
|
sum += diff*diff;
|
|
diff = x2[i] - center[i];
|
|
sum += diff*diff;
|
|
diff = x3[i] - center[i];
|
|
sum += diff*diff;
|
|
}
|
|
|
|
if ( (sum /= 3.0) > VTK_DOUBLE_MAX )
|
|
{
|
|
return VTK_DOUBLE_MAX;
|
|
}
|
|
else
|
|
{
|
|
return sum;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Given a 2D point x[2], determine the barycentric coordinates of the point.
|
|
// Barycentric coordinates are a natural coordinate system for simplices that
|
|
// express a position as a linear combination of the vertices. For a
|
|
// triangle, there are three barycentric coordinates (because there are
|
|
// fourthree vertices), and the sum of the coordinates must equal 1. If a
|
|
// point x is inside a simplex, then all three coordinates will be strictly
|
|
// positive. If two coordinates are zero (so the third =1), then the
|
|
// point x is on a vertex. If one coordinates are zero, the point x is on an
|
|
// edge. In this method, you must specify the vertex coordinates x1->x3.
|
|
// Returns 0 if triangle is degenerate.
|
|
int vtkTriangle::BarycentricCoords(double x[2], double x1[2], double x2[2],
|
|
double x3[2], double bcoords[3])
|
|
{
|
|
double *A[3], p[3], a1[3], a2[3], a3[3];
|
|
int i;
|
|
|
|
// Homogenize the variables; load into arrays.
|
|
//
|
|
a1[0] = x1[0]; a1[1] = x2[0]; a1[2] = x3[0];
|
|
a2[0] = x1[1]; a2[1] = x2[1]; a2[2] = x3[1];
|
|
a3[0] = 1.0; a3[1] = 1.0; a3[2] = 1.0;
|
|
p[0] = x[0]; p[1] = x[1]; p[2] = 1.0;
|
|
|
|
// Now solve system of equations for barycentric coordinates
|
|
//
|
|
A[0] = a1;
|
|
A[1] = a2;
|
|
A[2] = a3;
|
|
|
|
if ( vtkMath::SolveLinearSystem(A,p,3) )
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
bcoords[i] = p[i];
|
|
}
|
|
return 1;
|
|
}
|
|
else
|
|
{
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Project triangle defined in 3D to 2D coordinates. Returns 0 if degenerate
|
|
// triangle; non-zero value otherwise. Input points are x1->x3; output 2D
|
|
// points are v1->v3.
|
|
int vtkTriangle::ProjectTo2D(double x1[3], double x2[3], double x3[3],
|
|
double v1[2], double v2[2], double v3[2])
|
|
{
|
|
double n[3], v21[3], v31[3], v[3], xLen;
|
|
|
|
// Get normal for triangle
|
|
vtkTriangle::ComputeNormal (x1, x2, x3, n);
|
|
|
|
for (int i=0; i < 3; i++)
|
|
{
|
|
v21[i] = x2[i] - x1[i];
|
|
v31[i] = x3[i] - x1[i];
|
|
}
|
|
|
|
if ( (xLen=vtkMath::Normalize(v21)) <= 0.0 )
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
// The first point is at (0,0); the next at (xLen,0); compute the other
|
|
// point relative to the first two.
|
|
v1[0] = v1[1] = 0.0;
|
|
v2[0] = xLen; v2[1] = 0.0;
|
|
|
|
vtkMath::Cross(n,v21,v);
|
|
|
|
v3[0] = vtkMath::Dot(v31,v21);
|
|
v3[1] = vtkMath::Dot(v31,v);
|
|
|
|
return 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Support triangle clipping. Note that the table defines triangles (three ids
|
|
// at a time define a triangle, -1 ends the list). Numbers in the list >= 100
|
|
// correspond to already existing vertices; otherwise the numbers refer to edge
|
|
// ids.
|
|
typedef int TRIANGLE_EDGE_LIST;
|
|
typedef struct {
|
|
TRIANGLE_EDGE_LIST edges[7];
|
|
} TRIANGLE_CASES;
|
|
|
|
static TRIANGLE_CASES triangleCases[] = {
|
|
{{-1, -1, -1, -1, -1, -1, -1}}, // 0
|
|
{{0, 2, 100, -1, -1, -1, -1}}, // 1
|
|
{{1, 0, 101, -1, -1, -1, -1}}, // 2
|
|
{{1, 2, 100, 1, 100, 101, -1}}, // 3
|
|
{{2, 1, 102, -1, -1, -1, -1}}, // 4
|
|
{{0, 1, 102, 102, 100, 0, -1}}, // 5
|
|
{{0, 101, 2, 2, 101, 102, -1}}, // 6
|
|
{{100, 101, 102, -1, -1, -1, -1}} // 7
|
|
};
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Clip this triangle using scalar value provided. Like contouring, except
|
|
// that it cuts the triangle to produce other triangles.
|
|
void vtkTriangle::Clip(double value, vtkDataArray *cellScalars,
|
|
vtkPointLocator *locator, vtkCellArray *tris,
|
|
vtkPointData *inPd, vtkPointData *outPd,
|
|
vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd,
|
|
int insideOut)
|
|
{
|
|
static int CASE_MASK[3] = {1,2,4};
|
|
TRIANGLE_CASES *triangleCase;
|
|
TRIANGLE_EDGE_LIST *edge;
|
|
int i, j, index, *vert;
|
|
int e1, e2, newCellId;
|
|
vtkIdType pts[3];
|
|
int vertexId;
|
|
double t, x1[3], x2[3], x[3], deltaScalar;
|
|
|
|
// Build the case table
|
|
if ( insideOut )
|
|
{
|
|
for ( i=0, index = 0; i < 3; i++)
|
|
{
|
|
if (cellScalars->GetComponent(i,0) <= value)
|
|
{
|
|
index |= CASE_MASK[i];
|
|
}
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for ( i=0, index = 0; i < 3; i++)
|
|
{
|
|
if (cellScalars->GetComponent(i,0) > value)
|
|
{
|
|
index |= CASE_MASK[i];
|
|
}
|
|
}
|
|
}
|
|
|
|
// Select the case based on the index and get the list of edges for this case
|
|
triangleCase = triangleCases + index;
|
|
edge = triangleCase->edges;
|
|
|
|
// generate each triangle
|
|
for ( ; edge[0] > -1; edge += 3 )
|
|
{
|
|
for (i=0; i<3; i++) // insert triangle
|
|
{
|
|
// vertex exists, and need not be interpolated
|
|
if (edge[i] >= 100)
|
|
{
|
|
vertexId = edge[i] - 100;
|
|
this->Points->GetPoint(vertexId, x);
|
|
if ( locator->InsertUniquePoint(x, pts[i]) )
|
|
{
|
|
outPd->CopyData(inPd,this->PointIds->GetId(vertexId),pts[i]);
|
|
}
|
|
}
|
|
|
|
else //new vertex, interpolate
|
|
{
|
|
vert = edges[edge[i]];
|
|
|
|
// calculate a preferred interpolation direction
|
|
deltaScalar = (cellScalars->GetComponent(vert[1],0) -
|
|
cellScalars->GetComponent(vert[0],0));
|
|
if (deltaScalar > 0)
|
|
{
|
|
e1 = vert[0]; e2 = vert[1];
|
|
}
|
|
else
|
|
{
|
|
e1 = vert[1]; e2 = vert[0];
|
|
deltaScalar = -deltaScalar;
|
|
}
|
|
|
|
// linear interpolation
|
|
if (deltaScalar == 0.0)
|
|
{
|
|
t = 0.0;
|
|
}
|
|
else
|
|
{
|
|
t = (value - cellScalars->GetComponent(e1,0)) / deltaScalar;
|
|
}
|
|
|
|
this->Points->GetPoint(e1, x1);
|
|
this->Points->GetPoint(e2, x2);
|
|
|
|
for (j=0; j<3; j++)
|
|
{
|
|
x[j] = x1[j] + t * (x2[j] - x1[j]);
|
|
}
|
|
if ( locator->InsertUniquePoint(x, pts[i]) )
|
|
{
|
|
vtkIdType p1 = this->PointIds->GetId(e1);
|
|
vtkIdType p2 = this->PointIds->GetId(e2);
|
|
outPd->InterpolateEdge(inPd,pts[i],p1,p2,t);
|
|
}
|
|
}
|
|
}
|
|
// check for degenerate tri's
|
|
if (pts[0] == pts[1] || pts[0] == pts[2] || pts[1] == pts[2])
|
|
{
|
|
continue;
|
|
}
|
|
|
|
newCellId = tris->InsertNextCell(3,pts);
|
|
outCd->CopyData(inCd,cellId,newCellId);
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// Given a point x, determine whether it is inside (within the
|
|
// tolerance squared, tol2) the triangle defined by the three
|
|
// coordinate values p1, p2, p3. Method is via comparing dot products.
|
|
// (Note: in current implementation the tolerance only works in the
|
|
// neighborhood of the three vertices of the triangle.
|
|
int vtkTriangle::PointInTriangle(double x[3], double p1[3], double p2[3],
|
|
double p3[3], double tol2)
|
|
{
|
|
double x1[3], x2[3], x3[3], v13[3], v21[3], v32[3];
|
|
double n1[3], n2[3], n3[3];
|
|
int i;
|
|
|
|
// Compute appropriate vectors
|
|
//
|
|
for (i=0; i<3; i++)
|
|
{
|
|
x1[i] = x[i] - p1[i];
|
|
x2[i] = x[i] - p2[i];
|
|
x3[i] = x[i] - p3[i];
|
|
v13[i] = p1[i] - p3[i];
|
|
v21[i] = p2[i] - p1[i];
|
|
v32[i] = p3[i] - p2[i];
|
|
}
|
|
|
|
// See whether intersection point is within tolerance of a vertex.
|
|
//
|
|
if ( (x1[0]*x1[0] + x1[1]*x1[1] + x1[2]*x1[2]) <= tol2 ||
|
|
(x2[0]*x2[0] + x2[1]*x2[1] + x2[2]*x2[2]) <= tol2 ||
|
|
(x3[0]*x3[0] + x3[1]*x3[1] + x3[2]*x3[2]) <= tol2 )
|
|
{
|
|
return 1;
|
|
}
|
|
|
|
// If not near a vertex, check whether point is inside of triangular face.
|
|
//
|
|
// Obtain normal off of triangular face
|
|
//
|
|
vtkMath::Cross (x1, v13, n1);
|
|
vtkMath::Cross (x2, v21, n2);
|
|
vtkMath::Cross (x3, v32, n3);
|
|
|
|
// Check whether normals go in same direction
|
|
//
|
|
if ( (vtkMath::Dot(n1,n2) >= 0.0) && (vtkMath::Dot(n2,n3) >= 0.0) )
|
|
{
|
|
return 1;
|
|
}
|
|
else
|
|
{
|
|
return 0;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
double vtkTriangle::GetParametricDistance(double pcoords[3])
|
|
{
|
|
int i;
|
|
double pDist, pDistMax=0.0;
|
|
double pc[3];
|
|
|
|
pc[0] = pcoords[0];
|
|
pc[1] = pcoords[1];
|
|
pc[2] = 1.0 - pcoords[0] - pcoords[1];
|
|
|
|
for (i=0; i<3; i++)
|
|
{
|
|
if ( pc[i] < 0.0 )
|
|
{
|
|
pDist = -pc[i];
|
|
}
|
|
else if ( pc[i] > 1.0 )
|
|
{
|
|
pDist = pc[i] - 1.0;
|
|
}
|
|
else //inside the cell in the parametric direction
|
|
{
|
|
pDist = 0.0;
|
|
}
|
|
if ( pDist > pDistMax )
|
|
{
|
|
pDistMax = pDist;
|
|
}
|
|
}
|
|
|
|
return pDistMax;
|
|
}
|
|
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkTriangle::ComputeQuadric(double x1[3], double x2[3], double x3[3],
|
|
double quadric[4][4])
|
|
{
|
|
double crossX1X2[3], crossX2X3[3], crossX3X1[3];
|
|
double determinantABC;
|
|
double ABCx[3][3];
|
|
double n[4];
|
|
int i, j;
|
|
|
|
for (i = 0; i < 3; i++)
|
|
{
|
|
ABCx[0][i] = x1[i];
|
|
ABCx[1][i] = x2[i];
|
|
ABCx[2][i] = x3[i];
|
|
}
|
|
|
|
vtkMath::Cross(x1, x2, crossX1X2);
|
|
vtkMath::Cross(x2, x3, crossX2X3);
|
|
vtkMath::Cross(x3, x1, crossX3X1);
|
|
determinantABC = vtkMath::Determinant3x3(ABCx);
|
|
|
|
n[0] = crossX1X2[0] + crossX2X3[0] + crossX3X1[0];
|
|
n[1] = crossX1X2[1] + crossX2X3[1] + crossX3X1[1];
|
|
n[2] = crossX1X2[2] + crossX2X3[2] + crossX3X1[2];
|
|
n[3] = -determinantABC;
|
|
|
|
for (i = 0; i < 4; i++)
|
|
{
|
|
for (j = 0; j < 4; j++)
|
|
{
|
|
quadric[i][j] = n[i] * n[j];
|
|
}
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkTriangle::ComputeQuadric(double x1[3], double x2[3], double x3[3],
|
|
vtkQuadric *quadric)
|
|
{
|
|
double quadricMatrix[4][4];
|
|
|
|
ComputeQuadric(x1, x2, x3, quadricMatrix);
|
|
quadric->SetCoefficients(quadricMatrix[0][0], quadricMatrix[1][1],
|
|
quadricMatrix[2][2], 2*quadricMatrix[0][1],
|
|
2*quadricMatrix[1][2], 2*quadricMatrix[0][2],
|
|
2*quadricMatrix[0][3], 2*quadricMatrix[1][3],
|
|
2*quadricMatrix[2][3], quadricMatrix[3][3]);
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
static double vtkTriangleCellPCoords[9] =
|
|
{0.0,0.0,0.0, 1.0,0.0,0.0, 0.0,1.0,0.0};
|
|
double *vtkTriangle::GetParametricCoords()
|
|
{
|
|
return vtkTriangleCellPCoords;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkTriangle::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
os << indent << "Line:\n";
|
|
this->Line->PrintSelf(os,indent.GetNextIndent());
|
|
}
|
|
|