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

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

#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkPolyData.h"
#include "vtkBitArray.h"
#include "vtkHyperTreeGrid.h"
#include "vtkDataSetAttributes.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkRenderer.h"
#include "vtkCamera.h"

#include "vtkIncrementalPointLocator.h"
#include "vtkMergePoints.h"

#include "vtkHyperTreeGridNonOrientedGeometryCursor.h"
#include "vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight.h"

static const unsigned int VonNeumannCursors3D[] = { 0, 1, 2, 4, 5, 6 };
static const unsigned int VonNeumannOrientations3D[]  = { 2, 1, 0, 0, 1, 2 };
static const unsigned int VonNeumannOffsets3D[]  = { 0, 0, 0, 1, 1, 1 };

vtkStandardNewMacro(vtkAdaptiveDataSetSurfaceFilter);

//-----------------------------------------------------------------------------
vtkAdaptiveDataSetSurfaceFilter::vtkAdaptiveDataSetSurfaceFilter()
{
  this->InData = nullptr;
  this->OutData = nullptr;
  this->Points = nullptr;
  this->Cells = nullptr;

  // Default dimension is 0
  this->Dimension = 0;

  // Default orientation is 0
  this->Orientation = 0;

  this->Renderer = nullptr;

  this->LevelMax = -1;

  this->ViewPointDepend = true;

  this->ParallelProjection = false;
  this->LastRendererSize[0] = 0;
  this->LastRendererSize[1] = 0;
  this->LastCameraFocalPoint[0] = 0.0;
  this->LastCameraFocalPoint[1] = 0.0;
  this->LastCameraFocalPoint[2] = 0.0;
  this->LastCameraParallelScale = 0;

  this->Scale = 1;

  this->CircleSelection = true;
  this->BBSelection = false;
  this->FixedLevelMax = -1;
  this->DynamicDecimateLevelMax = 0;

  // Default Locator is 0
  this->Merging = false;
  this->Locator = nullptr;
}

//-----------------------------------------------------------------------------
vtkAdaptiveDataSetSurfaceFilter::~vtkAdaptiveDataSetSurfaceFilter()
{
}

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

  if( this->InData )
  {
    os << indent << "InData:\n";
    this->InData->PrintSelf( os, indent.GetNextIndent() );
  }
  else
  {
    os << indent << "InData: ( none )\n";
  }

  if( this->OutData )
  {
    os << indent << "OutData:\n";
    this->OutData->PrintSelf( os, indent.GetNextIndent() );
  }
  else
  {
    os << indent << "OutData: ( none )\n";
  }

  if( this->Points )
  {
    os << indent << "Points:\n";
    this->Points->PrintSelf( os, indent.GetNextIndent() );
  }
  else
  {
    os << indent << "Points: ( none )\n";
  }

  if( this->Cells )
  {
    os << indent << "Cells:\n";
    this->Cells->PrintSelf( os, indent.GetNextIndent() );
  }
  else
  {
    os << indent << "Cells: ( none )\n";
  }

  os << indent << "Dimension: " << this->Dimension << endl;
  os << indent << "Orientation: " << this->Orientation << endl;
  os << indent << "Axis1: " << this->Axis1 << endl;
  os << indent << "Axis2: " << this->Axis2 << endl;
  os << indent << "Radius: " << this->Radius << endl;
  os << indent << "LevelMax: " << this->LevelMax << endl;
  os << indent << "ViewPointDepend: " << this->ViewPointDepend << endl;
  os << indent << "ParallelProjection: " << this->ParallelProjection << endl;
  os << indent << "Scale: " << this->Scale << endl;
  os << indent << "FixedLevelMax: " << this->FixedLevelMax << endl;
  os << indent << "DynamicDecimateLevelMax: " << this->DynamicDecimateLevelMax << endl;
  os << indent << "LastCameraParallelScale: " << this->LastCameraParallelScale << endl;
  os << indent << "LastRendererSize: "
     << this->LastRendererSize[0] << ", "
     << this->LastRendererSize[1] << endl;
  os << indent << "LastCameraFocalPoint: "
     << this->LastCameraFocalPoint[0] << ", "
     << this->LastCameraFocalPoint[1] << ", "
     << this->LastCameraFocalPoint[2] << endl;
}

