Difference between revisions of "Point-in-polygon"

From aemwiki
Jump to: navigation, search
(New page: ; Language: C++ ; Purpose: Determines whether or not a point is in a polygon. Jumps in complex arguments are preserved. ; Author: James R. Craig ; GW Gourmet Functions Used: logZp1oZm1...)
 
Line 35: Line 35:
 
  }
 
  }
 
  }
 
  }
  //if not close, use faster algorithm
+
  //if not close, use faster ray-tracing algorithm
 
  //adapted from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
 
  //adapted from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
 
  if (!close)
 
  if (!close)
Line 51: Line 51:
 
  return in;
 
  return in;
 
  }
 
  }
  //if close, use complex algortithm
+
  //if close, use complex winding algorithm
 +
        //this is consistent with the jump calculation for line elements
 
  //sum ln (Z-1)/(Z+1) should be zero if outside
 
  //sum ln (Z-1)/(Z+1) should be zero if outside
 
  else
 
  else

Revision as of 13:08, 7 December 2007

Language
C++
Purpose
Determines whether or not a point is in a polygon. Jumps in complex arguments are preserved.
Author
James R. Craig
GW Gourmet Functions Used
logZp1oZm1
#include <math.h>
#include <complex>
#include <iostream>
typedef complex<double> cmplex;
const double CLOSE_TO_LINE=1.02;// ellipse radius around line which determines closeness
const double REALSMALL=1e-10; 	//almost zero
/**************************************************************************
IsInsidePolygon:
Returns if a point z is in polygon defined by N vertices stored in 
array ze[]. The last polygon point is NOT repeated.
Adapted from Paul Bourke <http://astronomy.swin.edu.au/pbourke/>
Uses function logZm1oZp1(Z), available at www.analyticelements.org 
"Groundwater Gourmet"
Author: James R. Craig
-------------------------------------------------------------------------*/
bool IsInsidePolygon(const cmplex &z, const cmplex *zv, const int N)
{
	double y(z.imag());
	double x(z.real()); 
	int  i;
	bool close(false); 
	if (abs(zv[0]-zv[N])>REALSMALL){
		cout<<"PolyIsInside : Bad polygon"<<endl;
		exit(0);
	} 
	//check if point is close to boundary of polygon
	for (i=0; i<N; i++){
		if (((abs(z-zv[i])+abs(z-zv[i+1]))/(abs(zv[i+1]-zv[i])))<CLOSE_TO_LINE){
			close=true;break;
		}
	}
	//if not close, use faster ray-tracing algorithm
	//adapted from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
	if (!close)
	{
		bool	 in(false);
		int 	 j;
		for (i=0, j=N-1; i<N; j=i++) {
			if ((((zv[i].imag()<=y) && (y<zv[j].imag())) ||
					 ((zv[j].imag()<=y) && (y<zv[i].imag()))) &&
				 (x<(zv[j].real()-zv[i].real())*(y-zv[i].imag())/
				 (zv[j].imag()-zv[i].imag())+zv[i].real())){
				in=!in;
			}
		}
		return in;
	}
	//if close, use complex winding algorithm
       //this is consistent with the jump calculation for line elements
	//sum ln (Z-1)/(Z+1) should be zero if outside
	else
	{
		cmplex Z;
		double sum(0);
		for (int i=0;i<N; i++)
		{
			Z=(z-0.5*(zv[i]+zv[i+1]))/(0.5*(zv[i+1]-zv[i]));
			sum+=logZm1oZp1(Z).imag();
		}
		return (fabs(sum)<REALSMALL ? false : true);
	}
}