/*****************************************************************************/
/*                                                                           */
/*            Copyright (C)           1999-2025         M. Held              */
/*                                                                           */
/* This code is not in the public domain. All rights reserved! Please make   */
/* sure to read the full copyright statement contained in "README.txt" or in */
/* the "main" file of this code, such as "main.cc".                          */
/*                                                                           */
/*****************************************************************************/
/*                                                                           */
/* Written by:  Martin Held                                                  */
/*                                                                           */
/* E-Mail:      held@cs.sbg.ac.at                                            */
/* Fax Mail:    (+43 662) 8044-172                                           */
/* Voice Mail:  (+43 662) 8044-6304                                          */
/* Snail Mail:  Martin Held                                                  */
/*              FB Informatik                                                */
/*              Universitaet Salzburg                                        */
/*              A-5020 Salzburg, Austria                                     */
/*                                                                           */
/*****************************************************************************/

/*                                                                           */
/* get standard libraries                                                    */
/*                                                                           */
#include <math.h>
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <vector>
#include <set>
#include <map>
#include <queue>
#include <algorithm>

/*                                                                           */
/* get my header files                                                       */
/*                                                                           */
#include "fpkernel.h"
#include "vronivector.h"
#include "vroni_object.h"
#include "numerics.h"
#include "roots.h"
#include "intersections.h"


#ifdef GRAPHICS
#define HandleSplitSeg(I, J, K) \
{\
    K = SplitSeg(I, J); \
    if (graphics) { \
       AddToCurBuffer(I, SEG, AlertColor); \
       AddToCurBuffer(K, SEG, AlertColor); \
       AddToCurBuffer(J, PNT, AlertColor); \
    }\
}
#else
#define HandleSplitSeg(I, J, K) \
{\
    K = SplitSeg(I, J); \
}
#endif

#ifdef GRAPHICS
#define HandleSplitArc(I, J, K) \
{\
    K = SplitArc(I, J); \
    if (graphics) { \
       AddToCurBuffer(I, ARC, AlertColor); \
       AddToCurBuffer(K, ARC, AlertColor); \
       AddToCurBuffer(J, PNT, AlertColor); \
    }\
}
#else
#define HandleSplitArc(I, J, K) \
{\
    K = SplitArc(I, J); \
}
#endif


#ifdef DBG_GRAPHICS
void SketchGrid(grid* grid);
#endif



void vroniObject::ResetIntersectionStatus(void)
{
   check_completed = true;
   eps_int = ZERO;

   return;
}

   
vr_bool vroniObject::IsCheckComplete(void)
{
   return check_completed;
}


vr_bool vroniObject::ThresholdReached(void)
{
   return (eps_int > THRESHOLD);
}


vr_bool vroniObject::CheckIntersections(void)
{
   vr_bool int_found;

   int_found = ComputeAllIntersections(eps_int);

   if (check_completed  &&  (eps_int <= THRESHOLD)) {
      eps_int *= 10.0;
   }

   return  int_found;
}


double vroniObject::GetIntersectionThreshold(void)
{
   return eps_int;
}

vr_bool vroniObject::AddSegIntersection(seg_entry* seg, coord intersection) 
{
   return seg->intersections->insert(intersection).second;
}


#ifdef GENUINE_ARCS
vr_bool vroniObject::AddArcIntersection(arc_entry* arc, coord intersection) 
{
   arc_intersection ai;
   int i1;
   coord center, start;
   double alpha;

   i1 = arc->num;
   center = GetArcCenter(i1);
   start = GetArcStartCoord(i1);

   alpha = atan2(intersection.y - center.y, intersection.x - center.x) - 
           atan2(start.y - center.y, start.x - center.x);
   if (alpha < 0)
      alpha += M_2PI;

   ai.c = intersection;
   ai.alpha = alpha;

   return arc->intersections.insert(ai).second;
}
#endif


grid* vroniObject::InitialiseGrid(int num_elements, coord bb_min, coord bb_max) 
{
   grid* g;
   int i, j;
   int n_x, n_y;
   double d_x, d_y;
   double min_x, min_y;

   n_x = (int)(sqrt(machine_double(num_elements))) + 3;
   n_y = n_x;
   d_x = (2 * ZERO_MAX + bb_max.x - bb_min.x)/(n_x - 2);
   d_y = (2 * ZERO_MAX + bb_max.y - bb_min.y)/(n_y - 2);
   min_x = bb_min.x - d_x;
   min_y = bb_min.y - d_y;

   g = new grid;

   g->d_x = d_x;
   g->d_y = d_y;
   g->n_x = n_x;
   g->n_y = n_y;

   g->x_min = new double[n_x];
   g->x_max = new double[n_x];
   g->y_min = new double[n_y];
   g->y_max = new double[n_y];

   for (i = 0; i < n_x; i++) {
      g->x_min[i] = (i == 0 ? min_x : g->x_max[i-1]); //min_x + i*d_x
      g->x_max[i] = min_x + (i+1)*d_x;
   }

   for (j = 0; j < n_y; j++) {
      g->y_min[j] = (j == 0 ? min_y : g->y_max[j-1]); //min_y + j*d_y
      g->y_max[j] = min_y + (j+1)*d_y;
   }

   g->cells = new grid_cell*[n_x];
   for (i = 0; i < n_x; i++) {
      g->cells[i] = new grid_cell[n_y];
      for (j = 0; j < n_y; j++) {
         g->cells[i][j].x_min = g->x_min[i];
         g->cells[i][j].x_max = g->x_max[i];
         g->cells[i][j].y_min = g->y_min[j];
         g->cells[i][j].y_max = g->y_max[j];
         g->cells[i][j].i = i;
         g->cells[i][j].j = j;
      }
   }

   return g;
}

int vroniObject::FindCellX(grid* g, double_arg x) {
   return REAL_TO_INT((x - g->x_min[0])/g->d_x);
}

int vroniObject::FindCellY(grid* g, double_arg y) {
   return REAL_TO_INT((y - g->y_min[0])/g->d_y);
}


vr_bool vroniObject::CheckSegIntersectsGridVerSeg(int seg, double_arg x, 
                                                  double_arg y1, double_arg y2, 
                                                  double_arg eps) 
{
   double a, b, c;
   coord p, q;
   double bbx_min, bbx_max, bby_min, bby_max;

   assert (y1 < y2);

   GetSegEqnData(seg, &a, &b, &c);

   p    = GetSegStartCoord(seg);
   q    = GetSegEndCoord(seg);

   bbx_min = std::min(p.x, q.x);
   bbx_max = std::max(p.x, q.x);
   bby_min = std::min(p.y, q.y);
   bby_max = std::max(p.y, q.y);

   //check whether the bounding box of the segment contains at least part of the grid line segment
   if (lt (x - bbx_min, eps) || lt(bbx_max - x, eps)) return false;
   if (lt (y2 - bby_min, eps) || lt(bby_max - y1, eps)) return false;

   //check whether the segment crosses the grid line at the x-coordinate x at an y-coordinate between y1 and y2
   if (eq(b, eps)) {
      return eq(a*x + c, eps);
   } else if (gt(b, eps)) {
      return le(a*x + b*y1 + c, eps) && ge(a*x + b*y2 + c, eps);
   } else {
      return ge(a*x + b*y1 + c, eps) && le(a*x + b*y2 + c, eps);
   }
}

