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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。* ]) b& {$ C; v

' O2 {  M/ m3 h' }6 U. X% a　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。* y  U! F" C; T

* j7 u8 d/ Z- K5 l. ~6 A9 v　　首先定义点结构如下：8 J  A6 o! A6 n& \8 [- U

/ X! v7 T- E7 V& w+ k以下是引用片段：8 v7 L4 D' J+ `0 V; D$ n. {* I( V
　　/* Vertex structure */
5 P; h4 B! H" f5 T# C　　typedef struct
% M6 w, U0 q/ g　　{ 8 M2 ^# n7 ^6 N3 R& G3 T0 u# D
　　double x, y; & Z, A. u6 M8 m" F7 S" S
　　} vertex_t;
0 W' Q/ e! p( _6 L& p/ j. o2 N; H: e- V/ V+ l6 v5 ]

* C! o8 r: p+ T! S# J) f　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
; @: q  t: u4 O# t( k/ }5 J6 b% b* J0 x% o* w
以下是引用片段：# k, O+ A- V4 g2 ]7 |1 N2 k
　　/* Vertex list structure – polygon */
- v* n- ^5 S1 T: C7 Y  G　　typedef struct . b! {  d( l, Q" n5 [& [/ U( O
　　{ - a* [7 V" E$ C# G# o1 j/ k7 @) J
　　int num_vertices; /* Number of vertices in list */
! U5 s) J4 \/ E2 k; [4 K　　vertex_t *vertex; /* Vertex array pointer */ * n# Z8 A8 T  U% Q
　　} vertexlist_t;
7 X0 n( z" [& ~; t; F2 U; ?7 X  Q) y7 T3 r
5 B& ]5 d$ W: ?$ E- D
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
6 ~% P- ~! g$ D- W7 C5 U& x* [9 o6 ~  }* m8 O0 @. A
以下是引用片段：; v9 l, o: x2 b! e9 C0 `8 z1 q
　　/* bounding rectangle type */
0 y% H, m, U+ d! M　　typedef struct 4 j, C9 _+ f5 i+ D  c
　　{
( S/ X) w; g% h. i  `+ a. W* N8 _　　double min_x, min_y, max_x, max_y; * j3 c4 j* `6 l* ~
　　} rect_t; , _+ `$ p# G! m8 }/ t1 _0 ]0 E2 d! S) D
　　/* gets extent of vertices */
! ~6 a5 t" [  p, f+ j- |　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ ( i$ S& [# x/ c, s3 a$ n
　　rect_t* rc /* out extent*/ )
# P5 ]$ M7 }$ u! |　　{
4 K' d: x% S% ?7 V3 R9 m9 f: s1 F' h& q　　int i;
8 l% m$ b: @/ x* ]; L( J6 z　　if (np > 0){
) C: H! v' U- `4 D' q　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
, I" v' f. D* s5 @1 g7 m　　}else{ ) X) [/ D. {2 c; n+ w
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
. H. h" w; X. B1 i, N　　} 7 v) Z* g8 Q1 r' V% R
　　for(i=1; i  3 k) B7 U: u; H; \5 J
　　{ , o4 J* N! `# T5 K! M
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
, J8 r! W9 ]& l　　if(vl.y < rc->min_y) rc->min_y = vl.y;
; c, }& O8 J$ P1 W5 A; i9 B) p　　if(vl.x > rc->max_x) rc->max_x = vl.x;   Y2 u' I3 I2 M$ O/ \+ m9 y
　　if(vl.y > rc->max_y) rc->max_y = vl.y;
+ y; V2 W: P  T　　} , d$ T* I; A  W5 o/ ~1 l7 V
　　} ! E& j6 _% T* J" G! O
+ K8 x7 q2 l. M7 @6 B$ Z) A
" G, t. X! K/ {2 ^) z* z+ X, v  L: ^
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。% |0 q5 ^, d' l; V& h% j: j

' _" h% h8 `8 L& u. F* r2 }　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：7 N( @- N& D1 \! C/ O

% F! U  t; t+ }, x2 `　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;) \! C' R9 |/ _8 L6 i0 `
: E0 a* ?* n& j9 N
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
! H; W$ W, G# L4 q6 O3 t; `* {9 J& W2 j; m# |$ S9 w$ T
以下是引用片段：
. h1 Z3 ]5 i% u8 @) i5 E4 E　　/* p, q is on the same of line l */ 4 Q3 v. \7 d# S: O0 H
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ , H' O' ^5 }* C7 r
　　const vertex_t* p, : v' b$ {+ T( v' ~, Y3 C
　　const vertex_t* q)
' f: Y$ X) w3 N/ V. e' S　　{
/ n$ n- s- }7 i9 u! `* k5 O　　double dx = l_end->x - l_start->x;
( k' z/ e. P6 @! m: s　　double dy = l_end->y - l_start->y; ' O+ U& t8 r1 ?6 I' m5 y/ i
　　double dx1= p->x - l_start->x; : Z& ^" `7 i- K+ J: ?! G% i; c9 C
　　double dy1= p->y - l_start->y;
: S; ]' ~6 {% S3 J, Y$ `0 S7 y% W1 B　　double dx2= q->x - l_end->x; . |: `' X6 S7 R, S( G: A
　　double dy2= q->y - l_end->y; 3 K7 O+ R4 ]$ [! f$ y
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
# T( X, Y" F( K" C! f' x+ ~7 X. G　　} . v+ F2 g$ n0 u9 L" c, u  f
　　/* 2 line segments (s1, s2) are intersect? */ 7 t: Z# c/ ^& g7 ?0 y
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, " t5 j, t0 H0 @9 _5 d& y2 b
　　const vertex_t* s2_start, const vertex_t* s2_end)
( w' L) t1 H! e6 e9 {. D* }　　{
1 Z% r9 k3 C7 ]; R8 e: F　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
( E/ r  T; {" J: C7 d　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; : s- x( C# j, v, y
　　}
! i/ R! r' W# R3 W
* N7 g- t) h8 _* Z; P! }1 l7 D% ]; ^& S; h+ y
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：0 d4 g  L5 ^- m% w8 x" N( d1 |

