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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。8 ?9 N) B! e( ?" c! g+ P

) G- n- F& }3 Y: S- a( S　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。: X6 l7 E5 N0 [+ V; I
3 |, {! j1 v0 h8 ]$ H# q
　　首先定义点结构如下：( l/ N& j) L% D
, z! ]6 }7 P* x4 k' _' S& }
以下是引用片段：
3 z0 a( Y/ L. B; |% o: m, |　　/* Vertex structure */
/ J+ f& @( y+ d( s1 ~　　typedef struct 1 p% M; m$ }2 ~. A" P' f
　　{ " r2 f  T: y+ B7 \
　　double x, y; , U. W6 k6 S: O$ y' P' K: d
　　} vertex_t;
+ t% A0 R& f0 s" ^& D- G+ M' R/ ]. k# C
& C" V- V* P6 i/ A- Y
　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：* k+ d( n4 J+ M# s3 V, L/ p
+ q) q( B+ t& p! U$ |7 C3 h& N5 |
以下是引用片段：; s2 T' H/ S# d! I& _+ k
　　/* Vertex list structure – polygon */
4 S2 x8 ^' t7 X% E" E　　typedef struct % T. P/ j1 m/ Q8 B& I* D
　　{
" O, F0 _! b: Q, `　　int num_vertices; /* Number of vertices in list */
! c3 v# V  a* e　　vertex_t *vertex; /* Vertex array pointer */
! p& w; e. f7 [& e/ j　　} vertexlist_t; 2 x$ l- w7 E  n" B1 [8 q4 T, t4 F

" n* u0 H) ]/ E2 k* `& R
1 V, N) {* x+ C: s5 Q+ J0 R! `1 X　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
& m- x$ j# H" S- E) f3 `4 Q$ p7 ?$ U$ g8 S8 m0 ^- x
以下是引用片段：. E! x2 f4 V. f' s6 n0 ~
　　/* bounding rectangle type */ , d8 o% J. w; I
　　typedef struct ' ^' L, N' T% @
　　{ * [7 V! ~9 D  Z( z4 b! B2 X- r& j7 `
　　double min_x, min_y, max_x, max_y;
+ S$ M3 a0 L2 O9 ]# q0 C　　} rect_t;
! y3 U$ k9 G  ~" o% P　　/* gets extent of vertices */
5 Q8 O! }* b0 Y+ {　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
. I% r' U0 O; p9 B. W1 e/ F　　rect_t* rc /* out extent*/ ) 9 X% O: L  [( h8 h) S" _, m
　　{
7 z/ Z) c" Z1 E5 [　　int i;
5 m, y1 w) W; h2 C9 J- {* k; s5 P　　if (np > 0){
8 c- k: v7 _& ]$ e5 ^　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; + U3 X/ b7 k4 l  s" C4 @) h
　　}else{
8 }) q" m! Z: f6 y7 a　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
) }( O$ L7 S1 a3 w! _- j* ?　　}
7 o# t; ^8 H% @! T　　for(i=1; i  + f% ~: X+ o) u. ~1 o$ g7 @8 D
　　{ ' h) E9 U- o4 U  J6 @& I  i0 F! q$ Q
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
7 O. I, A& h# m, w　　if(vl.y < rc->min_y) rc->min_y = vl.y; 7 H: S2 y( S' I: j! A/ C. k6 o( l
　　if(vl.x > rc->max_x) rc->max_x = vl.x;
( _! S0 e+ P+ h% V7 J6 F　　if(vl.y > rc->max_y) rc->max_y = vl.y;
& Q2 k- @/ p- H1 [　　} 9 f* U1 l+ p: J2 P/ R9 a
　　} 3 n* \# n" H# }2 o" `1 u
/ S1 F' Z7 Y. J! V6 F& [

. h; h/ d- R+ Y, G7 m　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
; t( c$ ~1 B' R; H6 M
4 Q- `5 u. X' X6 X- u9 _. H' c0 ]3 E　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
: ?( k3 R3 ]. L( ~. K0 \$ k; ?, F2 s, F( i
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;  w* R: H, t- [. {' g

2 ]" k5 W& ?3 X6 {　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;# Q" h5 {7 _; L) c5 o

1 P7 Q3 S* f. ?4 N$ [以下是引用片段：9 l/ C0 n4 k1 K7 w' d& o7 `; Z
　　/* p, q is on the same of line l */
# E& A  n; l; W' X5 Y% V4 Z　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
/ t3 T' N, O& z( z& i( c; K　　const vertex_t* p, ( }7 Z* t" o% o$ @* g, T
　　const vertex_t* q) . z" K2 s2 M7 z! G& _  \) ~
　　{
0 ]; Y! U. J/ T( C9 f8 H8 Z　　double dx = l_end->x - l_start->x;
# O5 c+ L; q0 _  S& S" h1 a. }　　double dy = l_end->y - l_start->y; " j/ g, `3 m2 u$ a
　　double dx1= p->x - l_start->x; 4 N) @. y( O9 |8 T: [# y
　　double dy1= p->y - l_start->y;
2 A! u1 ^3 ?6 N7 p　　double dx2= q->x - l_end->x; 4 U% e. A( T6 J; n! B$ T' x3 I  X
　　double dy2= q->y - l_end->y; ( K2 K$ C$ }( \+ [
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
. q; Z0 K% Y: u9 X8 h　　} 4 Q5 Q$ x+ @" w& X+ Y. ]
　　/* 2 line segments (s1, s2) are intersect? */
8 F5 Y/ T, ]- P. O7 l; t9 `6 [　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, / O, V, o# x9 ^+ g5 P/ W  G
　　const vertex_t* s2_start, const vertex_t* s2_end)
: c" `% Z  p, u* E2 z0 d7 |　　{ ( o& K! d: r7 B3 x; D. P# s
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && 1 P8 _' C7 l. h% b, u+ o
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
2 Z! C  U# W1 R$ V　　}
. ]1 t2 L" v9 e  Z2 ~+ |5 C1 C6 q2 J% I) I* o5 B: X# K* E
" F/ U: O# e, X7 m8 E1 }5 d
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
4 `# I7 F0 j4 r  b4 I+ @1 z* R! _7 w0 _
以下是引用片段：4 s  `) n5 @2 Y: k
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
3 n1 V) w6 B: F2 _3 L  W9 C  J　　const vertex_t* v) / k: i* D- E* o0 k9 G9 ]
　　{
7 i7 h) u7 I, Q. Y2 h　　int i, j, k1, k2, c; ! f* w' `' A1 s. q
　　rect_t rc;
8 a' h: E, `% X/ s% M4 H　　vertex_t w; . x( B: m, k/ H; M( }- ~( ?+ ~
　　if (np < 3) * e. x' O0 Y5 i4 q& N$ t3 q
　　return 0;
( D8 s( z: s' @" R) ]  m　　vertices_get_extent(vl, np, &rc);
; @3 F# S; ]& I　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
( _* y2 e; Y6 }& [, p' |　　return 0; 6 _6 D' M2 t' ^
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ & G! M, a# b+ y) R, E
　　w.x = rc.max_x + DBL_EPSILON; 6 Y( V5 Y( r0 c: i
　　w.y = v->y; 0 s' `, C% o! T6 Q8 a3 S
　　c = 0; /* Intersection points counter */
0 F: S  }, }4 I% I　　for(i=0; i  
+ U( D5 `( G* b/ F% K% J% M　　{ : a9 U$ ~8 O9 V! I6 k+ T- G8 H* n) R
　　j = (i+1) % np; 9 V5 D; `! n; w% A( B4 W
　　if(is_intersect(vl+i, vl+j, v, &w))
  g5 `/ y) z$ ]1 y9 r; {" C　　{
( \# d" ]2 J0 k8 u& p　　C++; 7 R1 A/ h: V- B! m! Y# G7 n' i; t
　　} ) k3 y4 u1 r# b0 {
　　else if(vl.y==w.y)
1 v1 \6 x/ W' p5 `' e　　{ " z2 W) _, h* O' t4 {
　　k1 = (np+i-1)%np; 0 g: m: S& M! k/ {  p' b- K- |3 ?. s
　　while(k1!=i && vl[k1].y==w.y) % Z" l  @& {+ Y
　　k1 = (np+k1-1)%np;
$ Y2 M, A0 x! {　　k2 = (i+1)%np;
' Z$ q# I4 ^5 ?) O0 A3 W' q　　while(k2!=i && vl[k2].y==w.y) $ |' Q- p1 W3 P7 K5 q; u" ]) U
　　k2 = (k2+1)%np;
2 H  h2 S3 v; v- @/ }　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) ; ^% l8 m0 U  p" E
　　C++;
% g# d0 |, |) K2 m　　if(k2 <= i)
( w* u& c! R: g  s5 m# X　　break; ! t( \, D" r! U, T+ }7 t
　　i = k2;
$ {4 w/ @( {7 `4 k, d　　}
9 |: d& Q7 z! k1 b  j% u　　} ( q3 e3 Q) m8 p9 {5 u4 O" \$ X3 t
　　return c%2;
  \& c; _# y' }/ y$ v. N' @　　}
$ F3 Q& Y/ C1 y# r
2 f% L+ x1 u; G8 t. g6 {5 R
. I  ^/ s4 l4 w8 }- F, C& V1 J; L$ j　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。