vtkBandedPolyDataContourFilter.cxx 26.7 KB
Newer Older
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkBandedPolyDataContourFilter.cxx

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9

10 11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
13 14 15

=========================================================================*/
#include "vtkBandedPolyDataContourFilter.h"
16

17 18
#include "vtkCellArray.h"
#include "vtkCellData.h"
19
#include "vtkEdgeTable.h"
20
#include "vtkExecutive.h"
21
#include "vtkFloatArray.h"
22 23
#include "vtkInformation.h"
#include "vtkInformationVector.h"
24
#include "vtkObjectFactory.h"
25
#include "vtkPointData.h"
26 27
#include "vtkPolyData.h"
#include "vtkTriangleStrip.h"
28
#include "vtkDoubleArray.h"
29

30
#include <float.h>
31

Brad King's avatar
Brad King committed
32
vtkStandardNewMacro(vtkBandedPolyDataContourFilter);
33

34
//------------------------------------------------------------------------------
35 36 37
// Construct object.
vtkBandedPolyDataContourFilter::vtkBandedPolyDataContourFilter()
{
38
  this->ContourValues = vtkContourValues::New();
Will Schroeder's avatar
Will Schroeder committed
39
  this->Clipping = 0;
Will Schroeder's avatar
Will Schroeder committed
40
  this->ScalarMode = VTK_SCALAR_MODE_INDEX;
41 42 43 44 45 46

  this->SetNumberOfOutputPorts(2);

  vtkPolyData *output2 = vtkPolyData::New();
  this->GetExecutive()->SetOutputData(1, output2);
  output2->Delete();
47
  this->ClipTolerance = FLT_EPSILON;
48
  this->InternalClipTolerance = FLT_EPSILON;
49
  this->GenerateContourEdges = 0;
50 51
}

52
//------------------------------------------------------------------------------
53 54
vtkBandedPolyDataContourFilter::~vtkBandedPolyDataContourFilter()
{
55 56 57
  this->ContourValues->Delete();
}

58
//------------------------------------------------------------------------------
Ken Martin's avatar
Ken Martin committed
59
int vtkBandedPolyDataContourFilter::ComputeScalarIndex(double val)
60
{
61

62 63 64 65 66 67 68
  for (int i=0; i < (this->NumberOfClipValues-1); i++)
    {
    if ( val >= this->ClipValues[i] && val < this->ClipValues[i+1]  )
      {
      return i;
      }
    }
69
  return this->NumberOfClipValues - 1;
70

71 72
}

73
//------------------------------------------------------------------------------
Ken Martin's avatar
Ken Martin committed
74
int vtkBandedPolyDataContourFilter::IsContourValue(double val)
Will Schroeder's avatar
Will Schroeder committed
75 76 77 78 79 80 81 82 83 84 85 86 87 88
{
  int i;

  // Check to see whether a vertex is an intersection point.
  for ( i=0; i < this->NumberOfClipValues; i++)
    {
    if ( val == this->ClipValues[i] )
      {
      return 1;
      }
    }
  return 0;
}

89
//------------------------------------------------------------------------------
Will Schroeder's avatar
Will Schroeder committed
90 91
// Return a flag that indicates that the ordering of vertices along the
// edge is not from v1->v2, where v1 < v2.
92
int vtkBandedPolyDataContourFilter::ClipEdge(int v1, int v2,
93
                                             vtkPoints *newPts,
94 95
                                             vtkDataArray *inScalars,
                                             vtkDoubleArray *outScalars,
96 97 98
                                             vtkPointData *inPD,
                                             vtkPointData *outPD)
{
99
  double x[3], t, sNew;
Ken Martin's avatar
Ken Martin committed
100
  double x1[3], x2[3];
101
  int ptId;
Will Schroeder's avatar
Will Schroeder committed
102
  int reverse = (v1 < v2 ? 0 : 1);
103

104 105
  newPts->GetPoint(v1, x1);
  newPts->GetPoint(v2, x2);
106

107 108
  double s1 = inScalars->GetTuple1(v1);
  double s2 = inScalars->GetTuple1(v2);
109

110
  if ( s1 <= s2 )
111
    {
112 113
    int idx1 = this->ComputeScalarIndex(s1);
    int idx2 = this->ComputeScalarIndex(s2);
114

115 116 117 118 119 120 121 122
    for (int i=1; i < (idx2-idx1+1); i++)
      {
      t = (this->ClipValues[idx1+i] - s1) / (s2 - s1);
      x[0] = x1[0] + t*(x2[0]-x1[0]);
      x[1] = x1[1] + t*(x2[1]-x1[1]);
      x[2] = x1[2] + t*(x2[2]-x1[2]);
      ptId = newPts->InsertNextPoint(x);
      outPD->InterpolateEdge(inPD,ptId,v1,v2,t);
123 124
      // We cannot use directly s1 + t*(s2-s1) as is causes rounding error
      sNew = this->ClipValues[idx1+i];
125
      outScalars->InsertTuple1(ptId,sNew);
126
      }
Will Schroeder's avatar
Will Schroeder committed
127
    return reverse;
128
    }
129 130
  else
    {
131 132
    int idx1 = this->ComputeScalarIndex(s1);
    int idx2 = this->ComputeScalarIndex(s2);
133

134 135 136 137 138 139 140 141
    for (int i=1; i < (idx1-idx2+1); i++)
      {
      t = (this->ClipValues[idx2+i] - s1) / (s2 - s1);
      x[0] = x1[0] + t*(x2[0]-x1[0]);
      x[1] = x1[1] + t*(x2[1]-x1[1]);
      x[2] = x1[2] + t*(x2[2]-x1[2]);
      ptId = newPts->InsertNextPoint(x);
      outPD->InterpolateEdge(inPD,ptId,v1,v2,t);
142 143
      // We cannot use directly s1 + t*(s2-s1) as is causes rounding error
      sNew = this->ClipValues[idx2+i];
144
      outScalars->InsertTuple1(ptId,sNew);
145
      }
Will Schroeder's avatar
Will Schroeder committed
146
    return ((reverse+1) % 2);
147
    }
148 149 150
}


