/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
#ifndef itkVoronoiDiagram2DGenerator_hxx
#define itkVoronoiDiagram2DGenerator_hxx

#include "itkIntTypes.h"
#include <algorithm>
#include "vnl/vnl_sample.h"
#include "itkMath.h"
#include "itkMakeUniqueForOverwrite.h"

namespace itk
{
const double     NUMERIC_TOLERENCE = 1.0e-10;
constexpr double DIFF_TOLERENCE = 0.001;

template <typename TCoordRepType>
VoronoiDiagram2DGenerator<TCoordRepType>::VoronoiDiagram2DGenerator()
  : m_OutputVD(Self::GetOutput())
  , m_BottomSite(nullptr)

{
  m_VorBoundary.Fill(0.0);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::SetRandomSeeds(int num)
{
  PointType curr;

  m_Seeds.clear();
  double ymax{ m_VorBoundary[1] };
  double xmax{ m_VorBoundary[0] };
  for (int i = 0; i < num; ++i)
  {
    curr[0] = (CoordRepType)(vnl_sample_uniform(0, xmax));
    curr[1] = (CoordRepType)(vnl_sample_uniform(0, ymax));
    m_Seeds.push_back(curr);
  }
  m_NumberOfSeeds = num;
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::SetSeeds(int num, SeedsIterator begin)
{
  m_Seeds.clear();
  SeedsIterator ii(begin);
  for (int i = 0; i < num; ++i)
  {
    m_Seeds.push_back(*ii++);
  }
  m_NumberOfSeeds = num;
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::SetBoundary(PointType vorsize)
{
  m_VorBoundary[0] = vorsize[0];
  m_VorBoundary[1] = vorsize[1];
  m_OutputVD->SetBoundary(vorsize);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::SetOrigin(PointType vorsize)
{
  m_Pxmin = vorsize[0];
  m_Pymin = vorsize[1];
  m_OutputVD->SetOrigin(vorsize);
}

template <typename TCoordRepType>
bool
VoronoiDiagram2DGenerator<TCoordRepType>::comp(PointType arg1, PointType arg2)
{
  if (arg1[1] < arg2[1])
  {
    return true;
  }
  else if (arg1[1] > arg2[1])
  {
    return false;
  }
  // arg1[1] == arg2[1]
  else if (arg1[0] < arg2[0])
  {
    return true;
  }
  else if (arg1[0] > arg2[0])
  {
    return false;
  }
  // arg1[0] == arg2[0]
  else
  {
    return false;
  }
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::SortSeeds()
{
  std::sort(m_Seeds.begin(), m_Seeds.end(), comp);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::AddSeeds(int num, SeedsIterator begin)
{
  auto ii(begin);

  for (int i = 0; i < num; ++i)
  {
    m_Seeds.push_back(*ii++);
  }
  m_NumberOfSeeds += num;
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::AddOneSeed(PointType inputSeed)
{
  m_Seeds.push_back(inputSeed);
  ++m_NumberOfSeeds;
}

template <typename TCoordRepType>
auto
VoronoiDiagram2DGenerator<TCoordRepType>::GetSeed(int SeedID) -> PointType
{
  PointType answer;

  answer[0] = m_Seeds[SeedID][0];
  answer[1] = m_Seeds[SeedID][1];
  return answer;
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::GenerateData()
{
  SortSeeds();
  m_OutputVD->SetSeeds(m_NumberOfSeeds, m_Seeds.begin());
  this->GenerateVDFortune();
  this->ConstructDiagram();
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::UpdateDiagram()
{
  this->GenerateData();
}

template <typename TCoordRepType>
bool
VoronoiDiagram2DGenerator<TCoordRepType>::differentPoint(PointType p1, PointType p2)
{
  double diffx = p1[0] - p2[0];
  double diffy = p1[1] - p2[1];

  return ((diffx < -DIFF_TOLERENCE) || (diffx > DIFF_TOLERENCE) || (diffy < -DIFF_TOLERENCE) ||
          (diffy > DIFF_TOLERENCE));
}

template <typename TCoordRepType>
bool
VoronoiDiagram2DGenerator<TCoordRepType>::almostsame(CoordRepType p1, CoordRepType p2)
{
  double diff = p1 - p2;
  bool   save;

  save = ((diff < -DIFF_TOLERENCE) || (diff > DIFF_TOLERENCE));
  return (!save);
}

template <typename TCoordRepType>
unsigned char
VoronoiDiagram2DGenerator<TCoordRepType>::Pointonbnd(int VertID)
{
  PointType currVert = m_OutputVD->GetVertex(VertID);

  if (almostsame(currVert[0], m_Pxmin))
  {
    return 1;
  }
  else if (almostsame(currVert[1], m_Pymax))
  {
    return 2;
  }
  else if (almostsame(currVert[0], m_Pxmax))
  {
    return 3;
  }
  else if (almostsame(currVert[1], m_Pymin))
  {
    return 4;
  }
  else
  {
    return 0;
  }
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::ConstructDiagram()
{
  const auto rawEdges = make_unique_for_overwrite<EdgeInfoDQ[]>(m_NumberOfSeeds);

  m_OutputVD->Reset();

  EdgeInfo currentPtID;
  int      edges = m_OutputVD->EdgeListSize();

  EdgeInfo LRsites;
  for (int i = 0; i < edges; ++i)
  {
    currentPtID = m_OutputVD->GetEdgeEnd(i);
    LRsites = m_OutputVD->GetLine(m_OutputVD->GetEdgeLineID(i));
    rawEdges[LRsites[0]].push_back(currentPtID);
    rawEdges[LRsites[1]].push_back(currentPtID);
    m_OutputVD->AddCellNeighbor(LRsites);
  }

  PointType corner[4];
  int       cornerID[4];

  corner[0][0] = m_Pxmin;
  corner[0][1] = m_Pymin;
  cornerID[0] = m_Nvert;
  ++m_Nvert;
  m_OutputVD->AddVert(corner[0]);
  corner[1][0] = m_Pxmin;
  corner[1][1] = m_Pymax;
  cornerID[1] = m_Nvert;
  ++m_Nvert;
  m_OutputVD->AddVert(corner[1]);
  corner[2][0] = m_Pxmax;
  corner[2][1] = m_Pymax;
  cornerID[2] = m_Nvert;
  ++m_Nvert;
  m_OutputVD->AddVert(corner[2]);
  corner[3][0] = m_Pxmax;
  corner[3][1] = m_Pymin;
  cornerID[3] = m_Nvert;
  ++m_Nvert;
  m_OutputVD->AddVert(corner[3]);

  std::list<EdgeInfo>                    buildEdges;
  typename std::list<EdgeInfo>::iterator BEiter;
  EdgeInfo                               curr;
  EdgeInfo                               curr1;
  EdgeInfo                               curr2;

  unsigned char               frontbnd;
  unsigned char               backbnd;
  std::vector<IdentifierType> cellPoints;
  for (unsigned int i = 0; i < m_NumberOfSeeds; ++i)
  {
    buildEdges.clear();
    curr = rawEdges[i].front();
    rawEdges[i].pop_front();
    buildEdges.push_back(curr);
    EdgeInfo front = curr;
    EdgeInfo back = curr;
    while (!(rawEdges[i].empty()))
    {
      curr = rawEdges[i].front();
      rawEdges[i].pop_front();
      frontbnd = Pointonbnd(front[0]);
      backbnd = Pointonbnd(back[1]);
      if (curr[0] == back[1])
      {
        buildEdges.push_back(curr);
        back = curr;
      }
      else if (curr[1] == front[0])
      {
        buildEdges.push_front(curr);
        front = curr;
      }
      else if (curr[1] == back[1])
      {
        curr1[1] = curr[0];
        curr1[0] = curr[1];
        buildEdges.push_back(curr1);
        back = curr1;
      }
      else if (curr[0] == front[0])
      {
        curr1[0] = curr[1];
        curr1[1] = curr[0];
        buildEdges.push_front(curr1);
        front = curr1;
      }
      else if ((frontbnd != 0) || (backbnd != 0))
      {
        unsigned char cfrontbnd = Pointonbnd(curr[0]);
        unsigned char cbackbnd = Pointonbnd(curr[1]);

        if ((cfrontbnd == backbnd) && (backbnd))
        {
          curr1[0] = back[1];
          curr1[1] = curr[0];
          buildEdges.push_back(curr1);
          buildEdges.push_back(curr);
          back = curr;
        }
        else if ((cbackbnd == frontbnd) && (frontbnd))
        {
          curr1[0] = curr[1];
          curr1[1] = front[0];
          buildEdges.push_front(curr1);
          buildEdges.push_front(curr);
          front = curr;
        }
        else if ((cfrontbnd == frontbnd) && (frontbnd))
        {
          curr1[0] = curr[0];
          curr1[1] = front[0];
          buildEdges.push_front(curr1);
          curr1[0] = curr[1];
          curr1[1] = curr[0];
          buildEdges.push_front(curr1);
          front = curr1;
        }
        else if ((cbackbnd == backbnd) && (backbnd))
        {
          curr1[0] = back[1];
          curr1[1] = curr[1];
          buildEdges.push_back(curr1);
          curr1[0] = curr[1];
          curr1[1] = curr[0];
          buildEdges.push_back(curr1);
          back = curr1;
        }
        else
        {
          rawEdges[i].push_back(curr);
        }
      }
      else
      {
        rawEdges[i].push_back(curr);
      }
    }

    curr = buildEdges.front();
    curr1 = buildEdges.back();
    if (curr[0] != curr1[1])
    {
      frontbnd = Pointonbnd(curr[0]);
      backbnd = Pointonbnd(curr1[1]);
      if ((frontbnd != 0) && (backbnd != 0))
      {
        if (frontbnd == backbnd)
        {
          curr2[0] = curr1[1];
          curr2[1] = curr[0];
          buildEdges.push_back(curr2);
        }
        else if ((frontbnd == backbnd + 1) || (frontbnd == backbnd - 3))
        {
          curr2[0] = cornerID[frontbnd - 1];
          curr2[1] = curr[0];
          buildEdges.push_front(curr2);
          curr2[1] = curr2[0];
          curr2[0] = curr1[1];
          buildEdges.push_front(curr2);
        }
        else if ((frontbnd == backbnd - 1) || (frontbnd == backbnd + 3))
        {
          curr2[0] = cornerID[backbnd - 1];
          curr2[1] = curr[0];
          buildEdges.push_front(curr2);
          curr2[1] = curr2[0];
          curr2[0] = curr1[1];
          buildEdges.push_front(curr2);
        }
        else
        {
          itkDebugMacro("Numerical problem 1" << curr[0] << ' ' << curr1[1]);
        }
      }
    }

    EdgeInfo pp;

    m_OutputVD->ClearRegion(i);

    for (BEiter = buildEdges.begin(); BEiter != buildEdges.end(); ++BEiter)
    {
      pp = *BEiter;
      m_OutputVD->VoronoiRegionAddPointId(i, pp[0]);
    }
    m_OutputVD->BuildEdge(i);
  }
  m_OutputVD->InsertCells();
}

template <typename TCoordRepType>
bool
VoronoiDiagram2DGenerator<TCoordRepType>::right_of(FortuneHalfEdge * el, PointType * p)
{
  FortuneEdge * e = el->m_Edge;
  FortuneSite * topsite = e->m_Reg[1];

  bool right_of_site = (((*p)[0]) > (topsite->m_Coord[0]));

  if (right_of_site && (!(el->m_RorL)))
  {
    return (true);
  }
  if ((!right_of_site) && (el->m_RorL))
  {
    return (false);
  }
  bool above;
  bool fast;
  if (e->m_A == 1.0)
  {
    double dyp = ((*p)[1]) - (topsite->m_Coord[1]);
    double dxp = ((*p)[0]) - (topsite->m_Coord[0]);
    fast = false;
    if (((!right_of_site) && ((e->m_B) < 0.0)) || (right_of_site && ((e->m_B) >= 0.0)))
    {
      above = (dyp >= (e->m_B) * dxp);
      fast = above;
    }
    else
    {
      above = ((((*p)[0]) + ((*p)[1]) * (e->m_B)) > e->m_C);
      if (e->m_B < 0.0)
      {
        above = !above;
      }
      if (!above)
      {
        fast = true;
      }
    }
    if (!fast)
    {
      double dxs = topsite->m_Coord[0] - ((e->m_Reg[0])->m_Coord[0]);
      above = (((e->m_B) * (dxp * dxp - dyp * dyp)) < (dxs * dyp * (1.0 + 2.0 * dxp / dxs + (e->m_B) * (e->m_B))));
      if ((e->m_B) < 0.0)
      {
        above = !above;
      }
    }
  }
  else
  { // e->m_B == 1.0
    double y1 = (e->m_C) - (e->m_A) * ((*p)[0]);
    double t1 = ((*p)[1]) - y1;
    double t2 = ((*p)[0]) - topsite->m_Coord[0];
    double t3 = y1 - topsite->m_Coord[1];
    above = ((t1 * t1) > (t2 * t2 + t3 * t3));
  }
  return (el->m_RorL ? (!above) : above);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::createHalfEdge(FortuneHalfEdge * task, FortuneEdge * e, bool pm)
{
  task->m_Edge = e;
  task->m_RorL = pm;
  task->m_Next = nullptr;
  task->m_Vert = nullptr;
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::PQshowMin(PointType * answer)
{
  while ((m_PQHash[m_PQmin].m_Next) == nullptr)
  {
    m_PQmin += 1;
  }
  (*answer)[0] = m_PQHash[m_PQmin].m_Next->m_Vert->m_Coord[0];
  (*answer)[1] = m_PQHash[m_PQmin].m_Next->m_Ystar;
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::deletePQ(FortuneHalfEdge * task)
{
  FortuneHalfEdge * last;

  if ((task->m_Vert) != nullptr)
  {
    last = &(m_PQHash[PQbucket(task)]);
    while ((last->m_Next) != task)
    {
      last = last->m_Next;
    }
    last->m_Next = (task->m_Next);
    --m_PQcount;
    task->m_Vert = nullptr;
  }
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::deleteEdgeList(FortuneHalfEdge * task)
{
  (task->m_Left)->m_Right = task->m_Right;
  (task->m_Right)->m_Left = task->m_Left;
  task->m_Edge = &(m_DELETED);
}

template <typename TCoordRepType>
int
VoronoiDiagram2DGenerator<TCoordRepType>::PQbucket(FortuneHalfEdge * task)
{
  int bucket;

  bucket = static_cast<int>((task->m_Ystar - m_Pymin) / m_Deltay * m_PQhashsize);
  if (bucket < 0)
  {
    bucket = 0;
  }
  if (bucket >= static_cast<int>(m_PQhashsize))
  {
    bucket = m_PQhashsize - 1;
  }
  if (bucket < static_cast<int>(m_PQmin))
  {
    m_PQmin = bucket;
  }
  return (bucket);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::insertPQ(FortuneHalfEdge * he, FortuneSite * v, double offset)
{
  he->m_Vert = v;
  he->m_Ystar = (v->m_Coord[1]) + offset;
  FortuneHalfEdge * last = &(m_PQHash[PQbucket(he)]);
  FortuneHalfEdge * enext;

  while (((enext = (last->m_Next)) != nullptr) &&
         (((he->m_Ystar) > (enext->m_Ystar)) ||
          ((Math::ExactlyEquals((he->m_Ystar), (enext->m_Ystar))) && ((v->m_Coord[0]) > (enext->m_Vert->m_Coord[0])))))
  {
    last = enext;
  }
  he->m_Next = (last->m_Next);
  last->m_Next = he;
  m_PQcount += 1;
}

template <typename TCoordRepType>
double
VoronoiDiagram2DGenerator<TCoordRepType>::dist(FortuneSite * s1, FortuneSite * s2)
{
  double dx = (s1->m_Coord[0]) - (s2->m_Coord[0]);
  double dy = (s1->m_Coord[1]) - (s2->m_Coord[1]);

  return (std::sqrt(dx * dx + dy * dy));
}

template <typename TCoordRepType>
auto
VoronoiDiagram2DGenerator<TCoordRepType>::ELgethash(int b) -> FortuneHalfEdge *
{
  if ((b < 0) || (b >= static_cast<int>(m_ELhashsize)))
  {
    return (nullptr);
  }
  FortuneHalfEdge * he = m_ELHash[b];
  if (he == nullptr)
  {
    return (he);
  }
  if (he->m_Edge == nullptr)
  {
    return (he);
  }
  if ((he->m_Edge) != (&m_DELETED))
  {
    return (he);
  }
  m_ELHash[b] = nullptr;

  return (nullptr);
}

template <typename TCoordRepType>
auto
VoronoiDiagram2DGenerator<TCoordRepType>::findLeftHE(PointType * p) -> FortuneHalfEdge *
{
  int  i;
  auto bucket = static_cast<int>((((*p)[0]) - m_Pxmin) / m_Deltax * m_ELhashsize);

  if (bucket < 0)
  {
    bucket = 0;
  }
  if (bucket >= static_cast<int>(m_ELhashsize))
  {
    bucket = static_cast<int>(m_ELhashsize) - 1;
  }
  FortuneHalfEdge * he = ELgethash(bucket);
  if (he == nullptr)
  {
    for (i = 1; true; ++i)
    {
      if ((he = ELgethash(bucket - i)) != nullptr)
      {
        break;
      }
      if ((he = ELgethash(bucket + i)) != nullptr)
      {
        break;
      }
    }
  }

  if ((he == (&m_ELleftend)) || ((he != (&m_ELrightend)) && right_of(he, p)))
  {
    do
    {
      he = he->m_Right;
    } while ((he != (&m_ELrightend)) && (right_of(he, p)));
    he = he->m_Left;
  }
  else
  {
    do
    {
      he = he->m_Left;
    } while ((he != (&m_ELleftend)) && (!right_of(he, p)));
  }

  if ((bucket > 0) && (bucket < static_cast<int>(m_ELhashsize) - 1))
  {
    m_ELHash[bucket] = he;
  }
  return (he);
}

template <typename TCoordRepType>
auto
VoronoiDiagram2DGenerator<TCoordRepType>::getRightReg(FortuneHalfEdge * he) -> FortuneSite *
{
  if ((he->m_Edge) == nullptr)
  {
    return (m_BottomSite);
  }
  else if (he->m_RorL)
  {
    return (he->m_Edge->m_Reg[0]);
  }
  else
  {
    return (he->m_Edge->m_Reg[1]);
  }
}

template <typename TCoordRepType>
auto
VoronoiDiagram2DGenerator<TCoordRepType>::getLeftReg(FortuneHalfEdge * he) -> FortuneSite *
{
  if ((he->m_Edge) == nullptr)
  {
    return (m_BottomSite);
  }
  else if (he->m_RorL)
  {
    return (he->m_Edge->m_Reg[1]);
  }
  else
  {
    return (he->m_Edge->m_Reg[0]);
  }
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::insertEdgeList(FortuneHalfEdge * lbase, FortuneHalfEdge * lnew)
{
  lnew->m_Left = lbase;
  lnew->m_Right = lbase->m_Right;
  (lbase->m_Right)->m_Left = lnew;
  lbase->m_Right = lnew;
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::bisect(FortuneEdge * answer, FortuneSite * s1, FortuneSite * s2)
{
  answer->m_Reg[0] = s1;
  answer->m_Reg[1] = s2;
  answer->m_Ep[0] = nullptr;
  answer->m_Ep[1] = nullptr;

  double dx = (s2->m_Coord[0]) - (s1->m_Coord[0]);
  double dy = (s2->m_Coord[1]) - (s1->m_Coord[1]);
  double adx = (dx > 0) ? dx : -dx;
  double ady = (dy > 0) ? dy : -dy;

  answer->m_C = (s1->m_Coord[0]) * dx + (s1->m_Coord[1]) * dy + (dx * dx + dy * dy) * 0.5;
  if (adx > ady)
  {
    answer->m_A = 1.0;
    answer->m_B = dy / dx;
    answer->m_C /= dx;
  }
  else
  {
    answer->m_A = dx / dy;
    answer->m_B = 1.0;
    answer->m_C /= dy;
  }
  answer->m_Edgenbr = m_Nedges;
  ++m_Nedges;
  Point<int, 2> outline;
  outline[0] = answer->m_Reg[0]->m_Sitenbr;
  outline[1] = answer->m_Reg[1]->m_Sitenbr;
  m_OutputVD->AddLine(outline);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::intersect(FortuneSite * newV, FortuneHalfEdge * el1, FortuneHalfEdge * el2)
{
  FortuneEdge *     e1 = el1->m_Edge;
  FortuneEdge *     e2 = el2->m_Edge;
  FortuneHalfEdge * saveHE;
  FortuneEdge *     saveE;

  if (e1 == nullptr)
  {
    newV->m_Sitenbr = -1;
    return;
  }
  if (e2 == nullptr)
  {
    newV->m_Sitenbr = -2;
    return;
  }
  if ((e1->m_Reg[1]) == (e2->m_Reg[1]))
  {
    newV->m_Sitenbr = -3;
    return;
  }

  double d = (e1->m_A) * (e2->m_B) - (e1->m_B) * (e2->m_A);

  if ((d > -NUMERIC_TOLERENCE) && (d < NUMERIC_TOLERENCE))
  {
    newV->m_Sitenbr = -4;
    return;
  }

  double xmeet = ((e1->m_C) * (e2->m_B) - (e2->m_C) * (e1->m_B)) / d;
  double ymeet = ((e2->m_C) * (e1->m_A) - (e1->m_C) * (e2->m_A)) / d;

  if (comp(e1->m_Reg[1]->m_Coord, e2->m_Reg[1]->m_Coord))
  {
    saveHE = el1;
    saveE = e1;
  }
  else
  {
    saveHE = el2;
    saveE = e2;
  }

  bool right_of_site = (xmeet >= (saveE->m_Reg[1]->m_Coord[0]));
  if ((right_of_site && (!(saveHE->m_RorL))) || ((!right_of_site) && (saveHE->m_RorL)))
  {
    newV->m_Sitenbr = -4;
    return;
  }

  newV->m_Coord[0] = xmeet;
  newV->m_Coord[1] = ymeet;
  newV->m_Sitenbr = -5;
}

template <typename TCoordRepType>
auto
VoronoiDiagram2DGenerator<TCoordRepType>::getPQmin() -> FortuneHalfEdge *
{
  FortuneHalfEdge * curr = m_PQHash[m_PQmin].m_Next;

  m_PQHash[m_PQmin].m_Next = curr->m_Next;
  --m_PQcount;
  return (curr);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::clip_line(FortuneEdge * task)
{
  // Clip line
  FortuneSite * s1;
  FortuneSite * s2;
  double        x1, y1, x2, y2;

  if (((task->m_A) == 1.0) && ((task->m_B) >= 0.0))
  {
    s1 = task->m_Ep[1];
    s2 = task->m_Ep[0];
  }
  else
  {
    s2 = task->m_Ep[1];
    s1 = task->m_Ep[0];
  }

  int id1;
  int id2;
  if ((task->m_A) == 1.0)
  {
    if ((s1 != nullptr) && ((s1->m_Coord[1]) > m_Pymin))
    {
      y1 = s1->m_Coord[1];
      if (y1 > m_Pymax)
      {
        return;
      }
      x1 = s1->m_Coord[0];
      id1 = s1->m_Sitenbr;
    }
    else
    {
      y1 = m_Pymin;
      x1 = (task->m_C) - (task->m_B) * y1;
      id1 = -1;
    }

    if ((s2 != nullptr) && ((s2->m_Coord[1]) < m_Pymax))
    {
      y2 = s2->m_Coord[1];
      if (y2 < m_Pymin)
      {
        return;
      }
      x2 = s2->m_Coord[0];
      id2 = s2->m_Sitenbr;
    }
    else
    {
      y2 = m_Pymax;
      x2 = (task->m_C) - (task->m_B) * y2;
      id2 = -1;
    }

    if ((x1 > m_Pxmax) && (x2 > m_Pxmax))
    {
      return;
    }
    if ((x1 < m_Pxmin) && (x2 < m_Pxmin))
    {
      return;
    }
    if (x1 > m_Pxmax)
    {
      x1 = m_Pxmax;
      y1 = ((task->m_C) - x1) / (task->m_B);
      id1 = -1;
    }
    if (x1 < m_Pxmin)
    {
      x1 = m_Pxmin;
      y1 = ((task->m_C) - x1) / (task->m_B);
      id1 = -1;
    }
    if (x2 > m_Pxmax)
    {
      x2 = m_Pxmax;
      y2 = ((task->m_C) - x2) / (task->m_B);
      id2 = -1;
    }
    if (x2 < m_Pxmin)
    {
      x2 = m_Pxmin;
      y2 = ((task->m_C) - x2) / (task->m_B);
      id2 = -1;
    }
  }
  else
  {
    if ((s1 != nullptr) && ((s1->m_Coord[0]) > m_Pxmin))
    {
      x1 = s1->m_Coord[0];
      if (x1 > m_Pxmax)
      {
        return;
      }
      y1 = s1->m_Coord[1];
      id1 = s1->m_Sitenbr;
    }
    else
    {
      x1 = m_Pxmin;
      y1 = (task->m_C) - (task->m_A) * x1;
      id1 = -1;
    }
    if ((s2 != nullptr) && ((s2->m_Coord[0]) < m_Pxmax))
    {
      x2 = s2->m_Coord[0];
      if (x2 < m_Pxmin)
      {
        return;
      }
      y2 = s2->m_Coord[1];
      id2 = s2->m_Sitenbr;
    }
    else
    {
      x2 = m_Pxmax;
      y2 = (task->m_C) - (task->m_A) * x2;
      id2 = -1;
    }
    if ((y1 > m_Pymax) && (y2 > m_Pymax))
    {
      return;
    }
    if ((y1 < m_Pymin) && (y2 < m_Pymin))
    {
      return;
    }
    if (y1 > m_Pymax)
    {
      y1 = m_Pymax;
      x1 = ((task->m_C) - y1) / (task->m_A);
      id1 = -1;
    }
    if (y1 < m_Pymin)
    {
      y1 = m_Pymin;
      x1 = ((task->m_C) - y1) / (task->m_A);
      id1 = -1;
    }
    if (y2 > m_Pymax)
    {
      y2 = m_Pymax;
      x2 = ((task->m_C) - y2) / (task->m_A);
      id2 = -1;
    }
    if (y2 < m_Pymin)
    {
      y2 = m_Pymin;
      x2 = ((task->m_C) - y2) / (task->m_A);
      id2 = -1;
    }
  }

  VoronoiEdge newInfo;
  newInfo.m_Left[0] = x1;
  newInfo.m_Left[1] = y1;
  newInfo.m_Right[0] = x2;
  newInfo.m_Right[1] = y2;
  newInfo.m_LineID = task->m_Edgenbr;

  if (id1 > -1)
  {
    newInfo.m_LeftID = id1;
  }
  else
  {
    newInfo.m_LeftID = m_Nvert;
    ++m_Nvert;
    PointType newv;
    newv[0] = x1;
    newv[1] = y1;
    m_OutputVD->AddVert(newv);
  }

  if (id2 > -1)
  {
    newInfo.m_RightID = id2;
  }
  else
  {
    newInfo.m_RightID = m_Nvert;
    ++m_Nvert;
    PointType newv;
    newv[0] = x2;
    newv[1] = y2;
    m_OutputVD->AddVert(newv);
  }
  m_OutputVD->AddEdge(newInfo);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::makeEndPoint(FortuneEdge * task, bool lr, FortuneSite * ends)
{
  task->m_Ep[lr] = ends;
  if ((task->m_Ep[1 - lr]) == nullptr)
  {
    return;
  }

  clip_line(task);
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::GenerateVDFortune()
{
  unsigned int i;

  // Build seed sites
  m_SeedSites.resize(m_NumberOfSeeds);
  for (i = 0; i < m_NumberOfSeeds; ++i)
  {
    m_SeedSites[i].m_Coord = m_Seeds[i];
    m_SeedSites[i].m_Sitenbr = i;
  }
  // Initialize boundary
  m_Pxmax = m_VorBoundary[0];
  m_Pymax = m_VorBoundary[1];

  m_Deltay = m_Pymax - m_Pymin;
  m_Deltax = m_Pxmax - m_Pxmin;
  m_SqrtNSites = std::sqrt(static_cast<float>(m_NumberOfSeeds + 4));

  // Initialize outputLists
  m_Nedges = 0;
  m_Nvert = 0;
  m_OutputVD->LineListClear();
  m_OutputVD->EdgeListClear();
  m_OutputVD->VertexListClear();

  // Initialize the hash table for circle event and point event
  m_PQcount = 0;
  m_PQmin = 0;
  m_PQhashsize = static_cast<int>(4 * m_SqrtNSites);
  m_PQHash.resize(m_PQhashsize);
  for (i = 0; i < m_PQhashsize; ++i)
  {
    m_PQHash[i].m_Next = nullptr;
  }
  m_ELhashsize = static_cast<int>(2 * m_SqrtNSites);
  m_ELHash.resize(m_ELhashsize);
  for (i = 0; i < m_ELhashsize; ++i)
  {
    m_ELHash[i] = nullptr;
  }
  createHalfEdge(&(m_ELleftend), nullptr, 0);
  createHalfEdge(&(m_ELrightend), nullptr, 0);
  m_ELleftend.m_Left = nullptr;
  m_ELleftend.m_Right = &(m_ELrightend);
  m_ELrightend.m_Left = &(m_ELleftend);
  m_ELrightend.m_Right = nullptr;
  m_ELHash[0] = &(m_ELleftend);
  m_ELHash[m_ELhashsize - 1] = &(m_ELrightend);

  m_BottomSite = &(m_SeedSites[0]);
  FortuneSite * currentSite = &(m_SeedSites[1]);

  PointType currentCircle;
  currentCircle.Fill(0.0);
  FortuneHalfEdge * leftHalfEdge;
  FortuneHalfEdge * rightHalfEdge;
  FortuneHalfEdge * left2HalfEdge;
  FortuneHalfEdge * right2HalfEdge;
  FortuneHalfEdge * newHE;
  FortuneSite *     findSite;
  FortuneSite *     topSite;
  FortuneSite *     saveSite;
  FortuneEdge *     newEdge;
  FortuneSite *     meetSite;
  FortuneSite *     newVert;
  bool              saveBool;

  std::vector<FortuneHalfEdge> HEpool;
  std::vector<FortuneEdge>     Edgepool;
  std::vector<FortuneSite>     Sitepool;
  HEpool.resize(5 * m_NumberOfSeeds);
  Edgepool.resize(5 * m_NumberOfSeeds);
  Sitepool.resize(5 * m_NumberOfSeeds);

  int HEid = 0;
  int Edgeid = 0;
  int Siteid = 0;

  i = 2;
  bool ok = true;
  while (ok)
  {
    if (m_PQcount != 0)
    {
      PQshowMin(&currentCircle);
    }
    if ((i <= m_NumberOfSeeds) && ((m_PQcount == 0) || comp(currentSite->m_Coord, currentCircle)))
    {
      // Handling site event
      leftHalfEdge = findLeftHE(&(currentSite->m_Coord));
      rightHalfEdge = leftHalfEdge->m_Right;

      findSite = getRightReg(leftHalfEdge);

      bisect(&(Edgepool[Edgeid]), findSite, currentSite);
      newEdge = &(Edgepool[Edgeid]);
      ++Edgeid;

      createHalfEdge(&(HEpool[HEid]), newEdge, 0);
      newHE = &(HEpool[HEid]);
      ++HEid;

      insertEdgeList(leftHalfEdge, newHE);

      intersect(&(Sitepool[Siteid]), leftHalfEdge, newHE);
      meetSite = &(Sitepool[Siteid]);

      if ((meetSite->m_Sitenbr) == -5)
      {
        deletePQ(leftHalfEdge);
        insertPQ(leftHalfEdge, meetSite, dist(meetSite, currentSite));
        ++Siteid;
      }

      leftHalfEdge = newHE;
      createHalfEdge(&(HEpool[HEid]), newEdge, 1);
      newHE = &(HEpool[HEid]);
      ++HEid;

      insertEdgeList(leftHalfEdge, newHE);

      intersect(&(Sitepool[Siteid]), newHE, rightHalfEdge);
      meetSite = &(Sitepool[Siteid]);
      if ((meetSite->m_Sitenbr) == -5)
      {
        ++Siteid;
        insertPQ(newHE, meetSite, dist(meetSite, currentSite));
      }
      if (i < m_SeedSites.size())
      {
        currentSite = &(m_SeedSites[i]);
      }
      ++i;
    }
    else if (m_PQcount != 0)
    {
      // Handling circle event

      leftHalfEdge = getPQmin();
      left2HalfEdge = leftHalfEdge->m_Left;
      rightHalfEdge = leftHalfEdge->m_Right;
      right2HalfEdge = rightHalfEdge->m_Right;
      findSite = getLeftReg(leftHalfEdge);
      topSite = getRightReg(rightHalfEdge);

      newVert = leftHalfEdge->m_Vert;
      newVert->m_Sitenbr = m_Nvert;
      ++m_Nvert;
      m_OutputVD->AddVert(newVert->m_Coord);

      makeEndPoint(leftHalfEdge->m_Edge, leftHalfEdge->m_RorL, newVert);
      makeEndPoint(rightHalfEdge->m_Edge, rightHalfEdge->m_RorL, newVert);
      deleteEdgeList(leftHalfEdge);
      deletePQ(rightHalfEdge);
      deleteEdgeList(rightHalfEdge);

      saveBool = false;
      if ((findSite->m_Coord[1]) > (topSite->m_Coord[1]))
      {
        saveSite = findSite;
        findSite = topSite;
        topSite = saveSite;
        saveBool = true;
      }

      bisect(&(Edgepool[Edgeid]), findSite, topSite);
      newEdge = &(Edgepool[Edgeid]);
      ++Edgeid;

      createHalfEdge(&(HEpool[HEid]), newEdge, saveBool);
      newHE = &(HEpool[HEid]);
      ++HEid;

      insertEdgeList(left2HalfEdge, newHE);
      makeEndPoint(newEdge, 1 - saveBool, newVert);

      intersect(&(Sitepool[Siteid]), left2HalfEdge, newHE);
      meetSite = &(Sitepool[Siteid]);

      if ((meetSite->m_Sitenbr) == -5)
      {
        ++Siteid;
        deletePQ(left2HalfEdge);
        insertPQ(left2HalfEdge, meetSite, dist(meetSite, findSite));
      }

      intersect(&(Sitepool[Siteid]), newHE, right2HalfEdge);
      meetSite = &(Sitepool[Siteid]);
      if ((meetSite->m_Sitenbr) == -5)
      {
        ++Siteid;
        insertPQ(newHE, meetSite, dist(meetSite, findSite));
      }
    }
    else
    {
      ok = false;
    }
  }

  for ((leftHalfEdge = m_ELleftend.m_Right); (leftHalfEdge != (&m_ELrightend));
       (leftHalfEdge = (leftHalfEdge->m_Right)))
  {
    newEdge = leftHalfEdge->m_Edge;
    clip_line(newEdge);
  }
}

template <typename TCoordRepType>
void
VoronoiDiagram2DGenerator<TCoordRepType>::PrintSelf(std::ostream & os, Indent indent) const
{
  Superclass::PrintSelf(os, indent);

  os << indent << "Number Of Seeds: " << m_NumberOfSeeds << std::endl;
  os << indent << "VorBoundary: " << static_cast<typename NumericTraits<PointType>::PrintType>(m_VorBoundary)
     << std::endl;
  os << indent << "OutputVD: " << static_cast<typename NumericTraits<OutputType>::PrintType>(m_OutputVD) << std::endl;

  os << indent << "Pxmin: " << m_Pxmin << std::endl;
  os << indent << "Pxmax: " << m_Pxmax << std::endl;
  os << indent << "Pymin: " << m_Pymin << std::endl;
  os << indent << "Pymax: " << m_Pymax << std::endl;
  os << indent << "Deltax: " << m_Deltax << std::endl;
  os << indent << "Deltay: " << m_Deltay << std::endl;
  os << indent << "SqrtNSites: " << m_SqrtNSites << std::endl;

  os << indent << "PQcount: " << m_PQcount << std::endl;
  os << indent << "PQmin: " << m_PQmin << std::endl;
  os << indent << "PQhashsize: " << m_PQhashsize << std::endl;
  os << indent << "Nedges: " << m_Nedges << std::endl;
  os << indent << "Nvert: " << m_Nvert << std::endl;
  os << indent << "BottomSite: " << m_BottomSite << std::endl;

  os << indent << "ELhashsize: " << m_ELhashsize << std::endl;

  os << indent << "ELHash: " << std::endl;
  for (unsigned int i = 0; i < m_ELHash.size(); ++i)
  {
    os << indent << '[' << i << "]: " << m_ELHash[i] << std::endl;
  }
}
} // namespace itk

#endif