C语言中显示 点在多边形内 算法

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
4 [' q( j+ ]. N5 Z4 d: K$ u) I9 o0 ?# ~6 Z$ ?& D
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。' X" R6 P$ G& C+ y

( @6 y7 o5 u6 R- N! y4 a' o　　首先定义点结构如下：+ w3 o% ^: @) Q

2 ^2 a$ G. b* Q: u3 _- c以下是引用片段：  Y4 L1 L- F, p3 F" A9 h! ?$ n
　　/* Vertex structure */
$ L2 T( s2 U4 h2 I1 K　　typedef struct 8 H/ Z# X/ F5 \) J: H' C: W
　　{ 3 y2 I: Y* f- N! e
　　double x, y; : J/ S3 x  z% e5 A# K
　　} vertex_t;
; s4 X. m5 R2 a* j! p2 f, v( a* j* H/ a  @4 T0 c
" `! V: m- j9 ]' _2 d/ H. B  u$ {
　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：- Z, O$ T/ c, }

8 [0 v$ n: y7 U/ _" r* x; Z1 q  g: o以下是引用片段：. Y# k- q4 d- I7 t
　　/* Vertex list structure – polygon */ 9 d/ X; Q/ E6 J6 f/ N( e6 i
　　typedef struct 7 F& r. I% s( A$ Y7 z; e2 H
　　{
) q. ^6 [* O& p. B$ U/ O- {　　int num_vertices; /* Number of vertices in list */ # {$ w. `7 O& ^: |4 {
　　vertex_t *vertex; /* Vertex array pointer */
6 K. g$ ]9 x! v, [( ?4 I9 I　　} vertexlist_t;
7 m" B. \$ |' F  N( H; Q# C& G8 A2 M1 ]# J* N" I! g

  v3 Z- }' T' a& G; d　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
' B8 h+ v" O% w+ A5 n
9 `+ _$ r/ _2 D  B% {1 L以下是引用片段：5 z! T( R- A$ c; ?6 `
　　/* bounding rectangle type */
/ g$ u7 r  P0 C& S　　typedef struct . j! V9 U, K! E( F9 B" z
　　{ / L1 S# u5 Z3 V4 J0 s) w3 X2 A; Q
　　double min_x, min_y, max_x, max_y; " [, E# S7 b" B  V7 Y3 L: k
　　} rect_t; : J$ u4 o+ E& g
　　/* gets extent of vertices */   j5 u1 \# \. k2 L
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ + n" K+ i! C- W7 r- }' Z# r1 \
　　rect_t* rc /* out extent*/ ) 1 g! n7 s9 d& F3 ]% ]* `! F
　　{ % q, }; b6 Q3 o. Q5 s% j: L8 J: {
　　int i; * d9 D8 F/ Y9 A
　　if (np > 0){ ! r1 i  \) [  R* i6 k
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
& D/ h; _9 q: g% _% K! o# Z　　}else{
3 V; o2 l% {; U* d  h　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ + X/ s) U' @' j* O. y1 \
　　} * P+ [1 a: R( K+ D- g
　　for(i=1; i  / C+ B/ ^1 G. J6 h/ f
　　{ $ G' `8 G. b0 B% P0 O% }) O
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
5 |6 Z& x4 _$ n* D, X: ]3 ?- i　　if(vl.y < rc->min_y) rc->min_y = vl.y;
! s; E! I0 [5 ?% n　　if(vl.x > rc->max_x) rc->max_x = vl.x; " L1 s, L$ N6 e) s( _
　　if(vl.y > rc->max_y) rc->max_y = vl.y; ' z5 \9 v- J6 G, F, ~
　　} 4 P) v9 y+ _# c: w( X
　　} ( m) {% D0 f6 J% z, D1 D, |: [6 C
# d  h* f& p- k5 o( Q
& Y3 P- s5 P! c3 f" [% G  l/ Y" q% s
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
. ~1 g; @( J0 V7 `8 C! o& T: k( G) i. U
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
) T; L" o; {% H; Z) T& e# O0 d; Q5 ^% A( k& l8 x1 j
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
. Z' j# K- X0 v/ \7 z% x1 C1 V5 E' A8 U3 U% b' k7 A
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;7 E2 o) O! g; k, Z0 H
( C2 Q) H; ]+ z. q$ c
以下是引用片段：
: v  M7 c) O: S0 O# [: }2 G7 N0 c　　/* p, q is on the same of line l */ ) O5 T5 p2 a' [4 u
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
7 I; `2 y8 A/ y: B　　const vertex_t* p,
0 A* \* _% V; _. Y/ b　　const vertex_t* q)
; ]7 L% G" i7 F5 A3 C$ p' r, m　　{
: {4 s7 k+ _$ A' d　　double dx = l_end->x - l_start->x; / x( u& R7 {7 ]8 D$ D( d3 Q6 L
　　double dy = l_end->y - l_start->y; : n6 x  y8 K& X/ I4 B$ c: @
　　double dx1= p->x - l_start->x;
& P% y+ |- W; I3 i( K" @* e　　double dy1= p->y - l_start->y;
0 E! F1 y* U# m- S( C% o, Q　　double dx2= q->x - l_end->x;
, O) e- ?( h" y5 ^* l5 A& X9 C　　double dy2= q->y - l_end->y;
$ V% H: Q; L) ^, i7 w4 C0 N　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
* \5 b& ~# D- r/ N　　}
0 [5 H) {" A& N+ y　　/* 2 line segments (s1, s2) are intersect? */
# d  y9 R# e, B6 h7 F3 k8 O　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
$ _2 o" R  I2 j( P# s　　const vertex_t* s2_start, const vertex_t* s2_end) 4 m$ \+ Y3 y9 J* v9 j( ^( ^1 z
　　{
# n: s) d8 p& J　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && : D+ ^2 {6 ]& L8 @6 b
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
5 d, O& @, x4 u. D　　} $ j/ W; G7 T5 }4 E# K0 U/ R
+ I6 V8 ?+ c2 H
9 x+ d* Z5 h8 M' }
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
! k7 [& F9 d3 p: |5 E* b1 p) @/ M7 J: p
以下是引用片段：
. T9 g  P7 b2 P; H) u9 w6 H　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ / W0 ?& z* W8 K: U( j
　　const vertex_t* v)
1 f- y# z# M* U; d6 u, t　　{
& A" @0 F2 g6 ?( d: K: c# L/ c　　int i, j, k1, k2, c; * J& I3 D9 @. ?$ c
　　rect_t rc; 4 h( u! l) {  p2 h0 W& d% j* n7 I
　　vertex_t w;
- j6 Y9 p2 H" i6 T　　if (np < 3)
1 t* L2 r& P$ p$ ]- H　　return 0;
0 ~' W: ~( [& v　　vertices_get_extent(vl, np, &rc);
: d! ?( ?2 u2 J　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
1 O) u6 c! ]0 h$ h" A　　return 0;
9 [2 ^6 }4 {+ z4 x0 e" M　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
7 z0 X6 j1 k3 e% D( A0 a8 u　　w.x = rc.max_x + DBL_EPSILON;
' ?8 n) c; g+ f. C9 N* f" a7 u, ^# ^$ d　　w.y = v->y;
5 r+ V# Z: K% a+ d* Z& s% X) H- y　　c = 0; /* Intersection points counter */ 1 o" w" g4 Z# j$ b& ^* W
　　for(i=0; i  & ^# h3 X6 t: c. U
　　{ ; `" l: J! G. A) Q- x  R/ r
　　j = (i+1) % np; $ P4 ^9 E( u1 y, Q1 J
　　if(is_intersect(vl+i, vl+j, v, &w)) . B5 L' Q1 ~- W  [
　　{ 9 U% y" |8 X$ t1 G4 g& P
　　C++; ' C4 h7 R" m) B6 i& z' S$ ?9 j
　　}
3 m( N0 v9 s0 {　　else if(vl.y==w.y) . |! G1 n2 p; v; B5 D4 E
　　{ 2 V0 [' ]* O0 j, n. A" B, o. M. k
　　k1 = (np+i-1)%np; 5 x$ [5 p& S) g) F9 F
　　while(k1!=i && vl[k1].y==w.y) ' i4 g: ?" ]: @* J( ~  B
　　k1 = (np+k1-1)%np;
* ^5 \8 W" b, B7 @9 Q* u　　k2 = (i+1)%np;
' Q% R! w9 ]1 {+ p  J7 p　　while(k2!=i && vl[k2].y==w.y)
; o) B6 [7 w5 Y　　k2 = (k2+1)%np; ; w: c# m  X+ ]5 u: g* L
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) & z6 t7 j" j, r  D+ y% Q! ^4 v
　　C++;
1 G! p% C2 J. A, h8 A6 l: `' v　　if(k2 <= i) ' f+ d+ D: u  M; ^; ?! E) _
　　break; % H* }  B9 ~4 q# s( m3 H
　　i = k2;
. w. [; Z4 W% s. X, s$ i: Y　　} 5 T7 P$ O2 I2 L- s
　　}
: u9 [8 K& M) t; S　　return c%2;
8 h* K) e! w0 t+ W　　} 6 z, g# O+ t% \7 J4 j

: K* S3 A4 g, m1 k) c
: i% A; U7 [# x　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。