#include "vtkCMBMeshTerrainWithArcs.h"

#include "vtkAppendPolyData.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCellLocator.h"
#include "vtkCleanPolyData.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkGenericCell.h"
#include "vtkIncrementalOctreePointLocator.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"
#include "vtkPolygon.h"

#include "smtk/extension/vtk/meshing/vtkCMBPrepareForTriangleMesher.h"
#include "smtk/extension/vtk/meshing/vtkCMBTriangleMesher.h"

#include "vtkMultiBlockDataSet.h"
#include "vtkXMLPolyDataWriter.h"

#include <iostream>

vtkStandardNewMacro(vtkCMBMeshTerrainWithArcs);

namespace
{
template <class T>
void zeroZValues(T* t, int size)
{
  for (int i = 0; i < size; i += 3, t += 3)
  {
    t[2] = 0.0;
  }
}

struct PInfo
{
  int NumberOfCells;
  int Id;
  double Height;
  PInfo()
    : NumberOfCells(0)
    , Id(0)
    , Height(0.0)
  {
  }
};
};

class vtkCMBMeshTerrainWithArcs::vtkCmbInternalTerrainInfo
{
public:
  vtkCmbInternalTerrainInfo(
    vtkInformationVector* input, const double& radius, const double& groundElevation);
  ~vtkCmbInternalTerrainInfo()
  {
    if (this->Locator)
    {
      this->Locator->Delete();
    }
  }

  template <class T>
  T getElevation(T* point);
  void setRadius(const double& r) { Radius = r; }
  void setGroundElevation(const double& ge) { GroundElevation = ge; }
protected:
  double Radius;
  double GroundElevation;
  std::map<vtkIdType, double> IdToElevation;
  vtkIncrementalOctreePointLocator* Locator;
};

class vtkCMBMeshTerrainWithArcs::vtkCmbPolygonInfo
{
public:
  std::vector<PInfo> Info;
  int GroundIndex;
};

vtkCMBMeshTerrainWithArcs::vtkCmbInternalTerrainInfo::vtkCmbInternalTerrainInfo(
  vtkInformationVector* input, const double& radius, const double& groundElevation)
  : Radius(radius)
  , GroundElevation(groundElevation)
{

  vtkIdType numPoints = 0;
  int numInputs = input->GetNumberOfInformationObjects();
  vtkPolyData* pd = NULL;
  for (int idx = 0; idx < numInputs; ++idx)
  {
    pd = vtkPolyData::GetData(input, idx);
    if (pd)
    {
      numPoints += pd->GetNumberOfPoints();
    }
  }

  vtkPoints* inputPoints = NULL;
  vtkPoints* points = vtkPoints::New();
  points->SetNumberOfPoints(numPoints);

  vtkIdType index = 0;
  double p[3];
  for (int idx = 0; idx < numInputs; ++idx)
  {
    pd = vtkPolyData::GetData(input, idx);
    if (pd)
    {
      inputPoints = pd->GetPoints();
      vtkIdType size = inputPoints->GetNumberOfPoints();
      for (vtkIdType i = 0; i < size; ++i)
      {

        inputPoints->GetPoint(i, p);

        double elev = p[2];
        p[2] = 0.0;
        points->InsertPoint(index, p);
        this->IdToElevation.insert(std::pair<vtkIdType, double>(index, elev));
        ++index;
      }
    }
  }

  vtkPolyData* pointSet = vtkPolyData::New();
  pointSet->SetPoints(points);
  points->FastDelete();

  this->Locator = vtkIncrementalOctreePointLocator::New();
  this->Locator->AutomaticOn();
  this->Locator->SetTolerance(0.0);
  this->Locator->SetDataSet(pointSet);
  this->Locator->BuildLocator();
  pointSet->Delete();
}

template <class T>
T vtkCMBMeshTerrainWithArcs::vtkCmbInternalTerrainInfo::getElevation(T* point)
{
  double dpoint[3];
  dpoint[0] = static_cast<double>(point[0]);
  dpoint[1] = static_cast<double>(point[1]);
  dpoint[2] = static_cast<double>(point[2]);
  vtkIdList* ids = vtkIdList::New();
  this->Locator->FindPointsWithinRadius(this->Radius, dpoint, ids);
  double sum = 0;
  std::map<vtkIdType, double>::const_iterator it;
  vtkIdType size = ids->GetNumberOfIds();
  for (vtkIdType i = 0; i < size; ++i)
  {

    it = this->IdToElevation.find(ids->GetId(i));
    if (it != this->IdToElevation.end())
    {
      sum += it->second;
    }
  }
  ids->Delete();

  T elev = static_cast<T>((size == 0) ? this->GroundElevation : sum / size);
  return elev;
}

vtkCMBMeshTerrainWithArcs::vtkCMBMeshTerrainWithArcs()
{
  this->PolygonInfo = new vtkCmbPolygonInfo();
  this->TerrainInfo = NULL;
  this->SetNumberOfInputPorts(2);
  this->ElevationRadius = 1.0;
  this->VOIBounds[0] = 0.0;
  this->VOIBounds[1] = 1.0;
  this->VOIBounds[2] = 0.0;
  this->VOIBounds[3] = 1.0;
  this->VOIBounds[4] = 0.0;
  this->VOIBounds[5] = 1.0;

  this->NumberOfProgressSteps = 10;
  this->StepIncrement = 0.1;

  this->Mesher = vtkCMBTriangleMesher::New();
  this->MesherMaxArea = 0.125;
}

vtkCMBMeshTerrainWithArcs::~vtkCMBMeshTerrainWithArcs()
{
  delete this->PolygonInfo;
  delete this->TerrainInfo;
  this->Mesher->Delete();
}

int vtkCMBMeshTerrainWithArcs::FillInputPortInformation(int port, vtkInformation* info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
  info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
  if (port == 1)
  {
    info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
  }
  return 1;
}

void vtkCMBMeshTerrainWithArcs::RemoveInputConnections()
{
  this->SetInputConnection(0, NULL);
}

void vtkCMBMeshTerrainWithArcs::RemoveSourceConnections()
{
  this->SetInputConnection(1, NULL);
}

