vtkBandedPolyDataContourFilter.cxx 27.8 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

Sean McBride's avatar
Sean McBride committed
30
#include <cfloat>
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
  this->Component = 0;
42 43 44 45 46 47

  this->SetNumberOfOutputPorts(2);

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

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

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

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

72 73
}

74
//------------------------------------------------------------------------------
Ken Martin's avatar
Ken Martin committed
75
int vtkBandedPolyDataContourFilter::IsContourValue(double val)
Will Schroeder's avatar
Will Schroeder committed
76 77 78 79 80
{
  int i;

  // Check to see whether a vertex is an intersection point.
  for ( i=0; i < this->NumberOfClipValues; i++)
81
  {
Will Schroeder's avatar
Will Schroeder committed
82
    if ( val == this->ClipValues[i] )
83
    {
Will Schroeder's avatar
Will Schroeder committed
84 85
      return 1;
    }
86
  }
Will Schroeder's avatar
Will Schroeder committed
87 88 89
  return 0;
}

90
//------------------------------------------------------------------------------
91 92 93 94 95 96
// Interpolate the input scalars and create intermediate points between
// v1 and v2 at the contour values.
// The point ids are returned in the edgePts array, arranged from v1 to v2 if
// v1<v2 or vice-versa.
// The input array edgePts must be large enough to hold the point ids.
// Return the number of intersection points created in edgePts.
97
int vtkBandedPolyDataContourFilter::ClipEdge(int v1, int v2,
98
                                             vtkPoints *newPts,
99 100
                                             vtkDataArray *inScalars,
                                             vtkDoubleArray *outScalars,
101
                                             vtkPointData *inPD,
102 103
                                             vtkPointData *outPD,
                                             vtkIdType edgePts[] )
104
{
105 106
  double low  = inScalars->GetTuple(v1)[this->Component];
  double high = inScalars->GetTuple(v2)[this->Component];
107

108 109
  int lowIdx  = this->ComputeScalarIndex(low);
  int highIdx = this->ComputeScalarIndex(high);
110

111

112
  if ( lowIdx == highIdx )
113
  {
114
    return 0;
115
  }
116 117 118 119 120 121 122

  double x[3], x1[3], x2[3];
  newPts->GetPoint(v1, x1);
  newPts->GetPoint(v2, x2);

  bool reverse = (v1 > v2);
  if ( low > high )
123
  {
124 125 126 127 128 129
    std::swap(low,high);
    std::swap(lowIdx,highIdx);
    std::swap(x1[0],x2[0]);
    std::swap(x1[1],x2[1]);
    std::swap(x1[2],x2[2]);
    reverse = !reverse;
130
  }
131

132 133
  int count = highIdx - lowIdx;
  for ( int i = 0; i < count; ++i )
134
  {
135 136 137 138 139 140 141 142 143 144 145
    int idx = lowIdx + 1 + i;
    double t = (this->ClipValues[idx] - low) / (high - low);
    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]);
    vtkIdType ptId = newPts->InsertNextPoint(x);
    outPD->InterpolateEdge(inPD,ptId,v1,v2,t);
    // We cannot use s1 + t*(s2-s1) as this causes a rounding error
    outScalars->InsertTuple(ptId,&this->ClipValues[idx]);
    int pos = reverse ? count-i-1 : i;
    edgePts[ pos ] = ptId;
146
  }
147
  return count;
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 156 157 158
  double v1 = *static_cast<const double*>(val1);
  double v2 = *static_cast<const double*>(val2);

  if ( v1 < v2 )
159
  {
160
    return -1;
161
  }
162
  else if ( v1 > v2 )
163
  {
164
    return 1;
165
  }
166
  else
167
  {
168
    return 0;
169
  }
170
}
171 172
}

