/* Richard A. DeVenezia
 * www.devenezia.com
 * 2/3/04
 *
 * Problem based on material originally posted to SAS-L March 27, 2003
 * Re: Sorting polygon boundary segments
 *
 * Fake data generator updated 2/3/04 to create polygons without overlapping edges.
 */

/* return values */
%global DONT_INTERSECT DO_INTERSECT PARALLEL;
%let DONT_INTERSECT = 0;
%let DO_INTERSECT   = 1;
%let PARALLEL       = 2;

/*
 * A macro that generates DATA Step code that determines the
 * intersection point of two line segments
 */

%macro segment_intersection (x1,y1,x2,y2,x3,y3,x4,y4,xvar,yvar,rcvar);
  /* Adapted from http://www.graphicsgems.org/gemsiii/insectc.c
   * (Antonio 1992) Antonio, Franklin,
   *  "Faster Line Segment Intersection,"
   *  Graphics Gems III (David Kirk, ed.),
   *  Academic Press, pp. 199-202, 1992.
   */

  /* Faster Line Segment Intersection   */
  /* Franklin Antonio                   */

  %local EndBlock;
  %let EndBlock = _%substr(%sysfunc(ranuni(0),9.7),3);

  &xvar = .;
  &yvar = .;
  &rcvar = .;

  length Ax Bx Cx Ay By Cy d e f 8;
  length x1lo x1hi y1lo y1hi 8;

  **********************************;

  Ax = &x2 - &x1;
  Bx = &x3 - &x4;

  /* X bound box test*/
  if ( Ax < 0 ) then do;
    x1lo = &x2;
    x1hi = &x1;
  end;
  else do;
    x1hi = &x2;
    x1lo = &x1;
  end;

  if ( Bx > 0 ) then do;
    if (x1hi < &x4  OR  &x3 < x1lo) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;
  else do;
    if (x1hi < &x3  OR  &x4 < x1lo) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;

  **********************************;

  Ay = &y2 - &y1;
  By = &y3 - &y4;

  /* Y bound box test*/
  if ( Ay < 0 ) then do;
    y1lo = &y2;
    y1hi = &y1;
  end;
  else do;
    y1hi = &y2;
    y1lo = &y1;
  end;

  if ( By > 0 ) then do;
    if (y1hi < &y4  OR  &y3 < y1lo) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;
  else do;
    if (y1hi < &y3  OR  &y4 < y1lo) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;

  **********************************;

  Cx = &x1 - &x3;
  Cy = &y1 - &y3;

  d = By*Cx - Bx*Cy;          /* alpha numerator*/
  f = Ay*Bx - Ax*By;          /* both denominator*/

  /* alpha tests*/
  if (f>0) then do;
    if (d<0  OR  d>f) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;
  else do;
    if(d>0  OR  d<f) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;

  e = Ax*Cy - Ay*Cx;          /* beta numerator*/

  /* beta tests*/
  if (f>0) then do;
    if(e<0  OR  e>f) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;
  else do;
    if(e>0  OR  e<f) then do;
      &rcvar = &DONT_INTERSECT;
      goto &EndBlock;
    end;
  end;

  /*compute intersection coordinates*/

  if (f=0) then do;
    &rcvar = &PARALLEL;
    goto &EndBlock;
  end;

  &xvar = &x1 + d*Ax;       /* intersection x */
  &yvar = &y1 + d*Ay;       /* intersection y */

  &rcvar = &DO_INTERSECT;

&EndBlock.:

%mend;


/*
 * Generate some polygon data such that no segment of a polygon
 * crosses any other of the same polygon
 */

options mprint mtrace;

%let maxvertices=12;
%let npolys=6;

data segs (keep=polyid segid x1 y1 x2 y2);

  format polyid segid 4.;


  array x[&maxvertices] _temporary_;
  array y[&maxvertices] _temporary_;

  retain seed 2;

  do polyid = 1 to &npolys;
    do i = 1 to hbound(x); x[i]=.; y[i]=.; end;

    x[1] = ranuni(seed);
    y[1] = ranuni(seed);

    nvertices = 3 + int ((&maxvertices-2)*ranuni(seed));

    do i = 2 to nvertices;

      ntries = 0;
      good = 0;
      do until (good);

        x[i] = ranuni(seed);
        y[i] = ranuni(seed);

        tryagain = 0;

        do j = 1 to i-2 while (not tryagain);
          %segment_intersection (x[j],y[j], x[j+1],y[j+1],  x[i-1],y[i-1],x[i],y[i], xi,yi,rc)
          tryagain = (rc = &DO_INTERSECT);
        end;

        if (i = nvertices) and not tryagain then do;

          do j = 1 to i-1 while (not tryagain);
            %segment_intersection (x[j],y[j], x[j+1],y[j+1],  x[i],y[i],x[1],y[1], xi,yi,rc)

            tryagain = (rc = &DO_INTERSECT);

            if j=1 and tryagain and xi=x[1] and yi=y[1] then
              tryagain = 0;
            else
            if j=i-1 and tryagain and xi=x[i] and yi=y[i] then
              tryagain = 0;
          end;
        end;

        ntries + 1;
        good = not tryagain;
      end;

    end;

    segid = 0;
    do i = 2 to nvertices;
      segid + 1;
      x1 = x[i-1];
      y1 = y[i-1];
      x2 = x[i];
      y2 = y[i];
      output;
    end;
    segid+1;
    x1 = x2;
    y1 = y2;
    x2 = x[1];
    y2 = y[1];
    output;
    * a dummy seg, so plot shows a closed figure;
    segid=.;
    x1 = x[1];
    y1 = y[1];
    output;
  end;
