/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkFlashContour.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 "vtkFlashContour.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkDoubleArray.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMarchingCubesTriangleCases.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkMultiPieceDataSet.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkUnsignedCharArray.h"

vtkStandardNewMacro(vtkFlashContour);

// How do we find edge/corner neighbors and neighbors in different levels.
// We could keep neighbors of global blocks (even ones not loaded),
// and save the global to local map.
// !!!!!  I added the ability to interpolate one cell array from the input
// to a point array of the output.  The problems with generalizing this is that
// for shared regions we will have to interpolated between two blocks.
// We cannot use InterpolatePoint (whatever it is called in vtkDataArray).
// I would have to write a similar method that takes two FieldDatas as input.

static int vtkFlashIsoEdgeToPointsTable[12][2] = { { 0, 1 }, { 1, 3 }, { 2, 3 }, { 0, 2 }, { 4, 5 },
  { 5, 7 }, { 6, 7 }, { 4, 6 }, { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 } };
static int vtkFlashIsoEdgeToVTKPointsTable[12][2] = { { 0, 1 }, { 1, 2 }, { 3, 2 }, { 0, 3 },
  { 4, 5 }, { 5, 6 }, { 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 3, 7 }, { 2, 6 } };

//============================================================================
//----------------------------------------------------------------------------
// Description:
// Construct object with initial range (0,1) and single contour value
// of 0.0. ComputeNormal is on, ComputeGradients is off and ComputeScalars is on.
vtkFlashContour::vtkFlashContour()
{
  this->IsoValue = 100.0;
  this->PassAttribute = 0;
  this->PassArray = 0;
  this->CellArrayNameToProcess = 0;

  // Pipeline
  this->SetNumberOfOutputPorts(1);
}

//----------------------------------------------------------------------------
vtkFlashContour::~vtkFlashContour()
{
  this->SetCellArrayNameToProcess(0);
  this->SetPassAttribute(0);
}

//----------------------------------------------------------------------------
void vtkFlashContour::PrintSelf(ostream& os, vtkIndent indent)
{
  // TODO print state
  this->Superclass::PrintSelf(os, indent);

  os << indent << "IsoValue: " << this->IsoValue << endl;
  if (this->PassAttribute)
  {
    os << indent << "PassAttribute: " << this->PassAttribute << endl;
  }
}

//----------------------------------------------------------------------------
int vtkFlashContour::FillInputPortInformation(int /*port*/, vtkInformation* info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");

  return 1;
}

//----------------------------------------------------------------------------
int vtkFlashContour::FillOutputPortInformation(int port, vtkInformation* info)
{
  switch (port)
  {
    case 0:
      info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkMultiBlockDataSet");
      break;
    default:
      vtkErrorMacro("Invalid output port.");
      break;
  }

  return 1;
}

// The placement of neighbors in the field-data global neighbor array
// generated by the flash reader.
#define XMIN 0
#define XMAX 1
#define YMIN 2
#define YMAX 3
#define ZMIN 4
#define ZMAX 5

//----------------------------------------------------------------------------
void vtkFlashContour::PropogateNeighbors(int neighbors[3][3][3], int x, int y, int z)
{
  if (neighbors[x][y][z] != -1)
  {
    return;
  }

  // It would be nice to loop over faces here.
  int nx, ny, nz;
  nx = x - 1;
  ny = y;
  nz = z;
  if (nx >= 0)
  {
    int blockId = neighbors[nx][ny][nz];
    if (blockId >= 0)
    {
      blockId = this->GlobalNeighborArray[blockId * 6 + XMAX];
      if (blockId >= 0)
      {
        neighbors[x][y][z] = blockId;
        return;
      }
    }
  }
  nx = x + 1;
  ny = y;
  nz = z;
  if (nx <= 2)
  {
    int blockId = neighbors[nx][ny][nz];
    if (blockId >= 0)
    {
      blockId = this->GlobalNeighborArray[blockId * 6 + XMIN];
      if (blockId >= 0)
      {
        neighbors[x][y][z] = blockId;
        return;
      }
    }
  }
  nx = 1;
  ny = y - 1;
  nz = z;
  if (ny >= 0)
  {
    int blockId = neighbors[nx][ny][nz];
    if (blockId >= 0)
    {
      blockId = this->GlobalNeighborArray[blockId * 6 + YMAX];
      if (blockId >= 0)
      {
        neighbors[x][y][z] = blockId;
        return;
      }
    }
  }
  nx = x;
  ny = y + 1;
  nz = z;
  if (ny <= 2)
  {
    int blockId = neighbors[nx][ny][nz];
    if (blockId >= 0)
    {
      blockId = this->GlobalNeighborArray[blockId * 6 + YMIN];
      if (blockId >= 0)
      {
        neighbors[x][y][z] = blockId;
        return;
      }
    }
  }
  nx = x;
  ny = y;
  nz = z - 1;
  if (nz >= 0)
  {
    int blockId = neighbors[nx][ny][nz];
    if (blockId >= 0)
    {
      blockId = this->GlobalNeighborArray[blockId * 6 + ZMAX];
      if (blockId >= 0)
      {
        neighbors[x][y][z] = blockId;
        return;
      }
    }
  }
  nx = x;
  ny = y;
  nz = z + 1;
  if (nx <= 2)
  {
    int blockId = neighbors[nx][ny][nz];
    if (blockId >= 0)
    {
      blockId = this->GlobalNeighborArray[blockId * 6 + ZMIN];
      if (blockId >= 0)
      {
        neighbors[x][y][z] = blockId;
        return;
      }
    }
  }
}

