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

  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 "vtkHyperTreeGridCursor.h"
#include "vtkDataSetAttributes.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkRenderer.h"
#include "vtkCamera.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->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;
}

//-----------------------------------------------------------------------------
vtkAdaptiveDataSetSurfaceFilter::~vtkAdaptiveDataSetSurfaceFilter() = default;

//----------------------------------------------------------------------------
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 << "ParallelProjection: " << this->ParallelProjection << endl;
  os << indent << "Scale: " << this->Scale << 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->ParallelProjection && this->Renderer )
  {
    // Generate planes XY, XZ o YZ
    unsigned int* gridSize = input->GetGridSize();
    if ( gridSize[0] == 1 )
    {
      this->Axis1 = 1;
      this->Axis2 = 2;
    }
    else if ( gridSize[1] == 1 )
    {
      this->Axis1 = 0;
      this->Axis2 = 2;
    }
    else if ( gridSize[2] == 1 )
    {
      this->Axis1 = 0;
      this->Axis2 = 1;
    }

    // Compute the window size
    double windowSize = std::min(
      ( double ) this->LastRendererSize[0]/gridSize[this->Axis1],
      ( double ) this->LastRendererSize[1]/gridSize[this->Axis2] );

    // Compute the zoom of the camera
    vtkCamera* cam = this->Renderer->GetActiveCamera();
    double height = cam->GetParallelScale() * 2;
    double bounds[6];
    input->GetBounds( bounds );
    double zoom1 = ( bounds[2 * this->Axis1 + 1] - bounds[2 * this->Axis1] ) / height;
    double zoom2 = ( bounds[2 * this->Axis2 + 1] - bounds[2 * this->Axis2] ) / height;
    double zoom = std::max( zoom1, zoom2 );

    // Compute how many levels of the tree we should process
    int f = input->GetBranchFactor();
    this->LevelMax = ( log( windowSize ) + log( zoom / this->Scale ) ) / log( f ) + 1;
    if ( this->LevelMax < 0 )
    {
      this->LevelMax = 0;
    }
    double ratio = ( ( double ) this->LastRendererSize[0] ) / this->LastRendererSize[1];
    this->Radius = cam->GetParallelScale() * sqrt( 1 + ratio * ratio );
  }
  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 )
{
  // Create storage for corners of leaf cells
  this->Points = vtkPoints::New();

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

  // Retrieve material mask
  vtkBitArray* mask = input->HasMaterialMask() ? input->GetMaterialMask() : nullptr;

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

  // Iterate over all hyper trees
  vtkIdType index;
  vtkHyperTreeGrid::vtkHyperTreeGridIterator it;
  input->InitializeTreeIterator( it );
  while ( it.GetNextTree( index ) )
  {
    // Initialize new cursor at root of current tree
    vtkHyperTreeGridCursor* cursor;
    if ( this->Dimension == 3 )
    {
      // In 3 dimensions, von Neumann neighborhood information is needed
      cursor = input->NewVonNeumannSuperCursor( index );
    } // if ( this->Dimension == 3 )
    else
    {
      // Otherwise, geometric properties of the cells suffice
      cursor = input->NewGeometricCursor( index );
    } // else

    // If this is not a ghost tree
    if ( ! ghost || ! ghost->GetTuple1( cursor->GetGlobalNodeIndex() ) )
    {
      // Build geometry recursively
      this->RecursivelyProcessTree( cursor, mask, 0 );
    }

    // Clean up
    cursor->Delete();
  } // it

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

  this->Points->Delete();
  this->Cells->Delete();
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::RecursivelyProcessTree( vtkHyperTreeGridCursor* cursor,
                                                              vtkBitArray* mask,
                                                              int level )
{
  // Retrieve input grid
  vtkHyperTreeGrid* input = cursor->GetGrid();

  if ( this->Dimension == 3 )
  {
    // Create geometry output if cursor is at leaf
    if ( cursor->IsLeaf() )
    {
      this->ProcessLeaf3D( cursor, mask );
    } // if ( cursor->IsLeaf() )
    else
    {
      // Cursor is not at leaf, recurse to all children
      int numChildren = input->GetNumberOfChildren();
      for ( int child = 0; child < numChildren; ++ child )
      {
        // Create child cursor from parent
        vtkHyperTreeGridCursor* childCursor = cursor->Clone();
        childCursor->ToChild( child );

        // Recurse
        this->RecursivelyProcessTree( childCursor, mask, level+1 );

        // Clean up
        childCursor->Delete();
        childCursor = nullptr;
      } // child
    } // else
  } // if ( this->Dimension == 3 )
  else
  {
    bool insideBB = ( this->LevelMax == -1 );
    if( ! insideBB )
    {
      // Check if the current node of the tree is going to be rendered
      double half = std::max( cursor->GetSize()[this->Axis1] / 2,
                              cursor->GetSize()[this->Axis2] / 2 );
      insideBB = (
        pow( cursor->GetOrigin()[this->Axis1] + half - this->LastCameraFocalPoint[this->Axis1], 2 ) +
        pow( cursor->GetOrigin()[this->Axis2] + half - this->LastCameraFocalPoint[this->Axis2], 2 ) ) <
        pow( this->Radius + half * sqrt(2.), 2 );
    }
    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, mask );
        }
        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 = input->GetNumberOfChildren();
        for ( int child = 0; child < numChildren; ++ child )
        {
          // Create child cursor from parent
          vtkHyperTreeGridCursor* childCursor = cursor->Clone();
          childCursor->ToChild( child );

          // Recurse
          this->RecursivelyProcessTree( childCursor, mask, level+1 );

          // Clean up
          childCursor->Delete();
          childCursor = nullptr;
        } // child
      } // else
    } // if( insideBB )
  } // else
}

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

  // First endpoint is at origin of cursor
  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( vtkHyperTreeGridCursor* cursor,
                                                     vtkBitArray* mask )

{
  // 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 ( ! mask  || ! mask->GetValue( id ) )
  {
    // Insert face into 2D geometry depending on orientation
    this->AddFace( id, cursor->GetOrigin(), cursor->GetSize(), 0, this->Orientation );
  }
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::ProcessLeaf3D( vtkHyperTreeGridCursor* superCursor,
                                                     vtkBitArray* mask )
{
  // Cell at super cursor center is a leaf, retrieve its global index, level, and mask
  vtkIdType id = superCursor->GetGlobalNodeIndex();
  unsigned level = superCursor->GetLevel();
  int masked = mask ? mask->GetValue( id ) : 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
    vtkHyperTreeGridCursor* cursorN = superCursor->GetCursor( VonNeumannCursors3D[c] );

    // Retrieve tree, leaf flag, and mask of neighbor cursor
    vtkHyperTree* treeN = cursorN->GetTree();
    bool leafN = cursorN->IsLeaf();
    vtkIdType idN = cursorN->GetGlobalNodeIndex();
    int maskedN  = mask ? mask->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 && cursorN->GetLevel() < level && ! maskedN ) )
    {
      // Generate face with corresponding normal and offset
      this->AddFace( id, superCursor->GetOrigin(), superCursor->GetSize(),
                     VonNeumannOffsets3D[c], VonNeumannOrientations3D[c] );
    }
  } // c
}

//----------------------------------------------------------------------------
void vtkAdaptiveDataSetSurfaceFilter::AddFace( vtkIdType inId,
                                               double* origin,
                                               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 ( 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();
}
