//****************************************************************************
// File: geometr.cc
// Operations such as determine the minimal bounding rectangle of a polygon.
//****************************************************************************

#include <math.h>
#include <stdio.h>
#include <std.h>

#include "geometr.h"


void geometr::len_vector(point v, coord *length)
// Calculates the length of a vector.
// 2D implementation only!
{
   *length = (coord)sqrt((double)v[0]*(double)v[0]+
                           (double)v[1]*(double)v[1]);
}


void geometr::inprod(point pt1, point vect, coord *c)
// Inproduct of point and vector, result is returned in c;
{
    int d;

    *c = (double)pt1[0] * (double)vect[0];
    for (d = 1; d < DIM; d++)
        *c += (double)pt1[d] * (double)vect[d];
}


boolean geometr::point_left_plane(point pt, point n_vector, coord *c)
// Determines whether a point is to the left of a hyperplane.
{
    coord c1;

    inprod(pt, n_vector, &c1);
    return (c1 <= *c + EPS);
}


boolean geometr::point_right_plane(point pt, point n_vector, coord *c)
// Determines whether a point is to the right of a hyperplane.
{
    coord c1;

    inprod(pt, n_vector, &c1);
    return (c1 >= *c - EPS);
}


boolean geometr::rect_left_plane(point rect[2], point n_vec, coord *lc)
// Determines whether a hyperrectangle is to the left of a hyperplane.
{
    point	corner;
    boolean	pt_left;
    int		nr[DIM], d;

    for (d = 0; d < DIM; d++) nr[d] = 0;
    do {
        for (d = 0; d < DIM; d++) corner[d] = rect[nr[d]][d];
        pt_left = point_left_plane(corner, n_vec, lc);
        if (!pt_left) {
            // Next corner of the hyperrectangle.
            for (d = 0; (d < DIM) && (nr[d] == 1); d++) nr[d] = 0;
            if (d < DIM) nr[d] = 1;
        }
    } while (!pt_left && (d < DIM));
    return pt_left;
}


boolean geometr::rect_right_plane(point rect[2], point n_vec, coord *rc)
// Determines whether a hyperrectangle is to the right of a hyperplane.
{
    point	corner;
    boolean	pt_right;
    int		nr[DIM], d;

    for (d = 0; d < DIM; d++) nr[d] = 0;
    do {
        for (d = 0; d < DIM; d++) corner[d] = rect[nr[d]][d];
        pt_right = point_right_plane(corner, n_vec, rc);
        if (!pt_right) {
            // Next corner of the hyperrectangle.
            for (d = 0; (d < DIM) && (nr[d] == 1); d++) nr[d] = 0;
            if (d < DIM) nr[d] = 1;
        }
    } while (!pt_right && (d < DIM));
    return pt_right;
}


boolean geometr::sphere_left_plane(point pt, coord *rad, point norm, coord *c)
// Determines whether a hypersphere is to the left of a hyperplane.
{
    coord c1, len_n;

    inprod(pt, norm, &c1);
    len_vector(norm, &len_n);
    c1 -= (*rad) * len_n;
    return (c1 <= *c + EPS);
}


boolean geometr::sphere_right_plane(point pt, coord *rad, point norm, coord *c)
// Determines whether a hypersphere is to the right of a hyperplane.
{
    coord c1, len_n;

    inprod(pt, norm, &c1);
    len_vector(norm, &len_n);
    c1 += (*rad) * len_n;
    return (c1 >= *c - EPS);
}


boolean geometr::pol_left_plane(int nop, point *pts, point norm, coord *c)
{
   boolean left_side = false;

   for (int i=0; !left_side && (i<nop); i++)
      left_side = point_left_plane(pts[i], norm, c);
   return left_side;
}


boolean geometr::pol_right_plane(int nop, point *pts, point norm, coord *c)
{
   boolean right_side = false;

   for (int i=0; !right_side && (i<nop); i++)
      right_side = point_right_plane(pts[i], norm, c);
   return right_side;
}