//----------------------------------------------------------------------------
int vtkFlashContour::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the data set which we are to process
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkMultiBlockDataSet* mbdsInput =
    vtkMultiBlockDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkDataArray* da = mbdsInput->GetFieldData()->GetArray("GlobalToLocalMap");
  vtkIntArray* globalToLocalMapArray = vtkIntArray::SafeDownCast(da);
  if (globalToLocalMapArray == 0)
  {
    vtkErrorMacro("Missing block map.");
    return 0;
  }
  this->GlobalToLocalMap = (int*)(globalToLocalMapArray->GetVoidPointer(0));
  this->NumberOfGlobalBlocks = globalToLocalMapArray->GetNumberOfTuples();

  da = mbdsInput->GetFieldData()->GetArray("BlockChildren");
  vtkIntArray* childrenIdArray = vtkIntArray::SafeDownCast(da);
  if (childrenIdArray == 0)
  {
    vtkErrorMacro("Missing children array.");
    return 0;
  }
  this->GlobalChildrenArray = (int*)(childrenIdArray->GetVoidPointer(0));

  da = mbdsInput->GetFieldData()->GetArray("BlockNeighbors");
  vtkIntArray* neighborIdArray = vtkIntArray::SafeDownCast(da);
  if (neighborIdArray == 0)
  {
    vtkErrorMacro("Missing children array.");
    return 0;
  }
  this->GlobalNeighborArray = (int*)(neighborIdArray->GetVoidPointer(0));

  da = mbdsInput->GetFieldData()->GetArray("BlockLevel");
  vtkIntArray* levelArray = vtkIntArray::SafeDownCast(da);
  if (levelArray == 0)
  {
    vtkErrorMacro("Missing level array.");
    return 0;
  }
  this->GlobalLevelArray = (int*)(levelArray->GetVoidPointer(0));

  // Get the outputs
  // 0
  vtkInformation* outInfo;
  outInfo = outputVector->GetInformationObject(0);
  vtkMultiBlockDataSet* mbdsOutput0 =
    vtkMultiBlockDataSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  mbdsOutput0->SetNumberOfBlocks(1);
  vtkMultiPieceDataSet* mpds = vtkMultiPieceDataSet::New();
  mbdsOutput0->SetBlock(0, mpds);

  mpds->SetNumberOfPieces(0);

  if (mbdsInput == 0)
  {
    // Do not deal with rectilinear grid
    vtkErrorMacro("This filter requires a vtkMultiBlockDataSet on its input.");
    return 0;
  }

  // This is alot to go through to get the name of the array to process.
  vtkInformationVector* inArrayVec = this->GetInformation()->Get(INPUT_ARRAYS_TO_PROCESS());
  if (!inArrayVec)
  {
    vtkErrorMacro("Problem finding array to process");
    return 0;
  }
  vtkInformation* inArrayInfo = inArrayVec->GetInformationObject(0);
  if (!inArrayInfo)
  {
    vtkErrorMacro("Problem getting name of array to process.");
    return 0;
  }
  if (!inArrayInfo->Has(vtkDataObject::FIELD_NAME()))
  {
    vtkErrorMacro("Missing field name.");
    return 0;
  }
  const char* arrayNameToProcess = inArrayInfo->Get(vtkDataObject::FIELD_NAME());
  this->SetCellArrayNameToProcess(arrayNameToProcess);

  this->Mesh = vtkPolyData::New();
  this->Points = vtkPoints::New();
  this->Faces = vtkCellArray::New();
  this->Mesh->SetPoints(this->Points);
  this->Mesh->SetPolys(this->Faces);
  mpds->SetPiece(0, this->Mesh);

  this->BlockIdCellArray = vtkIntArray::New();
  this->BlockIdCellArray->SetName("GlobalBlockId");
  this->LevelCellArray = vtkUnsignedCharArray::New();
  this->LevelCellArray->SetName("Level");
  this->RemainingDepthCellArray = vtkUnsignedCharArray::New();
  this->RemainingDepthCellArray->SetName("HiddenLevels");
  this->Mesh->GetCellData()->AddArray(this->BlockIdCellArray);
  this->Mesh->GetCellData()->AddArray(this->LevelCellArray);
  this->Mesh->GetCellData()->AddArray(this->RemainingDepthCellArray);

  // For debugging.
  // this->BlockIdCellArray = vtkIntArray::New();
  // this->BlockIdCellArray->SetName("BlockIds");
  // mesh->GetCellData()->AddArray(this->BlockIdCellArray);

  this->PassArray = 0;
  if (this->PassAttribute && mbdsInput->GetNumberOfBlocks() > 0)
  {
    // Find the array we are supposed to pass.
    // pain empty blocks.
    int blockId = 0;
    int numBlocks = mbdsInput->GetNumberOfBlocks();
    vtkImageData* block = 0;
    while (block == 0 && blockId < numBlocks)
    {
      block = vtkImageData::SafeDownCast(mbdsInput->GetBlock(blockId));
      ++blockId;
    }
    if (block)
    {
      da = block->GetCellData()->GetArray(this->PassAttribute);
      if (da == 0 || da->GetDataType() != VTK_DOUBLE)
      {
        vtkErrorMacro("We can only pass double arrays.");
      }
      else
      {
        this->PassArray = vtkDoubleArray::New();
        this->PassArray->SetName(this->PassAttribute);
        this->Mesh->GetPointData()->AddArray(this->PassArray);
      }
    }
    else
    {
      vtkWarningMacro("Cannot find first block for attribute to pass.");
    }
  }

  // Find all roots and recurse on each.
  int* levelPtr = this->GlobalLevelArray;
  for (int i = 0; i < this->NumberOfGlobalBlocks; ++i)
  {
    if (*levelPtr++ == 1)
    { // Root node
      // Find neighbors
      // I am going to traverse the tree of blocks keeping the 26 neighbors in an array.
      int neighborhood[3][3][3];
      for (int x = 0; x < 3; ++x)
      {
        for (int y = 0; y < 3; ++y)
        {
          for (int z = 0; z < 3; ++z)
          {
            neighborhood[x][y][z] = -1;
          }
        }
      }
      neighborhood[1][1][1] = i;
      neighborhood[0][1][1] = this->GlobalNeighborArray[i * 6 + XMIN];
      neighborhood[2][1][1] = this->GlobalNeighborArray[i * 6 + XMAX];
      neighborhood[1][0][1] = this->GlobalNeighborArray[i * 6 + YMIN];
      neighborhood[1][2][1] = this->GlobalNeighborArray[i * 6 + YMAX];
      neighborhood[1][1][0] = this->GlobalNeighborArray[i * 6 + ZMIN];
      neighborhood[1][1][2] = this->GlobalNeighborArray[i * 6 + ZMAX];
      // Now we have to fill in edge and corner neighbors by traversing face neighbors.
      // It would be nice to generalize this but .....
      this->PropogateNeighbors(neighborhood, 0, 0, 1);
      this->PropogateNeighbors(neighborhood, 2, 0, 1);
      this->PropogateNeighbors(neighborhood, 0, 2, 1);
      this->PropogateNeighbors(neighborhood, 2, 2, 1);
      this->PropogateNeighbors(neighborhood, 0, 1, 0);
      this->PropogateNeighbors(neighborhood, 2, 1, 0);
      this->PropogateNeighbors(neighborhood, 0, 1, 2);
      this->PropogateNeighbors(neighborhood, 2, 1, 2);
      this->PropogateNeighbors(neighborhood, 2, 0, 1);
      this->PropogateNeighbors(neighborhood, 0, 2, 1);
      this->PropogateNeighbors(neighborhood, 2, 2, 1);
      // Corners have to be done after edges.
      this->PropogateNeighbors(neighborhood, 0, 0, 0);
      this->PropogateNeighbors(neighborhood, 2, 0, 0);
      this->PropogateNeighbors(neighborhood, 0, 2, 0);
      this->PropogateNeighbors(neighborhood, 2, 2, 0);
      this->PropogateNeighbors(neighborhood, 0, 0, 2);
      this->PropogateNeighbors(neighborhood, 2, 0, 2);
      this->PropogateNeighbors(neighborhood, 0, 2, 2);
      this->PropogateNeighbors(neighborhood, 2, 2, 2);

      this->RecurseTree(neighborhood, mbdsInput);
    }
  }

  this->Mesh->Delete();
  this->Points->Delete();
  this->Points = 0;
  this->Faces->Delete();
  this->Faces = 0;
  if (this->PassArray)
  {
    this->PassArray->Delete();
    this->PassArray = 0;
  }
  this->BlockIdCellArray->Delete();
  this->BlockIdCellArray = 0;
  ;
  this->LevelCellArray->Delete();
  this->LevelCellArray = 0;
  ;
  this->RemainingDepthCellArray->Delete();
  this->RemainingDepthCellArray = 0;

  mpds->Delete();

  return 1;
}

