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

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

#include "vtkAMRInterpolatedVelocityField.h"
#include "vtkAppendPolyData.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkCompositeDataIterator.h"
#include "vtkCompositeDataPipeline.h"
#include "vtkCompositeDataSet.h"
#include "vtkDataSetAttributes.h"
#include "vtkDoubleArray.h"
#include "vtkExecutive.h"
#include "vtkGenericCell.h"
#include "vtkIdList.h"
#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkCellLocator.h"
#include "vtkModifiedBSPTree.h"
#include "vtkInterpolatedVelocityField.h"
#include "vtkAbstractInterpolatedVelocityField.h"
#include "vtkCellLocatorInterpolatedVelocityField.h"
#include "vtkMath.h"
#include "vtkMultiBlockDataSet.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkOverlappingAMR.h"
#include "vtkPointData.h"
#include "vtkPointSet.h"
#include "vtkPolyData.h"
#include "vtkPolyDataCollection.h"
#include "vtkPolyLine.h"
#include "vtkRungeKutta2.h"
#include "vtkRungeKutta4.h"
#include "vtkRungeKutta45.h"
#include "vtkSmartPointer.h"
#include "vtkStreamTracer.h"

#include <algorithm>
#include <array>
#include <iostream>
#include <iterator>
#include <vector>

vtkObjectFactoryNewMacro(vtkEvenlySpacedStreamlines2D);
vtkCxxSetObjectMacro(vtkEvenlySpacedStreamlines2D,Integrator,
                     vtkInitialValueProblemSolver);
vtkCxxSetObjectMacro(vtkEvenlySpacedStreamlines2D,InterpolatorPrototype,
                     vtkAbstractInterpolatedVelocityField);

vtkEvenlySpacedStreamlines2D::vtkEvenlySpacedStreamlines2D()
{
  this->Integrator = vtkRungeKutta2::New();
  for(int i=0; i<3; i++)
  {
    this->StartPosition[i] = 0.0;
  }

  this->IntegrationStepUnit    = vtkStreamTracer::CELL_LENGTH_UNIT;
  this->InitialIntegrationStep = 0.5;
  this->ClosedLoopMaximumDistance = 1.0e-6;
  this->ClosedLoopMaximumDistanceArcLength = 1.0e-6;
  this->LoopAngle = 0.349066; // 20 degrees in radians
  this->MaximumNumberOfSteps = 2000;
  this->MinimumNumberOfLoopPoints = 4;
  this->DirectionStart = 0;
  // invalid integration direction so that we trigger a change the first time
  this->PreviousDirection = 0;

  this->TerminalSpeed        = 1.0E-12;

  this->ComputeVorticity = true;

  this->InterpolatorPrototype = nullptr;

  // by default process active point vectors
  this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
                               vtkDataSetAttributes::VECTORS);
  this->SeparatingDistance = 1;
  this->SeparatingDistanceArcLength = 1;
  this->SeparatingDistanceRatio = 0.5;
  this->SuperposedGrid = vtkImageData::New();
  this->Streamlines = vtkPolyDataCollection::New();
  // by default process active point vectors
  this->SetInputArrayToProcess(0,0,0,vtkDataObject::FIELD_ASSOCIATION_POINTS,
                               vtkDataSetAttributes::VECTORS);
}

vtkEvenlySpacedStreamlines2D::~vtkEvenlySpacedStreamlines2D()
{
  this->SetIntegrator(nullptr);
  this->SetInterpolatorPrototype(nullptr);
  this->SuperposedGrid->Delete();
  this->Streamlines->Delete();
}