vr_bool vroniObject::CheckSegIntersectsGridHorSeg(int seg, 
                                                  double_arg x1, double_arg x2,
                                                  double_arg y, double_arg eps) 
{
   double a, b, c;
   coord p, q;
   double bbx_min, bbx_max, bby_min, bby_max;

   assert (x1 < x2);

   GetSegEqnData(seg, &a, &b, &c);

   p    = GetSegStartCoord(seg);
   q    = GetSegEndCoord(seg);

   bbx_min = std::min(p.x, q.x);
   bbx_max = std::max(p.x, q.x);
   bby_min = std::min(p.y, q.y);
   bby_max = std::max(p.y, q.y);

   //check whether the bounding box of the segment contains at least part of the grid line segment
   if (lt (x2 - bbx_min, eps) || lt(bbx_max - x1, eps)) return false;
   if (lt (y - bby_min, eps) || lt(bby_max - y, eps)) return false;

   //check whether the segment crosses the grid line at the y-coordinate y at an x-coordinate between x1 and x2
   if (eq(a, eps)) {
      return eq(b*y + c, eps);
   } else if (gt(a, eps)) {
      return le((a*x1 + b*y + c), eps) && ge(a*x2 + b*y + c, eps);
   } else {
      return ge((a*x1 + b*y + c), eps) && le(a*x2 + b*y + c, eps);
   }
}

#ifdef GENUINE_ARCS
vr_bool vroniObject::CheckArcIntersectsGridVerSeg(int arc, double_arg x, 
                                                  double_arg y1, double_arg y2, 
                                                  double_arg eps)
{
   double r;
   coord c;
   double xa, y1a, y2a, yy;
   double yymin, yymax;

   assert (y1 < y2);

   c = GetArcCenter(arc);

   if (y1 < c.y && y2 > c.y) { //simplify handling of this special case
      return CheckArcIntersectsGridVerSeg(arc, x, y1, c.y, eps)
            || CheckArcIntersectsGridVerSeg(arc, x, c.y, y2, eps);
   }

   xa = x - c.x;
   y1a = y1 - c.y;
   y2a = y2 - c.y;
   r = GetArcRadius(arc);

   yy = r*r - xa*xa;

   yymin = std::min(y1a*y1a, y2a*y2a);
   yymax = std::max(y1a*y1a, y2a*y2a);

   if (le(yymin - yy, eps) && ge(yymax - yy, eps)) {
      //grid segment intersects (or is very near to) the full circle - now check only the arc

      coord p, q;
      double arc_xmin, arc_xmax;

      p = GetArcStartCoord(arc);
      q = GetArcEndCoord(arc);

      //Note: arcs are always stored CCW

      if (y1 >= c.y) {
         //we only consider the upper half because y1, y2 >= c.y
         assert (y2 >= c.y);

         //Note: we don't allow arcs that span an angle of 180 or more degrees
         if (p.y < c.y) {
            //arc starts in the lower half of the circle
            if (q.y < c.y)
               return false; //arc is completely in the lower half
            arc_xmax = c.x + r;
         } else {
            arc_xmax = p.x;
         }

         if (q.y < c.y) {
            arc_xmin = c.x - r;
         } else {
            arc_xmin = q.x;
         }
      } else {
         //we only consider the lower half because y1, y2 <= c.y
         assert (y2 <= c.y);

         if (p.y > c.y) {
            if (q.y > c.y)
               return false; //arc is completely in the upper half
            arc_xmin = c.x - r;
         } else {
            //Note: we don't allow arcs that span an angle of 180 or more degrees
            arc_xmin = p.x;
         }

         if (q.y > c.y) {
            arc_xmax = c.x + r;
         } else {
            arc_xmax = q.x;
         }
      }

      return le(arc_xmin - x, eps) && ge(arc_xmax - x, eps);
   } else {
      return false;
   }
}

vr_bool vroniObject::CheckArcIntersectsGridHorSeg(int arc, 
                                                  double_arg x1, double_arg x2, 
                                                  double_arg y, double_arg eps) 
{
   double r;
   coord c;
   double ya, x1a, x2a, xx;
   double xxmin, xxmax;

   assert (x1 < x2);

   c = GetArcCenter(arc);

   if (x1 < c.x && x2 > c.x) { //simplify handling of this special case
      return CheckArcIntersectsGridHorSeg(arc, x1, c.x, y, eps)
            || CheckArcIntersectsGridHorSeg(arc, c.x, x2, y, eps);
   }

   ya = y - c.y;
   x1a = x1 - c.x;
   x2a = x2 - c.x;
   r = GetArcRadius(arc);

   xx = r*r - ya*ya;

   xxmin = std::min(x1a*x1a, x2a*x2a);
   xxmax = std::max(x1a*x1a, x2a*x2a);

   if (le(xxmin - xx, eps) && ge(xxmax - xx, eps)) {
      //grid segment intersects (or is very near to) the full circle - now check only the arc

      coord p, q;
      double arc_ymin, arc_ymax;

      p = GetArcStartCoord(arc);
      q = GetArcEndCoord(arc);

      //Note: arcs are always stored CCW

      if (x1 >= c.x) {
         //we only consider the right half because x1, x2 >= c.x
         assert (x2 >= c.x);

         //Note: we don't allow arcs that span an angle of 180 or more degrees
         if (p.x < c.x) {
            //arc starts in the left half of the circle
            if (q.x < c.x)
               return false; //arc is completely in the left half
            arc_ymin = c.y - r;
         } else {
            arc_ymin = p.y;
         }

         if (q.x < c.x) {
            arc_ymax = c.y + r;
         } else {
            arc_ymax = q.y;
         }
      } else {
         //we only consider the left half because x1, x2 <= c.x
         assert (x2 <= c.x);

         if (p.x > c.x) {
            //arc starts in the right half of the circle
            if (q.x > c.x)
               return false; //arc is completely in the right half
            arc_ymax = c.y + r;
         } else {
            //Note: we don't allow arcs that span an angle of 180 or more degrees
            arc_ymax = p.y;
         }

         if (q.x > c.x) {
            arc_ymin = c.y - r;
         } else {
            arc_ymin = q.y;
         }
      }

      return le(arc_ymin - y, eps) && ge(arc_ymax - y, eps);
   } else {
      return false;
   }
}
#endif


