vtkFlyingEdges2D.cxx 31.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkFlyingEdges2D.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 "vtkFlyingEdges2D.h"

#include "vtkImageData.h"
#include "vtkCellArray.h"
19
#include "vtkImageTransform.h"
20 21 22 23 24 25 26
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkFloatArray.h"
#include "vtkStreamingDemandDrivenPipeline.h"
27
#include "vtkSMPTools.h"
28

Sean McBride's avatar
Sean McBride committed
29
#include <cmath>
30 31 32 33 34 35 36 37 38 39 40

vtkStandardNewMacro(vtkFlyingEdges2D);

// This templated class is the heart of the algorithm. Templated across
// scalar type T. vtkFlyingEdges2D populates the information in this class
// and then invokes ContourImage() to actually initiate executions.
template <class T>
class vtkFlyingEdges2DAlgorithm
{
public:
  // Edge case table values.
Will Schroeder's avatar
Will Schroeder committed
41
  enum EdgeClass {
42 43 44 45 46 47
    Below = 0, //below isovalue
    Above = 1, //above isovalue
    LeftAbove = 1, //left vertex is above isovalue
    RightAbove = 2, //right vertex is above isovalue
    BothAbove = 3 //entire edge is above isovalue
  };
48 49

  // Dealing with boundary situations when processing images.
Will Schroeder's avatar
Will Schroeder committed
50
  enum CellClass {
51 52 53 54
    Interior = 0,
    MinBoundary = 1,
    MaxBoundary = 2
  };
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
  // Edges to generate output line primitives (aka case table).
  static const unsigned char EdgeCases[16][5];

  // A table that lists pixel point ids as a function of edge ids (edge ids
  // for edge-based case table).
  static const unsigned char VertMap[4][2];

  // A table describing vertex offsets (in index space) from the pixel axes
  // origin for each of the four vertices of a pixel.
  static const unsigned char VertOffsets[4][2];

  // This table is used to accelerate the generation of output lines and
  // points. The EdgeUses array, a function of the pixel case number,
  // indicates which pixel edges intersect with the contour (i.e., require
  // interpolation). This array is filled in at instantiation during the case
  // table generation process.
  unsigned char EdgeUses[16][4];

  // Flags indicate whether a particular case requires pixel axes to be
  // processed. A cheap acceleration structure computed from the case
  // tables at the point of instantiation.
  unsigned char IncludesAxes[16];

  // Algorithm-derived data
  unsigned char *XCases;
  vtkIdType *EdgeMetaData;

  // Internal variables used by the various algorithm methods. Interfaces VTK
  // image data in a form more convenient to the algorithm.
  vtkIdType Dims[2];
85
  int K;
86 87 88 89 90 91 92 93 94 95 96
  int Axis0;
  int Min0;
  int Max0;
  int Inc0;
  int Axis1;
  int Min1;
  int Max1;
  int Inc1;
  int Axis2;

  // Output data. Threads write to partitioned memory.
97
  T         *Scalars;
98 99 100 101 102 103 104 105 106 107
  T         *NewScalars;
  vtkIdType *NewLines;
  float     *NewPoints;

  // Instantiate and initialize key data members.
  vtkFlyingEdges2DAlgorithm();

  // The three passes of the algorithm.
  void ProcessXEdge(double value, T* inPtr, vtkIdType row); //PASS 1
  void ProcessYEdges(vtkIdType row); //PASS 2
108
  void GenerateOutput(double value, T* inPtr, vtkIdType row); //PASS 4
109 110 111 112 113 114 115 116

  // Place holder for now in case fancy bit fiddling is needed later.
  void SetXEdge(unsigned char *ePtr, unsigned char edgeCase)
    {*ePtr = edgeCase;}

  // Given the two x-edge cases defining this pixel, return the pixel case
  // number.
  unsigned char GetEdgeCase(unsigned char *ePtr0, unsigned char *ePtr1)
117
  {
118
    return ( (*ePtr0) | ((*ePtr1)<<2) );
119
  }
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135

  // Return number of contouring primitives (line segments) for a particular case.
  unsigned char GetNumberOfPrimitives(unsigned char caseNum)
    { return this->EdgeCases[caseNum][0]; }

  // Return an array indicating which pixel edges intersect the contour.
  unsigned char *GetEdgeUses(unsigned char eCase)
    { return this->EdgeUses[eCase]; }

  // Indicate whether pixel axes need processing for this case.
  unsigned char CaseIncludesAxes(unsigned char eCase)
    { return this->IncludesAxes[eCase]; }

  // Count edge intersections near image boundaries.
  void CountBoundaryYInts(unsigned char loc, unsigned char *edgeUses,
                          vtkIdType *eMD)
136
  {
137
      if ( loc == 2 )
138
      { //+x boundary
139
        eMD[1] += edgeUses[3];
140 141
      }
  }
142 143 144 145