int vtkEvenlySpacedStreamlines2D::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  if (!this->SetupOutput(inInfo, outInfo))
  {
    return 0;
  }
  double bounds[6];
  vtkEvenlySpacedStreamlines2D::GetBounds(this->InputData, bounds);
  if (bounds[5] != bounds[4])
  {
    this->InputData->UnRegister(this);
    vtkErrorMacro(
      "vtkEvenlySpacedStreamlines2D does not support planes not alligned with XY.");
    return 0;
  }
  std::array<double, 3> v = {{bounds[1] - bounds[0], bounds[3] - bounds[2],
                              bounds[5] - bounds[4]}};
  double length = vtkMath::Norm(&v[0]);

  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  // compute the separating distance arc length
  double cellLength = 0;
  if (! this->ComputeCellLength(&cellLength))
  {
    this->InputData->UnRegister(this);
    return 0;
  }
  this->SeparatingDistanceArcLength = this->ConvertToLength(
    this->SeparatingDistance, this->IntegrationStepUnit, cellLength);
  this->ClosedLoopMaximumDistanceArcLength = this->ConvertToLength(
    this->ClosedLoopMaximumDistance, this->IntegrationStepUnit, cellLength);
  this->InitializeSuperposedGrid(bounds);
  auto streamTracer = vtkSmartPointer<vtkStreamTracer>::New();
  streamTracer->SetInputDataObject(this->InputData);
  streamTracer->SetMaximumPropagation(length);
  streamTracer->SetMaximumNumberOfSteps(this->MaximumNumberOfSteps);
  streamTracer->SetIntegrationDirection(vtkStreamTracer::BOTH);
  streamTracer->SetInputArrayToProcess(0, this->GetInputArrayInformation(0));
  streamTracer->SetStartPosition(this->StartPosition);
  streamTracer->SetTerminalSpeed(this->TerminalSpeed);
  streamTracer->SetInitialIntegrationStep(this->InitialIntegrationStep);
  streamTracer->SetIntegrationStepUnit(this->IntegrationStepUnit);
  streamTracer->SetIntegrator(this->Integrator);
  streamTracer->SetComputeVorticity(this->ComputeVorticity);
  streamTracer->SetInterpolatorPrototype(this->InterpolatorPrototype);
  // we end streamlines after one loop iteration
  streamTracer->AddCustomTerminationCallback(
    &vtkEvenlySpacedStreamlines2D::IsStreamlineLooping, this,
    vtkStreamTracer::FIXED_REASONS_FOR_TERMINATION_COUNT);
  streamTracer->Update();

  auto streamline = vtkSmartPointer<vtkPolyData>::New();
  streamline->ShallowCopy(streamTracer->GetOutput());
  this->AddToAllPoints(streamline);

  auto append = vtkSmartPointer<vtkAppendPolyData>::New();
  append->UserManagedInputsOn();
  append->SetNumberOfInputs(2);
  output->ShallowCopy(streamline);
  int currentSeedId = 1;
  int processedSeedId = 0;

  this->Streamlines->RemoveAllItems();
  this->Streamlines->AddItem(streamline);
  // we also end streamlines when they are close to other streamlines
  streamTracer->AddCustomTerminationCallback(
    &vtkEvenlySpacedStreamlines2D::IsStreamlineTooCloseToOthers, this,
    vtkStreamTracer::FIXED_REASONS_FOR_TERMINATION_COUNT + 1);

  const char* velocityName = this->GetInputArrayToProcessName();
  double deltaOne = this->SeparatingDistanceArcLength / 1000;
  double delta[3] = {deltaOne, deltaOne, deltaOne};
  while (this->Streamlines->GetNumberOfItems())
  {
    streamline = vtkPolyData::SafeDownCast(this->Streamlines->GetItemAsObject(0));
    vtkDataArray* normals = streamline->GetPointData()->GetArray("Normals");
    vtkDataArray* velocity = streamline->GetPointData()->GetArray(velocityName);
    for (vtkIdType pointId = 0; pointId < streamline->GetNumberOfPoints(); ++pointId)
    {
      // generate 2 new seeds for every streamline point
      double newSeedVector[3];
      vtkMath::Cross(normals->GetTuple(pointId), velocity->GetTuple(pointId),
                     newSeedVector);
      // floating point errors move newSeedVector out of XY plane.
      newSeedVector[2] = 0;
      vtkMath::Normalize(newSeedVector);
      vtkMath::MultiplyScalar(newSeedVector, this->SeparatingDistanceArcLength);
      double point[3];
      streamline->GetPoint(pointId, point);
      std::array<std::array<double, 3>, 2> newSeeds;
      vtkMath::Add(point, newSeedVector, &newSeeds[0][0]);
      vtkMath::Subtract(point, newSeedVector, &newSeeds[1][0]);

      for (auto newSeed: newSeeds)
      {
        if (vtkMath::PointIsWithinBounds(&newSeed[0], bounds, delta) &&
            ! this->ForEachCell(&newSeed[0],
                                &vtkEvenlySpacedStreamlines2D::IsTooClose<DISTANCE>))
        {
          streamTracer->SetStartPosition(&newSeed[0]);
          streamTracer->Update();
          auto newStreamline = vtkSmartPointer<vtkPolyData>::New();
          newStreamline->ShallowCopy(streamTracer->GetOutput());

          vtkIntArray* seedIds = vtkIntArray::SafeDownCast(
            newStreamline->GetCellData()->GetArray("SeedIds"));
          for (int cellId = 0; cellId < newStreamline->GetNumberOfCells(); ++cellId)
          {
            seedIds->SetValue(cellId, currentSeedId);
          }
          currentSeedId++;
          this->AddToAllPoints(newStreamline);
          append->SetInputDataByNumber(0, output);
          append->SetInputDataByNumber(1, newStreamline);
          append->Update();
          output->ShallowCopy(append->GetOutput());
          this->Streamlines->AddItem(newStreamline);
        }
      }
    }
    this->Streamlines->RemoveItem(0);
    ++processedSeedId;
  }
  this->InputData->UnRegister(this);
  return 1;
}

