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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
, \2 V& @5 w2 M6 U
. W* [5 K- q. }4 p3 B0 ^% I3 t　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
( e% i: Q: V. ?2 P9 g+ z$ ^9 a. Q% q. g, B; U, Z# J
　　首先定义点结构如下：3 E1 u. I+ ~# ?! f& E
. S  d, X+ `. L9 d& H& L
以下是引用片段：! T& |4 T2 P7 |$ |5 A, L4 u" `. q
　　/* Vertex structure */ 3 W9 U, E- z3 L4 z  L
　　typedef struct ; {3 |' g8 O. |- @4 S: A3 u; C
　　{
" _/ l! q, `# w! {6 K, [# Y　　double x, y;
8 r8 U  t1 w: \/ `; @+ l( k　　} vertex_t;
6 S* H7 y  d9 R: u& W& c: s
# f' z* g1 T4 \; G, \) @  v# U6 n/ I% H1 Y$ ~7 P0 N
　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
1 n) b) O0 c- J6 D0 S9 A! Z. G. Y, y$ h
以下是引用片段：
2 a. d6 E: _# ~　　/* Vertex list structure – polygon */ . _& @- L2 Z7 @$ S
　　typedef struct - M1 P5 B( Y/ g: r7 Z" u; {$ ^
　　{ 2 ^. C8 ]& v/ y0 a/ j8 N6 O2 `7 o. A5 D
　　int num_vertices; /* Number of vertices in list */
+ S) R. x  k  q' a" j% y　　vertex_t *vertex; /* Vertex array pointer */
$ o4 l- ]7 e" r/ B. x/ d" v9 q　　} vertexlist_t; 7 S4 I+ |/ s+ ]9 E) m
. W! `: m; \% P1 i3 F( P2 z4 A6 D
- B* y- l  M- e2 A! C
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：/ d  I; ^0 o7 ^0 s( ^# A

- {4 F! H; {' `: \) e9 Z6 X以下是引用片段：
" Q, Y$ |" q  d. d+ Z2 ^, J　　/* bounding rectangle type */ # E' P% r" s- G" r& @9 R# ]4 _0 }
　　typedef struct
& Y5 W" M1 z6 R2 C; `　　{ ) P6 x; s3 ~, ^, [3 d* x
　　double min_x, min_y, max_x, max_y; ' z, a$ a; R( o) Y; {! b
　　} rect_t; 0 ~& C* i. p7 c- y8 [
　　/* gets extent of vertices */
& Q% u9 s" t! n2 Q9 I- T# N　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ / H2 L; h# V1 R& M4 R' y& \9 u
　　rect_t* rc /* out extent*/ )
9 E( z; q" m. f! Y: K$ i% C　　{
8 \$ @% R* Z: \+ F7 A; i/ T' e　　int i; : l1 D4 b0 @' g9 x  q
　　if (np > 0){
& k7 [- _% j, O3 }7 o6 J/ M　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
$ g; h- L5 W+ @3 k' o　　}else{ % Y! N2 z, u0 ]( Z; l+ |
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ ! v( T: f0 {* N7 T( |
　　}
$ A, L+ U3 Z. _. N# |　　for(i=1; i  
( R: `0 g! j5 i/ H( G  q) g7 T　　{ ( L2 X- g: P% o- `* @
　　if(vl.x < rc->min_x) rc->min_x = vl.x; : v/ |( A# j7 w. r1 x/ K# D8 H+ y8 Q
　　if(vl.y < rc->min_y) rc->min_y = vl.y;
9 {5 z* @1 L9 B# b9 c3 L. T0 v2 l　　if(vl.x > rc->max_x) rc->max_x = vl.x; . s, A0 M# c; L1 v" P
　　if(vl.y > rc->max_y) rc->max_y = vl.y;
3 P' q2 H  A- y, K( z  _/ d　　}
# O1 p: F% B2 J# A2 k5 A2 V0 u0 {　　}
& O' _  S& ?0 F3 n. ], C0 \, `, \* ]4 J

8 k% f% G3 R  \( X: O* ?) I　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。" G. q3 V, Q' o& O2 ?& @  d
, O4 O1 S$ E2 S4 H2 c7 k
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：, W" V) ^7 K4 s6 m

. d# C* f1 @0 a4 d! L5 G# H　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;* _8 W4 f; N0 {- O3 {
& Y! y6 B/ K# p" j, {
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;4 p+ y9 f7 w4 n$ n# q% `
( o5 q1 _' n2 V5 M& c, u
以下是引用片段：; w, k9 i' D6 K% W! I* m
　　/* p, q is on the same of line l */ 3 J9 K8 ?9 p$ n$ x" d+ `
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ / S" A6 f" P2 H. y2 ?- t: |
　　const vertex_t* p, 2 U/ H7 |! V5 m- C3 G
　　const vertex_t* q)
% T7 T1 Y* U* @( U6 a　　{
5 o( E3 C# A0 Z  O/ j& ]( T& V　　double dx = l_end->x - l_start->x; 9 u# b+ J. r0 M! o. D
　　double dy = l_end->y - l_start->y; $ y3 K( q. A' S  }$ t
　　double dx1= p->x - l_start->x; % E" f. b+ Y' B2 x1 E; B0 R! w
　　double dy1= p->y - l_start->y;
7 m# H- V6 T8 Y4 U. B5 f% P" R8 u　　double dx2= q->x - l_end->x; " H& V7 S! u2 O- r
　　double dy2= q->y - l_end->y; % |# @7 P6 f! h; O; ~
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
* T/ r. m' U2 o2 ^* t) e, ?　　} : \4 C8 z& O' ~: k/ Z8 z. Q: F
　　/* 2 line segments (s1, s2) are intersect? */ ( b7 b3 H  Z6 L; }3 A
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, - n3 |0 G- V( i( |; _3 Y# r
　　const vertex_t* s2_start, const vertex_t* s2_end) ' e( c8 ?: S& `4 }
　　{
% m5 e6 ^: x4 ]- k; f$ _! R　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && 3 e+ i* \9 ^7 o3 {( x
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
% r! ]  i4 k& _0 F% _3 I. q- B　　}
% U, v/ \; Q- X) W7 i" O+ W# e, F# {

) T0 W% J9 ]3 i　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
* l/ f) {  j7 g0 s  S( u( \* I; U1 _3 l- i' C# z/ P. O  C
以下是引用片段：
; a7 j9 K0 R, P+ |　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ 4 ^! t! Y4 F  J, L* [3 F
　　const vertex_t* v) / p& c. e: r. R3 I  I- F6 Z
　　{
8 [5 b8 P. P4 ^　　int i, j, k1, k2, c; % q+ R+ W$ s  I) r/ ^4 I. Q& c
　　rect_t rc;
  h! ]. l' d0 \8 B+ _5 Z) H+ f+ L# \　　vertex_t w;
5 j+ @. f/ r- C! m5 w　　if (np < 3) 3 Q* {9 V; H, |, S9 L
　　return 0; ( d  S. Y" {1 Y+ H4 D! L, ]
　　vertices_get_extent(vl, np, &rc); & L; ^, ^$ g% P/ G0 l) X5 `
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
& H4 T' ]8 d$ x( e　　return 0; ; h4 H* {4 G, h# D9 A/ W' [" K
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ 6 O% A4 m- j- Z% _% }" m6 v' p
　　w.x = rc.max_x + DBL_EPSILON; 4 N$ ]. m# J* i8 E: J6 U8 z" r
　　w.y = v->y;
: F& e! w8 S' B1 A% X2 e7 E　　c = 0; /* Intersection points counter */
6 O2 \  ]3 }; x6 A4 j　　for(i=0; i  " v4 P2 t* G: u
　　{ ; ^% G1 M9 m$ m0 T' U
　　j = (i+1) % np; 5 t" t: n/ J0 U) ~8 y
　　if(is_intersect(vl+i, vl+j, v, &w))   S- `0 Y7 W, F, z! ~3 L% C* V7 {+ K
　　{
: |8 ?: B. f% \0 w; l$ F　　C++; 1 ~! U, j1 G% E. u
　　}
/ m  m6 D) ^- v" C9 l" A) F( R　　else if(vl.y==w.y)
9 j0 E& R, ^+ p　　{ : M5 ~1 Z3 G$ c. w2 P% @
　　k1 = (np+i-1)%np;
* S" g2 k- R3 `' R* t& d　　while(k1!=i && vl[k1].y==w.y) 1 t" }9 X, }8 v; H3 f" H% l
　　k1 = (np+k1-1)%np; 4 J3 L% D( y! Y7 i
　　k2 = (i+1)%np;
$ S1 \# X7 {/ L' B9 P　　while(k2!=i && vl[k2].y==w.y) 3 i) `8 h: n' Z4 j1 n# n4 ^
　　k2 = (k2+1)%np;
2 O1 v1 _# d# t6 [: O+ _　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) 7 B3 l8 A3 W: V: D! ]
　　C++;
3 j% R. z2 e2 q0 m) m　　if(k2 <= i)
' g. w/ D3 t) \; i　　break;
! [8 \$ |9 ~+ L  \6 m; ^; s6 R　　i = k2;
0 s1 o9 j% s+ G& [　　}
) Q& W0 h' ?: @- B+ G　　}
: E" t1 W4 u" x! E2 |- P& T* b; x- S　　return c%2; ) g8 `! T) r- G! r
　　} % d$ Q3 t/ i! ]) Z

/ r' D! k) G4 y' X  W" ^& A1 G2 |9 t7 P+ v$ d5 \/ v
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。