  // Produce the line segments for this pixel cell.
  void GenerateLines(unsigned char eCase, unsigned char numLines, vtkIdType *eIds,
                    vtkIdType &lineId)
146
  {
147 148 149
      vtkIdType *line;
      const unsigned char *edges = this->EdgeCases[eCase] + 1;
      for (int i=0; i < numLines; ++i, edges+=2)
150
      {
151 152 153 154
        line = this->NewLines + 3*lineId++;
        line[0] = 2;
        line[1] = eIds[edges[0]];
        line[2] = eIds[edges[1]];
155 156
      }
  }
157 158

  // Interpolate along a pixel axes edge.
159 160
  void InterpolateAxesEdge(double value, T *s0, int ijk0[3], T* s1,
                           int ijk1[3], vtkIdType vId)
161
  {
162 163
      double t = (value - *s0) / (*s1 - *s0);
      float *x = this->NewPoints + 3*vId;
164 165 166
      x[0] = ijk0[0] + t*(ijk1[0]-ijk0[0]) + this->Min0;
      x[1] = ijk0[1] + t*(ijk1[1]-ijk0[1]) + this->Min1;
      x[2] = this->K;
167
  }
168 169 170 171

  // Interpolate along an arbitrary edge, typically one that may be on the
  // volume boundary. This means careful computation of stuff requiring
  // neighborhood information (e.g., gradients).
172
  void InterpolateEdge(double value, T *s, int ijk[3], unsigned char edgeNum,
Will Schroeder's avatar
Will Schroeder committed
173
                       unsigned char edgeUses[4], vtkIdType *eIds);
174 175

  // Produce the output points on the pixel axes for this pixel cell.
176
  void GeneratePoints(double value, unsigned char loc, T *sPtr, int ijk[3],
Will Schroeder's avatar
Will Schroeder committed
177
                      unsigned char *edgeUses, vtkIdType *eIds);
178 179 180 181

  // Helper function to set up the point ids on pixel edges.
  unsigned char InitPixelIds(unsigned char *ePtr0, unsigned char *ePtr1,
                             vtkIdType *eMD0, vtkIdType *eMD1, vtkIdType *eIds)
182
  {
183 184 185 186 187 188
      unsigned char eCase = GetEdgeCase(ePtr0,ePtr1);
      eIds[0] = eMD0[0]; //x-edges
      eIds[1] = eMD1[0];
      eIds[2] = eMD0[1]; //y-edges
      eIds[3] = eIds[2] + this->EdgeUses[eCase][2];
      return eCase;
189
  }
190 191 192

  // Helper function to advance the point ids along pixel rows.
  void AdvancePixelIds(unsigned char eCase, vtkIdType *eIds)
193
  {
194 195 196 197
      eIds[0] += this->EdgeUses[eCase][0]; //x-edges
      eIds[1] += this->EdgeUses[eCase][1];
      eIds[2] += this->EdgeUses[eCase][2]; //y-edges
      eIds[3] = eIds[2] + this->EdgeUses[eCase][3];
198
  }
199 200 201

  // Threading integration via SMPTools
  template <class TT> class Pass1
202
  {
203 204 205 206 207 208
    public:
      Pass1(vtkFlyingEdges2DAlgorithm<TT> *algo, double value)
        { this->Algo = algo; this->Value = value;}
      vtkFlyingEdges2DAlgorithm<TT> *Algo;
      double Value;
      void  operator()(vtkIdType row, vtkIdType end)
209
      {
210 211
        TT *rowPtr = this->Algo->Scalars + row*this->Algo->Inc1;
        for ( ; row < end; ++row)
212
        {
213 214
          this->Algo->ProcessXEdge(this->Value, rowPtr, row);
          rowPtr += this->Algo->Inc1;
215 216 217
        }//for all rows in this batch
      }
  };
218
  template <class TT> class Pass2
219
  {
220 221 222 223 224
    public:
      Pass2(vtkFlyingEdges2DAlgorithm<TT> *algo)
        { this->Algo = algo;}
      vtkFlyingEdges2DAlgorithm<TT> *Algo;
      void  operator()(vtkIdType row, vtkIdType end)
225
      {
226
        for ( ; row < end; ++row)
227
        {
228
          this->Algo->ProcessYEdges(row);
229 230 231
        }//for all rows in this batch
      }
  };
232
  template <class TT> class Pass4
233
  {
234
    public:
235
      Pass4(vtkFlyingEdges2DAlgorithm<TT> *algo, double value)
236 237 238 239
        { this->Algo = algo; this->Value = value;}
      vtkFlyingEdges2DAlgorithm<TT> *Algo;
      double Value;
      void  operator()(vtkIdType row, vtkIdType end)
240
      {
241 242
        T *rowPtr = this->Algo->Scalars + row*this->Algo->Inc1;
        for ( ; row < end; ++row)
243
        {
244 245
          this->Algo->GenerateOutput(this->Value, rowPtr, row);
          rowPtr += this->Algo->Inc1;
246 247 248
        }//for all rows in this batch
      }
  };
249 250 251 252 253

