Difference between revisions of "Point-in-polygon"
Analytophile (talk | contribs) |
m |
||
(One intermediate revision by one other user not shown) | |||
Line 4: | Line 4: | ||
; GW Gourmet Functions Used: [[logZp1oZm1]] | ; GW Gourmet Functions Used: [[logZp1oZm1]] | ||
− | #include < | + | #include <cmath> |
#include <complex> | #include <complex> | ||
#include <iostream> | #include <iostream> | ||
− | typedef complex<double> cmplex; | + | typedef std::complex<double> cmplex; |
const double CLOSE_TO_LINE=1.02;// ellipse radius around line which determines closeness | const double CLOSE_TO_LINE=1.02;// ellipse radius around line which determines closeness | ||
const double REALSMALL=1e-10; //almost zero | const double REALSMALL=1e-10; //almost zero | ||
Line 21: | Line 21: | ||
bool IsInsidePolygon(const cmplex &z, const cmplex *zv, const int N) | 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; | |
− | bool close | ||
− | if (abs(zv[0]-zv[N])>REALSMALL){ | ||
− | cout<<"PolyIsInside : Bad polygon"<<endl; | ||
exit(0); | exit(0); | ||
} | } | ||
//check if point is close to boundary of polygon | //check if point is close to boundary of polygon | ||
− | for (i=0; i<N; i++){ | + | for (int i=0; i<N; i++){ |
− | if (((abs(z-zv[i])+abs(z-zv[i+1]))/ | + | 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; | close=true;break; | ||
} | } | ||
Line 37: | Line 34: | ||
//if not close, use faster ray-tracing 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> | ||
+ | double y = z.imag(); | ||
+ | double x = z.real(); | ||
if (!close) | if (!close) | ||
{ | { | ||
− | bool in | + | bool in = false; |
− | + | for (int i=0, j=N-1; i<N; j=i++) { | |
− | for (i=0, j=N-1; i<N; j=i++) { | ||
if ((((zv[i].imag()<=y) && (y<zv[j].imag())) || | if ((((zv[i].imag()<=y) && (y<zv[j].imag())) || | ||
((zv[j].imag()<=y) && (y<zv[i].imag()))) && | ((zv[j].imag()<=y) && (y<zv[i].imag()))) && | ||
Line 56: | Line 54: | ||
else | else | ||
{ | { | ||
− | + | double sum = 0; | |
− | double sum | ||
for (int i=0;i<N; i++) | for (int i=0;i<N; i++) | ||
{ | { | ||
− | Z=(z-0.5*(zv[i]+zv[i+1]))/(0.5*(zv[i+1]-zv[i])); | + | cmplex Z=(z-0.5*(zv[i]+zv[i+1]))/(0.5*(zv[i+1]-zv[i])); |
sum+=logZm1oZp1(Z).imag(); | sum+=logZm1oZp1(Z).imag(); | ||
} | } | ||
− | return | + | return std::abs(sum)>=REALSMALL; |
} | } | ||
} | } |
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; } }