151
//------------------------------------------------------------------------------
152
extern "C" {
Sean McBride's avatar
Sean McBride committed
153
static int vtkCompareClipValues(const void *val1, const void *val2)
154
{
155
  if ( *((double*)val1) < *((double*)val2) )
156 157 158
    {
    return (-1);
    }
159
  else if ( *((double*)val1) > *((double*)val2) )
160 161 162
    {
    return (1);
    }
163
  else
164 165 166
    {
    return (0);
    }
167
}
168 169
}

170
//------------------------------------------------------------------------------
Will Schroeder's avatar
Will Schroeder committed
171 172
inline int vtkBandedPolyDataContourFilter::InsertCell(vtkCellArray *cells,
                                                      int npts, vtkIdType *pts,
Ken Martin's avatar
Ken Martin committed
173
                                                      int cellId, double s,
Will Schroeder's avatar
Will Schroeder committed
174 175
                                                      vtkFloatArray *newS)
{
176
  int idx = this->ComputeScalarIndex(s+this->InternalClipTolerance);
Will Schroeder's avatar
Will Schroeder committed
177

178
  if ( !this->Clipping ||
179
       (idx >= this->ClipIndex[0] && idx < this->ClipIndex[1]) )
Will Schroeder's avatar
Will Schroeder committed
180 181 182 183 184 185 186 187 188
    {
    cells->InsertNextCell(npts,pts);

    if ( this->ScalarMode == VTK_SCALAR_MODE_INDEX )
      {
      newS->InsertTuple1(cellId++,idx);
      }
    else
      {
189
      newS->InsertTuple1(cellId++,this->ClipValues[idx]);
Will Schroeder's avatar
Will Schroeder committed
190 191 192 193 194 195
      }
    }
  return cellId;
}