  // Interface between VTK and templated functions
  static void ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
                           vtkDataArray *newScalars, vtkCellArray *newLines,
                           vtkImageData *input, int *updateExt);
254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
};

//----------------------------------------------------------------------------
// Specify the points that define each edge.
template <class T> const unsigned char vtkFlyingEdges2DAlgorithm<T>::
VertMap[4][2] = {{0,1}, {2,3}, {0,2}, {1,3}};

//----------------------------------------------------------------------------
// The offsets of each vertex (in index space) from the pixel axes origin.
template <class T> const unsigned char vtkFlyingEdges2DAlgorithm<T>::
VertOffsets[4][2] = {{0,0}, {1,0}, {0,1}, {1,1}};

// Initialize case array
template <class T> const unsigned char vtkFlyingEdges2DAlgorithm<T>::
EdgeCases[16][5] = {
    {0,0,0,0,0}, {1,0,2,0,0}, {1,3,0,0,0}, {1,3,2,0,0},
    {1,2,1,0,0}, {1,0,1,0,0}, {2,2,1,3,0}, {1,3,1,0,0},
    {1,1,3,0,0}, {2,0,2,3,1}, {1,1,0,0,0}, {1,1,2,0,0},
    {1,2,3,0,0}, {1,0,3,0,0}, {1,2,0,0,0}, {0,0,0,0,0}};

//----------------------------------------------------------------------------
// Instantiate and initialize key data members. Mostly we build some
// acceleration structures from the case table.
template <class T> vtkFlyingEdges2DAlgorithm<T>::
278 279
vtkFlyingEdges2DAlgorithm():XCases(nullptr),EdgeMetaData(nullptr),Scalars(nullptr),
                            NewScalars(nullptr),NewLines(nullptr),NewPoints(nullptr)
280 281 282 283 284 285
{
  int j, eCase, numLines;
  const unsigned char *edgeCase;

  // Initialize
  for (eCase=0; eCase<16; ++eCase)
286
  {
287
    for (j=0; j<4; ++j)
288
    {
289 290
      this->EdgeUses[eCase][j] = 0;
    }
291 292
    this->IncludesAxes[eCase] = 0;
  }
293 294 295

  // Populate
  for (eCase=0; eCase<16; ++eCase)
296
  {
297 298 299 300 301
    edgeCase = this->EdgeCases[eCase];
    numLines = *edgeCase++;

    // Mark edges that are used by this case.
    for (j=0; j < numLines*2; ++j) //just loop over all edges
302
    {
303
      this->EdgeUses[eCase][edgeCase[j]] = 1;
304
    }
305 306 307

    this->IncludesAxes[eCase] = this->EdgeUses[eCase][0] |
      this->EdgeUses[eCase][2];
308
  }
309 310 311 312 313 314
}

//----------------------------------------------------------------------------
// Interpolate a new point along a boundary edge. Make sure to consider
// proximity to boundary when computing gradients, etc.
template <class T> void vtkFlyingEdges2DAlgorithm<T>::
315
InterpolateEdge(double value, T *s, int ijk[3], unsigned char edgeNum,
Will Schroeder's avatar
Will Schroeder committed
316
                unsigned char edgeUses[12], vtkIdType *eIds)
317 318 319
{
  // if this edge is not used then get out
  if ( ! edgeUses[edgeNum] )
320
  {
321
    return;
322
  }
323 324 325 326 327

  // build the edge information
  const unsigned char *vertMap = this->VertMap[edgeNum];
  T *s0, *s1;
  float x0[3], x1[3];
Will Schroeder's avatar
Will Schroeder committed
328
  vtkIdType vId=eIds[edgeNum];
329 330 331

  const unsigned char *offsets = this->VertOffsets[vertMap[0]];
  s0 = s + offsets[0]*this->Inc0 + offsets[1]*this->Inc1;
332 333
  x0[0] = ijk[0] + offsets[0];
  x0[1] = ijk[1] + offsets[1];
334 335 336

  offsets = this->VertOffsets[vertMap[1]];
  s1 = s + offsets[0]*this->Inc0 + offsets[1]*this->Inc1;
337 338
  x1[0] = ijk[0] + offsets[0];
  x1[1] = ijk[1] + offsets[1];
339 340 341 342

  // Okay interpolate
  double t = (value - *s0) / (*s1 - *s0);
  float *xPtr = this->NewPoints + 3*vId;
343 344 345
  xPtr[0] = x0[0] + t*(x1[0]-x0[0]) + this->Min0;
  xPtr[1] = x0[1] + t*(x1[1]-x0[1]) + this->Min1;
  xPtr[2] = this->K;
346 347 348 349 350 351
}

//----------------------------------------------------------------------------
// Generate the output points and optionally normals, gradients and
// interpolate attributes.
template <class T> void vtkFlyingEdges2DAlgorithm<T>::
352
GeneratePoints(double value, unsigned char loc, T *sPtr, int ijk[3],
Will Schroeder's avatar
Will Schroeder committed
353
               unsigned char *edgeUses, vtkIdType *eIds)