boolean geometr::left_touched(point mid, coord *rad, point n_vec, coord *c)
{
    coord c1, len_n;

    inprod(mid, n_vec, &c1);
    len_vector(n_vec, &len_n);
    c1 += (*rad) * len_n;
    return (c1 == *c);
}


boolean geometr::right_touched(point mid, coord *rad, point n_vec, coord *c)
{
    coord c1, len_n;

    inprod(mid, n_vec, &c1);
    len_vector(n_vec, &len_n);
    c1 -= (*rad) * len_n;
    return (c1 == *c);
}


void geometr::rect_cpy(point rects[2], point rectd[2])
// Copy the first rectangle to the second rectangle.
{
   for (int d = 0; d < DIM; d++) {
      rectd[0][d] = rects[0][d];
      rectd[1][d] = rects[1][d];
   }
}


void geometr::pt_cpy(point pts, point ptd)
// Copy the first point to the second point.
{
   for (int d = 0; d < DIM; d++) {
      ptd[d] = pts[d];
   }
}


void geometr::cpy_to_pol(point rectangle[2], point *pol)
// Copy the rectangle in polyline format.
// 2D implementation only.
{
   pt_cpy(rectangle[0], pol[0]);
   pol[1][0] = rectangle[1][0];
   pol[1][1] = rectangle[0][1];
   pt_cpy(rectangle[1], pol[2]);
   pol[3][0] = rectangle[0][0];
   pol[3][1] = rectangle[1][1];
}


void geometr::to_rectangle(point *polygon, int nop, point rect[2])
// Determine the minimal bounding rectangle of a polygon or a polyline.
{
   coord min, max;

   for (int d = 0; d < DIM; d++) {
      min = max = polygon[0][d];
      for (int i = 1; i < nop; i++)
    if (polygon[i][d] < min)
       min = polygon[i][d];
    else if (polygon[i][d] > max)
       max = polygon[i][d];
      rect[0][d] = min;
      rect[1][d] = max;
   }
}


void geometr::calc_rect(point rects1[2], point rects2[2], point rectd[2])
// Determine the minimal bounding rectangle of two rectangles.
{
   for (int d = 0; d < DIM; d++) {
      rectd[0][d] = rects1[0][d] < rects2[0][d] ? rects1[0][d] : rects2[0][d];
      rectd[1][d] = rects1[1][d] > rects2[1][d] ? rects1[1][d] : rects2[1][d];
   }
}


coord geometr::calc_area(point rect[2])
// Calculate area of the rectangle.
{
   coord temp_area = rect[1][0] - rect[0][0];

   for (int d = 1; d < DIM; d++)
      temp_area *= rect[1][d] - rect[0][d];
   return temp_area;
}

coord geometr::distance(point p1, point p2)
// Distance between two points.
{
   double d=0.0;
   for (int i=0;i<DIM;i++) d += (double)sqr((p1[i]-p2[i]));
   d = sqrt(d);
   return((coord)d);
}


boolean geometr::cover(point rect1[2], point rect2[2])
// Calculate whether the first rectangle covers the second rectangle.
{
   boolean inside_rect1 = true;

   for (int d = 0; inside_rect1 && (d < DIM); d++)
      if (rect2[0][d] < (rect1[0][d]) || (rect2[1][d] > rect1[1][d]))
    inside_rect1 = false;
   return inside_rect1;
}

boolean geometr::cover(point pt1, coord r1, point pt2, coord r2)
// Calculate whether the first sphere covers the second sphere.
{
   if ((distance(pt1, pt2) -(r1-r2))> (RELEPS*r1)) return false;
      else return true;
}


boolean geometr::overlap(point rect1[2], point rect2[2])
// Calculate whether the first rectangle overlaps the second rectangle.
{
   boolean intersect = true;

   for (int d = 0; intersect && (d < DIM); d++)
      if ((rect1[1][d] < rect2[0][d]) || (rect1[0][d] > rect2[1][d]))
         intersect = false;
   return intersect;
}


