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 +#include #include @@ -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 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 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::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 *);