354 355
{
  // Create a slightly faster path for pixel axes interior to the image.
356
  int ijk1[3];
357
  if ( edgeUses[0] ) //x axes edge
358
  {
359 360 361
    ijk1[0] = ijk[0] + 1;
    ijk1[1] = ijk[1];
    this->InterpolateAxesEdge(value, sPtr, ijk, sPtr+this->Inc0, ijk1, eIds[0]);
362
  }
363
  if ( edgeUses[2] ) //y axes edge
364
  {
365 366 367
    ijk1[0] = ijk[0];
    ijk1[1] = ijk[1] + 1;
    this->InterpolateAxesEdge(value, sPtr, ijk, sPtr+this->Inc1, ijk1, eIds[2]);
368
  }
369 370 371 372 373 374

  // Otherwise do more general gyrations. These are boundary situations where
  // the pixel axes is not fully formed. These situations occur on the
  // +x,+y image boundaries. (The other cases are handled by the default:
  // case and are expected.)
  switch (loc)
375
  {
376
    case 2: //+x edge
377
      this->InterpolateEdge(value, sPtr, ijk, 3, edgeUses, eIds);
378 379 380
      break;

    case 8: //+y
381
      this->InterpolateEdge(value, sPtr, ijk, 1, edgeUses, eIds);
382 383 384
      break;

    case 10: //+x +y
385 386
      this->InterpolateEdge(value, sPtr, ijk, 1, edgeUses, eIds);
      this->InterpolateEdge(value, sPtr, ijk, 3, edgeUses, eIds);
387 388 389 390
      break;

    default: //interior, or -x,-y boundary
      return;
391
  }
392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
}

//----------------------------------------------------------------------------
// PASS 1: Process a single x-row (and all of the pixel edges that compose
// the row).  Start building cell contour case table, determine the number of
// intersections, figure out where intersections along row begin and end
// (computational trimming).
template <class T> void vtkFlyingEdges2DAlgorithm<T>::
ProcessXEdge(double value, T* inPtr, vtkIdType row)
{
  vtkIdType nxcells=this->Dims[0]-1;
  vtkIdType minInt=nxcells, maxInt = 0;
  unsigned char edgeCase;
  vtkIdType *eMD=this->EdgeMetaData+row*5;
  unsigned char *ePtr=this->XCases + row*nxcells;
  double s0, s1 = static_cast<double>(*inPtr);

  //run along the entire x-edge computing edge cases
  std::fill_n(eMD, 5, 0);
  for (vtkIdType i=0; i < nxcells; ++i, ++ePtr)
412
  {
413 414 415 416 417 418 419 420 421 422 423 424 425
    s0 = s1;
    s1 = static_cast<double>(*(inPtr + (i+1)*this->Inc0));

    edgeCase  = ( s0 < value ? vtkFlyingEdges2DAlgorithm::Below :
                  vtkFlyingEdges2DAlgorithm::LeftAbove );
    edgeCase |= ( s1 < value ? vtkFlyingEdges2DAlgorithm::Below :
                  vtkFlyingEdges2DAlgorithm::RightAbove );

    this->SetXEdge(ePtr, edgeCase);

    // if edge intersects contour
    if ( edgeCase == vtkFlyingEdges2DAlgorithm::LeftAbove ||
         edgeCase == vtkFlyingEdges2DAlgorithm::RightAbove )
426
    {
427 428 429
      eMD[0]++; //increment number of intersections along x-edge
      minInt = ( i < minInt ? i : minInt);
      maxInt = i + 1;
430 431
    }//if contour interacts with this x-edge
  }//for all x-cell edges along this x-edge
432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447

  // The beginning and ending of intersections along the edge is used for
  // computational trimming.
  eMD[3] = minInt; //where intersections start along x edge
  eMD[4] = maxInt; //where intersections end along x edge
}