boolean geometr::overlap(point pt1, coord r1, point pt2, coord r2)
// Calculate whether two spheres overlap.
{
// if (distance(pt1,pt2) > (r1+r2)) return false; else return true;
   if (((pt1[0]-pt2[0])*(pt1[0]-pt2[0])+(pt1[1]-pt2[1])*(pt1[1]-pt2[1])) > (r1+r2)*(r1+r2)) return false; else return true; //2D only!
}


line_situation geometr::evaluate_rect_line(point p1, point p2, point rect[2])
// 2D implementation only!
{
   double dist;
   point n, pt[4], r[2];
   line_situation ls=IN;

   cpy_to_pol(rect, pt);

   n[0]= ( p1[1] - p2[1]);
   n[1]= ( p2[0] - p1[0]);
   for (int i=0; ((i<4)&&(ls!=OUT)); i++) {
      dist = (n[0]*(p1[0]-pt[i][0])) + (n[1]*(p1[1]-pt[i][1]));
      if (dist<0) {if ((ls==IN)||(ls==LEFT)) ls=LEFT; else ls=OUT;}
         else {if ((ls==IN)||(ls==RIGHT)) ls=RIGHT; else ls=OUT;}
   }
   if (ls!=OUT) return ls;
   if (p1[1]<p2[1]) {r[0][1]=p1[1];r[1][1]=p2[1];}
      else {r[0][1]=p2[1];r[1][1]=p1[1];}
   if (p1[0]<p2[0]) {r[0][0]=p1[0];r[1][0]=p2[0];}
      else {r[0][0]=p2[0];r[1][0]=p1[0];}
   if (overlap(rect, r)) return IN; else return OUT;
}


line_situation geometr::evaluate_spl_line(point p1,point p2, point pc,coord r)
// 2D implementation only!
{
   double a, b, c, d, l1, l2, dist, len;
   point n;

   a = sqr((p2[0]-p1[0])) + sqr((p2[1]-p1[1]));
   len = sqrt(a);
   n[0]= ( p1[1] - p2[1])/len;
   n[1]= ( p2[0] - p1[0])/len;
   dist = (n[0]*(p1[0]-pc[0])) + (n[1]*(p1[1]-pc[1]));
   if (dist<0) {if (r< -dist) return LEFT;}
     else {if (r< dist) return RIGHT;}

   b = 2*(((p2[0]-p1[0])*(p1[0]-pc[0]))+((p2[1]-p1[1])*(p1[1]-pc[1])));
   c = sqr((p1[0]-pc[0])) + sqr((p1[1]-pc[1])) - sqr(r);
   d = sqr(b) - 4*a*c;
   if (d > 0) {
      l1 = (-b + sqrt(d))/(2*a);
      l2 = (-b - sqrt(d))/(2*a);
      if (((l1 >= 0)&&(l1 <= 1)) || ((l2 >= 0)&&(l2 <= 1))) return IN;
         else if (((l1 < 0)&&(l2 > 1)) || ((l2 < 0)&&(l1 > 1))) return IN;
            else return OUT;
   }
   fprintf(stderr, "Problem in geometr::evaluate_spl_line\n");
   abort();
   return IN;
}


line_situation geometr::evaluate_pt_line(point p1, point p2, point pc)
// 2D implementation only!
{
   double dist;
   point n;

   n[0]= ( p1[1] - p2[1]);
   n[1]= ( p2[0] - p1[0]);
   dist = (n[0]*(p1[0]-pc[0])) + (n[1]*(p1[1]-pc[1]));
   if (dist<0) return LEFT; else return RIGHT;
}


boolean geometr::overlap(int nop, point *pts, point pt, coord r)
// Calculates whether a convex polygon has overlap with a sphere.
{
   boolean out_intersect=false;
   for (int i=0; i<nop; i++)
      switch (evaluate_spl_line(pts[i],pts[(i+1==nop)?0:i+1], pt, r)) {
      case RIGHT: return false;
      case IN:    return true;
      case OUT:   out_intersect=true; break;
      case LEFT:  break;
      }

   if (out_intersect) return false; else return true;
}