void vtkCMBMeshTerrainWithArcs::DetermineNumberOfProgressSteps(const int& numInputs)
{

  this->NumberOfProgressSteps = 4 + numInputs * 4;
  this->CurrentProgressStep = 0;
  this->StepIncrement = 0.9 / this->NumberOfProgressSteps;
  this->SetProgressText("vtkCMBMeshTerrainWithArcs");

}

void vtkCMBMeshTerrainWithArcs::NextProgressStep()
{
  if (this->CurrentProgressStep < this->NumberOfProgressSteps)
  {
    this->UpdateProgress(this->StepIncrement * this->CurrentProgressStep);
    ++this->CurrentProgressStep;
  }
  else
  {
    vtkWarningMacro("Updating progress in vtkCMBMeshTerrainWithArcs more times than expected");
  }
}

int vtkCMBMeshTerrainWithArcs::RequestData(
  vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  bool validMesh = true;

  this->DetermineNumberOfProgressSteps(inputVector[1]->GetNumberOfInformationObjects());

  vtkNew<vtkPolyData> groundMesh;
  vtkPolyData* mesherResult = NULL;

  bool readyForGroundMeshing = this->PrepForMeshing(inputVector[1], groundMesh.GetPointer());

  if (readyForGroundMeshing)
  {

    this->Mesher->SetPreserveBoundaries(false);
    this->Mesher->SetPreserveEdgesAndNodes(true);
    this->Mesher->SetMaxAreaMode(vtkCMBTriangleMesher::RelativeToBoundsAndSegments);
    this->Mesher->SetMaxArea(0.015);
    this->Mesher->SetUseMinAngle(true);
    this->Mesher->SetInputData(groundMesh.GetPointer());
    this->Mesher->Update();
    this->NextProgressStep();

    this->MesherMaxArea = this->Mesher->GetComputedMaxArea();

    mesherResult = vtkPolyData::New();
    mesherResult->ShallowCopy(this->Mesher->GetOutput());
  }
  else
  {
    return false;
  }

  if (validMesh && mesherResult)
  {

    this->TerrainInfo =
      new vtkCmbInternalTerrainInfo(inputVector[0], this->ElevationRadius, this->VOIBounds[4]);

    vtkPolyData* finalMesh = vtkPolyData::New();
    finalMesh->DeepCopy(mesherResult);
    int numOfArrays = finalMesh->GetCellData()->GetNumberOfArrays();
    for (int j = numOfArrays - 1; j >= 0; --j)
    {
      finalMesh->GetCellData()->RemoveArray(j);
    }
    finalMesh->SetLines(NULL); 

    validMesh = this->GenerateGroundMesh(finalMesh);
    if (validMesh && this->PolygonInfo->Info.size() > 1)
    {

      this->GenerateExtrudedArcSets(mesherResult, finalMesh);
    }
    mesherResult->Delete();

    delete this->TerrainInfo;
    this->TerrainInfo = NULL;

    vtkPolyData* output = vtkPolyData::GetData(outputVector);
    output->ShallowCopy(finalMesh);
    finalMesh->Delete();
  }

  this->UpdateProgress(1.0);
  return validMesh;
}

bool vtkCMBMeshTerrainWithArcs::PrepForMeshing(vtkInformationVector* input, vtkPolyData* mesh)
{

  mesh->SetVerts(NULL);
  mesh->SetPolys(NULL);
  mesh->SetStrips(NULL);

  this->PolygonInfo->Info.clear();

  vtkCellArray* newCells = vtkCellArray::New();

  vtkCellArray* cells = NULL;
  int numInputs = input->GetNumberOfInformationObjects();
  for (int idx = 0; idx < numInputs; ++idx)
  {
    cells = NULL;
    vtkPolyData* pd = vtkPolyData::GetData(input, idx);
    cells = pd->GetLines();
    if (!cells)
    {
      continue;
    }

    if (cells->GetNumberOfCells() > 0)
    {
      PInfo info;
      int numCells = 0;
      vtkIdType npts, *pts;
      bool firstCell = true;
      for (cells->InitTraversal(); cells->GetNextCell(npts, pts);)
      {
        if (firstCell)
        {

          info.Height = pd->GetPoint(pts[0])[2];
          firstCell = false;
        }
        for (vtkIdType i = 0; i < (npts - 1); ++i, ++numCells)
        {
          newCells->InsertNextCell(2, pts + i);
        }
      }

      info.Id = static_cast<int>(this->PolygonInfo->Info.size()) + 1;
      info.NumberOfCells = numCells;
      this->PolygonInfo->Info.push_back(info);
      this->NextProgressStep();
    }
  }

  vtkNew<vtkPoints> points;
  const int numInputsToAppend = input->GetNumberOfInformationObjects();
  if (numInputsToAppend > 0)
  {
    vtkNew<vtkAppendPolyData> append;
    for (int i = 0; i < numInputsToAppend; ++i)
    {
      append->AddInputData(vtkPolyData::GetData(input, i));
    }
    append->Update();

    points->ShallowCopy(append->GetOutput()->GetPoints());
  }

  this->InsertGroundPlane(points.GetPointer(), newCells);

  vtkIdType size = newCells->GetNumberOfCells();
  mesh->SetLines(newCells);
  newCells->FastDelete();

  bool validMesh = true;

  validMesh &= this->AssignPolygonIds(mesh, size);
  this->NextProgressStep();

  validMesh &= this->FlattenMesh(points.GetPointer());
  this->NextProgressStep();

  mesh->SetPoints(points.GetPointer());

  return validMesh;
}

