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

ENH: Added contouring to cells.

parent 6b0a073f
......@@ -26,6 +26,8 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
#include "Object.hh"
#include "FPoints.hh"
#include "FScalars.hh"
#include "CellArr.hh"
#include "IdList.hh"
class vlCell : public vlObject
......@@ -38,11 +40,9 @@ public:
// Because these objects (cells and derived classes) are computational
// objects, and because they are used internally, do not use memory
// reference counting.
// Point coordinates for cell.
vlFloatPoints *GetPoints() {return &this->Points;};
// Point ids for cell.
vlIdList *GetPointIds() {return &this->PointIds;};
// Number of points in cell
int NumberOfPoints() {return this->PointIds.NumberOfIds();};
// Dimensionality of cell (0,1,2, or 3)
virtual int CellDimension() = 0;
......@@ -55,12 +55,24 @@ public:
// Determine global coordinate from subId and parametric coordinates
virtual void EvaluateLocation(int& subId, float pcoords[3], float x[3]) = 0;
// Generate contouring primitives
virtual void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *scalars) = 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]);
// Point coordinates for cell.
vlFloatPoints *GetPoints() {return &this->Points;};
// Point ids for cell.
vlIdList *GetPointIds() {return &this->PointIds;};
// left public for quick computational access
vlFloatPoints Points;
vlIdList PointIds;
......
......@@ -21,19 +21,25 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
#include "DS2PolyF.hh"
#define MAX_CONTOURS 256
class vlContourFilter : public vlDataSetToPolyFilter
{
public:
vlContourFilter(float value=0.0) {this->Value = value;};
vlContourFilter();
~vlContourFilter() {};
char *GetClassName() {return "vlContourFilter";};
vlSetMacro(Value,float);
vlGetMacro(Value,float);
void SetValue(int i, float value);
vlGetVectorMacro(Values,float);
void GenerateValues(int numContours, float range[2]);
protected:
void Execute();
float Value;
float Values[MAX_CONTOURS];
int NumberOfContours;
float Range[2];
};
#endif
......
......@@ -30,6 +30,9 @@ public:
char *GetClassName() {return "vlHexahedron";};
int CellDimension() {return 3;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
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]);
......
......@@ -30,6 +30,10 @@ public:
char *GetClassName() {return "vlLine";};
int CellDimension() {return 1;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -30,6 +30,9 @@ public:
char *GetClassName() {return "vlRectangle";};
int CellDimension() {return 2;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -37,6 +37,10 @@ public:
int GenerateNormals(vlPoints *, vlCellArray *, vlFloatNormals *);
int CellDimension() {return 1;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -30,6 +30,9 @@ public:
char *GetClassName() {return "vlPolyPoints";};
int CellDimension() {return 0;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -37,6 +37,9 @@ public:
void ComputeNormal(vlFloatPoints *p, float n[3]);
int CellDimension() {return 2;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points,vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -30,6 +30,9 @@ public:
char *GetClassName() {return "vlQuad";};
int CellDimension() {return 2;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
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]);
......
......@@ -30,6 +30,9 @@ public:
char *GetClassName() {return "vlTetra";};
int CellDimension() {return 3;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -30,6 +30,10 @@ public:
char *GetClassName() {return "vlTriangleStrip";};
int CellDimension() {return 2;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -30,6 +30,10 @@ public:
char *GetClassName() {return "vlTriangle";};
int CellDimension() {return 2;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -30,6 +30,9 @@ public:
char *GetClassName() {return "vlPoint";};
int CellDimension() {return 0;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
float EvaluatePosition(float x[3], int& subId, float pcoords[3]);
void EvaluateLocation(int& subId, float pcoords[3], float x[3]);
......
......@@ -30,6 +30,9 @@ public:
char *GetClassName() {return "vlBrick";};
int CellDimension() {return 3;};
void Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys, vlFloatScalars *s);
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]);
......
......@@ -17,40 +17,127 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
#include "FScalars.hh"
#include "Cell.hh"
vlContourFilter::vlContourFilter()
{
for (int i=0; i<MAX_CONTOURS; i++) this->Values[i] = 0.0;
this->NumberOfContours = 1;
this->Range[0] = 0.0;
this->Range[1] = 1.0;
}
void vlContourFilter::SetValue(int i, float value)
{
i = (i >= MAX_CONTOURS ? MAX_CONTOURS-1 : (i < 0 ? 0 : i) );
if ( this->Values[i] != value )
{
this->Modified();
this->Values[i] = value;
if ( i > NumberOfContours ) this->NumberOfContours = i;
if ( value < this->Range[0] ) this->Range[0] = value;
if ( value > this->Range[1] ) this->Range[1] = value;
}
}
void vlContourFilter::GenerateValues(int numContours, float range[2])
{
float val, incr;
int i;
numContours = (numContours >= MAX_CONTOURS ? MAX_CONTOURS-1 :
(numContours < 0 ? 0 : numContours) );
if ( numContours < 1 ) numContours = 1;
incr = (range[1] - range[0]) / numContours;
for (i=0, val=range[0]; i < numContours; i++, val+=incr)
{
this->SetValue(i,val); // don't modify object unless absolutely nec.
}
}
//
// General contouring filter. Handles arbitrary input.
//
void vlContourFilter::Execute()
{
static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
int cellId, i;
int index;
vlIdList *cellPts;
vlScalars *inScalars;
vlFloatScalars cellScalars(8);
vlFloatScalars cellScalars(MAX_CELL_SIZE);
vlCell *cell;
float *range, value;
vlFloatScalars *newScalars;
vlCellArray *newVerts, *newLines, *newPolys;
vlFloatPoints *newPoints;
//
// Check that value is within scalar range
// Initialize and check input
//
this->Initialize();
if ( ! (inScalars = this->Input->GetPointData()->GetScalars()) )
{
vlErrorMacro(<<"No scalar data to contour");
return;
}
range = inScalars->GetRange();
//
// Create objects to hold output of contour operation
//
// Loop over all cells
newPoints = new vlFloatPoints(1000,1000);
newVerts = new vlCellArray(1000,1000);
newLines = new vlCellArray(1000,10000);
newPolys = new vlCellArray(1000,10000);
newScalars = new vlFloatScalars(3000,30000);
//
for (cellId=0; cellId<Input->NumberOfCells(); cellId++)
// Loop over all contour values. For each contour value, loop over all cells.
//
for (i=0; i<this->NumberOfContours; i++)
{
cell = Input->GetCell(cellId);
cellPts = cell->GetPointIds();
inScalars->GetScalars(*cellPts,cellScalars);
value = this->Values[i];
for (cellId=0; cellId<Input->NumberOfCells(); cellId++)
{
cell = Input->GetCell(cellId);
cellPts = cell->GetPointIds();
inScalars->GetScalars(*cellPts,cellScalars);
// Build the case table
for ( i=0, index = 0; i < cellPts->NumberOfIds(); i++)
if (cellScalars.GetScalar(i) >= this->Value)
index |= CASE_MASK[i];
cell->Contour(value, &cellScalars, newPoints, newVerts, newLines, newPolys, newScalars);
// Use different primitives to generate
} // for all cells
} // for all contour values
//
// Update ourselves. Because we don't know upfront how many verts, lines,
// polys we've created, take care to reclaim memory.
//
if (newVerts->GetNumberOfCells())
{
newVerts->Squeeze();
this->SetVerts(newVerts);
}
else
{
delete newVerts;
}
if (newLines->GetNumberOfCells())
{
newLines->Squeeze();
this->SetLines(newLines);
}
else
{
delete newLines;
}
} // for all cells
if (newPolys->GetNumberOfCells())
{
newPolys->Squeeze();
this->SetPolys(newPolys);
}
else
{
delete newPolys;
}
this->PointData.SetScalars(newScalars);
}
......@@ -18,6 +18,7 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
#include <math.h>
#include "Hexa.hh"
#include "vlMath.hh"
#include "Brick.hh"
//
// Note: the ordering of the Points and PointIds is important. See text.
......@@ -219,3 +220,48 @@ void vlHexahedron::EvaluateLocation(int& subId, float pcoords[3], float x[3])
}
}
}
//
// Marching cubes case table
//
#include "MC_Cases.h"
void vlHexahedron::Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points,
vlCellArray *verts, vlCellArray *lines,
vlCellArray *polys, vlFloatScalars *scalars)
{
static int CASE_MASK[8] = {1,2,4,8,16,32,64,128};
TRIANGLE_CASES *triCase;
EDGE_LIST *edge;
int i, j, index, *vert;
static int edges[12][2] = { {0,1}, {1,2}, {2,3}, {3,0},
{4,5}, {5,6}, {6,7}, {7,4},
{0,4}, {1,5}, {3,7}, {2,6}};
int pts[3];
float t, *x1, *x2, x[3];
// Build the case table
for ( i=0, index = 0; i < 8; i++)
if (cellScalars->GetScalar(i) >= value)
index |= CASE_MASK[i];
triCase = triCases + index;
edge = triCase->edges;
for ( ; edge[0] > -1; edge += 3 )
{
for (i=0; i<3; i++) // insert triangle
{
vert = edges[edge[i]];
t = (value - cellScalars->GetScalar(vert[0])) /
(cellScalars->GetScalar(vert[1]) - cellScalars->GetScalar(vert[0]));
x1 = this->Points.GetPoint(vert[0]);
x2 = this->Points.GetPoint(vert[1]);
for (j=0; j<3; j++) x[j] = x1[j] + t * (x2[j] - x1[j]);
pts[i] = points->InsertNextPoint(x);
scalars->InsertNextScalar(value);
}
polys->InsertNextCell(3,pts);
}
}
......@@ -143,3 +143,53 @@ int vlLine::Intersection (float a1[3], float a2[3], float b1[3], float b2[3],
return NO_INTERSECTION;
}
}
//
// marching lines case table
//
typedef int VERT_LIST;
typedef struct {
VERT_LIST verts[2];
} LINE_CASES;
static LINE_CASES lineCases[]= {
{-1,-1},
{1,0},
{0,1},
{-1,-1}};
void vlLine::Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points,
vlCellArray *verts, vlCellArray *lines,
vlCellArray *polys, vlFloatScalars *scalars)
{
static int CASE_MASK[2] = {1,2};
int index, i;
LINE_CASES *lineCase;
VERT_LIST *vert;
float t, x[3], *x1, *x2;
int pts[1];
//
// Build the case table
//
for ( i=0, index = 0; i < 2; i++)
if (cellScalars->GetScalar(i) >= value)
index |= CASE_MASK[i];
lineCase = lineCases + index;
vert = lineCase->verts;
while ( vert[0] > -1 )
{
t = (value - cellScalars->GetScalar(vert[0])) /
(cellScalars->GetScalar(vert[1]) - cellScalars->GetScalar(vert[0]));
x1 = this->Points.GetPoint(vert[0]);
x2 = this->Points.GetPoint(vert[1]);
for (i=0; i<3; i++) x[i] = x1[i] + t * (x2[i] - x1[i]);
pts[0] = points->InsertNextPoint(x);
verts->InsertNextCell(1,pts);
scalars->InsertNextScalar(value);
}
}
......@@ -95,3 +95,68 @@ void vlRectangle::EvaluateLocation(int& subId, float pcoords[3], float x[3])
pcoords[1]*(pt3[i] - pt1[i]);
}
}
//
// Marching (convex) quadrilaterals
//
typedef int EDGE_LIST;
typedef struct {
EDGE_LIST edges[5];
} LINE_CASES;
static LINE_CASES lineCases[] = {
{-1, -1, -1, -1, -1},
{0, 3, -1, -1, -1},
{1, 0, -1, -1, -1},
{1, 3, -1, -1, -1},
{2, 1, -1, -1, -1},
{0, 3, 2, 1, -1},
{2, 0, -1, -1, -1},
{2, 3, -1, -1, -1},
{3, 2, -1, -1, -1},
{0, 2, -1, -1, -1},
{1, 0, 3, 2, -1},
{1, 2, -1, -1, -1},
{3, 1, -1, -1, -1},
{0, 1, -1, -1, -1},
{3, 0, -1, -1, -1},
{-1, -1, -1, -1, -1}
};
void vlRectangle::Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *scalars)
{
static int CASE_MASK[4] = {1,2,4,8};
LINE_CASES *lineCase;
EDGE_LIST *edge;
int i, j, index, *vert;
static int edges[4][2] = { {0,1}, {1,2}, {2,3}, {3,0} };
int pts[2];
float t, *x1, *x2, x[3];
// Build the case table
for ( i=0, index = 0; i < 4; i++)
if (cellScalars->GetScalar(i) >= 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]];
t = (value - cellScalars->GetScalar(vert[0])) /
(cellScalars->GetScalar(vert[1]) - cellScalars->GetScalar(vert[0]));
x1 = this->Points.GetPoint(vert[0]);
x2 = this->Points.GetPoint(vert[1]);
for (j=0; j<3; j++) x[j] = x1[j] + t * (x2[j] - x1[j]);
pts[i] = points->InsertNextPoint(x);
scalars->InsertNextScalar(value);
}
lines->InsertNextCell(2,pts);
}
}
......@@ -85,9 +85,13 @@ int vlPolyLine::GenerateNormals(vlPoints *pts, vlCellArray *lines, vlFloatNormal
return 1;
}
//
// eliminate constructor / destructor calls
//
static vlLine line;
float vlPolyLine::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
vlLine line;
float pc[3], dist2, minDist2;
int ignoreId, i;
......@@ -120,3 +124,26 @@ void vlPolyLine::EvaluateLocation(int& subId, float pcoords[3], float x[3])
x[i] = a1[i] + pcoords[0]*(a2[i] - a1[i]);
}
}
void vlPolyLine::Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *scalars)
{
int i;
vlFloatScalars lineScalars(2);
for ( i=0; i<this->Points.NumberOfPoints()-1; i++)
{
line.Points.SetPoint(0,this->Points.GetPoint(i));
line.Points.SetPoint(1,this->Points.GetPoint(i+1));
lineScalars.SetScalar(0,cellScalars->GetScalar(i));
lineScalars.SetScalar(1,cellScalars->GetScalar(i+1));
line.Contour(value, &lineScalars, points, verts,
lines, polys, scalars);
}
}
......@@ -20,7 +20,7 @@ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen 1993, 1994
float vlPolyPoints::EvaluatePosition(float x[3], int& subId, float pcoords[3])
{
int numPts;
int numPts=this->Points.NumberOfPoints();
float *X;
float dist2, minDist2;
int i;
......@@ -56,3 +56,22 @@ void vlPolyPoints::EvaluateLocation(int& subId, float pcoords[3], float x[3])
x[1] = X[1];
x[2] = X[2];
}
void vlPolyPoints::Contour(float value, vlFloatScalars *cellScalars,
vlFloatPoints *points, vlCellArray *verts,
vlCellArray *lines, vlCellArray *polys,
vlFloatScalars *scalars)
{
int i, pts[1];
for (i=0; i<this->Points.NumberOfPoints(); i++)
{
if ( value == cellScalars-><