- 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;
}
}