//----------------------------------------------------------------------------
unsigned char vtkFlashContour::ComputeBranchDepth(int globalBlockId)
{
  int* children = this->GlobalChildrenArray + (globalBlockId * 8);
  if (children[0] < 0)
  {
    return 0;
  }
  unsigned char max = 0;
  for (int i = 0; i < 8; ++i)
  {
    unsigned char tmp = this->ComputeBranchDepth(children[i]);
    if (tmp > max)
    {
      max = tmp;
    }
  }
  return max + 1;
}

//----------------------------------------------------------------------------
void vtkFlashContour::RecurseTree(int neighborhood[3][3][3], vtkMultiBlockDataSet* input)
{
  int parent = neighborhood[1][1][1];
  int* children = this->GlobalChildrenArray + (parent << 3);
  int* neighborChildren;
  // Check to see if the center of the neighborhood is a leaf.
  // Local map is set to -32 if node is below leaf.
  // Local map is set to -1 if node is not loaded but decendant is.
  // Assume all or none children. only check 0.
  if (children[0] > 0 && this->GlobalToLocalMap[children[0]] != -32)
  { // Not a leaf.  Traverse.
    for (int childIdx = 0; childIdx < 8; ++childIdx)
    {
      if (this->GlobalToLocalMap[children[childIdx]] == -32)
      {
        vtkErrorMacro("Partial node refinement.");
        continue;
      }
      int childNeighborhood[3][3][3];
      int neighbor;
      // This assumes children are indexed x, y, z.
      // Maybe we can find a symetry between cases and avoid a long case.
      // Names are taken from the childIdx == 0 case.
      int nx0, nx2, ny0, ny2, nz0, nz2;
      int cx0, cx1, cy0, cy1, cz0, cz1;
      cx0 = cy0 = cz0 = nx0 = ny0 = nz0 = 0;
      cx1 = cy1 = cz1 = 1;
      nx2 = ny2 = nz2 = 2;
      if (childIdx & 1)
      {
        nx0 = 2;
        nx2 = 0;
        cx0 = 1;
        cx1 = 0;
      }
      if (childIdx & 2)
      {
        ny0 = 2;
        ny2 = 0;
        cy0 = 1;
        cy1 = 0;
      }
      if (childIdx & 4)
      {
        nz0 = 2;
        nz2 = 0;
        cz0 = 1;
        cz1 = 0;
      }
      // Since children use a 1D 8 array rather than a 3D array like neighbors,
      // xyz children indexes need to be converted into bits.
      cy0 = cy0 << 1;
      cy1 = cy1 << 1;
      cz0 = cz0 << 2;
      cz1 = cz1 << 2;
      // Take care of siblings of chosen child.
      childNeighborhood[1][1][1] = children[childIdx];
      childNeighborhood[nx2][1][1] = children[cx1 | cy0 | cz0];     // (1) 1 0 0
      childNeighborhood[1][ny2][1] = children[cx0 | cy1 | cz0];     // (2) 0 1 0
      childNeighborhood[nx2][ny2][1] = children[cx1 | cy1 | cz0];   // (3) 1 1 0
      childNeighborhood[1][1][nz2] = children[cx0 | cy0 | cz1];     // (4) 0 0 1
      childNeighborhood[nx2][1][nz2] = children[cx1 | cy0 | cz1];   // (5) 1 0 1
      childNeighborhood[1][ny2][nz2] = children[cx0 | cy1 | cz1];   // (6) 0 1 1
      childNeighborhood[nx2][ny2][nz2] = children[cx1 | cy1 | cz1]; // (7) 1 1 1

      neighbor = neighborhood[nx0][1][1];
      neighborChildren = this->GlobalChildrenArray + (neighbor << 3);
      if (neighbor < 0 || neighborChildren[0] < 0 ||
        this->GlobalToLocalMap[neighborChildren[0]] == -32)
      {
        childNeighborhood[nx0][1][1] = neighbor;
        childNeighborhood[nx0][ny2][1] = neighbor;
        childNeighborhood[nx0][1][nz2] = neighbor;
        childNeighborhood[nx0][ny2][nz2] = neighbor;
      }
      else
      {
        childNeighborhood[nx0][1][1] = neighborChildren[cx1 | cy0 | cz0];     // (1) 1 0 0
        childNeighborhood[nx0][ny2][1] = neighborChildren[cx1 | cy1 | cz0];   // (3) 1 1 0
        childNeighborhood[nx0][1][nz2] = neighborChildren[cx1 | cy0 | cz1];   // (5) 1 0 1
        childNeighborhood[nx0][ny2][nz2] = neighborChildren[cx1 | cy1 | cz1]; // (7) 1 1 1
      }

      neighbor = neighborhood[1][ny0][1];
      neighborChildren = this->GlobalChildrenArray + (neighbor << 3);
      if (neighbor < 0 || neighborChildren[0] < 0 ||
        this->GlobalToLocalMap[neighborChildren[0]] == -32)
      {
        childNeighborhood[1][ny0][1] = neighbor;
        childNeighborhood[nx2][ny0][1] = neighbor;
        childNeighborhood[1][ny0][nz2] = neighbor;
        childNeighborhood[nx2][ny0][nz2] = neighbor;
      }
      else
      {
        childNeighborhood[1][ny0][1] = neighborChildren[cx0 | cy1 | cz0];     // (2) 0 1 0
        childNeighborhood[nx2][ny0][1] = neighborChildren[cx1 | cy1 | cz0];   // (3) 1 1 0
        childNeighborhood[1][ny0][nz2] = neighborChildren[cx0 | cy1 | cz1];   // (6) 0 1 1
        childNeighborhood[nx2][ny0][nz2] = neighborChildren[cx1 | cy1 | cz1]; // (7) 1 1 1
      }

      neighbor = neighborhood[1][1][nz0];
      neighborChildren = this->GlobalChildrenArray + (neighbor << 3);
      if (neighbor < 0 || neighborChildren[0] < 0 ||
        this->GlobalToLocalMap[neighborChildren[0]] == -32)
      {
        childNeighborhood[1][1][nz0] = neighbor;
        childNeighborhood[nx2][1][nz0] = neighbor;
        childNeighborhood[1][ny2][nz0] = neighbor;
        childNeighborhood[nx2][ny2][nz0] = neighbor;
      }
      else
      {
        childNeighborhood[1][1][nz0] = neighborChildren[cx0 | cy0 | cz1];     // (4) 0 0 1
        childNeighborhood[nx2][1][nz0] = neighborChildren[cx1 | cy0 | cz1];   // (5) 1 0 1
        childNeighborhood[1][ny2][nz0] = neighborChildren[cx0 | cy1 | cz1];   // (6) 0 1 1
        childNeighborhood[nx2][ny2][nz0] = neighborChildren[cx1 | cy1 | cz1]; // (7) 1 1 1
      }

      neighbor = neighborhood[nx0][ny0][1];
      neighborChildren = this->GlobalChildrenArray + (neighbor << 3);
      if (neighbor < 0 || neighborChildren[0] < 0 ||
        this->GlobalToLocalMap[neighborChildren[0]] == -32)
      {
        childNeighborhood[nx0][ny0][1] = neighbor;
        childNeighborhood[nx0][ny0][nz2] = neighbor;
      }
      else
      {
        childNeighborhood[nx0][ny0][1] = neighborChildren[cx1 | cy1 | cz0];   // (3) 1 1 0
        childNeighborhood[nx0][ny0][nz2] = neighborChildren[cx1 | cy1 | cz1]; // (7) 1 1 1
      }

      neighbor = neighborhood[nx0][1][nz0];
      neighborChildren = this->GlobalChildrenArray + (neighbor << 3);
      if (neighbor < 0 || neighborChildren[0] < 0 ||
        this->GlobalToLocalMap[neighborChildren[0]] == -32)
      {
        childNeighborhood[nx0][1][nz0] = neighbor;
        childNeighborhood[nx0][ny2][nz0] = neighbor;
      }
      else
      {
        childNeighborhood[nx0][1][nz0] = neighborChildren[cx1 | cy0 | cz1];   // (5) 1 0 1
        childNeighborhood[nx0][ny2][nz0] = neighborChildren[cx1 | cy1 | cz1]; // (7) 1 1 1
      }

      neighbor = neighborhood[1][ny0][nz0];
      neighborChildren = this->GlobalChildrenArray + (neighbor << 3);
      if (neighbor < 0 || neighborChildren[0] < 0 ||
        this->GlobalToLocalMap[neighborChildren[0]] == -32)
      {
        childNeighborhood[1][ny0][nz0] = neighbor;
        childNeighborhood[nx2][ny0][nz0] = neighbor;
      }
      else
      {
        childNeighborhood[1][ny0][nz0] = neighborChildren[cx0 | cy1 | cz1];   // (6) 0 1 1
        childNeighborhood[nx2][ny0][nz0] = neighborChildren[cx1 | cy1 | cz1]; // (7) 1 1 1
      }

      neighbor = neighborhood[nx0][ny0][nz0];
      neighborChildren = this->GlobalChildrenArray + (neighbor << 3);
      if (neighbor < 0 || neighborChildren[0] < 0 ||
        this->GlobalToLocalMap[neighborChildren[0]] == -32)
      {
        childNeighborhood[nx0][ny0][nz0] = neighbor;
      }
      else
      {
        childNeighborhood[nx0][ny0][nz0] = neighborChildren[cx1 | cy1 | cz1]; // (7) 1 1 1
      }
      this->RecurseTree(childNeighborhood, input);
    }
    return;
  }

  // Center of neighborhood is a leaf.
  // Contour the block and the shared regions it owns.
  int globalBlockId = neighborhood[1][1][1];
  vtkDataObject* block = input->GetBlock(this->GlobalToLocalMap[globalBlockId]);
  vtkImageData* image = vtkImageData::SafeDownCast(block);
  if (image)
  {
    this->CurrentLevel = this->GlobalLevelArray[globalBlockId];
    this->CurrentBlockId = globalBlockId;
    // Recursively find the maximum depth of the children branches (not loaded).
    this->RemainingDepth = this->ComputeBranchDepth(globalBlockId);

    this->ProcessBlock(image);
    // Now lets process the regions shared with neighbors.
    int r[3];
    for (r[2] = 0; r[2] < 3; ++r[2])
    {
      for (r[1] = 0; r[1] < 3; ++r[1])
      {
        for (r[0] = 0; r[0] < 3; ++r[0])
        {
          if (r[0] != 1 || r[1] != 1 || r[2] != 1)
          {
            this->ProcessNeighborhoodSharedRegion(neighborhood, r, input);
          }
        }
      }
    }
  }
}