//----------------------------------------------------------------------------
int vtkAdaptiveDataSetSurfaceFilter::RequestData( vtkInformation* 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
  vtkDataSet *input = vtkDataSet::SafeDownCast(
    inInfo->Get( vtkDataObject::DATA_OBJECT() ) );
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get( vtkDataObject::DATA_OBJECT() ) );

  int MeshType = input->GetDataObjectType();
  if ( MeshType != VTK_HYPER_TREE_GRID )
  {
    return this->Superclass::RequestData( request, inputVector, outputVector );
  }

  return this->DataSetExecute( input, output );
}

//----------------------------------------------------------------------------
int vtkAdaptiveDataSetSurfaceFilter::DataSetExecute( vtkDataSet* inputDS,
                                                     vtkPolyData* output )
{
  // Retrieve input grid
  vtkHyperTreeGrid* input = vtkHyperTreeGrid::SafeDownCast( inputDS );
  if ( ! input )
  {
    return vtkDataSetSurfaceFilter::DataSetExecute( inputDS, output );
  }

  // Retrieve useful grid parameters for speed of access
  this->Dimension = input->GetDimension();
  this->Orientation = input->GetOrientation();

  // Initialize output cell data
  this->InData = static_cast<vtkDataSetAttributes*>( input->GetPointData() );
  this->OutData = static_cast<vtkDataSetAttributes*>( output->GetCellData() );
  this->OutData->CopyAllocate( this->InData );
  if ( this->PassThroughCellIds )
  {
    this->OriginalCellIds = vtkIdTypeArray::New();
    this->OriginalCellIds->SetName( this->GetOriginalCellIdsName() );
    this->OriginalCellIds->SetNumberOfComponents( 1 );
    this->OutData->AddArray( this->OriginalCellIds );
  }

  // Init renderer information
  if ( this->ViewPointDepend && this->ParallelProjection && this->Renderer )
  {
    // Generate planes XY, XZ o YZ
    unsigned int* gridSize = input->GetGridSize();
    bool isInit = false;
    if ( gridSize[0] == 1 )
    {
      this->Axis1 = 1;
      this->Axis2 = 2;
      isInit = true;
    }
    else if ( gridSize[1] == 1 )
    {
      this->Axis1 = 0;
      this->Axis2 = 2;
      isInit = true;
    }
    else if ( gridSize[2] == 1 )
    {
      this->Axis1 = 0;
      this->Axis2 = 1;
      isInit = true;
    }

    // Compute the zoom of the camera
    vtkCamera* cam = this->Renderer->GetActiveCamera();

    // Compute the bounding box
    double bounds[6];
    input->GetBounds( bounds );

    // JB Recupere le branch factor
    int f = input->GetBranchFactor();

    // JB Le calcul qui suit a pour objet de determiner le niveau de parcour en profondeur utile pour l'affichage
    // JB en fonction de la distance
    if ( isInit )
    {
      // JB Taille Moyenne d'une maille du niveau 0 dans les coordonnees reelles suivant chaque direction
      double worldCellAverageScaleAxis1 = ( bounds[2 * this->Axis1 + 1] - bounds[2 * this->Axis1] ) / (double) ( gridSize[this->Axis1] ) / this->Scale;
      double worldCellAverageScaleAxis2 = ( bounds[2 * this->Axis2 + 1] - bounds[2 * this->Axis2] ) / (double) ( gridSize[this->Axis2] ) / this->Scale;

      // JB Taille de la fenetre dans les coordonnees reelles (GetParallelScale) suivant chaque direction
      double worldWindScaleAxis1 = cam->GetParallelScale() * this->LastRendererSize[0] / (double) ( this->LastRendererSize[1] );
      double worldWindScaleAxis2 = cam->GetParallelScale();

      // JB Taille de la fenetre en pixel, ecran, suivant chaque direction
      double windScaleAxis1 = this->LastRendererSize[0];
      double windScaleAxis2 = this->LastRendererSize[1];

      // JB Compute how many levels of the tree we should process by direction
      // JB 1) application du theoreme de Thales ; le rapport taille d'une maille de niveau 0 par la taille de la fenetre est identique si le
      // JB calcul se fait en coordonnees reelles ou en coordonnees ecran (pixel)
      // JB 2) la taille d'une maille de niveau L vaut la taille d'une maille de niveau 0 divisee par le facteur de raffinement eleve a la puissance L
      // JB 3) au final, le calcul qui suit a pour objet de determiner quand une maille fera un pixel
      double levelMaxiAxis1 = ( log ( windScaleAxis1 ) + log ( worldCellAverageScaleAxis1 ) - log ( worldWindScaleAxis1 ) ) / log ( f );
      double levelMaxiAxis2 = ( log ( windScaleAxis2 ) + log ( worldCellAverageScaleAxis2 ) - log ( worldWindScaleAxis2 ) ) / log ( f );

      // JB On opte pour le niveau le plus eleve
      this->LevelMax = std::ceil( std::max( levelMaxiAxis1, levelMaxiAxis2 ) );
    } else {
      // JB En 3D, par defaut, on prend tous les niveaux
      this->LevelMax = 65536;
    }

    // JB Par option, on peut reduire cette valeur... tres utile pour avoir un LOD leger.
    this->LevelMax -= this->DynamicDecimateLevelMax;
    if ( this->LevelMax < 0 )
    {
      this->LevelMax = 0;
    }

    // JB Par option, on peut fixer le niveau max independemment du calcul dynamique realise precedemment
    if ( this->FixedLevelMax >= 0 )
    {
      this->LevelMax = this->FixedLevelMax;
    }

    // JB Le calcul qui suit a pour objet de determiner le rayon du cercle dans les coordonnees reelles incluant la projection
    // JB de la fenetre. L'activation de CircleSelection permettra de ne produire que les mailles intersectant ce cercle centre
    // JB au camera focal point.
    // JB LastCameraFocalPoint retourne le centre de l'ecran dans les coordonnees reelles
    double ratio = this->LastRendererSize[0] / (double) ( this->LastRendererSize[1] );
    this->Radius = cam->GetParallelScale() * sqrt( 1 + pow(  ratio, 2 ) );

    // JB Le calcul qui suit a pour objet de determiner la boite englobante dans les coordonnees reelles (et sans tenir compte
    // JB d'un point de vue qui aurait tourne) incluant la projection de la fenetre. L'activation de BBSelection permettra de ne
    // JB produire que les mailles intersectant cette boite englobante.
    this->WindowBounds[0] = this->LastCameraFocalPoint[ 0 ] - cam->GetParallelScale() * ratio;
    this->WindowBounds[1] = this->LastCameraFocalPoint[ 0 ] + cam->GetParallelScale() * ratio;
    this->WindowBounds[2] = this->LastCameraFocalPoint[ 1 ] - cam->GetParallelScale();
    this->WindowBounds[3] = this->LastCameraFocalPoint[ 1 ] + cam->GetParallelScale();

#ifndef NDEBUG
    this->NbRejectByCircle = 0;
    this->NbRejectByBB = 0;

    std::cerr << "LevelMax        " << this->LevelMax << std::endl;
    std::cerr << "CircleSelection " << this->CircleSelection << std::endl;
    std::cerr << "Circle R        " << this->Radius << std::endl;
    std::cerr << "       CX       " << this->LastCameraFocalPoint[this->Axis1] << std::endl;
    std::cerr << "       CY       " << this->LastCameraFocalPoint[this->Axis2] << std::endl;
    std::cerr << "BBSelection     " << this->BBSelection << std::endl;
    std::cerr << "Bounds X        " << this->WindowBounds[0] << " : " << this->WindowBounds[1] << std::endl;
    std::cerr << "       Y        " << this->WindowBounds[2] << " : " << this->WindowBounds[3] << std::endl;
#endif
  }
  else
  {
    // Recurse all the tree
    this->LevelMax = -1;
  }

  // Extract geometry from hyper tree grid
  this->ProcessTrees( input, output );

  this->UpdateProgress( 1. );

  return 1;
}

