Commit c2337989 authored by Will Schroeder's avatar Will Schroeder
Browse files

ENH: Added computational guts to cells.

parent 01113b1c
......@@ -21,18 +21,34 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
#ifndef __vlCell_h
#define __vlCell_h
#define TOL 1.e-05 // Tolerance for geometric calculation
#include "Object.hh"
#include "FPoints.hh"
class vlCell : public vlObject
{
public:
vlCell() {};
vlCell():Points(0) {};
char *GetClassName() {return "vlCell";};
void SetPoints(vlFloatPoints *pts) {this->Points = pts;};
virtual float DistanceToPoint(float *x) = 0;
// return inside cell (!=0) if <= tolerance squared; evaluate parametric
// coordinates, sub-cell id (if cell is composite), and distance squared
// of point x[3] to cell.
virtual float EvaluatePosition(float x[3], int& subId, float pcoords[3]) = 0;
// virtual int EvaluatePosition(float x[3], float& tol2,
// int& subId, float pcoords[3], float& dist2) = 0;
// Determine global coordinate from subId and parametric coordinates
virtual void EvaluateLocation(int& subId, float pcoords[3], float x[3]) = 0;
// Compute cell bounding box (xmin,xmax,ymin,ymax,zmin,zmax)
float *GetBounds();
// Quick intersection of cell bounding box. Returns != 0 for hit.
char HitBBox(float bounds[6], float origin[3], float dir[3], float coord[3]);
protected:
vlFloatPoints *Points;
......
......@@ -29,8 +29,11 @@ public:
vlLine() {};
char *GetClassName() {return "vlLine";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
Intersection(float x[3], float xray[3], float x1[3], float x2[3],
float& u, float& v);
};
#endif
......
......@@ -36,7 +36,8 @@ public:
int GenerateNormals(vlPoints *, vlCellArray *, vlFloatNormals *);
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
};
......
/*=========================================================================
Program: Visualization Library
......@@ -29,7 +30,8 @@ public:
vlPolyPoints() {};
char *GetClassName() {return "vlPolyPoints";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
};
......
......@@ -32,10 +32,17 @@ public:
vlPolygon() {};
char *GetClassName() {return "vlPolygon";};
void ComputeNormal(vlPoints *p, int numPts, int *pts, float *n);
void ComputeNormal(float *v1, float *v2, float *v3, float *n);
void ComputeNormal(vlPoints *p, int numPts, int *pts, float n[3]);
void ComputeNormal(float v1[3], float v2[3], float v3[3], float n[3]);
void ComputeNormal(vlFloatPoints *p, float n[3]);
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
int ParameterizePolygon(float p0[3], float p10[3], float &l10,
float p20[3], float &l20, float n[3]);
int PointInPolygon(float bounds[6], float x[3], float n[3]);
};
#endif
......
......@@ -29,7 +29,10 @@ public:
vlQuad() {};
char *GetClassName() {return "vlQuad";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
void ShapeFunctions(float pcoords[3], float sf[4]);
void ShapeDerivs(float pcoords[3], float derivs[12]);
};
......
......@@ -29,7 +29,8 @@ public:
vlTetra() {};
char *GetClassName() {return "vlTetra";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
};
......
......@@ -29,7 +29,8 @@ public:
vlTriangleStrip() {};
char *GetClassName() {return "vlTriangleStrip";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
};
......
......@@ -29,7 +29,8 @@ public:
vlTriangle() {};
char *GetClassName() {return "vlTriangle";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
};
......
......@@ -29,7 +29,8 @@ public:
vlPoint() {};
char *GetClassName() {return "vlPoint";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
};
......
......@@ -29,7 +29,10 @@ public:
vlBrick() {};
char *GetClassName() {return "vlBrick";};
float DistanceToPoint(float *x);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
void ShapeFunctions(float pcoords[3], float sf[8]);
void ShapeDerivs(float pcoords[3], float derivs[24]);
};
......
......@@ -16,3 +16,109 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include "Cell.hh"
//
// Bounding box intersection modified from Graphics Gems Vol I.
// Note: the intersection ray is assumed normalized such that
// valid intersections can only occur between [0,1].
//
#define RIGHT 0
#define LEFT 1
#define MIDDLE 2
char vlCell::HitBBox (float bounds[6], float origin[3], float dir[3], float coord[3])
{
char inside=1;
char quadrant[3];
int i, whichPlane=0;
float maxT[3], candidatePlane[3];
//
// First find closest planes
//
for (i=0; i<3; i++)
{
if ( origin[i] < bounds[2*i] )
{
quadrant[i] = LEFT;
candidatePlane[i] = bounds[2*i];
inside = 0;
}
else if ( origin[i] > bounds[2*i+1] )
{
quadrant[i] = RIGHT;
candidatePlane[i] = bounds[2*i+1];
inside = 0;
}
else
{
quadrant[i] = MIDDLE;
}
}
//
// Check whether origin of ray is inside bbox
//
if (inside) return 1;
//
// Calculate parametric distances to plane
//
for (i=0; i<3; i++)
{
if ( quadrant[i] != MIDDLE && dir[i] != 0.0 )
{
maxT[i] = (candidatePlane[i]-origin[i]) / dir[i];
}
else
{
maxT[i] = -1.0;
}
}
//
// Find the largest parametric value of intersection
//
for (i=0; i<3; i++)
if ( maxT[whichPlane] < maxT[i] )
whichPlane = i;
//
// Check for valid intersection along line
//
if ( maxT[whichPlane] > 1.0 || maxT[whichPlane] < 0.0 )
return 0;
//
// Intersection point along line is okay. Check bbox.
//
for (i=0; i<3; i++)
{
if (whichPlane != i)
{
coord[i] = origin[i] + maxT[whichPlane]*dir[i];
if ( coord[i] < bounds[2*i] || coord[i] > bounds[2*i+1] )
return 0;
}
else
{
coord[i] = candidatePlane[i];
}
}
return 1;
}
float *vlCell::GetBounds ()
{
float *x;
int i, j;
static float bounds[6];
bounds[0] = bounds[2] = bounds[4] = LARGE_FLOAT;
bounds[1] = bounds[3] = bounds[5] = -LARGE_FLOAT;
for (i=0; i<this->Points->NumberOfPoints(); i++)
{
x = this->Points->GetPoint(i);
for (j=0; j<3; j++)
{
if ( x[j] < bounds[2*j] ) bounds[2*j] = x[j];
if ( x[j] > bounds[2*j+1] ) bounds[2*j+1] = x[j];
}
}
return bounds;
}
......@@ -16,9 +16,130 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include "Line.hh"
#include "vlMath.hh"
float vlLine::DistanceToPoint(float *x)
#define NO_INTERSECTION 1
#define INTERSECTION 2
#define ON_LINE 6
float vlLine::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
float *a1, *a2, a21[3], denom, num, *closestPoint;
int i, numPts;
vlMath math;
subId = 0;
pcoords[1] = pcoords[2] = 0.0;
a1 = this->Points->GetPoint(0);
a2 = this->Points->GetPoint(1);
//
// Determine appropriate vectors
//
for (i=0; i<3; i++) a21[i] = a2[i] - a1[i];
//
// Get parametric location
//
num = a21[0]*(x[0]-a1[0]) + a21[1]*(x[1]-a1[1]) + a21[2]*(x[2]-a1[2]);
denom = math.Dot(a21,a21);
if ( (denom = math.Dot(a21,a21)) < fabs(TOL*num) )
{
return LARGE_FLOAT;
}
else
{
pcoords[0] = num / denom;
}
//
// If parametric coordinate is within 0<=p<=1, then the point is closest to
// the line. Otherwise, it's closest to a point at the end of the line.
//
if ( pcoords[0] < 0.0 )
{
closestPoint = a1;
}
else if ( pcoords[0] > 1.0 )
{
closestPoint = a2;
}
else
{
closestPoint = a21;
for (i=0; i<3; i++) a21[i] = a1[i] + pcoords[0]*a21[i];
}
return math.Distance2BetweenPoints(closestPoint,x);
}
return 1.0;
void vlLine::EvaluateLocation(int& subId, float pcoords[3], float x[3])
{
int i;
float *a1 = this->Points->GetPoint(0);
float *a2 = this->Points->GetPoint(1);
for (i=0; i<3; i++)
{
x[i] = a1[i] + pcoords[0]*(a2[i] - a1[i]);
}
}
//
// Intersect two 3D lines
//
int vlLine::Intersection (float a1[3], float a2[3], float b1[3], float b2[3],
float& u, float& v)
{
float a21[3], b21[3], b1a1[3];
float sys[2][2], c[2], det;
int i;
vlMath math;
//
// Initialize
//
u = v = 0.0;
//
// Determine line vectors.
//
for (i=0; i<3; i++)
{
a21[i] = a2[i] - a1[i];
b21[i] = b2[i] - b1[i];
b1a1[i] = b1[i] - a1[i];
}
//
// Compute the system (least squares) matrix.
//
sys[0][0] = math.Dot ( a21, a21 );
sys[0][1] = -math.Dot ( a21, b21 );
sys[1][0] = sys[0][1];
sys[1][1] = math.Dot ( b21, b21 );
//
// Compute the least squares system constant term.
//
c[0] = math.Dot ( a21, b1a1 );
c[1] = -math.Dot ( b21, b1a1 );
//
// Solve the system of equations
//
if ( (det=math.Determinate2x2(sys[0],sys[1])) <= TOL )
{
return ON_LINE;
}
else
{
u = math.Determinate2x2(c,sys[1]) / det;
v = math.Determinate2x2(sys[0],c) / det;
}
//
// Check parametric coordinates for intersection.
//
if ( (0.0 <= u) && (u <= 1.0) && (0.0 <= v) && (v <= 1.0) )
{
return INTERSECTION;
}
else
{
return NO_INTERSECTION;
}
}
......@@ -17,6 +17,7 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include "PolyLine.hh"
#include "vlMath.hh"
#include "Line.hh"
int vlPolyLine::GenerateNormals(vlPoints *pts, vlCellArray *lines, vlFloatNormals *normals)
{
......@@ -84,8 +85,40 @@ int vlPolyLine::GenerateNormals(vlPoints *pts, vlCellArray *lines, vlFloatNormal
return 1;
}
float vlPolyLine::DistanceToPoint(float *x)
float vlPolyLine::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
vlLine line;
vlFloatPoints pts(2);
float pc[3], dist2, minDist2;
int ignoreId, i;
return 1.0;
pcoords[1] = pcoords[2] = 0.0;
line.SetPoints(&pts);
for (minDist2=LARGE_FLOAT,i=0; i<this->Points->NumberOfPoints()-1; i++)
{
pts.SetPoint(0,this->Points->GetPoint(i));
pts.SetPoint(1,this->Points->GetPoint(i+1));
dist2 = line.EvaluatePosition(x, ignoreId, pc);
if ( dist2 < minDist2 )
{
subId = i;
pcoords[0] = pc[0];
minDist2 = dist2;
}
}
return minDist2;
}
void vlPolyLine::EvaluateLocation(int& subId, float pcoords[3], float x[3])
{
int i;
float *a1 = this->Points->GetPoint(subId);
float *a2 = this->Points->GetPoint(subId+1);
for (i=0; i<3; i++)
{
x[i] = a1[i] + pcoords[0]*(a2[i] - a1[i]);
}
}
......@@ -16,9 +16,43 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include "PolyPts.hh"
#include "vlMath.hh"
float vlPolyPoints::DistanceToPoint(float *x)
float vlPolyPoints::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
int numPts;
float *X;
float dist2, minDist2;
int i;
vlMath math;
return 1.0;
for (minDist2=LARGE_FLOAT, i=0; i<numPts; i++)
{
X = this->Points->GetPoint(i);
dist2 = math.Distance2BetweenPoints(X,x);
if (dist2 < minDist2)
{
minDist2 = dist2;
subId = i;
}
}
if (minDist2 == 0.0)
{
pcoords[0] = 0.0;
}
else
{
pcoords[0] = -10.0;
}
return minDist2;
}
void vlPolyPoints::EvaluateLocation(int& subId, float pcoords[3], float x[3])
{
float *X = this->Points->GetPoint(subId);
x[0] = X[0];
x[1] = X[1];
x[2] = X[2];
}
......@@ -15,8 +15,17 @@ without the express written consent of the authors.
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
=========================================================================*/
#include <math.h>
#include "Polygon.hh"
#include "vlMath.hh"
#include "Line.hh"
#include "Plane.hh"
#define FAILURE 0
#define INTERSECTION 2
#define OUTSIDE 3
#define INSIDE 4
#define ON_LINE 6
void vlPolygon::ComputeNormal(vlPoints *p, int numPts, int *pts, float *n)
{
......@@ -80,8 +89,314 @@ void vlPolygon::ComputeNormal(float *v1, float *v2, float *v3, float *n)
}
}
float vlPolygon::DistanceToPoint(float *x)
void vlPolygon::ComputeNormal(vlFloatPoints *p, float *n)
{
int i, numPts;
float *v1, *v2, *v3;
float length;
float ax, ay, az;
float bx, by, bz;
//
// Because some polygon vertices are colinear, need to make sure
// first non-zero normal is found.
//
numPts = p->NumberOfPoints();
v1 = p->GetPoint(0);
v2 = p->GetPoint(1);
v3 = p->GetPoint(2);
for (i=0; i<numPts; i++)
{
ax = v2[0] - v1[0]; ay = v2[1] - v1[1]; az = v2[2] - v1[2];
bx = v3[0] - v1[0]; by = v3[1] - v1[1]; bz = v3[2] - v1[2];
n[0] = (ay * bz - az * by);
n[1] = (az * bx - ax * bz);
n[2] = (ax * by - ay * bx);
length = sqrt (n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
if (length != 0.0)
{
n[0] /= length;
n[1] /= length;
n[2] /= length;
return;
}
else
{
v1 = v2;
v2 = v3;
v3 = p->GetPoint((i+3)%numPts);
}
}
}
float vlPolygon::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
int i;
float p0[3], p10[3], l10, p20[3], l20, n[3];
float xproj[3], ray[3];
vlPlane plane;
vlMath math;
this->ParameterizePolygon(p0, p10, l10, p20, l20, n);
plane.ProjectPoint(x,p0,n,xproj);
for (i=0; i<3; i++) ray[i] = xproj[i] - p0[i];
pcoords[0] = math.Dot(ray,p10) / l10;
pcoords[1] = math.Dot(ray,p20) / l20;
if ( pcoords[0] >= 0.0 && pcoords[0] <= 1.0 &&
pcoords[1] >= 0.0 && pcoords[1] <= 1.0 &&
this->PointInPolygon(this->GetBounds(),xproj,n) == INSIDE )
{
return math.Distance2BetweenPoints(x,xproj);
}
//
// If here, point is outside of polygon, so need to find distance to boundary
//
vlLine line;
float pc[3], dist2, minDist2;
int ignoreId, numPts;
vlFloatPoints pts(2);
numPts = this->Points->NumberOfPoints();
line.SetPoints(&pts);
for (minDist2=LARGE_FLOAT,i=0; i<numPts - 1; i++)
{
pts.SetPoint(0,this->Points->GetPoint(i));
pts.SetPoint(1,this->Points->GetPoint(i+1));
dist2 = line.EvaluatePosition(x, ignoreId, pc);
if ( dist2 < minDist2 )
{