home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / P_APP.C < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  19.0 KB  |  605 lines
  1. /* 
  2.  * p_app.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * P(oly_)APP(roximation)
  11.  *
  12.  */
  13.  
  14.  
  15. /* Code fragments for computation of polygonal approximation and convex hull
  16.  * of a given boundary.  To use, first write a routine to supply the boundary,
  17.  * as decribed in struct Bdy.  To compute an approximation, choose an error
  18.  * tolerance (say 1.5 pixels), and call
  19.  *   poly_approx(&bdy,tolerance);
  20.  * which will enrich Bdy bdy with an approximation.  To compute the hull of 
  21.  * the approximation, call:
  22.  *   convex_hull(&bdy);
  23.  *
  24.  * Assumes the availability of a sorting function qsort() defined as follows:
  25.  *   qsort(base, nel, width, compar)
  26.  *   char *base;
  27.  *   int (*compar)();
  28.  *   
  29.  * The first argument is a pointer to the base of the data; the second is
  30.  * the number of elements; the third is the width of an element in bytes;
  31.  * the last is the name of the comparison routine to be called with two 
  32.  * arguments which are pointers to the elements being compared.  It should 
  33.  * be declared as:
  34.  *  compar(a, b)
  35.  *  char *a, *b;
  36.  * 
  37.  * The routine must return an integer less than, equal to, or greater than 0
  38.  * according as the first argument is to be considered less than, equal to,
  39.  * or greater than the second.
  40.  *
  41.  * Also assumes the availability of a square-root function sqrt(), taking a
  42.  * double argument and returing a double value.  This may be a fast 
  43.  * approximate implementation.  It is used rarely.
  44.  *
  45.  * Also assumes memory allocation functions malloc() and free(), and stdio.
  46.  */
  47. /* 
  48.  * A boundary is an ordered list of vertices, assumed to close.  `vn' counts
  49.  * the vertices.  We use 8-connected boundaries, but that is immaterial
  50.  * to the code.  For the convenience of poly_approx(), the 1st vertex must be
  51.  * duplicated at the end (but `vn' doesn't count the redundant vertex).
  52.  * Except for this case, the input boundary has no consecutive duplicate
  53.  * vertices, and horizontal and vertical runs may (or may not, as you wish)
  54.  * be compressed by omission of their interior colinear vertices.  The winding
  55.  * order of the boundary must be specified by the Bdy_ccw bit in ident (our
  56.  * convention is that counter-clockwise boundaries enclose black regions,
  57.  * and clockwise boundaries enclose holes, but this is immaterial here).
  58.  * `per' counts the no. of pixels on the 8-connected boundary.
  59.  * (definition in bdy_struct.h or sgl_stat.h)
  60.  *
  61.  * NOTE:
  62.  * Original code by H. S. Baird (Spring 1987).
  63.  * Substantially revised, speeded up, and improved by R J Elliott (Spring 1988).
  64.  * Further bug fixes and speedups by H. S. Baird (Fall 1989).
  65.  * Thorsten von Eiken suggests:  backup trick might be sharper 
  66.  *  if you backup only as long as the incremental area is falling fast;  
  67.  *  when it knees, stop.
  68.  */
  69.  
  70.  
  71. #include <stdio.h>
  72. #include <math.h>
  73. #include <malloc.h>
  74. #include <stdlib.h>
  75. #include <string.h>
  76. #include <search.h>
  77.  
  78. #include "ip.h"
  79.  
  80. #undef    DEBUG
  81.  
  82. #define Bdy_ccw    0000000001      /* winding order is counter-clockwise */
  83.  
  84. #define BACKUP_F (0.75)         /* break if distance from the last breakpoint turns
  85.                                  * * * around and shortens by more than BACKUP_F *
  86.                                  * * * effective tolerance, compared to the longest so 
  87.                                  * * * far
  88.                                  */
  89.  
  90. #define RAMER_F (3)             /* the number of backup points traversed
  91.                                  * * * after the length of the dist starts decreasing
  92.                                  * * * is RAMER_F * effective tolerance
  93.                                  */
  94.  
  95.  
  96. /*
  97.  * poly_approx()
  98.  *   DESCRIPTION:
  99.  *     (Re)Compute a polygonal (straight-line segment) approximation to a given
  100.  *     2D boundary.
  101.  *   ARGUMENTS:
  102.  *     bdp(Bdy *) pointer polygon boundary struct (see ph.h)
  103.  *     tol(double) tolerance in pixels
  104.  *   RETURN VALUE:
  105.  *     none
  106.  *   COMMENTS:
  107.  *     The breakpoint vertices are chosen from among the given 
  108.  *     vertices.
  109.  *     In vestigial cases an empty approximation may be returned; otherwise, the
  110.  *     approximation will have at least two vertices.  The original boundary need 
  111.  *     not be connected; but the 1st boundary vertex must be repeated after the 
  112.  *     last (simplifying this code).
  113.  *     Uses the fast merging method of Wall & Danielsson (CVGIP 28, pp 220-227,
  114.  *     1984) followed by local Ramer-like corrections.  W&D guarantees that the
  115.  *     ``included area'' between the original boundary and each approximation edge
  116.  *     is less than maxarea.  In an extension to the original algorithm 
  117.  *     description, we also ensure that the approximation never backs up:  
  118.  *     i.e. turns around and shortens the maximum distance seen so far from the 
  119.  *     last breakpoint by more than BACKUP_F times the tolerance.
  120.  *     Also, W&D tends to overshoot.  Each breakpoint is therefore adjusted by
  121.  *     backing up to a local maximum perpendicular distance from the line between
  122.  *     the prior (already-corrected) breakpoint and the next (not yet corrected)
  123.  *     breakpoint.  The interval looked back on stretches from the current 
  124.  *     breakpoint through at most BKUP*maxarea pixels of monotonically 
  125.  *     non-increasing perp. distance.  This is a kind of local Ramer correction.
  126.  *     Finally, the firstbreakpoint may be adjusted  or even deleted if in its
  127.  *     absence no breakpoint would occur.
  128.  *     If there is already an approximation, it is deleted and replaced.
  129.  *     Does not ensure that the approximation has no two consecutive identical
  130.  *     breakpoints (hard to do on the fly).
  131.  *     Absolutely fast implementation:  uses no sqrt() or transcendental fns.
  132.  *     Uses floating-point arithmetic sparingly.
  133.  */
  134.  
  135. void
  136. poly_approx (bdp, tol)
  137.      Bdy *bdp;
  138.      double tol;                /* tolerance, in pixels */
  139.  
  140. {
  141.   int rx, ry, il2, ll2;
  142.   Sp *or, *ip, *lp, *eip;
  143.   double f, a2;
  144.   double efftol;                /* effective tolerance (larger for half-pixel bdys) */
  145.   Sp *maxp;
  146.   int maxd, pdist;
  147.   int backcnt, cnt;             /* limit backup by this many pixels */
  148.  
  149. #ifdef DEBUG
  150.   printf ("\n...arrived in poly_appox()");
  151.   printf ("    tol: %lf", tol);
  152.   printf ("    bdp->an initialized to: %ld\n", bdp->an);
  153. #endif
  154.  
  155.   if (bdp->an > 0 && bdp->ap != NULL) {
  156.     /* free any prior approximation */
  157.     bdp->an = 0;
  158.     free (bdp->ap);
  159.     bdp->err = (float) -1.0;
  160.   };
  161.   if (bdp->vn < 2)
  162.     return;
  163.  
  164.   if ((bdp->ap = (Sp **) malloc ((bdp->vn + 1) * sizeof (Sp *))) == NULL)
  165.     abort1 ("can't alloc Bdy.ap[%d]", bdp->vn + 1,
  166.             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
  167.  
  168.   if (bdp->vn <= 4 || tol == (double) 0) {
  169.     /* don't approximate - make every vertex a breakpoint */
  170.  
  171. #ifdef DEBUG
  172.     printf ("...no approximation: every vertex is a breakpoint\n");
  173.     printf ("    bdp->vn =%ld tol = %lf", bdp->vn, tol);
  174. #endif
  175.  
  176.     for (bdp->an = 0, ip = bdp->v; bdp->an < bdp->vn; bdp->an++) {
  177.       bdp->ap[bdp->an] = ip++;
  178.     }
  179.     bdp->err = (float) 0;
  180.     return;
  181.   };
  182.  
  183.   efftol = tol;
  184.   a2 = efftol * efftol;
  185.   if ((backcnt = (int) (RAMER_F * efftol)) < 1)
  186.     backcnt = 1;
  187.   else if (backcnt > bdp->vn / 2)
  188.     backcnt = bdp->vn / 2;
  189.  
  190.   /* 1st breakpoint is (provisionally) the 1st vertex */
  191.   eip = (ip = bdp->ap[0] = bdp->v) + bdp->vn;
  192.   bdp->an = 1;
  193.   while (ip < eip) {
  194.     or = ip++;
  195.     for (f = 0.0, ll2 = 0, lp = ip++; ip <= eip; lp = ip++) {
  196.       /* next point relative to origin */
  197.       rx = ip->x - or->x;
  198.       ry = ip->y - or->y;
  199.       il2 = rx * rx + ry * ry;
  200.       if (il2 >= ll2)
  201.         ll2 = il2;
  202.       else {                    /* fast, approximate test for backing up */
  203.         double a, b;
  204.         a = (double) (ll2 - il2);
  205.         a = a * a;
  206.         b = (4.0 * BACKUP_F * BACKUP_F) * efftol * efftol * ll2;
  207.         if (a > b) {
  208.           break;
  209.         };
  210.       };
  211.       /* discrepancy area so far */
  212.       f += (double) (rx * (ip->y - lp->y) - ry * (ip->x - lp->x));
  213.       if ((f * f) > (a2 * il2)) {
  214.         break;
  215.       }
  216.     };
  217.     if (ip <= eip) {
  218.       /* backing up or exceeding threshold: a new breakpt */
  219.       bdp->ap[bdp->an++] = --ip;
  220.  
  221.       /* Attempt to adjust prior breakpt bdp->ap[bdp->an-2]
  222.        * by backing it up until its perpendicular distance 
  223.        * from the line between its neighboring breakpoints 
  224.        * is maximum.
  225.        * Only look back a distance bounded by 
  226.        * RAMER_F*efftol^2.
  227.        * Among those with equal perpendicular distance,
  228.        * choose the one farthest along the boundary. */
  229.       if (bdp->an > 2) {
  230.         or = bdp->ap[bdp->an - 3];
  231.         rx = ip->x - or->x;
  232.         ry = ip->y - or->y;
  233.         for (cnt = 0, maxd = 0, maxp = ip = bdp->ap[bdp->an - 2];
  234.              ip != or; ip--) {
  235.           pdist = rx * (ip->y - or->y) - ry * (ip->x - or->x);
  236.           pdist = abs (pdist);
  237.           if (pdist > maxd) {
  238.             maxd = pdist;
  239.             maxp = ip;
  240.             cnt = 0;
  241.           }
  242.           else if (pdist < maxd) {
  243.             if (++cnt >= backcnt)
  244.               break;
  245.           }
  246.           else
  247.             cnt = 0;
  248.         };
  249.         bdp->ap[bdp->an - 2] = maxp;
  250.       };
  251.       ip = bdp->ap[bdp->an - 1];
  252.     };
  253.   };
  254.   /* Have arrived back at 1st vertex, the 1st breakpoint */
  255.  
  256.   /* Try deleting the 1st breakpoint.  Look past it to test
  257.    * whether a breakpoint is called for:  if not, delete it. */
  258.   if (bdp->an >= 2) {
  259.     or = bdp->ap[bdp->an - 1];
  260.     if ((ip = or + 1) >= bdp->v + bdp->vn)
  261.       ip = bdp->v;
  262.     if ((eip = bdp->ap[1] + 1) >= bdp->v + bdp->vn)
  263.       eip = bdp->v;
  264.     f = 0.0;
  265.     ll2 = 0;
  266.     lp = ip++;
  267.     while (ip != eip) {
  268.       /* next point relative to origin */
  269.       rx = ip->x - or->x;
  270.       ry = ip->y - or->y;
  271.       il2 = rx * rx + ry * ry;
  272.       if (il2 >= ll2)
  273.         ll2 = il2;
  274.       else {                    /* fast, approximate test for backing up */
  275.         double a, b;
  276.         a = (double) (ll2 - il2);
  277.         a = a * a;
  278.         b = (4.0 * BACKUP_F * BACKUP_F) * efftol * efftol * ll2;
  279.         if (a > b)
  280.           break;
  281.       };
  282.       /* discrepancy area so far */
  283.       f += (double) (rx * (ip->y - lp->y) - ry * (ip->x - lp->x));
  284.       if ((f * f) > (a2 * il2))
  285.         break;
  286.       if (ip >= (bdp->v + bdp->vn)) {
  287.         lp = ip;
  288.         ip = bdp->v;
  289.       }
  290.       else
  291.         lp = ip++;
  292.     };
  293.     if (ip == eip) {
  294.       /*  No need for a breakpoint: delete 1st breakpoint */
  295.       bdp->ap[0] = bdp->ap[bdp->an - 1];
  296.       bdp->an -= 1;
  297.  
  298. #ifdef DEBUG
  299.       printf ("...upon closure: no need for a breakpt\n");
  300.       printf ("  after deleting 1st breakpt: an = %ld\n",
  301.               bdp->an);
  302. #endif
  303.  
  304.       if (bdp->an > 2) {        /* Try adjusting new 1st breakpoint */
  305.         or = bdp->ap[bdp->an - 1];
  306.         rx = bdp->ap[1]->x - or->x;
  307.         ry = bdp->ap[1]->y - or->y;
  308.         for (cnt = 0, maxd = 0, maxp = ip = bdp->ap[0]; ip != or; ip--) {
  309.           if (ip == bdp->v)
  310.             ip = bdp->v + bdp->vn;
  311.           pdist = rx * (ip->y - or->y) - ry * (ip->x - or->x);
  312.           pdist = abs (pdist);
  313.           if (pdist > maxd) {
  314.             maxd = pdist;
  315.             maxp = ip;
  316.             cnt = 0;
  317.           }
  318.           else if (pdist < maxd) {
  319.             if (++cnt >= backcnt)
  320.               break;
  321.           }
  322.           else
  323.             cnt = 0;
  324.         };
  325.         bdp->ap[0] = maxp;
  326.       };
  327.     }
  328.     else {                      /* Try adjusting last breakpoint */
  329.       or = bdp->ap[bdp->an - 2];
  330.       rx = bdp->ap[0]->x - or->x;
  331.       ry = bdp->ap[0]->y - or->y;
  332.       for (cnt = 0, maxd = 0, maxp = ip = bdp->ap[bdp->an - 1]; ip != or; ip--) {
  333.         if (ip == bdp->v)
  334.           ip = bdp->v + bdp->vn;
  335.         pdist = rx * (ip->y - or->y) - ry * (ip->x - or->x);
  336.         pdist = abs (pdist);
  337.         if (pdist > maxd) {
  338.           maxd = pdist;
  339.           maxp = ip;
  340.           cnt = 0;
  341.         }
  342.         else if (pdist < maxd) {
  343.           if (++cnt >= backcnt)
  344.             break;
  345.         }
  346.         else
  347.           cnt = 0;
  348.       };
  349.       bdp->ap[bdp->an - 1] = maxp;
  350.  
  351.       /* Try adjusting 1st breakpoint */
  352.       or = bdp->ap[bdp->an - 1];
  353.       rx = bdp->ap[1]->x - or->x;
  354.       ry = bdp->ap[1]->y - or->y;
  355.       for (cnt = 0, maxd = 0, maxp = ip = bdp->ap[0]; ip != or; ip--) {
  356.         if (ip == bdp->v)
  357.           ip = bdp->v + bdp->vn;
  358.         pdist = rx * (ip->y - or->y) - ry * (ip->x - or->x);
  359.         pdist = abs (pdist);
  360.         if (pdist > maxd) {
  361.           maxd = pdist;
  362.           maxp = ip;
  363.           cnt = 0;
  364.         }
  365.         else if (pdist < maxd) {
  366.           if (++cnt >= backcnt)
  367.             break;
  368.         }
  369.         else
  370.           cnt = 0;
  371.       };
  372.       bdp->ap[0] = maxp;
  373.     };
  374.   };
  375.  
  376.   /* Force approximation to have >= 2 breakpoints; also, check for and
  377.    * correct bug where the 1st or last breakpoint is off the end of v[] */
  378.   if (bdp->an >= 2) {
  379.     if (bdp->ap[0] == (bdp->v + bdp->vn))
  380.       bdp->ap[0] = bdp->v;
  381.     else if (bdp->an > 1) {
  382.       if (bdp->ap[bdp->an - 1] == (bdp->v + bdp->vn))
  383.         bdp->ap[bdp->an - 1] = bdp->v;
  384.     };
  385.     if ((bdp->ap = (Sp **) realloc (bdp->ap, bdp->an * sizeof (Sp *))) == NULL)
  386.       abort1 ("can't realloc Bdy.ap[%d]", bdp->an,
  387.               0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
  388.     bdp->err = (float) tol;
  389.   }
  390.   else {
  391.     if (bdp->ap != NULL) {
  392.       free (bdp->ap);
  393.       bdp->ap = NULL;
  394.     };
  395.     bdp->an = 0;
  396.   };
  397.  
  398. #ifdef DEBUG
  399.   printf ("...reached end of poly_approx():bdp->an = %ld", bdp->an);
  400. #endif
  401. }
  402.  
  403.  
  404. #define old_BACKUP_F (0.5)      /* break if distance from the last breakpoint turns
  405.                                  * * * around and shortens by more than BACKUP_F *
  406.                                  * * * tolerance compared to the longest seen */
  407. #define old_RAMER_F (5)         /* vary the number of backup points traversed
  408.                                  * * * after the length of the dist starts decreasing
  409.                                  * * * with RAMER_F * tolerance */
  410.  
  411.  
  412. void
  413. err (s, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12)
  414.      char *s;
  415.      int a1, a2, a3, a5, a6, a7, a8, a9, a10, a11, a12;
  416.      double a4;
  417.  
  418. {
  419.   char m[100];
  420.   strcpy (m, "polyapprox: ");
  421.   strcat (m, s);
  422.   strcat (m, "\n");
  423.   fprintf (stderr, m, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12);
  424. }
  425.  
  426. void
  427. abort1 (s, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12)
  428.      char *s;
  429.      double a4;
  430. {
  431.   static char m[100];
  432.   strcpy (m, s);
  433.   strcat (m, " - abort");
  434.   err (m, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12);
  435.   exit (1);
  436. }
  437.  
  438.  
  439.  
  440. /* 
  441.  *Globals used by convex_hull functions 
  442.  */
  443.  
  444. Sp HULL_c, HULL_o;              /* center of bdy bx and origin vector */
  445.  
  446. /* macros for use within winding order function */
  447. /* dot product */
  448. #define dot(a,b) ( (a).x*(b).x+(a).y*(b).y )
  449.  
  450. /* perp product */
  451. #define perp(a,b) (  -((a).y)*(b).x+(a).x*(b).y )
  452.  
  453. /* 
  454.  * return the quadrant no: 0,1,2,3 of pt `a' wrt to the X-Y coordinate system
  455.  * defined by aligning the positive X-axis along vector `o' (`a' is a vector
  456.  * wrt to `o') 
  457.  */
  458. #define sgn(a) ( ((a) > 0) ? +1 : -1)
  459. #define quad(a,o) (  ((dot((o),(a)))>=0) ?  ( (perp((o),(a))>=0) ? 0 : 3)  :  ( (perp((o),(a))>=0) ? 1 : 2) )
  460.  
  461. int
  462. winding_ccw (appp, bppp)
  463.      Sp ***appp, ***bppp;
  464. {
  465.   Sp a, b;
  466.   int qa, qb;
  467.  
  468. /* 
  469.  * return -1,0,1 as pt `a' is less than, equal to, or greater than pt `b' in
  470.  * counterclockwise winding order about an (implicit) center `HULL_c', where
  471.  * the vector `o' defines the `origin' direction. 
  472.  */
  473. #define wind(a,b,o) (  ((qa=quad((a),(o))) == (qb=quad((b),(o)))) ?  sgn(perp((a),(b)))  :  sgn(qb-qa) )
  474.  
  475.   a.x = (**appp)->x - HULL_c.x;
  476.   a.y = (**appp)->y - HULL_c.y;
  477.   b.x = (**bppp)->x - HULL_c.x;
  478.   b.y = (**bppp)->y - HULL_c.y;
  479.   return (wind (a, b, HULL_o));
  480. }
  481.  
  482.  
  483. int
  484. winding_cw (appp, bppp)
  485.      Sp ***appp, ***bppp;
  486. {
  487.   return (-winding_ccw (appp, bppp));
  488. }
  489.  
  490. /*
  491.  * convex_hull()
  492.  *   DESCRIPTION:
  493.  *     (Re)Compute the convex non-collinear hull of the polygonal approximation
  494.  *     of a given boundary.  Uses a sort in winding order about the center,
  495.  *     followed by a Graham scan.
  496.  *     The hull vertices are chosen from the vertices of the approximation.
  497.  *     Assumes that the first approximation point is on the hull.
  498.  *     Does not delete the first point if it turns out that it is on the line from
  499.  *     last to second (in short, if the first point is not on the non-collinear hull);
  500.  *   ARGUMENTS:
  501.  *     bdp(Bdy *) pointer polygon boundary struct (see ph.h)
  502.  *   RETURN VALUE:
  503.  *     none
  504.  */
  505.  
  506. void
  507. convex_hull (bdp)
  508.      Bdy *bdp;
  509. {
  510.   int pi;
  511.   Sp **pp, ***ppp, ***top, ***next;
  512.   static Sp or =
  513.   {0, 0};
  514.  
  515.   if (bdp->hn > 0 && bdp->hpp != NULL)
  516.     free (bdp->hpp);
  517.   if (bdp->an < 2) {
  518.     bdp->hn = 0;
  519.     bdp->hpp = NULL;
  520.     return;
  521.   };
  522.  
  523.   if ((bdp->hpp = (Sp ***) malloc ((int) (bdp->an * sizeof (Sp **)))) == NULL)
  524.     abort1 ("can't alloc bdp->hpp[%d]", (int) bdp->an,
  525.             0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
  526.  
  527.   /* start with the assumption that every vertex in the polygonal
  528.    * approximation is on the hull */
  529.   /* also, find "center" of approximation:  */
  530.   /* need some point guaranteed to fall inside the hull */
  531.   HULL_c.x = 0;
  532.   HULL_c.y = 0;
  533.   for (bdp->hn = 0, pp = bdp->ap; bdp->hn < (int) bdp->an; bdp->hn++) {
  534.     HULL_c.x += (*pp)->x;
  535.     HULL_c.y += (*pp)->y;
  536.     bdp->hpp[bdp->hn] = pp++;
  537.   };
  538.   if (bdp->an < 4)
  539.     return;                     /* too small */
  540.   HULL_c.x /= (short) bdp->an;
  541.   HULL_c.y /= (short) bdp->an;
  542.  
  543. #ifdef DEBUG
  544.   printf ("\n HULL_center = (%d,%d)\n", HULL_c.x, HULL_c.y);
  545. #endif
  546.  
  547. #define init_scan {top=bdp->hpp-1; next=top-1; bdp->hn=0;}
  548. #define push(pp) {next=top++; *top=(pp); bdp->hn++;}
  549. #define pop {top=next--; bdp->hn--;}
  550.  
  551. /* 
  552.  * interior(a,b,c) - T iff point *a is ``interior'' of the line from *b to *c, that
  553.  * is, in the closed half-plane to the left of b-->c if the winding
  554.  * order is counterclockwise; else in the right half-plane.  The use of closed
  555.  * half-planes implies that only non-collinear hull points are selected. 
  556.  */
  557. #define leftof(a,b,c) (  (bdp->ident&Bdy_ccw) ?  ((-((c)->y-(b)->y))*((a)->x-(b)->x)+((c)->x-(b)->x)*((a)->y-(b)->y)) <= 0  :  ((-((c)->y-(b)->y))*((a)->x-(b)->x)+((c)->x-(b)->x)*((a)->y-(b)->y)) >= 0 )
  558.  
  559.   /* define the ``origin'' vector from `c' to the 1st point */
  560.   HULL_o.x = (**bdp->hpp)->x - HULL_c.x;
  561.   HULL_o.y = (**bdp->hpp)->y - HULL_c.y;
  562.   if (HULL_o.x == 0 && HULL_o.y == 0)
  563.     return;                     /* too small */
  564.  
  565.   /* sort approximation pts about center in winding order (1st as 0) */
  566.   if (bdp->ident & Bdy_ccw)
  567. #if defined(LINUX)
  568.     qsort (bdp->hpp + 1, (size_t) bdp->hn - 1, sizeof (Sp ***), (__compar_fn_t) winding_ccw);
  569.   else
  570.     qsort (bdp->hpp + 1, (size_t) bdp->hn - 1, sizeof (Sp ***), (__compar_fn_t) winding_cw);
  571. #else
  572.     qsort (bdp->hpp + 1, (size_t) bdp->hn - 1, sizeof (Sp ***), winding_ccw);
  573.   else
  574.     qsort (bdp->hpp + 1, (size_t) bdp->hn - 1, sizeof (Sp ***), winding_cw);
  575. #endif
  576.  
  577.   /* Graham scan */
  578.   ppp = bdp->hpp;
  579.   init_scan;
  580.   push (*(ppp++));
  581.   push (*(ppp++));
  582.   for (pi = 2; pi < bdp->an; pi++, ppp++) {
  583.     while (bdp->hn > 1 && leftof (**top, **next, **ppp))
  584.       pop;
  585.     push (*ppp);
  586.   };
  587.   while (bdp->hn > 1 && leftof (**top, **next, **(bdp->hpp)))
  588.     pop;
  589.   if (bdp->hn < bdp->an) {
  590.     if ((bdp->hpp = (Sp ***) realloc (bdp->hpp, bdp->hn * sizeof (Sp **))) == NULL)
  591.       abort1 ("can't realloc bdp->hpp[%d]", bdp->hn,
  592.               0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
  593.   };
  594.  
  595. }
  596.  
  597. /*
  598.  * report coordinates of HULL center 
  599.  */
  600. Sp *
  601. reprt_hull_center ()
  602. {
  603.   return (&HULL_c);
  604. }
  605.