//-----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::ProcessTrees( vtkHyperTreeGrid* input,
                                                    vtkPolyData* output )
{
  if ( this->Points )
  {
    this->Points->Delete();
  }
  // Create storage for corners of leaf cells
  this->Points = vtkPoints::New();

  // Create storage for untructured leaf cells
  if ( this->Cells )
  {
    this->Cells->Delete();
  }
  this->Cells = vtkCellArray::New();

  //JB Initialize a Locator
  if ( this->Merging )
  {
    this->Locator = vtkMergePoints::New();
    this->Locator->InitPointInsertion ( this->Points, input->GetBounds() );
  }

  // Retrieve material mask
  this->MaterialMask = input->HasMaterialMask() ? input->GetMaterialMask() : 0;

  //
  vtkUnsignedCharArray* ghost = input->GetPointGhostArray();
  if ( ghost )
  {
    this->OutData->CopyFieldOff( vtkDataSetAttributes::GhostArrayName() );
  }

  // Iterate over all hyper trees
  if ( this->Dimension == 3 )
  {
    vtkIdType index;
    vtkHyperTreeGrid::vtkHyperTreeGridIterator it;
    input->InitializeTreeIterator( it );
    vtkNew<vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight> cursor;
    while ( it.GetNextTree( index ) )
    {
      // In 3 dimensions, von Neumann neighborhood information is needed
      input->InitializeNonOrientedVonNeumannSuperCursorLight( cursor, index );
      // If this is not a ghost tree
      if ( ! ghost || ! ghost->GetTuple1( cursor->GetGlobalNodeIndex() ) )
      {
        // Build geometry recursively
        this->RecursivelyProcessTree3D( cursor, 0 );
      }
    } // it
  } // if ( this->Dimension == 3 )
  else
  {
    vtkIdType index;
    vtkHyperTreeGrid::vtkHyperTreeGridIterator it;
    input->InitializeTreeIterator( it );
    vtkNew<vtkHyperTreeGridNonOrientedGeometryCursor> cursor;
    while ( it.GetNextTree( index ) )
    {
      // Otherwise, geometric properties of the cells suffice
      input->InitializeNonOrientedGeometryCursor( cursor, index );
      // If this is not a ghost tree
      if ( ! ghost || ! ghost->GetTuple1( cursor->GetGlobalNodeIndex() ) )
      {
        // Build geometry recursively
        this->RecursivelyProcessTreeNot3D( cursor, 0 );
      }
    } // it
  } // else

  // Set output geometry and topology
  output->SetPoints( this->Points );
  if ( this->Dimension == 1  )
  {
    output->SetLines( this->Cells );
  }
  else
  {
    output->SetPolys( this->Cells );
  }

#ifndef NDEBUG
  std::cerr << "vtkAdaptiveDataSetSurfaceFilter #Points            " << this->Points->GetNumberOfPoints() << std::endl;
  std::cerr << "                                #Cells             " << this->Cells->GetNumberOfCells() << std::endl;
  std::cerr << "                                #Type&Connectivity " << this->Cells->GetNumberOfConnectivityEntries() << std::endl;
  std::cerr << "                          Cells #NbRejectByBB      " << this->NbRejectByBB << std::endl;
  std::cerr << "                                #NbRejectByCircle  " << this->NbRejectByCircle << std::endl;
#endif
  std::cerr << "vtkAdaptiveDataSetSurfaceFilter #Points            " << this->Points->GetNumberOfPoints() << std::endl;
  std::cerr << "                                #Cells             " << this->Cells->GetNumberOfCells() << std::endl;
  std::cerr << "                                #Type&Connectivity " << this->Cells->GetNumberOfConnectivityEntries() << std::endl;

  this->Points->Delete();
  this->Points = NULL;
  this->Cells->Delete();
  this->Cells = NULL;

  if ( this->Locator )
  {
    this->Locator->UnRegister(this);
    this->Locator = nullptr;
  }
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::RecursivelyProcessTreeNot3D(
  vtkHyperTreeGridNonOrientedGeometryCursor* cursor,
  int level )
{
  bool insideBB = ( this->LevelMax == -1 );
  if( ! insideBB && ( this->CircleSelection || this->BBSelection ) )
  {
    double originAxis1 = cursor->GetOrigin()[ this->Axis1 ];
    double originAxis2 = cursor->GetOrigin()[ this->Axis2 ];
    double halfAxis1 = cursor->GetSize()[ this->Axis1 ] / 2;
    double halfAxis2 = cursor->GetSize()[ this->Axis2 ] / 2;
    if ( this->CircleSelection )
    {
      //JB On determine si la maille correspondant au current node of the tree
      //JB is going to be renedered.
      //JB Pour cela, on fait une premiere approximation en considerant la maille
      //JB carre qui l'englobe en conservant la meme origine et en fixant sa demi-largeur
      //JB a la valeur maximale entre les valeurs de demi-largeur et demi-longueur.
      double half = std::max( halfAxis1, halfAxis2 );
      //JB This cell must to be rendered si le centre de cette maille se trouve a moins de
      //JB Radius + half * sqrt(2) du camera focal point. Radius est le rayon minimal du cercle
      //JB centre sur la camera focal point couvrant la fenetre de rendu.
      //JB Le centre de la maille se trouve a Origin + half, par direction.
      //JB La comparaison se fait sur les distances au carre afin d'eviter le calcul
      //JB couteux de racines carres.
      insideBB = (
        pow( originAxis1 + half - this->LastCameraFocalPoint[this->Axis1], 2 ) +
        pow( originAxis2 + half - this->LastCameraFocalPoint[this->Axis2], 2 ) ) <
        // pow( this->Radius + half * sqrt(2.), 2 );
        pow( this->Radius + half * 1.414213562, 2 );
    } else {
      insideBB = true;
    }

    if ( insideBB && this->BBSelection )
    {
      //JB On determine si la maille correspondant au current node of the tree
      //JB is going to be renedered.
      //JB Pour cela, on verifie si la maille est dans une boite englobante correspondant a la
      //JB projection de l'ecran dans le monde du maillage.
      insideBB = ( ( originAxis1 + 2 * halfAxis1 >= this->WindowBounds[0] )
                && ( originAxis1 <= this->WindowBounds[1] )
                && ( originAxis2 + 2 * halfAxis2 >= this->WindowBounds[2] )
                && ( originAxis2 <= this->WindowBounds[3] ) );
#ifndef NDEBUG
      if ( ! insideBB )
      {
        this->NbRejectByBB ++;
      }
    } else {
        this->NbRejectByCircle ++;
#endif
    }
  }
  if( insideBB )
  {
    // We only process those nodes than are going to be rendered
    if ( cursor->IsLeaf() || ( this->LevelMax != -1 && level >= this->LevelMax ) )
    {
      if ( this->Dimension == 2 )
      {
        this->ProcessLeaf2D( cursor );
      }
      else
      {
        this->ProcessLeaf1D( cursor );
      } // else
    } // if ( cursor->IsLeaf() || ( this->LevelMax!=-1 && level >= this->LevelMax ) )
    else
    {
      // Cursor is not at leaf, recurse to all children
      int numChildren = cursor->GetNumberOfChildren();
      for ( int ichild = 0; ichild < numChildren; ++ ichild )
      {
        cursor->ToChild( ichild );
        // Recurse
        this->RecursivelyProcessTreeNot3D( cursor, level+1 );
        cursor->ToParent( );
      } // ichild
    } // else
  } // if( insideBB )
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::ProcessLeaf1D( vtkHyperTreeGridNonOrientedGeometryCursor* cursor )
{
  // In 1D the geometry is composed of edges, create storage for endpoint IDs
  vtkIdType id[2];

  // First endpoint is at origin of cursor
  const double* origin = cursor->GetOrigin();
  id[0] = this->Points->InsertNextPoint( origin );

  // Second endpoint is at origin of cursor plus its length
  double pt[3];
  memcpy( pt, origin, 3 * sizeof( double ) );
  switch ( this->Orientation )
  {
    case 3: // 1 + 2
      pt[2] += cursor->GetSize()[2];
      break;
    case 5: // 1 + 4
      pt[1] += cursor->GetSize()[1];
      break;
    case 6: // 2 + 4
      pt[0] += cursor->GetSize()[0];
      break;
  } // switch
  id[1] = this->Points->InsertNextPoint( pt );

  // Insert edge into 1D geometry
  this->Cells->InsertNextCell( 2, id );
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::ProcessLeaf2D( vtkHyperTreeGridNonOrientedGeometryCursor* cursor )

{
  // Cell at cursor center is a leaf, retrieve its global index
  vtkIdType id = cursor->GetGlobalNodeIndex();
  if ( id < 0 )
  {
    return;
  }

  // In 2D all unmasked faces are generated
  if ( ! this->MaterialMask  || ! this->MaterialMask->GetValue( id ) )
  {
    // Insert face into 2D geometry depending on orientation
    this->AddFace( id, cursor->GetOrigin(), cursor->GetSize(), 0, this->Orientation );
  }
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::RecursivelyProcessTree3D(
  vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight* cursor,
  int level )
{
  // Create geometry output if cursor is at leaf
  if ( cursor->IsLeaf() )
  {
    this->ProcessLeaf3D( cursor );
  } // if ( cursor->IsLeaf() )
  else
  {
    // Cursor is not at leaf, recurse to all children
    int numChildren = cursor->GetNumberOfChildren();
    for ( int ichild = 0; ichild < numChildren; ++ ichild )
    {
      cursor->ToChild( ichild );
      // Recurse
      this->RecursivelyProcessTree3D( cursor, level+1 );
      cursor->ToParent( );
    } // child
  } // else
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::ProcessLeaf3D(
  vtkHyperTreeGridNonOrientedVonNeumannSuperCursorLight* superCursor )
{
  // Cell at super cursor center is a leaf, retrieve its global index, level, and mask
  vtkIdType idcenter = superCursor->GetGlobalNodeIndex();
  unsigned level = superCursor->GetLevel();
  int masked = this->MaterialMask ? this->MaterialMask->GetValue( idcenter ) : 0;

  // Iterate over all cursors of Von Neumann neighborhood around center
  unsigned int nc = superCursor->GetNumberOfCursors() - 1;
  for ( unsigned int c = 0 ; c < nc; ++ c )
  {
    // Retrieve cursor to neighbor across face
    // Retrieve tree, leaf flag, and mask of neighbor cursor
    unsigned int levelN;
    bool leafN;
    vtkIdType idN;
    vtkHyperTree* treeN = superCursor->GetInformation( VonNeumannCursors3D[c], levelN, leafN, idN );

    int maskedN = 1;
    if ( treeN )
    {
      maskedN = this->MaterialMask ? this->MaterialMask->GetValue( idN ) : 0;
    }

    // In 3D masked and unmasked cells are handled differently:
    // . If cell is unmasked, and face neighbor is a masked leaf, or no such neighbor
    //   exists, then generate face.
    // . If cell is masked, and face neighbor exists and is an unmasked leaf, then
    //   generate face, breaking ties at same level. This ensures that faces between
    //   unmasked and masked cells will be generated once and only once.
    if ( ( ! masked && ( ! treeN || ( leafN && maskedN ) ) )
         ||
         ( masked && treeN && leafN && levelN < level && ! maskedN ) )
    {
      // Generate face with corresponding normal and offset
      this->AddFace( idcenter, superCursor->GetOrigin(), superCursor->GetSize(),
                     VonNeumannOffsets3D[c], VonNeumannOrientations3D[c] );
    }
  } // c
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::AddFace( vtkIdType inId,
                                               const double* origin,
                                               const double* size,
                                               int offset,
                                               unsigned int orientation )
{
  // Storage for point coordinates
  double pt[] = { 0., 0., 0. };

  // Storage for face vertex IDs
  vtkIdType ids[4];

  // First cell vertex is always at origin of cursor
  memcpy( pt, origin, 3 * sizeof( double ) );

  if ( this->Locator )
  {
    if ( offset )
    {
      // Offset point coordinate as needed
      pt[orientation] += size[orientation];
    }
    this->Locator->InsertUniquePoint( pt, ids[0]);
    // Create other face vertices depending on orientation
    unsigned int axis1 = orientation ? 0 : 1;
    unsigned int axis2 = orientation == 2 ? 1 : 2;
    pt[axis1] += size[axis1];
    this->Locator->InsertUniquePoint( pt, ids[1] );
    pt[axis2] += size[axis2];
    this->Locator->InsertUniquePoint( pt, ids[2] );
    pt[axis1] = origin[axis1];
    this->Locator->InsertUniquePoint( pt, ids[3] );
  } else {
    if ( offset )
    {
      // Offset point coordinate as needed
      pt[orientation] += size[orientation];
    }
    ids[0] = this->Points->InsertNextPoint( pt );

    // Create other face vertices depending on orientation
    unsigned int axis1 = orientation ? 0 : 1;
    unsigned int axis2 = orientation == 2 ? 1 : 2;
    pt[axis1] += size[axis1];
    ids[1] = this->Points->InsertNextPoint( pt );
    pt[axis2] += size[axis2];
    ids[2] = this->Points->InsertNextPoint( pt );
    pt[axis1] = origin[axis1];
    ids[3] = this->Points->InsertNextPoint( pt );
  }

  // Insert next face
  vtkIdType outId = this->Cells->InsertNextCell( 4, ids );

  // Copy face data from that of the cell from which it comes
  this->OutData->CopyData( this->InData, inId, outId );
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::SetRenderer( vtkRenderer* ren )
{
  if ( ren != this->Renderer )
  {
    this->Renderer = ren;
    this->Modified();
  }
}

//----------------------------------------------------------------------------
vtkMTimeType vtkAdaptiveDataSetSurfaceFilter::GetMTime()
{
  // Check for minimal changes
  if ( this->Renderer )
  {
    vtkCamera* cam = this->Renderer->GetActiveCamera();
    if ( cam )
    {
      // Check & Update parallel projection
      bool para = (cam->GetParallelProjection()!=0);

      if ( this->ParallelProjection != para )
      {
        this->ParallelProjection = para;
        this->Modified();
      }

      // Check & Update renderer size
      int* sz = this->Renderer->GetSize();
      if (   this->LastRendererSize[0] != sz[0]
                || this->LastRendererSize[1] != sz[1] )
      {
        this->LastRendererSize[0] = sz[0];
        this->LastRendererSize[1] = sz[1];
        this->Modified();
      }

      // Check & Update camera focal point
      double* fp = cam->GetFocalPoint();
      if ( this->LastCameraFocalPoint[0] != fp[0]
        || this->LastCameraFocalPoint[1] != fp[1]
        || this->LastCameraFocalPoint[2] != fp[2] )
      {
        this->LastCameraFocalPoint[0] = fp[0];
        this->LastCameraFocalPoint[1] = fp[1];
        this->LastCameraFocalPoint[2] = fp[2];
        this->Modified();
      }

      // Check & Update camera scale
      double scale = cam->GetParallelScale();
      if( this->LastCameraParallelScale != scale)
      {
        this->LastCameraParallelScale = scale;
        this->Modified();
      }
    } // if ( cam )
  } // if ( this->Renderer )
  return this->Superclass::GetMTime();
}
