/* 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 d0) then do; if(e<0 OR e>f) then do; &rcvar = &DONT_INTERSECT; goto &EndBlock; end; end; else do; if(e>0 OR e 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;