bool vtkCMBMeshTerrainWithArcs::AssignPolygonIds(vtkPolyData* mesh, const vtkIdType& ) const
{
  typedef vtkCMBPrepareForTriangleMesher vtkPrepareForMesher;
  vtkNew<vtkPrepareForMesher> mapInterface;
  mapInterface->SetPolyData(mesh); 
  mapInterface->SetNumberOfArcs(this->PolygonInfo->Info.size());
  mapInterface->SetNumberOfLoops(this->PolygonInfo->Info.size());
  mapInterface->InitializeNewMapInfo();

  size_t i = 0;
  vtkIdType index = 0;
  int numCells = 0;

  vtkIdType groundId = this->PolygonInfo->Info[this->PolygonInfo->GroundIndex].Id;
  int asId = 0;
  int endpointId = 0;
  for (i = 0; i < this->PolygonInfo->Info.size(); ++i)
  {
    numCells = this->PolygonInfo->Info[i].NumberOfCells;
    asId = this->PolygonInfo->Info[i].Id;
    if (asId == groundId)
    {

      vtkIdType loopId = mapInterface->AddLoop(groundId, -1);

      mapInterface->AddArc(index * 3, numCells * 3, groundId, loopId, -1, endpointId, endpointId);
      endpointId++;

      index += numCells;
    }
    else
    {

      vtkIdType loopId = mapInterface->AddLoop(-1, groundId);
      mapInterface->AddArc(index * 3, numCells * 3, asId, loopId, -1, endpointId, endpointId);
      endpointId++;

      index += numCells;
    }
  }
  mapInterface->FinalizeNewMapInfo(); 

  return true;
}

bool vtkCMBMeshTerrainWithArcs::FlattenMesh(vtkPoints* points) const
{
  bool valid = true;
  vtkDataArray* dataArray = points->GetData();
  int size = points->GetNumberOfPoints() * 3;
  if (dataArray->GetDataType() == VTK_FLOAT)
  {
    vtkFloatArray* floatArray = static_cast<vtkFloatArray*>(dataArray);
    float* pt = floatArray->GetPointer(0);
    zeroZValues(pt, size);
  }
  else if (dataArray->GetDataType() == VTK_DOUBLE)
  {
    vtkDoubleArray* doubleArray = static_cast<vtkDoubleArray*>(dataArray);
    double* pt = doubleArray->GetPointer(0);
    zeroZValues(pt, size);
  }
  else
  {
    valid = false;
  }
  return valid;
}

void vtkCMBMeshTerrainWithArcs::InsertGroundPlane(vtkPoints* points, vtkCellArray* newCells) const
{
  vtkIdType bl, br, tr, tl;
  bl = points->InsertNextPoint(this->VOIBounds[0], this->VOIBounds[2], 0); 
  br = points->InsertNextPoint(this->VOIBounds[1], this->VOIBounds[2], 0); 
  tr = points->InsertNextPoint(this->VOIBounds[1], this->VOIBounds[3], 0); 
  tl = points->InsertNextPoint(this->VOIBounds[0], this->VOIBounds[3], 0); 
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(br);
  newCells->InsertCellPoint(bl);
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(bl);
  newCells->InsertCellPoint(tl);
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(tl);
  newCells->InsertCellPoint(tr);
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(tr);
  newCells->InsertCellPoint(br);

  PInfo info;

  info.Id = static_cast<int>(this->PolygonInfo->Info.size()) + 1;
  info.NumberOfCells = 4;
  this->PolygonInfo->Info.push_back(info);
  this->PolygonInfo->GroundIndex = static_cast<int>(this->PolygonInfo->Info.size() - 1);
}

bool vtkCMBMeshTerrainWithArcs::GenerateGroundMesh(vtkPolyData* finalMesh)
{
  bool valid = this->ExtrudeMeshPoints(finalMesh->GetPoints());
  this->NextProgressStep();
  return valid;
}

bool vtkCMBMeshTerrainWithArcs::GenerateExtrudedArcSets(vtkPolyData* input, vtkPolyData* output)
{
  bool valid;
  for (size_t i = 0; i < this->PolygonInfo->Info.size(); i++)
  {

    if (static_cast<size_t>(this->PolygonInfo->GroundIndex) == i)
    {
      continue;
    }
    vtkNew<vtkPolyData> arcsetMesh;
    int id = this->PolygonInfo->Info[i].Id;
    this->CreateArcSetForMeshing(input, arcsetMesh.GetPointer(), id);
    this->NextProgressStep();

    this->Mesher->SetPreserveBoundaries(true);
    this->Mesher->SetMaxAreaMode(vtkCMBTriangleMesher::AbsoluteArea);
    this->Mesher->SetMaxArea(this->MesherMaxArea);
    this->Mesher->SetInputData(arcsetMesh.GetPointer());
    this->Mesher->Update();
    this->NextProgressStep();

    vtkNew<vtkPolyData> resultMesh;
    resultMesh->ShallowCopy(this->Mesher->GetOutput());

    int numOfArrays = resultMesh->GetCellData()->GetNumberOfArrays();
    for (int j = numOfArrays - 1; j >= 0; --j)
    {
      resultMesh->GetCellData()->RemoveArray(j);
    }

    valid = this->ExtrudeArcMesh(resultMesh.GetPointer(), static_cast<int>(i));
    this->NextProgressStep();
    if (valid)
    {

      vtkNew<vtkAppendPolyData> append;
      append->AddInputData(output);
      append->AddInputData(resultMesh.GetPointer());
      append->Update();
      vtkPolyData* result = append->GetOutput();
      output->ShallowCopy(result);
    }
  }

  return true;
}

void vtkCMBMeshTerrainWithArcs::CreateArcSetForMeshing(
  vtkPolyData* input, vtkPolyData* output, const int& index) const
{
  vtkCellData* cellData = input->GetCellData();
  vtkPoints* points = input->GetPoints();
  if (cellData == NULL || points == NULL)
  {
    return;
  }

  vtkPoints* newPoints = vtkPoints::New();
  newPoints->DeepCopy(points);
  output->SetPoints(newPoints);
  newPoints->Delete();

  vtkCellArray* newCells = vtkCellArray::New();
  newCells->Allocate(input->GetNumberOfLines());

  vtkIdTypeArray* elementIds = vtkIdTypeArray::SafeDownCast(cellData->GetArray("ElementIds"));
  vtkIdType numCells = input->GetNumberOfCells();

  for (vtkIdType i = 0; i < numCells; ++i)
  {
    if ((input->GetCellType(i) == VTK_LINE || input->GetCellType(i) == VTK_POLY_LINE) &&
      elementIds->GetValue(i) == index)
    {
      newCells->InsertNextCell(input->GetCell(i));
    }
  }

  output->SetLines(newCells);
  vtkIdType size = output->GetNumberOfCells();
  newCells->FastDelete();

  typedef vtkCMBPrepareForTriangleMesher vtkPrepareForMesher;
  vtkNew<vtkPrepareForMesher> mapInterface;
  mapInterface->SetPolyData(output); 
  mapInterface->SetNumberOfArcs(1);
  mapInterface->SetNumberOfLoops(1);
  mapInterface->InitializeNewMapInfo();

  vtkIdType loopId = mapInterface->AddLoop(index, -1);
  mapInterface->AddArc(0, size * 3, 1, loopId, -1, 1, 1);

  mapInterface->FinalizeNewMapInfo();
}