//----------------------------------------------------------------------------
void vtkFlashContour::ProcessBlock(vtkImageData* image)
{
  const double* spacing = image->GetSpacing();
  double blockOrigin[3];
  image->GetOrigin(blockOrigin);
  int dims[3];

  // Shift origin half a pixel for the dual grid.
  blockOrigin[0] += 0.5 * spacing[0];
  blockOrigin[1] += 0.5 * spacing[1];
  blockOrigin[2] += 0.5 * spacing[2];

  // How to do the dual grid contour?
  // Simple marching cubes (marching corners for dual :)
  vtkDataArray* da = image->GetCellData()->GetArray(this->CellArrayNameToProcess);
  if (da->GetDataType() != VTK_DOUBLE)
  {
    vtkErrorMacro("Expecting doubles");
    return;
  }
  void* ptr = da->GetVoidPointer(0);
  double* dPtr = (double*)(ptr);
  // For passing / interpolating one double array.
  double* pPtr = 0;
  if (this->PassArray)
  {
    da = image->GetCellData()->GetArray(this->PassAttribute);
    if (da->GetDataType() != VTK_DOUBLE)
    {
      vtkErrorMacro("Expecting doubles");
      return;
    }
    ptr = da->GetVoidPointer(0);
    pPtr = (double*)(ptr);
  }

  double origin[3];

  image->GetDimensions(dims);

  // Change point dimensions to cell dimensions.
  dims[0] -= 1;
  dims[1] -= 1;
  dims[2] -= 1;
  // Compute offset to 8 neighbor cells (corners of dual cell).
  int yInc = dims[0];
  int zInc = yInc * dims[1];
  double cornerValues[8];
  double passValues[8];
  int cornerOffsets[8];
  cornerOffsets[0] = 0;
  cornerOffsets[1] = 1;
  cornerOffsets[2] = 1 + yInc;
  cornerOffsets[3] = yInc;
  cornerOffsets[4] = zInc;
  cornerOffsets[5] = 1 + zInc;
  cornerOffsets[6] = 1 + yInc + zInc;
  cornerOffsets[7] = yInc + zInc;

  // Change cell dims to dual cell dimensions.
  dims[0] -= 1;
  dims[1] -= 1;
  dims[2] -= 1;
  // Loop through the dual cells.
  origin[2] = blockOrigin[2];
  for (int z = 0; z < dims[2]; ++z)
  {
    origin[1] = blockOrigin[1];
    for (int y = 0; y < dims[1]; ++y)
    {
      origin[0] = blockOrigin[0];
      for (int x = 0; x < dims[0]; ++x)
      {
        cornerValues[0] = dPtr[cornerOffsets[0]];
        cornerValues[1] = dPtr[cornerOffsets[1]];
        cornerValues[2] = dPtr[cornerOffsets[2]];
        cornerValues[3] = dPtr[cornerOffsets[3]];
        cornerValues[4] = dPtr[cornerOffsets[4]];
        cornerValues[5] = dPtr[cornerOffsets[5]];
        cornerValues[6] = dPtr[cornerOffsets[6]];
        cornerValues[7] = dPtr[cornerOffsets[7]];
        if (pPtr)
        {
          passValues[0] = pPtr[cornerOffsets[0]];
          passValues[1] = pPtr[cornerOffsets[1]];
          passValues[2] = pPtr[cornerOffsets[2]];
          passValues[3] = pPtr[cornerOffsets[3]];
          passValues[4] = pPtr[cornerOffsets[4]];
          passValues[5] = pPtr[cornerOffsets[5]];
          passValues[6] = pPtr[cornerOffsets[6]];
          passValues[7] = pPtr[cornerOffsets[7]];
        }

        // I adding interpolation of attributes after the fact.
        // I need ids of the corner cells (dual points).
        this->ProcessCell(origin, spacing, cornerValues, passValues);
        ++dPtr;
        if (pPtr)
        {
          ++pPtr;
        }
        origin[0] += spacing[0];
      }
      // Because we are looping over cells and pointer is for corners,
      // We need to skip the last (unshared) corner (max faces).
      // Skip the last corner in X row.
      ++dPtr;
      if (pPtr)
      {
        ++pPtr;
      }
      origin[1] += spacing[1];
    }
    // Skip the last row of Y
    dPtr += yInc;
    if (pPtr)
    {
      pPtr += yInc;
    }
    origin[2] += spacing[2];
  }
  // The last z face is skip automatically.
}