void vroniObject::AddPointToGrid(grid* g, int i1, double_arg eps) 
{
   int i, j;
   point_entry* entry;
   coord p;

   assert(InPntsList(i1));

   p   = GetPntCoords(i1);
   entry = new point_entry;
   entry->num = i1;

   for (i = FindCellX(g, p.x - eps); i <= FindCellX(g, p.x + eps); i++) {
      for (j = FindCellY(g, p.y - eps); j <= FindCellY(g, p.y + eps); j++) {
         assert (!(i < 0 || j < 0 || i >= g->n_x || j >= g->n_y));

         entry->cells.insert(&(g->cells[i][j]));
         g->cells[i][j].points.insert(entry);
      }
   }

   g->points[i1] = entry;
}

void vroniObject::AddSegToGrid(grid* g, int i1, double_arg eps) 
{
   int i, j;
   seg_entry* entry;
   coord p, q;
   std::queue<std::pair<int,int> > queue;
   double x, y;

   assert(InSegsList(i1));

   entry = new seg_entry(comp_eps);
   entry->num = i1;

   p    = GetSegStartCoord(i1);
   q    = GetSegEndCoord(i1);

   /*
    * we basically do a BFS from the seg's start vertex;
    * two cells are connected if and only if the seg intersects the grid line segment between them
    */

   //First find the start cell(s) -> can be up to four cells (if the start vertex is at a corner of a grid cell)
   for (i = FindCellX(g, p.x - eps); i <= FindCellX(g, p.x + eps); i++) {
      for (j = FindCellY(g, p.y - eps); j <= FindCellY(g, p.y + eps); j++) {
         assert (!(i < 0 || j < 0 || i >= g->n_x || j >= g->n_y));
         queue.push(std::make_pair(i, j));
      }
   }

   while (!queue.empty()) {
      std::pair<int,int> cell;

      cell = queue.front(); queue.pop();
      i = cell.first;
      j = cell.second;

      if (!(i < 0 || j < 0 || i >= g->n_x || j >= g->n_y)) {

         if (!g->cells[i][j].segs.insert(entry).second)
            continue; //we have already visited this cell
         
         x = (q.x >= p.x ? g->x_max[i] : g->x_min[i]);
         y = (q.y >= p.y ? g->y_max[j] : g->y_min[j]);
         
         if (CheckSegIntersectsGridVerSeg(i1, x, g->y_min[j], g->y_max[j], eps)) {
            queue.push(std::make_pair(q.x >= p.x ? i+1 : i-1, j));
         }
         if (CheckSegIntersectsGridHorSeg(i1, g->x_min[i], g->x_max[i], y, eps)) {
            queue.push(std::make_pair(i, q.y >= p.y ? j+1 : j-1));
         }
      }
   }

   g->segs[i1] = entry;
}


#ifdef GENUINE_ARCS
void vroniObject::AddArcToGrid(grid* g, int i1, double_arg eps) 
{
   int i, j;
   arc_entry* entry;
   coord p;
   std::queue<std::pair<int,int> > queue;

   assert(InArcsList(i1));

   entry = new arc_entry;
   entry->num = i1;

   p = GetArcStartCoord(i1);

   /*
    * we basically do a BFS from the arc's start vertex;
    * two cells are connected if and only if the arc intersects the grid line segment between them
    */

   //First find the start cell(s) -> can be up to four cells (if the start vertex is at a corner of a grid cell)
   for (i = FindCellX(g, p.x - eps); i <= FindCellX(g, p.x + eps); i++) {
      for (j = FindCellY(g, p.y - eps); j <= FindCellY(g, p.y + eps); j++) {
         assert (!(i < 0 || j < 0 || i >= g->n_x || j >= g->n_y));

         queue.push(std::make_pair(i, j));
      }
   }

   while (!queue.empty()) {
      std::pair<int,int> cell;

      cell = queue.front(); queue.pop();
      i = cell.first;
      j = cell.second;

      if (!(i < 0 || j < 0 || i >= g->n_x || j >= g->n_y)) {

         if (!g->cells[i][j].arcs.insert(entry).second)
            continue; //we have already visited this cell
         
         if (CheckArcIntersectsGridVerSeg(i1, g->x_min[i], g->y_min[j], g->y_max[j], eps))
            queue.push(std::make_pair(i-1, j));
         if (CheckArcIntersectsGridHorSeg(i1, g->x_min[i], g->x_max[i], g->y_min[j], eps))
            queue.push(std::make_pair(i, j-1));
         if (CheckArcIntersectsGridVerSeg(i1, g->x_max[i], g->y_min[j], g->y_max[j], eps))
            queue.push(std::make_pair(i+1, j));
         if (CheckArcIntersectsGridHorSeg(i1, g->x_min[i], g->x_max[i], g->y_max[j], eps))
            queue.push(std::make_pair(i, j+1));
      }
   }
    
   g->arcs[i1] = entry;
}
#endif


void vroniObject::FillGrid(grid* g, double_arg eps)
{
   int i;

   for (i = 2; i <= num_pnts - 3; i++) {
      if (!IsDeletedPnt(i)) {
         AddPointToGrid(g, i, eps);
      }
   }

   for (i = 0; i <= num_segs - 1; i++)
      AddSegToGrid(g, i, eps);

#ifdef GENUINE_ARCS
   for (i = 0; i <= num_arcs - 1; i++)
      AddArcToGrid(g, i, eps);
#endif
}


void vroniObject::DeleteGrid(grid* grid) 
{
   std::map<int,point_entry*>::iterator it1;
   std::map<int,seg_entry*>::iterator it2;
#ifdef GENUINE_ARCS
   std::map<int,arc_entry*>::iterator it3;
#endif
   int i;

   for (it1 = grid->points.begin(); it1 != grid->points.end(); it1++)
      delete (*it1).second;

   for (it2 = grid->segs.begin(); it2 != grid->segs.end(); it2++)
      delete (*it2).second;

#ifdef GENUINE_ARCS
   for (it3 = grid->arcs.begin(); it3 != grid->arcs.end(); it3++)
      delete (*it3).second;
#endif

   delete [] grid->x_min;
   delete [] grid->x_max;
   delete [] grid->y_min;
   delete [] grid->y_max;

   for (i = 0; i < grid->n_x; i++)
      delete [] grid->cells[i];

   delete [] grid->cells;

   delete grid;
}


#ifdef DBG_GRAPHICS

static void SketchCell(grid* grid, int i, int j)
{
   coord p, q;
   grid_cell* c = &grid->cells[i][j];

   int color = !c->segs.empty()
#ifdef GENUINE_ARCS
         || !c->arcs.empty()
#endif
         || !c->points.empty() ? GridColor : NoColor;

   p.x = c->x_min;
   q.x = c->x_max;
   p.y = c->y_min;
   q.y = c->y_min;
   AddEdgeToBuffer(p.x, p.y, q.x, q.y, color);
   p.y = c->y_max;
   q.y = c->y_max;
   AddEdgeToBuffer(p.x, p.y, q.x, q.y, color);
   p.x = c->x_min;
   q.x = c->x_min;
   p.y = c->y_min;
   q.y = c->y_max;
   AddEdgeToBuffer(p.x, p.y, q.x, q.y, color);
   p.x = c->x_max;
   q.x = c->x_max;
   AddEdgeToBuffer(p.x, p.y, q.x, q.y, color);
}


