Difference between revisions of "Point-in-polygon"

Latest revision as of 04:05, 15 February 2011

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 <cmath>
#include <complex>
#include <iostream>
typedef std::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)
{
	bool close = false; 
	if (std::abs(zv[0]-zv[N])>REALSMALL){
		std::cout<<"PolyIsInside : Bad polygon"<<std::endl;
		exit(0);
	} 
	//check if point is close to boundary of polygon
	for (int i=0; i<N; i++){
		if (((std::abs(z-zv[i])+std::abs(z-zv[i+1]))/std::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>
	double y = z.imag();
	double x = z.real(); 
	if (!close)
	{
		bool	 in = false;
		for (int 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
	{
		double sum = 0;
		for (int i=0;i<N; i++)
		{
			cmplex Z=(z-0.5*(zv[i]+zv[i+1]))/(0.5*(zv[i+1]-zv[i]));
			sum+=logZm1oZp1(Z).imag();
		}
		return std::abs(sum)>=REALSMALL;
	}
}