//----------------------------------------------------------------------------
// Assume the same level: easy.
void vtkFlashContour::ProcessNeighborhoodSharedRegion(
  int neighborhood[3][3][3], int r[3], vtkMultiBlockDataSet* input)
{
  int regionDims[3];   // dual cell dimensions of region
  double* ptrs[8];     // Pointer to corner scalars
  double* aptrs[8];    // Pointer to attribute to pass
  double corners[32];  // world location of corners
  double spacings[32]; // spacing for each corner block.
  int levelDiff[8];    // level of corner relative to center block
  int incs[3];         // all blocks have the same increments

  // Assumptions that all blocks have the same increments is only
  // to minimize the arguments to the ProcessSharedRegion method.
  int block1GlobalId = neighborhood[1][1][1];
  int level1 = this->GlobalLevelArray[block1GlobalId];
  vtkDataObject* block1 = input->GetBlock(this->GlobalToLocalMap[block1GlobalId]);
  vtkImageData* image1 = vtkImageData::SafeDownCast(block1);
  const int* dims1 = image1->GetDimensions();
  const double* spacing1 = image1->GetSpacing();
  const double* origin1 = image1->GetOrigin();
  // Compute increments for cell array (cell array is one less than point).
  incs[0] = 1;
  incs[1] = (dims1[0] - 1);
  incs[2] = incs[1] * (dims1[1] - 1);

  int cornerNeighborhoodAxisIndex[2][3];
  int cornerDualPointAxisIndex[2][3];

  // Find the dimensions of the dual cell shared region we need to process.
  for (int i = 0; i < 3; ++i)
  {
    // Now we need to find neighbor index and index (offset) for each corner.
    switch (r[i])
    {
      case 0:
        cornerNeighborhoodAxisIndex[0][i] = 0;
        cornerNeighborhoodAxisIndex[1][i] = 1;
        cornerDualPointAxisIndex[0][i] = dims1[i] - 2;
        cornerDualPointAxisIndex[1][i] = 0;
        regionDims[i] = 1; // Collapse axis to get face/edge/corner region dimensions.
        break;
      case 1:
        cornerNeighborhoodAxisIndex[0][i] = 1;
        cornerNeighborhoodAxisIndex[1][i] = 1;
        cornerDualPointAxisIndex[0][i] = 0;
        cornerDualPointAxisIndex[1][i] = 1;
        regionDims[i] =
          dims1[i] - 2; // -1 to convert point to cell, -1 to convert cell to dual cell.
        break;
      case 2:
        cornerNeighborhoodAxisIndex[0][i] = 1;
        cornerNeighborhoodAxisIndex[1][i] = 2;
        cornerDualPointAxisIndex[0][i] = dims1[i] - 2;
        cornerDualPointAxisIndex[1][i] = 0;
        regionDims[i] = 1; // Collapse axis to get face/edge/corner region dimensions.
        break;
      default:
        vtkErrorMacro("Bad neighbor index.");
    }
  }
  // Note:  using index 1 and 2 (i.e. image1 and image2) is a bit confusing.
  // image1 is the center block of the neighborhood (owner of the region).
  // image2 is for the block found for the corner in the corner loop.
  // Each corner can have a different block (not just two blocks).
  // Block2 may (and is usually) the same as block1.

  // Now loop over corners and fill in information we need to march.
  for (int cornerId = 0; cornerId < 8; ++cornerId)
  {
    int cx = (cornerId & 1) ? 1 : 0;
    int cy = (cornerId & 2) ? 1 : 0;
    int cz = (cornerId & 4) ? 1 : 0;
    int block2GlobalId =
      neighborhood[cornerNeighborhoodAxisIndex[cx][0]][cornerNeighborhoodAxisIndex[cy][1]]
                  [cornerNeighborhoodAxisIndex[cz][2]];
    if (block2GlobalId < 0)
    { // Missing neighbor. Skip this shared region.
      return;
    }
    // If the corner block has children then the children have priority for
    // ownership of this region and we should not process it.
    int* block2Children = this->GlobalChildrenArray + (block2GlobalId << 3);
    if (this->GlobalToLocalMap[block2Children[0]] >= 0)
    {
      return;
    }

    int level2 = this->GlobalLevelArray[block2GlobalId];
    // Skip the region if the center block (1) does not own it.
    if (level2 > level1 || (level2 == level1 && block2GlobalId < block1GlobalId))
    { // block2 dominates block1 in ownership.
      return;
    }
    levelDiff[cornerId] = level1 - level2;
    // Get the image for block2
    int block2LocalId = this->GlobalToLocalMap[block2GlobalId];
    if (block2LocalId < 0)
    { // Missing neighbor. Skip this shared region.
      return;
    }
    vtkDataObject* block2 = input->GetBlock(block2LocalId);
    vtkImageData* image2 = vtkImageData::SafeDownCast(block2);
    if (image2 == 0)
    {
      return;
    }
    // Sanity check. All blocks must have the same dimensions.
    int* dims2 = image2->GetDimensions();
    if (dims1[0] != dims2[0] || dims1[1] != dims2[1] || dims1[2] != dims2[2])
    {
      vtkErrorMacro("Neighbor dimensions do not match.");
      return;
    }
    // Origin of block that contains the corner.
    const double* origin2 = image2->GetOrigin();
    // Spacing of block that contains corner.
    double* spacing2 = spacings + (cornerId << 2);
    ;
    image2->GetSpacing(spacing2);
    // Find world location for dual point / corner.
    double* cornerPoint = corners + (cornerId << 2);
    // Find index of cell in block for corner.
    int dualPoint2Index[3];
    if (levelDiff[cornerId] == 0)
    { // Same level, so just use precomputed index.
      dualPoint2Index[0] = cornerDualPointAxisIndex[cx][0];
      dualPoint2Index[1] = cornerDualPointAxisIndex[cy][1];
      dualPoint2Index[2] = cornerDualPointAxisIndex[cz][2];
      // + 0.5 to move to the center of the cell (dual point).
      cornerPoint[0] = origin2[0] + spacing2[0] * (dualPoint2Index[0] + 0.5);
      cornerPoint[1] = origin2[1] + spacing2[1] * (dualPoint2Index[1] + 0.5);
      cornerPoint[2] = origin2[2] + spacing2[2] * (dualPoint2Index[2] + 0.5);
    }
    else
    { // Translate into the coordinates of the lower level block.
      // First fing the world location of the corner assuming the same level.
      // Note: we use dualPoint2Index to temporarily hold dualPoint1Index.
      dualPoint2Index[0] =
        (cornerNeighborhoodAxisIndex[cx][0] - 1) * (dims1[0] - 1) + cornerDualPointAxisIndex[cx][0];
      dualPoint2Index[1] =
        (cornerNeighborhoodAxisIndex[cy][1] - 1) * (dims1[1] - 1) + cornerDualPointAxisIndex[cy][1];
      dualPoint2Index[2] =
        (cornerNeighborhoodAxisIndex[cz][2] - 1) * (dims1[2] - 1) + cornerDualPointAxisIndex[cz][2];
      cornerPoint[0] = origin1[0] + spacing1[0] * (dualPoint2Index[0] + 0.5);
      cornerPoint[1] = origin1[1] + spacing1[1] * (dualPoint2Index[1] + 0.5);
      cornerPoint[2] = origin1[2] + spacing1[2] * (dualPoint2Index[2] + 0.5);
      // Now find the index of the block2 cell containing the point.
      dualPoint2Index[0] = (int)((cornerPoint[0] - origin2[0]) / spacing2[0]);
      dualPoint2Index[1] = (int)((cornerPoint[1] - origin2[1]) / spacing2[1]);
      dualPoint2Index[2] = (int)((cornerPoint[2] - origin2[2]) / spacing2[2]);
      // Now compute the correct (degenerate) dual point location for lower level.
      cornerPoint[0] = origin2[0] + spacing2[0] * (dualPoint2Index[0] + 0.5);
      cornerPoint[1] = origin2[1] + spacing2[1] * (dualPoint2Index[1] + 0.5);
      cornerPoint[2] = origin2[2] + spacing2[2] * (dualPoint2Index[2] + 0.5);
    }
    // Get the pointer for the corner.
    vtkDataArray* da = image2->GetCellData()->GetArray(this->CellArrayNameToProcess);
    if (da->GetDataType() != VTK_DOUBLE)
    {
      vtkErrorMacro("Expecting doubles");
      return;
    }
    ptrs[cornerId] = (double*)(da->GetVoidPointer(0));
    // Move pointer to the correct position.
    ptrs[cornerId] +=
      incs[0] * dualPoint2Index[0] + incs[1] * dualPoint2Index[1] + incs[2] * dualPoint2Index[2];
    // For passing attributes.
    if (this->PassArray)
    {
      da = image2->GetCellData()->GetArray(this->PassAttribute);
      if (da->GetDataType() != VTK_DOUBLE)
      {
        vtkErrorMacro("Expecting doubles");
        return;
      }
      aptrs[cornerId] = (double*)(da->GetVoidPointer(0));
      // Move pointer to the correct position.
      aptrs[cornerId] +=
        incs[0] * dualPoint2Index[0] + incs[1] * dualPoint2Index[1] + incs[2] * dualPoint2Index[2];
    }
  }
  // Now that we have all of the information for the starting cell corners
  // Contour the region.
  this->ProcessSharedRegion(regionDims, ptrs, incs, corners, spacings, levelDiff, aptrs);
}

