Method Contour

Common/DataModel/vtkPixel.cxx:310–367

static int edges[4][2] = { { 0, 1 }, { 1, 3 }, { 2, 3 }, { 0, 2 } };

void vtkPixel::Contour(double value, vtkDataArray* cellScalars, vtkIncrementalPointLocator* locator,
  vtkCellArray* vtkNotUsed(verts), vtkCellArray* lines, vtkCellArray* vtkNotUsed(polys),
  vtkPointData* inPd, vtkPointData* outPd, vtkCellData* inCd, vtkIdType cellId, vtkCellData* outCd)
{
  static const int CASE_MASK[4] = { 1, 2, 8, 4 }; // note differenceom quad!
  vtkMarchingSquaresLineCases* lineCase;
  int* edge;
  int i, j, index, *vert;
  int newCellId;
  vtkIdType pts[2];
  double t, x1[3], x2[3], x[3];

  // Build the case table
  for (i = 0, index = 0; i < 4; i++)
  {
    if (cellScalars->GetComponent(i, 0) >= value)
    {
      index |= CASE_MASK[i];
    }
  }

  lineCase = vtkMarchingSquaresLineCases::GetCases() + index;
  edge = lineCase->edges;

  for (; edge[0] > -1; edge += 2)
  {
    for (i = 0; i < 2; i++) // insert line
    {
      vert = edges[edge[i]];
      t = (value - cellScalars->GetComponent(vert[0], 0)) /
        (cellScalars->GetComponent(vert[1], 0) - cellScalars->GetComponent(vert[0], 0));
      this->Points->GetPoint(vert[0], x1);
      this->Points->GetPoint(vert[1], x2);
      for (j = 0; j < 3; j++)
      {
        x[j] = x1[j] + t * (x2[j] - x1[j]);
      }
      if (locator->InsertUniquePoint(x, pts[i]))
      {
        if (outPd)
        {
          int p1 = this->PointIds->GetId(vert[0]);
          int p2 = this->PointIds->GetId(vert[1]);
          outPd->InterpolateEdge(inPd, pts[i], p1, p2, t);
        }
      }
    }
    // check for degenerate line
    if (pts[0] != pts[1])
    {
      newCellId = lines->InsertNextCell(2, pts);
      if (outCd)
      {
        outCd->CopyData(inCd, cellId, newCellId);
      }
    }
  }
}