//----------------------------------------------------------------------------
// PASS 2: Process the y-cell edges (that form the cell axes) along a single
// x-row.  Continue building cell contour case table, and determine the
// number of cell y-edge intersections. Use computational trimming to reduce
// work.
template <class T> void vtkFlyingEdges2DAlgorithm<T>::
ProcessYEdges(vtkIdType row)
{
  // Grab the two edge cases bounding this pixel x-row.
Will Schroeder's avatar
Will Schroeder committed
448
  unsigned char *ePtr0, *ePtr1, ec0, ec1, xInts=1;
449 450 451 452 453 454 455 456 457 458 459
  ePtr0 = this->XCases + row*(this->Dims[0]-1);
  ePtr1 = ePtr0 + this->Dims[0]-1;

  // And metadata
  vtkIdType *eMD0 = this->EdgeMetaData + row*5;
  vtkIdType *eMD1 = this->EdgeMetaData + (row+1)*5;

  // Determine whether this row of x-cells needs processing. If there are no
  // x-edge intersections, and there is no y-edge (any y-edge) then the
  // row is contour free.
  if ( (eMD0[0] | eMD1[0]) == 0 ) //any x-ints?
460
  {
461
    if ( *ePtr0 == *ePtr1 )
462
    {
Will Schroeder's avatar
Will Schroeder committed
463
      return; //there are no x- or y-ints, thus no contour, skip pixel row
464
    }
Will Schroeder's avatar
Will Schroeder committed
465
    else
466
    {
Will Schroeder's avatar
Will Schroeder committed
467
      xInts = 0; //there are y- edge ints however
468
    }
469
  }
470 471 472 473 474 475 476 477 478 479

  // Determine proximity to the boundary of the image. This information is used
  // to count edge intersections in boundary situations.
  unsigned char loc, yLoc;
  yLoc = ( (row >= (this->Dims[1]-2) ? MaxBoundary : Interior) << 2);

  // The trim y-edges may need adjustment if the contour travels between
  // the top and bottom rows of x-edges (without intersecting x-edges).
  vtkIdType xL = ( (eMD0[3] < eMD1[3]) ? eMD0[3] : eMD1[3]);
  vtkIdType xR = ( (eMD0[4] > eMD1[4]) ? eMD0[4] : eMD1[4]);
Will Schroeder's avatar
Will Schroeder committed
480
  if ( xInts )
481
  {
Will Schroeder's avatar
Will Schroeder committed
482
    if ( xL > 0 )
483
    {
Will Schroeder's avatar
Will Schroeder committed
484 485 486
      ec0 = *(ePtr0 + xL);
      ec1 = *(ePtr1 + xL);
      if ( (ec0 & 0x1) != (ec1 & 0x1) )
487
      {
Will Schroeder's avatar
Will Schroeder committed
488
        xL = eMD0[3] = 0; //reset left trim
489
      }
490
    }
Will Schroeder's avatar
Will Schroeder committed
491
    if ( xR < (this->Dims[0]-1) )
492
    {
Will Schroeder's avatar
Will Schroeder committed
493 494 495
      ec0 = *(ePtr0 + xR);
      ec1 = *(ePtr1 + xR);
      if ( (ec0 & 0x2) != (ec1 & 0x2) )
496
      {
Will Schroeder's avatar
Will Schroeder committed
497
        xR = eMD0[4] = this->Dims[0] - 1; //reset right trim
498 499
      }
    }
500
  }
Will Schroeder's avatar
Will Schroeder committed
501
  else //contour cuts through without intersecting x-edges, reset trim edges
502
  {
Will Schroeder's avatar
Will Schroeder committed
503 504
    xL = eMD0[3] = 0;
    xR = eMD0[4] = this->Dims[0]-1;
505
  }
506 507 508 509 510 511 512 513

  // Okay run along the x-pixels and count the number of
  // y-intersections. Here we are just checking y edges that make up the
  // pixel axes. Also check the number of primitives generated.
  unsigned char *edgeUses, eCase, numLines;
  ePtr0 += xL; ePtr1 += xL;
  vtkIdType i;
  for (i=xL; i < xR; ++i) //run along the trimmed x-pixels
514
  {
515 516
    eCase = this->GetEdgeCase(ePtr0,ePtr1);
    if ( (numLines=this->GetNumberOfPrimitives(eCase)) > 0 )
517
    {
518 519 520 521 522 523 524 525 526 527
      // Okay let's increment the triangle count.
      eMD0[2] += numLines;

      // Count the number of y-points to be generated. Pass# 1 counted
      // the number of x-intersections along the x-edges. Now we count all
      // intersections on the y-pixel axes.
      edgeUses = this->GetEdgeUses(eCase);
      eMD0[1] += edgeUses[2]; //y-pixel axes edge always counted
      loc = yLoc | (i >= (this->Dims[0]-2) ? MaxBoundary : Interior);
      if ( loc != 0 )
528
      {
529
        this->CountBoundaryYInts(loc,edgeUses,eMD0);
530 531
      }
    }//if cell contains contour
532 533 534

    // advance the two pointers along pixel row
    ePtr0++; ePtr1++;
535
  }//for all pixels along this x-edge
536 537 538
}