bool vtkCMBMeshTerrainWithArcs::ExtrudeArcMesh(vtkPolyData* mesh, const int& index) const
{

  if (mesh == NULL)
  {
    return false;
  }
  vtkPoints* points = mesh->GetPoints();
  if (points == NULL)
  {
    return false;
  }

  double height = this->PolygonInfo->Info[index].Height;
  this->UnFlattenMesh(points, height);

  mesh->BuildCells();

  std::map<vtkIdType, vtkIdType> SidesBottomToTop;
  vtkIdType meshId, newGroundId, npts;
  vtkIdList* pIds = vtkIdList::New();
  double p[3];

  vtkCellArray* lines = mesh->GetLines();
  lines->InitTraversal();
  while (lines->GetNextCell(pIds))
  {
    npts = pIds->GetNumberOfIds();
    for (vtkIdType i = npts - 1; i >= 0; --i)
    {
      meshId = pIds->GetId(i);
      if (SidesBottomToTop.find(meshId) == SidesBottomToTop.end())
      {
        points->GetPoint(meshId, p);
        p[2] = 0; 
        p[2] = this->TerrainInfo->getElevation(p);
        newGroundId = points->InsertNextPoint(p);
      }
      else
      {
        newGroundId = SidesBottomToTop.find(meshId)->second;
      }
      SidesBottomToTop.insert(std::pair<vtkIdType, vtkIdType>(meshId, newGroundId));
      pIds->InsertNextId(newGroundId);
    }
    mesh->InsertNextCell(VTK_POLYGON, pIds);
    pIds->Reset();
  }

  pIds->Delete();

  mesh->SetLines(NULL);

  return true;
}

bool vtkCMBMeshTerrainWithArcs::UnFlattenMesh(vtkPoints* points, const double& height) const
{

  vtkDataArray* dataArray = points->GetData();
  vtkIdType size = points->GetNumberOfPoints();
  if (dataArray->GetDataType() == VTK_FLOAT)
  {
    vtkFloatArray* floatArray = static_cast<vtkFloatArray*>(dataArray);
    float* pos = floatArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = height;
    }
  }
  else
  {
    vtkDoubleArray* doubleArray = static_cast<vtkDoubleArray*>(dataArray);
    double* pos = doubleArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = height;
    }
  }
  return true;
}

bool vtkCMBMeshTerrainWithArcs::ExtrudeMeshPoints(vtkPoints* points) const
{

  vtkDataArray* dataArray = points->GetData();
  vtkIdType size = points->GetNumberOfPoints();
  if (dataArray->GetDataType() == VTK_FLOAT)
  {
    vtkFloatArray* floatArray = static_cast<vtkFloatArray*>(dataArray);
    float* pos = floatArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = this->TerrainInfo->getElevation(pos);
    }
  }
  else
  {
    vtkDoubleArray* doubleArray = static_cast<vtkDoubleArray*>(dataArray);
    double* pos = doubleArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = this->TerrainInfo->getElevation(pos);
    }
  }
  return true;
}

void vtkCMBMeshTerrainWithArcs::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "NumberOfProgressSteps: " << this->NumberOfProgressSteps << std::endl;
  os << indent << "CurrentProgressStep: " << this->CurrentProgressStep << std::endl;
  os << indent << "StepIncrement: " << this->StepIncrement << std::endl;
  os << indent << "VOIBounds: " << this->VOIBounds[0] << ", " << this->VOIBounds[1] << ", "
     << this->VOIBounds[2] << ", " << this->VOIBounds[3] << ", " << this->VOIBounds[4] << ", "
     << this->VOIBounds[5] << ", " << std::endl;
  os << indent << "ElevationRadius: " << this->ElevationRadius << std::endl;
  os << indent << "MesherMaxArea: " << this->MesherMaxArea << std::endl;
  os << indent << "Mesher: " << this->Mesher << std::endl;
}
//=========================================================================
//  Copyright (c) Kitware, Inc.
//  All rights reserved.
//  See LICENSE.txt 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 "vtkCMBMeshTerrainWithArcs.h"

#include "vtkAppendPolyData.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCellLocator.h"
#include "vtkCleanPolyData.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkGenericCell.h"
#include "vtkIncrementalOctreePointLocator.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"
#include "vtkPolygon.h"

#include "smtk/extension/vtk/meshing/vtkCMBPrepareForTriangleMesher.h"
#include "smtk/extension/vtk/meshing/vtkCMBTriangleMesher.h"

#include "vtkMultiBlockDataSet.h"
#include "vtkXMLPolyDataWriter.h"

#include <iostream>

vtkStandardNewMacro(vtkCMBMeshTerrainWithArcs);

namespace
{
template <class T>
void zeroZValues(T* t, int size)
{
  for (int i = 0; i < size; i += 3, t += 3)
  {
    t[2] = 0.0;
  }
}

struct PInfo
{
  int NumberOfCells;
  int Id;
  double Height;
  PInfo()
    : NumberOfCells(0)
    , Id(0)
    , Height(0.0)
  {
  }
};
};