boolean geometr::overlap(int nop, point *pts, point rect[2]) 
// Calculates whether a convex polygon has overlap with a rectangle.
{
   boolean out_intersect=false;
   for (int i=0; i<nop; i++) {
      switch (evaluate_rect_line(pts[i],pts[(i+1==nop)?0:i+1], rect)) {
      case RIGHT: return false;
      case IN:    return true;
      case OUT:   out_intersect=true; break;
      case LEFT:  break;
      }
   }

   if (out_intersect) return false; else return true;
}


boolean geometr::point_inside(point rect[2], point pt)
// Calculate whether the point lies inside the rectangle.
{
   boolean between = true;

   for (int d = 0; between && (d < DIM); d++)
      if ((pt[d] < rect[0][d]) || (rect[1][d] < pt[d]))
    between = false;
   return between;
}


boolean geometr::point_inside(point pt, point mid_pt, coord radius)
// Calculate whether the point lies inside the sphere.
{
//    if ((distance(pt,mid_pt)-radius) <= (RELEPS*radius)) return true;
//       else return false;
//  Fast 2D-implementation:
    if ((((pt[0]-mid_pt[0])*(pt[0]-mid_pt[0]))+((pt[1]-mid_pt[1])*(pt[1]-mid_pt[1]))) <= ((1+RELEPS)*radius)*((1+RELEPS)*radius)) return true;
       else return false;
}


boolean geometr::point_inside(point pt, int nop, point *pts)
// Calculate whether point lies in convex polygon.
{
   for (int i=0; i<nop; i++)
      switch (evaluate_pt_line(pts[i],pts[(i+1==nop)?0:i+1], pt)) {
      case RIGHT: return false;
      }
   return true;
}


boolean geometr::overlap(point rect[2], point pt, coord r)
// Calculate whether there is overlap between a rectangle and a sphere.
// 2D implementation only !
{
   boolean overlap = true, cornertest = false;
   point corner;

   for (int d = 0; overlap && (d < DIM); d++)
      if ((pt[d] < (rect[0][d]-r)) || (pt[d] > (rect[1][d]+r)))
         overlap = false;

   if (overlap && (pt[0] < rect[0][0])) {
      corner[0] = rect[0][0];
      if (pt[1] < rect[0][1]) {
         cornertest = true;
         corner[1] = rect[0][1];
      } else if (pt[1] > rect[1][1]) {
         cornertest = true;
         corner[1] = rect[1][1];
      }
   } else if (overlap && (pt[0] > rect[1][0])) {   
      corner[0] = rect[1][0];
      if (pt[1] < rect[0][1]) {
         cornertest = true;
         corner[1] = rect[0][1];
      } else if (pt[1] > rect[1][1]) {
         cornertest = true;
         corner[1] = rect[1][1];
      }
   }

   if (cornertest) {
      if (((pt[0]-corner[0])*(pt[0]-corner[0])+
           (pt[1]-corner[1])*(pt[1]-corner[1])) > (r*r)) overlap = false;
   }

   return overlap;
}


boolean geometr::point_equal(point pt1, point pt2)
// Test whether two points are equal.
{
    boolean equal = true;

    for (int d = 0; (d < DIM) && equal; d++)
        if (pt1[d] != pt2[d])
            equal = false;
    return equal;
}


boolean geometr::sphere_equal(point m_pt1, coord rad1, point m_pt2, coord rad2)
// Test whether two spheres are equal.
{
    return (point_equal(m_pt1, m_pt2) && (rad1 == rad2));
}


boolean geometr::rect_equal(point rect1[2], point rect2[2])
// Test whether two rectangles are equal.
{
   boolean equal = true;

   for (int d = 0; equal && (d < DIM); d++)
      if (rect1[0][d] != rect2[0][d] || rect1[1][d] != rect2[1][d])
         equal = false;
   return equal;
}


void geometr::diff_vect(point p1, point p2, point p3)
{
    for (int d = 0; d < DIM; d++)
        p3[d] = (coord) ((double)p1[d] - (double)p2[d]);
}