//only for debugging purposes
static void SketchGrid(grid* grid)
{
   int i, j;
   for (i = 0;  i < grid->n_x;  ++i)
      for (j = 0;  j < grid->n_y;  ++j)
         SketchCell(grid, i, j);
}
#endif


void vroniObject::SplitAtIntersections(grid* grid, double_arg eps) 
{
   int i1, i2, i3;
   std::set<coord,comp_seg_iscoord>::iterator it;
   std::map<int,seg_entry*>::iterator it2;
   vr_bool reverse;
#ifdef GENUINE_ARCS
   std::set<arc_intersection,comp_arc_iscoord>::iterator it1;
   std::map<int,arc_entry*>::iterator it3;
#endif

   comp_eps = eps;

   for (it2 = grid->segs.begin(); it2 != grid->segs.end(); it2++) {
      seg_entry* sege = (*it2).second;
      i1 = sege->num;
      coord s, e;

      s = GetSegStartCoord(i1);
      e = GetSegEndCoord(i1);

      reverse = lt(e.x - s.x, comp_eps) || (eq(e.x - s.x, comp_eps) && e.y < s.y);

      for (it = sege->intersections->begin(); it != sege->intersections->end(); it++) {
         coord c = *it;
#ifdef EXT_APPL_PNTS
        // MODIF: al nuovo punto aggiunto dall'intersezione posso assegnare soltanto il loop originario, visto che non deriva
        // da alcun punto della curva di partenza
         i3 = StorePnt(c.x, c.y, {segs[i1].ext_appl.first, -1});
#else
         i3 = StorePnt(c.x, c.y);
#endif

         HandleSplitSeg(i1, i3, i2);
         if (!reverse)
            i1 = i2;
      }
   }

#ifdef GENUINE_ARCS
   for (it3 = grid->arcs.begin(); it3 != grid->arcs.end(); it3++) {
      arc_entry* arce = (*it3).second;
      i1 = arce->num;

      for (it1 = arce->intersections.begin(); it1 != arce->intersections.end(); it1++) {
         coord c = (*it1).c;
#ifdef EXT_APPL_PNTS
        // MODIF: al nuovo punto aggiunto dall'intersezione posso assegnare soltanto il loop originario, visto che non deriva
        // da alcun punto della curva di partenza
         i3 = StorePnt(c.x, c.y, {arcs[i1].ext_appl.first, -1});
#else
         i3 = StorePnt(c.x, c.y);
#endif

         HandleSplitArc(i1, i3, i1);
      }
   }
#endif
}


vr_bool vroniObject::CheckPntPnt(int i1, int i2, double_arg eps)
{
   int i3;
   double dist;
   coord p1, p2;

   assert(InPntsList(i1));
   assert(InPntsList(i2));

   if (i1 == i2)  return  false;

   if (i2 < i1)  VroniSwap(i1, i2, i3);

   p1   = GetPntCoords(i1);
   p2   = GetPntCoords(i2);
   dist = PntPntDist(p1, p2);

   if (eq(dist, eps)) {
      IntWarning("pnt", i1, "pnt", i2);
      return  true;
   }

   return  false;
}



vr_bool vroniObject::CheckPntOnSeg(int seg, int pnt, double_arg eps)
{
   double a, b, c, lgth, dist_q;
   coord  p, q;

   if (!IsSegStartPnt(seg, pnt) && !IsSegEndPnt(seg, pnt)) {
      GetSegEqnData(seg, &a, &b, &c);
      q    = GetPntCoords(pnt);
      p    = GetSegStartCoord(seg);
      lgth = GetSegLgth(seg);
      dist_q = PntSegDist(p, q, a, b, c, lgth, TINY);

      if (eq(dist_q, eps)) {
         IntWarning("seg", seg, "pnt", pnt);
         return true;
      }
   }

   return  false;
}



#ifdef GENUINE_ARCS
vr_bool vroniObject::CheckPntOnArc(int arc, int pnt, double_arg eps)
{
   double r, dist;
   coord  p3, p, cp, ns, ne;

   if (!IsArcStartPnt(arc, pnt) && !IsArcEndPnt(arc, pnt)) {
      p3 = GetArcCenter(arc);
      r  = GetArcRadius(arc);
      p  = GetPntCoords(pnt);
      ns = GetArcStartNormal(arc);
      ne = GetArcEndNormal(arc);
      cp = VecSub(p, p3);
      dist = PntArcDist(ns, ne, cp, r, TINY);

      if (eq(dist, eps)) {
         IntWarning("arc", arc, "pnt", pnt);
         return  true;
      }
   }

   return  false;
}
#endif




/*                                                                           */
/* check (and compute) the intersection of seg  i1  and seg  i2,  and split  */
/* the sites accordingly.                                                    */
/*                                                                           */
/* note: (1) emphasis is placed on reliability rather than speed!            */
/*       (2) an iterative solver is used if the analytical computation of    */
/*           the points of intersection does not yield accurate results.     */
/*       (3) it is assumed that a pnt-on-seg test has been carried out prior */
/*           to calling this routine.                                        */
/*                                                                           */
vr_bool vroniObject::GetSegSegIntersection(int i1, int i2, coord* result, 
                                           double_arg eps, vr_bool *ident_segs)
{
   int s_p, s_q, s_u, s_v, i3, i5;
   double a1, b1, c1, lgth, dist_p, dist_q, dist_u, dist_v;
   double a2, b2, c2, t;
   coord dir, p, q, u, v, w;
   int ident = 0;
#ifndef NDEBUG 
   vr_bool success;
#endif

   if (i1 == i2)  return  false;
   *ident_segs = false;

   if (i2 < i1) VroniSwap(i1, i2, i3);

   /*                                                                        */
   /* compute the sidedness of the end points  u,v  of seg 2 w.r.t. the      */
   /* supporting line of seg 1. note that the segs cannot overlap as we      */
   /* assume that a pnt-on-seg test has been carried out prior to calling    */
   /* this function. thus, we may stop the test once one vertex lies on the  */
   /* supporting line of the other seg.                                      */
   /*                                                                        */
   GetSegEqnData(i1, &a1, &b1, &c1);
   p    = GetSegStartCoord(i1);
   q    = GetSegEndCoord(i1);
   lgth = GetSegLgth(i1);
   u    = GetSegStartCoord(i2);
   v    = GetSegEndCoord(i2);

   i5 = Get1stVtx(i2, SEG);
   if (IsSegStartPnt(i1, i5) || IsSegEndPnt(i1, i5)) {
      s_u    = 0;
      ident  = 1;
   }
   else {
      dist_u = PntLineDist(a1, b1, c1, u);
      s_u    = SignEps(dist_u, eps);
   }

   i5 = Get2ndVtx(i2, SEG);
   if (IsSegStartPnt(i1, i5) || IsSegEndPnt(i1, i5)) {
      s_v    = 0;
      ++ident;
   }
   else {
      dist_v = PntLineDist(a1, b1, c1, v);
      s_v    = SignEps(dist_v, eps);
   }


   if (ident == 2) {
      //we don't allow duplicate segs
      *ident_segs = true;
      return false;
   }

   if (s_u * s_v == -1) {
      /*                                                                     */
      /* compute the sidedness of the end points  p,q  of seg 1 w.r.t. the   */
      /* supporting line of seg 2.                                           */
      /*                                                                     */
      GetSegEqnData(i2, &a2, &b2, &c2);
      dist_p = PntLineDist(a2, b2, c2, p);
      s_p    = SignEps(dist_p, eps);
      dist_q = PntLineDist(a2, b2, c2, q);
      s_q    = SignEps(dist_q, eps);

      if (s_p * s_q == -1) {
         /*                                                                  */
         /* this is a genuine intersection. let's compute it as accurately   */
         /* as possible.                                                     */
         /*                                                                  */
         dir = GetSegDirection(i1);
         t = lgth * (Abs(dist_p) / (Abs(dist_p) + Abs(dist_q)));
         w = RayPnt(p, dir, t);

#ifndef NDEBUG 
         success = IterativeIntersection(SEG, SEG, a2, b2, c2, p, 0.0, p,
                                         0.0, lgth, dir, 0.0, &t, &w);
         assert(success);
#else
         (void) IterativeIntersection(SEG, SEG, a2, b2, c2, p, 0.0, p,
                                      0.0, lgth, dir, 0.0, &t, &w);
#endif

         IntWarning("seg", i1, "seg", i2);

         *result = w;
         return  true;
      }
   }

   return  false;
}