196
//------------------------------------------------------------------------------
197
// Create filled contours for polydata
198 199 200 201
int vtkBandedPolyDataContourFilter::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
202
{
203 204 205 206
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

207
  // get the input and output
208 209 210 211 212
  vtkPolyData *input = vtkPolyData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

213
  vtkPointData *pd = input->GetPointData();
214
  vtkPointData *outPD = output->GetPointData();
215
  vtkCellData *outCD = output->GetCellData();
216
  vtkPoints *inPts = input->GetPoints();
217
  vtkDataArray *inScalars = pd->GetScalars();
Will Schroeder's avatar
Will Schroeder committed
218
  int abort=0;
219
  vtkPoints *newPts;
220 221 222 223
  int i, j, idx=0;
  vtkIdType npts = 0;
  vtkIdType cellId=0;
  vtkIdType *pts = 0;
Will Schroeder's avatar
Will Schroeder committed
224 225
  int numEdgePts, numNewPts, maxCellSize;
  vtkIdType v, vR, *intPts;
226 227 228
  int intsIdx, reverse;
  vtkIdType intLoc;
  vtkIdType numIntPts, intsInc;
229
  vtkIdType numPts, numCells, estimatedSize;
230

Will Schroeder's avatar
Will Schroeder committed
231
  vtkDebugMacro(<<"Executing banded contour filter");
232 233 234

  //  Check input
  //
235

236
  numCells = input->GetNumberOfCells();
237
  if ( !inPts || (numPts=inPts->GetNumberOfPoints()) < 1 ||
238
       !inScalars || numCells < 1 )
239 240
    {
    vtkErrorMacro(<<"No input data!");
241
    return 1;
242 243
    }

244 245 246 247
  // Set up supplemental data structures for processing edge/generating
  // intersections. First we sort the contour values into an ascending
  // list of clip values including the extreme min/max values.
  this->NumberOfClipValues = this->ContourValues->GetNumberOfContours() + 2;
Ken Martin's avatar
Ken Martin committed
248
  this->ClipValues = new double[this->NumberOfClipValues];
Ken Martin's avatar
Ken Martin committed
249
  double range[2];
250 251 252
  inScalars->GetRange(range);

  // base clip tolerance on overall input scalar range
253
  this->InternalClipTolerance = this->ClipTolerance*(range[1] - range[0]);
254

255
  this->ClipValues[0] =
256 257 258
    (range[0]<this->ContourValues->GetValue(0))?
    (range[0]):
    (this->ContourValues->GetValue(0));
259 260

  this->ClipValues[this->NumberOfClipValues - 1] =
261
    (range[1]>this->ContourValues->GetValue(this->NumberOfClipValues-2-1))?
262
    (range[1]):
263
    (this->ContourValues->GetValue(this->NumberOfClipValues-2-1));
264

265
  for ( i=1; i<this->NumberOfClipValues-1; i++)
266
    {
267
    this->ClipValues[i] = this->ContourValues->GetValue(i-1);
268
    }
Will Schroeder's avatar
Will Schroeder committed
269

270
  qsort(this->ClipValues, this->NumberOfClipValues, sizeof(double),
Will Schroeder's avatar
Will Schroeder committed
271 272
        vtkCompareClipValues);

273 274 275
  // toss out values which are too close together, currently within FLT_EPSILON%
  // of each other based on full scalar range, but could define temporary based
  // on percentage of scalar range...
276 277
  for ( i=0; i<(this->NumberOfClipValues-1); i++)
    {
278
    if ( (this->ClipValues[i] + this->InternalClipTolerance) >= this->ClipValues[i+1] )
279 280 281
      {
      for (j=i+1; j<(this->NumberOfClipValues-2); j++)
        {
282
        this->ClipValues[j] = this->ClipValues[j+1];
283 284 285 286
        }
      this->NumberOfClipValues--;
      }
    }
287

288
  this->ClipIndex[0] =
289 290
    this->ComputeScalarIndex(this->ContourValues->GetValue(0));
  this->ClipIndex[1] = this->ComputeScalarIndex(
Will Schroeder's avatar
Will Schroeder committed
291 292
    this->ContourValues->GetValue(this->ContourValues->GetNumberOfContours()-1));

293 294 295
  //
  // Estimate allocation size, stolen from vtkContourGrid...
  //
296
  estimatedSize=static_cast<vtkIdType>(pow(static_cast<double>(numCells),.9));
297 298 299 300 301 302 303
  estimatedSize *= this->NumberOfClipValues;
  estimatedSize = estimatedSize / 1024 * 1024; // multiple of 1024
  if (estimatedSize < 1024)
    {
    estimatedSize = 1024;
    }

304
  // The original set of points and point data are copied. Later on
305 306
  // intersection points due to clipping will be created.
  newPts = vtkPoints::New();
307

308 309 310
  // Note: since we use the output scalars in the execution of the algorithm,
  // the output point scalars MUST BE double or bad things happen due to
  // numerical precision issues.
311
  newPts->Allocate(estimatedSize,estimatedSize);
312
  outPD->CopyScalarsOff();
313
  outPD->InterpolateAllocate(pd,3*numPts,numPts);
314 315 316
  vtkDoubleArray *outScalars = vtkDoubleArray::New();
  outScalars->Allocate(3*numPts,numPts);
  outPD->SetScalars(outScalars);
Berk Geveci's avatar
Berk Geveci committed
317
  outScalars->Delete();
318

319
  for (i=0; i<numPts; i++)
320
    {
321 322
    newPts->InsertPoint(i,inPts->GetPoint(i));
    outPD->CopyData(pd, i, i);
323
    outScalars->InsertTuple1(i, inScalars->GetTuple1(i));
324 325
    }

326 327 328
  // These are the new cell scalars
  vtkFloatArray *newScalars = vtkFloatArray::New();
  newScalars->Allocate(numCells*5,numCells);
329
  newScalars->SetName("Scalars");
330

Will Schroeder's avatar
Will Schroeder committed
331 332 333 334
  // Used to keep track of intersections
  vtkEdgeTable *edgeTable = vtkEdgeTable::New();
  vtkCellArray *intList = vtkCellArray::New(); //intersection point ids

335 336
  // All vertices are filled and passed through; poly-vertices are broken
  // into single vertices. Cell data per vertex is set.
337
  //
338
  if ( input->GetVerts()->GetNumberOfCells() > 0 )
339
    {
340 341 342
    vtkCellArray *verts = input->GetVerts();
    vtkCellArray *newVerts = vtkCellArray::New();
    newVerts->Allocate(verts->GetSize());
343
    for ( verts->InitTraversal(); verts->GetNextCell(npts,pts) && !abort;
344
          abort=this->GetAbortExecute() )
345
      {
346
      for (i=0; i<npts; i++)
347
        {
348
        newVerts->InsertNextCell(1,pts+i);
349
        idx = this->ComputeScalarIndex(inScalars->GetTuple1(pts[i]));
350
        newScalars->InsertTuple1(cellId++,idx);
351 352
        }
      }
353 354 355
    output->SetVerts(newVerts);
    newVerts->Delete();
    }
356
  this->UpdateProgress(0.05);
357

358 359 360 361 362
  // Lines are chopped into line segments.
  //
  if ( input->GetLines()->GetNumberOfCells() > 0 )
    {
    vtkCellArray *lines = input->GetLines();
Will Schroeder's avatar
Will Schroeder committed
363 364 365 366 367

    maxCellSize = lines->GetMaxCellSize();
    maxCellSize *= (1 + this->NumberOfClipValues);

    vtkIdType *fullLine = new vtkIdType [maxCellSize];
368 369
    vtkCellArray *newLines = vtkCellArray::New();
    newLines->Allocate(lines->GetSize());
Will Schroeder's avatar
Will Schroeder committed
370 371 372
    edgeTable->InitEdgeInsertion(numPts,1); //store attributes on edge

    //start by generating intersection points
373
    for ( lines->InitTraversal(); lines->GetNextCell(npts,pts) && !abort;
374
          abort=this->GetAbortExecute() )
375
      {
376
      for (i=0; i<(npts-1); i++)
377
        {
Will Schroeder's avatar
Will Schroeder committed
378
        numNewPts = newPts->GetNumberOfPoints();
379 380
        reverse = this->ClipEdge(pts[i],pts[i+1],newPts,inScalars,outScalars,
                                 pd,outPD);
Will Schroeder's avatar
Will Schroeder committed
381 382
        numEdgePts = newPts->GetNumberOfPoints() - numNewPts;
        if ( numEdgePts > 0 ) //there is an intersection
383
          {
Will Schroeder's avatar
Will Schroeder committed
384
          if ( !reverse )
385 386
            {
            for (j=0; j<numEdgePts; j++) fullLine[j] = numNewPts + j;
Will Schroeder's avatar
Will Schroeder committed
387 388
            }
          else
389 390 391
            {
            for (j=0; j<numEdgePts; j++)
              fullLine[j] = numNewPts + numEdgePts - j - 1;
Will Schroeder's avatar
Will Schroeder committed
392 393 394 395
            }
          intList->InsertNextCell(numEdgePts,fullLine);
          edgeTable->InsertEdge(pts[i],pts[i+1], //associate ints with edge
                                intList->GetInsertLocation(numEdgePts));
396
          }
Will Schroeder's avatar
Will Schroeder committed
397
        else //no intersection points along the edge
398
          {
Will Schroeder's avatar
Will Schroeder committed
399
          edgeTable->InsertEdge(pts[i],pts[i+1],-1); //-1 means no points
400
          }
Will Schroeder's avatar
Will Schroeder committed
401 402 403 404
        }//for all line segments in this line
      }

    //now create line segments
405
    for ( lines->InitTraversal(); lines->GetNextCell(npts,pts) && !abort;
406
          abort=this->GetAbortExecute() )
Will Schroeder's avatar
Will Schroeder committed
407 408 409 410 411 412 413
      {
      for (i=0; i<(npts-1); i++)
        {
        v = pts[i];
        vR = pts[i+1];

        newLines->InsertNextCell(2);
414

415
        newScalars->InsertTuple1(cellId++,
416
                    this->ComputeScalarIndex(outScalars->GetTuple1(v)));
Will Schroeder's avatar
Will Schroeder committed
417 418 419
        newLines->InsertCellPoint(v);

        if ( (intLoc=edgeTable->IsEdge(v,vR)) != -1 )
420
          {
Will Schroeder's avatar
Will Schroeder committed
421 422 423 424 425 426 427 428
          intList->GetCell(intLoc,numIntPts,intPts);
          if ( v < vR ) {intsIdx = 0; intsInc=1;} //order of the edge
          else {intsIdx=numIntPts-1; intsInc=(-1);}

          for ( ; intsIdx >= 0 && intsIdx < numIntPts; intsIdx += intsInc )
            {
            newLines->InsertCellPoint(intPts[intsIdx]);
            newLines->InsertNextCell(2);
429 430

            newScalars->InsertTuple1(cellId++, this->ComputeScalarIndex(
Will Schroeder's avatar
Will Schroeder committed
431 432 433
              outScalars->GetTuple1(intPts[intsIdx])));
            newLines->InsertCellPoint(intPts[intsIdx]);
            }
434
          }
Will Schroeder's avatar
Will Schroeder committed
435
        newLines->InsertCellPoint(vR);
436
        }
437
      }
Will Schroeder's avatar
Will Schroeder committed
438 439 440

    delete [] fullLine;

441 442
    output->SetLines(newLines);
    newLines->Delete();
443
    }
444
  this->UpdateProgress(0.1);
445

446
  // Polygons are assumed convex and chopped into filled, convex polygons.
447
  // Triangle strips are treated similarly.
448
  //
Sean McBride's avatar
Sean McBride committed
449 450
  vtkIdType numPolys = input->GetPolys()->GetNumberOfCells();
  vtkIdType numStrips = input->GetStrips()->GetNumberOfCells();
Will Schroeder's avatar
Will Schroeder committed
451
  if ( numPolys > 0 || numStrips > 0 )
452 453 454 455 456 457
    {
    // Set up processing. We are going to store an ordered list of
    // intersections along each edge (ordered from smallest point id
    // to largest). These will later be connected into convex polygons
    // which represent a filled region in the cell.
    //
Will Schroeder's avatar
Will Schroeder committed
458
    edgeTable->InitEdgeInsertion(numPts,1); //store attributes on edge
Will Schroeder's avatar
Will Schroeder committed
459
    intList->Reset();
Will Schroeder's avatar
Will Schroeder committed
460 461

    vtkCellArray *polys = input->GetPolys();
Will Schroeder's avatar
Will Schroeder committed
462
    vtkCellArray *tmpPolys = NULL;
Will Schroeder's avatar
Will Schroeder committed
463

464
    // If contour edges requested, set things up.
465
    vtkCellArray *contourEdges=0;
466 467 468 469 470 471 472 473 474
    if ( this->GenerateContourEdges )
      {
      contourEdges = vtkCellArray::New();
      contourEdges->Allocate(numCells);
      this->GetContourEdgesOutput()->SetLines(contourEdges);
      contourEdges->Delete();
      this->GetContourEdgesOutput()->SetPoints(newPts);
      }

Will Schroeder's avatar
Will Schroeder committed
475 476
    // Set up structures for processing polygons
    maxCellSize = polys->GetMaxCellSize();
477 478 479 480
    if( maxCellSize == 0 )
      {
      maxCellSize = input->GetStrips()->GetMaxCellSize();
      }
Will Schroeder's avatar
Will Schroeder committed
481 482 483
    maxCellSize *= (1 + this->NumberOfClipValues);

    vtkIdType *newPolygon = new vtkIdType [maxCellSize];
Ken Martin's avatar
Ken Martin committed
484
    double *s = new double [maxCellSize]; //scalars at vertices
Will Schroeder's avatar
Will Schroeder committed
485 486 487 488 489 490
    int *isContourValue = new int [maxCellSize];
    int *isOriginalVertex = new int [maxCellSize];
    vtkIdType *fullPoly = new vtkIdType [maxCellSize];

    // Lump strips and polygons together.
    // Decompose strips into triangles.
491
    if ( numStrips > 0 )
Will Schroeder's avatar
Will Schroeder committed
492 493 494 495 496 497 498
      {
      vtkCellArray *strips = input->GetStrips();
      tmpPolys = vtkCellArray::New();
      if ( numPolys > 0 )
        {
        tmpPolys->DeepCopy(polys);
        }
499
      else
Will Schroeder's avatar
Will Schroeder committed
500 501 502 503 504 505 506 507 508
        {
        tmpPolys->Allocate(polys->EstimateSize(numStrips,5));
        }
      for ( strips->InitTraversal(); strips->GetNextCell(npts,pts); )
        {
        vtkTriangleStrip::DecomposeStrip(npts, pts, tmpPolys);
        }
      polys = tmpPolys;
      }
509

Will Schroeder's avatar
Will Schroeder committed
510 511
    // Process polygons to produce edge intersections.------------------------
    //
512 513 514
    numPolys = polys->GetNumberOfCells();
    vtkIdType updateCount = numPolys/20 + 1;
    vtkIdType count=0;
515
    for ( polys->InitTraversal(); polys->GetNextCell(npts,pts) && !abort;
516
      abort=this->GetAbortExecute() )
Will Schroeder's avatar
Will Schroeder committed
517
      {
518 519
      if  ( ! (++count % updateCount) )
        {
520
        this->UpdateProgress(0.1 + 0.45*(static_cast<double>(count)/numPolys));
521 522
        }

Will Schroeder's avatar
Will Schroeder committed
523 524 525 526 527 528 529
      for (i=0; i<npts; i++)
        {
        v = pts[i];
        vR = pts[(i+1) % npts];
        if ( edgeTable->IsEdge(v,vR) == -1 )
          {
          numNewPts = newPts->GetNumberOfPoints();
530
          reverse = this->ClipEdge(v,vR,newPts,inScalars,outScalars,pd,outPD);
Will Schroeder's avatar
Will Schroeder committed
531 532 533
          numEdgePts = newPts->GetNumberOfPoints() - numNewPts;
          if ( numEdgePts > 0 )
            {
Will Schroeder's avatar
Will Schroeder committed
534
            if ( !reverse )
535 536
              {
              for (j=0; j<numEdgePts; j++) fullPoly[j] = numNewPts + j;
Will Schroeder's avatar
Will Schroeder committed
537 538
              }
            else
539 540 541
              {
              for (j=0; j<numEdgePts; j++)
                fullPoly[j] = numNewPts + numEdgePts - j - 1;
Will Schroeder's avatar
Will Schroeder committed
542
              }
Will Schroeder's avatar
Will Schroeder committed
543 544 545
            intList->InsertNextCell(numEdgePts,fullPoly);
            edgeTable->InsertEdge(v,vR, //associate ints with edge
                                  intList->GetInsertLocation(numEdgePts));
Will Schroeder's avatar
Will Schroeder committed
546 547 548 549 550 551 552 553
            }
          else //no intersection points along the edge
            {
            edgeTable->InsertEdge(v,vR,-1); //-1 means no points
            }
          }//if edge not processed yet
        }
      }//for all polygons
554

555
    // Process polygons to produce output triangles------------------------
Will Schroeder's avatar
Will Schroeder committed
556
    //
Will Schroeder's avatar
Will Schroeder committed
557 558
    vtkCellArray *newPolys = vtkCellArray::New();
    newPolys->Allocate(polys->GetSize());
Will Schroeder's avatar
Will Schroeder committed
559
    int intersectionPoint;
560 561
    int mL, mR, m2L, m2R;
    int numPointsToAdd, numLeftPointsToAdd, numRightPointsToAdd;
Will Schroeder's avatar
Will Schroeder committed
562
    int numPolyPoints, numFullPts;
563
    count = 0;
564
    for ( polys->InitTraversal(); polys->GetNextCell(npts,pts) && !abort;
565
          abort=this->GetAbortExecute() )
Will Schroeder's avatar
Will Schroeder committed
566
      {
567 568
      if  ( ! (++count % updateCount) )
        {
569 570
        this->UpdateProgress(0.55 +
                             0.45*(static_cast<double>(count)/numPolys));
571 572
        }

Will Schroeder's avatar
Will Schroeder committed
573
      //Create a new polygon that includes all the points including the
Will Schroeder's avatar
Will Schroeder committed
574
      //intersection vertices. This hugely simplifies the logic of the
Will Schroeder's avatar
Will Schroeder committed
575 576
      //code.
      for ( intersectionPoint=0, numFullPts=0, i=0; i<npts; i++)
Will Schroeder's avatar
Will Schroeder committed
577
        {
Will Schroeder's avatar
Will Schroeder committed
578 579
        v = pts[i];
        vR = pts[(i+1)%npts];
580

Will Schroeder's avatar
Will Schroeder committed
581
        s[numFullPts] = outScalars->GetTuple1(v);
Will Schroeder's avatar
Will Schroeder committed
582
        isContourValue[numFullPts] = this->IsContourValue(s[numFullPts]);
Will Schroeder's avatar
Will Schroeder committed
583
        isOriginalVertex[numFullPts] = 1;
Will Schroeder's avatar
Will Schroeder committed
584
        fullPoly[numFullPts++] = v;
Will Schroeder's avatar
Will Schroeder committed
585

Will Schroeder's avatar
Will Schroeder committed
586
        //see whether intersection points need to be added.
Will Schroeder's avatar
Will Schroeder committed
587
        if ( (intLoc=edgeTable->IsEdge(v,vR)) != -1 )
Will Schroeder's avatar
Will Schroeder committed
588 589 590
          {
          intersectionPoint = 1;
          intList->GetCell(intLoc,numIntPts,intPts);
Will Schroeder's avatar
Will Schroeder committed
591
          if ( v < vR ) {intsIdx = 0; intsInc=1;} //order of the edge
Will Schroeder's avatar
Will Schroeder committed
592
          else {intsIdx=numIntPts-1; intsInc=(-1);}
Will Schroeder's avatar
Will Schroeder committed
593 594 595 596 597 598 599 600 601
          for ( ; intsIdx >= 0 && intsIdx < numIntPts; intsIdx += intsInc )
            {
            s[numFullPts] = outScalars->GetTuple1(intPts[intsIdx]);
            isContourValue[numFullPts] = 1;
            isOriginalVertex[numFullPts] = 0;
            fullPoly[numFullPts++] = intPts[intsIdx];
            }
          }
        } //for all points and edges
602

Will Schroeder's avatar
Will Schroeder committed
603 604
      //Very important: have to find the right starting vertex. The vertex
      //needs to be one where the contour values increase in both directions.
Will Schroeder's avatar
Will Schroeder committed
605
      //Really should check whether the vertex is convex.
Ken Martin's avatar
Ken Martin committed
606
      double minValue=VTK_DOUBLE_MAX;
Will Schroeder's avatar
Will Schroeder committed
607
      for ( i=0; i<numFullPts; i++)
Will Schroeder's avatar
Will Schroeder committed
608
        {
Will Schroeder's avatar
Will Schroeder committed
609 610
        if ( isOriginalVertex[i] )
          {
611
          if ( s[i] < minValue && s[i] <= s[(i+numFullPts-1)%numFullPts] &&
Will Schroeder's avatar
Will Schroeder committed
612
               s[i] <= s[(i+1)%numFullPts] )
Will Schroeder's avatar
Will Schroeder committed
613 614
            {
            idx = i;
Will Schroeder's avatar
Will Schroeder committed
615
            minValue = s[i];
Will Schroeder's avatar
Will Schroeder committed
616 617
            }
          }
Will Schroeder's avatar
Will Schroeder committed
618
        }
Will Schroeder's avatar
Will Schroeder committed
619

620
      //Trivial output - completely in a contour band or a triangle
621
      if ( ! intersectionPoint || numFullPts == 3 )
622
        {
Will Schroeder's avatar
Will Schroeder committed
623
        cellId = this->InsertCell(newPolys,npts,pts,cellId,s[idx],newScalars);
624 625 626
        continue;
        }

627 628 629 630 631 632 633 634 635 636 637 638 639 640 641
      //Produce contour edges if requested
      if ( this->GenerateContourEdges )
        {
        for (i=0; i < numFullPts; i++)
          {
          if ( isContourValue[i] && isContourValue[(i+1)%numFullPts] &&
               s[i] == s[(i+1)%numFullPts] )
            {
            contourEdges->InsertNextCell(2);
            contourEdges->InsertCellPoint(fullPoly[i]);
            contourEdges->InsertCellPoint(fullPoly[(i+1)%numFullPts]);
            }
          }
        }

Will Schroeder's avatar
Will Schroeder committed
642
      //Find the first intersection points in the polygons starting
Will Schroeder's avatar
Will Schroeder committed
643
      //from this vertex and build a polygon.
644
      numPointsToAdd = 1;
Will Schroeder's avatar
Will Schroeder committed
645
      for ( mR=idx, intersectionPoint=0; !intersectionPoint; )
Will Schroeder's avatar
Will Schroeder committed
646
        {
647
        numPointsToAdd++;
Will Schroeder's avatar
Will Schroeder committed
648
        mR = (mR + 1) % numFullPts;
649
        if ( isContourValue[mR] && s[mR] != s[idx] ) intersectionPoint = 1;
Will Schroeder's avatar
Will Schroeder committed
650
        }
Will Schroeder's avatar
Will Schroeder committed
651
      for ( mL=idx, intersectionPoint=0; !intersectionPoint; )
Will Schroeder's avatar
Will Schroeder committed
652
        {
653
        numPointsToAdd++;
Will Schroeder's avatar
Will Schroeder committed
654
        mL = (mL + numFullPts - 1) % numFullPts;
Will Schroeder's avatar
Will Schroeder committed
655
        if ( isContourValue[mL] && s[mL] != s[idx] ) intersectionPoint = 1;
Will Schroeder's avatar
Will Schroeder committed
656
        }
657
      for ( numPolyPoints=0, i=0; i<numPointsToAdd; i++)
Will Schroeder's avatar
Will Schroeder committed
658 659 660
        {
        newPolygon[numPolyPoints++] = fullPoly[(mL+i)%numFullPts];
        }
661 662 663 664 665
      if(numPolyPoints >= 3)
        {
        cellId = this->InsertCell(newPolys,numPolyPoints,newPolygon,
                                  cellId,s[idx],newScalars);
        }
666 667 668 669 670 671
      if ( this->GenerateContourEdges )
        {
        contourEdges->InsertNextCell(2);
        contourEdges->InsertCellPoint(fullPoly[mR]);
        contourEdges->InsertCellPoint(fullPoly[mL]);
        }
Will Schroeder's avatar
Will Schroeder committed
672 673 674

      //We've got an edge (mL,mR) that marks the edge of the region not yet
      //clipped. We move this edge forward from intersection point to
675
      //intersection point.
Will Schroeder's avatar
Will Schroeder committed
676 677 678 679
      m2R = mR;
      m2L = mL;
      while ( m2R != m2L )
        {
680
        numPointsToAdd = (mL > mR ? mL-mR+1 : numFullPts-(mR-mL)+1);
681
        if ( numPointsToAdd == 3 )
Will Schroeder's avatar
Will Schroeder committed
682
          {//just a triangle left
683
          for (i=0; i<numPointsToAdd; i++)
Will Schroeder's avatar
Will Schroeder committed
684
            {
Will Schroeder's avatar
Will Schroeder committed
685
            newPolygon[i] = fullPoly[(mR+i)%numFullPts];
Will Schroeder's avatar
Will Schroeder committed
686
            }
Will Schroeder's avatar
Will Schroeder committed
687 688
          cellId = this->InsertCell(newPolys,numPointsToAdd,newPolygon,
                                    cellId,s[mR],newScalars);
689 690 691 692 693 694
          if ( this->GenerateContourEdges )
            {
            contourEdges->InsertNextCell(2);
            contourEdges->InsertCellPoint(fullPoly[mR]);
            contourEdges->InsertCellPoint(fullPoly[mL]);
            }
Will Schroeder's avatar
Will Schroeder committed
695
          break;
Will Schroeder's avatar
Will Schroeder committed
696
          }
Will Schroeder's avatar
Will Schroeder committed
697
        else //find the next intersection points
Will Schroeder's avatar
Will Schroeder committed
698
          {
699 700
          numLeftPointsToAdd = 0;
          numRightPointsToAdd = 0;
701
          for ( intersectionPoint=0;
702
                !intersectionPoint && ((m2R+1)%numFullPts) != m2L; )
Will Schroeder's avatar
Will Schroeder committed
703
            {
704
            numRightPointsToAdd++;
Will Schroeder's avatar
Will Schroeder committed
705
            m2R = (m2R + 1) % numFullPts;
706
            if ( isContourValue[m2R] ) intersectionPoint = 1;
Will Schroeder's avatar
Will Schroeder committed
707
            }
708
          for ( intersectionPoint=0;
709
                !intersectionPoint && ((m2L+numFullPts-1)%numFullPts) != m2R; )
Will Schroeder's avatar
Will Schroeder committed
710
            {
711
            numLeftPointsToAdd++;
Will Schroeder's avatar
Will Schroeder committed
712
            m2L = (m2L + numFullPts - 1) % numFullPts;
713
            if ( isContourValue[m2L] ) intersectionPoint = 1;
Will Schroeder's avatar
Will Schroeder committed
714 715
            }

716 717
          //specify the polygon vertices. From m2L to mL, then mR to m2R.
          for ( numPolyPoints=0, i=0; i<numLeftPointsToAdd; i++)
Will Schroeder's avatar
Will Schroeder committed
718 719 720 721 722
            {
            newPolygon[numPolyPoints++] = fullPoly[(m2L+i)%numFullPts];
            }
          newPolygon[numPolyPoints++] = fullPoly[mL];
          newPolygon[numPolyPoints++] = fullPoly[mR];
723
          for ( i=1; i<=numRightPointsToAdd; i++)
Will Schroeder's avatar
Will Schroeder committed
724 725 726
            {
            newPolygon[numPolyPoints++] = fullPoly[(mR+i)%numFullPts];
            }
Will Schroeder's avatar
Will Schroeder committed
727

Will Schroeder's avatar
Will Schroeder committed
728
          //add the polygon
729 730 731 732
          if(numPolyPoints < 3)
            {
            break;
            }
Will Schroeder's avatar
Will Schroeder committed
733 734
          cellId = this->InsertCell(newPolys,numPolyPoints,newPolygon,
                                    cellId,s[mR],newScalars);
735 736 737 738 739 740
          if ( this->GenerateContourEdges )
            {
            contourEdges->InsertNextCell(2);
            contourEdges->InsertCellPoint(fullPoly[mR]);
            contourEdges->InsertCellPoint(fullPoly[mL]);
            }
Will Schroeder's avatar
Will Schroeder committed
741 742 743 744 745
          mL = m2L;
          mR = m2R;
          }//add a polygon
        }//while still removing polygons
      }//for all polygons
Will Schroeder's avatar
Will Schroeder committed
746 747 748 749 750

    delete [] s;
    delete [] newPolygon;
    delete [] isContourValue;
    delete [] isOriginalVertex;
Will Schroeder's avatar
Will Schroeder committed
751
    delete [] fullPoly;
752

Will Schroeder's avatar
Will Schroeder committed
753 754
    output->SetPolys(newPolys);
    newPolys->Delete();
Will Schroeder's avatar
Will Schroeder committed
755 756
    if ( tmpPolys ) {tmpPolys->Delete(); }
    }//for all polygons (and strips) in input
757

Will Schroeder's avatar
Will Schroeder committed
758
  vtkDebugMacro(<<"Created " << cellId << " total cells\n");
759
  vtkDebugMacro(<<"Created " << output->GetVerts()->GetNumberOfCells()
760
                << " verts\n");
761
  vtkDebugMacro(<<"Created " << output->GetLines()->GetNumberOfCells()
762
                << " lines\n");
763
  vtkDebugMacro(<<"Created " << output->GetPolys()->GetNumberOfCells()
764
                << " polys\n");
765
  vtkDebugMacro(<<"Created " << output->GetStrips()->GetNumberOfCells()
766 767 768
                << " strips\n");

  //  Update ourselves and release temporary memory
769
  //
Will Schroeder's avatar
Will Schroeder committed
770 771 772
  delete [] this->ClipValues;
  intList->Delete();
  edgeTable->Delete();
773

774 775
  output->SetPoints(newPts);
  newPts->Delete();
776

777 778
  int arrayIdx = outCD->AddArray(newScalars);
  outCD->SetActiveAttribute(arrayIdx, vtkDataSetAttributes::SCALARS);
779

780
  newScalars->Delete();
781

782
  output->Squeeze();
783 784

  return 1;
785 786
}