class vtkCMBMeshTerrainWithArcs::vtkCmbInternalTerrainInfo
{
public:
  vtkCmbInternalTerrainInfo(
    vtkInformationVector* input, const double& radius, const double& groundElevation);
  ~vtkCmbInternalTerrainInfo()
  {
    if (this->Locator)
    {
      this->Locator->Delete();
    }
  }
  //returns the average elevation for a circle given the radius
  template <class T>
  T getElevation(T* point);
  void setRadius(const double& r) { Radius = r; }
  void setGroundElevation(const double& ge) { GroundElevation = ge; }
protected:
  double Radius;
  double GroundElevation;
  std::map<vtkIdType, double> IdToElevation;
  vtkIncrementalOctreePointLocator* Locator;
};

class vtkCMBMeshTerrainWithArcs::vtkCmbPolygonInfo
{
public:
  std::vector<PInfo> Info;
  int GroundIndex;
};

vtkCMBMeshTerrainWithArcs::vtkCmbInternalTerrainInfo::vtkCmbInternalTerrainInfo(
  vtkInformationVector* input, const double& radius, const double& groundElevation)
  : Radius(radius)
  , GroundElevation(groundElevation)
{
  //1. Go through all the input data objects and collect the size of all the point dataset
  //2. Create a super set of points while removing all the z values from the points
  //and storing them in the map.
  //3. Create locator of the resulting 2D point set

  vtkIdType numPoints = 0;
  int numInputs = input->GetNumberOfInformationObjects();
  vtkPolyData* pd = NULL;
  for (int idx = 0; idx < numInputs; ++idx)
  {
    pd = vtkPolyData::GetData(input, idx);
    if (pd)
    {
      numPoints += pd->GetNumberOfPoints();
    }
  }

  //second iteration is building the point set and elevation mapping
  vtkPoints* inputPoints = NULL;
  vtkPoints* points = vtkPoints::New();
  points->SetNumberOfPoints(numPoints);

  vtkIdType index = 0;
  double p[3];
  for (int idx = 0; idx < numInputs; ++idx)
  {
    pd = vtkPolyData::GetData(input, idx);
    if (pd)
    {
      inputPoints = pd->GetPoints();
      vtkIdType size = inputPoints->GetNumberOfPoints();
      for (vtkIdType i = 0; i < size; ++i)
      {
        //get the point
        inputPoints->GetPoint(i, p);

        //store the z value & flatten
        double elev = p[2];
        p[2] = 0.0;
        points->InsertPoint(index, p);
        this->IdToElevation.insert(std::pair<vtkIdType, double>(index, elev));
        ++index;
      }
    }
  }

  vtkPolyData* pointSet = vtkPolyData::New();
  pointSet->SetPoints(points);
  points->FastDelete();

  this->Locator = vtkIncrementalOctreePointLocator::New();
  this->Locator->AutomaticOn();
  this->Locator->SetTolerance(0.0);
  this->Locator->SetDataSet(pointSet);
  this->Locator->BuildLocator();
  pointSet->Delete();
}

template <class T>
T vtkCMBMeshTerrainWithArcs::vtkCmbInternalTerrainInfo::getElevation(T* point)
{
  double dpoint[3];
  dpoint[0] = static_cast<double>(point[0]);
  dpoint[1] = static_cast<double>(point[1]);
  dpoint[2] = static_cast<double>(point[2]);
  vtkIdList* ids = vtkIdList::New();
  this->Locator->FindPointsWithinRadius(this->Radius, dpoint, ids);
  double sum = 0;
  std::map<vtkIdType, double>::const_iterator it;
  vtkIdType size = ids->GetNumberOfIds();
  for (vtkIdType i = 0; i < size; ++i)
  {
    //average the elevation
    it = this->IdToElevation.find(ids->GetId(i));
    if (it != this->IdToElevation.end())
    {
      sum += it->second;
    }
  }
  ids->Delete();

  //handle the zero size use case
  T elev = static_cast<T>((size == 0) ? this->GroundElevation : sum / size);
  return elev;
}

vtkCMBMeshTerrainWithArcs::vtkCMBMeshTerrainWithArcs()
{
  this->PolygonInfo = new vtkCmbPolygonInfo();
  this->TerrainInfo = NULL;
  this->SetNumberOfInputPorts(2);
  this->ElevationRadius = 1.0;
  this->VOIBounds[0] = 0.0;
  this->VOIBounds[1] = 1.0;
  this->VOIBounds[2] = 0.0;
  this->VOIBounds[3] = 1.0;
  this->VOIBounds[4] = 0.0;
  this->VOIBounds[5] = 1.0;

  this->NumberOfProgressSteps = 10;
  this->StepIncrement = 0.1;

  this->Mesher = vtkCMBTriangleMesher::New();
  this->MesherMaxArea = 0.125;
}

vtkCMBMeshTerrainWithArcs::~vtkCMBMeshTerrainWithArcs()
{
  delete this->PolygonInfo;
  delete this->TerrainInfo;
  this->Mesher->Delete();
}

int vtkCMBMeshTerrainWithArcs::FillInputPortInformation(int port, vtkInformation* info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
  info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
  if (port == 1)
  {
    info->Set(vtkAlgorithm::INPUT_IS_OPTIONAL(), 1);
  }
  return 1;
}

void vtkCMBMeshTerrainWithArcs::RemoveInputConnections()
{
  this->SetInputConnection(0, NULL);
}

void vtkCMBMeshTerrainWithArcs::RemoveSourceConnections()
{
  this->SetInputConnection(1, NULL);
}

void vtkCMBMeshTerrainWithArcs::DetermineNumberOfProgressSteps(const int& numInputs)
{
  //the number of progress steps is:
  // prepForMeshing calls numInputs + 2 times
  // requestData 1
  // GenerateGroundMesh 1
  // GenerateExtrudedArcSets numInputs * 3

  this->NumberOfProgressSteps = 4 + numInputs * 4;
  this->CurrentProgressStep = 0;
  this->StepIncrement = 0.9 / this->NumberOfProgressSteps;
  this->SetProgressText("vtkCMBMeshTerrainWithArcs");

  //the end of request data is where we set progress to 1.0
}

void vtkCMBMeshTerrainWithArcs::NextProgressStep()
{
  if (this->CurrentProgressStep < this->NumberOfProgressSteps)
  {
    this->UpdateProgress(this->StepIncrement * this->CurrentProgressStep);
    ++this->CurrentProgressStep;
  }
  else
  {
    vtkWarningMacro("Updating progress in vtkCMBMeshTerrainWithArcs more times than expected");
  }
}