#ifdef GENUINE_ARCS
/*                                                                           */
/* check (and compute) the intersection of seg  i1   and arc  i2,  and split */
/* the sites accordingly.                                                    */
/*                                                                           */
/* note: (1) emphasis is placed on reliability rather than speed!            */
/*       (2) an iterative solver is used if the analytical computation of    */
/*           the points of intersection does not yield accurate results.     */
/*       (3) it is assumed that a pnt-on-seg and pnt-on-arc test has been    */
/*           carried out.                                                    */
/*                                                                           */
int vroniObject::GetArcArcIntersection(int i1, int i2, 
                                       coord* res1, coord* res2,
                                       double_arg eps, vr_bool *ident_arcs)
{
   double R1, R2, d, alpha, alpha1, alpha2, t1, t2, f;
   coord q1, q2, q3, q4, C1, C2, v, Ns1, Ne1, Ns2, Ne2, shared_point;
   vr_bool success, shared = false;
   int i3;
   int tmp, intersections = 0;
   double lo, hi, a1, a2, d1, d2;
   int i, MAX_ITERATIONS = 300;


   if (i1 == i2) return  0;
   *ident_arcs = false;
   

   R1 = GetArcRadius(i1);
   assert(gt(R1, eps));
   R2 = GetArcRadius(i2);
   assert(gt(R2, eps));

   if ((R1 < R2)  ||  ((R1 == R2)  &&  (i2 < i1)))  {
      VroniSwap(R1, R2, d);
      VroniSwap(i1, i2, tmp);
   }

   /*                                                                        */
   /* first check whether the arcs share one or two endpoints                */
   /*                                                                        */
   i3 = Get1stVtxArc(i1);
   shared_point = GetPntCoords(i3);
   if (IsArcStartPnt(i2, i3) || IsArcEndPnt(i2, i3)) {
      shared = true;
   }

   i3 = Get2ndVtxArc(i1);
   if (IsArcStartPnt(i2, i3) || IsArcEndPnt(i2, i3)) {
      if (shared) {
         /*                                                                 */
         /* the arcs share both end points                                  */
         /*                                                                 */
         C1 = GetArcCenter(i1);
         C2 = GetArcCenter(i2);
         d  = PntPntDist(C1, C2);
         if (eq(d, eps)) {
            /*                                                               */
            /* the arcs are either identical or complementary. since we do   */
            /* not allow arcs that span an angle of 180 or more degrees, the */
            /* arcs have to be identical.                                    */
            /*                                                               */
            *ident_arcs = true;
            return 0;
         }
         else {
            /*                                                               */
            /* there can't be a non-trivial intersection                     */
            /*                                                               */
            return 0;
         }
      }
      shared = true;
      shared_point = GetPntCoords(i3);
   }

   C1 = GetArcCenter(i1);
   C2 = GetArcCenter(i2);
   v = VecSub(C2, C1);
   d = VecLen(v);

   Ns1 = GetArcStartNormal(i1);
   Ne1 = GetArcEndNormal(i1);
   Ns2 = GetArcStartNormal(i2);
   Ne2 = GetArcEndNormal(i2);

   if (eq(d, eps)) {
      /*                                                                     */
      /* the two arcs have the same center. a partial overlap of the arcs    */
      /* cannot occur as we assume that a pnt-on-arc test has been carried   */
      /* out                                                                 */
      /*                                                                     */
      return  0;
   }

   if (lt(d - R1 + R2, eps)   ||  gt(d - R1 - R2 , eps)) { 
      /*                                                                     */
      /* circle 2 is completely within circle 1   or                         */
      /* the circles don't overlap                                           */
      /*                                                                     */
      return 0;
   }

   alpha = atan2(v.y, v.x);
   alpha1 = alpha - M_PI;
   alpha2 = alpha + M_PI;
   t1 = alpha - M_PI_2;
   t2 = alpha + M_PI_2;

   if (eq(d - R1 - R2, eps)) { 
      /*                                                                     */
      /* the circles touch at one point (from the outside)                   */
      /*                                                                     */
      if ((d - R1 - R2) >= 0.0) {
         /*                                                                  */
         /* in this case we don't really need an iterative search, as we     */
         /* know the circles intersect exactly once and the analytic result  */
         /* is accurate enough                                               */
         /*                                                                  */
         f = (R1 - R2 + d)/2/d;
         q1 = VecAdd(C1, VecMult(f, v));

         q3 = VecSub(q1, C1);
         q4 = VecSub(q1, C2);

         if ((IsInArcConeStrict(Ns1, Ne1, q3) == 0)  && 
             (IsInArcConeStrict(Ns2, Ne2, q4) == 0)  && 
             (!shared || ne(PntPntDist(q1, shared_point), TINY))) {
            *res1 = q1;
            return 1;
         } 
         else {
            return 0;
         }
      }

      /*                                                                     */
      /* calculate the point on circle 2 that should be nearest to circle 1  */
      /*                                                                     */
      q1 = CirclePnt(C2, R2, alpha1);
      d1 = PntCircleDist(C1, R1, q1);

      if (d1 >= 0) {
         /*                                                                  */
         /* The calculated nearest point is not inside circle 1 (this might  */
         /* be due to numerical inaccuracies or because no such point        */
         /* exists): thus we can't do a iterative search using bisection...  */
         /* Use a ternary search to improve the calculation of the farthest  */
         /* point                                                            */
         /*                                                                  */
         lo = alpha1 - M_PI_2;
         hi = alpha1 + M_PI_2;

         for (i = 0; i < MAX_ITERATIONS; i++) {
            a1 = lo + (hi - lo)/3;
            a2 = lo + (hi - lo)*2/3;
            if (a1 == lo || a2 == hi) 
               /* we have reached machine precision, so we can stop          */
               break;

            q1 = CirclePnt(C2, R2, a1);
            q2 = CirclePnt(C2, R2, a2);

            d1 = PntCircleDist(C1, R1, q1);
            d2 = PntCircleDist(C1, R1, q2);

            if (d1 < d2)
               hi = a2;
            else
               lo = a1;
         }

         if (d1 >= 0) { 
            /*                                                               */
            /* nearest point on circle 2 is on or outside circle 1 (but very */
            /* near), so we return it as an intersection                     */
            /*                                                               */
            assert (eq(d1, 2*eps));

            q3 = VecSub(q1, C1);
            q4 = VecSub(q1, C2);

            if ((IsInArcConeStrict(Ns1, Ne1, q3) == 0)  && 
                (IsInArcConeStrict(Ns2, Ne2, q4) == 0)  && 
                (!shared || ne(PntPntDist(q1, shared_point), TINY))) {
               *res1 = q1;
               return 1;
            } 
            else {
               return 0;
            }
         } 
         else { 
            /*                                                               */
            /* we found a point on circle 2 that lies inside circle 1, so we */
            /* can do the usual iterative search (bisection)                 */
            /*                                                               */
            alpha = a1;
            alpha1 = alpha - M_PI;
            alpha2 = alpha + M_PI;
            t1 = a1;
            t2 = a1;
         }
      }
   } 
   else if (eq(R2 - R1 + d, eps)) { 
      /*                                                                     */
      /* circle 2's outside touches circle 1 from the inside                 */
      /*                                                                     */
      /* calculate the point on circle 2 that should be farthest from the    */
      /* center of circle 1                                                  */
      /*                                                                     */
      q1 = CirclePnt(C2, R2, alpha);
      d1 = PntCircleDist(C1, R1, q1);

      if (d1 <= 0) {
         /*                                                                  */
         /* the calculated farthest point is not outside circle 1 (this      */
         /* might be due to numerical inaccuracies or because no such point  */
         /* exists): thus, we can't do a iterative search using bisection... */
         /* Use a ternary search to improve the calculation of the farthest  */
         /* point                                                            */
         /*                                                                  */
         lo = alpha - M_PI_2;
         hi = alpha + M_PI_2;

         for (i = 0; i < MAX_ITERATIONS; i++) {
            a1 = lo + (hi - lo)/3;
            a2 = lo + (hi - lo)*2/3;
            if (a1 == lo || a2 == hi) 
               /* we have reached machine precision, so we can stop          */
               break;

            q1 = CirclePnt(C2, R2, a1);
            q2 = CirclePnt(C2, R2, a2);

            d1 = PntCircleDist(C1, R1, q1);
            d2 = PntCircleDist(C1, R1, q2);

            if (d1 > d2)
               hi = a2;
            else
               lo = a1;
         }

         if (d1 <= 0) { 
            /*                                                               */
            /* farthest point on circle 2 is on or inside circle 1 (but      */
            /* very near), so we return it as an intersection                */
            /*                                                               */
            assert (eq(d1, 2*eps));

            q3 = VecSub(q1, C1);
            q4 = VecSub(q1, C2);

            if ((IsInArcConeStrict(Ns1, Ne1, q3) == 0)  && 
                (IsInArcConeStrict(Ns2, Ne2, q4) == 0)  && 
                (!shared || ne(PntPntDist(q1, shared_point), TINY))) {
               *res1 = q1;
               return 1;
            } 
            else {
               return 0;
            }
         } 
         else { 
            /*                                                               */
            /* we found a point on circle 2 that lies outside circle 1, so   */
            /* we can do the usual iterative search (bisection)              */
            /*                                                               */
            alpha = a1;
            alpha1 = alpha - M_PI;
            alpha2 = alpha + M_PI;
            t1 = a1;
            t2 = a1;
         }
      }
   }

   success = IterativeIntersection(ARC, ARC, 0.0, 0.0, 0.0, C1, R1, C2,
         alpha1, alpha, v, R2, &t1, &q1);

   q3 = VecSub(q1, C1);
   q4 = VecSub(q1, C2);
   if ((success && IsInArcConeStrict(Ns1, Ne1, q3) == 0)  && 
       (IsInArcConeStrict(Ns2, Ne2, q4) == 0)  && 
       (!shared || ne(PntPntDist(q1, shared_point), TINY))) {
      *res1 = q1;
      intersections++;
   }

   success = IterativeIntersection(ARC, ARC, 0.0, 0.0, 0.0, C1, R1, C2,
                                   alpha, alpha2, v, R2, &t2, &q2);

   q3 = VecSub(q2, C1);
   q4 = VecSub(q2, C2);
   if (success  && 
       (IsInArcConeStrict(Ns1, Ne1, q3) == 0)  && 
       (IsInArcConeStrict(Ns2, Ne2, q4) == 0)  && 
       (!shared || ne(PntPntDist(q2, shared_point), TINY))) {
      if (intersections == 1) {
         if (ne(PntPntDist(q1, q2), eps)) {
            *res2 = q2;
            intersections++;
         }
      } 
      else {
         *res1 = q2;
         intersections++;
      }
   }

   if (intersections > 0)
      IntWarning("arc", i1, "arc", i2);

   return intersections;
}


/*                                                                           */
/* check (and compute) the intersection of seg  i1   and arc  i2,  and split */
/* the sites accordingly.                                                    */
/*                                                                           */
/* note: (1) emphasis is placed on reliability rather than speed!            */
/*       (2) an iterative solver is used if the analytical computation of    */
/*           the points of intersection does not yield accurate results.     */
/*       (3) it is assumed that a pnt-on-seg and pnt-on-arc test has been    */
/*           carried out.                                                    */
/*                                                                           */
int vroniObject::GetSegArcIntersection(int seg, int arc, coord* res1, 
                                       coord* res2, double_arg eps)
{
   int j1, j2;
   double lgth, d, R, u1, u2, t, dist, a, b, c;
   coord C, S, E, v1, v2, Ns, Ne, q1, q2, q3;
   vr_bool shared, success;
   coord shared_point;
   int intersections;

   /*                                                                        */
   /* first check whether the seg and the arc share one or two end points    */
   /*                                                                        */
   intersections = 0;
   shared = false;
   j1 = Get1stVtx(seg, SEG);
   shared_point = GetSegStartCoord(seg);
   if (IsArcStartPnt(arc, j1)  ||  IsArcEndPnt(arc, j1)) {
      shared = true;
   }

   j2 = Get2ndVtx(seg, SEG);
   if (IsArcStartPnt(arc, j2)  ||  IsArcEndPnt(arc, j2)) {
      if (shared) {
         /*                                                                  */
         /* the arc and the seg share both end points. thus, only trivial    */
         /* intersections exist.                                             */
         /*                                                                  */
         return 0;
      }
      shared = true;
      shared_point = GetSegEndCoord(seg);
   }

   /*                                                                        */
   /* compute the parameter  d  of the projection  p1 + v1 * d  of the arc's */
   /* center on the infinite line.                                           */
   /* d is important to determine the interval for the iterative computation */
   /*                                                                        */
   C     = GetArcCenter(arc);
   S     = GetSegStartCoord(seg);
   E     = GetSegEndCoord(seg);
   lgth  = GetSegLgth(seg);
   v1    = GetSegDirection(seg);
   Ns    = GetArcStartNormal(arc);
   Ne    = GetArcEndNormal(arc);

   assert(eq(VecLen(v1) - 1.0, eps));

   v2 = VecSub(C, S);
   d  = VecDotProd(v1, v2);

   GetSegEqnData(seg, &a, &b, &c);
   R  = GetArcRadius(arc);
   dist = AbsPntLineDist(a, b, c, C);

   if (lt(R - dist, eps)) { //line is fully outside of circle
      return 0;
   }

   /*                                                                        */
   /* check whether the line is tangential                                   */
   /*                                                                        */
   /* calculate the point on the segment that is nearest to the circle's     */
   /* center                                                                 */
   /*                                                                        */
   q1 = VecAdd(S, VecMult(d, v1));
   if (PntCircleDist(C, R, q1) >= 0) {
      /*                                                                     */
      /* the nearest point lies on the circle (or slightly outside)          */
      /* -> in this case the iterative solver won't work, because we don't   */
      /*    know a point on the segment that is inside the circle            */
      /* -> assume the line is tangential and use the nearest point as       */
      /*    a candidate for the intersection                                 */
      /*                                                                     */
	   q3 = VecSub(q1, C);
	   if ((IsInArcConeStrict(Ns, Ne, q3) == 0)  && 
           (IsInSegConeStrict(S, q1, a, b, lgth) == 0)  && 
           (!shared || ne(PntPntDist(q1, shared_point), TINY))) {
	      *res1 = q1;
          //          IntWarning("seg", seg, "arc", arc);
	      return 1;
	   } 
       else {
	      return 0;
	   }
   } 

   /*                                                                        */
   /* else we can use the iterative solver                                   */
   /*                                                                        */
   u1 = d - R - eps;
   u2 = d + R + eps;

   t = (u1 + d)/2;
   success = IterativeIntersection(ARC, SEG, 0.0, 0.0, 0.0,
                                            C, R, S, u1, d, v1, 0.0,
                                            &t, &q1);

   q3 = VecSub(q1, C);

   if (success  &&  
       (IsInArcConeStrict(Ns, Ne, q3) == 0)  && 
       (IsInSegConeStrict(S, q1, a, b, lgth) == 0)  && 
       (!shared || ne(PntPntDist(q1, shared_point), TINY))) {
      *res1 = q1;
      intersections++;
   }

   t = (u2 + d)/2;
   success = IterativeIntersection(ARC, SEG, 0.0, 0.0, 0.0,
                                            C, R, S, d, u2, v1, 0.0,
                                            &t, &q2);

   q3 = VecSub(q2, C);

   if (success  && 
       (IsInArcConeStrict(Ns, Ne, q3) == 0)  && 
       (IsInSegConeStrict(S, q2, a, b, lgth) == 0)  && 
       (!shared || ne(PntPntDist(q2, shared_point), TINY))) {
      if (intersections == 1) {
         if (ne(PntPntDist(q1, q2), eps)) {
            *res2 = q2;
            intersections++;
         }
      } 
      else {
         *res1 = q2;
         intersections++;
      }
   }

   if (intersections > 0)
      IntWarning("seg", seg, "arc", arc);

   return intersections;
}
#endif


vr_bool vroniObject::CheckIntersectionsLocally(int i1, t_site t1, 
                                               int i2, t_site t2,
                                               int i3, t_site t3)
{
   coord res1, res2;
   vr_bool duplicate;

   if (t1 == SEG) {
      if (t2 == SEG) {
         if (GetSegSegIntersection(i1, i2, &res1, eps_int, &duplicate)) {
            SetInvalidSites(i1, SEG, i2, SEG, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
      if (t3 == SEG) {
         if (GetSegSegIntersection(i1, i3, &res1, eps_int, &duplicate)) {
            SetInvalidSites(i1, SEG, i3, SEG, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
#ifdef GENUINE_ARCS
      else if (t2 == ARC) {
         if (GetSegArcIntersection(i1, i2, &res1, &res2, eps_int)) {
            SetInvalidSites(i1, SEG, i2, ARC, eps_int); 
            if (IsInvalidData(false))  return  true;
         }
      }
      else if (t3 == ARC) {
         if (GetSegArcIntersection(i1, i3, &res1, &res2, eps_int)) {
            SetInvalidSites(i1, SEG, i3, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
#endif
   }
#ifdef GENUINE_ARCS
   else if (t1 == ARC) {
      if (t2 == SEG) {
         if (GetSegArcIntersection(i2, i1, &res1, &res2, eps_int)) {
            SetInvalidSites(i1, SEG, i2, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
      else if (t2 == ARC) {
         if (GetArcArcIntersection(i1, i2, &res1, &res2, eps_int, 
                                   &duplicate)) {
            SetInvalidSites(i1, ARC, i2, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
      if (t3 == SEG) {
         if (GetSegArcIntersection(i3, i1, &res1, &res2, eps_int)) {
            SetInvalidSites(i1, SEG, i3, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
      else if (t3 == ARC) {
         if (GetArcArcIntersection(i1, i3, &res1, &res2, eps_int, 
                                   &duplicate)) {
            SetInvalidSites(i1, ARC, i3, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
   }
#endif
   if (t2 == SEG) {
      if (t3 == SEG) {
         if (GetSegSegIntersection(i2, i3, &res1, eps_int, &duplicate)) {
            SetInvalidSites(i2, SEG, i3, SEG, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
#ifdef GENUINE_ARCS
      else if (t3 == ARC) {
         if (GetSegArcIntersection(i2, i3, &res1, &res2, eps_int)) {
            SetInvalidSites(i2, SEG, i3, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
#endif
   }
#ifdef GENUINE_ARCS
   else if (t2 == ARC) {
      if (t3 == SEG) {
         if (GetSegArcIntersection(i3, i2, &res1, &res2, eps_int)) {
            SetInvalidSites(i2, SEG, i3, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
      else if (t3 == ARC) {
         if (GetArcArcIntersection(i2, i3, &res1, &res2, eps_int, 
                                   &duplicate)) {
            SetInvalidSites(i2, ARC, i3, ARC, eps_int);
            if (IsInvalidData(false))  return  true;
         }
      }
   }
#endif

   return  false;
}


/*                                                                           */
/* this function finds all intersections among                               */
/*                                                                           */
/*      pnts[2,num_pnts-3],                                                  */
/*      segs[0,num_segs-1],                                                  */
/*      arcs[0,num_arcs-1].                                                  */
/*                                                                           */
/*                                                                           */
/* a pnt and a seg or arc are reported as intersecting if the pnt lies on    */
/* the seg or arc and if it is neither an end point of the seg or arc nor a  */
/* supporting point of the seg or arc.                                       */
/*                                                                           */
/* note: this routine does not delete duplicates. (all duplicates will be    */
/* discarded within the next call to CleanData().                            */
/*                                                                           */
vr_bool vroniObject::ComputeAllIntersections(double_arg eps)
{
   int i, j;
   coord res1, res2;
#ifdef GENUINE_ARCS
   int cnt;
#endif
   int intersections = 0;
   grid* grid;
   vr_bool ident = false, identical_sites = false;

   VD_Warning("ComputeAllIntersections() - full check for intersections!");

   comp_eps = eps;

   grid = InitialiseGrid(num_pnts + num_segs + num_arcs, bb_min, bb_max);
   FillGrid(grid, eps);

   if (check_completed) {
      check_completed = false;
      if (verbose)
         FP_printf("...checking point intersections with precision %2.1e", 
                FP_PRNTARG(eps));

      for (i = 0; i < grid->n_x; i++) {
         for (j = 0; j < grid->n_y; j++) {
            grid_cell* cell;
            std::set<point_entry*>::iterator it1, it2;
            std::set<seg_entry*>::iterator it3;
#ifdef GENUINE_ARCS
            std::set<arc_entry*>::iterator it4;
            arc_entry *arc;
#endif
            point_entry *p1, *p2;
            seg_entry *seg;

            cell = &grid->cells[i][j];
            for (it1 = cell->points.begin(); it1 != cell->points.end(); it1++){
               p1 = *it1;
               if (IsDeletedPnt(p1->num)) continue;

               //pnt - pnt
               it2 = it1;
               it2++;
               for (; it2 != cell->points.end(); it2++) {
                  p2 = *it2;
                  if (IsDeletedPnt(p2->num)) continue;

                  if (CheckPntPnt(p1->num, p2->num, eps)) {
                     MergePoints(p1->num, p2->num);
                     intersections++;
                  }
               }

               //pnt - seg
               for (it3 = cell->segs.begin(); it3 != cell->segs.end(); it3++) {
                  seg = *it3;

                  if (CheckPntOnSeg(seg->num, p1->num, eps)) {
                     if (AddSegIntersection(seg, GetPntCoords(p1->num)))
                        intersections++;
                  }
               }

   #ifdef GENUINE_ARCS
               //pnt - arc
               for (it4 = cell->arcs.begin(); it4 != cell->arcs.end(); it4++) {
                  arc = *it4;

                  if (CheckPntOnArc(arc->num, p1->num, eps)) {
                     if (AddArcIntersection(arc, GetPntCoords(p1->num)))
                        intersections++;
                  }
               }
   #endif
            }
         }
      }
      if (verbose) {
         printf(" -- done!\n");
         printf("   %d point intersections found!\n", intersections);
      }

      if (intersections > 0) {
         SplitAtIntersections(grid, eps);
         DeleteGrid(grid);
         return true;
      }
   }

   if (verbose)
      FP_printf("...checking seg/arc intersections with precision %2.1e", 
             FP_PRNTARG(eps));

   for (i = 0; i < grid->n_x; i++) {
      for (j = 0; j < grid->n_y; j++) {
         std::set<seg_entry*>::iterator it1, it2;
         seg_entry *seg1, *seg2;
#ifdef GENUINE_ARCS
         std::set<arc_entry*>::iterator it3, it4;
         arc_entry *arc1, *arc2;
#endif
         grid_cell* cell;

         cell = &grid->cells[i][j];

         for (it1 = cell->segs.begin(); it1 != cell->segs.end(); it1++) {
            seg1 = *it1;

            //seg - seg
            it2 = it1;
            it2++;
            for (; it2 != cell->segs.end(); it2++) {
               seg2 = *it2;

               if (GetSegSegIntersection(seg1->num, seg2->num, &res1, eps,
                                         &ident)) {
                  if (AddSegIntersection(seg1, res1) | 
                      AddSegIntersection(seg2, res1))
                     intersections++;
               }
               else if (ident) {
                  ident = false;
                  identical_sites = true;
               }
            }


#ifdef GENUINE_ARCS
            //seg - arc
            for (it3 = cell->arcs.begin(); it3 != cell->arcs.end(); it3++) {
               arc1 = *it3;

               if ((cnt = GetSegArcIntersection(seg1->num, arc1->num, 
                                                &res1, &res2, eps)) > 0) {
                  if ( AddSegIntersection(seg1, res1) | 
                       AddArcIntersection(arc1, res1))
                     intersections++;
                  if (cnt > 1) {
                     assert(cnt == 2);
                     if (AddSegIntersection(seg1, res2) | 
                         AddArcIntersection(arc1, res2))
                        intersections++;
                  }
               }
            }
#endif
         }

#ifdef GENUINE_ARCS
         //arc - arc
         for (it3 = cell->arcs.begin(); it3 != cell->arcs.end(); it3++) {
            arc1 = *it3;

            it4 = it3;
            it4++;
            for (; it4 != cell->arcs.end(); it4++) {
               arc2 = *it4;

               if ((cnt = GetArcArcIntersection(arc1->num, arc2->num, &res1, 
                                                &res2, eps, &ident)) > 0){
                  if (AddArcIntersection(arc1, res1) | 
                      AddArcIntersection(arc2, res1))
                     intersections++;

                  if (cnt > 1) {
                     assert(cnt == 2);
                     if (AddArcIntersection(arc1, res2) | 
                         AddArcIntersection(arc2, res2))
                        intersections++;
                  }
               }
               else if (ident)  {
                  ident = false;
                  identical_sites = true;
               }
            }
         }
#endif
      }
   }

   SplitAtIntersections(grid, eps);
   DeleteGrid(grid);

   check_completed = true;

   if (verbose) {
      printf(" -- done!\n");
      if ((intersections > 0)  &&  identical_sites)
         printf("   %d seg/arc intersections and at least one duplicate!\n", 
                intersections);
      else if (intersections > 0)
         printf("   %d seg/arc intersections!\n", 
                intersections);
      else if (identical_sites)
         printf("   at least one duplicate seg/arc!\n");
      else 
         printf("   no seg/arc intersections and no duplicates!\n");
   }
   
   return ((intersections > 0)  ||  identical_sites);
}