//----------------------------------------------------------------------------
// cornerPtr and cornerPoints get modified.
void vtkFlashContour::ProcessSharedRegion(int regionDims[3], double* cornerPtrs[8], int incs[3],
  double cornerPoints[32], double cornerSpacings[32], int cornerLevelDiffs[8], double* passPtrs[8])
{
  // Skip schedule for lower levels.
  // The 2's have not effect when levelDiff = 0.
  // When levelDiff > 0, they cause the max corner on an axis to move to
  // the next pointer/location one step ahead of the min corner.
  int levelCountX[8] = { 1, 2, 1, 2, 1, 2, 1, 2 };
  int levelCountY[8] = { 1, 1, 2, 2, 1, 1, 2, 2 };
  int levelCountZ[8] = { 1, 1, 1, 1, 2, 2, 2, 2 };

  // We have to keep y and x values for all cornerPtrs and cornerPoints.
  // There should be no problems overwriting arrays passed by using
  // cornerPtrs for cornerPtrsZ and cornerPoints for cornerPointsZ.
  double* cornerPtrsY[8];
  double* cornerPtrsX[8];
  double cornerPointsY[32];
  double cornerPointsX[32];
  double* passPtrsY[8];
  double* passPtrsX[8];

  // Loop over the region.
  for (int z = 0; z < regionDims[2]; ++z)
  {
    memcpy(cornerPtrsY, cornerPtrs, 8 * sizeof(double*));
    memcpy(cornerPointsY, cornerPoints, 32 * sizeof(double));
    if (this->PassArray)
    {
      memcpy(passPtrsY, passPtrs, 8 * sizeof(double*));
    }
    for (int y = 0; y < regionDims[1]; ++y)
    {
      memcpy(cornerPtrsX, cornerPtrsY, 8 * sizeof(double*));
      memcpy(cornerPointsX, cornerPointsY, 32 * sizeof(double));
      if (this->PassArray)
      {
        memcpy(passPtrsX, passPtrsY, 8 * sizeof(double*));
      }
      for (int x = 0; x < regionDims[0]; ++x)
      {
        this->ProcessDegenerateCell(cornerPointsX, cornerPtrsX, passPtrsX);
        // Increment x corners
        for (int i = 0; i < 8; ++i)
        {
          if (++levelCountX[i] > (1 << cornerLevelDiffs[i]))
          { // Increment
            cornerPtrsX[i] += incs[0];
            cornerPointsX[i << 2] += cornerSpacings[i << 2];
            levelCountX[i] = 1;
            if (this->PassArray)
            {
              passPtrsX[i] += incs[0];
            }
          }
        }
      }
      // Increment y corners
      for (int i = 0; i < 8; ++i)
      {
        if (++levelCountY[i] > (1 << cornerLevelDiffs[i]))
        { // Increment
          if (this->PassArray)
          {
            passPtrsY[i] += incs[1];
          }
          cornerPtrsY[i] += incs[1];
          cornerPointsY[(i << 2) | 1] += cornerSpacings[(i << 2) | 1];
          levelCountY[i] = 1;
        }
        // This is the most convient place to reset the levelCounts
        // after processing a row.
        levelCountX[i] = (i & 1) ? 2 : 1;
      }
    }
    // Increment corners
    for (int i = 0; i < 8; ++i)
    {
      if (++levelCountZ[i] > (1 << cornerLevelDiffs[i]))
      { // Increment
        if (this->PassArray)
        {
          passPtrs[i] += incs[2];
        }
        cornerPtrs[i] += incs[2];
        cornerPoints[(i << 2) | 2] += cornerSpacings[(i << 2) | 2];
        levelCountZ[i] = 1;
      }
      // This is the most convient place to reset the levelCounts
      // after processing a plane.
      levelCountY[i] = (i & 2) ? 2 : 1;
    }
  }
}