173
//------------------------------------------------------------------------------
Will Schroeder's avatar
Will Schroeder committed
174 175
inline int vtkBandedPolyDataContourFilter::InsertCell(vtkCellArray *cells,
                                                      int npts, vtkIdType *pts,
Ken Martin's avatar
Ken Martin committed
176
                                                      int cellId, double s,
Will Schroeder's avatar
Will Schroeder committed
177
                                                      vtkFloatArray *newS)
178 179 180
{
  int idx = this->ComputeClippedIndex( s );
  if ( idx < 0 )
181
  {
182
    return cellId;
183
  }
184 185 186 187 188 189 190 191 192 193 194 195 196
  cells->InsertNextCell(npts,pts);
  return InsertNextScalar( newS, cellId, idx );
}

//------------------------------------------------------------------------------
inline int vtkBandedPolyDataContourFilter::InsertLine(vtkCellArray *cells,
                                                      vtkIdType pt1,
                                                      vtkIdType pt2,
                                                      int cellId, double s,
                                                      vtkFloatArray *newS)
{
  int idx = this->ComputeClippedIndex( s );
  if ( idx < 0 )
197
  {
198
    return cellId;
199
  }
200 201 202 203 204 205 206 207
  cells->InsertNextCell(2);
  cells->InsertCellPoint(pt1);
  cells->InsertCellPoint(pt2);
  return InsertNextScalar( newS, cellId, idx );
}

//------------------------------------------------------------------------------
int vtkBandedPolyDataContourFilter::ComputeClippedIndex(double s )
Will Schroeder's avatar
Will Schroeder committed
208
{
209
  int idx = this->ComputeScalarIndex(s+this->InternalClipTolerance);
Will Schroeder's avatar
Will Schroeder committed
210

211
  if ( !this->Clipping ||
212
       (idx >= this->ClipIndex[0] && idx < this->ClipIndex[1]) )
213
  {
214
    return idx;
215
  }
216 217
  return -1;
}
Will Schroeder's avatar
Will Schroeder committed
218

219 220 221 222 223 224
//------------------------------------------------------------------------------
int vtkBandedPolyDataContourFilter::InsertNextScalar( vtkFloatArray* scalars,
                                                      int            cellId,
                                                      int            idx )
{
  if ( idx < 0 )
225
  {
226
    return cellId;
227
  }
228 229

  if ( this->ScalarMode == VTK_SCALAR_MODE_INDEX )
230
  {
231 232
    double value = idx;
    scalars->InsertTuple(cellId++,&value);
233
  }
234
  else
235
  {
236
    scalars->InsertTuple(cellId++,&this->ClipValues[idx]);
237
  }
Will Schroeder's avatar
Will Schroeder committed
238 239 240
  return cellId;
}

