Permalink
Switch branches/tags
Nothing to show
Find file
Fetching contributors…
Cannot retrieve contributors at this time
890 lines (694 sloc) 21.3 KB
/*=========================================================================
Program: Visualization Toolkit
Module: vtkPlanesIntersection.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.
=========================================================================*/
/*----------------------------------------------------------------------------
Copyright (c) Sandia Corporation
See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/
#include "vtkMath.h"
#include "vtkPlanesIntersection.h"
#include "vtkPointsProjectedHull.h"
#include "vtkFloatArray.h"
#include "vtkCell.h"
#include "vtkObjectFactory.h"
vtkStandardNewMacro(vtkPlanesIntersection);
// Experiment shows that we get plane equation values on the
// order of 10e-6 when the point is actually on the plane
#define VTK_SMALL_DOUBLE (10e-5)
const int Inside = 0;
const int Outside = 1;
const int Straddle = 2;
const int Xdim=0; // don't change these three values
const int Ydim=1;
const int Zdim=2;
vtkPlanesIntersection::vtkPlanesIntersection()
{
this->Plane = NULL;
this->regionPts = NULL;
}
vtkPlanesIntersection::~vtkPlanesIntersection()
{
if (this->regionPts)
{
this->regionPts->Delete();
this->regionPts = NULL;
}
if (this->Plane)
{
delete [] this->Plane;
this->Plane = NULL;
}
}
void vtkPlanesIntersection::SetRegionVertices(vtkPoints *v)
{
int i;
if (this->regionPts) this->regionPts->Delete();
this->regionPts = vtkPointsProjectedHull::New();
if (v->GetDataType() == VTK_DOUBLE)
{
this->regionPts->DeepCopy(v);
}
else
{
this->regionPts->SetDataTypeToDouble();
int npts = v->GetNumberOfPoints();
this->regionPts->SetNumberOfPoints(npts);
double *pt;
for (i=0; i<npts; i++)
{
pt = v->GetPoint(i);
regionPts->SetPoint(i, pt[0], pt[1], pt[2]);
}
}
}
void vtkPlanesIntersection::SetRegionVertices(double *v, int nvertices)
{
int i;
if (this->regionPts) this->regionPts->Delete();
this->regionPts = vtkPointsProjectedHull::New();
this->regionPts->SetDataTypeToDouble();
this->regionPts->SetNumberOfPoints(nvertices);
for (i=0; i<nvertices; i++)
{
this->regionPts->SetPoint(i, v + (i*3));
}
}
int vtkPlanesIntersection::GetRegionVertices(double *v, int nvertices)
{
int i;
if (this->regionPts == NULL) this->ComputeRegionVertices();
int npts = this->regionPts->GetNumberOfPoints();
if (npts > nvertices) npts = nvertices;
for (i=0; i<npts; i++)
{
this->regionPts->GetPoint(i, v + i*3);
}
return npts;
}
int vtkPlanesIntersection::GetNumRegionVertices()
{
if (this->regionPts == NULL) this->ComputeRegionVertices();
return this->regionPts->GetNumberOfPoints();
}
//---------------------------------------------------------------------
// Determine whether the axis aligned box provided intersects
// the convex region bounded by the planes.
//---------------------------------------------------------------------
int vtkPlanesIntersection::IntersectsRegion(vtkPoints *R)
{
int plane;
int allInside=0;
int nplanes = this->GetNumberOfPlanes();
if (nplanes < 4)
{
vtkErrorMacro("invalid region - less than 4 planes");
return 0;
}
if (this->regionPts == NULL)
{
this->ComputeRegionVertices();
if (this->regionPts->GetNumberOfPoints() < 4)
{
vtkErrorMacro("Invalid region: zero-volume intersection");
return 0;
}
}
if (R->GetNumberOfPoints() < 8)
{
vtkErrorMacro("invalid box");
return 0;
}
int *where = new int[nplanes];
int intersects = -1;
// Here's the algorithm from Graphics Gems IV, page 81,
//
// R is an axis aligned box (could represent a region in a spatial
// spatial partitioning of a volume of data).
//
// P is a set of planes defining a convex region in space (could be
// a view frustum).
//
// The question is does P intersect R. We expect to be doing the
// calculation for one P and many Rs.
// You may wonder why we don't do what vtkClipPolyData does, which
// computes the following on every point of it's PolyData input:
//
// for each point in the input
// for each plane defining the convex region
// evaluate plane eq to determine if outside, inside or on plane
//
// For each cell, if some points are inside and some outside, then
// vtkClipPolyData decides it straddles the region and clips it. If
// every point is inside, it tosses it.
//
// The reason is that the Graphics Gems algorithm is faster in some
// cases (we may only need to evaluate one vertex of the box). And
// also because if the convex region passes through the box without
// including any vertices of the box, all box vertices will be
// "outside" and the algorithm will fail. vtkClipPolyData assumes
// cells are very small relative to the clip region. In general
// the axis-aligned box may be a large portion of world coordinate
// space, and the convex region a view frustum representing a
// small portion of the screen.
// 1. If R does not intersect P's bounding box, return 0.
if (this->IntersectsBoundingBox(R) == 0)
{
intersects = 0;
}
// 2. If P's bounding box is entirely inside R, return 1.
else if (this->EnclosesBoundingBox(R) == 1)
{
intersects = 1;
}
// 3. For each face plane F of P
//
// Suppose the plane equation is negative inside P and
// positive outside P. Choose the vertex (n) of R which is
// most in the direction of the negative pointing normal of
// the plane. The opposite vertex (p) is most in the
// direction of the positive pointing normal. (This is
// a very quick calculation.)
//
// If n is on the positive side of the plane, R is
// completely outside of P, so return 0.
//
// If n and p are both on the negative side, then R is on
// the "inside" of F. Keep track to see if all R is inside
// all planes defining the region.
else
{
if (this->Plane == NULL) this->SetPlaneEquations();
allInside = 1;
for (plane=0; plane < nplanes; plane++)
{
where[plane] = this->EvaluateFacePlane(plane, R);
if (allInside && (where[plane] != Inside))
{
allInside = 0;
}
if (where[plane] == Outside)
{
intersects = 0;
break;
}
}
}
if (intersects == -1)
{
// 4. If n and p were "inside" all faces, R is inside P
// so return 1.
if ( allInside)
{
intersects = 1;
}
// 5. For each of three orthographic projections (X, Y and Z)
//
// Compute the equations of the edge lines of P in those views.
//
// If R's projection lies outside any of these lines (using 2D
// version of n & p tests), return 0.
else if ((this->IntersectsProjection(R, Xdim) == 0) ||
(this->IntersectsProjection(R, Ydim) == 0) ||
(this->IntersectsProjection(R, Zdim) == 0) )
{
}
else
{
// 6. Return 1.
intersects = 1;
}
}
delete [] where;
return (intersects==1);
}
// a static convenience function - since we have all the machinery
// in this class, we can compute whether an arbitrary polygon intersects
// an axis aligned box
//
// it is assumed "pts" represents a planar polygon
//
int vtkPlanesIntersection::PolygonIntersectsBBox(double bounds[6], vtkPoints *pts)
{
// a bogus vtkPlanesIntersection object containing only one plane
vtkPlanesIntersection *pi = vtkPlanesIntersection::New();
pi->SetRegionVertices(pts);
vtkPoints *Box = vtkPoints::New();
Box->SetNumberOfPoints(8);
Box->SetPoint(0, bounds[0], bounds[2], bounds[4]);
Box->SetPoint(1, bounds[1], bounds[2], bounds[4]);
Box->SetPoint(2, bounds[1], bounds[3], bounds[4]);
Box->SetPoint(3, bounds[0], bounds[3], bounds[4]);
Box->SetPoint(4, bounds[0], bounds[2], bounds[5]);
Box->SetPoint(5, bounds[1], bounds[2], bounds[5]);
Box->SetPoint(6, bounds[1], bounds[3], bounds[5]);
Box->SetPoint(7, bounds[0], bounds[3], bounds[5]);
int intersects = -1;
// 1. Does Box intersect the polygon's bounding box?
if (pi->IntersectsBoundingBox(Box) == 0)
{
intersects = 0;
}
// 2. If so, does Box entirely contain the polygon's bounding box?
else if (pi->EnclosesBoundingBox(Box) == 1)
{
intersects = 1;
}
if (intersects == -1)
{
// 3. If not, determine whether the Box intersects the plane of the polygon
vtkPoints *origin = vtkPoints::New();
origin->SetNumberOfPoints(1);
origin->SetPoint(0, pts->GetPoint(0));
vtkFloatArray *normal = vtkFloatArray::New();
normal->SetNumberOfComponents(3);
normal->SetNumberOfTuples(1);
// find 3 points that are not co-linear and compute a normal
double nvec[3], p0[3], p1[3], pp[3];
int npts = pts->GetNumberOfPoints();
pts->GetPoint(0, p0);
pts->GetPoint(1, p1);
for (int p = 2; p < npts; p++)
{
pts->GetPoint(p, pp);
vtkPlanesIntersection::ComputeNormal(p0, p1, pp, nvec);
if (vtkPlanesIntersection::GoodNormal(nvec))
{
break;
}
}
normal->SetTuple(0, nvec);
pi->SetPoints(origin);
pi->SetNormals(normal);
origin->Delete();
normal->Delete();
pi->SetPlaneEquations();
int where = pi->EvaluateFacePlane(0, Box);
if (where != Straddle)
{
intersects = 0;
}
}
if (intersects == -1)
{
// 4. The Box intersects the plane of the polygon.
//
// For each of three orthographic projections (X, Y and Z),
// compute the equations of the edge lines of the polygon in those views.
//
// If Box's projection lies outside any of these projections, they
// don't intersect in 3D. Otherwise they do intersect in 3D.
//
// KDM: I'm pretty sure the above statement is untrue. I can think of a
// situation where all 3 projections intersect, but the 3D intersection
// does not. However, if the two intersect in 3D, then they will
// intersect in the 3 2D projections. Since I'm not worried about
// false positives, I'm not going to fix this right now.
if ((pi->IntersectsProjection(Box, Xdim) == 0) ||
(pi->IntersectsProjection(Box, Ydim) == 0) ||
(pi->IntersectsProjection(Box, Zdim) == 0) )
{
intersects = 0;
}
else
{
intersects = 1;
}
}
Box->Delete();
pi->Delete();
return intersects;
}
//---------------------------------------------------------------------
// Some convenience functions that build a vtkPlanesIntersection object
// out of a convex region.
//---------------------------------------------------------------------
// a static convenience function that converts a 3D cell into a
// vtkPlanesIntersection object
vtkPlanesIntersection *vtkPlanesIntersection::Convert3DCell(vtkCell *cell)
{
int i;
int nfaces = cell->GetNumberOfFaces();
vtkPoints *origins = vtkPoints::New();
origins->SetNumberOfPoints(nfaces);
vtkFloatArray *normals = vtkFloatArray::New();
normals->SetNumberOfComponents(3);
normals->SetNumberOfTuples(nfaces);
double inside[3] = {0.0, 0.0, 0.0};
for (i=0; i < nfaces; i++)
{
vtkCell *face = cell->GetFace(i);
vtkPoints *facePts = face->GetPoints();
int npts = facePts->GetNumberOfPoints();
double p0[3], p1[3], pp[3], n[3];
facePts->GetPoint(0, p0);
facePts->GetPoint(1, p1);
for (int p = 2; p < npts; p++)
{
facePts->GetPoint(p, pp);
vtkPlanesIntersection::ComputeNormal(pp, p1, p0, n);
if (vtkPlanesIntersection::GoodNormal(n))
{
break;
}
}
origins->SetPoint(i, pp);
normals->SetTuple(i, n);
inside[0] += p1[0];
inside[1] += p1[1];
inside[2] += p1[2];
}
inside[0] /= static_cast<double>(nfaces);
inside[1] /= static_cast<double>(nfaces);
inside[2] /= static_cast<double>(nfaces);
// ensure that all normals are outward pointing
for (i=0; i < nfaces; i++)
{
double ns[3], xs[3];
double n[3], x[3], p[4];
normals->GetTuple(i, ns);
origins->GetPoint(i, xs);
n[0] = ns[0];
x[0] = xs[0];
n[1] = ns[1];
x[1] = xs[1];
n[2] = ns[2];
x[2] = xs[2];
double outside[3];
outside[0] = x[0] + n[0];
outside[1] = x[1] + n[1];
outside[2] = x[2] + n[2];
vtkPlanesIntersection::PlaneEquation(n, x, p);
double insideVal = vtkPlanesIntersection::EvaluatePlaneEquation(inside, p);
double normalDirection =
vtkPlanesIntersection::EvaluatePlaneEquation(outside, p);
int sameSide =
( (insideVal < 0) && (normalDirection < 0)) ||
( (insideVal > 0) && (normalDirection > 0));
if (sameSide)
{
ns[0] = -ns[0];
ns[1] = -ns[1];
ns[2] = -ns[2];
normals->SetTuple(i, ns);
}
}
vtkPlanesIntersection *pi = vtkPlanesIntersection::New();
pi->SetPoints(origins);
pi->SetNormals(normals);
origins->Delete();
normals->Delete();
pi->SetRegionVertices(cell->GetPoints());
return pi;
}
//--------------------------------------------------------------------------
void vtkPlanesIntersection::ComputeNormal(double *p1, double *p2, double *p3,
double normal[3])
{
double v1[3], v2[3];
v1[0] = p1[0] - p2[0]; v1[1] = p1[1] - p2[1]; v1[2] = p1[2] - p2[2];
v2[0] = p3[0] - p2[0]; v2[1] = p3[1] - p2[1]; v2[2] = p3[2] - p2[2];
vtkMath::Cross(v1, v2, normal);
return;
}
int vtkPlanesIntersection::GoodNormal(double *n)
{
if ( (n[0] < VTK_SMALL_DOUBLE) || (n[0] > VTK_SMALL_DOUBLE) ||
(n[1] < VTK_SMALL_DOUBLE) || (n[1] > VTK_SMALL_DOUBLE) ||
(n[2] < VTK_SMALL_DOUBLE) || (n[2] > VTK_SMALL_DOUBLE) )
{
return 1;
}
else
{
return 0;
}
}
double vtkPlanesIntersection::EvaluatePlaneEquation(double *x, double *p)
{
return (x[0]*p[0] + x[1]*p[1] + x[2]*p[2] + p[3]);
}
void vtkPlanesIntersection::PlaneEquation(double *n, double *x, double *p)
{
p[0] = n[0];
p[1] = n[1];
p[2] = n[2];
p[3] = -(n[0]*x[0] + n[1]*x[1] + n[2]*x[2]);
}
// The plane equations ***********************************************
void vtkPlanesIntersection::SetPlaneEquations()
{
int i;
int nplanes = this->GetNumberOfPlanes();
// vtkPlanes stores normals & pts instead of
// plane equation coefficients
if (this->Plane) delete [] this->Plane;
this->Plane = new double[nplanes*4];
for (i=0; i<nplanes; i++)
{
double n[3], x[3];
this->Points->GetPoint(i, x);
this->Normals->GetTuple(i, n);
double nd[3], xd[3];
nd[0] = n[0]; xd[0] = x[0];
nd[1] = n[1]; xd[1] = x[1];
nd[2] = n[2]; xd[2] = x[2];
double *p = this->Plane + (i*4);
vtkPlanesIntersection::PlaneEquation(nd, xd, p);
}
}
// Compute region vertices if not set explicity ********************
void vtkPlanesIntersection::ComputeRegionVertices()
{
double M[3][3];
double rhs[3];
double testv[3];
int i, j, k;
int nplanes = this->GetNumberOfPlanes();
if (this->regionPts) this->regionPts->Delete();
this->regionPts = vtkPointsProjectedHull::New();
if (nplanes <= 3)
{
vtkErrorMacro( <<
"vtkPlanesIntersection::ComputeRegionVertices invalid region");
return;
}
if (this->Plane == NULL)
{
this->SetPlaneEquations();
}
// This is an expensive process. Better if vertices are
// set in SetRegionVertices(). We're testing every triple of
// planes to see if they intersect in a point that is
// not "outside" any plane.
int nvertices=0;
for (i=0; i < nplanes; i++)
{
for (j=i+1; j < nplanes; j++)
{
for (k=j+1; k < nplanes; k++)
{
this->planesMatrix(i, j, k, M);
int notInvertible = this->Invert3x3(M);
if (notInvertible) continue;
this->planesRHS(i, j, k, rhs);
vtkMath::Multiply3x3(M, rhs, testv);
if (duplicate(testv)) continue;
int outside = this->outsideRegion(testv);
if (!outside)
{
this->regionPts->InsertPoint(nvertices, testv);
nvertices++;
}
}
}
}
}
int vtkPlanesIntersection::duplicate(double testv[3]) const
{
int i;
double pt[3];
int npts = this->regionPts->GetNumberOfPoints();
for (i=0; i<npts; i++)
{
this->regionPts->GetPoint(i, pt);
if ( (pt[0] == testv[0]) && (pt[1] == testv[1]) && (pt[2] == testv[2]))
{
return 1;
}
}
return 0;
}
void vtkPlanesIntersection::planesMatrix(int p1, int p2, int p3, double M[3][3]) const
{
int i;
for (i=0; i<3; i++)
{
M[0][i] = this->Plane[p1*4 + i];
M[1][i] = this->Plane[p2*4 + i];
M[2][i] = this->Plane[p3*4 + i];
}
}
void vtkPlanesIntersection::planesRHS(int p1, int p2, int p3, double r[3]) const
{
r[0] = -(this->Plane[p1*4 + 3]);
r[1] = -(this->Plane[p2*4 + 3]);
r[2] = -(this->Plane[p3*4 + 3]);
}
int vtkPlanesIntersection::outsideRegion(double testv[3])
{
int i;
int outside = 0;
int nplanes = this->GetNumberOfPlanes();
for (i=0; i<nplanes; i++)
{
int row=i*4;
double fx =
vtkPlanesIntersection::EvaluatePlaneEquation(testv, this->Plane + row);
if (fx > VTK_SMALL_DOUBLE)
{
outside = 1;
break;
}
}
return outside;
}
int vtkPlanesIntersection::Invert3x3(double M[3][3])
{
int i, j;
double temp[3][3];
double det = vtkMath::Determinant3x3(M);
if ( (det > -VTK_SMALL_DOUBLE) && (det < VTK_SMALL_DOUBLE)) return -1;
vtkMath::Invert3x3(M, temp);
for (i=0; i<3; i++)
{
for (j=0; j<3; j++)
{
M[i][j] = temp[i][j];
}
}
return 0;
}
// Region / box intersection tests *******************************
int vtkPlanesIntersection::IntersectsBoundingBox(vtkPoints *R)
{
double BoxBounds[6], RegionBounds[6];
R->GetBounds(BoxBounds);
this->regionPts->GetBounds(RegionBounds);
if ((BoxBounds[1] < RegionBounds[0]) ||
(BoxBounds[0] > RegionBounds[1]) ||
(BoxBounds[3] < RegionBounds[2]) ||
(BoxBounds[2] > RegionBounds[3]) ||
(BoxBounds[5] < RegionBounds[4]) ||
(BoxBounds[4] > RegionBounds[5]))
{
return 0;
}
return 1;
}
int vtkPlanesIntersection::EnclosesBoundingBox(vtkPoints *R)
{
double BoxBounds[6], RegionBounds[6];
R->GetBounds(BoxBounds);
this->regionPts->GetBounds(RegionBounds);
if ((BoxBounds[0] > RegionBounds[0]) ||
(BoxBounds[1] < RegionBounds[1]) ||
(BoxBounds[2] > RegionBounds[2]) ||
(BoxBounds[3] < RegionBounds[3]) ||
(BoxBounds[4] > RegionBounds[4]) ||
(BoxBounds[5] < RegionBounds[5]))
{
return 0;
}
return 1;
}
int vtkPlanesIntersection::EvaluateFacePlane(int plane, vtkPoints *R)
{
int i;
double n[3], bounds[6];
double withN[3], oppositeN[3];
R->GetBounds(bounds);
this->Normals->GetTuple(plane, n);
// Find vertex of R most in direction of normal, and find
// oppposite vertex
for (i=0; i<3; i++)
{
if (n[i] < 0)
{
withN[i] = bounds[i*2];
oppositeN[i] = bounds[i*2 + 1];
}
else
{
withN[i] = bounds[i*2 + 1];
oppositeN[i] = bounds[i*2];
}
}
// Determine whether R is in negative half plane ("inside" frustum),
// positive half plane, or whether it straddles the plane.
// The normal points in direction of positive half plane.
double *p = this->Plane + (plane * 4);
double negVal =
vtkPlanesIntersection::EvaluatePlaneEquation(oppositeN, p);
if (negVal > 0)
{
return Outside;
}
double posVal =
vtkPlanesIntersection::EvaluatePlaneEquation(withN, p);
if (posVal < 0)
{
return Inside;
}
else return Straddle;
}
int vtkPlanesIntersection::IntersectsProjection(vtkPoints *R, int dir)
{
int intersects = 0;
switch (dir)
{
case Xdim:
intersects = this->regionPts->RectangleIntersectionX(R);
break;
case Ydim:
intersects = this->regionPts->RectangleIntersectionY(R);
break;
case Zdim:
intersects = this->regionPts->RectangleIntersectionZ(R);
break;
}
return intersects;
}
void vtkPlanesIntersection::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Plane: " << this->Plane << endl;
os << indent << "regionPts: " << this->regionPts << endl;
int i, npts;
if (this->Points)
{
npts = this->Points->GetNumberOfPoints();
for (i=0; i<npts; i++)
{
double *pt = this->Points->GetPoint(i);
double *n = this->Normals->GetTuple(i);
os << indent << "Origin " << pt[0] << " " << pt[1] << " " << pt[2] << " " ;
os << indent << "Normal " << n[0] << " " << n[1] << " " << n[2] << endl;
}
}
if (this->regionPts)
{
npts = this->regionPts->GetNumberOfPoints();
for (i=0; i<npts; i++)
{
double *pt = this->regionPts->GetPoint(i);
os << indent << "Vertex " << pt[0] << " " << pt[1] << " " << pt[2] << endl;
}
}
}