int vtkCMBMeshTerrainWithArcs::RequestData(
  vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  bool validMesh = true;

  this->DetermineNumberOfProgressSteps(inputVector[1]->GetNumberOfInformationObjects());

  //setup the input mesh for the mesher
  vtkNew<vtkPolyData> groundMesh;
  vtkPolyData* mesherResult = NULL;

  //we need to add 3 properties to the polydata to
  //correctly prep it for meshing
  bool readyForGroundMeshing = this->PrepForMeshing(inputVector[1], groundMesh.GetPointer());

  //now lets mesh it
  if (readyForGroundMeshing)
  {
    //first call to triangle will mesh only the ground plane with
    //all the arc sets as holes
    this->Mesher->SetPreserveBoundaries(false);
    this->Mesher->SetPreserveEdgesAndNodes(true);
    this->Mesher->SetMaxAreaMode(vtkCMBTriangleMesher::RelativeToBoundsAndSegments);
    this->Mesher->SetMaxArea(0.015);
    this->Mesher->SetUseMinAngle(true);
    this->Mesher->SetInputData(groundMesh.GetPointer());
    this->Mesher->Update();
    this->NextProgressStep();

    //Do not get the Max area, as that is just a parameter to computing
    //the actual max area depending on the max area mode.
    this->MesherMaxArea = this->Mesher->GetComputedMaxArea();

    //update the ground mesh
    mesherResult = vtkPolyData::New();
    mesherResult->ShallowCopy(this->Mesher->GetOutput());
  }
  else
  {
    return false;
  }

  if (validMesh && mesherResult)
  {
    //we now have to call triangle with each arc sets updated
    //boundary that the ground plane generated.

    //Construct the TerrainInfo first
    this->TerrainInfo =
      new vtkCmbInternalTerrainInfo(inputVector[0], this->ElevationRadius, this->VOIBounds[4]);

    vtkPolyData* finalMesh = vtkPolyData::New();
    finalMesh->DeepCopy(mesherResult);
    int numOfArrays = finalMesh->GetCellData()->GetNumberOfArrays();
    for (int j = numOfArrays - 1; j >= 0; --j)
    {
      finalMesh->GetCellData()->RemoveArray(j);
    }
    finalMesh->SetLines(NULL); //don't want lines in the output

    validMesh = this->GenerateGroundMesh(finalMesh);
    if (validMesh && this->PolygonInfo->Info.size() > 1)
    {
      //At this point mesherResult and finalMesh are the same except
      //final mesh doesn't have lines, arrays, and has extruded points
      this->GenerateExtrudedArcSets(mesherResult, finalMesh);
    }
    mesherResult->Delete();

    delete this->TerrainInfo;
    this->TerrainInfo = NULL;

    vtkPolyData* output = vtkPolyData::GetData(outputVector);
    output->ShallowCopy(finalMesh);
    finalMesh->Delete();
  }

  //force the progress to finish no matter what.
  this->UpdateProgress(1.0);
  return validMesh;
}

bool vtkCMBMeshTerrainWithArcs::PrepForMeshing(vtkInformationVector* input, vtkPolyData* mesh)
{
  //make sure we have no verts, triangles, and strips
  //reduces possible memory foot print
  mesh->SetVerts(NULL);
  mesh->SetPolys(NULL);
  mesh->SetStrips(NULL);

  //clear the information each time so we don't
  //have information of the previous modified run
  this->PolygonInfo->Info.clear();

  vtkCellArray* newCells = vtkCellArray::New();

  vtkCellArray* cells = NULL;
  int numInputs = input->GetNumberOfInformationObjects();
  for (int idx = 0; idx < numInputs; ++idx)
  {
    cells = NULL;
    vtkPolyData* pd = vtkPolyData::GetData(input, idx);
    cells = pd->GetLines();
    if (!cells)
    {
      continue;
    }

    if (cells->GetNumberOfCells() > 0)
    {
      PInfo info;
      int numCells = 0;
      vtkIdType npts, *pts;
      bool firstCell = true;
      for (cells->InitTraversal(); cells->GetNextCell(npts, pts);)
      {
        if (firstCell)
        {
          //in the first phase  we can presume the contour has a constant z value
          //across all points
          info.Height = pd->GetPoint(pts[0])[2];
          firstCell = false;
        }
        for (vtkIdType i = 0; i < (npts - 1); ++i, ++numCells)
        {
          newCells->InsertNextCell(2, pts + i);
        }
      }

      //all Ids have to be 1 based
      info.Id = static_cast<int>(this->PolygonInfo->Info.size()) + 1;
      info.NumberOfCells = numCells;
      this->PolygonInfo->Info.push_back(info);
      this->NextProgressStep();
    }
  }

  //combine all the arcs together into a single polydata
  //so that we can get the combined point set
  vtkNew<vtkPoints> points;
  const int numInputsToAppend = input->GetNumberOfInformationObjects();
  if (numInputsToAppend > 0)
  {
    vtkNew<vtkAppendPolyData> append;
    for (int i = 0; i < numInputsToAppend; ++i)
    {
      append->AddInputData(vtkPolyData::GetData(input, i));
    }
    append->Update();
    //inc ref count so it doesn't go to null when append scope deletes
    points->ShallowCopy(append->GetOutput()->GetPoints());
  }

  //Insert the ground plane
  this->InsertGroundPlane(points.GetPointer(), newCells);

  vtkIdType size = newCells->GetNumberOfCells();
  mesh->SetLines(newCells);
  newCells->FastDelete();

  //bitwise and each method call to get if the mesh is valid
  bool validMesh = true;
  //Setup the ids in the cell data
  validMesh &= this->AssignPolygonIds(mesh, size);
  this->NextProgressStep();

  //when the meshing happens it will require all the points have 0 zvalues
  validMesh &= this->FlattenMesh(points.GetPointer());
  this->NextProgressStep();

  mesh->SetPoints(points.GetPointer());

  return validMesh;
}