787
//------------------------------------------------------------------------------
788 789
vtkPolyData *vtkBandedPolyDataContourFilter::GetContourEdgesOutput()
{
790 791 792 793 794 795 796
  if (this->GetNumberOfOutputPorts() < 2)
    {
    return NULL;
    }

  return vtkPolyData::SafeDownCast(
    this->GetExecutive()->GetOutputData(1));
797
}
798

799
//------------------------------------------------------------------------------
800 801
unsigned long int vtkBandedPolyDataContourFilter::GetMTime()
{
802
  unsigned long mTime=this->Superclass::GetMTime();
803 804
  unsigned long time;

805 806
  time = this->ContourValues->GetMTime();
  mTime = ( time > mTime ? time : mTime );
807

808 809 810
  return mTime;
}

811
//------------------------------------------------------------------------------
812 813
void vtkBandedPolyDataContourFilter::PrintSelf(ostream& os, vtkIndent indent)
{
Brad King's avatar
Brad King committed
814
  this->Superclass::PrintSelf(os,indent);
815

816
  os << indent << "Generate Contour Edges: "
817 818
     << (this->GenerateContourEdges ? "On\n" : "Off\n");

819
  this->ContourValues->PrintSelf(os,indent.GetNextIndent());
Will Schroeder's avatar
Will Schroeder committed
820
  os << indent << "Clipping: " << (this->Clipping ? "On\n" : "Off\n");
821

Will Schroeder's avatar
Will Schroeder committed
822 823 824 825 826 827 828 829 830
  os << indent << "Scalar Mode: ";
  if ( this->ScalarMode == VTK_SCALAR_MODE_INDEX )
    {
    os << "INDEX\n";
    }
  else
    {
    os << "VALUE\n";
    }
831 832

  os << indent << "Clip Tolerance: " << this->ClipTolerance << "\n";
833
}