void geometr::intersect(point rectangle[2], point n, coord c,
                        double_point seg[2])
// Determine intersection of a rectangle and a line.
{
   if (n[0] < n[1]) {
      seg[0][1] = (double) rectangle[0][1];
      seg[0][0] = (double) (c - n[1] * seg[0][1]) / n[0];
      if (seg[0][0] < (double) rectangle[0][0]) {
         seg[0][0] = (double) rectangle[0][0];
         seg[0][1] = (double) (c - n[0] * seg[0][0]) / n[1];
      }
      seg[1][1] = (double) rectangle[1][1];
      seg[1][0] = (double) (c - n[1] * seg[1][1]) / n[0];
      if (seg[1][0] > (double) rectangle[1][0]) {
         seg[1][0] = rectangle[1][0];
         seg[1][1] = (double)(c - n[0] * seg[1][0]) / n[1];
      }
   }
   else {
      seg[0][0] = rectangle[0][0];
      seg[0][1] = (double) (c - n[0] * seg[0][0]) / n[1];
      if (seg[0][1] < (double) rectangle[0][1]) {
         seg[0][1] = rectangle[0][1];
         seg[0][0] = (double)(c - n[1] * seg[0][1]) / n[0];
      }
      seg[1][0] = rectangle[1][0];
      seg[1][1] = (double) (c - n[0] * seg[1][0]) / n[1];
      if (seg[1][1] > (double) rectangle[1][1]) {
         seg[1][1] = rectangle[1][1];
         seg[1][0] = (double)(c - n[1] * seg[1][1]) / n[0];
      }
   }
}


void geometr::intersect(point b, point e, point n, coord *c, point pt)
{
   double labda;

   labda = ((double)*c-(double)b[0]*n[0]-(double)b[1]*n[1])/
           (((double)e[0]-b[0])*(double)n[0]+((double)e[1]-b[1])*(double)n[1]);
   pt[0] = (coord)(labda*((double)e[0]-b[0])+(double)b[0]);
   pt[1] = (coord)(labda*((double)e[1]-b[1])+(double)b[1]);
}


boolean geometr::intersect(point n, coord *c, int nop, point *pol, side_type s,
                           int &new_nop, point* &new_pol, point seg[2])
// Determine intersection of a polygon and a line.
{
   point	pt;
   int		i, j;
   int		intersecting = 0;
   side_type    prev_time;

   if ((new_pol = new point[nop+2]) == NULL) {
      printf("geometr::intersect: Out of memory.\n");
      exit(1);
   }
   new_nop = 0;
   if ((s == left) ?
        point_left_plane(pol[0], n, c) : point_right_plane(pol[0], n, c)) {
      pt_cpy(pol[0], new_pol[new_nop]);
      new_nop++;
      prev_time = s;
   }
   else {	// Wrong side.
      if (s == left)
         prev_time = right;
      else
         prev_time = left;
   }
   for (i = 1; i < nop; i++) {
      if ((s == left) ?
          point_left_plane(pol[i], n, c) : point_right_plane(pol[i], n, c)) {
         // Right side.
         if (s != prev_time) {	// Line intersects line segment.
            intersect(pol[i-1], pol[i], n, c, pt);
            pt_cpy(pt, seg[intersecting]);
            intersecting++;
            pt_cpy(pt, new_pol[new_nop]);
            new_nop++;
         }
         pt_cpy(pol[i], new_pol[new_nop]);
         new_nop++;
         prev_time = s;
      }
      else {	// Wrong side.
         if (s == prev_time) {	// Line intersects line segment.
            intersect(pol[i-1], pol[i], n, c, pt);
            pt_cpy(pt, seg[intersecting]);
            intersecting++;
            pt_cpy(pt, new_pol[new_nop]);
            new_nop++;
            if (s == left)
               prev_time = right;
            else
               prev_time = left;
         }
      }
   }
   if ((s == left) ?
        point_left_plane(pol[0], n, c) :
        point_right_plane(pol[0], n, c)) {	// Check seg(0,nop-1)
      // Right side.
      if (s != prev_time) {	// Line intersects line segment.
         intersect(pol[nop-1], pol[0], n, c, pt);
         pt_cpy(pt, seg[intersecting]);
         intersecting++;
         pt_cpy(pt, new_pol[new_nop]);
         new_nop++;
      }
   }
   else {	// Wrong side
      if (s == prev_time) {	// Line intersects line segment.
         intersect(pol[nop-1], pol[0], n, c, pt);
         pt_cpy(pt, seg[intersecting]);
         intersecting++;
         pt_cpy(pt, new_pol[new_nop]);
         new_nop++;
         if (s == left)
            prev_time = right;
         else
            prev_time = left;
      }
   }
   if (intersecting > 2) {
      printf("number of intersecting points = %d\n", intersecting);
      exit(1);
   }
   return (intersecting > 0) ? true : false;
}