int vtkEvenlySpacedStreamlines2D::ComputeCellLength(
  double* cellLength)
{
  vtkAbstractInterpolatedVelocityField* func;
  int maxCellSize = 0;
  if (this->CheckInputs(func, &maxCellSize) != VTK_OK)
  {
      if(func)
      {
        func->Delete();
      }
      return 0;
  }
  vtkDataSet* input;
  auto cell = vtkSmartPointer<vtkGenericCell>::New();
  double velocity[3];
  // access the start postion
  if (!func->FunctionValues(this->StartPosition, velocity))
  {
    func->Delete();
    return 0;
  }
  // Make sure we use the dataset found by the vtkAbstractInterpolatedVelocityField
  input = func->GetLastDataSet();
  input->GetCell(func->GetLastCellId(), cell);
  *cellLength = sqrt(static_cast<double>(cell->GetLength2()));
  func->Delete();
  return 1;
}

int vtkEvenlySpacedStreamlines2D::FillInputPortInformation(
  int port, vtkInformation *info)
{
  if (port == 0)
  {
    info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataObject");
  }
  return 1;
}

bool vtkEvenlySpacedStreamlines2D::IsStreamlineTooCloseToOthers(
  void* clientdata,
  vtkPoints* points, vtkDataArray* velocity, int direction)
{
  (void)velocity;
  (void) direction;
  vtkEvenlySpacedStreamlines2D* This =
    static_cast<vtkEvenlySpacedStreamlines2D*>(clientdata);
  vtkIdType count = points->GetNumberOfPoints();
  double point[3];
  points->GetPoint(count - 1, point);
  return This->ForEachCell(
    point, &vtkEvenlySpacedStreamlines2D::IsTooClose<DISTANCE_RATIO>);
}

