text: /*=========================================================================
  
    Program:   Visualization Toolkit
    Module:    $RCSfile: vtkHyperOctreeContourFilter.cxx,v $
  
    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 "vtkHyperOctreeContourFilter.h"
  
  #include "vtkHyperOctree.h"
  #include "vtkCellArray.h"
  #include "vtkCellData.h"
  #include "vtkClipVolume.h"
  #include "vtkExecutive.h"
  #include "vtkDoubleArray.h"
  #include "vtkGenericCell.h"
  #include "vtkImageData.h"
  #include "vtkImplicitFunction.h"
  #include "vtkInformation.h"
  #include "vtkInformationVector.h"
  #include "vtkIntArray.h"
  #include "vtkMergePoints.h"
  #include "vtkObjectFactory.h"
  #include "vtkPointData.h"
  #include "vtkUnsignedCharArray.h"
  #include "vtkUnstructuredGrid.h"
  #include "vtkHyperOctreeCursor.h"
  #include "vtkVoxel.h"
  #include "vtkPixel.h"
  #include "vtkLine.h"
  #include "vtkTetra.h"
  #include "vtkPolygon.h"
  #include "vtkOrderedTriangulator.h"
  #include "vtkStreamingDemandDrivenPipeline.h"
  #include <math.h>
  #include <assert.h>
  #include <vtkstd/set>
  #include "vtkBitArray.h"
  
  #include "vtkHyperOctreePointsGrabber.h"
  
  
  //----------------------------------------------------------------------------
  void vtkHyperOctreeContourFilter::PrintSelf(ostream& os, vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
    
    this->ContourValues->PrintSelf(os,indent.GetNextIndent());
  
    if ( this->Locator )
      {
      os << indent << "Locator: " << this->Locator << "\n";
      }
    else
      {
      os << indent << "Locator: (none)\n";
      }
  }
  
  
  class vtkHyperOctreeIdSet // Pimpl idiom
  {
  public:
    vtkstd::set<vtkIdType> Set;
  };
  
  class vtkHyperOctreeContourPointsGrabber; 
  
  class vtkHyperOctreeContourPointsGrabber : public vtkHyperOctreePointsGrabber
  {
  public:
    static vtkHyperOctreeContourPointsGrabber *New();
    
    vtkTypeRevisionMacro(vtkHyperOctreeContourPointsGrabber,vtkHyperOctreePointsGrabber);
    
    void PrintSelf(ostream& os, vtkIndent indent);
    
    // Description:
    // Set the dimension of the hyperoctree.
    // \pre valid_dim: (dim==2 || dim==3)
    // \post is_set: GetDimension()==dim
    virtual void SetDimension(int dim);
    
    // Description:
    // Initialize the points insertion scheme.
    // Actually, it is just a trick to initialize the IdSet from the filter.
    // The IdSet class cannot be shared with the filter because it is a Pimpl.
    // It is used by clip,cut and contour filters to build the points
    // that lie on an hyperoctant.
    // \pre only_in_3d: GetDimension()==3
    virtual void InitPointInsertion();
    
    // Description:
    // Insert a point, assuming the point is unique and does not require a
    // locator. Tt does not mean it does not use a locator. It just mean that
    // some implementation may skip the use of a locator.
    virtual void InsertPoint(vtkIdType ptId,
                             double pt[3],
                             double pcoords[3],
                             int ijk[3]);
    
    // Description:
    // Insert a point using a locator.
    virtual void InsertPointWithMerge(vtkIdType ptId,
                                      double pt[3],
                                      double pcoords[3],
                                      int ijk[3]);
    
    // Description:
    // Insert a point in the quadtree case.
    virtual void InsertPoint2D(double pt[3],
                               int ijk[3]);
    
    // Description:
    // Return the ordered triangulator.
    vtkOrderedTriangulator *GetTriangulator();
    
    // Description:
    // Return the polygon.
    vtkPolygon *GetPolygon();
    
    // Description:
    // Init the locator.
    void InitLocator(vtkPoints *pts,
                     double bounds[6]);
    
    // Description:
    // Return the id in the locator after a call to InsertPoint*().
    vtkIdType GetLastPtId();
    
    void SetFilter(vtkHyperOctreeContourFilter *filter);
    
  protected:
    // Constructor with default bounds (0,1, 0,1, 0,1).
    vtkHyperOctreeContourPointsGrabber();
    ~vtkHyperOctreeContourPointsGrabber();
    
    vtkHyperOctreeContourFilter *Filter;
    vtkOrderedTriangulator *Triangulator;
    vtkPolygon *Polygon;
    vtkMergePoints *Locator;
    vtkIdType LastPtId;
    vtkHyperOctreeIdSet *IdSet;
  private:
    vtkHyperOctreeContourPointsGrabber(const vtkHyperOctreeContourPointsGrabber&);  // Not implemented.
    void operator=(const vtkHyperOctreeContourPointsGrabber&);    // Not implemented.
  };
    
  vtkCxxRevisionMacro(vtkHyperOctreeContourFilter, "$Revision: 1.4 $");
  vtkStandardNewMacro(vtkHyperOctreeContourFilter);
  
  //----------------------------------------------------------------------------
  // Construct with user-specified implicit function; value
  // set to 0.0; and generate cut scalars turned off.
  vtkHyperOctreeContourFilter::vtkHyperOctreeContourFilter()
  {
    this->ContourValues = vtkContourValues::New();
    
    this->Locator = NULL;
    
    this->SetNumberOfOutputPorts(1);
  
    // by default process active cell scalars
    // This is points because octree returns dual grid.
    this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
                                 vtkDataSetAttributes::SCALARS);
    
    this->Input=0;
    this->Output=0;
   
    this->NewVerts=0;
    this->NewLines=0;
    this->NewPolys=0;
    
    this->InCD=0;
    this->OutCD=0;
    this->OutPD=0;
    this->Triangulator=0;
    
    this->Tetra=0;
    this->TetScalars=0;
    this->PointScalars=0;
    this->CellScalars=0;
    
    this->Cursor=0;
    this->NeighborCursor=0;
    this->Sibling=0;
    
    this->InScalars=0;
    this->Grabber=0;
    this->Polygon=0;
    this->SortBy=VTK_SORT_BY_VALUE;
    
    this->Line=0;
  }
  
  //----------------------------------------------------------------------------
  vtkHyperOctreeContourFilter::~vtkHyperOctreeContourFilter()
  {
    this->ContourValues->Delete();
    if ( this->Locator )
      {
      this->Locator->UnRegister(this);
      this->Locator = NULL;
      }
  }
  
  //----------------------------------------------------------------------------
  // Overload standard modified time function. If Cut function is modified,
  // then this object is modified as well.
  unsigned long vtkHyperOctreeContourFilter::GetMTime()
  {
    unsigned long mTime=this->Superclass::GetMTime();
    unsigned long contourValuesMTime=this->ContourValues->GetMTime();
    unsigned long time;
    
    mTime = ( contourValuesMTime > mTime ? contourValuesMTime : mTime );
    
    if ( this->Locator != NULL )
      {
      time = this->Locator->GetMTime();
      mTime = ( time > mTime ? time : mTime );
      }
  
    return mTime;
  }
  
  //----------------------------------------------------------------------------
  //
  // Cut through data generating surface.
  //
  int vtkHyperOctreeContourFilter::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
    this->Input = vtkHyperOctree::SafeDownCast(
      inInfo->Get(vtkDataObject::DATA_OBJECT()));
    
    if(this->Input->GetNumberOfLevels()==1)
      {
      // just the root. There is absolutely no chance
      // to get an isosurface here.
      this->Input=0;
      return 1;
      }
    
    this->InScalars=this->GetInputArrayToProcess(0,inputVector);
    if(this->InScalars==0)
      {
      vtkDebugMacro(<<"No data to contour");
      this->Input=0;
      return 1;
      }
    
    int numContours=this->ContourValues->GetNumberOfContours();
    if(numContours==0)
      {
      vtkDebugMacro(<<"No contour");
      this->Input=0;
      return 1;
      }
    
    double *values=this->ContourValues->GetValues();
    
    // If all the contour values are out of the range of the input scalar
    // there is no chance to get a contour, just exit.
   
    double range[2];
    this->InScalars->GetRange(range);
    int i=0;
    int allOut=1;
    while(allOut && i<numContours)
      {
      allOut=(values[i]<range[0]) || (values[i]>range[1]);
      ++i;
      }
    if(allOut)
      {
      // empty output
      this->Input=0;
      return 1;
      }
    
    
    this->Output=vtkPolyData::SafeDownCast(
      outInfo->Get(vtkDataObject::DATA_OBJECT()));
    
    vtkIdType numPts=this->Input->GetMaxNumberOfPoints(0);
    vtkIdType numCells = this->Input->GetNumberOfLeaves();
    
    vtkPoints *newPoints = vtkPoints::New();
    newPoints->Allocate(numPts,numPts/2);
    
    // allocate the output and associated helper classes
    vtkIdType estimatedSize = numCells;
    estimatedSize = estimatedSize / 1024 * 1024; //multiple of 1024
    if (estimatedSize < 1024)
      {
      estimatedSize = 1024;
      }
    
   
    this->NewVerts=vtkCellArray::New();
    this->NewVerts->Allocate(estimatedSize,estimatedSize/2);
    this->NewLines = vtkCellArray::New();
    this->NewLines->Allocate(estimatedSize,estimatedSize/2);
    this->NewPolys = vtkCellArray::New();
    this->NewPolys->Allocate(estimatedSize,estimatedSize/2);
    
    // locator used to merge potentially duplicate points
    if ( this->Locator == NULL )
      {
      this->CreateDefaultLocator();
      }
    
    this->Locator->InitPointInsertion (newPoints, this->Input->GetBounds());
    
   
    this->InCD=(vtkCellData*)(this->Input->GetLeafData());
    // Scalars are added to this, so we need to make a copy.
    this->InPD = vtkPointData::New();
    // Since the dataset API returns the dual, cell and point data are switched.
    this->InPD->ShallowCopy(this->Input->GetCellData());
    
    this->OutCD=this->Output->GetCellData();
    this->OutCD->CopyAllocate(this->InCD,estimatedSize,estimatedSize/2);
      
    this->OutPD=this->Output->GetPointData();
    this->OutPD->CopyAllocate(this->Input->GetPointData(),estimatedSize,estimatedSize/2);
    this->OutPD->CopyScalarsOn();
    
    static double bounds[6]={0,1,0,1,0,1};
    
    vtkPoints *originalPoints=0;
    
    switch(this->Input->GetDimension())
      {
      case 3:
        this->Tetra=vtkTetra::New();
        this->TetScalars=vtkDoubleArray::New();
        this->TetScalars->SetNumberOfComponents(1);
        this->TetScalars->SetNumberOfTuples(4);
        this->Grabber=vtkHyperOctreeContourPointsGrabber::New();
        this->Grabber->SetFilter(this);
        this->Grabber->SetDimension(3);
        originalPoints=vtkPoints::New();
        originalPoints->Allocate(numPts,numPts/2);
        this->Grabber->InitLocator(originalPoints,bounds);
        this->Triangulator=this->Grabber->GetTriangulator();
        break;
      case 2:
        this->Grabber=vtkHyperOctreeContourPointsGrabber::New();
        this->Grabber->SetFilter(this);
        this->Grabber->SetDimension(2);
        originalPoints=vtkPoints::New();
        originalPoints->Allocate(numPts,numPts/2);
        
        this->Grabber->InitLocator(originalPoints,this->Input->GetBounds());
        this->Polygon=this->Grabber->GetPolygon();
        break;
      case 1:
        this->Line=vtkLine::New();
        break;
      default:
        assert("check: impossible case" && 0); // do nothing
        break;
      }
    this->PointScalars=vtkDoubleArray::New();
    this->PointScalars->SetName(this->InScalars->GetName());
    this->PointScalars->Allocate(estimatedSize);
    
    this->TotalCounter=0;
    this->TemplateCounter=0;
    
    int j=0;
    while(j<65536)
      {
      this->CellTypeCounter[j]=0; // up-to-65536 points per octant
      ++j;
      }
    
    this->Cursor=this->Input->NewCellCursor();
    this->NeighborCursor=this->Input->NewCellCursor();
    this->Sibling=this->Input->NewCellCursor();
    
    this->PointScalars->SetNumberOfComponents(1);
    this->PointScalars->Allocate(estimatedSize);
    this->Cursor->ToRoot();
    
    
    this->CellScalars=vtkDoubleArray::New();
    
    // let's go
    if(this->Input->GetDimension()==1)
      {
      // got the first leaf on the left side.
      this->Cursor->ToRoot();
      while(!this->Cursor->CurrentIsLeaf())
        {
        this->Cursor->ToChild(0);
        }
      this->LeftValue=this->InScalars->GetTuple1(this->Cursor->GetLeafId());
      this->LeftCoord=this->Input->GetOrigin()[0];
      this->Cursor->ToRoot();
      this->ContourNode1D();
      }
    else
      {
      this->ContourNode();
      }
    
    this->CellScalars->UnRegister(this);
    this->CellScalars=0;
    
    this->PointScalars->UnRegister(this);
    this->PointScalars=0;
    
    this->InCD=0;
    
    this->Cursor->UnRegister(this);
    this->Cursor=0;
    this->NeighborCursor->UnRegister(this);
    this->NeighborCursor=0;
    this->Sibling->UnRegister(this);
    this->Sibling=0;
   
  //  cout<<"ClipHyperOctree: "<<this->TemplateCounter<<" templates over "<<this->TotalCounter<<" octants, ratio="<<(this->TemplateCounter/static_cast<double>(this->TotalCounter))<<endl;
    
    j=0;
    while(j<65536)
      {
      if(this->CellTypeCounter[j]>0)
        {
  //      cout<<this->CellTypeCounter[j]<<" with "<<j+1<<"points"<<endl;
        }
      ++j;
      }
    
    switch(this->Input->GetDimension())
      {
      case 3:
        this->Tetra->UnRegister(this);
        this->Tetra=0;
        this->TetScalars->UnRegister(this);
        this->TetScalars=0;
        this->Triangulator=0;
        this->Grabber->UnRegister(this);
        this->Grabber=0;
        originalPoints->UnRegister(this);
        break;
      case 2:
        this->Polygon=0;
        this->Grabber->UnRegister(this);
        this->Grabber=0;
        originalPoints->UnRegister(this);
        break;
      case 1:
        this->Line->UnRegister(this);
        this->Line=0;
        break;
      default:
        assert("check: impossible case" && 0);
        break;
      }
    
    this->OutPD=0;
    this->Input=0;
    
    this->Output->SetPoints(newPoints);
    newPoints->Delete();
    
    if (this->NewVerts->GetNumberOfCells()>0)
      {
      this->Output->SetVerts(this->NewVerts);
      }
    this->NewVerts->Delete();
    this->NewVerts=0;
  
    if (this->NewLines->GetNumberOfCells()>0)
      {
      this->Output->SetLines(this->NewLines);
      }
    this->NewLines->Delete();
    this->NewLines=0;
    
    if (this->NewPolys->GetNumberOfCells()>0)
      {
      this->Output->SetPolys(this->NewPolys);
      }
    this->NewPolys->Delete();
    this->NewPolys=0;
    
    this->OutCD=0;
    this->InPD->Delete();
    this->InPD = 0;
    
    this->Locator->Initialize();//release any extra memory
    this->Output->Squeeze();
    this->Output=0;
    
    assert("post: input_is_null" && this->Input==0);
    assert("post: output_is_null" && this->Output==0);
    assert("post: incd_is_null" && this->InCD==0);
    assert("post: outpd_is_null" && this->OutPD==0);
    assert("post: outcd_is_null" && this->OutCD==0);
    
    return 1;
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperOctreeContourFilter::ContourNode1D()
  {
    if(!this->Cursor->CurrentIsLeaf())
      {
      int child=0;
      while(child<2)
        {
        this->Cursor->ToChild(child);
        this->ContourNode1D();
        this->Cursor->ToParent();
        ++child;
        }
      }
    else
      {
      vtkIdType cellId=this->Cursor->GetLeafId();
      double cellValue=this->InScalars->GetTuple1(cellId);
      int level=this->Cursor->GetCurrentLevel();
      
      int target[3];
      target[1]=0;
      target[2]=0;
      target[0]=this->Cursor->GetIndex(0);
      double rightValue;
    
      if(target[0]>=(1<<(level-1)))
        {
        rightValue=cellValue;
        }
      else
        {
        this->NeighborCursor->MoveToNode(target,level);
        rightValue=(cellValue+this->InScalars->GetTuple1(this->NeighborCursor->GetLeafId()))*0.5;
        }
      
      // build a line and contour it.
      double pt[3];
      pt[1]=this->Input->GetOrigin()[1];
      pt[2]=this->Input->GetOrigin()[2];
      // build the point
      pt[0]=this->LeftCoord;
      this->Line->GetPoints()->SetPoint(0,pt);
      pt[0]+=this->Input->GetSize()[0]/(1<<level);
      this->LeftCoord=pt[0]; // for the next step
      this->Line->GetPoints()->SetPoint(1,pt);
      
      this->CellScalars->SetNumberOfComponents(1);
      this->CellScalars->SetNumberOfTuples(2);
      
      this->CellScalars->SetValue(0,this->LeftValue);
      this->CellScalars->SetValue(1,rightValue);
        
      // It was a bad idea to modify the input PointData.
      // We made a copy.
      //vtkPointData *inPD=this->Input->GetPointData(); // void
      this->InPD->SetScalars(this->PointScalars); // void
      
      if ( this->SortBy == VTK_SORT_BY_CELL )
        {
        double value = this->ContourValues->GetValue(this->Iter);
        this->Line->Contour(value, this->CellScalars, this->Locator, 
                            this->NewVerts,this->NewLines,this->NewPolys,
                            this->InPD,this->OutPD,this->InCD,cellId,
                            this->OutCD);
        }
      else
        {
        // VTK_SORT_BY_VALUE
        int iter=0;
        int numContours=this->ContourValues->GetNumberOfContours();
        while(iter<numContours)
          {
          double value = this->ContourValues->GetValue(iter);
          this->Line->Contour(value, this->CellScalars, this->Locator, 
                              this->NewVerts,this->NewLines,
                              this->NewPolys,this->InPD,this->OutPD,
                              this->InCD,cellId,this->OutCD);
            
            ++iter;
            }
          }
      
      
      // initialize the left value for the next leaf.
      this->LeftValue=rightValue;
      }
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // Do the recursive contour of the node pointed by Cursor.
  void vtkHyperOctreeContourFilter::ContourNode()
  {
    if(!this->Cursor->CurrentIsLeaf())
      {
      int child=0;
      int c=this->Cursor->GetNumberOfChildren();
      while(child<c)
        {
        this->Cursor->ToChild(child);
        this->ContourNode();
        this->Cursor->ToParent();
        ++child;
        }
      }
    else
      {
      // some parent=>have sibling=>some sibling may have children
      // => those children may create points on some face of cursor
      // => difficult case
      //
      // Even worst, if the siblings don't have children, the
      // sibling of the parent may have children that create points
      // on some face.
      //
      // Even if there is no children, the neighbor cell tessellation
      // has to be compatible with the current cell tessellation.
      // In any case, we need the ordered triangulator.
      
        
      // Add the points of the current leaf
      // use the bounds
      int i;
      int j;
      int k;
      
      int x;
      int y;
      int z;
      
      // resolution in point along each axis.
      int resolution=(1<<(this->Input->GetNumberOfLevels()-1))+1;
  
      int level=this->Cursor->GetCurrentLevel();
      int deltaLevel=this->Input->GetNumberOfLevels()-1-level;
      assert("check: positive_deltaLevel" && deltaLevel>=0);
      
      double ratio=1.0/(resolution-1);
      
      double pt[3];
      double pcoords[3];
      int pijk[3];
      int qijk[3];
      
      int numContours=0; // initialized to removed warnings
      
      if ( this->SortBy == VTK_SORT_BY_VALUE )
        {
        numContours=this->ContourValues->GetNumberOfContours();
        }
      // index of the node
      if(this->Input->GetDimension()==3)
        {
        vtkIdType nbpts=this->Input->GetMaxNumberOfPointsOnBoundary(level);
        double pbounds[6]={0,1,0,1,0,1};
        double *size=this->Input->GetSize();
        double *origin=this->Input->GetOrigin();
        
        this->Triangulator->InitTriangulation(pbounds,nbpts);
        
        this->Triangulator->PreSortedOff();
        this->Grabber->InitPointInsertion();
        
        i=this->Cursor->GetIndex(0);
        j=this->Cursor->GetIndex(1);
        k=this->Cursor->GetIndex(2);
        z=0;
        pijk[2]=k;
        while(z<2)
          {
          pijk[1]=j;
          y=0;
          while(y<2)
            {
            pijk[0]=i;
            x=0;
            while(x<2)
              {
              // Get some parametric coords in [0,1]
              // [0,1] covers the whole dataset axis.
              int coord=0;
              while(coord<3)
                {
                qijk[coord]=pijk[coord]<<deltaLevel;
                pcoords[coord]=qijk[coord]*ratio;
                pt[coord]=pcoords[coord]*size[coord]+origin[coord];
                ++coord;
                }
              
              assert("check: in_bounds" && pt[0]>=this->Input->GetBounds()[0] && pt[0]<=this->Input->GetBounds()[1] && pt[1]>=this->Input->GetBounds()[2] && pt[1]<=this->Input->GetBounds()[3] && pt[2]>=this->Input->GetBounds()[4] && pt[2]<=this->Input->GetBounds()[5]);
              
              this->Grabber->InsertPoint(0,pt,pcoords,qijk);
              ++x;
              ++pijk[0];
              }
            ++y;
            ++pijk[1];
            }
          ++z;
          ++pijk[2];
          }
        }
      
      int lastLevelLeaf=level>=(this->Input->GetNumberOfLevels()-1);
      
      if(this->Input->GetDimension()==3)
        {
        if(!lastLevelLeaf)
          {
          // Ok, now ask my parent if I have sibling with children on my
          // faces and even worst, if my parent has sibling with children
          // that have children on my face, or if the parent of my parent
          // has sibling with children that have children, that have children
          // on my face, until I reach the root...
          
          // list the 3 faces of the parent, the current node is lying on.
          int faces[3];
          
          int child=this->Cursor->GetChildIndex();
          
          faces[0]=(child&1)==1; // false: -x, true: +x
          faces[1]=(child&2)==2; // false: -y, true: +y
          faces[2]=(child&4)==4; // false: -z, true: +z
          
          // sibling on faces that are not on a parent face
          int siblings[3];
          
          int inc=1;
          i=0;
          while(i<3)
            {
            if(faces[i])
              {
              siblings[i]=child-inc;
              }
            else
              {
              siblings[i]=child+inc;
              }
            ++i;
            inc<<=1;
            }
          
          this->Sibling->ToSameNode(this->Cursor);
          this->Sibling->ToParent();
          // ask the 3 sibling, one on each face of the current node
          i=0; 
          int faceOffset=0;
          while(i<3)
            {
            this->Sibling->ToChild(siblings[i]);
            assert("check: we are not visiting ourselves" && this->Sibling->GetChildIndex()!=child);
            if(!this->Sibling->CurrentIsLeaf())
              {
              assert("check: if the sibling is not a leaf we cannot be at the last level" && level<(this->Input->GetNumberOfLevels()-1));
              
              // get the points of this sibling on some given face
              int siblingFace=faceOffset;
              if(faces[i])
                {
                ++siblingFace;
                }
              this->Input->GetPointsOnFace(this->Sibling,siblingFace,level,
                                           this->Grabber);
              }
            this->Sibling->ToParent();
            ++i;
            faceOffset+=2;
            }
          
          // Get points on faces shared with the parent node.
          this->Input->GetPointsOnParentFaces(faces,level,this->Cursor,
                                              this->Grabber);
          
          // Get the points from the edge-only neighbors.
          int childIndices[3];
          childIndices[2]=(child&4)>>2;
          childIndices[1]=(child&2)>>1;
          childIndices[0]=child&1;
          
          assert("check valid_range_c2" && childIndices[2]>=0 &&
                 childIndices[2]<=1);
          assert("check valid_range_c1" && childIndices[1]>=0 &&
                 childIndices[1]<=1);
          assert("check valid_range_c0" && childIndices[0]>=0 &&
                 childIndices[0]<=1);
          
          // First the edges aligned on X axis
          int axis=0;
          int a=2;
          int b=1;
          
          this->Sibling->ToSameNode(this->Cursor);
          this->Sibling->ToParent();
          
          while(axis<3)
            {
            k=0;
            while(k<2)
              {
              j=0;
              while(j<2)
                {
                if(k!=childIndices[a]&&j!=childIndices[b])
                  {
                  this->Sibling->ToChild((k<<a)+(j<<b)+
                                         (childIndices[axis]<<axis));
                  if(!this->Sibling->CurrentIsLeaf())
                    {
                    this->Input->GetPointsOnEdge(this->Sibling,level,axis,!k,
                                                 !j,this->Grabber);
                    }
                  this->Sibling->ToParent();
                  }
                else
                  {
                  this->Input->GetPointsOnParentEdge(this->Cursor,level,axis,k,j,
                                                     this->Grabber);
                  }
                ++j;
                }
              ++k;
              }
            ++axis;
            ++a;
            if(a>2)
              {
              a-=3;
              }
            ++b;
            if(b>2)
              {
              b-=3;
              }
            }
          } // if not leaf at last level
        }
      else
        {
        // this->Input->GetDimension()==2
        // counter- clockwise direction matters here.
        
        int edges[2];
        int child=this->Cursor->GetChildIndex();
        double *size=this->Input->GetSize();
        double *origin=this->Input->GetOrigin();
        
        this->Polygon->GetPointIds()->SetNumberOfIds(0);
        this->Polygon->GetPoints()->SetNumberOfPoints(0);
        
        if(!lastLevelLeaf)
          {
          this->Sibling->ToSameNode(this->Cursor);
          this->Sibling->ToParent();
          
          // list the 2 edges of the parent, the current node is lying on.
          edges[0]=(child&1)==1; // false: -x, true: +x
          edges[1]=(child&2)==2; // false: -y, true: +y
          }
        else
          {
          edges[0]=0;
          edges[1]=0;
          }
        
        i=this->Cursor->GetIndex(0);
        j=this->Cursor->GetIndex(1);
      
        
        // Insert vertex (xmin,ymin)
        pijk[0]=i;
        pijk[1]=j;
        
        qijk[0]=pijk[0]<<deltaLevel;
        qijk[1]=pijk[1]<<deltaLevel;
        
        pcoords[0]=qijk[0]*ratio;
        pcoords[1]=qijk[1]*ratio;
        
        pt[0]=pcoords[0]*size[0]+origin[0];
        pt[1]=pcoords[1]*size[1]+origin[1];
        pt[2]=origin[2];
        
        this->Grabber->InsertPoint2D(pt,qijk);
        
        if(!lastLevelLeaf)
          {
          // Process edge (-y)
          if(edges[1])
            {
            // sibling
            this->Sibling->ToChild(child-2);
            if(!this->Sibling->CurrentIsLeaf())
              {
              this->Input->GetPointsOnEdge2D(this->Sibling,3,level,
                                             this->Grabber); // 3==+y
              }
            this->Sibling->ToParent();
            }
          else
            {
            // parent
            this->Input->GetPointsOnParentEdge2D(this->Cursor,2,level,
                                                 this->Grabber); // 2==-y
            }
          }
        
        // Insert vertex (xmax,ymin)
        pijk[0]=i+1;
        qijk[0]=pijk[0]<<deltaLevel;
        pcoords[0]=qijk[0]*ratio;
        pt[0]=pcoords[0]*size[0]+origin[0];
        
        this->Grabber->InsertPoint2D(pt,qijk);
        
        if(!lastLevelLeaf)
          {
          // Process edge (+x)
          if(edges[0])
            {
            // parent
            this->Input->GetPointsOnParentEdge2D(this->Cursor,1,level,
                                                 this->Grabber); // 1==+x
            }
          else
            {
            // sibling
            this->Sibling->ToChild(child+1);
            if(!this->Sibling->CurrentIsLeaf())
              {
              this->Input->GetPointsOnEdge2D(this->Sibling,0,level,
                                             this->Grabber); //0==-x
              }
            this->Sibling->ToParent();
            }
          }
        
        // Insert vertex (xmax,ymax)
        pijk[1]=j+1;
        qijk[1]=pijk[1]<<deltaLevel;
        pcoords[1]=qijk[1]*ratio;
        pt[1]=pcoords[1]*size[1]+origin[1];
        this->Grabber->InsertPoint2D(pt,qijk);
      
        if(!lastLevelLeaf)
          {
          // Process edge (+y)
          if(edges[1])
            {
            // parent
            this->Input->GetPointsOnParentEdge2D(this->Cursor,3,level,
                                                 this->Grabber); // 3==+y
            }
          else
            {
            // sibling
            this->Sibling->ToChild(child+2);
            if(!this->Sibling->CurrentIsLeaf())
              {
              this->Input->GetPointsOnEdge2D(this->Sibling,2,level,
                                             this->Grabber); //2==-y
              }
            this->Sibling->ToParent();
            }
          }
        
        
        // Insert vertex (xmin,ymax)
        pijk[0]=i;
        qijk[0]=pijk[0]<<deltaLevel;
        pcoords[0]=qijk[0]*ratio;
        pt[0]=pcoords[0]*size[0]+origin[0];
        
        this->Grabber->InsertPoint2D(pt,qijk);
        
        if(!lastLevelLeaf)
          {
          // Process edge (-x)
          if(edges[0])
            {
            // sibling
            this->Sibling->ToChild(child-1);
            if(!this->Sibling->CurrentIsLeaf())
              {
              this->Input->GetPointsOnEdge2D(this->Sibling,1,level,
                                             this->Grabber); // 1==+x
              }
            this->Sibling->ToParent();
            }
          else
            {
            // parent
            this->Input->GetPointsOnParentEdge2D(this->Cursor,0,level,
                                                 this->Grabber); // 0==-x
            }
          }
        }
      
      if(this->Input->GetDimension()==3)
        {     
        int c=this->Triangulator->GetNumberOfPoints();
        
        
        this->CellScalars->SetNumberOfComponents(1);
        this->CellScalars->SetNumberOfTuples(c);
        
        
        // very important: we have to build the scalars value on the cell
        // BEFORE a call to Triangulate().
        i=0;
        while(i<c)
          {
          vtkIdType ptId=this->Triangulator->GetPointId(i);
          this->CellScalars->InsertValue(i,this->PointScalars->GetValue(ptId));
          ++i;
          }
        
        if(c==8)
          {
          // only the vertices of a voxel: fast path.
          this->Triangulator->UseTemplatesOn();
          this->Triangulator->TemplateTriangulate(VTK_VOXEL,8,12);
          ++this->TotalCounter;
          ++this->TemplateCounter;
          }
        else
          {
          // slow path
          this->Triangulator->UseTemplatesOff();
          this->Triangulator->Triangulate();
          ++this->TotalCounter;
          if(this->Triangulator->GetNumberOfPoints()<=65536)
            {
            this->CellTypeCounter[this->Triangulator->GetNumberOfPoints()-1]++;
            }
          }
        
        // perform contour
        
        vtkIdType cellId=this->Cursor->GetLeafId();
             
        // I made a copy of the input point data so it is ok modify.
        this->InPD->SetScalars(this->PointScalars); // otherwise, it crashes
         
        if ( this->SortBy == VTK_SORT_BY_CELL )
          {
          double value = this->ContourValues->GetValue(this->Iter);
          this->Triangulator->InitTetraTraversal();
          int numTetras=0; // debug
          int done=this->Triangulator->GetNextTetra(0,this->Tetra,
                                                    this->CellScalars,
                                                    this->TetScalars)==0;
          while(!done)
            {
            ++numTetras;
            this->Tetra->Contour(value, this->TetScalars, this->Locator, 
                                 this->NewVerts,this->NewLines,this->NewPolys,
                                 this->InPD,this->OutPD,this->InCD,cellId,
                                 this->OutCD);
            done=this->Triangulator->GetNextTetra(0,this->Tetra,
                                                  this->CellScalars,
                                                  this->TetScalars)==0;
            } //while
          }
        else
          {
          // VTK_SORT_BY_VALUE
          this->Triangulator->InitTetraTraversal();
          int numTetras=0; // debug
          int done=this->Triangulator->GetNextTetra(0,this->Tetra,
                                                    this->CellScalars,
                                                    this->TetScalars)==0;
          while(!done)
            {
            ++numTetras;
            int iter=0;
            while(iter<numContours)
              {
              double value = this->ContourValues->GetValue(iter);
              this->Tetra->Contour(value, this->TetScalars, this->Locator, 
                                   this->NewVerts,this->NewLines,
                                   this->NewPolys,this->InPD,this->OutPD,
                                   this->InCD,cellId,this->OutCD);
              ++iter;
              }
            done=this->Triangulator->GetNextTetra(0,this->Tetra,
                                                  this->CellScalars,
                                                  this->TetScalars)==0;
            } //while
          }
        }
      else
        {
        // this->Input->GetDimension()==2
       
        // perform contour
        
        vtkIdType cellId=this->Cursor->GetLeafId();
        
        int c=this->Polygon->GetPoints()->GetNumberOfPoints();
       
        this->CellScalars->SetNumberOfComponents(1);
        this->CellScalars->SetNumberOfTuples(c);
        
        i=0;
        while(i<c)
          {
          vtkIdType ptId=this->Polygon->GetPointId(i);
          this->CellScalars->SetValue(i,this->PointScalars->GetValue(ptId));
          ++i;
          }
        
        // I made a copy of the input point data, so it is OK to modify.
        this->InPD->SetScalars(this->PointScalars); // void
        
        if ( this->SortBy == VTK_SORT_BY_CELL )
          {
          double value = this->ContourValues->GetValue(this->Iter);
          this->Polygon->Contour(value, this->CellScalars, this->Locator, 
                                 this->NewVerts,this->NewLines,this->NewPolys,
                                 this->InPD,this->OutPD,this->InCD,cellId,
                                 this->OutCD);
          }
        else
          {
          // VTK_SORT_BY_VALUE
          int iter=0;
          while(iter<numContours)
            {
            double value = this->ContourValues->GetValue(iter);
            this->Polygon->Contour(value, this->CellScalars, this->Locator, 
                                   this->NewVerts,this->NewLines,
                                   this->NewPolys,this->InPD,this->OutPD,
                                   this->InCD,cellId,this->OutCD);
            
            ++iter;
            }
          }
        }
      }
  }
  
  //----------------------------------------------------------------------------
  // Description:
  // (i,j,k) are point coordinates at last level.
  // \pre all_set: this->Input->GetDimension()==2 implies ptIndices[2]==0
  double vtkHyperOctreeContourFilter::ComputePointValue(int ptIndices[3])
  {
    // this->Input->GetDimension()==2 implies ptIndices[2]==0
    // A implies B: !A||B
    assert("pre: all_set" &&
           (this->Input->GetDimension()!=2 || ptIndices[2]==0));
    
    int target[3];
    double result=0;
    double nb=0;
    
    int kmax;
    int jmax;
    int level=this->Input->GetNumberOfLevels();
    int maxCellIdx=1<<(level-1);
    
    switch(this->Input->GetDimension())
      {
      case 3:
        kmax=2;
        jmax=2;
        break;
      case 2:
        kmax=1;
        jmax=2;
        break;
      default:
        kmax=1;
        jmax=1;
        break;
      }
    
    int k=0;
    while(k<kmax)
      {
      target[2]=ptIndices[2]-k;
      int j=0;
      while(j<jmax)
        {
        target[1]=ptIndices[1]-j;
        int i=0;
        while(i<2)
          {
          target[0]=ptIndices[0]-i;
          if(target[0]>=0 && target[0]<maxCellIdx && target[1]>=0 && target[1]<maxCellIdx &&target[2]>=0 && target[2]<maxCellIdx )
            {
            this->NeighborCursor->MoveToNode(target,level-1);
            result+=this->InScalars->GetTuple1(this->NeighborCursor->GetLeafId());
            ++nb;
            }
          ++i;
          }
        ++j;
        }
      ++k;
      }
    
    if(nb>1)
      {
      result/=nb;
      }
    return result;
  }
  
  //----------------------------------------------------------------------------
  // Specify a spatial locator for merging points. By default, 
  // an instance of vtkMergePoints is used.
  void vtkHyperOctreeContourFilter::SetLocator(vtkPointLocator *locator)
  {
    if ( this->Locator == locator)
      {
      return;
      }
    
    if ( this->Locator )
      {
      this->Locator->UnRegister(this);
      this->Locator = NULL;
      }
  
    if ( locator )
      {
      locator->Register(this);
      }
  
    this->Locator = locator;
    this->Modified();
  }
  
  //----------------------------------------------------------------------------
  void vtkHyperOctreeContourFilter::CreateDefaultLocator()
  {
    if ( this->Locator == NULL )
      {
      this->Locator = vtkMergePoints::New();
      this->Locator->Register(this);
      this->Locator->Delete();
      }
  }
  
  //----------------------------------------------------------------------------
  int vtkHyperOctreeContourFilter::RequestUpdateExtent(
    vtkInformation *,
    vtkInformationVector **inputVector,
    vtkInformationVector *)
  {
    vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
    inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);
    return 1;
  }
  
  //----------------------------------------------------------------------------
  int vtkHyperOctreeContourFilter::FillInputPortInformation(int,
                                                            vtkInformation *info)
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkHyperOctree");
    return 1;
  }
  
  
  vtkCxxRevisionMacro(vtkHyperOctreeContourPointsGrabber, "$Revision: 1.4 $");
  vtkStandardNewMacro(vtkHyperOctreeContourPointsGrabber);
  
  //-----------------------------------------------------------------------------
  // Description:
  // Default constructor.
  vtkHyperOctreeContourPointsGrabber::vtkHyperOctreeContourPointsGrabber()
  {
    this->Triangulator=vtkOrderedTriangulator::New();
    this->Locator=vtkMergePoints::New();
    this->Polygon=0;
    this->Dimension=3;
    this->IdSet=new vtkHyperOctreeIdSet;
    this->Filter=0;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Destructor.
  vtkHyperOctreeContourPointsGrabber::~vtkHyperOctreeContourPointsGrabber()
  {
    if(this->Triangulator!=0)
      {
      this->Triangulator->UnRegister(this);
      delete this->IdSet;
      }
    if(this->Polygon!=0)
      {
      this->Polygon->UnRegister(this);
      }
    this->Locator->UnRegister(this);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Set the dimension of the hyperoctree.
  // \pre valid_dim: (dim==2 || dim==3)
  // \post is_set: GetDimension()==dim
  void vtkHyperOctreeContourPointsGrabber::SetDimension(int dim)
  {
    assert("pre: valid_dim" && (dim==2 || dim==3));
    if(dim!=this->Dimension)
      {
      if(dim==3)
        {
        this->Polygon->UnRegister(this);
        this->Polygon=0;
        this->Triangulator=vtkOrderedTriangulator::New();
        this->IdSet=new vtkHyperOctreeIdSet;
        }
      else
        {
         this->Triangulator->UnRegister(this);
         this->Triangulator=0;
         delete this->IdSet;
         this->Polygon=vtkPolygon::New();
        }
      this->Dimension=dim;
      }
    assert("post: is_set" && GetDimension()==dim);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Initialize the points insertion scheme.
  // Actually, it is just a trick to initialize the IdSet from the filter.
  // The IdSet class cannot be shared with the filter because it is a Pimpl.
  // It is used by clip,cut and contour filters to build the points
  // that lie on an hyperoctant.
  // \pre only_in_3d: GetDimension()==3
  void vtkHyperOctreeContourPointsGrabber::InitPointInsertion()
  {
    assert("pre: only_in_3d" && this->GetDimension()==3);
    this->IdSet->Set.clear();
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Insert a point, assuming the point is unique and does not require a
  // locator. Tt does not mean it does not use a locator. It just mean that
  // some implementation may skip the use of a locator.
  void vtkHyperOctreeContourPointsGrabber::InsertPoint(vtkIdType vtkNotUsed(ptId),
                                                       double pt[3],
                                                       double pcoords[3],
                                                       int ijk[3])
  {
    if(this->Locator->InsertUniquePoint(pcoords,this->LastPtId))
      {
      double value=this->Filter->ComputePointValue(ijk);
      this->Filter->PointScalars->InsertValue(this->LastPtId,value);
      }
    this->Triangulator->InsertPoint(this->LastPtId,pt,pcoords,0);
  }
  
  
  //-----------------------------------------------------------------------------
  // Description:
  // Insert a point using a locator.
  void vtkHyperOctreeContourPointsGrabber::InsertPointWithMerge(
    vtkIdType vtkNotUsed(ptId),
    double pt[3],
    double pcoords[3],
    int ijk[3])
  {
    if(this->Locator->InsertUniquePoint(pcoords,this->LastPtId))
      {
      double value=this->Filter->ComputePointValue(ijk);
      this->Filter->PointScalars->InsertValue(this->LastPtId,value);
      }
    if(this->IdSet->Set.find(this->LastPtId)==this->IdSet->Set.end()) // not find
      {
      this->IdSet->Set.insert(this->LastPtId);
      this->Triangulator->InsertPoint(this->LastPtId,pt,pcoords,0);
      }
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Insert a point in the quadtree case.
  void vtkHyperOctreeContourPointsGrabber::InsertPoint2D(double pt[3],
                                                         int ijk[3])
  {
    if(this->Locator->InsertUniquePoint(pt,this->LastPtId))
      {
      ijk[2]=0;
      double value=this->Filter->ComputePointValue(ijk);
      this->Filter->PointScalars->InsertValue(this->LastPtId,value);
      }
      
    this->Polygon->GetPointIds()->InsertNextId(this->LastPtId);
    this->Polygon->GetPoints()->InsertNextPoint(pt);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Return the ordered triangulator.
  vtkOrderedTriangulator *vtkHyperOctreeContourPointsGrabber::GetTriangulator()
  {
    return this->Triangulator;
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Return the polygon.
  vtkPolygon *vtkHyperOctreeContourPointsGrabber::GetPolygon()
  {
    return this->Polygon;
  }
  
  //-----------------------------------------------------------------------------
  void vtkHyperOctreeContourPointsGrabber::PrintSelf(ostream& os,
                                                     vtkIndent indent)
  {
    this->Superclass::PrintSelf(os,indent);
  }
  
  //-----------------------------------------------------------------------------
  // Description:
  // Return the id in the locator after a call to InsertPoint*().
  vtkIdType vtkHyperOctreeContourPointsGrabber::GetLastPtId()
  {
    return this->LastPtId;
  }
  //-----------------------------------------------------------------------------
  // Description:
  // Init the bounds of the locator.
  void vtkHyperOctreeContourPointsGrabber::InitLocator(vtkPoints *pts,
                                                       double bounds[6])
  {
    this->Locator->InitPointInsertion(pts,bounds);
  }
  //-----------------------------------------------------------------------------
  void vtkHyperOctreeContourPointsGrabber::SetFilter(
    vtkHyperOctreeContourFilter *filter)
  {
    this->Filter=filter;
  }
  