// Utility functions for calculation minimal bounding sphere
// of a set of points (2D only).
// Copy the first point to the second double_point. **************************
void pt_cpy2d(point pts, double_point ptd)
{
   for(int d = 0; d < DIM; d++) ptd[d] = (double) pts[d];
}

// Copy the first double_point to the second point.
void pt_cpy2s(double_point pts, point ptd)
{
   for(int d = 0; d < DIM; d++) ptd[d] = (coord) pts[d];
}

// Copy the first double_point to the second double_point.
void pt_cpyd(double_point pts, double_point ptd)
{
   for(int d = 0; d < DIM; d++) ptd[d] = pts[d];
}

void InitSolution(point *pts, int nop, double_point solution[DIM+2])
{
   // Copy the first four points.
   for (int i = 0; (i < DIM + 2) && (i < nop); i++)
      pt_cpy2d(pts[i], solution[i]);
   for (int j = i; j < DIM + 2; j++)
      pt_cpy2d(pts[nop - 1], solution[j]);
}

boolean InSphere(double_point c, double *r, // sphere
                 double_point curr,         // point
                 double *rr)                // distance between M & P
{
   *rr = 0.0;
   for (int d = 0; d < DIM; d++) 
      *rr += (c[d] - curr[d])*(c[d] - curr[d]);
   *rr = sqrt(*rr);
   return ((*rr - *r) <= ((*r)* RELEPS)) ? true : false;
}

void Try2(double_point p1, double_point p2, double_point c, double *r)
{
   c[0] = (p1[0] + p2[0]) / 2.0;
   c[1] = (p1[1] + p2[1]) / 2.0;
   *r = sqrt((c[0]-p1[0])*(c[0]-p1[0]) + (c[1]-p2[1])*(c[1]-p2[1]));
//if (*r == 0.0) printf("geometr Try2: danger (double points)\n");
}

void Try3(double_point p1, double_point p2, double_point p3,    // Points
          double_point c, double *r)                     // Sphere
{
   double teller, noemer;

   noemer = 2.0 * ((p3[0]-p1[0])*(p2[1]-p1[1]) + (p3[1]-p1[1])*(p1[0]-p2[0]));
   teller = (p3[0]-p1[0]) * (p3[0]-p2[0]) + (p3[1]-p1[1]) *(p3[1]-p2[1]);
   if (noemer != 0.0) {
      c[0] = (teller*(p2[1]-p1[1]))/noemer + (p1[0]+p2[0])/ 2.0;
      c[1] = (teller*(p1[0]-p2[0]))/noemer + (p1[1]+p2[1])/ 2.0;
      *r = sqrt((c[0]-p3[0])*(c[0]-p3[0]) + (c[1]- p3[1])*(c[1]- p3[1]));
   }
   else {
      pt_cpyd(p2, c);
//printf("geometr Try3: danger (double points)\n");
      *r = 0.0;
   }
}