//----------------------------------------------------------------------------
void vtkFlashContour::ProcessDegenerateCell(
  double cornerPoints[32], double* cornerPtrs[8], double* passPtrs[8])
{
  int cubeCase = 0;
  double cornerValues[8];
  double passValues[8];

  cornerValues[0] = *cornerPtrs[0];
  cornerValues[1] = *cornerPtrs[1];
  cornerValues[2] = *cornerPtrs[3];
  cornerValues[3] = *cornerPtrs[2];
  cornerValues[4] = *cornerPtrs[4];
  cornerValues[5] = *cornerPtrs[5];
  cornerValues[6] = *cornerPtrs[7];
  cornerValues[7] = *cornerPtrs[6];

  cubeCase = 0;
  if (cornerValues[0] > this->IsoValue)
  {
    cubeCase += 1;
  }
  if (cornerValues[1] > this->IsoValue)
  {
    cubeCase += 2;
  }
  if (cornerValues[2] > this->IsoValue)
  {
    cubeCase += 4;
  }
  if (cornerValues[3] > this->IsoValue)
  {
    cubeCase += 8;
  }
  if (cornerValues[4] > this->IsoValue)
  {
    cubeCase += 16;
  }
  if (cornerValues[5] > this->IsoValue)
  {
    cubeCase += 32;
  }
  if (cornerValues[6] > this->IsoValue)
  {
    cubeCase += 64;
  }
  if (cornerValues[7] > this->IsoValue)
  {
    cubeCase += 128;
  }

  // I am trying to exit as quick as possible if there is
  // no surface to generate.  I could also check that the index
  // is not on boundary.
  if (cubeCase == 0 || cubeCase == 255)
  {
    return;
  }

  if (this->PassArray && passPtrs)
  {
    passValues[0] = *passPtrs[0];
    passValues[1] = *passPtrs[1];
    passValues[2] = *passPtrs[3];
    passValues[3] = *passPtrs[2];
    passValues[4] = *passPtrs[4];
    passValues[5] = *passPtrs[5];
    passValues[6] = *passPtrs[7];
    passValues[7] = *passPtrs[6];
  }

  this->ProcessCellFinal(cornerPoints, cornerValues, cubeCase, passValues);
}

