corkum@csri.toronto.edu (Brent Thomas Corkum) (10/22/90)
This routine and sample program computes whether a 2D polygon is clockwise
or counter-clockwise. The function takes vertices passed in an array
(x1,y1,x2,y2,...xn,yn where n=#vertices) and the number of vertices.
Of course for a 3D polygon, clockwise or counter-clockwise depends
on which way you look at it.
Brent Corkum
corkum@csri.toronto.edu
This function is also available via anonymous ftp at
boulder.civ.toronto.edu (128.100.14.11)
int the
pub/CCW_OR_CW
directory.
P.S. Feel free to use and distribute this function. Any use or misuse
of this function is the responsibility of the user. If you
use this function, the header must remain in place to acknowledge
my contribution.
/*---------------------------cut here--------------------------------*/
#include <stdio.h>
#include <math.h>
#define ABS(x) ((x)<0.0 ? -(x) : (x))
/*
compile with:
cc <file.c> -lm
*/
main()
{
float pl[100];
int numv,i,flag;
printf("Enter #numvertices: ");
scanf("%d",&numv);
for(i=1;i<=numv;++i){
printf("Enter Vertex %d: ",i);
scanf("%f%f",&pl[2*(i-1)],&pl[2*(i-1)+1]);
}
flag=Isclockwise(pl,numv);
if(flag==1) printf("polygon is clockwise\n");
else if(flag==0) printf("polygon is counter-clockwise\n");
else printf("invalid polygon\n");
}
/*
Function by :
Brent Corkum
Data Visualization Lab
Civil Engineering
University of Toronto
Toronto Ontario Canada
corkum@csri.toronto.edu or
corkum@boulder.civ.toronto.edu (128.100.14.11)
That determines whether a 2d polygons vertices are orderred clockwise
or counter-clockwise. The polygon is passed as an array of vertices
(x1,y1,x2,y2,...,xn,yn where n=numv=#vertices or edges in the closed polygon).
The function returns 1 if the polygon is clockwise and 0 if the polygon
is counterclockwise. The function returns -1 if numv<3.
*/
int Isclockwise(pl,numv)
float *pl;
int numv;
{
int index,pindex,aindex,i,flag;
float maxv,eps,dx,dy,adx,ady,xi,yi,xj,yj,xh,yh,xtest,trend;
if(numv<3)return(-1);
/* calculate an epsilon for floating point comparisons */
eps = 1e-5 * ABS(pl[0]);
if(eps<1e-5)eps=1e-5;
/* filter out last point if equal to first point */
if(pl[0]==pl[2*(numv-1)] && pl[1]==pl[2*(numv-1)+1])--numv;
/* compute the vertex with the minimum x value, if there's more than one
compute the vertex that has the minimum x and y value (min vertex=P) */
index=0;
for(i=1;i<numv;++i){
dx=pl[2*i]-pl[2*index];dy=pl[2*i+1]-pl[2*index+1];
adx = ABS(dx); ady = ABS(dy);
if(adx>eps){
if(dx<0.0)index=i;
}
else{
if(dy<0.0)index=i;
}
}
/* compute the previous vertex, making sure that the previous vertex does
not equal the min vertex (Pp)*/
pindex=index;flag=1;
do{
--pindex;
if(pindex<0)pindex=numv-1;
dx=pl[2*pindex]-pl[2*index];dy=pl[2*pindex+1]-pl[2*index+1];
adx = ABS(dx); ady = ABS(dy);
if(adx>eps || ady>eps)flag=0;
else if(pindex==index)flag=0;
}while(flag);
/* compute the vertex coming after the min vertex, making sure that it
doesn't equal the min vertex (Pa)*/
aindex=index;flag=1;
do{
++aindex;
if(aindex==numv)aindex=0;
dx=pl[2*aindex]-pl[2*index];dy=pl[2*aindex+1]-pl[2*index+1];
adx = ABS(dx); ady = ABS(dy);
if(adx>eps || ady>eps)flag=0;
else if(aindex==index)flag=0;
}while(flag);
/* compute whether angle made by Pp -> P ->Pa is convex or concave,
the angle is measured counterclockwise from Pp->P */
if(aindex != index && pindex != index){
xh = pl[2*pindex] ; yh = pl[2*pindex+1];
xi = pl[2*index] ; yi = pl[2*index+1];
xj = pl[2*aindex] ; yj = pl[2*aindex+1];
/* compute azimuth of Pp->P */
dx = xi-xh ; dy = yi-yh;
if(dx==0.0 && dy==0.0) trend = 0.0;
else{
trend = atan2(dx,dy);
if(trend < 0.0) trend += 6.283185308; /* += 2*PI */
}
/* rotate x value of Pp->P to positive y axis(x=0) and x of
P-Pa by same angle */
dx = xj-xi ; dy = yj-yi;
xtest = dx*cos(trend) - dy*sin(trend);
/* if xtest is >= 0.0 then you know that the angle between Pp->P and
P->Pa is CONVEX and therefore the polygon is clockwise */
if(xtest>=0.0)return(1);
else return(0);
}
return(-1);
}td@alice.att.com (Tom Duff) (10/23/90)
Brent Corkum of Toronto posted a 165 line solution to somebody's
request for a polygon orientation (cw/ccw) test. A much simpler
method is to compute the polygon's area by
summing trapezoid areas. If the sum is positive, the vertices
are passed clockwise, otherwise they're passed counter-clockwise.
If a polygon is not simple (i.e. if its boundary intersects
itself) it may be impossible to assign a meaningful orientation.
In that case this method fails, as it ought to. It does not deliver
a diagnostic in that case, as simplicity testing is a rather more
complicated problem than computing area -- increasing the size
of the code by a factor of 20 (just a guess, 50 may be more realistic)
for the sake of an error check would be silly.
Also, another poster (who will remain nameless) submitted a hilarious
proposal for finding the orientation of a polygon in 3-space. In
3 dimensions, a polygon has no orientation: if it appears clockwise viewed
from one side, then it's counter-clockwise from the other.
Here, then, is a simpler orientation test, built to Corkum's
specification. The arguments are a pointer to an array of coordinates
arranged (x0, y0, x1, y1, ..., xn, yn), and a count of vertices. The
return value is 0 is the points are passed clockwise, 1 if counter-clockwise
and -1 if no determination can be made. This works only for simple polygons.
Isclockwise(float *pv, int nv){
float *p1, *p2, *endp=pv+2*nv;
float area=0; /* acually twice the area */
for(p1=pv, p2=endp-1; p1!=endp; p2=p1, p1+=2)
area+=(p2[0]-p1[0])*(p2[1]+p1[1]);
return area<0 ? 1 : area>2 ? 0 : -1;
}td@alice.att.com (Tom Duff) (10/23/90)
My last posting contained a small error. Perhaps this one
does too. The correct text of the orientation-testing routine
is (I think):
Isclockwise(float *pv, int nv){
float *p1, *p2, *endp=pv+2*nv;
float area=0; /* acually twice the area */
for(p1=pv, p2=endp-1; p1!=endp; p2=p1, p1+=2)
area+=(p2[0]-p1[0])*(p2[1]+p1[1]);
return area<0 ? 1 : area>0 ? 0 : -1;
}davidle@microsoft.UUCP (David LEVINE) (10/25/90)
In article <1990Oct22.114039.8799@jarvis.csri.toronto.edu> corkum@csri.toronto.edu (Brent Thomas Corkum) writes: >This routine and sample program computes whether a 2D polygon is clockwise >or counter-clockwise. The function takes vertices passed in an array >(x1,y1,x2,y2,...xn,yn where n=#vertices) and the number of vertices. >Of course for a 3D polygon, clockwise or counter-clockwise depends I believe that it would be much easier to use the simple formula which computes the AREA of the polygon and then look at the resulting sign. For a clockwise polygon, the area is: sum ( ( x(i) - x(i-1) ) * ( ( y(i) + y(i-1) ) ) / 2 Where x(i),y(i) and x(i-1),y(i-1) are adjacent pairs of points going clockwise around the polygon. If the sum comes out negative then the polygon was counter-clockwise. Of course, if you only want the area then don't bother to divide by two. This is much easier to code, although it is possible that for very long, thin or otherwise strange polygons, some special case code might be better. Similarly, this formula is not guarenteed to work with polygons that have crossed edges (I don't know what is defined to be clockwise then, anyway). I posted the derivation of the area formula sometime last year. DLL