void CurrSphere(double_point sol[DIM+2], double_point c, double *s, 
                int *notused)
{
   double d; /* dummy, not really used */
   double hs;
   double_point hc;
   boolean found;

   if (DIM != 2) {
      fprintf(stderr, "Only 2 DIM case is implemented\n");
      exit(1);
   }
   Try2(sol[0], sol[1], c, s);
   *notused = 2;
   if (InSphere(c, s, sol[2], &d) && InSphere(c, s, sol[3], &d))
      return;
   Try2(sol[0], sol[2], c, s);
   *notused = 1;
   if (InSphere(c, s, sol[1], &d) && InSphere(c, s, sol[3], &d))
      return;
   Try2(sol[0], sol[3], c, s);
   *notused = 1;
   if (InSphere(c, s, sol[1], &d) && InSphere(c, s, sol[2], &d))
      return;
   Try2(sol[1], sol[2], c, s);
   *notused = 0;
   if (InSphere(c, s, sol[0], &d) && InSphere(c, s, sol[3], &d))
      return;
   Try2(sol[1], sol[3], c, s);
   *notused = 0;
   if (InSphere(c, s, sol[0], &d) && InSphere(c, s, sol[2], &d))
      return;
   Try2(sol[2], sol[3], c, s);
   *notused = 0;
   if (InSphere(c, s, sol[0], &d) && InSphere(c, s, sol[1], &d))
      return;
   found = false;
// printf("Try2 has not produced solution, continue with Try3\n");
   Try3(sol[0], sol[1], sol[2], hc, &hs);
   if (InSphere(hc, &hs, sol[3], &d)) {
      *s = hs;
      pt_cpyd(hc, c);
      *notused = 3;
      found = true;
   }
   Try3(sol[0], sol[1], sol[3], hc, &hs);
   if (InSphere(hc, &hs, sol[2], &d)) {
      if (!found || (hs < *s)) {
         *s = hs;
         pt_cpyd(hc, c);
         *notused = 2;
         found = true;
      }
   }
   Try3(sol[0], sol[2], sol[3], hc, &hs);
   if (InSphere(hc, &hs, sol[1], &d)) {
      if (!found || (hs < *s)) {
         *s = hs;
         pt_cpyd(hc, c);
         *notused = 1;
         found = true;
      }
   }
   Try3(sol[1], sol[2], sol[3], hc, &hs);
   if (InSphere(hc, &hs, sol[0], &d)) {
      if (!found || (hs < *s)) {
         *s = hs;
         pt_cpyd(hc, c);
         *notused = 0;
         found = true;
      }
   }
   if (!found) {
      fprintf(stderr, "CurrSphere: No sphere found !!!\n");
      exit(1);
   }
}

boolean SphereOK(point *pts, int nop,      // points
                 double_point c, double *r,   // sphere
                 double_point outside)             // point outside the sphere
{
   boolean     OK;
   double rr, max;
   double_point dummy;

   max = 0.0;
   OK = true;
   for (int i = 0 ; i < nop; i++) {
      pt_cpy2d(pts[i], dummy);
      if (!InSphere(c, r, dummy, &rr)) {
         OK = false;
         if (rr > max) {
            pt_cpyd(dummy, outside);
            max = rr;
         }; 
      }
   }
   return (OK);
}


void geometr::to_sphere(point *pts, int nop, point mid_pt, coord &radius)
// Determine the minimal bounding sphere of a polygon or a polyline.
{
   double_point  solution[DIM+2], outside, mid;
   int          notused;
   double  newr, oldr = 0.0;         // radius

   if (DIM != 2) {
      fprintf(stderr, "Only for 2d case implemented.\n");
      exit(1);
   }
   if (nop <= 0) {
      fprintf(stderr, "Object must contain at least 1 point.\n");
      exit(1);
   }
   InitSolution(pts, nop, solution);
   CurrSphere(solution, mid, &newr, &notused);
   while (!SphereOK(pts,nop,mid,&newr,outside)&&
      ((newr-oldr) > (RELEPS*newr))) {
      oldr = newr;
      pt_cpyd(outside, solution[notused]);
      CurrSphere(solution, mid, &newr, &notused);
   }
   radius = (coord)(newr);
   pt_cpy2s(mid, mid_pt);
}