home *** CD-ROM | disk | FTP | other *** search
/ Resource Library: Graphics / graphics-16000.iso / general / raytrace / renaisnc / delaunay.lha / Delaunay / .‾GrahamHull.c next >
C/C++ Source or Header  |  1992-01-04  |  3KB  |  127 lines
  1. /*
  2.  File: GrahamHull.c
  3.  Authors: K.R. Sloan
  4.  Last Modified: 5 January 1987
  5.  Purpose: calculates the convex hull of a set of <x,y> points,
  6.           using the Graham Scan method.
  7. */
  8. #include <malloc.h>
  9. #include "DCEL.h"
  10. #include "GrahamHull.h"
  11.  
  12. typedef struct HeapKeyRecord HeapKey;
  13. struct HeapKeyRecord
  14.         {  
  15.          REAL A;
  16.          int P;
  17.         };
  18.  
  19. static void Heapify(V, Root, Last)
  20.  HeapKey V[];
  21.  int Root, Last;
  22.  {
  23.   int r,s;
  24.   HeapKey t;
  25.  
  26.   t = V[Root]; r = Root;
  27.   for (s=r;r==s;)
  28.    {
  29.     s = r+r;
  30.     if (s <  Last) if (V[s].A < V[s+1].A) s = s+1;
  31.     if (s <= Last) if (t.A    < V[s].A  ) { V[r] = V[s]; r = s; }
  32.    }
  33.  
  34.   V[r] = t;
  35.  }
  36.  
  37. void HeapSort(V, LastElement)
  38.  HeapKey V[];
  39.  int LastElement;
  40.  {
  41.   int i;
  42.   HeapKey t;
  43.  
  44.   for (i = LastElement / 2; i > 1; i--) Heapify(V,i,LastElement);
  45.   for (i = LastElement; i > 1; i--)
  46.     { /* 2..i form a Heap */
  47.      Heapify(V, 1,i);
  48.      t = V[1]; V[1] = V[i]; V[i] = t;
  49.     }
  50.   }
  51.  
  52. PointSetPointer GrahamHull(Points, LastPoint)
  53.  Point Points[];
  54.  int LastPoint;
  55.  {
  56.   HeapKey *Heap;
  57.   int i;
  58.   PointSetPointer Hull, H, X, FirstHullPoint;
  59.   REAL Ax, Ay, Bx, By, Cx, Cy, Ox, Oy;
  60.  
  61.   Heap = (HeapKey *) malloc( (LastPoint+1) * sizeof(HeapKey) );
  62.  
  63.   Ax = Points[0].x; Ay = Points[0].y;
  64.   Bx = Points[1].x; By = Points[1].y;
  65.   for (i=2;;i++)
  66.    {
  67.     Cx  =  Points[i].x; Cy  =  Points[i].y;
  68.     if (((Ax-Bx)*(Ay-Cy) - (Ay-By)*(Ax-Cx)) != 0.0) break; /* not colinear */
  69.    }
  70.  
  71.   Ox  =  (Ax+Bx+Cx) / 3.0; Oy  =  (Ay+By+Cy) / 3.0;   /* an internal point */
  72.  
  73.   /* Points[LastPoint] is reserved for just such emergencies */
  74.   Points[LastPoint].x  =  Ox; Points[LastPoint].y  =  Oy;
  75.   for (i=1; i<=LastPoint; i++)
  76.    {
  77.     Heap[i].P  =  i-1; /* off by one to help HeapSort */
  78.     Heap[i].A  =  FowlerAngle(LastPoint, i-1);
  79.    }
  80.  
  81.   HeapSort(Heap, LastPoint);
  82.  
  83.   Hull = NewPoint();
  84.   Hull->P  =  Heap[1].P;
  85.   Hull->Previous  =  Hull; Hull->Next  =  Hull; /* circular */
  86.  
  87.   for (i=2; i<=LastPoint; i++)
  88.    {
  89.     H = NewPoint();
  90.     H->P  =  Heap[i].P; H->Previous  =  Hull->Previous; H->Next  =  Hull;
  91.     Hull->Previous->Next  =  H; Hull->Previous  =  H;
  92.    }  
  93.   /*
  94.    All points are now in a doubly-linked, circular list, sorted
  95.    by angle about O
  96.   */
  97.   H = Hull;
  98.   FirstHullPoint = Hull;
  99.   for(;;)
  100.    {
  101.     H  =  H->Next;
  102.     if (Points[H->P].x < Points[FirstHullPoint->P].x) FirstHullPoint  =  H;
  103.     if (H == Hull) break;
  104.    }
  105.  
  106.  
  107.  /*
  108.   FirstHullPoint is on the hull - now find the others
  109.  */  
  110.   for(H=FirstHullPoint; H->Next != FirstHullPoint; )
  111.    {
  112.     if (LeftTurn(H->P, H->Next->P, H->Next->Next->P))
  113.      H = H->Next;
  114.     else 
  115.      { /* next point is internal */
  116.       X  =  H->Next;
  117.       H->Next  =  X->Next;
  118.       X->Next->Previous  =  H;
  119.       DisposePoint(X);
  120.       if (H != FirstHullPoint)
  121.        H  =  H->Previous; /* back up, maybe */
  122.      }  /* next point is internal */
  123.    }
  124.   free(Heap);
  125.   return (FirstHullPoint);
  126.  }
  127.