bool vtkEvenlySpacedStreamlines2D::IsStreamlineLooping(
  void* clientdata,
  vtkPoints* points, vtkDataArray* velocity, int direction)
{
  vtkEvenlySpacedStreamlines2D* This =
    static_cast<vtkEvenlySpacedStreamlines2D*>(clientdata);
  vtkIdType p0 = points->GetNumberOfPoints() - 1;

  // reinitialize when changing direction
  if (direction != This->PreviousDirection)
  {
    This->InitializePoints(This->CurrentPoints);
    This->InitializeMinPointIds();
    This->PreviousDirection = direction;
    This->DirectionStart = p0;
  }

  double p0Point[3];
  points->GetPoint(p0, p0Point);
  int ijk[3] = {0, 0, 0};
  ijk[0] = floor(p0Point[0] / This->SeparatingDistanceArcLength);
  ijk[1] = floor(p0Point[1] / This->SeparatingDistanceArcLength);
  vtkIdType cellId = This->SuperposedGrid->ComputeCellId(&ijk[0]);

  bool retVal = This->ForEachCell(
    p0Point, &vtkEvenlySpacedStreamlines2D::IsLooping, points, velocity, direction);

  // add the point to the list
  This->CurrentPoints[cellId].push_back(p0);
  if (p0 < This->MinPointIds[cellId])
  {
    This->MinPointIds[cellId] = p0;
  }
  return retVal;
}

template<typename CellCheckerType>
bool vtkEvenlySpacedStreamlines2D::ForEachCell(
  double* point, CellCheckerType checker,
  vtkPoints* points, vtkDataArray* velocity, int direction)
{
  // point current cell
  int ijk[3] = {0, 0, 0};
  ijk[0] = floor(point[0] / this->SeparatingDistanceArcLength);
  ijk[1] = floor(point[1] / this->SeparatingDistanceArcLength);
  vtkIdType cellId = this->SuperposedGrid->ComputeCellId(&ijk[0]);
  if ((this->*checker)(point, cellId, points, velocity, direction))
  {
    return true;
  }
  // and check cells around the current cell
  std::array<std::array<int, 3>, 8> around =
  {{
      {{ijk[0] - 1, ijk[1] + 1, ijk[2]}},
      {{ijk[0],     ijk[1] + 1, ijk[2]}},
      {{ijk[0] + 1, ijk[1] + 1, ijk[2]}},
      {{ijk[0] - 1, ijk[1],     ijk[2]}},
      {{ijk[0] + 1, ijk[1],     ijk[2]}},
      {{ijk[0] - 1, ijk[1] - 1, ijk[2]}},
      {{ijk[0],     ijk[1] - 1, ijk[2]}},
      {{ijk[0] + 1, ijk[1] - 1, ijk[2]}}
  }};
  int extent[6];
  this->SuperposedGrid->GetExtent(extent);
  for (auto cellPos: around)
  {
    cellId = this->SuperposedGrid->ComputeCellId(&cellPos[0]);
    if (cellPos[0] >= extent[0] && cellPos[0] < extent[1] &&
        cellPos[1] >= extent[2] && cellPos[1] < extent[3] &&
        (this->*checker)(point, cellId, points, velocity, direction))
    {
      return true;
    }
  }
  return false;
}