//----------------------------------------------------------------------------
539 540
// PASS 4: Process the x-row cells to generate output primitives, including
// point coordinates and line segments. This is the fourth pass of the
541 542 543 544 545 546 547 548
// algorithm.
template <class T> void vtkFlyingEdges2DAlgorithm<T>::
GenerateOutput(double value, T* rowPtr, vtkIdType row)
{
  vtkIdType *eMD0 = this->EdgeMetaData + row*5;
  vtkIdType *eMD1 = this->EdgeMetaData + (row+1)*5;
  // Return if there is nothing to do (i.e., no lines to generate)
  if ( eMD0[2] == eMD1[2] )
549
  {
550
    return;
551
  }
552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579

  //Get the trim edges and prepare to generate
  vtkIdType i;
  vtkIdType xL = ( (eMD0[3] < eMD1[3]) ? eMD0[3] : eMD1[3]);
  vtkIdType xR = ( (eMD0[4] > eMD1[4]) ? eMD0[4] : eMD1[4]);

  // Grab the two edge cases bounding this pixel x-row. Begin at left trim edge.
  unsigned char *ePtr0, *ePtr1;
  ePtr0 = this->XCases + row*(this->Dims[0]-1) + xL;
  ePtr1 = ePtr0 + this->Dims[0]-1;

  // Traverse all pixels in this row, those containing the contour are
  // further identified for processing, meaning generating points and
  // triangles. Begin by setting up point ids on pixel edges.
  vtkIdType lineId = eMD0[2];
  vtkIdType eIds[4]; //the ids of generated points
  unsigned char *edgeUses, numLines;
  unsigned char eCase = this->InitPixelIds(ePtr0,ePtr1,eMD0,eMD1,eIds);

  // Determine the proximity to the boundary of volume. This information is
  // used to generate edge intersections.
  unsigned char loc, yLoc;
  yLoc = ( (row >= (this->Dims[1]-2) ? MaxBoundary : Interior) << 2);

  // Run along pixels in x-row direction and generate output primitives. Note
  // that active pixel axes edges are interpolated to produce points and
  // possibly interpolate attribute data.
  T *sPtr;
580 581 582
  int ijk[3];
  ijk[1] = row;
  ijk[2] = this->K;
583
  for (i=xL; i < xR; ++i)
584
  {
585
    if ( (numLines=this->GetNumberOfPrimitives(eCase)) > 0 )
586
    {
587 588 589 590 591 592 593
      // Start by generating triangles for this case
      this->GenerateLines(eCase,numLines,eIds,lineId);

      // Now generate point(s) along pixel axes if needed. Remember to take
      // boundary into account.
      loc = yLoc | (i >= (this->Dims[0]-2) ? MaxBoundary : Interior);
      if ( this->CaseIncludesAxes(eCase) || loc != Interior )
594
      {
595
        sPtr = rowPtr + i*this->Inc0;
596
        ijk[0] = i;
597
        edgeUses = this->GetEdgeUses(eCase);
598
        this->GeneratePoints(value, loc, sPtr, ijk, edgeUses, eIds);
599
      }
600 601

      this->AdvancePixelIds(eCase,eIds);
602
    }
603 604 605 606

    // advance along pixel row
    ePtr0++; ePtr1++;
    eCase = GetEdgeCase(ePtr0,ePtr1);
607
  } //for all non-trimmed cells along this x-edge
608 609 610 611 612 613 614
}

//----------------------------------------------------------------------------
// Contouring filter specialized for images. This templated function interfaces the
// vtkFlyingEdges2D class with the templated algorithm class. It also invokes
// the three passes of the Flying Edges algorithm.
//
615 616 617 618
template <class T> void vtkFlyingEdges2DAlgorithm<T>::
ContourImage(vtkFlyingEdges2D *self, T *scalars, vtkPoints *newPts,
             vtkDataArray *newScalars, vtkCellArray *newLines,
             vtkImageData *input, int *updateExt)