bool vtkCMBMeshTerrainWithArcs::AssignPolygonIds(vtkPolyData* mesh, const vtkIdType& /*size*/) const
{
  typedef vtkCMBPrepareForTriangleMesher vtkPrepareForMesher;
  vtkNew<vtkPrepareForMesher> mapInterface;
  mapInterface->SetPolyData(mesh); //Add necessary field data to the mesh
  mapInterface->SetNumberOfArcs(this->PolygonInfo->Info.size());
  mapInterface->SetNumberOfLoops(this->PolygonInfo->Info.size());
  mapInterface->InitializeNewMapInfo();

  size_t i = 0;
  vtkIdType index = 0;
  int numCells = 0;
  //set up the cell arrays
  vtkIdType groundId = this->PolygonInfo->Info[this->PolygonInfo->GroundIndex].Id;
  int asId = 0;
  int endpointId = 0;
  for (i = 0; i < this->PolygonInfo->Info.size(); ++i)
  {
    numCells = this->PolygonInfo->Info[i].NumberOfCells;
    asId = this->PolygonInfo->Info[i].Id;
    if (asId == groundId)
    {
      //In this program a loop is just one arc
      vtkIdType loopId = mapInterface->AddLoop(groundId, -1);
      //endpoints don't mean anything, just make sure they are unique
      //arcs will never be part of more than one loop in this code
      mapInterface->AddArc(index * 3, numCells * 3, groundId, loopId, -1, endpointId, endpointId);
      endpointId++;

      index += numCells;
    }
    else
    {
      //This loop will be a hole of the ground
      vtkIdType loopId = mapInterface->AddLoop(-1, groundId);
      mapInterface->AddArc(index * 3, numCells * 3, asId, loopId, -1, endpointId, endpointId);
      endpointId++;

      index += numCells;
    }
  }
  mapInterface->FinalizeNewMapInfo(); //add the field data

  return true;
}

bool vtkCMBMeshTerrainWithArcs::FlattenMesh(vtkPoints* points) const
{
  bool valid = true;
  vtkDataArray* dataArray = points->GetData();
  int size = points->GetNumberOfPoints() * 3;
  if (dataArray->GetDataType() == VTK_FLOAT)
  {
    vtkFloatArray* floatArray = static_cast<vtkFloatArray*>(dataArray);
    float* pt = floatArray->GetPointer(0);
    zeroZValues(pt, size);
  }
  else if (dataArray->GetDataType() == VTK_DOUBLE)
  {
    vtkDoubleArray* doubleArray = static_cast<vtkDoubleArray*>(dataArray);
    double* pt = doubleArray->GetPointer(0);
    zeroZValues(pt, size);
  }
  else
  {
    valid = false;
  }
  return valid;
}

void vtkCMBMeshTerrainWithArcs::InsertGroundPlane(vtkPoints* points, vtkCellArray* newCells) const
{
  vtkIdType bl, br, tr, tl;
  bl = points->InsertNextPoint(this->VOIBounds[0], this->VOIBounds[2], 0); //bottom left
  br = points->InsertNextPoint(this->VOIBounds[1], this->VOIBounds[2], 0); //bottom right
  tr = points->InsertNextPoint(this->VOIBounds[1], this->VOIBounds[3], 0); //top right
  tl = points->InsertNextPoint(this->VOIBounds[0], this->VOIBounds[3], 0); //top left
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(br);
  newCells->InsertCellPoint(bl);
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(bl);
  newCells->InsertCellPoint(tl);
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(tl);
  newCells->InsertCellPoint(tr);
  newCells->InsertNextCell(2);
  newCells->InsertCellPoint(tr);
  newCells->InsertCellPoint(br);

  PInfo info;
  //all Ids have to be 1 based
  info.Id = static_cast<int>(this->PolygonInfo->Info.size()) + 1;
  info.NumberOfCells = 4;
  this->PolygonInfo->Info.push_back(info);
  this->PolygonInfo->GroundIndex = static_cast<int>(this->PolygonInfo->Info.size() - 1);
}

bool vtkCMBMeshTerrainWithArcs::GenerateGroundMesh(vtkPolyData* finalMesh)
{
  bool valid = this->ExtrudeMeshPoints(finalMesh->GetPoints());
  this->NextProgressStep();
  return valid;
}

bool vtkCMBMeshTerrainWithArcs::GenerateExtrudedArcSets(vtkPolyData* input, vtkPolyData* output)
{
  bool valid;
  for (size_t i = 0; i < this->PolygonInfo->Info.size(); i++)
  {
    //skip the ground plane
    if (static_cast<size_t>(this->PolygonInfo->GroundIndex) == i)
    {
      continue;
    }
    vtkNew<vtkPolyData> arcsetMesh;
    int id = this->PolygonInfo->Info[i].Id;
    this->CreateArcSetForMeshing(input, arcsetMesh.GetPointer(), id);
    this->NextProgressStep();

    this->Mesher->SetPreserveBoundaries(true);
    this->Mesher->SetMaxAreaMode(vtkCMBTriangleMesher::AbsoluteArea);
    this->Mesher->SetMaxArea(this->MesherMaxArea);
    this->Mesher->SetInputData(arcsetMesh.GetPointer());
    this->Mesher->Update();
    this->NextProgressStep();

    vtkNew<vtkPolyData> resultMesh;
    resultMesh->ShallowCopy(this->Mesher->GetOutput());

    //remove all the cell properties since they are uneeded
    int numOfArrays = resultMesh->GetCellData()->GetNumberOfArrays();
    for (int j = numOfArrays - 1; j >= 0; --j)
    {
      resultMesh->GetCellData()->RemoveArray(j);
    }

    valid = this->ExtrudeArcMesh(resultMesh.GetPointer(), static_cast<int>(i));
    this->NextProgressStep();
    if (valid)
    {
      //now add this mesh on the final mesh
      vtkNew<vtkAppendPolyData> append;
      append->AddInputData(output);
      append->AddInputData(resultMesh.GetPointer());
      append->Update();
      vtkPolyData* result = append->GetOutput();
      output->ShallowCopy(result);
    }
  }

  return true;
}