bool vtkEvenlySpacedStreamlines2D::IsLooping(
  double* point, vtkIdType cellId,
  vtkPoints* points, vtkDataArray* velocity, int direction)
{
  (void)point;
  // do we have enough points to form a loop
  vtkIdType p0 = points->GetNumberOfPoints() - 1;
  vtkIdType minLoopPoints = std::max(vtkIdType(3), this->MinimumNumberOfLoopPoints);
  if (this->CurrentPoints[cellId].size() > 0 &&
      p0 - this->MinPointIds[cellId] + 1 >= minLoopPoints)
  {
    vtkIdType p1 = p0 - 1;
    double testDistance2 = this->SeparatingDistanceArcLength * this->SeparatingDistanceArcLength *
      this->SeparatingDistanceRatio * this->SeparatingDistanceRatio;
    double maxDistance2 =
      this->ClosedLoopMaximumDistanceArcLength * this->ClosedLoopMaximumDistanceArcLength;
    for (vtkIdType q: this->CurrentPoints[cellId])
    {
      // do we have enough points to form a loop
      if (p0 - q + 1 < minLoopPoints)
      {
        continue;
      }
      double p0Point[3];
      points->GetPoint(p0, p0Point);
      double qPoint[3];
      points->GetPoint(q, qPoint);
      double distance2 = vtkMath::Distance2BetweenPoints(p0Point, qPoint);
      if (distance2 <= maxDistance2)
      {
        // closed loop
        return true;
      }
      if (distance2 >= testDistance2)
      {
        // we might loop but points are too far.
        continue;
      }
      double p1Point[3];
      points->GetPoint(p1, p1Point);
      double v1[3];
      vtkMath::Subtract(p0Point, p1Point, v1);
      vtkMath::MultiplyScalar(v1, direction);
      double* qVector = velocity->GetTuple(q);
      if (vtkMath::Dot(qVector, v1) < cos(this->LoopAngle))
      {
        // qVector makes a large angle with p0p1
        continue;
      }
      double u0[3], u1[3];
      vtkMath::Subtract(p0Point, qPoint, u0);
      vtkMath::MultiplyScalar(u0, direction);
      vtkMath::Subtract(p1Point, qPoint, u1);
      vtkMath::MultiplyScalar(u1, direction);
      if (vtkMath::Dot(u0, v1) >= 0 && vtkMath::Dot(u1, v1) >=0)
      {
        // we found a "proponent point" See Liu et al.
        continue;
      }
      // the algorithm in Liu at al. has another test that determines if the
      // loop is closed or spiraling. We don't care about that so we skip it.
      return true;
    }
  }
  return false;
}

template<int distanceType>
bool vtkEvenlySpacedStreamlines2D::IsTooClose(
  double* point, vtkIdType cellId,
  vtkPoints* points, vtkDataArray* velocity, int direction)
{
  (void) points;
  (void) velocity;
  (void) direction;
  double testDistance2 = this->SeparatingDistanceArcLength * this->SeparatingDistanceArcLength;
  if (distanceType == DISTANCE_RATIO)
  {
    testDistance2 *= (this->SeparatingDistanceRatio * this->SeparatingDistanceRatio);
  }
  for (auto cellPoint : this->AllPoints[cellId])
  {
    double distance2 = vtkMath::Distance2BetweenPoints(point, &cellPoint[0]);
    if (distance2 < testDistance2)
    {
      return true;
    }
  }
  return false;
}

int vtkEvenlySpacedStreamlines2D::GetIntegratorType()
{
  if (!this->Integrator)
  {
    return vtkStreamTracer::NONE;
  }
  if (!strcmp(this->Integrator->GetClassName(), "vtkRungeKutta2"))
  {
    return vtkStreamTracer::RUNGE_KUTTA2;
  }
  if (!strcmp(this->Integrator->GetClassName(), "vtkRungeKutta4"))
  {
    return vtkStreamTracer::RUNGE_KUTTA4;
  }
  return vtkStreamTracer::UNKNOWN;
}

void vtkEvenlySpacedStreamlines2D::SetInterpolatorTypeToDataSetPointLocator()
{
  this->SetInterpolatorType(
    static_cast<int>(vtkStreamTracer::INTERPOLATOR_WITH_DATASET_POINT_LOCATOR));
}

void vtkEvenlySpacedStreamlines2D::SetInterpolatorTypeToCellLocator()
{
  this->SetInterpolatorType
    (  static_cast<int> ( vtkStreamTracer::INTERPOLATOR_WITH_CELL_LOCATOR )  );
}