( {1 s7 w0 k3 f9 Y. n- A以下是引用片段：5 [1 H7 S1 \& @- x4 [# j8 {
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
% \8 D* r0 z+ ~9 r$ K9 Y　　const vertex_t* v) ' G, N9 B6 [: u: `4 Z
　　{ 0 }& h/ U: q- Y' J, u' X- ?9 R
　　int i, j, k1, k2, c; , z. U' I, Q9 P* @5 i
　　rect_t rc; 6 t+ }& s' ^3 n1 l8 M5 R- D( t
　　vertex_t w; 6 S! G/ l( s& V. V9 k
　　if (np < 3)
+ e. q+ Z9 d) v' }$ d# o　　return 0; ! t, V9 r, s8 u
　　vertices_get_extent(vl, np, &rc);
, l! s5 J8 v/ d, ^* L9 G　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) - D9 ?) A/ a1 i! a+ l; J
　　return 0; 0 V" O% ?' U* X0 r
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
& P  T. U; P) F3 {9 Z　　w.x = rc.max_x + DBL_EPSILON;
. @% b1 f/ x* \+ e1 {$ D　　w.y = v->y; 9 C7 S  \1 X, }2 f; T5 C" s
　　c = 0; /* Intersection points counter */ " f1 _  R! ]" Q0 w8 j: j0 v
　　for(i=0; i  
& |) H0 j9 X; ~, Z: b* Z3 V　　{
; ]# z. ~* l5 W: U　　j = (i+1) % np; ; `6 z% T4 t2 S( r, L/ g
　　if(is_intersect(vl+i, vl+j, v, &w)) " n0 i! L+ f+ o8 G+ U+ E
　　{ " {: d) [" b0 J9 j$ ]3 d
　　C++; + d- N! Y: w' a/ J- l. M; \
　　}
1 R4 d! O- A) x9 d6 ]  g3 m: Q　　else if(vl.y==w.y)
" j8 d' ^- ?5 A% H5 o　　{ $ Y4 c0 J6 w+ n: [" ?5 e
　　k1 = (np+i-1)%np; . P6 X. X% D& A! d3 t
　　while(k1!=i && vl[k1].y==w.y) ' h1 o0 q) e" ?! @
　　k1 = (np+k1-1)%np;
  R# p) j! r8 I( ^1 Z0 g　　k2 = (i+1)%np; ! `$ c! ]) N0 R! e% h; Y
　　while(k2!=i && vl[k2].y==w.y)
/ e! b, }3 t; f7 M9 L　　k2 = (k2+1)%np;
0 H  e. J& Y& Q　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) 8 g6 R7 B# c, x9 u: H
　　C++; # h- |! A8 \3 N' _
　　if(k2 <= i) 7 |; Q) C2 N5 ?. ]4 x/ _1 X; u
　　break; 2 y5 D: R1 A6 q7 Z' I; s% Q
　　i = k2;
7 u% P9 v" @; x3 i) T, c　　}
4 w, V6 W5 Q4 T* [! B' q  F　　} ' m( j" u8 O' N0 U! W
　　return c%2; . w! y3 z# |- s5 C
　　} 8 j- E$ W9 o# W. z
! s. l, w( u: U. A1 A
2 v9 N8 _, e1 u9 m; W
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。