Commit 88f0ff89 authored by Will Schroeder's avatar Will Schroeder

ENH:Added new filters for filling holes in meshes, and for decimating polylines

parent ff94b0be
package require vtk
package require vtkinteraction
package require vtktesting
# Create the RenderWindow, Renderer and both Actors
#
vtkRenderer ren1
vtkRenderWindow renWin
renWin AddRenderer ren1
vtkRenderWindowInteractor iren
iren SetRenderWindow renWin
# Create some data, remove some polygons
vtkPolyData pd
vtkPoints pts
vtkCellArray polys
pd SetPoints pts
pd SetPolys polys
set xRes 10
set yRes 10
set xPtsRes [expr $xRes + 1]
set yPtsRes [expr $yRes + 1]
#insert points
for {set j 0} {$j < $yPtsRes} {incr j} {
for {set i 0} {$i < $xPtsRes} {incr i} {
pts InsertNextPoint $i $j 0.0
}
}
#insert cells
for {set j 1} {$j <= $yRes} {incr j} {
for {set i 1} {$i <= $xRes} {incr i} {
set cellId [expr $i -1 + $yRes*($j-1)]
if { $cellId != 48 &&\
$cellId != 12 && $cellId != 13 && $cellId != 23 &&\
$cellId != 60 &&\
$cellId != 83 && $cellId != 72 &&\
$cellId != 76 && $cellId != 77 && $cellId != 78 && $cellId != 87 } {
polys InsertNextCell 4
polys InsertCellPoint [expr $i - 1 + (($j-1)*$yPtsRes)]
polys InsertCellPoint [expr $i + (($j-1)*$yPtsRes)]
polys InsertCellPoint [expr $i + ($j*$yPtsRes)]
polys InsertCellPoint [expr $i - 1 + ($j*$yPtsRes)]
}
}
}
# Fill the holes
vtkFillHolesFilter fill
fill SetInput pd
fill SetHoleSize 20.0
# Mapping and actor
vtkPolyDataMapper map
map SetInput pd
# map SetInputConnection [fill GetOutputPort]
vtkActor actor
actor SetMapper map
[actor GetProperty] SetColor 1 0 0
# Add the actors to the renderer, set the background and size
#
ren1 AddActor actor
ren1 SetBackground 1 1 1
renWin SetSize 500 500
renWin Initialize
# render the image
#
iren AddObserver UserEvent {wm deiconify .vtkInteract}
# prevent the tk window from showing up then start the event loop
wm withdraw .
/*=========================================================================
Program: Visualization Toolkit
Module: vtkDecimatePolylineFilter.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkDecimatePolylineFilter.h"
#include "vtkDoubleArray.h"
#include "vtkLine.h"
#include "vtkMath.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkLine.h"
#include "vtkPolyData.h"
#include "vtkPoints.h"
#include "vtkCellArray.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include <vtkstd/list>
#include <vtkstd/vector>
#include <vtkstd/queue>
vtkCxxRevisionMacro(vtkDecimatePolylineFilter, "1.1");
vtkStandardNewMacro(vtkDecimatePolylineFilter);
//------------------------------------------------------------------
// STL list for implementing the algorithm
struct vtkPLineVertex
{
vtkIdType PtId;
double Error;
vtkPLineVertex(vtkIdType id, double error) : PtId(id), Error(error) {}
};
typedef vtkstd::list<vtkPLineVertex> PLineVertexListType;
typedef vtkstd::list<vtkPLineVertex>::iterator PLineVertexListIterator;
class CompareError {
public:
int operator()( const PLineVertexListIterator &x, const PLineVertexListIterator &y )
{ return x->Error > y->Error; }
};
typedef vtkstd::priority_queue<PLineVertexListIterator,vtkstd::vector<PLineVertexListIterator>,CompareError> PQueueType;
//---------------------------------------------------------------------
// Create object with specified reduction of 90%.
vtkDecimatePolylineFilter::vtkDecimatePolylineFilter()
{
this->TargetReduction = 0.90;
}
//---------------------------------------------------------------------
vtkDecimatePolylineFilter::~vtkDecimatePolylineFilter()
{
}
//---------------------------------------------------------------------
// Reduce the number of points in a set of polylines
//
int vtkDecimatePolylineFilter::RequestData(
vtkInformation *vtkNotUsed(request),
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
// get the info objects
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
vtkInformation *outInfo = outputVector->GetInformationObject(0);
// get the input and ouptut
vtkPolyData *input = vtkPolyData::SafeDownCast(
inInfo->Get(vtkDataObject::DATA_OBJECT()));
vtkPolyData *output = vtkPolyData::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
vtkCellArray *inputLines = input->GetLines();
vtkPoints * inputPoints = input->GetPoints();
vtkDebugMacro("Decimating polylines");
if (!inputLines || !inputPoints)
{
return 1;
}
vtkIdType numLines = inputLines->GetNumberOfCells();
vtkIdType numPts = inputPoints->GetNumberOfPoints();
if ( numLines < 1 || numPts < 1 )
{
return 1;
}
// Allocate memory and prepare for data processing
vtkPoints *newPts = vtkPoints::New();
vtkCellArray *newLines = vtkCellArray::New();
newLines->Allocate(numLines,2);
vtkPointData *inPD = input->GetPointData();
vtkPointData *outPD = output->GetPointData();
vtkCellData *inCD = input->GetCellData();
vtkCellData *outCD = output->GetCellData();
outPD->CopyAllocate(inPD);
outCD->CopyAllocate(inCD);
// Loop over all polylines, decimating each independently
vtkIdType npts, *pts, i, cellId, newId;
double x[3], x1[3], x2[3], error, len, dist;
PLineVertexListType VertList;
PLineVertexListIterator vertIterator;
PQueueType PQueue;
for (cellId=0, inputLines->InitTraversal(); inputLines->GetNextCell(npts,pts); cellId++)
{
if ( npts < 3 )
{
newId = newLines->InsertNextCell(npts,pts);
outCD->CopyData(inCD,cellId,newId);
for (i=0; i < npts; i++)
{
newId = newPts->InsertNextPoint(inputPoints->GetPoint(pts[i]));
outPD->CopyData(inPD,pts[i],newId);
}
continue; //skip the rest
}
// Start off by loading data into data structures
VertList.clear();
for (i=0; i < npts; i++)
{
if ( i == 0 || i == (npts-1) )
{
error = VTK_LARGE_FLOAT;
}
else
{
inputPoints->GetPoint(pts[i-1],x1);
inputPoints->GetPoint(pts[i],x);
inputPoints->GetPoint(pts[i+1],x2);
if ( (len = sqrt(vtkMath::Distance2BetweenPoints(x1,x2))) <= 0.0 )
{
error = 0.0;
}
else
{
dist = sqrt(vtkLine::DistanceToLine(x,x1,x2));
error = dist/len;
}
}
vertIterator = VertList.insert(VertList.end(),vtkPLineVertex(pts[i],error));
PQueue.push(vertIterator);
}//for all points in polyline
// Now process structures, deleting points until the decimation target is met.
vtkIdType currentNumPts = npts;
while ( 1.0 - (static_cast<double>(currentNumPts)/static_cast<double>(npts)) < this->TargetReduction &&
currentNumPts > 2 )
{
vertIterator = PQueue.top();
PQueue.pop();
currentNumPts--;
VertList.erase(vertIterator);
}
// What's left over is now spit out as a new polyline
newId = newLines->InsertNextCell(currentNumPts);
outCD->CopyData(inCD,cellId,newId);
for ( vertIterator = VertList.begin(); vertIterator != VertList.end(); ++ vertIterator )
{
newId = newPts->InsertNextPoint(inputPoints->GetPoint(vertIterator->PtId));
newLines->InsertCellPoint(newId);
outPD->CopyData(inPD,vertIterator->PtId,newId);
}
// Clean up in preparation for the next line
while ( !PQueue.empty() ) //empty the queue in preparation for the next line
{
PQueue.pop();
}
}//loop over all polylines
// Create output and clean up
output->SetPoints( newPts );
output->SetLines( newLines );
newLines->Delete();
newPts->Delete();
return 1;
}
//---------------------------------------------------------------------
void vtkDecimatePolylineFilter::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Target Reduction: " << this->TargetReduction << "\n";
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkDecimatePolylineFilter.h
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
// .NAME vtkDecimatePolylineFilter - reduce the number of lines in a polyline
// .SECTION Description
// vtkDecimatePolylineFilter is a filter to reduce the number of lines in a
// polyline. The algorithm functions by evaluating an error metric for each
// vertex (i.e., the distance of the vertex to a line defined from the two
// vertices on either side of the vertex). Then, these vertices are placed
// into a priority queue, and those with larger errors are deleted first.
// The decimation continues until the target reduction is reached.
//
// .SECTION Caveats
// This algorithm is a very simple implementation that overlooks some
// potential complexities. First, if a vertex is multiply connected,
// meaning that it is used by multiple polylines, then the extra
// topological constraints are ignored. Second, the error is not updated
// as vertices are deleted (similar to iteratively computing a quadric
// error metric). Thus, once calculated, the error is used to determine
// which vertices are removed. This can produce less than optimal results.
//
// .SECTION See Also
// vtkDecimate vtkDecimateProp vtkQuadricClustering vtkQuadricDecimation
#ifndef __vtkDecimatePolylineFilter_h
#define __vtkDecimatePolylineFilter_h
#include "vtkPolyDataAlgorithm.h"
class vtkPriorityQueue;
class VTK_GRAPHICS_EXPORT vtkDecimatePolylineFilter : public vtkPolyDataAlgorithm
{
public:
// Description:
// Standard methods for type information and printing.
vtkTypeRevisionMacro(vtkDecimatePolylineFilter,vtkPolyDataAlgorithm);
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Instantiate this object with a target reduction of 0.90.
static vtkDecimatePolylineFilter *New();
// Description:
// Specify the desired reduction in the total number of polygons (e.g., if
// TargetReduction is set to 0.9, this filter will try to reduce the data set
// to 10% of its original size).
vtkSetClampMacro(TargetReduction,double,0.0,1.0);
vtkGetMacro(TargetReduction,double);
protected:
vtkDecimatePolylineFilter();
~vtkDecimatePolylineFilter();
int RequestData(vtkInformation *, vtkInformationVector **, vtkInformationVector *);
double TargetReduction;
private:
vtkDecimatePolylineFilter(const vtkDecimatePolylineFilter&); // Not implemented.
void operator=(const vtkDecimatePolylineFilter&); // Not implemented.
};
#endif
/*=========================================================================
Program: Visualization Toolkit
Module: vtkFillHolesFilter.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
#include "vtkFillHolesFilter.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDoubleArray.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkTriangleStrip.h"
#include "vtkPolygon.h"
#include "vtkSphere.h"
#include "vtkMath.h"
vtkCxxRevisionMacro(vtkFillHolesFilter, "1.1");
vtkStandardNewMacro(vtkFillHolesFilter);
//------------------------------------------------------------------------
vtkFillHolesFilter::vtkFillHolesFilter()
{
this->HoleSize = 1.0;
}
//------------------------------------------------------------------------
vtkFillHolesFilter::~vtkFillHolesFilter()
{
}
//------------------------------------------------------------------------
int vtkFillHolesFilter::RequestData(
vtkInformation *vtkNotUsed(request),
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
// get the info objects
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
vtkInformation *outInfo = outputVector->GetInformationObject(0);
// get the input and ouptut
vtkPolyData *input = vtkPolyData::SafeDownCast(
inInfo->Get(vtkDataObject::DATA_OBJECT()));
vtkPolyData *output = vtkPolyData::SafeDownCast(
outInfo->Get(vtkDataObject::DATA_OBJECT()));
vtkPointData *pd=input->GetPointData(), *outPD=output->GetPointData();
vtkDebugMacro(<<"Executing hole fill operation");
// check the input, build data structures as necessary
vtkIdType numPts, npts, *pts;
vtkPoints *inPts=input->GetPoints();
vtkIdType numPolys = input->GetNumberOfPolys();
vtkIdType numStrips = input->GetNumberOfStrips();
if ( (numPts=input->GetNumberOfPoints()) < 1 || !inPts ||
(numPolys < 1 && numStrips < 1) )
{
vtkDebugMacro(<<"No input data!");
return 1;
}
vtkPolyData *Mesh = vtkPolyData::New();
Mesh->SetPoints(inPts);
vtkCellArray *newPolys, *inPolys=input->GetPolys(), *inStrips=input->GetStrips();
if ( numStrips > 0 )
{
newPolys = vtkCellArray::New();
if ( numPolys > 0 )
{
newPolys->DeepCopy(inPolys);
}
else
{
newPolys->Allocate(newPolys->EstimateSize(numStrips,5));
}
inStrips = input->GetStrips();
for ( inStrips->InitTraversal(); inStrips->GetNextCell(npts,pts); )
{
vtkTriangleStrip::DecomposeStrip(npts, pts, newPolys);
}
Mesh->SetPolys(newPolys);
newPolys->Delete();
}
else
{
newPolys = inPolys;
Mesh->SetPolys(newPolys);
}
Mesh->BuildLinks();
// Allocate storage for lines/points (arbitrary allocation sizes)
//
vtkPolyData *Lines = vtkPolyData::New();
vtkCellArray *newLines = vtkCellArray::New();
newLines->Allocate(numPts/10);
Lines->SetLines(newLines);
Lines->SetPoints(inPts);
// grab all free edges and place them into a temporary polydata
int abort=0;
vtkIdType cellId, p1, p2, numNei, i, newId, numCells=newPolys->GetNumberOfCells();
vtkIdType progressInterval=numCells/20+1;
vtkIdList *neighbors = vtkIdList::New();
neighbors->Allocate(VTK_CELL_SIZE);
for (cellId=0, newPolys->InitTraversal();
newPolys->GetNextCell(npts,pts) && !abort; cellId++)
{
if ( ! (cellId % progressInterval) ) //manage progress / early abort
{
this->UpdateProgress ((double)cellId / numCells);
abort = this->GetAbortExecute();
}
for (i=0; i < npts; i++)
{
p1 = pts[i];
p2 = pts[(i+1)%npts];
Mesh->GetCellEdgeNeighbors(cellId,p1,p2, neighbors);
numNei = neighbors->GetNumberOfIds();
if ( numNei < 1 )
{
newId = newLines->InsertNextCell(2);
newLines->InsertCellPoint(p1);
newLines->InsertCellPoint(p2);
}
}
}
// Track all free edges and see whether polygons can be built from them.
// For each polygon of appropriate HoleSize, triangulate the hole and
// add to the output list of cells
vtkIdType numHolesFilled=0;
numCells = newLines->GetNumberOfCells();
vtkCellArray *newCells = NULL;
if ( numCells > 3 ) //only do the work if there are free edges
{
double sphere[4];
vtkIdType startId, neiId, currentCellId, hints[2]; hints[0]=0; hints[1]=0;
vtkPolygon *polygon=vtkPolygon::New();
polygon->Points->SetDataTypeToDouble();
vtkIdList *endId = vtkIdList::New();
endId->SetNumberOfIds(1);
char *visited = new char [numCells];
memset(visited, 0, numCells);
Lines->BuildLinks(); //build the neighbor data structure
newCells = vtkCellArray::New();
newCells->DeepCopy(inPolys);
for (cellId=0; cellId < numCells && !abort; cellId++)
{
if ( ! visited[cellId] )
{
visited[cellId] = 1;
// Setup the polygon
Lines->GetCellPoints(cellId, npts, pts);
startId = pts[0];
polygon->PointIds->Reset();
polygon->Points->Reset();
polygon->PointIds->InsertId(0,pts[0]);
polygon->Points->InsertPoint(0,inPts->GetPoint(pts[0]));
// Work around the loop and terminate when the loop ends
endId->SetId(0,pts[1]);
int valid = 1;
currentCellId = cellId;
while ( startId != endId->GetId(0) && valid )
{
polygon->PointIds->InsertNextId(endId->GetId(0));
polygon->Points->InsertNextPoint(inPts->GetPoint(endId->GetId(0)));
Lines->GetCellNeighbors(currentCellId,endId,neighbors);
if ( neighbors->GetNumberOfIds() == 0 )
{
valid = 0;
}
else if ( neighbors->GetNumberOfIds() > 1 )
{
//have to logically split this vertex
valid = 0;
}
else
{
neiId = neighbors->GetId(0);
visited[neiId] = 1;
Lines->GetCellPoints(neiId,npts,pts);
endId->SetId( 0, (pts[0] != endId->GetId(0) ? pts[0] : pts[1] ) );
currentCellId = neiId;
}
}//while loop connected
// Evaluate the size of the loop and see if it is small enough
if ( valid )
{
vtkSphere::ComputeBoundingSphere(static_cast<vtkDoubleArray*>(polygon->Points->GetData())->GetPointer(0),
polygon->PointIds->GetNumberOfIds(),sphere,hints);
if ( sphere[3] <= this->HoleSize )
{
// Now triangulate the loop and pass to the output
numHolesFilled++;
polygon->NonDegenerateTriangulate(neighbors);
for ( i=0; i < neighbors->GetNumberOfIds(); i+=3 )
{
newCells->InsertNextCell(3);
newCells->InsertCellPoint(polygon->PointIds->GetId(neighbors->GetId(i)));
newCells->InsertCellPoint(polygon->PointIds->GetId(neighbors->GetId(i+1)));
newCells->InsertCellPoint(polygon->PointIds->GetId(neighbors->GetId(i+2)));
}
}//if hole small enough
}//if a valid loop
}//if not yet visited a line
}//for all lines
polygon->Delete();
endId->Delete();
delete [] visited;
}//if loops present in the input
// Clean up
neighbors->Delete();
Lines->Delete();
// No new points are created, so the points and point data can be passed
// through to the output.
output->SetPoints(inPts);
outPD->PassData(pd);
// New cells are created, so currently we do not pass the cell data.
// It would be pretty easy to extend the existing cell data and mark
// the new cells with special data values.
output->SetVerts(input->GetVerts());
output->SetLines(input->GetLines());
if ( newCells )
{
output->SetPolys(newCells);
newCells->Delete();
}
output->SetStrips(input->GetStrips());
Mesh->Delete();
newLines->Delete();
return 1;
}
//------------------------------------------------------------------------
void vtkFillHolesFilter::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os,indent);
os << indent << "Hole Size: " << this->HoleSize << "\n";
}
/*=========================================================================
Program: Visualization Toolkit
Module: vtkFillHolesFilter.h
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notice for more information.
=========================================================================*/
// .NAME vtkFillHolesFilter - identify and fill holes in meshes
// .SECTION Description
// vtkFillHolesFilter is a filter that identifies and fills holes in
// input vtkPolyData meshes. Holes are identified by locating
// boundary edges, linking them together into loops, and then
// triangulating the resulting loops. Note that you can specify
// an approximate limit to the size of the hole that can be filled.
//
// .SECTION Caveats
// Note that any mesh with boundary edges by definition has a
// topological hole. This even includes a reactangular grid
// (e.g., the output of vtkPlaneSource). In such situations, if
// the outer hole is filled, retriangulation of the hole will cause
// geometric overlap of the mesh. This can be prevented by using
// the hole size instance variable to prevent the larger holes
// from being triangulated.
//
// Note this filter only operates on polygons and triangle strips.
// Vertices and polylines are passed through untouched.
#ifndef __vtkFillHolesFilter_h
#define __vtkFillHolesFilter_h
#include "vtkPolyDataAlgorithm.h"
class vtkAbstractTransform;
class VTK_GRAPHICS_EXPORT vtkFillHolesFilter : public vtkPolyDataAlgorithm
{
public:
// Description:
// Standard methods for instantiation, type information and printing.
static vtkFillHolesFilter *New();
vtkTypeRevisionMacro(vtkFillHolesFilter,vtkPolyDataAlgorithm);
void PrintSelf(ostream& os, vtkIndent indent);
// Description:
// Specify the maximum hole size to fill. This is represented as a radius
// to the bounding circumsphere containing the hole. Note that this is an
// approximate area; the actual area cannot be computed without first