//----------------------------------------------------------------------------
void vtkFlashContour::ProcessCell(
  const double* origin, const double* spacing, const double* cornerValues, const double* passValues)
{
  int cubeCase = 0;

  cubeCase = 0;
  if (cornerValues[0] > this->IsoValue)
  {
    cubeCase += 1;
  }
  if (cornerValues[1] > this->IsoValue)
  {
    cubeCase += 2;
  }
  if (cornerValues[2] > this->IsoValue)
  {
    cubeCase += 4;
  }
  if (cornerValues[3] > this->IsoValue)
  {
    cubeCase += 8;
  }
  if (cornerValues[4] > this->IsoValue)
  {
    cubeCase += 16;
  }
  if (cornerValues[5] > this->IsoValue)
  {
    cubeCase += 32;
  }
  if (cornerValues[6] > this->IsoValue)
  {
    cubeCase += 64;
  }
  if (cornerValues[7] > this->IsoValue)
  {
    cubeCase += 128;
  }

  // I am trying to exit as quick as possible if there is
  // no surface to generate.  I could also check that the index
  // is not on boundary.
  if (cubeCase == 0 || cubeCase == 255)
  {
    return;
  }

  double cornerPoints[32]; // 4 is easier to optimize than 3.
  // Loop over the corners.
  for (int c = 0; c < 8; ++c)
  {
    int px, py, pz; // Corner global xyz index.
    // CornerIndex
    px = (c & 1) ? 1 : 0;
    py = (c & 2) ? 1 : 0;
    pz = (c & 4) ? 1 : 0;
    cornerPoints[c << 2] = origin[0] + spacing[0] * ((double)(px));
    cornerPoints[(c << 2) | 1] = origin[1] + spacing[1] * ((double)(py));
    cornerPoints[(c << 2) | 2] = origin[2] + spacing[2] * ((double)(pz));
  }

  this->ProcessCellFinal(cornerPoints, cornerValues, cubeCase, passValues);
}

//----------------------------------------------------------------------------
// It appears that cornerValues use VTK indexing scheme but
// cornerPoints does not.
void vtkFlashContour::ProcessCellFinal(const double cornerPoints[32], const double cornerValues[8],
  int cubeCase, const double passValues[8])
{
  vtkIdType pointIds[6];
  vtkMarchingCubesTriangleCases *triCase, *triCases;
  EDGE_LIST* edge;
  double k, v0, v1;
  triCases = vtkMarchingCubesTriangleCases::GetCases();

  // We have the points, now contour the cell.
  // Get edges.
  triCase = triCases + cubeCase;
  edge = triCase->edges;
  double pt[3];

  // loop over triangles
  while (*edge > -1)
  {
    // I want to avoid adding degenerate triangles.
    // Maybe the best way to do this is to have a point locator
    // merge points first.
    // Create brute force locator for a block, and resuse it.
    // Only permanently keep locator for edges shared between two blocks.
    for (int ii = 0; ii < 3; ++ii, ++edge) // insert triangle
    {
      vtkIdType ptId = -1; // = this->BlockLocator->GetEdgePointer(x,y,z,*edge);

      if (ptId == -1)
      {
        // Compute the interpolation factor.
        v0 = cornerValues[vtkFlashIsoEdgeToVTKPointsTable[*edge][0]];
        v1 = cornerValues[vtkFlashIsoEdgeToVTKPointsTable[*edge][1]];
        k = (this->IsoValue - v0) / (v1 - v0);
        // Add the point to the output and get the index of the point.
        int pt1Idx = (vtkFlashIsoEdgeToPointsTable[*edge][0] << 2);
        int pt2Idx = (vtkFlashIsoEdgeToPointsTable[*edge][1] << 2);
        // I wonder if this is any faster than incrementing a pointer.
        pt[0] = cornerPoints[pt1Idx] + k * (cornerPoints[pt2Idx] - cornerPoints[pt1Idx]);
        pt[1] =
          cornerPoints[pt1Idx | 1] + k * (cornerPoints[pt2Idx | 1] - cornerPoints[pt1Idx | 1]);
        pt[2] =
          cornerPoints[pt1Idx | 2] + k * (cornerPoints[pt2Idx | 2] - cornerPoints[pt1Idx | 2]);
        ptId = this->Points->InsertNextPoint(pt);

        if (this->PassArray)
        {
          double p0;
          double p1;
          p0 = passValues[vtkFlashIsoEdgeToVTKPointsTable[*edge][0]];
          p1 = passValues[vtkFlashIsoEdgeToVTKPointsTable[*edge][1]];
          double value = p0 + k * (p1 - p0);
          this->PassArray->InsertNextValue(value);
        }
      }
      pointIds[ii] = ptId;
    }
    if (pointIds[0] != pointIds[1] && pointIds[0] != pointIds[2] && pointIds[1] != pointIds[2])
    {
      this->Faces->InsertNextCell(3, pointIds);
      this->BlockIdCellArray->InsertNextValue(this->CurrentBlockId);
      this->LevelCellArray->InsertNextValue(this->CurrentLevel);
      this->RemainingDepthCellArray->InsertNextValue(this->RemainingDepth);
    }
  }
}
