Attached Files | 14042-polydatatoimagestencil-robustness-3.diff [^] (21,608 bytes) 2014-09-15 11:50 [Show Content] [Hide Content]diff --git a/Imaging/Stencil/vtkPolyDataToImageStencil.cxx b/Imaging/Stencil/vtkPolyDataToImageStencil.cxx
index b7cf8c8..99a4ea0 100644
--- a/Imaging/Stencil/vtkPolyDataToImageStencil.cxx
+++ b/Imaging/Stencil/vtkPolyDataToImageStencil.cxx
@@ -53,7 +53,6 @@ POSSIBILITY OF SUCH DAMAGES.
#include "vtkCellArray.h"
#include "vtkDoubleArray.h"
#include "vtkSignedCharArray.h"
-#include "vtkMergePoints.h"
#include "vtkPoints.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
@@ -64,6 +63,8 @@ POSSIBILITY OF SUCH DAMAGES.
#include "vtkInformationVector.h"
#include "vtkStreamingDemandDrivenPipeline.h"
+#include <map>
+#include <utility>
#include <math.h>
@@ -109,10 +110,191 @@ void vtkPolyDataToImageStencil::PrintSelf(ostream& os,
}
//----------------------------------------------------------------------------
+// A helper class to quickly locate an edge, given the endpoint ids.
+// It uses an stl map rather than a table partitioning scheme, since
+// we have no idea how many entries there will be when we start. So
+// the performance is approximately log(n).
+namespace {
+
+// A Node in a linked list that contains information about one edge
+class EdgeLocatorNode
+{
+public:
+ EdgeLocatorNode() :
+ ptId0(-1), ptId1(-1), edgeId(-1), next(0) {}
+
+ // Free the list that this node is the head of
+ void FreeList() {
+ EdgeLocatorNode *ptr = this->next;
+ while (ptr)
+ {
+ EdgeLocatorNode *tmp = ptr;
+ ptr = ptr->next;
+ tmp->next = 0;
+ delete tmp;
+ }
+ }
+
+ vtkIdType ptId0;
+ vtkIdType ptId1;
+ vtkIdType edgeId;
+ EdgeLocatorNode *next;
+};
+
+// The EdgeLocator class itself, for keeping track of edges
+class EdgeLocator
+{
+private:
+ typedef std::map<vtkIdType, EdgeLocatorNode> MapType;
+ MapType EdgeMap;
+
+public:
+ EdgeLocator() : EdgeMap() {}
+ ~EdgeLocator() { this->Initialize(); }
+
+ // Description:
+ // Initialize the locator.
+ void Initialize();
+
+ // Description:
+ // If edge (i0, i1) is not in the list, then it will be added and
+ // a pointer for storing the new edgeId will be returned.
+ // If edge (i0, i1) is in the list, then edgeId will be set to the
+ // stored value and a null pointer will be returned.
+ vtkIdType *InsertUniqueEdge(vtkIdType i0, vtkIdType i1, vtkIdType &edgeId);
+
+ // Description:
+ // A helper function for interpolating a new point along an edge. It
+ // stores the index of the interpolated point in "i", and returns 1 if
+ // a new point was added to the locator. The values i0, i1, v0, v1 are
+ // the edge endpoints and scalar values, respectively.
+ int InterpolateEdge(
+ vtkPoints *inPoints, vtkPoints *outPoints,
+ vtkIdType i0, vtkIdType i1, double v0, double v1,
+ vtkIdType &i);
+};
+
+void EdgeLocator::Initialize()
+{
+ for (MapType::iterator i = this->EdgeMap.begin();
+ i != this->EdgeMap.end();
+ ++i)
+ {
+ i->second.FreeList();
+ }
+ this->EdgeMap.clear();
+}
+
+vtkIdType *EdgeLocator::InsertUniqueEdge(
+ vtkIdType i0, vtkIdType i1, vtkIdType &edgeId)
+{
+ // Ensure consistent ordering of edge
+ if (i1 < i0)
+ {
+ vtkIdType tmp = i0;
+ i0 = i1;
+ i1 = tmp;
+ }
+
+ // Generate a integer key, try to make it unique
+#ifdef VTK_USE_64BIT_IDS
+ vtkIdType key = ((i1 << 32) ^ i0);
+#else
+ vtkIdType key = ((i1 << 16) ^ i0);
+#endif
+
+ EdgeLocatorNode *node = &this->EdgeMap[key];
+
+ if (node->ptId1 < 0)
+ {
+ // Didn't find key, so add a new edge entry
+ node->ptId0 = i0;
+ node->ptId1 = i1;
+ return &node->edgeId;
+ }
+
+ // Search through the list for i0 and i1
+ if (node->ptId0 == i0 && node->ptId1 == i1)
+ {
+ edgeId = node->edgeId;
+ return 0;
+ }
+
+ int i = 1;
+ while (node->next != 0)
+ {
+ i++;
+ node = node->next;
+
+ if (node->ptId0 == i0 && node->ptId1 == i1)
+ {
+ edgeId = node->edgeId;
+ return 0;
+ }
+ }
+
+ // No entry for i1, so make one and return
+ node->next = new EdgeLocatorNode;
+ node = node->next;
+ node->ptId0 = i0;
+ node->ptId1 = i1;
+ node->edgeId = this->EdgeMap.size()-1;
+ return &node->edgeId;
+}
+
+int EdgeLocator::InterpolateEdge(
+ vtkPoints *points, vtkPoints *outPoints,
+ vtkIdType i0, vtkIdType i1, double v0, double v1,
+ vtkIdType &i)
+{
+ // This swap guarantees that exactly the same point is computed
+ // for both line directions, as long as the endpoints are the same.
+ if (v1 > 0)
+ {
+ vtkIdType tmpi = i0;
+ i0 = i1;
+ i1 = tmpi;
+
+ double tmp = v0;
+ v0 = v1;
+ v1 = tmp;
+ }
+
+ // After the above swap, i0 will be kept, and i1 will be clipped
+
+ // Check to see if this point has already been computed
+ vtkIdType *iptr = this->InsertUniqueEdge(i0, i1, i);
+ if (iptr == 0)
+ {
+ return 0;
+ }
+
+ // Get the edge and interpolate the new point
+ double p0[3], p1[3], p[3];
+ points->GetPoint(i0, p0);
+ points->GetPoint(i1, p1);
+
+ double f = v0/(v0 - v1);
+ double s = 1.0 - f;
+ double t = 1.0 - s;
+
+ p[0] = s*p0[0] + t*p1[0];
+ p[1] = s*p0[1] + t*p1[1];
+ p[2] = s*p0[2] + t*p1[2];
+
+ // Add the point, store the new index in the locator
+ i = outPoints->InsertNextPoint(p);
+ *iptr = i;
+
+ return 1;
+}
+
+} // end anonymous namespace
+
+//----------------------------------------------------------------------------
// Select contours within slice z
void vtkPolyDataToImageStencil::PolyDataSelector(
- vtkPolyData *input, vtkPolyData *output, double z, double thickness,
- vtkMergePoints *locator)
+ vtkPolyData *input, vtkPolyData *output, double z, double thickness)
{
vtkPoints *points = input->GetPoints();
vtkCellArray *lines = input->GetLines();
@@ -124,11 +306,8 @@ void vtkPolyDataToImageStencil::PolyDataSelector(
double minz = z - 0.5*thickness;
double maxz = z + 0.5*thickness;
- double bounds[6];
- input->GetBounds(bounds);
- bounds[4] = minz;
- bounds[5] = maxz;
- locator->InitPointInsertion(newPoints, bounds);
+ // use a map to avoid adding duplicate points
+ std::map<vtkIdType, vtkIdType> pointLocator;
vtkIdType loc = 0;
vtkIdType numCells = lines->GetNumberOfCells();
@@ -155,10 +334,21 @@ void vtkPolyDataToImageStencil::PolyDataSelector(
newLines->InsertNextCell(npts);
for (i = 0; i < npts; i++)
{
- vtkIdType ptId;
- double point[3];
- points->GetPoint(ptIds[i], point);
- locator->InsertUniquePoint(point, ptId);
+ vtkIdType oldId = ptIds[i];
+ std::map<vtkIdType, vtkIdType>::iterator iter =
+ pointLocator.lower_bound(oldId);
+ vtkIdType ptId = 0;
+ if (iter == pointLocator.end() || iter->first != oldId)
+ {
+ double point[3];
+ points->GetPoint(oldId, point);
+ ptId = newPoints->InsertNextPoint(point);
+ pointLocator.insert(iter, std::make_pair(oldId, ptId));
+ }
+ else
+ {
+ ptId = iter->second;
+ }
newLines->InsertCellPoint(ptId);
}
}
@@ -172,75 +362,100 @@ void vtkPolyDataToImageStencil::PolyDataSelector(
//----------------------------------------------------------------------------
// This method was taken from vtkCutter and slightly modified
void vtkPolyDataToImageStencil::PolyDataCutter(
- vtkPolyData *input, vtkPolyData *output, double z,
- vtkMergePoints *locator)
+ vtkPolyData *input, vtkPolyData *output, double z)
{
- vtkCellData *inCD = input->GetCellData();
- vtkCellData *outCD = output->GetCellData();
- vtkDoubleArray *cellScalars = vtkDoubleArray::New();
-
- // For the new points and lines
+ vtkPoints *points = input->GetPoints();
+ vtkCellArray *inputPolys = input->GetPolys();
+ vtkCellArray *inputStrips = input->GetStrips();
vtkPoints *newPoints = vtkPoints::New();
newPoints->Allocate(333);
vtkCellArray *newLines = vtkCellArray::New();
newLines->Allocate(1000);
- // No verts or polys are expected
- vtkCellArray *newVerts = vtkCellArray::New();
- vtkCellArray *newPolys = vtkCellArray::New();
-
- // Allocate space for the cell data
- outCD->CopyAllocate(inCD, 1000);
+ // An edge locator to avoid point duplication while clipping
+ EdgeLocator edgeLocator;
- // locator used to merge potentially duplicate points
- locator->InitPointInsertion(newPoints, input->GetBounds());
+ // Go through all cells and clip them.
+ vtkIdType numPolys = input->GetNumberOfPolys();
+ vtkIdType numStrips = input->GetNumberOfStrips();
+ vtkIdType numCells = numPolys + numStrips;
- // Compute some information for progress methods
- vtkGenericCell *cell = vtkGenericCell::New();
-
- // Loop over all cells; get scalar values for all cell points
- // and process each cell.
- vtkIdType numCells = input->GetNumberOfCells();
+ vtkIdType loc = 0;
+ vtkCellArray *cellArray = inputPolys;
for (vtkIdType cellId = 0; cellId < numCells; cellId++)
{
- input->GetCell(cellId, cell);
- vtkPoints *cellPts = cell->GetPoints();
- vtkIdList *cellIds = cell->GetPointIds();
+ // switch to strips when polys are done
+ if (cellId == numPolys)
+ {
+ loc = 0;
+ cellArray = inputStrips;
+ }
- if (cell->GetCellDimension() == 2)
+ vtkIdType npts, *ptIds;
+ cellArray->GetCell(loc, npts, ptIds);
+ loc += npts + 1;
+
+ vtkIdType numSubCells = 1;
+ if (cellArray == inputStrips)
+ {
+ numSubCells = npts - 2;
+ npts = 3;
+ }
+
+ for (vtkIdType subId = 0; subId < numSubCells; subId++)
{
- vtkIdType numCellPts = cellPts->GetNumberOfPoints();
- cellScalars->SetNumberOfTuples(numCellPts);
- for (vtkIdType i = 0; i < numCellPts; i++)
+ vtkIdType i1 = ptIds[npts-1];
+ double point[3];
+ points->GetPoint(i1, point);
+ double v1 = point[2] - z;
+ bool c1 = (v1 > 0);
+ bool odd = ((subId & 1) != 0);
+
+ // To store the ids of the contour line
+ vtkIdType linePts[2];
+ linePts[0] = 0;
+ linePts[1] = 0;
+
+ for (vtkIdType i = 0; i < npts; i++)
+ {
+ // Save previous point info
+ vtkIdType i0 = i1;
+ double v0 = v1;
+ bool c0 = c1;
+
+ // Generate new point info
+ i1 = ptIds[i];
+ points->GetPoint(i1, point);
+ v1 = point[2] - z;
+ c1 = (v1 > 0);
+
+ // If at least one edge end point wasn't clipped
+ if ( (c0 | c1) )
+ {
+ // If only one end was clipped, interpolate new point
+ if ( (c0 ^ c1) )
+ {
+ edgeLocator.InterpolateEdge(
+ points, newPoints, i0, i1, v0, v1, linePts[c0 ^ odd]);
+ }
+ }
+ }
+
+ // Insert the contour line if one was created
+ if (linePts[0] != linePts[1])
{
- // scalar value is distance from the specified z plane
- cellScalars->SetValue(i, input->GetPoint(cellIds->GetId(i))[2]);
+ newLines->InsertNextCell(2, linePts);
}
- cell->Contour(z, cellScalars, locator,
- newVerts, newLines, newPolys, NULL, NULL,
- inCD, cellId, outCD);
+ // Increment to get to the next triangle, if cell is a strip
+ ptIds++;
}
}
- // Update ourselves. Because we don't know upfront how many verts, lines,
- // polys we've created, take care to reclaim memory.
- cell->Delete();
- cellScalars->Delete();
-
output->SetPoints(newPoints);
+ output->SetLines(newLines);
newPoints->Delete();
-
- if (newLines->GetNumberOfCells())
- {
- output->SetLines(newLines);
- }
newLines->Delete();
- newVerts->Delete();
- newPolys->Delete();
-
- //release any extra memory
- locator->Initialize();
}
//----------------------------------------------------------------------------
@@ -279,9 +494,6 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
// get the input data
vtkPolyData *input = this->GetInput();
- // the locator to use with the data
- vtkMergePoints *locator = vtkMergePoints::New();
-
// the output produced by cutting the polydata with the Z plane
vtkPolyData *slice = vtkPolyData::New();
@@ -312,12 +524,12 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
// Step 1: Cut the data into slices
if (input->GetNumberOfPolys() > 0 || input->GetNumberOfStrips() > 0)
{
- this->PolyDataCutter(input, slice, z, locator);
+ this->PolyDataCutter(input, slice, z);
}
else
{
// if no polys, select polylines instead
- this->PolyDataSelector(input, slice, z, spacing[2], locator);
+ this->PolyDataSelector(input, slice, z, spacing[2]);
}
if (!slice->GetNumberOfLines())
@@ -360,7 +572,7 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
lines->InitTraversal();
vtkIdType npts;
vtkIdType *pointIds;
- while( lines->GetNextCell(npts, pointIds) )
+ while (lines->GetNextCell(npts, pointIds))
{
int isClosed = 0;
if (pointIds[0] == pointIds[npts-1])
@@ -395,82 +607,114 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
if (numberOfNeighbors == 1)
{
// store the loose end
- looseEndIdList->InsertNextId( i );
- looseEndNeighborList->InsertNextId( neighborId );
+ looseEndIdList->InsertNextId(i);
+ looseEndNeighborList->InsertNextId(neighborId);
}
// mark inflection points
- inflectionPointList->InsertNextValue( bottomPoint | topPoint );
+ inflectionPointList->InsertNextValue(bottomPoint | topPoint);
}
+ // join any loose ends
while (looseEndIdList->GetNumberOfIds() >= 2)
{
- // first loose end point in the list
- vtkIdType firstLooseEndId = looseEndIdList->GetId(0);
- vtkIdType neighborId = looseEndNeighborList->GetId(0);
- double firstLooseEnd[3];
- slice->GetPoint( firstLooseEndId, firstLooseEnd );
- double neighbor[3];
- slice->GetPoint( neighborId, neighbor);
-
- // second loose end in the list
- vtkIdType secondLooseEndId = looseEndIdList->GetId(1);
- double secondLooseEnd[3];
- slice->GetPoint( secondLooseEndId, secondLooseEnd );
+ vtkIdType n = looseEndIdList->GetNumberOfIds();
- // search for the loose end closest to the first one
+ // search for the two closest loose ends
double maxval = -VTK_FLOAT_MAX;
+ vtkIdType firstIndex = 0;
+ vtkIdType secondIndex = 1;
+ bool isCoincident = false;
- for(vtkIdType j = 1; j < looseEndIdList->GetNumberOfIds(); j++)
+ for (vtkIdType i = 0; i < n && !isCoincident; i++)
{
- vtkIdType currentLooseEndId = looseEndIdList->GetId( j );
- if (currentLooseEndId != neighborId)
+ // first loose end
+ vtkIdType firstLooseEndId = looseEndIdList->GetId(i);
+ vtkIdType neighborId = looseEndNeighborList->GetId(i);
+ double firstLooseEnd[3];
+ slice->GetPoint(firstLooseEndId, firstLooseEnd);
+ double neighbor[3];
+ slice->GetPoint(neighborId, neighbor);
+
+ for (vtkIdType j = i+1; j < n && !isCoincident; j++)
{
- double currentLooseEnd[3];
- slice->GetPoint( currentLooseEndId, currentLooseEnd );
-
- // When connecting loose ends, use dot product to favor
- // continuing in same direction as the line already
- // connected to the loose end, but also favour short
- // distances by dividing dotprod by square of distance.
- double v1[2], v2[2];
- v1[0] = firstLooseEnd[0] - neighbor[0];
- v1[1] = firstLooseEnd[1] - neighbor[1];
- v2[0] = currentLooseEnd[0] - firstLooseEnd[0];
- v2[1] = currentLooseEnd[1] - firstLooseEnd[1];
- double dotprod = v1[0]*v2[0] + v1[1]*v2[1];
- double distance2 = v2[0]*v2[0] + v2[1]*v2[1];
-
- if (dotprod > maxval*distance2 && distance2 > 0.0)
+ vtkIdType secondLooseEndId = looseEndIdList->GetId(j);
+ if (secondLooseEndId != neighborId)
{
- maxval = dotprod/distance2;
- secondLooseEndId = currentLooseEndId;
+ double currentLooseEnd[3];
+ slice->GetPoint(secondLooseEndId, currentLooseEnd);
+
+ // When connecting loose ends, use dot product to favor
+ // continuing in same direction as the line already
+ // connected to the loose end, but also favour short
+ // distances by dividing dotprod by square of distance.
+ double v1[2], v2[2];
+ v1[0] = firstLooseEnd[0] - neighbor[0];
+ v1[1] = firstLooseEnd[1] - neighbor[1];
+ v2[0] = currentLooseEnd[0] - firstLooseEnd[0];
+ v2[1] = currentLooseEnd[1] - firstLooseEnd[1];
+ double dotprod = v1[0]*v2[0] + v1[1]*v2[1];
+ double distance2 = v2[0]*v2[0] + v2[1]*v2[1];
+
+ if (distance2 == 0)
+ {
+ firstIndex = i;
+ secondIndex = j;
+ isCoincident = true;
+ }
+ else if (dotprod > maxval*distance2)
+ {
+ firstIndex = i;
+ secondIndex = j;
+ maxval = dotprod/distance2;
+ }
}
}
}
- // create a new line segment by connecting these two points
- looseEndIdList->DeleteId( firstLooseEndId );
- looseEndIdList->DeleteId( secondLooseEndId );
- looseEndNeighborList->DeleteId( firstLooseEndId );
- looseEndNeighborList->DeleteId( secondLooseEndId );
-
- lines->InsertNextCell( 2 );
- lines->InsertCellPoint( firstLooseEndId );
- lines->InsertCellPoint( secondLooseEndId );
+ // get info about the two loose ends and their neighbors
+ vtkIdType firstLooseEndId = looseEndIdList->GetId(firstIndex);
+ vtkIdType neighborId = looseEndNeighborList->GetId(firstIndex);
+ double firstLooseEnd[3];
+ slice->GetPoint(firstLooseEndId, firstLooseEnd);
+ double neighbor[3];
+ slice->GetPoint(neighborId, neighbor);
- // check if the new vertices are vertical inflection points
- slice->GetPoint( secondLooseEndId, secondLooseEnd );
- vtkIdType secondNeighborId = looseEndNeighborList->GetId(0);
+ vtkIdType secondLooseEndId = looseEndIdList->GetId(secondIndex);
+ vtkIdType secondNeighborId = looseEndNeighborList->GetId(secondIndex);
+ double secondLooseEnd[3];
+ slice->GetPoint(secondLooseEndId, secondLooseEnd);
double secondNeighbor[3];
- slice->GetPoint( secondNeighborId, secondNeighbor);
+ slice->GetPoint(secondNeighborId, secondNeighbor);
- inflectionPointList->SetValue(firstLooseEndId,
- ((firstLooseEnd[1] - neighbor[1])*
- (secondLooseEnd[1] - firstLooseEnd[1]) <= 0));
+ // remove these loose ends from the list
+ looseEndIdList->DeleteId(firstLooseEndId);
+ looseEndIdList->DeleteId(secondLooseEndId);
+ looseEndNeighborList->DeleteId(firstLooseEndId);
+ looseEndNeighborList->DeleteId(secondLooseEndId);
- inflectionPointList->SetValue(secondLooseEndId,
- ((secondLooseEnd[1] - firstLooseEnd[1])*
- (secondNeighbor[1] - secondLooseEnd[1]) <= 0));
+ if (isCoincident)
+ {
+ // check if the new vertex is a vertical inflection points
+ inflectionPointList->SetValue(firstLooseEndId,
+ ((firstLooseEnd[1] - neighbor[1])*
+ (secondNeighbor[1] - firstLooseEnd[1]) <= 0));
+ }
+ else
+ {
+ // create a new line segment by connecting these two points
+ lines->InsertNextCell(2);
+ lines->InsertCellPoint(firstLooseEndId);
+ lines->InsertCellPoint(secondLooseEndId);
+
+ // check if the new vertices are vertical inflection points
+ inflectionPointList->SetValue(firstLooseEndId,
+ ((firstLooseEnd[1] - neighbor[1])*
+ (secondLooseEnd[1] - firstLooseEnd[1]) <= 0));
+
+ inflectionPointList->SetValue(secondLooseEndId,
+ ((secondLooseEnd[1] - firstLooseEnd[1])*
+ (secondNeighbor[1] - secondLooseEnd[1]) <= 0));
+ }
}
// Step 3: Go through all the line segments for this slice,
@@ -480,7 +724,7 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
vtkIdType *pts = 0;
vtkIdType npts = 0;
- while ( lines->GetNextCell(npts, pts) )
+ while (lines->GetNextCell(npts, pts))
{
for (vtkIdType j = 1; j < npts; j++)
{
@@ -507,7 +751,6 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
}
slice->Delete();
- locator->Delete();
}
//----------------------------------------------------------------------------
diff --git a/Imaging/Stencil/vtkPolyDataToImageStencil.h b/Imaging/Stencil/vtkPolyDataToImageStencil.h
index 1b79281..751130f 100644
--- a/Imaging/Stencil/vtkPolyDataToImageStencil.h
+++ b/Imaging/Stencil/vtkPolyDataToImageStencil.h
@@ -94,11 +94,10 @@ protected:
int extent[6], int threadId);
static void PolyDataCutter(vtkPolyData *input, vtkPolyData *output,
- double z, vtkMergePoints *locator);
+ double z);
static void PolyDataSelector(vtkPolyData *input, vtkPolyData *output,
- double z, double thickness,
- vtkMergePoints *locator);
+ double z, double thickness);
virtual int RequestData(vtkInformation *, vtkInformationVector **,
vtkInformationVector *);
|