vtkFillHolesFilter.cxx 8.38 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
/*=========================================================================

  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"

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);

53
  // get the input and output
54 55 56 57 58 59 60 61
  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");
62

63 64 65 66 67
  // check the input, build data structures as necessary
  vtkIdType numPts, npts, *pts;
  vtkPoints *inPts=input->GetPoints();
  vtkIdType numPolys = input->GetNumberOfPolys();
  vtkIdType numStrips = input->GetNumberOfStrips();
68
  if ( (numPts=input->GetNumberOfPoints()) < 1 || !inPts ||
69
       (numPolys < 1 && numStrips < 1) )
70
  {
71 72
    vtkDebugMacro(<<"No input data!");
    return 1;
73
  }
74 75 76

  vtkPolyData *Mesh = vtkPolyData::New();
  Mesh->SetPoints(inPts);
77
  vtkCellArray *newPolys, *inPolys=input->GetPolys();
78
  if ( numStrips > 0 )
79
  {
80 81
    newPolys = vtkCellArray::New();
    if ( numPolys > 0 )
82
    {
83
      newPolys->DeepCopy(inPolys);
84
    }
85
    else
86
    {
87
      newPolys->Allocate(newPolys->EstimateSize(numStrips,5));
88
    }
89
    vtkCellArray *inStrips = input->GetStrips();
90
    for ( inStrips->InitTraversal(); inStrips->GetNextCell(npts,pts); )
91
    {
92
      vtkTriangleStrip::DecomposeStrip(npts, pts, newPolys);
93
    }
94 95
    Mesh->SetPolys(newPolys);
    newPolys->Delete();
96
  }
97
  else
98
  {
99 100
    newPolys = inPolys;
    Mesh->SetPolys(newPolys);
101
  }
102 103 104 105 106 107 108 109 110
  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);
111

112 113
  // grab all free edges and place them into a temporary polydata
  int abort=0;
114
  vtkIdType cellId, p1, p2, numNei, i, numCells=newPolys->GetNumberOfCells();
115 116 117
  vtkIdType progressInterval=numCells/20+1;
  vtkIdList *neighbors = vtkIdList::New();
  neighbors->Allocate(VTK_CELL_SIZE);
118
  for (cellId=0, newPolys->InitTraversal();
119
       newPolys->GetNextCell(npts,pts) && !abort; cellId++)
120
  {
121
    if ( ! (cellId % progressInterval) ) //manage progress / early abort
122
    {
123
      this->UpdateProgress (static_cast<double>(cellId) / numCells);
124
      abort = this->GetAbortExecute();
125
    }
126

127
    for (i=0; i < npts; i++)
128
    {
129 130 131 132 133 134 135
      p1 = pts[i];
      p2 = pts[(i+1)%npts];

      Mesh->GetCellEdgeNeighbors(cellId,p1,p2, neighbors);
      numNei = neighbors->GetNumberOfIds();

      if ( numNei < 1 )
136
      {
137
        newLines->InsertNextCell(2);
138 139 140 141
        newLines->InsertCellPoint(p1);
        newLines->InsertCellPoint(p2);
      }
    }
142
  }
143

144 145 146 147 148
  // 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();
149
  vtkCellArray *newCells = nullptr;
150
  if ( numCells >= 3 ) //only do the work if there are free edges
151
  {
152 153 154 155 156 157 158 159 160 161 162 163 164
    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++)
165
    {
166
      if ( ! visited[cellId] )
167
      {
168 169 170 171 172 173 174 175 176 177 178 179 180 181
        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 )
182
        {
183 184 185 186
          polygon->PointIds->InsertNextId(endId->GetId(0));
          polygon->Points->InsertNextPoint(inPts->GetPoint(endId->GetId(0)));
          Lines->GetCellNeighbors(currentCellId,endId,neighbors);
          if ( neighbors->GetNumberOfIds() == 0 )
187
          {
188
            valid = 0;
189
          }
190
          else if ( neighbors->GetNumberOfIds() > 1 )
191
          {
192 193
            //have to logically split this vertex
            valid = 0;
194
          }
195
          else
196
          {
197 198 199 200 201
            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;
202 203
          }
        }//while loop connected
204 205 206

        // Evaluate the size of the loop and see if it is small enough
        if ( valid )
207
        {
208 209 210
          vtkSphere::ComputeBoundingSphere(static_cast<vtkDoubleArray*>(polygon->Points->GetData())->GetPointer(0),
                                           polygon->PointIds->GetNumberOfIds(),sphere,hints);
          if ( sphere[3] <= this->HoleSize )
211
          {
212 213 214 215
            // Now triangulate the loop and pass to the output
            numHolesFilled++;
            polygon->NonDegenerateTriangulate(neighbors);
            for ( i=0; i < neighbors->GetNumberOfIds(); i+=3 )
216
            {
217 218 219 220
              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)));
221 222 223 224 225
            }
          }//if hole small enough
        }//if a valid loop
      }//if not yet visited a line
    }//for all lines
226 227 228
    polygon->Delete();
    endId->Delete();
    delete [] visited;
229
  }//if loops present in the input
230

231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
  // 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 )
246
  {
247 248
    output->SetPolys(newCells);
    newCells->Delete();
249
  }
250
  else
251
  {
252
    output->SetPolys(inPolys);
253
  }
254 255 256 257 258 259 260 261 262 263 264 265 266 267
  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";
}