/*******************************************************************************
+
+  LEDA 3.5
+
+  closest_pair.c
+
+  This file is part of the LEDA research version (LEDA-R) that can be 
+  used free of charge in academic research and teaching. Any commercial
+  use of this software requires a license which is distributed by the
+  LEDA Software GmbH, Postfach 151101, 66041 Saarbruecken, FRG
+  (fax +49 681 31104).
+
+  Copyright (c) 1991-1997  by  Max-Planck-Institut fuer Informatik
+  Im Stadtwald, 66123 Saarbruecken, Germany     
+  All rights reserved.
+ 
*******************************************************************************/
//------------------------------------------------------------------------------
// C. Uhrig & S. N"aher (1995)
//------------------------------------------------------------------------------

#include<LEDA/array.h>
#include <math.h>
#include <assert.h>

static int cmp_y(const POINT& a, const POINT& b)
{ return compare(a.ycoord(),b.ycoord()); }

inline COORD square(COORD x) { return x*x; }


static void x_search(POINT* first, POINT* last, COORD& dist, 
                                                POINT&  r1,
                                                POINT&  r2)
{ 
  for (POINT* p = first; p < last; p++)
  { POINT* q = p+1;
    while (q <= last && square(p->xdist(*q)) < dist)
    { COORD d = p->sqr_dist(*q);
      if (d < dist)
      { dist = d;
        r1 = *p;
        r2 = *q;
       }
      q++;
     }
   }

 }


static void y_search(POINT* first, POINT* last, COORD& dist,
                                                POINT&  r1,
                                                POINT&  r2)
{ 
  for (POINT* p = first; p < last; p++)
  { POINT* q = p+1;
    while (q <= last && square(p->ydist(*q)) < dist)
    { COORD d = p->sqr_dist(*q);
      if (d < dist)
      { dist = d;
        r1 = *p;
        r2 = *q;
       }
      q++;
     }
   }

 }




static void cl_pair(POINT* start, POINT* stop, POINT& r1,
                                               POINT& r2, 
                                               COORD& dist, 
                                               array<POINT>& Ay)
{ if (dist == 0) return;

  int n = stop-start+1;

  if (n < 30) 
     x_search(start,stop,dist,r1,r2);
  else
     { POINT* mid = start + n/2 - 1;

       cl_pair(start, mid, r1, r2, dist, Ay); 
       cl_pair(mid+1, stop, r1, r2, dist, Ay);
   
       POINT* ystart = &Ay[0];
       POINT* pp = ystart;
       POINT* p;

       for(p = mid-1; p >= start; p--)
       { if (square(mid->xdist(*p)) >= dist) break;
          *pp++ = *p;
        }
   
       for(p = mid+1; p <= stop; p++)
       { if (square(mid->xdist(*p)) >= dist) break;
          *pp++ = *p;
        }
     
       int m = pp - ystart;

       if (m > 1) 
       { Ay.sort(cmp_y,0,m-1);
         y_search(ystart,pp,dist,r1,r2);
        }
     }
}



COORD CLOSEST_PAIR(list<POINT>& L, POINT& r1, POINT& r2)
{
  int n = L.length();

  if (n < 2)
    error_handler(1, "CLOSEST PAIR: input has less than two POINTs.");

  array<POINT> Ax(n);  // array of input POINTs sorted by x-coordinates
  array<POINT> Ay(n);  // auxiliary array

  // copy input POINTs to Ax and sort

  POINT* pp = &Ax[0];
  POINT p;
  forall(p,L) *pp++ = p;

  Ax.sort();

  COORD dist = MAXINT;

  cl_pair(&Ax[0],&Ax[n-1],r1,r2,dist,Ay);

  return dist;

}

