本文介绍了两种判断点是否在多边形内的算法：射线交点法和绕数法。射线交点法通过统计射线与多边形边的交点数量判定，而绕数法则计算多边形绕点的次数。针对特定测试点，不同算法可能得出不同结果，例如在某些实现中，点落在边上的情况可能会导致误判。
判断一点P是否在多边形内部的算法一般有两种：射线交点法(ray crossing number method)，绕数法(windig number method)
射线交点法(ray crossing number method)
以从此点为原点引一射线，然后统计此射线与多边形各边交点的数量。如果与射线相交的边的数量是偶数，则此点落在多变形外，否则落在多边形内。
绕数法(windig number method)
统计多边形绕点P的次数(绕数)。如果这个绕数是0，则点在多边形内，否则点在多边形外。
[1]给出了射线交点法的C++语言的一个实现，根据[1]给出的代码，将点坐标数据类型由int转换为double
#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <vector>
static double ESP = 0.0001;
double max( double a, double b )
{
return (a > b)? a:b;
}
double min( double a, double b )
{
return (a < b)? a:b;
}
class Grid
{
public:
Grid()
{
_X = 0.0;
_Y = 0.0;
_Z = 0.0;
}
Grid( double X, double Y, double Z )
{
_X = X;
_Y = Y;
_Z = Z;
}
;
virtual ~Grid()
{
}
;
public:
bool isSame(Grid theGrid )
{
bool result = false;
if (fabs(_X - theGrid._X) > ESP)
{
return result;
}
if (fabs(_Y - theGrid._Y) > ESP)
{
return result;
}
if (fabs(_Z - theGrid._Z) > ESP)
{
return result;
}
result = true;
return result;
}
void show()
{
printf("[%f, %f, %f]\n", _X, _Y, _Z);
}
public:
double _X;
double _Y;
double _Z;
};
// http://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/
// Define Infinite (Using INT_MAX caused overflow problems)
#define INF 10000
double max( double a, double b )
{
return (a > b)? a:b;
}
double min( double a, double b )
{
return (a < b)? a:b;
}
//struct Point
//{
// int x;
// int y;
//};
// http://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/
// Given three colinear points p, q, r, the function checks if
// point q lies on line segment 'pr'
bool onSegment(Grid p, Grid q, Grid r)
{
if (q._X <= max(p._X, r._X) && q._X >= min(p._X, r._X) &&
q._Y <= max(p._Y, r._Y) && q._Y >= min(p._Y, r._Y))
return true;
return false;
}
// http://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/
// To find orientation of ordered triplet (p, q, r).
// The function returns following values
// 0 --> p, q and r are colinear
// 1 --> Clockwise
// 2 --> Counterclockwise
int orientation(Grid p, Grid q, Grid r)
{
int val = (q._Y - p._Y) * (r._X - q._X) -
(q._X - p._X) * (r._Y - q._Y);
if (val == 0) return 0; // colinear
return (val > 0)? 1: 2; // clock or counterclock wise
}
// http://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/
// The function that returns true if line segment 'p1q1'
// and 'p2q2' intersect.
bool doIntersect(Grid p1, Grid q1, Grid p2, Grid q2)
{
// Find the four orientations needed for general and
// special cases
int o1 = orientation(p1, q1, p2);
int o2 = orientation(p1, q1, q2);
int o3 = orientation(p2, q2, p1);
int o4 = orientation(p2, q2, q1);
// General case
if (o1 != o2 && o3 != o4)
return true;
// Special Cases
// p1, q1 and p2 are colinear and p2 lies on segment p1q1
if (o1 == 0 && onSegment(p1, p2, q1)) return true;
// p1, q1 and p2 are colinear and q2 lies on segment p1q1
if (o2 == 0 && onSegment(p1, q2, q1)) return true;
// p2, q2 and p1 are colinear and p1 lies on segment p2q2
if (o3 == 0 && onSegment(p2, p1, q2)) return true;
// p2, q2 and q1 are colinear and q1 lies on segment p2q2
if (o4 == 0 && onSegment(p2, q1, q2)) return true;
return false; // Doesn't fall in any of the above cases
}
// http://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/
// Returns true if the point p lies inside the polygon[] with n vertices
bool isInside(std::vector<Grid> polygon, Grid p)
{
int n = polygon.size();
// There must be at least 3 vertices in polygon[]
if (n < 3) return false;
// Create a point for line segment from p to infinite
Grid extreme = Grid(INF, p._Y, 0);
// Count intersections of the above line with sides of polygon
int count = 0, i = 0;
do
{
int next = (i+1)%n;
// Check if the line segment from 'p' to 'extreme' intersects
// with the line segment from 'polygon[i]' to 'polygon[next]'
if (doIntersect(polygon[i], polygon[next], p, extreme))
{
// If the point 'p' is colinear with line segment 'i-next',
// then check if it lies on segment. If it lies, return true,
// otherwise false
if (orientation(polygon[i], p, polygon[next]) == 0)
return onSegment(polygon[i], p, polygon[next]);
count++;
}
i = next;
} while (i != 0);
// Return true if count is odd, false otherwise
return count&1; // Same as (count%2 == 1)
}
int main( int argc, char* argv )
{
Grid P1( 10.000000, 20.000000, 10.000000 );
Grid P2( -10.000000, 20.000000, 10.000000 );
Grid P3( -10.000000, 10.000000, -10.00 );
Grid P4( 10.000000, -10.000000, 10.0 );
std::vector<Grid> polygon1 ;
polygon1.push_back( P1 );
polygon1.push_back( P2 );
polygon1.push_back( P3 );
polygon1.push_back( P4 ) ;
Grid section2P1( -10.000000, 10.000000, 10.000000);
isInside( polygon1, section2P1)? printf("Inside.\n"): printf("Outside.\n");
}
上述代码的算例中，测试点落在多边形的边上，[1]算法给出的结果是inside，
对于判定一点是否在多边形内，[2][3] 给出了绕数法的C++语言代码， 但对 这个算例给出的结果是outside。
[4]分别给出了射线交点法和绕数法的C++语言代码， 但对 这个算例给出的结果都是outside。下面是[4]的代码
#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <vector>
#include <map>
static double ESP = 0.01;
static double dotMultiply(double* v1, double* v2)
{
double result = v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
return result;
}
static void crossMultiply(double* v1, double* v2, double* result)
{
result[0] = v1[1] * v2[2] - v1[2] * v2[1];
result[1] = v1[2] * v2[0] - v1[0] * v2[2];
result[2] = v1[0] * v2[1] - v1[1] * v2[0];
return;
}
class Grid
{
public:
Grid()
{
_ID = -1;
_CoordID = 0;
_X = 0.0;
_Y = 0.0;
_Z = 0.0;
}
Grid( double X, double Y, double Z )
{
_ID = -1;
_CoordID = 0;
_X = X;
_Y = Y;
_Z = Z;
}
;
virtual ~Grid()
{
}
;
public:
bool isSame(Grid theGrid, double error)
{
bool result = false;
if (fabs(_X - theGrid._X) > error)
{
return result;
}
if (fabs(_Y - theGrid._Y) > error)
{
return result;
}
if (fabs(_Z - theGrid._Z) > error)
{
return result;
}
if (_CoordID != theGrid._CoordID)
{
return result;
}
result = true;
return result;
}
void show()
{
printf("Grid %d: [%f, %f, %f]\n", _ID, _X, _Y, _Z);
}
public:
int _ID;
int _CoordID;
double _X;
double _Y;
double _Z;
};
//
// http://geomalgorithms.com/a03-_inclusion.html
//
// Copyright 2000 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
// a Point is defined by its coordinates {int x, y;}
//===================================================================
// isLeft(): tests if a point is Left|On|Right of an infinite line.
// Input: three points P0, P1, and P2
// Return: >0 for P2 left of the line through P0 and P1
// =0 for P2 on the line
// <0 for P2 right of the line
// See: Algorithm 1 "Area of Triangles and Polygons"
inline int
isLeft( Grid P0, Grid P1, Grid P2 )
{
return ( (P1._X - P0._X) * (P2._Y - P0._Y)
- (P2._X - P0._X) * (P1._Y - P0._Y) );
}
//===================================================================
// cn_PnPoly(): crossing number test for a point in a polygon
// Input: P = a point,
// V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
// Return: 0 = outside, 1 = inside
// This code is patterned after [Franklin, 2000]
int cn_PnPoly( Grid P, std::vector<Grid> V, int n )
{
int cn = 0; // the crossing number counter
// loop through all edges of the polygon
for (int i=0; i<n; i++) { // edge from V[i] to V[i+1]
if (((V[i]._Y <= P._Y) && (V[i+1]._Y > P._Y)) // an upward crossing
|| ((V[i]._Y > P._Y) && (V[i+1]._Y <= P._Y))) { // a downward crossing
// compute the actual edge-ray intersect x-coordinate
float vt = (float)(P._Y - V[i]._Y) / (V[i+1]._Y - V[i]._Y);
if (P._X < V[i]._X + vt * (V[i+1]._X - V[i]._X)) // P._X < intersect
++cn; // a valid crossing of y=P._Y right of P._X
}
}
return (cn&1); // 0 if even (out), and 1 if odd (in)
}
//===================================================================
// wn_PnPoly(): winding number test for a point in a polygon
// Input: P = a point,
// V[] = vertex points of a polygon V[n+1] with V[n]=V[0]
// Return: wn = the winding number (=0 only when P is outside)
int wn_PnPoly( Grid P, std::vector<Grid> V, int n )
{
int wn = 0; // the winding number counter
// loop through all edges of the polygon
for (int i=0; i<n; i++) { // edge from V[i] to V[i+1]
if (V[i]._Y <= P._Y) { // start y <= P._Y
if (V[i+1]._Y > P._Y) // an upward crossing
if (isLeft( V[i], V[i+1], P) > 0) // P left of edge
++wn; // have a valid up intersect
}
else { // start y > P._Y (no test needed)
if (V[i+1]._Y <= P._Y) // a downward crossing
if (isLeft( V[i], V[i+1], P) < 0) // P right of edge
--wn; // have a valid down intersect
}
}
return wn;
}
//===================================================================
int main( int argc, char* argv )
{
//
Grid section1P3( -10.000000, 10.000000, -10.00 );
Grid section1P4( 10.000000, -10.000000, 10.0 );
Grid section1P1( 10.000000, 20.000000, 10.000000 );
Grid section1P2( -10.000000, 20.000000, 10.000000 );
std::vector<Grid> prism1BaseVertex ;
prism1BaseVertex.push_back( section1P1 );
prism1BaseVertex.push_back( section1P2 );
prism1BaseVertex.push_back( section1P3 );
prism1BaseVertex.push_back( section1P4 ) ;
Grid section2P1( -10.000000, 10.000000, 10.000000);
int r = wn_PnPoly( section2P1, prism1BaseVertex, prism1BaseVertex.size() );
printf("winding number result = %d.\n", r );
r = cn_PnPoly( section2P1, prism1BaseVertex, prism1BaseVertex.size() );
printf("crosing number result = %d.\n", r );
}
Tarjan’s Algorithm to find Strongly Connected Components[4]
// A C++ program to find strongly connected components in a given
// directed graph using Tarjan's algorithm (single DFS)
#include<iostream>
#include <list>
#include <stack>
#define NIL -1
using namespace std;
double max( double a, double b )
{
return (a > b)? a:b;
}
double min( double a, double b )
{
return (a < b)? a:b;
}
// A class that represents an directed graph
class Graph
{
int V; // No. of vertices
list<int> *adj; // A dynamic array of adjacency lists
// A Recursive DFS based function used by SCC()
void SCCUtil(int u, int disc[], int low[],
stack<int> *st, bool stackMember[]);
public:
Graph(int V); // Constructor
void addEdge(int v, int w); // function to add an edge to graph
void SCC(); // prints strongly connected components
};
Graph::Graph(int V)
{
this->V = V;
adj = new list<int>[V];
}
void Graph::addEdge(int v, int w)
{
adj[v].push_back(w);
}
// A recursive function that finds and prints strongly connected
// components using DFS traversal
// u --> The vertex to be visited next
// disc[] --> Stores discovery times of visited vertices
// low[] -- >> earliest visited vertex (the vertex with minimum
// discovery time) that can be reached from subtree
// rooted with current vertex
// *st -- >> To store all the connected ancestors (could be part
// of SCC)
// stackMember[] --> bit/index array for faster check whether
// a node is in stack
void Graph::SCCUtil(int u, int disc[], int low[], stack<int> *st,
bool stackMember[])
{
// A static variable is used for simplicity, we can avoid use
// of static variable by passing a pointer.
static int time = 0;
// Initialize discovery time and low value
disc[u] = low[u] = ++time;
st->push(u);
stackMember[u] = true;
// Go through all vertices adjacent to this
list<int>::iterator i;
for (i = adj[u].begin(); i != adj[u].end(); ++i)
{
int v = *i; // v is current adjacent of 'u'
// If v is not visited yet, then recur for it
if (disc[v] == -1)
{
SCCUtil(v, disc, low, st, stackMember);
// Check if the subtree rooted with 'v' has a
// connection to one of the ancestors of 'u'
// Case 1 (per above discussion on Disc and Low value)
low[u] = min(low[u], low[v]);
}
// Update low value of 'u' only of 'v' is still in stack
// (i.e. it's a back edge, not cross edge).
// Case 2 (per above discussion on Disc and Low value)
else if (stackMember[v] == true)
low[u] = min(low[u], disc[v]);
}
// head node found, pop the stack and print an SCC
int w = 0; // To store stack extracted vertices
if (low[u] == disc[u])
{
while (st->top() != u)
{
w = (int) st->top();
cout << w << " ";
stackMember[w] = false;
st->pop();
}
w = (int) st->top();
cout << w << "\n";
stackMember[w] = false;
st->pop();
}
}
// The function to do DFS traversal. It uses SCCUtil()
void Graph::SCC()
{
int *disc = new int[V];
int *low = new int[V];
bool *stackMember = new bool[V];
stack<int> *st = new stack<int>();
// Initialize disc and low, and stackMember arrays
for (int i = 0; i < V; i++)
{
disc[i] = NIL;
low[i] = NIL;
stackMember[i] = false;
}
// Call the recursive helper function to find strongly
// connected components in DFS tree with vertex 'i'
for (int i = 0; i < V; i++)
if (disc[i] == NIL)
SCCUtil(i, disc, low, st, stackMember);
}
// Driver program to test above function
int main()
{
cout << "\nSCCs in first graph \n";
Graph g1(5);
g1.addEdge(1, 0);
g1.addEdge(0, 2);
g1.addEdge(2, 1);
g1.addEdge(0, 3);
g1.addEdge(3, 4);
g1.SCC();
cout << "\nSCCs in second graph \n";
Graph g2(4);
g2.addEdge(0, 1);
g2.addEdge(1, 2);
g2.addEdge(2, 3);
g2.SCC();
cout << "\nSCCs in third graph \n";
Graph g3(7);
g3.addEdge(0, 1);
g3.addEdge(1, 2);
g3.addEdge(2, 0);
g3.addEdge(1, 3);
g3.addEdge(1, 4);
g3.addEdge(1, 6);
g3.addEdge(3, 5);
g3.addEdge(4, 5);
g3.SCC();
cout << "\nSCCs in fourth graph \n";
Graph g4(11);
g4.addEdge(0,1);g4.addEdge(0,3);
g4.addEdge(1,2);g4.addEdge(1,4);
g4.addEdge(2,0);g4.addEdge(2,6);
g4.addEdge(3,2);
g4.addEdge(4,5);g4.addEdge(4,6);
g4.addEdge(5,6);g4.addEdge(5,7);g4.addEdge(5,8);g4.addEdge(5,9);
g4.addEdge(6,4);
g4.addEdge(7,9);
g4.addEdge(8,9);
g4.addEdge(9,8);
g4.SCC();
cout << "\nSCCs in fifth graph \n";
Graph g5(5);
g5.addEdge(0,1);
g5.addEdge(1,2);
g5.addEdge(2,3);
g5.addEdge(2,4);
g5.addEdge(3,0);
g5.addEdge(4,2);
g5.SCC();
return 0;
}
Triangle
Triangle is a 2D quality mesh generator and Delaunay triangulator.
下载到源程序后，如果是Windows操作系统，还需要在triangle.h之前做些配置，如定义以下几个宏[8]：
#define REAL double
#define ANSI_DECLARATORS
#include "triangle.h"
#undef REAL
在triangle.c中定义宏：#define NO_TIMER。有了上面的宏定义，可以编译出一个triangle.exe程序了。如果要将triangle用在自己的程序中，还需要定义#define TRILIBRARY。更多宏定义可以参考源程序。
Doubly connected edge list(DCEL)
https://en.wikipedia.org/wiki/Doubly_connected_edge_list
http://www.holmes3d.net/graphics/dcel/
http://www.cnblogs.com/lihao102/archive/2013/04/14/3020238.html
http://www.sintef.no/projectweb/geometry-toolkits/
http://what-when-how.com/advanced-methods-in-computer-graphics/mesh-processing-advanced-methods-in-computer-graphics-part-2/
http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/model/
http://imr.sandia.gov/papers/imr22/IMR22_25_Dyedov.pdf
http://cg.cs.tsinghua.edu.cn/people/~Yongjin/CAD13-exact-geodesic.pdf
https://www.researchgate.net/publication/255983160_AHF_Array-based_half-facet_data_structures_for_mixed-dimensional_and_non-manifold_meshes
http://www.industrial-geometry.at/events/strobl2006/2006-06steiner.pdf
http://staff.ustc.edu.cn/~lgliu/Courses/DGP_2012_spring-summer/Courseware/DGP04_MeshDataStructure.pdf
https://www2.cs.arizona.edu/classes/cs437/fall07/Lecture2.pdf
http://www.iue.tuwien.ac.at/phd/halama/node110.html
https://github.com/PeddleSpam/PolygonGraph
Reference
[1] http://www.geeksforgeeks.org/how-to-check-if-a-given-point-lies-inside-a-polygon/
http://www.dcs.gla.ac.uk/~pat/52233/slides/Geometry1x1.pdf
[2] https://stackoverflow.com/questions/11716268/point-in-polygon-algorithm
[3] http://alienryderflex.com/polygon/
[4] http://geomalgorithms.com/a03-_inclusion.html
[5] Tarjan’s Algorithm to find Strongly Connected Components
http://www.geeksforgeeks.org/tarjan-algorithm-find-strongly-connected-components/
[6] https://en.wikipedia.org/wiki/Tarjan%27s_strongly_connected_components_algorithm
[7 http://www.ics.uci.edu/~eppstein/161/960220.html#sca
[8] http://www.cnblogs.com/opencascade/p/3632705.html
扫一扫
1411
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
