[vtkusers] Subtracting polydata

[Johnny Norris]

Hi,

I'm trying to subtract vtkPolyData B from vtkPolyData A to create
vtkPolyData C; that is, C will contain all cells in A that are not also in B.
To say that I'm having a little trouble with this would be an understatement.

Right now I have a routine that, for each cell X in B, it searches for points
in A within a certain tolerance of the points in X.  Then it calls
A->GetCellNeighbors() to find all the cells in A that use these points.
Finally, for each cell Y returned by GetCellNeighbors(), it makes sure that
X and Y are the same cell type, have the same number of points, and that the
bounds of Y fall completely within the bounds of B.  If all of these tests
pass, then I call A->DeleteCell() to delete this cell.

Anyway, it isn't working well.  Sometimes polygons are subtracted, and
sometimes they aren't.  I'm hoping that someone out there has already done
something like this, and is willing to share their code with me.

Here's my SubtractPolyData() routine, in case anyone can see something
obviously wrong with it.  ds is the polydata I subtracting from (A), locator
is a point locator initialized with ds (using SetDataSet() and
BuildLocator()), and del is the polydata I'm subtracting out (B).

Any help would be appreciated,
John
-- 
John Norris
Research Programmer
Center for Simulation of Advanced Rockets
http://www.uiuc.edu/ph/www/jnorris
-------------- next part --------------
static void SubtractPolyData(vtkPolyData *ds, vtkPointLocator *locator,
                             vtkPolyData *del)
{
   vtkPoints *pts = ds->GetPoints();
   vtkPoints *dPts = del->GetPoints();
   vtkCell *cell;
   vtkIdList *ids = vtkIdList::New();
   vtkIdList *ptIds = vtkIdList::New();
   vtkIdList *dPtIds = vtkIdList::New();
   vtkIdList *cellIds = vtkIdList::New();
   vtkIdType dCellId, numDCells = del->GetNumberOfCells(), dPtId, numDPts;
   vtkIdType cellId, ptId, numPts, type, dType;
   float dP[3], radius = locator->GetTolerance();
   float dBounds[6], cBounds[6];
   int i, j, found;

   // Stretch the bounds by the given tolerance.
   del->GetBounds(dBounds);
   dBounds[0] -= radius;
   dBounds[1] += radius;
   dBounds[2] -= radius;
   dBounds[3] += radius;
   dBounds[4] -= radius;
   dBounds[5] += radius;

   for (dCellId=0; dCellId<numDCells; ++dCellId) {
      del->GetCellPoints(dCellId, dPtIds);
      numDPts = dPtIds->GetNumberOfIds();
      ptIds->Initialize();
      ptIds->SetNumberOfIds(numDPts);
      found = 1;
      for (i=0; i<numDPts; ++i) {
         dPtId = dPtIds->GetId(i);
         dPts->GetPoint(dPtId, dP);
         locator->FindPointsWithinRadius(radius, dP, ids);
         numPts = ids->GetNumberOfIds();
         if (numPts == 0) {
            found = 0;
            break;
         }
         locator->FindClosestNPoints(numPts, dP, ids);
         for (j=0; j<numPts; ++j) {
            ptId = ids->GetId(j);
            if (ptIds->IsId(ptId) == -1) {
               ptIds->SetId(i, ptId);
               break;
            }
         }
         if (j == numPts) {
            found = 0;
            break;
         }
      }

      if (found) {
         dType = del->GetCellType(dCellId);
         ds->GetCellNeighbors(-1, ptIds, cellIds);
         for (i=0; i<cellIds->GetNumberOfIds(); ++i) {
            cellId = cellIds->GetId(i);
            cell = ds->GetCell(cellId);
            cell->GetBounds(cBounds);
            if (dType == ds->GetCellType(cellId)
                && cell->GetNumberOfPoints() == numDPts
                && IsContainedIn(cBounds, dBounds)) {
               continue;
            }
            ds->DeleteCell(cellId);
         }
      }
   }

   cellIds->Delete();
   dPtIds->Delete();
   ptIds->Delete();
   ids->Delete();
}

inline int IsContainedIn(float *little, float *big)
{
   int i;
   for (i=0; i<6; i+=2)
      if (little[i] < big[i] || little[i+1] > big[i+1])
         return 0;

   return 1;
}