void vtkEvenlySpacedStreamlines2D::SetInterpolatorType( int interpType )
{
  if ( interpType == vtkStreamTracer::INTERPOLATOR_WITH_CELL_LOCATOR )
  {
    // create an interpolator equipped with a cell locator
    vtkSmartPointer< vtkCellLocatorInterpolatedVelocityField > cellLoc =
    vtkSmartPointer< vtkCellLocatorInterpolatedVelocityField >::New();

    // specify the type of the cell locator attached to the interpolator
    vtkSmartPointer< vtkModifiedBSPTree > cellLocType =
    vtkSmartPointer< vtkModifiedBSPTree >::New();
    cellLoc->SetCellLocatorPrototype( cellLocType );

    this->SetInterpolatorPrototype( cellLoc );
  }
  else
  {
    // create an interpolator equipped with a point locator (by default)
    vtkSmartPointer< vtkInterpolatedVelocityField > pntLoc =
    vtkSmartPointer< vtkInterpolatedVelocityField >::New();
    this->SetInterpolatorPrototype( pntLoc );
  }
}

void vtkEvenlySpacedStreamlines2D::SetIntegratorType(int type)
{
  vtkInitialValueProblemSolver* ivp=nullptr;
  switch (type)
  {
    case vtkStreamTracer::RUNGE_KUTTA2:
      ivp = vtkRungeKutta2::New();
      break;
    case vtkStreamTracer::RUNGE_KUTTA4:
      ivp = vtkRungeKutta4::New();
      break;
    default:
      vtkWarningMacro("Unrecognized integrator type. Keeping old one.");
      break;
  }
  if (ivp)
  {
    this->SetIntegrator(ivp);
    ivp->Delete();
  }
}

void vtkEvenlySpacedStreamlines2D::SetIntegrationStepUnit( int unit )
{
  if ( unit != vtkStreamTracer::LENGTH_UNIT &&
       unit != vtkStreamTracer::CELL_LENGTH_UNIT )
  {
    unit = vtkStreamTracer::CELL_LENGTH_UNIT;
  }

  if ( unit == this->IntegrationStepUnit )
  {
    return;
  }

  this->IntegrationStepUnit = unit;
  this->Modified();
}

double vtkEvenlySpacedStreamlines2D::ConvertToLength(
  double interval, int unit, double cellLength )
{
  double retVal = 0.0;
  if ( unit == vtkStreamTracer::LENGTH_UNIT )
  {
    retVal = interval;
  }
  else if ( unit == vtkStreamTracer::CELL_LENGTH_UNIT )
  {
    retVal = interval * cellLength;
  }
  return retVal;
}

int vtkEvenlySpacedStreamlines2D::SetupOutput(vtkInformation* inInfo,
                                              vtkInformation* outInfo)
{
  int piece=outInfo->Get(
    vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER());
  int numPieces =
    outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES());

  vtkDataObject* input = inInfo->Get(vtkDataObject::DATA_OBJECT());

  vtkCompositeDataSet *hdInput = vtkCompositeDataSet::SafeDownCast(input);
  vtkDataSet* dsInput = vtkDataSet::SafeDownCast(input);
  if (hdInput)
  {
    this->InputData = hdInput;
    hdInput->Register(this);
    return 1;
  }
  else if (dsInput)
  {
    auto mb = vtkSmartPointer<vtkMultiBlockDataSet>::New();
    mb->SetNumberOfBlocks(numPieces);
    mb->SetBlock(piece, dsInput);
    this->InputData = mb;
    mb->Register(this);
    return 1;
  }
  else
  {
    vtkErrorMacro("This filter cannot handle input of type: "
                  << (input?input->GetClassName():"(none)"));
    return 0;
  }
}