void vtkCMBMeshTerrainWithArcs::CreateArcSetForMeshing(
  vtkPolyData* input, vtkPolyData* output, const int& index) const
{
  vtkCellData* cellData = input->GetCellData();
  vtkPoints* points = input->GetPoints();
  if (cellData == NULL || points == NULL)
  {
    return;
  }

  vtkPoints* newPoints = vtkPoints::New();
  newPoints->DeepCopy(points);
  output->SetPoints(newPoints);
  newPoints->Delete();

  vtkCellArray* newCells = vtkCellArray::New();
  newCells->Allocate(input->GetNumberOfLines());
  //ElementIds will be the same as the arcSetIds
  vtkIdTypeArray* elementIds = vtkIdTypeArray::SafeDownCast(cellData->GetArray("ElementIds"));
  vtkIdType numCells = input->GetNumberOfCells();

  for (vtkIdType i = 0; i < numCells; ++i)
  {
    if ((input->GetCellType(i) == VTK_LINE || input->GetCellType(i) == VTK_POLY_LINE) &&
      elementIds->GetValue(i) == index)
    {
      newCells->InsertNextCell(input->GetCell(i));
    }
  }

  output->SetLines(newCells);
  vtkIdType size = output->GetNumberOfCells();
  newCells->FastDelete();

  typedef vtkCMBPrepareForTriangleMesher vtkPrepareForMesher;
  vtkNew<vtkPrepareForMesher> mapInterface;
  mapInterface->SetPolyData(output); //Add necessary field data to the mesh
  mapInterface->SetNumberOfArcs(1);
  mapInterface->SetNumberOfLoops(1);
  mapInterface->InitializeNewMapInfo();

  //There is only 1 arc and only 1 loop
  vtkIdType loopId = mapInterface->AddLoop(index, -1);
  mapInterface->AddArc(0, size * 3, 1, loopId, -1, 1, 1);

  mapInterface->FinalizeNewMapInfo();
}

bool vtkCMBMeshTerrainWithArcs::ExtrudeArcMesh(vtkPolyData* mesh, const int& index) const
{
  //note this method modifies the mesh
  if (mesh == NULL)
  {
    return false;
  }
  vtkPoints* points = mesh->GetPoints();
  if (points == NULL)
  {
    return false;
  }

  //restore this arc set to the correct height
  double height = this->PolygonInfo->Info[index].Height;
  this->UnFlattenMesh(points, height);

  //build up the cell links, so we can insert new cells
  mesh->BuildCells();

  //now convert the lines into polygons that form the sides of the
  //extrusion
  std::map<vtkIdType, vtkIdType> SidesBottomToTop;
  vtkIdType meshId, newGroundId, npts;
  vtkIdList* pIds = vtkIdList::New();
  double p[3];

  vtkCellArray* lines = mesh->GetLines();
  lines->InitTraversal();
  while (lines->GetNextCell(pIds))
  {
    npts = pIds->GetNumberOfIds();
    for (vtkIdType i = npts - 1; i >= 0; --i)
    {
      meshId = pIds->GetId(i);
      if (SidesBottomToTop.find(meshId) == SidesBottomToTop.end())
      {
        points->GetPoint(meshId, p);
        p[2] = 0; //set it to zero so get the right info
        p[2] = this->TerrainInfo->getElevation(p);
        newGroundId = points->InsertNextPoint(p);
      }
      else
      {
        newGroundId = SidesBottomToTop.find(meshId)->second;
      }
      SidesBottomToTop.insert(std::pair<vtkIdType, vtkIdType>(meshId, newGroundId));
      pIds->InsertNextId(newGroundId);
    }
    mesh->InsertNextCell(VTK_POLYGON, pIds);
    pIds->Reset();
  }

  pIds->Delete();

  //now the we are done iterating the lines, delete them as they are not needed
  mesh->SetLines(NULL);

  return true;
}

bool vtkCMBMeshTerrainWithArcs::UnFlattenMesh(vtkPoints* points, const double& height) const
{
  //get the point data void pointer
  vtkDataArray* dataArray = points->GetData();
  vtkIdType size = points->GetNumberOfPoints();
  if (dataArray->GetDataType() == VTK_FLOAT)
  {
    vtkFloatArray* floatArray = static_cast<vtkFloatArray*>(dataArray);
    float* pos = floatArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = height;
    }
  }
  else
  {
    vtkDoubleArray* doubleArray = static_cast<vtkDoubleArray*>(dataArray);
    double* pos = doubleArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = height;
    }
  }
  return true;
}

bool vtkCMBMeshTerrainWithArcs::ExtrudeMeshPoints(vtkPoints* points) const
{
  //get the point data void pointer
  vtkDataArray* dataArray = points->GetData();
  vtkIdType size = points->GetNumberOfPoints();
  if (dataArray->GetDataType() == VTK_FLOAT)
  {
    vtkFloatArray* floatArray = static_cast<vtkFloatArray*>(dataArray);
    float* pos = floatArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = this->TerrainInfo->getElevation(pos);
    }
  }
  else
  {
    vtkDoubleArray* doubleArray = static_cast<vtkDoubleArray*>(dataArray);
    double* pos = doubleArray->GetPointer(0);
    for (vtkIdType i = 0; i < size; ++i, pos += 3)
    {
      pos[2] = this->TerrainInfo->getElevation(pos);
    }
  }
  return true;
}

void vtkCMBMeshTerrainWithArcs::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "NumberOfProgressSteps: " << this->NumberOfProgressSteps << std::endl;
  os << indent << "CurrentProgressStep: " << this->CurrentProgressStep << std::endl;
  os << indent << "StepIncrement: " << this->StepIncrement << std::endl;
  os << indent << "VOIBounds: " << this->VOIBounds[0] << ", " << this->VOIBounds[1] << ", "
     << this->VOIBounds[2] << ", " << this->VOIBounds[3] << ", " << this->VOIBounds[4] << ", "
     << this->VOIBounds[5] << ", " << std::endl;
  os << indent << "ElevationRadius: " << this->ElevationRadius << std::endl;
  os << indent << "MesherMaxArea: " << this->MesherMaxArea << std::endl;
  os << indent << "Mesher: " << this->Mesher << std::endl;
}
