// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-FileCopyrightText: Copyright 2008 Sandia Corporation
// SPDX-License-Identifier: LicenseRef-BSD-3-Clause-Sandia-NVIDIA-USGov

#include "vtkCoincidentPoints.h"

#include "vtkIdList.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPoints.h"
#include "vtkSmartPointer.h"

#include <map>
#include <set>
#include <vector>

//------------------------------------------------------------------------------
// vtkCoincidentPoints::implementation

VTK_ABI_NAMESPACE_BEGIN
class vtkCoincidentPoints::implementation
{
public:
  implementation() { this->TraversalIterator = this->CoordMap.end(); }

  ~implementation() { this->CoordMap.clear(); }

  struct Coord
  {
    double coord[3];

    Coord()
    {
      this->coord[0] = -1.0;
      this->coord[1] = -1.0;
      this->coord[2] = -1.0;
    }
    Coord(const Coord& src)
    {
      this->coord[0] = src.coord[0];
      this->coord[1] = src.coord[1];
      this->coord[2] = src.coord[2];
    }
    Coord(const double src[3])
    {
      this->coord[0] = src[0];
      this->coord[1] = src[1];
      this->coord[2] = src[2];
    }
    Coord(const double x, const double y, const double z)
    {
      this->coord[0] = x;
      this->coord[1] = y;
      this->coord[2] = z;
    }

    inline bool operator<(const Coord& other) const
    {
      return this->coord[0] < other.coord[0] ||
        (this->coord[0] == other.coord[0] &&
          (this->coord[1] < other.coord[1] ||
            (this->coord[1] == other.coord[1] && this->coord[2] < other.coord[2])));
    }
  };

  typedef std::map<Coord, vtkSmartPointer<vtkIdList>>::iterator MapCoordIter;

  vtkCoincidentPoints* Self;

  std::map<Coord, vtkSmartPointer<vtkIdList>> CoordMap;
  std::map<vtkIdType, vtkIdType> CoincidenceMap;
  MapCoordIter TraversalIterator;
};

//------------------------------------------------------------------------------
// vtkCoincidentPoints

vtkStandardNewMacro(vtkCoincidentPoints);

vtkCoincidentPoints::vtkCoincidentPoints()
{
  this->Implementation = new implementation();
  this->Implementation->Self = this;
}

vtkCoincidentPoints::~vtkCoincidentPoints()
{
  delete this->Implementation;
}

void vtkCoincidentPoints::Clear()
{
  this->Implementation->CoordMap.clear();
  this->Implementation->CoincidenceMap.clear();
}

void vtkCoincidentPoints::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  // os << indent << "MaximumDepth: " << this->MaximumDepth << "\n";
}

void vtkCoincidentPoints::AddPoint(vtkIdType Id, const double point[3])
{
  implementation::Coord coord(point);
  implementation::MapCoordIter mapIter = this->Implementation->CoordMap.find(coord);
  if (mapIter == this->Implementation->CoordMap.end())
  {
    vtkSmartPointer<vtkIdList> idSet = vtkSmartPointer<vtkIdList>::New();
    idSet->InsertNextId(Id);
    this->Implementation->CoordMap[coord] = idSet;
  }
  else
  {
    (*mapIter).second->InsertNextId(Id);
  }
}

vtkIdList* vtkCoincidentPoints::GetCoincidentPointIds(const double point[3])
{
  implementation::Coord coord(point);
  implementation::MapCoordIter mapIter = this->Implementation->CoordMap.find(coord);
  if (mapIter == this->Implementation->CoordMap.end())
  {
    return nullptr;
  }

  if ((*mapIter).second->GetNumberOfIds() > 1)
  {
    return (*mapIter).second;
  }
  else
  {
    // No Coincident Points
    return nullptr;
  }
}

void vtkCoincidentPoints::RemoveNonCoincidentPoints()
{
  implementation::MapCoordIter mapIter = this->Implementation->CoordMap.begin();
  while (mapIter != this->Implementation->CoordMap.end())
  {
    if ((*mapIter).second->GetNumberOfIds() <= 1)
    {
      this->Implementation->CoordMap.erase(mapIter++);
    }
    else
    {
      ++mapIter;
    }
  }
}

vtkIdList* vtkCoincidentPoints::GetNextCoincidentPointIds()
{
  vtkIdList* rvalue = nullptr;
  if (this->Implementation->TraversalIterator != this->Implementation->CoordMap.end())
  {
    rvalue = (*this->Implementation->TraversalIterator).second;
    ++this->Implementation->TraversalIterator;
    return rvalue;
  }

  return rvalue;
}

void vtkCoincidentPoints::InitTraversal()
{
  this->Implementation->TraversalIterator = this->Implementation->CoordMap.begin();
}

//------------------------------------------------------------------------------
// vtkSpiralkVertices - calculate points at a regular interval along a parametric
// spiral.
void vtkCoincidentPoints::SpiralPoints(vtkIdType num, vtkPoints* offsets)
{
  int maxIter = 10;
  double pi = vtkMath::Pi();
  double a = 1 / (4 * pi * pi);
  offsets->Initialize();
  offsets->SetNumberOfPoints(num);

  for (vtkIdType i = 0; i < num; i++)
  {
    double d = 2.0 * i / sqrt(3.0);
    // We are looking for points at regular intervals along the parametric spiral
    // x = t*cos(2*pi*t)
    // y = t*sin(2*pi*t)
    // We cannot solve this equation exactly, so we use newton's method.
    // Using an Excel trendline, we find that
    // t = 0.553*d^0.502
    // is an excellent starting point./g
    double t = 0.553 * pow(d, 0.502);
    for (int iter = 0; iter < maxIter; iter++)
    {
      double r = sqrt(t * t + a * a);
      double f = pi * (t * r + a * a * log(t + r)) - d;
      double df = 2 * pi * r;
      t = t - f / df;
    }
    double x = t * cos(2 * pi * t);
    double y = t * sin(2 * pi * t);
    offsets->SetPoint(i, x, y, 0);
  }
}
VTK_ABI_NAMESPACE_END
