Point-in-polygon
From aemwiki
- 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; } }