241
//------------------------------------------------------------------------------
242
// Create filled contours for polydata
243 244 245 246
int vtkBandedPolyDataContourFilter::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
247
{
248 249 250 251
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

252
  // get the input and output
253 254 255 256 257
  vtkPolyData *input = vtkPolyData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

258
  vtkPointData *pd = input->GetPointData();
259
  vtkPointData *outPD = output->GetPointData();
260
  vtkCellData *outCD = output->GetCellData();
261
  vtkPoints *inPts = input->GetPoints();
262
  vtkDataArray *inScalars = pd->GetScalars();
Will Schroeder's avatar
Will Schroeder committed
263
  int abort=0;
264
  vtkPoints *newPts;
265
  int i, j, idx = 0;
266 267
  vtkIdType npts = 0;
  vtkIdType cellId=0;
268
  vtkIdType *pts = nullptr;
269
  int numEdgePts, maxCellSize;
Will Schroeder's avatar
Will Schroeder committed
270
  vtkIdType v, vR, *intPts;
271
  int intsIdx;
272 273
  vtkIdType intLoc;
  vtkIdType numIntPts, intsInc;
274
  vtkIdType numPts, numCells, estimatedSize;
275

Will Schroeder's avatar
Will Schroeder committed
276
  vtkDebugMacro(<<"Executing banded contour filter");
277 278 279

  //  Check input
  //
280

281
  numCells = input->GetNumberOfCells();
282
  if ( !inPts || (numPts=inPts->GetNumberOfPoints()) < 1 ||
283
       !inScalars || numCells < 1 )
284
  {
285
    vtkErrorMacro(<<"No input data!");
286
    return 1;
287
  }
288

289
  if ( inScalars->GetNumberOfComponents() < this->Component + 1 )
290
  {
291 292 293
    vtkErrorMacro( << "Input scalars expected to have "<<this->Component+1
                   << " components" );
    return 0;
294
  }
295

296 297 298 299
  // 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
300
  this->ClipValues = new double[this->NumberOfClipValues];
Ken Martin's avatar
Ken Martin committed
301
  double range[2];
302 303 304
  inScalars->GetRange(range);

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

307
  this->ClipValues[0] =
308 309 310
    (range[0]<this->ContourValues->GetValue(0))?
    (range[0]):
    (this->ContourValues->GetValue(0));
311 312

  this->ClipValues[this->NumberOfClipValues - 1] =
313
    (range[1]>this->ContourValues->GetValue(this->NumberOfClipValues-2-1))?
314
    (range[1]):
315
    (this->ContourValues->GetValue(this->NumberOfClipValues-2-1));
316

317
  for ( i=1; i<this->NumberOfClipValues-1; i++)
318
  {
319
    this->ClipValues[i] = this->ContourValues->GetValue(i-1);
320
  }
Will Schroeder's avatar
Will Schroeder committed
321

322
  qsort(this->ClipValues, this->NumberOfClipValues, sizeof(double),
Will Schroeder's avatar
Will Schroeder committed
323 324
        vtkCompareClipValues);

325 326 327
  // 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...
328
  for ( i=0; i<(this->NumberOfClipValues-1); i++)
329
  {
330
    if ( (this->ClipValues[i] + this->InternalClipTolerance) >= this->ClipValues[i+1] )
331
    {
332
      for (j=i+1; j<(this->NumberOfClipValues-2); j++)
333
      {
334
        this->ClipValues[j] = this->ClipValues[j+1];
335
      }
336
      this->NumberOfClipValues--;
337
    }
338
  }
339

340
  this->ClipIndex[0] =
341 342
    this->ComputeScalarIndex(this->ContourValues->GetValue(0));
  this->ClipIndex[1] = this->ComputeScalarIndex(
Will Schroeder's avatar
Will Schroeder committed
343 344
    this->ContourValues->GetValue(this->ContourValues->GetNumberOfContours()-1));

345 346 347
  //
  // Estimate allocation size, stolen from vtkContourGrid...
  //
348
  estimatedSize=static_cast<vtkIdType>(pow(static_cast<double>(numCells),.9));
349 350 351
  estimatedSize *= this->NumberOfClipValues;
  estimatedSize = estimatedSize / 1024 * 1024; // multiple of 1024
  if (estimatedSize < 1024)
352
  {
353
    estimatedSize = 1024;
354
  }
355

356
  // The original set of points and point data are copied. Later on
357 358
  // intersection points due to clipping will be created.
  newPts = vtkPoints::New();
359

360 361 362
  // 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.
363
  newPts->Allocate(estimatedSize,estimatedSize);
364
  outPD->CopyScalarsOff();
365
  outPD->InterpolateAllocate(pd,3*numPts,numPts);
366 367 368
  vtkDoubleArray *outScalars = vtkDoubleArray::New();
  outScalars->Allocate(3*numPts,numPts);
  outPD->SetScalars(outScalars);
Berk Geveci's avatar
Berk Geveci committed
369
  outScalars->Delete();
370

371
  for (i=0; i<numPts; i++)
372
  {
373 374
    newPts->InsertPoint(i,inPts->GetPoint(i));
    outPD->CopyData(pd, i, i);
375 376
    double value = inScalars->GetTuple(i)[this->Component];
    outScalars->InsertTuple(i, &value);
377
  }
378

379 380 381
  // These are the new cell scalars
  vtkFloatArray *newScalars = vtkFloatArray::New();
  newScalars->Allocate(numCells*5,numCells);
382
  newScalars->SetName("Scalars");
383

Will Schroeder's avatar
Will Schroeder committed
384 385 386 387
  // Used to keep track of intersections
  vtkEdgeTable *edgeTable = vtkEdgeTable::New();
  vtkCellArray *intList = vtkCellArray::New(); //intersection point ids

388 389
  // All vertices are filled and passed through; poly-vertices are broken
  // into single vertices. Cell data per vertex is set.
390
  //
391
  if ( input->GetVerts()->GetNumberOfCells() > 0 )
392
  {
393 394 395
    vtkCellArray *verts = input->GetVerts();
    vtkCellArray *newVerts = vtkCellArray::New();
    newVerts->Allocate(verts->GetSize());
396
    for ( verts->InitTraversal(); verts->GetNextCell(npts,pts) && !abort;
397
          abort=this->GetAbortExecute() )
398
    {
399
      for (i=0; i<npts; i++)
400
      {
401 402 403 404
        cellId = this->InsertCell( newVerts,
                                   1,pts+i,cellId,
                                   inScalars->GetTuple(pts[i])[this->Component],
                                   newScalars);
405
      }
406
    }
407 408
    output->SetVerts(newVerts);
    newVerts->Delete();
409
  }
410
  this->UpdateProgress(0.05);
411

412 413 414
  // Lines are chopped into line segments.
  //
  if ( input->GetLines()->GetNumberOfCells() > 0 )
415
  {
416
    vtkCellArray *lines = input->GetLines();
Will Schroeder's avatar
Will Schroeder committed
417 418 419 420 421

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

    vtkIdType *fullLine = new vtkIdType [maxCellSize];
422 423
    vtkCellArray *newLines = vtkCellArray::New();
    newLines->Allocate(lines->GetSize());
Will Schroeder's avatar
Will Schroeder committed
424 425 426
    edgeTable->InitEdgeInsertion(numPts,1); //store attributes on edge

    //start by generating intersection points
427
    for ( lines->InitTraversal(); lines->GetNextCell(npts,pts) && !abort;
428
          abort=this->GetAbortExecute() )
429
    {
430
      for (i=0; i<(npts-1); i++)
431
      {
432 433
        numEdgePts = this->ClipEdge(pts[i],pts[i+1],newPts,inScalars,outScalars,
                                    pd,outPD, fullLine );
Will Schroeder's avatar
Will Schroeder committed
434
        if ( numEdgePts > 0 ) //there is an intersection
435
        {
Will Schroeder's avatar
Will Schroeder committed
436 437 438
          intList->InsertNextCell(numEdgePts,fullLine);
          edgeTable->InsertEdge(pts[i],pts[i+1], //associate ints with edge
                                intList->GetInsertLocation(numEdgePts));
439
        }
Will Schroeder's avatar
Will Schroeder committed
440
        else //no intersection points along the edge
441
        {
Will Schroeder's avatar
Will Schroeder committed
442
          edgeTable->InsertEdge(pts[i],pts[i+1],-1); //-1 means no points
443 444 445
        }
      }//for all line segments in this line
    }
Will Schroeder's avatar
Will Schroeder committed
446 447

    //now create line segments
448
    for ( lines->InitTraversal(); lines->GetNextCell(npts,pts) && !abort;
449
          abort=this->GetAbortExecute() )
450
    {
Will Schroeder's avatar
Will Schroeder committed
451
      for (i=0; i<(npts-1); i++)
452
      {
Will Schroeder's avatar
Will Schroeder committed
453 454
        v = pts[i];
        vR = pts[i+1];
455
        bool reverse = (v > vR );
Will Schroeder's avatar
Will Schroeder committed
456

457 458 459
        double s1 = inScalars->GetTuple(v)[this->Component];
        double s2 = inScalars->GetTuple(vR)[this->Component];
        bool increasing = ( s2 > s1 );
Will Schroeder's avatar
Will Schroeder committed
460

461
        vtkIdType p1=v;
Will Schroeder's avatar
Will Schroeder committed
462
        if ( (intLoc=edgeTable->IsEdge(v,vR)) != -1 )
463
        {
Will Schroeder's avatar
Will Schroeder committed
464
          intList->GetCell(intLoc,numIntPts,intPts);
465 466 467
          int incr;
          int k;
          if ( !reverse )
468
          {
469 470
            k = 0;
            incr = 1;
471
          }
472
          else
473
          {
474 475
            k = numIntPts-1;
            incr = -1;
476
          }
477
          for ( int n = 0; n < numIntPts; ++n, k += incr )
478
          {
479 480 481 482 483
            vtkIdType p2 = intPts[k];
            double value = outScalars->GetTuple(increasing ? p1 : p2)[0];
            cellId = this->InsertLine(newLines, p1, p2, cellId,
                                      value, newScalars);
            p1 = p2;
484
          }
485 486 487
          double value = outScalars->GetTuple(increasing ? p1 : vR)[0];
          cellId = this->InsertLine(newLines, p1, vR, cellId,
                                    value, newScalars);
488
        }
489
        else
490
        {
491 492 493
          double value = outScalars->GetTuple(vR)[0];
          cellId = this->InsertLine(newLines, v, vR, cellId,
                                    value, newScalars);
494
        }
495
      }
496
    }
Will Schroeder's avatar
Will Schroeder committed
497 498 499

    delete [] fullLine;

500 501
    output->SetLines(newLines);
    newLines->Delete();
502
  }
503
  this->UpdateProgress(0.1);
504

505
  // Polygons are assumed convex and chopped into filled, convex polygons.
506
  // Triangle strips are treated similarly.
507
  //
Sean McBride's avatar
Sean McBride committed
508 509
  vtkIdType numPolys = input->GetPolys()->GetNumberOfCells();
  vtkIdType numStrips = input->GetStrips()->GetNumberOfCells();
Will Schroeder's avatar
Will Schroeder committed
510
  if ( numPolys > 0 || numStrips > 0 )
511
  {
512 513 514 515 516
    // 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
517
    edgeTable->InitEdgeInsertion(numPts,1); //store attributes on edge
Will Schroeder's avatar
Will Schroeder committed
518
    intList->Reset();
Will Schroeder's avatar
Will Schroeder committed
519 520

    vtkCellArray *polys = input->GetPolys();
521
    vtkCellArray *tmpPolys = nullptr;
Will Schroeder's avatar
Will Schroeder committed
522

523
    // If contour edges requested, set things up.
524
    vtkCellArray *contourEdges=nullptr;
525
    if ( this->GenerateContourEdges )
526
    {
527 528 529 530 531
      contourEdges = vtkCellArray::New();
      contourEdges->Allocate(numCells);
      this->GetContourEdgesOutput()->SetLines(contourEdges);
      contourEdges->Delete();
      this->GetContourEdgesOutput()->SetPoints(newPts);
532
    }
533

Will Schroeder's avatar
Will Schroeder committed
534 535
    // Set up structures for processing polygons
    maxCellSize = polys->GetMaxCellSize();
536
    if( maxCellSize == 0 )
537
    {
538
      maxCellSize = input->GetStrips()->GetMaxCellSize();
539
    }
Will Schroeder's avatar
Will Schroeder committed
540 541 542
    maxCellSize *= (1 + this->NumberOfClipValues);

    vtkIdType *newPolygon = new vtkIdType [maxCellSize];
Ken Martin's avatar
Ken Martin committed
543
    double *s = new double [maxCellSize]; //scalars at vertices
Will Schroeder's avatar
Will Schroeder committed
544 545 546 547 548 549
    int *isContourValue = new int [maxCellSize];
    int *isOriginalVertex = new int [maxCellSize];
    vtkIdType *fullPoly = new vtkIdType [maxCellSize];

    // Lump strips and polygons together.
    // Decompose strips into triangles.
550
    if ( numStrips > 0 )
551
    {
Will Schroeder's avatar
Will Schroeder committed
552 553 554
      vtkCellArray *strips = input->GetStrips();
      tmpPolys = vtkCellArray::New();
      if ( numPolys > 0 )
555
      {
Will Schroeder's avatar
Will Schroeder committed
556
        tmpPolys->DeepCopy(polys);
557
      }
558
      else
559
      {
Will Schroeder's avatar
Will Schroeder committed
560
        tmpPolys->Allocate(polys->EstimateSize(numStrips,5));
561
      }
Will Schroeder's avatar
Will Schroeder committed
562
      for ( strips->InitTraversal(); strips->GetNextCell(npts,pts); )
563
      {
Will Schroeder's avatar
Will Schroeder committed
564 565
        vtkTriangleStrip::DecomposeStrip(npts, pts, tmpPolys);
      }
566 567
      polys = tmpPolys;
    }
568

Will Schroeder's avatar
Will Schroeder committed
569 570
    // Process polygons to produce edge intersections.------------------------
    //
571 572 573
    numPolys = polys->GetNumberOfCells();
    vtkIdType updateCount = numPolys/20 + 1;
    vtkIdType count=0;
574
    for ( polys->InitTraversal(); polys->GetNextCell(npts,pts) && !abort;
575
      abort=this->GetAbortExecute() )
576
    {
577
      if  ( ! (++count % updateCount) )
578
      {
579
        this->UpdateProgress(0.1 + 0.45*(static_cast<double>(count)/numPolys));
580
      }
581

Will Schroeder's avatar
Will Schroeder committed
582
      for (i=0; i<npts; i++)
583
      {
Will Schroeder's avatar
Will Schroeder committed
584 585 586
        v = pts[i];
        vR = pts[(i+1) % npts];
        if ( edgeTable->IsEdge(v,vR) == -1 )
587
        {
588 589
          numEdgePts = this->ClipEdge(v,vR,newPts,inScalars,outScalars,
                                      pd,outPD,fullPoly);
Will Schroeder's avatar
Will Schroeder committed
590
          if ( numEdgePts > 0 )
591
          {
Will Schroeder's avatar
Will Schroeder committed
592 593 594
            intList->InsertNextCell(numEdgePts,fullPoly);
            edgeTable->InsertEdge(v,vR, //associate ints with edge
                                  intList->GetInsertLocation(numEdgePts));
595
          }
Will Schroeder's avatar
Will Schroeder committed
596
          else //no intersection points along the edge
597
          {
Will Schroeder's avatar
Will Schroeder committed
598
            edgeTable->InsertEdge(v,vR,-1); //-1 means no points
599 600 601 602
          }
        }//if edge not processed yet
      }
    }//for all polygons
603

604
    // Process polygons to produce output triangles------------------------
Will Schroeder's avatar
Will Schroeder committed
605
    //
Will Schroeder's avatar
Will Schroeder committed
606 607
    vtkCellArray *newPolys = vtkCellArray::New();
    newPolys->Allocate(polys->GetSize());
Will Schroeder's avatar
Will Schroeder committed
608
    int intersectionPoint;
609 610
    int mL, mR, m2L, m2R;
    int numPointsToAdd, numLeftPointsToAdd, numRightPointsToAdd;
Will Schroeder's avatar
Will Schroeder committed
611
    int numPolyPoints, numFullPts;
612
    count = 0;
613
    for ( polys->InitTraversal(); polys->GetNextCell(npts,pts) && !abort;
614
          abort=this->GetAbortExecute() )
615
    {
616
      if  ( ! (++count % updateCount) )
617
      {
618 619
        this->UpdateProgress(0.55 +
                             0.45*(static_cast<double>(count)/numPolys));
620
      }
621

Will Schroeder's avatar
Will Schroeder committed
622
      //Create a new polygon that includes all the points including the
Will Schroeder's avatar
Will Schroeder committed
623
      //intersection vertices. This hugely simplifies the logic of the
Will Schroeder's avatar
Will Schroeder committed
624 625
      //code.
      for ( intersectionPoint=0, numFullPts=0, i=0; i<npts; i++)
626
      {
Will Schroeder's avatar
Will Schroeder committed
627 628
        v = pts[i];
        vR = pts[(i+1)%npts];
629

630
        s[numFullPts] = outScalars->GetTuple(v)[0];
Will Schroeder's avatar
Will Schroeder committed
631
        isContourValue[numFullPts] = this->IsContourValue(s[numFullPts]);
Will Schroeder's avatar
Will Schroeder committed
632
        isOriginalVertex[numFullPts] = 1;
Will Schroeder's avatar
Will Schroeder committed
633
        fullPoly[numFullPts++] = v;
Will Schroeder's avatar
Will Schroeder committed
634

Will Schroeder's avatar
Will Schroeder committed
635
        //see whether intersection points need to be added.
Will Schroeder's avatar
Will Schroeder committed
636
        if ( (intLoc=edgeTable->IsEdge(v,vR)) != -1 )
637
        {
Will Schroeder's avatar
Will Schroeder committed
638 639
          intersectionPoint = 1;
          intList->GetCell(intLoc,numIntPts,intPts);
Will Schroeder's avatar
Will Schroeder committed
640
          if ( v < vR ) {intsIdx = 0; intsInc=1;} //order of the edge
Will Schroeder's avatar
Will Schroeder committed
641
          else {intsIdx=numIntPts-1; intsInc=(-1);}
Will Schroeder's avatar
Will Schroeder committed
642
          for ( ; intsIdx >= 0 && intsIdx < numIntPts; intsIdx += intsInc )
643
          {
644
            s[numFullPts] = outScalars->GetTuple(intPts[intsIdx])[0];
Will Schroeder's avatar
Will Schroeder committed
645 646 647 648
            isContourValue[numFullPts] = 1;
            isOriginalVertex[numFullPts] = 0;
            fullPoly[numFullPts++] = intPts[intsIdx];
          }
649 650
        }
      } //for all points and edges
651

Will Schroeder's avatar
Will Schroeder committed
652 653
      //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
654
      //Really should check whether the vertex is convex.
Ken Martin's avatar
Ken Martin committed
655
      double minValue=VTK_DOUBLE_MAX;
Will Schroeder's avatar
Will Schroeder committed
656
      for ( i=0; i<numFullPts; i++)
657
      {
Will Schroeder's avatar
Will Schroeder committed
658
        if ( isOriginalVertex[i] )
659
        {
660
          if ( s[i] < minValue && s[i] <= s[(i+numFullPts-1)%numFullPts] &&
Will Schroeder's avatar
Will Schroeder committed
661
               s[i] <= s[(i+1)%numFullPts] )
662
          {
Will Schroeder's avatar
Will Schroeder committed
663
            idx = i;
Will Schroeder's avatar
Will Schroeder committed
664
            minValue = s[i];
Will Schroeder's avatar
Will Schroeder committed
665
          }
Will Schroeder's avatar
Will Schroeder committed
666
        }
667
      }
Will Schroeder's avatar
Will Schroeder committed
668

669
      //Trivial output - completely in a contour band or a triangle
670
      if ( ! intersectionPoint || numFullPts == 3 )
671
      {
Will Schroeder's avatar
Will Schroeder committed
672
        cellId = this->InsertCell(newPolys,npts,pts,cellId,s[idx],newScalars);
673
        continue;
674
      }
675

676 677
      //Produce contour edges if requested
      if ( this->GenerateContourEdges )
678
      {
679
        for (i=0; i < numFullPts; i++)
680
        {
681 682
          if ( isContourValue[i] && isContourValue[(i+1)%numFullPts] &&
               s[i] == s[(i+1)%numFullPts] )
683
          {
684 685 686 687 688
            contourEdges->InsertNextCell(2);
            contourEdges->InsertCellPoint(fullPoly[i]);
            contourEdges->InsertCellPoint(fullPoly[(i+1)%numFullPts]);
          }
        }
689
      }
690

Will Schroeder's avatar
Will Schroeder committed
691
      //Find the first intersection points in the polygons starting
Will Schroeder's avatar
Will Schroeder committed
692
      //from this vertex and build a polygon.
693
      numPointsToAdd = 1;
Will Schroeder's avatar
Will Schroeder committed
694
      for ( mR=idx, intersectionPoint=0; !intersectionPoint; )
695
      {
696
        numPointsToAdd++;
Will Schroeder's avatar
Will Schroeder committed
697
        mR = (mR + 1) % numFullPts;
698
        if ( isContourValue[mR] && s[mR] != s[idx] ) intersectionPoint = 1;
699
      }
Will Schroeder's avatar
Will Schroeder committed
700
      for ( mL=idx, intersectionPoint=0; !intersectionPoint; )
701
      {
702
        numPointsToAdd++;
Will Schroeder's avatar
Will Schroeder committed
703
        mL = (mL + numFullPts - 1) % numFullPts;
Will Schroeder's avatar
Will Schroeder committed
704
        if ( isContourValue[mL] && s[mL] != s[idx] ) intersectionPoint = 1;
705
      }
706
      for ( numPolyPoints=0, i=0; i<numPointsToAdd; i++)
707
      {
Will Schroeder's avatar
Will Schroeder committed
708
        newPolygon[numPolyPoints++] = fullPoly[(mL+i)%numFullPts];
709
      }
710
      if(numPolyPoints >= 3)
711
      {
712 713
        cellId = this->InsertCell(newPolys,numPolyPoints,newPolygon,
                                  cellId,s[idx],newScalars);
714
      }
715
      if ( this->GenerateContourEdges )
716
      {
717 718 719
        contourEdges->InsertNextCell(2);
        contourEdges->InsertCellPoint(fullPoly[mR]);
        contourEdges->InsertCellPoint(fullPoly[mL]);
720
      }
Will Schroeder's avatar
Will Schroeder committed
721 722 723

      //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
724
      //intersection point.
Will Schroeder's avatar
Will Schroeder committed
725 726 727
      m2R = mR;
      m2L = mL;
      while ( m2R != m2L )
728
      {
729
        numPointsToAdd = (mL > mR ? mL-mR+1 : numFullPts-(mR-mL)+1);
730
        if ( numPointsToAdd == 3 )
731
        {//just a triangle left
732
          for (i=0; i<numPointsToAdd; i++)
733
          {
Will Schroeder's avatar
Will Schroeder committed
734
            newPolygon[i] = fullPoly[(mR+i)%numFullPts];
735
          }
Will Schroeder's avatar
Will Schroeder committed
736 737
          cellId = this->InsertCell(newPolys,numPointsToAdd,newPolygon,
                                    cellId,s[mR],newScalars);
738
          if ( this->GenerateContourEdges )
739
          {
740 741 742
            contourEdges->InsertNextCell(2);
            contourEdges->InsertCellPoint(fullPoly[mR]);
            contourEdges->InsertCellPoint(fullPoly[mL]);
Will Schroeder's avatar
Will Schroeder committed
743
          }
744 745
          break;
        }
Will Schroeder's avatar
Will Schroeder committed
746
        else //find the next intersection points
747
        {
748 749
          numLeftPointsToAdd = 0;
          numRightPointsToAdd = 0;
750
          for ( intersectionPoint=0;
751
                !intersectionPoint && ((m2R+1)%numFullPts) != m2L; )
752
          {
753
            numRightPointsToAdd++;
Will Schroeder's avatar
Will Schroeder committed
754
            m2R = (m2R + 1) % numFullPts;
Mathieu Malaterre's avatar <