Commit 480a27dd authored by T.J. Corona's avatar T.J. Corona

Fix bugs in vtkMath, vtkLine, and add test for vtkLine static methods.

parent e04f7da4
......@@ -27,6 +27,7 @@
#include "vtkDataArray.h"
#include <cassert>
#include <cmath>
#include <limits>
#include <vector>
#include "vtkBoxMuellerRandomSequence.h"
......@@ -378,7 +379,9 @@ int vtkMath::SolveLinearSystem(double **A, double *x, int size)
det = vtkMath::Determinant2x2(A[0][0], A[0][1], A[1][0], A[1][1]);
if (det == 0.0)
static const double eps = 256*std::numeric_limits<double>::epsilon();
if (std::fabs(det) < eps)
{
// Unable to solve linear system
return 0;
......
......@@ -58,6 +58,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
UnitTestCells.cxx
UnitTestImplicitDataSet.cxx
UnitTestImplicitVolume.cxx
UnitTestLine.cxx
UnitTestPlanesIntersection.cxx
)
vtk_add_test_cxx(${vtk-module}CxxTests data_tests
......
This diff is collapsed.
......@@ -54,7 +54,7 @@ int vtkLine::EvaluatePosition(double x[3], double* closestPoint,
this->Points->GetPoint(0, a1);
this->Points->GetPoint(1, a2);
// DistanceToLine sets pcoords[0] to a value t, 0 <= t <= 1
// DistanceToLine sets pcoords[0] to a value t
dist2 = this->DistanceToLine(x,a1,a2,pcoords[0],closestPoint);
// pcoords[0] == t, need weights to be 1-t and t
......@@ -90,10 +90,11 @@ void vtkLine::EvaluateLocation(int& vtkNotUsed(subId), double pcoords[3],
}
//----------------------------------------------------------------------------
// Performs intersection of two finite 3D lines. An intersection is found if
// the projection of the two lines onto the plane perpendicular to the cross
// product of the two lines intersect. The parameters (u,v) are the
// parametric coordinates of the lines at the position of closest approach.
// Performs intersection of the projection of two finite 3D lines onto a 2D
// plane. An intersection is found if the projection of the two lines onto
// the plane perpendicular to the cross product of the two lines intersect.
// The parameters (u,v) are the parametric coordinates of the lines at the
// position of closest approach.
int vtkLine::Intersection (double a1[3], double a2[3],
double b1[3], double b2[3],
double& u, double& v)
......@@ -142,7 +143,7 @@ int vtkLine::Intersection (double a1[3], double a2[3],
{
minDist = dist;
*(uv1[i]) = t;
*(uv2[i]) = static_cast<double>(i%2);
*(uv2[i]) = static_cast<double>(i%2); // the corresponding extremum
}
}
return VTK_ON_LINE;
......@@ -164,6 +165,41 @@ int vtkLine::Intersection (double a1[3], double a2[3],
}
}
//----------------------------------------------------------------------------
int vtkLine::Intersection3D(double a1[3], double a2[3],
double b1[3], double b2[3],
double& u, double& v)
{
// Description:
// Performs intersection of two finite 3D lines. An intersection is found if
// the projection of the two lines onto the plane perpendicular to the cross
// product of the two lines intersect, and if the distance between the
// closest points of approach are within a relative tolerance. The parameters
// (u,v) are the parametric coordinates of the lines at the position of
// closest approach.
int projectedIntersection = vtkLine::Intersection(a1,a2,b1,b2,u,v);
if (projectedIntersection == VTK_YES_INTERSECTION)
{
double a_i,b_i;
double lenA = 0.;
double lenB = 0.;
double dist = 0.;
for (unsigned i=0;i<3;i++)
{
a_i = a1[i] + (a2[i]-a1[i])*u;
b_i = b1[i] + (b2[i]-b1[i])*v;
lenA += (a2[i]-a1[i])*(a2[i]-a1[i]);
lenB += (b2[i]-b1[i])*(b2[i]-b1[i]);
dist += (a_i-b_i)*(a_i-b_i);
}
if (dist > 1.e-6*(lenA > lenB ? lenA : lenB))
return VTK_NO_INTERSECTION;
}
return projectedIntersection;
}
//----------------------------------------------------------------------------
int vtkLine::CellBoundary(int vtkNotUsed(subId), double pcoords[3],
vtkIdList *pts)
......@@ -350,108 +386,33 @@ double vtkLine::DistanceBetweenLineSegments(
if (D < 1e-6)
{
// the lines are almost parallel. Where on the line is the closest point.
// First check if the lines overlap. If they do, then the distance is
// just the distance between the infinite lines..
double dist;
dist = vtkLine::DistanceToLine( l0, m0, m1, t2, closestPt2 );
if (t2 >= 0.0 && t2 <= 1.0)
{
t1 = 0.0;
closestPt1[0] = l0[0];
closestPt1[1] = l0[1];
closestPt1[2] = l0[2];
return dist;
}
dist = vtkLine::DistanceToLine( l1, m0, m1, t2, closestPt2 );
if (t2 >= 0.0 && t2 <= 1.0)
{
t1 = 1.0;
closestPt1[0] = l1[0];
closestPt1[1] = l1[1];
closestPt1[2] = l1[2];
return dist;
}
dist = vtkLine::DistanceToLine( m0, l0, l1, t1, closestPt1 );
if (t1 >= 0.0 && t1 <= 1.0)
{
t1 = 0.0;
closestPt2[0] = m0[0];
closestPt2[1] = m0[1];
closestPt2[2] = m0[2];
return dist;
}
dist = vtkLine::DistanceToLine( m1, l0, l1, t1, closestPt1 );
if (t1 >= 0.0 && t1 <= 1.0)
{
t1 = 1.0;
closestPt2[0] = m1[0];
closestPt2[1] = m1[1];
closestPt2[2] = m1[2];
return dist;
}
// The lines don't overlap.... The shortest distance is the min of the
// distance between the end points.
const double d1 = vtkMath::Distance2BetweenPoints( l0, m0 );
const double d2 = vtkMath::Distance2BetweenPoints( l0, m1 );
const double d3 = vtkMath::Distance2BetweenPoints( l1, m0 );
const double d4 = vtkMath::Distance2BetweenPoints( l1, m1 );
if (d1 <= d2 && d1 <= d3 && d1 <= d4)
{
t1 = t2 = 0.0;
closestPt1[0] = l0[0];
closestPt1[1] = l0[1];
closestPt1[2] = l0[2];
closestPt2[0] = m0[0];
closestPt2[1] = m0[1];
closestPt2[2] = m0[2];
return d1;
}
if (d2 <= d1 && d2 <= d3 && d2 <= d4)
{
t1 = 0.0; t2 = 1.0;
closestPt1[0] = l0[0];
closestPt1[1] = l0[1];
closestPt1[2] = l0[2];
closestPt2[0] = m1[0];
closestPt2[1] = m1[1];
closestPt2[2] = m1[2];
return d2;
}
if (d3 <= d1 && d3 <= d2 && d3 <= d4)
{
t1 = 1.0; t2 = 0.0;
closestPt1[0] = l1[0];
closestPt1[1] = l1[1];
closestPt1[2] = l1[2];
closestPt2[0] = m0[0];
closestPt2[1] = m0[1];
closestPt2[2] = m0[2];
return d3;
}
if (d4 <= d1 && d4 <= d2 && d4 <= d3)
// The lines are colinear. Therefore, one of the four endpoints is the
// point of closest approach
double minDist = VTK_DOUBLE_MAX;
double* p[4] = {l0,l1,m0,m1};
double* a1[4] = {m0,m0,l0,l0};
double* a2[4] = {m1,m1,l1,l1};
double* uv1[4] = {&t2,&t2,&t1,&t1};
double* uv2[4] = {&t1,&t1,&t2,&t2};
double* pn1[4] = {closestPt2,closestPt2,closestPt1,closestPt1};
double* pn2[4] = {closestPt1,closestPt1,closestPt2,closestPt2};
double dist,t,pn_[3];
for (unsigned i=0;i<4;i++)
{
t1 = t2 = 1.0;
closestPt1[0] = l1[0];
closestPt1[1] = l1[1];
closestPt1[2] = l1[2];
closestPt2[0] = m1[0];
closestPt2[1] = m1[1];
closestPt2[2] = m1[2];
return d4;
dist = vtkLine::DistanceToLine(p[i],a1[i],a2[i],t,pn_);
if (dist < minDist)
{
minDist = dist;
*(uv1[i]) = (t < 0. ? 0. : (t > 1. ? 1. : t));
*(uv2[i]) = static_cast<double>(i%2); // the corresponding extremum
for (unsigned j=0;j<3;j++)
{
pn1[i][j] = pn_[j];
pn2[i][j] = p[i][j];
}
}
}
// Can never get here.
return 0.0;
return minDist;
}
// The lines aren't parallel.
......
......@@ -75,15 +75,32 @@ public:
// Description:
// Performs intersection of two finite 3D lines. An intersection is found if
// the projection of the two lines onto the plane perpendicular to the cross
// product of the two lines intersect. The parameters (u,v) are the
// parametric coordinates of the lines at the position of closest approach.
// Performs intersection of the projection of two finite 3D lines onto a 2D
// plane. An intersection is found if the projection of the two lines onto
// the plane perpendicular to the cross product of the two lines intersect.
// The parameters (u,v) are the parametric coordinates of the lines at the
// position of closest approach.
static int Intersection(double p1[3], double p2[3],
double x1[3], double x2[3],
double& u, double& v);
// Description:
// Performs intersection of two finite 3D lines. An intersection is found if
// the projection of the two lines onto the plane perpendicular to the cross
// product of the two lines intersect, and if the distance between the
// closest points of approach are within a relative tolerance. The parameters
// (u,v) are the parametric coordinates of the lines at the position of
// closest approach.
//
// NOTE: "Unlike Intersection(), which determines whether the projections of
// two lines onto a plane intersect, Intersection3D() determines whether the
// lines themselves in 3D space intersect, within a tolerance.
static int Intersection3D(double p1[3], double p2[3],
double x1[3], double x2[3],
double& u, double& v);
// Description:
// Compute the distance of a point x to a finite line (p1,p2). The method
// computes the parametric coordinate t and the point location on the
......@@ -163,5 +180,3 @@ inline int vtkLine::GetParametricCenter(double pcoords[3])
}
#endif
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment