/*****************************************************************************/
/*                                                                           */
/*           Copyright (C)         Martin Held             1999-2025         */
/*                                                                           */
/* This code is provided at no charge to you for non-commercial purposes     */
/* only and only for use internal to your institution. You may use this      */
/* code for academic evaluation and applications, following standard         */
/* rules of academic conduct including crediting the author(s) and the       */
/* copyright holder in any publications.  This code is not in the public     */
/* domain, and may not be duplicated, altered, sold or re-distributed        */
/* without obtaining the prior written consent of the copyright holder.      */
/* No part of this code may be used for or included in any commercial        */
/* application unless your institution obtains a commercial license.         */
/* Please contact the author, Martin Held (held@cs.sbg.ac.at), for           */
/* commercial licensing terms.                                               */
/*                                                                           */
/* This code is provided as is, and you use it at your own risk. The author  */
/* does not accept any responsibility, to the extent permitted by applicable */
/* law, for the consequences of using it or for its usefulness for any       */
/* particular application.                                                   */
/*                                                                           */
/*****************************************************************************/
/*                                                                           */
/* 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                                     */
/*                                                                           */
/*****************************************************************************/
/*                                                                           */
/* Please read the README.txt and README_data.txt files before compiling or  */
/* using this code!                                                          */
/*                                                                           */
/*****************************************************************************/
/*                                                                           */
/* get standard libraries                                                    */
/*                                                                           */
#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <assert.h>
#include <stdlib.h>

/*                                                                           */
/* get my header files                                                       */
/*                                                                           */
#include "fpkernel.h"
#include "vronivector.h"
#include "vroni_object.h"
#include "defs.h"
#include "random.h"
#include "numerics.h"


/*
 * please turn those define's only on if you know what you do! (they are
 * here only for my debugging purposes!!
 *
#define TRACE_VD
#define INIT_FORCED
#define TRACE_INPUT
 *
 */


/*                                                                           */
/* this subroutine drives the computation of the voronoi diagram. the VD is  */
/* computed in an area whose width equals (2*bound+1) times the width of the */
/* bounding box of the input data, and whose height equals (2*bound+1) times */
/* the height of the bounding box of the input data.                         */
/*                                                                           */
/* if ComputeVD() were to be called from a user's application program then   */
/* it is the responsibility of this application program to initialize all    */
/* arrays appropriately. see ReadPolygon(), ReadPoly() or ReadSites() in the */
/* file io.c. It might be best to interface to this subroutine only via the  */
/* API-functions provided; see main().                                       */
/*                                                                           */
void vroniObject::ComputeVD(vr_bool save_data, vr_bool new_input, 
                            int step_size, vr_bool time, int bound, 
                            int sample, int approx, char output_file[], 
                            vr_bool discard_duplicate_sites, int p_step_size, 
                            vr_bool *finished, int a_step_size,
                            vr_bool pnts_only, vr_bool write_vd_dt,
                            char vd_dt_file[], 
                            vr_bool left_vd, vr_bool right_vd,
                            vr_bool clean_up)
{
#ifdef GRAPHICS
   int i;
#endif
   
  int j;
#ifdef TRACE
   vr_bool trace_write_debug = false;
#endif
#ifdef TRACE_INPUT
   vr_bool trace_pnt = true, trace_seg = true, trace_arc = true;
#endif
#ifdef VRONI_INFO
   double eps;
#endif

   /*                                                                        */
   /* make sure that the compile-time options selected by the user are       */
   /* reasonable.                                                            */
   /*                                                                        */
#ifdef INVERSE
#ifndef SORTED
   throw std::runtime_error("VRONI error: ComputeVD() - Inverse sorting requires -DSORTED and -DINVERSE!");
#endif
#endif

#ifdef RANDOM
#ifdef SORTED
   throw std::runtime_error("VRONI error: ComputeVD() - Do not use -DRANDOM and -DSORTED together!");
#else /* else if not SORTED */
#ifdef ORDERED
   throw std::runtime_error("VRONI error: ComputeVD() - Do not use -DRANDOM and -DORDERED together!");
#else /* else if not ORDERED */
#ifdef FORCED
   throw std::runtime_error("VRONI error: ComputeVD() - Do not use -DRANDOM and -DFORCED together!");
#endif /* endif FORCED */
#endif /* endif ORDERED */
#endif /* endif SORTED */
#else /* else if not RANDOM */
#ifdef SORTED
#ifdef ORDERED
   throw std::runtime_error("VRONI error: ComputeVD() - Do not use -DSORTED and -DORDERED together!");
#else /* else if not ORDERED */
#ifdef FORCED
   throw std::runtime_error("VRONI error: ComputeVD() - Do not use -DSORTED and -DFORCED together!");
#endif /* endif FORCED */
#endif /* endif ORDERED */
#else /* else if not SORTED */
#ifdef ORDERED
#ifdef FORCED
   throw std::runtime_error("VRONI error: ComputeVD() - Do not use -DORDERED and -DFORCED together!");
#endif /* endif FORCED */
#endif /* endif ORDERED */
#endif /* endif SORTED */
#endif /* endif RANDOM */
#ifndef RANDOM
#ifndef SORTED
#ifndef ORDERED
#ifndef FORCED
   VD_Warning("ComputeVD() - it is recommended to use -DRANDOM!");
#endif
#endif
#endif
#endif
#ifdef INIT_FORCED
#ifndef RANDOM
   VD_Warning("ComputeVD() - please use -INIT_FORCED only with -DRANDOM!");
#endif
#endif

   /*                                                                        */
   /* actual start of code for inserting sites into a VD                     */
   /*                                                                        */
#ifdef GRAPHICS
   if (graphics) ResetCurBuffer();  /* only if OpenGL display is enabled     */
#else
   step_size   = INT_MAX;
   p_step_size = INT_MAX;
   a_step_size = INT_MAX;
   left_vd     = false;
   right_vd    = false;
#endif

   /*                                                                        */
   /* check whether there seems to be input data; note that every seg or arc */
   /* results in two pnts to handled!                                        */
   /*                                                                        */
   if (num_pnts <= 0) {
      VD_Warning("ComputeVD() - no sites!");
      *finished   = true;
      initialized = false;
      return;
   }
   else {
      *finished = false;
   }

   /*                                                                        */
   /* branch depending on whether this call of ComputeVD() is for            */
   /*  - entirely new data, with no pre-existing VD and no appropriate       */
   /*    preprocessing done yet;                                             */
   /*  - for another round of insertions of sites into an already existing   */
   /*    VD, with preprocessing already done;                                */
   /*  - a restart of the code on data that had already been subjected to    */
   /*    preprocessing once; this may happen if the input data suffers from  */
   /*    deficiencies such as self-intersections.                            */
   /*                                                                        */
   if (restart) {
      if (vr_incr) {
         if (vr_file_handle != ((FILE*)NULL))  fclose(vr_file_handle);
      }
      restart = false;
      if (restart_due_to_intersection) {
         restart_due_to_intersection = false;
         restart_without_intersection = 0;
         ++intersection_counter;
      }
      else {
         ++restart_counter;
         ++restart_without_intersection;
      }
      if (restart_without_intersection > MAX_RESTART) {
         throw std::runtime_error("VRONI error: ComputeVD() - aborting computation -- too many restarts!");
      }
#ifdef VRONI_INFO
      else if (!quiet) {
         printf("warning in ComputeVD() - restart # %d of VD computation!\n", 
                restart_counter + intersection_counter);
      }
#endif
      ExtApplComputeRestarts;
#ifdef GRAPHICS
      ResetBufferData();
#endif
      ResetVDData();
      ResetVDConstructionData();
      ResetFlags();
#ifdef VRONI_INFO
      InitRelaxationMemory(TINY);
#endif

      sample                  = 0;
      new_input               = true;
      discard_duplicate_sites = true;

#ifdef INIT_FORCED
      fprintf(fcd_output, "\n");
      fclose(fcd_output);
#endif

      /*                                                                     */
      /* let the user call some application-specific function                */
      /*                                                                     */
      ExtApplFuncRestart;
   }
   else if (new_input) {
#ifdef VRONI_INFO
      InitRelaxationMemory(TINY);
#endif
      removed_pnts     = removed_segs = removed_arcs = 0;
      restart_counter  = intersection_counter = 0;
      restart_without_intersection = 0;
      bbox_added       = false;
      io_troubles      = false;

      cpu_time_pre     = 0.0;
      cpu_time_pnt     = 0.0;
      cpu_time_seg     = 0.0;
      cpu_time_arc     = 0.0;
      cpu_time_cln     = 0.0;

      /*                                                                     */
      /* we reset the random number generator in order to guarantee that the */
      /* randomized incremental insertion will yield the same result when    */
      /* applying it repeatedly to the same data set.                        */
      /*                                                                     */
      InitRandom(vr_seed); 
      
#ifdef INIT_FORCED
      /* the next three lines might cause troubles on your machine...        */
      system("rm forced_order_pnts.txt");
      system("rm forced_order_segs.txt");
      system("rm forced_order_arcs.txt");
#endif
   }
   else if (!initialized) {
      VD_Warning("ComputeVD() - no input!?");
      *finished = true;
      return;
   }

   if (new_input) {
      /*                                                                     */
      /* preprocessing phase of the VD algorithm                             */
      /*                                                                     */
      VD_Info("...preprocessing -- ");

      if (time)  (void) elapsed();

      num_pnts_deleted = 0;
      initialized   = true;
      new_input     = false;
      restart       = false;
      troubles      = false;
      desperate     = false;
      numerical     = false;
      written       = false;
      done_pnts     = false;
      done_segs     = false;
      done_arcs     = false;
      first_pnts    = true;
      first_segs    = true;
#ifdef GENUINE_ARCS
      first_arcs    = true;
#endif
      incremental   = false;
      
      /*                                                                     */
      /* compute the bounding box of the points pnts[0,...,num_pnts-1].      */
      /*                                                                     */
      if (!bbox_initialized) SetBoundingBox();
      if ((num_arcs > 0) && bbox_initialized && 
          (!bbox_added  || (num_critical_arcs > 0)))  ExtendBoundingBox();

      /*                                                                     */
      /* scale the data if scaling is requested.                             */
      /*                                                                     */
      if (bbox_initialized) {
         if (scale_data && !data_scaled) {
            ScaleData();
#if defined(OGL_GRAPHICS)
            if (graphics  &&  data_scaled)  UpdateScaleData();
#endif         
         }
      }

      /*                                                                     */
      /* enlarge the bounding box according to "bound" by adding four dummy  */
      /* points.                                                             */
      /*                                                                     */
      if (bbox_initialized && !bbox_added) AddDummyCorners(bound);

      /*                                                                     */
      /* if requested, replace the segs/arcs by a dense sample of points,    */
      /* note that we also need to adjust the bounding box if arcs are       */
      /* sampled.                                                            */
      /*                                                                     */
      if ((sample > 0)  &&  ((num_segs > 0)  ||  (num_arcs > 0)))
         SampleData(sample);

      /*                                                                     */
      /* if requested, approximate the circular arcs. note that we also need */
      /* to adjust the bounding box.                                         */
      /*                                                                     */

      if (pnts_only) {
         /*                                                                  */
         /* we'll only compute the VD and DT of points.                      */
         /*                                                                  */
         num_segs      = 0;
         num_arcs      = 0;
         isolated_pnts = true;
      }

#ifndef GENUINE_ARCS
      if (num_arcs > 0) 
         throw std::runtime_error("VRONI error: ComputeVD() - genuine circular arcs not supported!\n");
#endif

#ifdef OTHER
      if (cgal  &&  (restart_counter == 0)) {
         VD_Info("...writing input data for CGAL's VD code --");
         WriteCgalSites("cgal_data.cin");
         VD_Info(" done.\n");
         printf("note that VRONI's CPU timing is incorrect due to this file I/O!\n");
      }
#endif

      /*                                                                     */
      /* sort the vertices in lexicographic order, discard duplicate         */
      /* vertices, and remove zero-length segs/arcs.                         */
      /*                                                                     */
      if (clean_up) discard_duplicate_sites = true;
      CleanData(&removed_pnts, &removed_segs, &removed_arcs, 
                discard_duplicate_sites);
      pnts_deleted = false;

      /*                                                                     */
      /* preprocess the circular arcs                                        */
      /*                                                                     */
      SetPreprocessingThreshold(GetIntersectionThreshold());
      num_critical_arcs = 0;
#ifdef GENUINE_ARCS
      if (num_arcs > 0)  PreprocessArcs(&num_critical_arcs);
#endif
      done_arcs = (num_arcs == 0);

      /*                                                                     */
      /* preprocess the line segments                                        */
      /*                                                                     */
      num_critical_segs = 0;
      if (num_segs > 0)  PreprocessSegments(&num_critical_segs);
      if (restart)  return;
      done_segs = (num_segs == 0);

      /*                                                                     */
      /* build a grid for speeding up the nearest-neighbor search            */
      /*                                                                     */
      BuildGrid();

      if (vr_incr  &&  pnts_only)  vr_file_handle = OpenFileVD(vr_file, "w");

      if (time)  cpu_time_pre += elapsed();

      /*                                                                     */
      /* let the user call some application-specific function                */
      /*                                                                     */
      ExtApplFuncNewInput;

#ifndef NDEBUG
      CheckBoundingBox();
#endif

#ifdef GRAPHICS
      if (graphics  && (p_step_size < INT_MAX))  return;
#endif

      VD_Info("done!\n");
   }
   
   if (!done_pnts) {
      /***********************************************************************/
      /*                                                                     */
      /* compute the Voronoi diagram of the points/vertices.                 */
      /*                                                                     */
      /***********************************************************************/
      if (first_pnts) {
         /*                                                                  */
         /* pnt-specific preprocessing: compute the initial VD and           */
         /* initialize the insertion order of the pnts                       */
         /*                                                                  */
         VD_Info("...computing VD of points -- ");
         if (time)  (void)elapsed();
         if ((num_segs == 0)  &&  (num_arcs == 0))  pnts_only = true;
         first_pnts = false;
         ComputeInitialVD();
         InitSiteOrder(PNT);
         n_pnts = 3;
         max_n_pnts = num_pnts - 2;
         assert(max_n_pnts >= 3);
#ifdef INIT_FORCED
         fcd_output = OpenFileVD("forced_order_pnts.txt", "a");
         fprintf(fcd_output, "restart: %d\n", restart_counter);
#endif
#ifdef TRACE_INPUT
         if (trace_pnt) {
            ClipWriteSiteData(2, PNT, "cwsd.site");
            // printf("2\n%20.16f %20.16f\n", pnts[2].p.x, pnts[2].p.y);
         }
#endif
#ifdef GRAPHICS
         draw_segs = false;
         draw_arcs = false;
         if (graphics) {
            if (p_step_size < INT_MAX)  {
               draw_pnts   = true;
               incremental = true;
               ResetVDBuffer();
               AddVDToBuffer(false, false); 
               AddDelToBuffer();
               return;
            }
            if ((step_size < INT_MAX)  ||  (a_step_size < INT_MAX))  {
               draw_pnts   = true;
            }
         }
#endif

         if (vr_incr  &&  pnts_only) {
            assert(vr_file_handle != ((FILE*)NULL));
            WriteVoronoiRegion(2, vr_file_handle);
         }

         /*                                                                  */
         /* let the user call some application-specific function             */
         /*                                                                  */
         ExtApplFuncFirstPnts;
      }
      
      /*                                                                     */
      /* insert the points  pnts[3,...,num_pnts-3]  into the VD according to */
      /* the insertion order selected.                                       */
      /* recall that   pnts[0,1,2,num_pnts-2,num_pnts-1]  have been inserted */
      /* into the VD by ComputeInitialVD().                                  */
      /*                                                                     */
#ifdef GRAPHICS
      i = 0;
      while ((i < p_step_size)  &&  (n_pnts < max_n_pnts)) {
         ++i;
#else
      while (n_pnts < max_n_pnts) {
#endif
         j = GetNextSite();
#ifdef INIT_FORCED
         fprintf(fcd_output, "%d\n", j);
#endif
#ifdef TRACE_INPUT
         if (trace_pnt) {
            ClipWriteSiteData(j, PNT, "cwsd.site");
            // printf("2\n%20.16f %20.16f\n", pnts[j].p.x, pnts[j].p.y);
         }
#endif
         InsertPntIntoVD(j);
         if (restart)  return; 

         /*                                                                  */
         /* output the newly computed Voronoi region of pnts[j], if required */
         /*                                                                  */
         if (vr_incr  &&  pnts_only) {
            assert(vr_file_handle != ((FILE*)NULL));
            WriteVoronoiRegion(j, vr_file_handle);
         }

         ++n_pnts;
      }

      if (n_pnts == max_n_pnts) {
         /*                                                                  */
         /* all pnts have been inserted into the VD                          */
         /*                                                                  */
         FreeGrid();
         FreeTree();
         FreeReverseIndex();

#ifdef INIT_FORCED
         fprintf(fcd_output, "\n");
         fclose(fcd_output);
#endif
#ifdef GRAPHICS
         if (graphics)  {
            if (p_step_size == INT_MAX)  ResetCurBuffer();
            ResetVDBuffer();
            if ((num_segs == 0) || (p_step_size < INT_MAX)  ||  
                (step_size < INT_MAX))  AddVDToBuffer(false, false);
            AddDelToBuffer();
         }
#endif
         /*                                                                  */
         /* if more than just the VD/DT of points is sought                  */
         /*                                                                  */
         if (!done_segs  ||  !done_arcs) {
            TakeSquareOfNodes();
            InsertDegreeTwoNodes();
         }
         done_pnts = true;

         if (vr_incr  &&  pnts_only) {
            assert(vr_file_handle != ((FILE*)NULL));
            fclose(vr_file_handle);
         }

         if (time)  cpu_time_pnt += elapsed();
         VD_Info("done!\n");
#ifndef NDEBUG
         if (!CheckVD(false)) {
            VD_Warning("ComputeVD() - topological problem in point VD!");
            troubles = true;
         }
#endif

         /*                                                                  */
         /* let the user call some application-specific function             */
         /*                                                                  */
         ExtApplFuncDonePnts;

#ifdef GRAPHICS
         if (graphics  &&  
             ((p_step_size < INT_MAX)  ||  (step_size < INT_MAX)))
            return;
#endif
      }
#ifdef GRAPHICS
      else if (graphics) {
         ResetVDBuffer();
         AddVDToBuffer(false, false);
         AddDelToBuffer();

#ifdef TRACE
         if (grid_used) DrawGrid();
         else           DrawTree();
#endif
         return;
      }
#endif
   }

   if (done_pnts  &&  !done_segs) {
      /***********************************************************************/
      /*                                                                     */
      /* insert the straight-line segments into the VD                       */
      /*                                                                     */
      /***********************************************************************/
      if (first_segs) {
         /*                                                                  */
         /* seg-specific preprocessing: initialize the insertion order of    */
         /* the segs                                                         */
         /*                                                                  */
         VD_Info("...computing VD of segments -- ");
         if (time)  (void)elapsed();
         first_segs = false;
         InitSiteOrder(SEG);
         n_segs = 0;
#ifdef INIT_FORCED
         fcd_output = OpenFileVD("forced_order_segs.txt", "a");
         fprintf(fcd_output, "restart: %d\n", restart_counter);
#endif
#ifdef GRAPHICS
         draw_segs = true;
         if (graphics  &&  (step_size < INT_MAX))
            incremental = true;
         else 
            incremental = false;
#endif
         /*                                                                  */
         /* let the user call some application-specific function             */
         /*                                                                  */
         ExtApplFuncFirstSegs;

         /*if (restart_counter == 5)  SetDebugFlagVD(true);*/
      }

      /*                                                                     */
      /* insert the segments  segs[0,...,num_segs-1]  into the VD in the     */
      /* order selected.                                                     */
      /*                                                                     */
#ifdef GRAPHICS
      i = 0;
      while ((i < step_size)  &&  (n_segs < num_segs)) {
         ++i;
#else
      while (n_segs < num_segs) {
#endif
#ifdef TRACE
         if (trace_write_debug)  WriteProcessedSites("debug.site");
#endif
         j = GetNextSite();
#ifdef INIT_FORCED
         fprintf(fcd_output, "%d\n", j);
#endif
#ifdef TRACE_INPUT
         if (trace_seg) { 
            //printf("0\n%20.16f %20.16f %20.16f %20.16f\n", 
            //       pnts[segs[j].i1].p.x, pnts[segs[j].i1].p.y,
            //       pnts[segs[j].i2].p.x, pnts[segs[j].i2].p.y);
            ClipWriteSiteData(j, SEG, "cwsd.site");
         }
#endif
         InsertSegIntoVD(j);

         if (restart) {
#ifdef TRACE_INPUT
            ClipWriteSiteData(j, UNKNOWN, "cwsd.site");
#endif
#ifdef TRACE_INPUT
            if (trace_seg) WriteSiteData(j, SEG, "restart\n");
#ifdef GRAPHICS
            if (graphics) AddToCurBuffer(j, SEG, AlertColor);
#endif
#endif
            return; 
         }
#ifdef VRONI_INFO
         ConfirmRelaxationMemory();
#endif

         ++n_segs;

#ifndef NDEBUG
         if (!CheckVD(false)) {
            fprintf(stderr, "CheckVD - topological problem in VD! ***\n");
            troubles = true;
         }
#ifdef VRONI_INFO
         else {
            if (!quiet) printf("VD OK after %d segs\n", n_segs);
         }
#endif
#endif
      }

      if (n_segs == num_segs) {
         /*                                                                  */
         /* all segs have been inserted into the VD                          */
         /*                                                                  */
         done_segs = true;
         if (time)  cpu_time_seg += elapsed();
         VD_Info("done!\n");
#ifdef INIT_FORCED
         fprintf(fcd_output, "\n");
         fclose(fcd_output);
#endif
#ifndef NDEBUG
         if (!done_arcs) {
            if (!CheckVD(false)) {
               VD_Warning("ComputeVD() - topological problem in seg VD!");
               troubles = true;
            }
         }
#endif
         /*                                                                  */
         /* let the user call some application-specific function             */
         /*                                                                  */
         ExtApplFuncDoneSegs;
      }

#ifdef GRAPHICS
      if (graphics  &&  (step_size < INT_MAX))  {
         ResetVDBuffer();
         AddVDToBuffer(false, false);
         return;
      }
#endif
   }

	
#ifdef GENUINE_ARCS
   if (done_pnts  &&  done_segs  &&  !done_arcs) {

      /***********************************************************************/
      /*                                                                     */
      /* insert the circular arcs                                            */
      /*                                                                     */
      /***********************************************************************/
      if (first_arcs) {
         /*                                                                  */
         /* arc-specific preprocessing: initialize the insertion order of    */
         /* the arcs                                                         */
         /*                                                                  */
         VD_Info("...computing VD of circular arcs -- ");
         if (time)  (void)elapsed();
         first_arcs = false;
         ResetVDConstructionData();
         InitSiteOrder(ARC);
         n_arcs = 0;
#ifdef INIT_FORCED
         fcd_output = OpenFileVD("forced_order_arcs.txt", "a");
         fprintf(fcd_output, "restart: %d\n", restart_counter);
#endif
#ifdef GRAPHICS
         draw_arcs = true;
         if (graphics  &&  (a_step_size < INT_MAX))
            incremental = true;
         else 
            incremental = false;
#endif
         /*                                                                  */
         /* let the user call some application-specific function             */
         /*                                                                  */
         ExtApplFuncFirstArcs;

#ifdef GRAPHICS
         if (graphics  &&  (a_step_size < INT_MAX))  {
            ResetVDBuffer();
            AddVDToBuffer(false, false);
            return;
         }
#endif
      }

      /*                                                                     */
      /* insert the circular arcs  arcs[0,...,num_arcs-1]  into the VD,      */
      /* either in the order they were stored or in random order.            */
      /*                                                                     */
#ifdef GRAPHICS
      i = 0;
      while ((i < a_step_size)  &&  (n_arcs < num_arcs)) {
         ++i;
#else
      while (n_arcs < num_arcs) {
#endif
#ifdef TRACE
         if (trace_write_debug)  WriteProcessedSites("debug.site");
#endif
         j = GetNextSite();
#ifdef INIT_FORCED
         fprintf(fcd_output, "%d\n", j);
#endif
#ifdef TRACE_INPUT 
         if (trace_arc) {
            // printf("arc %d\n", j);
            // printf("1\n%20.16f %20.16f %20.16f %20.16f %20.16f %20.16f\n", 
            //       pnts[arcs[j].i1].p.x, pnts[arcs[j].i1].p.y,
            //       pnts[arcs[j].i2].p.x, pnts[arcs[j].i2].p.y,
            //       arcs[j].c.x, arcs[j].c.y);
            ClipWriteSiteData(j, ARC, "cwsd.site");
         }
#endif
         InsertArcIntoVD(j);

         if (restart) {
#ifdef TRACE_INPUT
            ClipWriteSiteData(j, UNKNOWN, "cwsd.site");
#endif
#ifdef TRACE_INPUT
            if (trace_arc) WriteSiteData(j, ARC, "restart\n");
#ifdef GRAPHICS
            if (graphics)  AddToCurBuffer(j, ARC, AlertColor);
#endif
#endif
            return;
         }
#ifdef VRONI_INFO
         ConfirmRelaxationMemory();
#endif

#ifndef NDEBUG
         if (!CheckVD(false)) {
            VD_Warning("CheckVD - topological problem in VD");
            troubles = true;
         }
#ifdef VRONI_INFO
         else {
            if (!quiet) printf("VD OK after %d arcs\n", n_arcs + 1);
         }
#endif
#endif
         ++n_arcs;
      }

      if (n_arcs == num_arcs) {
         /*                                                                  */
         /* all arcs have been inserted into the VD                          */
         /*                                                                  */
         done_arcs = true;
         if (time)  cpu_time_arc += elapsed();
         VD_Info("done!\n");

         /*                                                                  */
         /* let the user call some application-specific function             */
         /*                                                                  */
         ExtApplFuncDoneArcs;

#ifdef INIT_FORCED
         fprintf(fcd_output, "\n");
         fclose(fcd_output);
#endif
      }

#ifdef GRAPHICS
      if (graphics  &&  (a_step_size < INT_MAX)) {
         ResetVDBuffer();
         AddVDToBuffer(false, false);
         return;
      }
#endif
   }
#else
   if (!done_arcs) {
      throw std::runtime_error("VRONI error: ComputeVD() - circular arcs cannot be handled without approximation!");
   }
#endif

   if (done_pnts  &&  done_segs  &&  done_arcs) {

      if (desperate  &&  (restart_counter < MAX_RESTART)) {
         restart = true;
         return;
      }

      /*                                                                     */
      /* do we have isolated points?                                         */
      /*                                                                     */
      if ((num_segs > 0)  ||  (num_arcs > 0)) {
         if (!GetDataSorted()) {
            CheckIsolatedPnts();
         }
      }
      else {
         isolated_pnts = true;
      }

#ifndef NDEBUG
      if (!CheckVD(false)) {
         VD_Warning("ComputeVD() - topological problem(s) in VD!");
      }
      else {
         VD_Info("\nComputeVD() - VD is topologically correct!");
      }
#endif

      SetVDStatusTrue();

#ifdef GRAPHICS
      if (graphics  &&  ((num_segs > 0)  ||  (num_arcs > 0)))  {
         draw_pnts   = false;
         incremental = false;
         ResetVDBuffer();
         if (step_size == INT_MAX)  ResetCurBuffer();
         AddVDToBuffer(left_vd, right_vd);
      }
#endif

      /*                                                                     */
      /* let the user call some application-specific function                */
      /*                                                                     */
      ExtApplFuncDoneVD;
      
      if (!written) {
         if (save_data) {
            VD_Info("\n...writing the output data -- ");
            WriteSites(output_file);
            written = true;
            VD_Info("done!\n");
         }
         if (pnts_only  &&  write_vd_dt) {
            VD_Info("\n...writing VD/DT data -- ");
            WriteVDDT(vd_dt_file);
            VD_Info("done!\n");
         }
         written = true;
      }

#ifdef VRONI_INFO
      if (!quiet) {
         if (IsCheckComplete()  ||  (intersection_counter > 0)) {
            eps = GetIntersectionThreshold();
            if (IsCheckComplete())  eps /= 10.0;
         }
         else {
            eps = 0.0;
         }
         printf("\nComputeVD() - intersection threshold used: %2.1e\n", eps);
#ifdef VRONI_DBG_WARN
         printf("ComputeVD() - maximum relaxation used:     %2.1e\n", 
                ReportRelaxationMemory());
#endif
         if ((restart_counter + intersection_counter) > 0) {
            printf("\nComputeVD() - problems with the input data caused %d restart(s)!",
                   restart_counter + intersection_counter);
         }
         else if (io_troubles) {
            VD_Warning("ComputeVD() - problems with input data occurred!");
         }
         if (desperate) {
            VD_Warning("ComputeVD() - VD computed; desperate mode was used!");
         }
         else if (troubles) {
            VD_Warning("ComputeVD() - VD computed; problems occurred!?");
         }
         else if (numerical) {
            VD_Warning("ComputeVD() - VD computed; numerical instabilities!?");
         }
         else {
            VD_Info("\nComputeVD() - VD computed successfully!\n");
         }
         printf("\n");
      }
#endif
   }

   *finished   = true;
   initialized = false;

   return;
}