619 620 621 622 623 624 625 626 627 628 629
{
  double value, *values = self->GetValues();
  int numContours = self->GetNumberOfContours();
  vtkIdType vidx, row, *eMD;
  vtkIdType numOutXPts, numOutYPts, numOutLines, numXPts=0, numYPts=0, numLines=0;
  vtkIdType startXPts, startYPts, startLines; startXPts=startYPts=startLines=0;

  // The update extent may be different than the extent of the image.
  // The only problem with using the update extent is that one or two
  // sources enlarge the update extent.  This behavior is slated to be
  // eliminated.
630 631
  vtkIdType incs[3];
  input->GetIncrements(incs);
632 633 634 635 636 637
  int *ext = input->GetExtent();

  // Figure out which 2D plane the image lies in. Capture information for
  // subsequent processing.
  vtkFlyingEdges2DAlgorithm<T> algo;
  if (updateExt[4] == updateExt[5])
638
  { // z collapsed
639 640 641 642 643 644 645 646
    algo.Axis0 = 0;
    algo.Min0 = updateExt[0];
    algo.Max0 = updateExt[1];
    algo.Inc0 = incs[0];
    algo.Axis1 = 1;
    algo.Min1 = updateExt[2];
    algo.Max1 = updateExt[3];
    algo.Inc1 = incs[1];
647
    algo.K = updateExt[4];
648
    algo.Axis2 = 2;
649
  }
650
  else if (updateExt[2] == updateExt[3])
651
  { // y collapsed
652 653 654 655 656 657 658 659
    algo.Axis0 = 0;
    algo.Min0 = updateExt[0];
    algo.Max0 = updateExt[1];
    algo.Inc0 = incs[0];
    algo.Axis1 = 2;
    algo.Min1 = updateExt[4];
    algo.Max1 = updateExt[5];
    algo.Inc1 = incs[2];
660
    algo.K = updateExt[2];
661
    algo.Axis2 = 1;
662
  }
663
  else if (updateExt[0] == updateExt[1])
664
  { // x collapsed
665 666 667 668 669 670 671 672
    algo.Axis0 = 1;
    algo.Min0 = updateExt[2];
    algo.Max0 = updateExt[3];
    algo.Inc0 = incs[1];
    algo.Axis1 = 2;
    algo.Min1 = updateExt[4];
    algo.Max1 = updateExt[5];
    algo.Inc1 = incs[2];
673
    algo.K = updateExt[0];
674
    algo.Axis2 = 0;
675
  }
676
  else
677
  {
678 679
    vtkGenericWarningMacro("Expecting 2D data.");
    return;
680
  }
681 682 683 684 685 686 687 688 689 690 691 692 693 694 695

  // Now allocate working arrays. The XCases array tracks case# for each cell.
  algo.Dims[0] = algo.Max0 - algo.Min0 + 1;
  algo.Dims[1] = algo.Max1 - algo.Min1 + 1;
  algo.XCases = new unsigned char [(algo.Dims[0]-1)*algo.Dims[1]];

  // Also allocate the characterization (metadata) array for the x edges.
  // This array tracks the number of intersections along each x-row, y-row;
  // as well as num line primitives, and the xMin_i and xMax_i (minimum
  // index of first intersection, maximum index of intersection for row i,
  // so-called trim edges used for computational trimming).
  algo.EdgeMetaData = new vtkIdType [algo.Dims[1]*5];

  // Compute the starting location for scalar data.  We may be operating
  // on a part of the image.
696 697 698
  algo.Scalars = scalars + incs[0]*(updateExt[0]-ext[0]) +
    incs[1]*(updateExt[2]-ext[2]) + incs[2]*(updateExt[4]-ext[4]) +
    self->GetArrayComponent();
699 700 701

  // The algorithm is separated into multiple passes. The first pass
  // computes intersections on row edges, counting the number of intersected edges
702
  // as it progresses. It also keeps track of the generated edge cases and
703 704 705 706 707 708
  // other incidental information about intersections along rows. The second
  // pass generates polylines from the cases and intersection information.
  // In the final and third pass output points and lines are generated.

  // Loop across each contour value. This encompasses all three passes.
  for (vidx = 0; vidx < numContours; vidx++)
709
  {
710 711 712 713 714
    value = values[vidx];

    // PASS 1: Traverse all rows generating intersection points and building
    // the case table. Also accumulate information necessary for later allocation.
    // For example the number of output points is computed.
715
    Pass1<T> pass1(&algo,value);
716
    vtkSMPTools::For(0,algo.Dims[1], pass1);
717 718 719 720

    // PASS 2: Traverse all rows and process cell y edges. Continue building
    // case table from y contributions (using computational trimming to reduce
    // work) and keep track of cell y intersections.
721
    Pass2<T> pass2(&algo);
722
    vtkSMPTools::For(0,algo.Dims[1]-1, pass2);
723 724 725 726 727 728 729 730 731 732 733

    // PASS 3: Now allocate and generate output. First we have to update the
    // x-Edge meta data to partition the output into separate pieces so
    // independent threads can write into separate memory partititions. Once
    // allocation is complete, process on a row by row basis and produce
    // output points, line primitives, and interpolate point attribute data
    // (if necessary).
    numOutXPts = startXPts;
    numOutYPts = startYPts;
    numOutLines = startLines;
    for (row=0; row < algo.Dims[1]; ++row)
734
    {
735 736 737 738 739 740 741 742 743 744
      eMD = algo.EdgeMetaData + row*5;
      numXPts = eMD[0];
      numYPts = eMD[1];
      numLines = eMD[2];
      eMD[0] = numOutXPts + numOutYPts;
      eMD[1] = eMD[0] + numXPts;
      eMD[2] = numOutLines;
      numOutXPts += numXPts;
      numOutYPts += numYPts;
      numOutLines += numLines;
745
    }
746 747

    // Output can now be allocated.
Will Schroeder's avatar
Will Schroeder committed
748 749
    vtkIdType totalPts = numOutXPts + numOutYPts;
    if ( totalPts > 0 )
750
    {
Will Schroeder's avatar
Will Schroeder committed
751 752 753 754 755
      newPts->GetData()->WriteVoidPointer(0,3*totalPts);
      algo.NewPoints = static_cast<float*>(newPts->GetVoidPointer(0));
      newLines->WritePointer(numOutLines,3*numOutLines);
      algo.NewLines = static_cast<vtkIdType*>(newLines->GetPointer());
      if (newScalars)
756
      {
757 758 759
        vtkIdType numPrevPts = newScalars->GetNumberOfTuples();
        vtkIdType numNewPts = totalPts - numPrevPts;
        newScalars->WriteVoidPointer(0,totalPts);
Will Schroeder's avatar
Will Schroeder committed
760 761
        algo.NewScalars = static_cast<T*>(newScalars->GetVoidPointer(0));
        T TValue = static_cast<T>(value);
762
        std::fill_n(algo.NewScalars+numPrevPts, numNewPts, TValue);
763
      }
764

Will Schroeder's avatar
Will Schroeder committed
765 766 767
      // PASS 4: Now process each x-row and produce the output primitives.
      Pass4<T> pass4(&algo,value);
      vtkSMPTools::For(0,algo.Dims[1]-1, pass4);
768
    }//if output generated
769 770 771 772 773

    // Handle multiple contours
    startXPts = numOutXPts;
    startYPts = numOutYPts;
    startLines = numOutLines;
774
  }// for all contour values
775 776 777 778 779 780

  // Clean up and return
  delete [] algo.XCases;
  delete [] algo.EdgeMetaData;
}