run;

/*
 * Draw the polygons
 */

ods listing;

symbol1 i=join;
axis1 order=(0 to 1) label=none minor=none value=none major=none;
proc gplot data=segs;
  plot y1*x1=polyid / haxis=axis1 vaxis=axis1 noaxis;
run;
quit;

/*
 * disorder the vertex information (from the perspective of polyid, segid)
 */

proc sql;
  create table rawsegs as
  select * from segs
  where segid
  order by ranuni(1234);
quit;

/*
 * output the segments as a series of points.
 * original polyid and segid is retained so results of algorithm can be
 * compared to original segment configuration
 * the P1 point of the segment is stored with myseginfo=1
 * the P2 point of the segment is stored with myseginfo=2
 */

data labelled_vertices (keep=polyid segid mysegid myseginfo x y);
  set rawsegs;
  mysegid+1;
  myseginfo=1;
  x = x1;
  y = y1;
  OUTPUT;
  myseginfo=2;
  x = x2;
  y = y2;
  OUTPUT;
  format myseg: 4.;
run;

/*
 * sort by x,y.
 * Presuming all vertices of all polygons are unique
 * then one can infer two consecutive rows with the same x and y
 * are endpoints of two different segments in the same polygon
 */

proc sql;
  create table ordered_vertices as
  select * from labelled_vertices
  order by x, y, myseginfo;
quit;

/*
 * process the ordered vertices
 *
 * Consider segment B (Bx1,By1 -> Bx2,By2) connected
 *       to segment A (Ax1,Ay1 -> Ax2,Ay2).
 *
 * From the data assumptions we know that Ax1,Ay1 = Bx2,By2
 *
 * If we have some x,y and determine it is Ax1,Ay1 then
 * we can determine segment B by locating Bx2,By2.
 */

%let nsegments = %eval(&sqlobs/2);

data polygons;
  array key [&nsegments] _temporary_;

  if _n_ = 1 then do;

    * populate an inverse mapping.;

    do while (not eov);
      set ordered_vertices end=eov;
      i+1;
      if myseginfo = 2 then continue;
      key[mysegid] = i;
    end;
  end;

  remain = &nsegments;

*put remain=;

  * Note: The nature of the construction of the data ensures
  * the first obs (nextpolyobs) will be a P1 (e.g. myseginfo=1);

  thepolyid = 0;
  nextpolyobs = 1;

  do until (remain=0);

    vertexobs = nextpolyobs;
    thepolyid + 1;
    thesegid = 0;
    done = 0;

    * http://astronomy.swin.edu.au/~pbourke/geometry/polyarea/ ;
    * http://astronomy.swin.edu.au/~pbourke/geometry/clockwise/ ;
    area = 0;

    do until (vertexobs<0);

      set ordered_vertices point=vertexobs;
      done = key [mysegid+0] < 0;

      if not done then do;
        thesegid + 1;

        x1 = x;
        y1 = y;

        * mark the P1 point of segment B as used;
        key [mysegid+0] = -key [mysegid+0];

        * the next obs in the data set is the P2 point of segment B ;
        vertexobs + 1;

        * read obs to get mysegid of segment B;
        set ordered_vertices(keep=mysegid) point=vertexobs;

        * inverse map mysegid to get P1 of segment B;

        vertexobs = abs(key [mysegid+0]) ;

        set ordered_vertices(keep=x y rename=(x=x2 y=y2)) point=vertexobs;
        OUTPUT;

        area + (x1*y2 - x2*y1);

        vertexobs = key [mysegid+0] ;

        * If the P1 vertex had been processed earlier, then vertexobs will
        * be <0 (marked as used) and cause the polygon do until loop to end;

        remain +(-1);
      end;
    end;

    area = area/2;
    if area > 0
      then orientation = 'CCW';
      else orientation = 'CW ';

    * dummy record so a closed figure can be graphed;
    thesegid=.;
    mysegid=.;
    segid=.;
    x1=x2;
    y1=y2;
    x2=.;
    y2=.;
    output;

    orientation = '';

    * search for start of next poly;
    nextpolyobs = 1;
    if remain then
    do until (key [mysegid+0] > 0);
      nextpolyobs+1;
      set ordered_vertices point=nextpolyobs;
    end;

  end;

  STOP;

  keep polyid segid thepolyid thesegid mysegid x1 y1 x2 y2 orientation;
  format the: 4.;
run;

/*
 * plot the results of the polygon discovery
 */

data anno;
  set polygons (rename=(x1=x y1=y));
  retain xsys ysys '2';
  length function $8 text style $40;
  retain style "'Arial'";
  function='move'; output;

  size=1.8;
  position = '+';
  function='label'; text=put(thesegid,2.-L); output;
  function='label'; size=2; position='>'; text = ' ' || orientation; output;


  keep xsys ysys x y function text style size position;
run;

proc gplot data=polygons anno=anno;
  plot y1*x1=thepolyid / haxis=axis1 vaxis=axis1 noaxis;
run;
quit;

options nocenter linesize=100 pagesize=100;
proc print data=polygons;
  by thepolyid;
  id thepolyid;
run;