int vtkEvenlySpacedStreamlines2D::CheckInputs(vtkAbstractInterpolatedVelocityField*& func,
                                              int* maxCellSize)
{
  if (!this->InputData)
  {
    return VTK_ERROR;
  }

  vtkOverlappingAMR* amrData = vtkOverlappingAMR::SafeDownCast(this->InputData);

  vtkSmartPointer<vtkCompositeDataIterator> iter;
  iter.TakeReference(this->InputData->NewIterator());

  vtkDataSet* input0 =nullptr;
  iter->GoToFirstItem();
  while (!iter->IsDoneWithTraversal() && input0==nullptr)
  {
    input0 = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
    iter->GoToNextItem();
  }
  if(!input0)
  {
    return VTK_ERROR;
  }

  int vecType(0);
  vtkDataArray *vectors = this->GetInputArrayToProcess(0,input0,vecType);
  if (!vectors)
  {
    return VTK_ERROR;
  }

  // Set the function set to be integrated
  if ( !this->InterpolatorPrototype )
  {
    if(amrData)
    {
      func = vtkAMRInterpolatedVelocityField::New();
    }
    else
    {
      func = vtkInterpolatedVelocityField::New();
    }
  }
  else
  {
    if(amrData && vtkAMRInterpolatedVelocityField::SafeDownCast(this->InterpolatorPrototype)==nullptr)
    {
      this->InterpolatorPrototype = vtkAMRInterpolatedVelocityField::New();
    }
    func = this->InterpolatorPrototype->NewInstance();
    func->CopyParameters(this->InterpolatorPrototype);
  }

  if(vtkAMRInterpolatedVelocityField::SafeDownCast(func))
  {
    assert(amrData);
    vtkAMRInterpolatedVelocityField::SafeDownCast(func)->SetAMRData(amrData);
    if(maxCellSize)
    {
      *maxCellSize = 8;
    }
  }
  else if(vtkCompositeInterpolatedVelocityField::SafeDownCast(func))
  {
    iter->GoToFirstItem();
    while (!iter->IsDoneWithTraversal())
    {
      vtkDataSet* inp = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
      if (inp)
      {
        int cellSize = inp->GetMaxCellSize();
        if ( cellSize > *maxCellSize )
        {
          *maxCellSize = cellSize;
        }
        vtkCompositeInterpolatedVelocityField::SafeDownCast(func)->AddDataSet(inp);
      }
      iter->GoToNextItem();
    }
  }
  else
  {
    assert(false);
  }

  const char *vecName = vectors->GetName();
  func->SelectVectors(vecType,vecName);
  return VTK_OK;
}

void vtkEvenlySpacedStreamlines2D::InitializeSuperposedGrid(double* bounds)
{
  this->SuperposedGrid->SetExtent(
    floor(bounds[0]/this->SeparatingDistanceArcLength),
    ceil(bounds[1]/this->SeparatingDistanceArcLength),
    floor(bounds[2]/this->SeparatingDistanceArcLength),
    ceil(bounds[3]/this->SeparatingDistanceArcLength), 0, 0);
  this->SuperposedGrid->SetSpacing(
    this->SeparatingDistanceArcLength, this->SeparatingDistanceArcLength, this->SeparatingDistanceArcLength);
  this->InitializePoints(this->AllPoints);
  this->InitializePoints(this->CurrentPoints);
}

template<typename T>
void vtkEvenlySpacedStreamlines2D::InitializePoints(T& points)
{
  points.resize(this->SuperposedGrid->GetNumberOfCells());
  for (std::size_t i = 0; i < points.size(); ++i)
  {
    points[i].clear();
  }
}

void vtkEvenlySpacedStreamlines2D::InitializeMinPointIds()
{
  this->MinPointIds.resize(this->SuperposedGrid->GetNumberOfCells());
  for (std::size_t i = 0; i < this->MinPointIds.size(); ++i)
  {
    this->MinPointIds[i] = std::numeric_limits<vtkIdType>::max();
  }
}