781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804
//----------------------------------------------------------------------------
// Here is the VTK class proper.
// Construct object with initial contour value of 0.0.
vtkFlyingEdges2D::vtkFlyingEdges2D()
{
  this->ContourValues = vtkContourValues::New();
  this->ComputeScalars = 1;
  this->ArrayComponent = 0;

  // by default process active point scalars
  this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
                               vtkDataSetAttributes::SCALARS);
}

//----------------------------------------------------------------------------
vtkFlyingEdges2D::~vtkFlyingEdges2D()
{
  this->ContourValues->Delete();
}

//----------------------------------------------------------------------------
// Description:
// Overload standard modified time function. If contour values are modified,
// then this object is modified as well.
Bill Lorensen's avatar
Bill Lorensen committed
805
vtkMTimeType vtkFlyingEdges2D::GetMTime()
806
{
Bill Lorensen's avatar
Bill Lorensen committed
807 808
  vtkMTimeType mTime=this->Superclass::GetMTime();
  vtkMTimeType mTime2=this->ContourValues->GetMTime();
809 810 811 812

  return ( mTime2 > mTime ? mTime2 : mTime );
}

813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835
//----------------------------------------------------------------------------
// Contouring filter specialized for images (or slices from images)
//
int vtkFlyingEdges2D::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 output
  vtkImageData *input = vtkImageData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkDataArray  *inScalars;

  vtkDebugMacro(<< "Executing 2D Flying Edges");

  int *ext =
    inInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT());
  inScalars = this->GetInputArrayToProcess(0,inputVector);
836
  if ( inScalars == nullptr )
837
  {
838 839
    vtkErrorMacro(<<"Scalars must be defined for contouring");
    return 1;
840
  }
841 842 843

  int numComps = inScalars->GetNumberOfComponents();
  if (this->ArrayComponent >= numComps)
844
  {
845 846 847
    vtkErrorMacro("Scalars have " << numComps << " components. "
                  "ArrayComponent must be smaller than " << numComps);
    return 1;
848
  }
849 850 851 852 853 854

  // Create necessary objects to hold output. We will defer the
  // actual allocation to a later point.
  vtkCellArray *newLines = vtkCellArray::New();
  vtkPoints *newPts = vtkPoints::New();
  newPts->SetDataTypeToFloat();
855
  vtkDataArray *newScalars = nullptr;
856 857

  if (this->ComputeScalars)
858
  {
859 860 861
    newScalars = inScalars->NewInstance();
    newScalars->SetNumberOfComponents(1);
    newScalars->SetName(inScalars->GetName());
862
  }
863 864 865 866

  // Check data type and execute appropriate function
  void *scalars = inScalars->GetVoidPointer(0);
  switch (inScalars->GetDataType())
867
  {
868 869
    vtkTemplateMacro(vtkFlyingEdges2DAlgorithm<VTK_TT>::
                     ContourImage(this,(VTK_TT *)scalars, newPts,
870
                      newScalars, newLines, input, ext));
871
  }//switch
872 873 874 875 876 877 878 879 880 881 882 883 884 885

  vtkDebugMacro(<<"Created: "
                << newPts->GetNumberOfPoints() << " points, "
                << newLines->GetNumberOfCells() << " lines");

  // Update ourselves.  Because we don't know up front how many lines
  // we've created, take care to reclaim memory.
  output->SetPoints(newPts);
  newPts->Delete();

  output->SetLines(newLines);
  newLines->Delete();

  if (newScalars)
886
  {
887 888 889
    int idx = output->GetPointData()->AddArray(newScalars);
    output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
    newScalars->Delete();
890
  }
891

892 893
  vtkImageTransform::TransformPointSet(input, output);

894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913
  return 1;
}

//----------------------------------------------------------------------------
int vtkFlyingEdges2D::FillInputPortInformation(int, vtkInformation *info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkImageData");
  return 1;
}

//----------------------------------------------------------------------------
void vtkFlyingEdges2D::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

  this->ContourValues->PrintSelf(os,indent.GetNextIndent());

  os << indent << "Compute Scalars: " << (this->ComputeScalars ? "On\n" : "Off\n");
  os << indent << "ArrayComponent: " << this->ArrayComponent << endl;
}