void vtkEvenlySpacedStreamlines2D::AddToAllPoints(vtkPolyData* streamline)
{
  vtkPoints* points = streamline->GetPoints();
  if (points)
  {
    for (vtkIdType i = 0; i < points->GetNumberOfPoints(); ++i)
    {
      double point[3];
      points->GetPoint(i, point);
      int ijk[3] = {0, 0, 0};
      ijk[0] = floor(point[0] / this->SeparatingDistanceArcLength);
      ijk[1] = floor(point[1] / this->SeparatingDistanceArcLength);
      vtkIdType cellId = this->SuperposedGrid->ComputeCellId(ijk);
      this->AllPoints[cellId].push_back({{point[0], point[1], point[2]}});
    }
  }
}


void vtkEvenlySpacedStreamlines2D::GetBounds(
  vtkCompositeDataSet* cds, double bounds[6])
{
  if(vtkOverlappingAMR::SafeDownCast(cds))
  {
    vtkOverlappingAMR* amr =vtkOverlappingAMR::SafeDownCast(cds);
    amr->GetBounds(bounds);
  }
  else
  {
    // initialize bounds
    for (int i: {0, 2, 4})
    {
      bounds[i] = std::numeric_limits<double>::max();
    }
    for (int i: {1, 3, 5})
    {
      bounds[i] = - std::numeric_limits<double>::max();
    }
    // go over all datasets in the composite data and find min,max
    // for components of all bounds
    vtkSmartPointer<vtkCompositeDataIterator> iter;
    iter.TakeReference (cds->NewIterator());
    iter->GoToFirstItem();
    while (!iter->IsDoneWithTraversal())
    {
      vtkDataSet* input = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
      if (input)
      {
        double b[6];
        input->GetBounds(b);
        for (int i: {0, 2, 4})
        {
          if (b[i] < bounds[i])
          {
            bounds[i] = b[i];
          }
        }
        for (int i: {1, 3, 5})
        {
          if (b[i] > bounds[i])
          {
            bounds[i] = b[i];
          }
        }
      }
      iter->GoToNextItem();
    }
  }
}

const char* vtkEvenlySpacedStreamlines2D::GetInputArrayToProcessName()
{
  vtkSmartPointer<vtkCompositeDataIterator> iter;
  iter.TakeReference(this->InputData->NewIterator());

  vtkDataSet* input0 =nullptr;
  iter->GoToFirstItem();
  while (!iter->IsDoneWithTraversal() && input0==nullptr)
  {
    input0 = vtkDataSet::SafeDownCast(iter->GetCurrentDataObject());
    iter->GoToNextItem();
  }
  if(!input0)
  {
    return "";
  }
  int vecType(0);
  vtkDataArray *vectors = this->GetInputArrayToProcess(0,input0,vecType);
  if (vectors)
  {
     return vectors->GetName();
  }
  else
  {
    vtkErrorMacro(
      "vtkEvenlySpacedStreamlines2D::SetInputArrayToProcess was not called");
    return nullptr;
  }
}


void vtkEvenlySpacedStreamlines2D::SetIntegratorTypeToRungeKutta2()
{this->SetIntegratorType(vtkStreamTracer::RUNGE_KUTTA2);}

void vtkEvenlySpacedStreamlines2D::SetIntegratorTypeToRungeKutta4()
{this->SetIntegratorType(vtkStreamTracer::RUNGE_KUTTA4);}


void vtkEvenlySpacedStreamlines2D::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  os << indent << "Start position: "
     << this->StartPosition[0] << " "
     << this->StartPosition[1] << " "
     << this->StartPosition[2] << endl;
  os << indent << "Terminal speed: " << this->TerminalSpeed << endl;

  os << indent << "Integration step unit: "
     << ( ( this->IntegrationStepUnit == vtkStreamTracer::LENGTH_UNIT )
          ? "length." : "cell length." ) << endl;

  os << indent << "Initial integration step: "
     << this->InitialIntegrationStep << endl;
  os << indent << "Separation distance: "
     << this->SeparatingDistance << endl;

  os << indent << "Integrator: " << this->Integrator << endl;
  os << indent << "Vorticity computation: "
     << (this->ComputeVorticity ? " On" : " Off") << endl;
}
