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

zw2004

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

+ n7 C2 D& W7 X( P6 O　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
' f4 z' a8 Z1 x* b' L6 T1 a, Z
' P& f4 R+ A# D　　首先定义点结构如下：
2 s0 p6 P$ W( j; h; M! E$ S, a. B, D5 {6 m5 E/ Y) M# q) b  j" q% n
以下是引用片段：
9 C6 n9 {0 d7 F% @　　/* Vertex structure */ 9 O' O- A; X. C* A" d. k
　　typedef struct # p# ~" y$ _8 n9 A! Q$ t) S( `
　　{
. h3 W7 v' U% z" A' S# _4 m　　double x, y;
3 u2 P0 f" y7 y& m  a9 Z; Q3 B　　} vertex_t; 9 M0 \, d: ]5 U8 j3 W. @
/ Y7 R" `# }. w  Q6 v/ Q" R

$ U/ ?) q7 ?9 y# Y3 U1 e　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：4 J- k+ g. o% }7 ^5 v8 L
" x* V) U7 y- X8 x
以下是引用片段：: y/ a- u8 T3 t
　　/* Vertex list structure – polygon */
+ s% @% y! X3 N! e　　typedef struct ( O0 b( ~- m! `& ]8 h
　　{ $ p  u0 W' n* L# z! s
　　int num_vertices; /* Number of vertices in list */
: c( x) N' y- Y( U  l7 g1 |7 T9 g3 W　　vertex_t *vertex; /* Vertex array pointer */ . Y+ t8 t* ]! _8 ^, C$ C! K# l. D
　　} vertexlist_t; 2 o- n  U( u% e
2 Y+ m2 Z/ q# {+ L
1 Z( L9 b/ ?8 Q2 f% ?
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：. ^( O. N3 t1 f5 o9 E" ^$ Z
" j$ K  d) Z1 J1 H% X( ^4 V( [
以下是引用片段：
& R3 P/ P) X, u; f　　/* bounding rectangle type */
$ D' ~; w) j8 g; X' ^2 ~0 h　　typedef struct
. f$ k0 x4 u1 ~* q' s  u　　{ $ P1 H9 r4 Q# z) |8 P( p
　　double min_x, min_y, max_x, max_y; $ O  r  F0 L1 i2 Z
　　} rect_t; 4 F  [( z8 v9 K0 v! w" F) A
　　/* gets extent of vertices */ - C  F" N; y" `" O; K+ K/ k
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
& y& X% e5 u5 u: T* v　　rect_t* rc /* out extent*/ )
3 i1 u& g5 `2 j- s3 X7 d6 t9 s& ]　　{ * k# K- I) O7 }$ v
　　int i;
$ W) E! M2 s0 H; o　　if (np > 0){ ' X7 G+ W) ~) M) h# z$ G  H6 S/ B
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; 4 ~% o# {2 W" I. X$ q$ ]
　　}else{ ; L: X* b5 X( [, Y" V. z0 B: h, f
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
3 U. L" @$ _% X- D. ]　　} 4 J$ v6 t! r" _% o5 |- [8 U/ J
　　for(i=1; i  
9 O# T5 e  Z$ F4 k; @6 R" \　　{
6 f8 h, ]: V+ y% u9 Y9 b! S　　if(vl.x < rc->min_x) rc->min_x = vl.x;   O9 Q5 {: N& w& ]2 c  I# U. W
　　if(vl.y < rc->min_y) rc->min_y = vl.y;
' x4 M& y4 ]; O9 v　　if(vl.x > rc->max_x) rc->max_x = vl.x;
3 I2 B- P- v6 p& _% O$ D3 x% g　　if(vl.y > rc->max_y) rc->max_y = vl.y;
6 l( I6 _: }6 D7 z+ F( h　　}
" _4 s" g3 \- J" a# S8 C　　} * j/ v( r$ j: G! U2 l' g( E9 U

4 X  E$ x8 J5 u$ _
9 m. p/ a. u( }9 g; [　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。) `% M+ {" R% P4 l

+ u* M3 K: \  M% U　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：! u, w4 \( @. f, i- p+ l
# A6 k; `' a0 V( m" y1 f- G
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;/ _( X2 ?4 }$ |* ]7 l' W' X
* z7 x" ]" U; D5 A1 Y
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;  H2 a3 l' L  F3 _

% K6 h; i# h8 f; i! ^  G% @& i以下是引用片段：
/ h3 q2 Y" j4 b/ o# x- ~　　/* p, q is on the same of line l */ ' r4 Y3 L2 m) {# e  [- s; ]
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
# e: c; O8 W  `% o5 K' ~　　const vertex_t* p, ) g! _- I% N2 m7 a3 M3 C
　　const vertex_t* q)
& p$ V1 }' y1 W. w9 i# b' w　　{ 8 _" H4 c1 q4 v$ O/ h' i
　　double dx = l_end->x - l_start->x; * V$ X# u* l/ k/ B$ S
　　double dy = l_end->y - l_start->y;
, V9 s- a, |" |6 S; \2 E& W5 r　　double dx1= p->x - l_start->x;
* T. D! A! J1 O- f0 l6 C  S　　double dy1= p->y - l_start->y;
: M5 r" j. D% S. R4 V& Z+ h　　double dx2= q->x - l_end->x;
! L* c% g( V0 b* \　　double dy2= q->y - l_end->y; 0 M' I1 W8 ]  l4 u! f
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); * z+ z/ p0 d4 _" L
　　} ( n3 q+ }/ G, T" ~: S. m0 R8 S' w
　　/* 2 line segments (s1, s2) are intersect? */   h9 G- R- W+ n4 `- g
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, 8 C' O+ u  {& X; a$ {5 C
　　const vertex_t* s2_start, const vertex_t* s2_end) 3 c! i$ J  ?" Y2 W+ Q
　　{ + S6 ?' I+ L" v# d% g3 _
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && ; d1 `2 w5 B/ s9 K* K
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; 6 }1 V& ?. i. I" |3 O/ o8 G
　　}
. X6 f) f8 e. c+ H* z$ K2 h" j* T0 _9 q5 o. B

' {6 y7 b* f: U8 z. I* G/ `; J　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：# W  l0 ~5 {& M
" E, \$ I: l" A2 X1 U
以下是引用片段：; u* `% ?/ F: i/ Q4 K
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ 7 p$ V; ?# S- e4 k2 f
　　const vertex_t* v)   r. a. K; m' K
　　{
' x$ c6 v! w4 K; }9 l　　int i, j, k1, k2, c; & d$ l. f5 I& n( W
　　rect_t rc; 1 e( `$ N- J( M9 S2 I' A
　　vertex_t w;
! Q  V4 v6 `4 b8 m: Q5 J　　if (np < 3) + f8 |* w. k: ]$ L2 D' w, P
　　return 0; 5 ?) C* X" C; t0 T( m6 K
　　vertices_get_extent(vl, np, &rc);
) L6 N) X8 z2 ]: z% M1 C! j! h. l, [5 S　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
- E* o4 V+ O# [1 |' J- C% K/ ?　　return 0; 0 ?; Q1 m' |7 z; ^+ v
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
3 R. |/ B2 @3 \( B5 \　　w.x = rc.max_x + DBL_EPSILON;
" S1 t% x) f2 Y' {" x1 z　　w.y = v->y;
9 \, ]) M1 w' m7 ?% u! o　　c = 0; /* Intersection points counter */
! h3 T9 |: u) E& L. a* l　　for(i=0; i  
2 l$ b: ?' O3 F  }- T, R" i　　{
) }; q" s; {7 X; P  ]5 C9 w! Q! ~  h　　j = (i+1) % np;
% v$ n  E0 k. U0 d# e9 g2 I$ `　　if(is_intersect(vl+i, vl+j, v, &w)) ! X% L6 @; u' ^7 y- t$ D
　　{
+ q  x1 O6 \+ W4 a% F: J: r9 q　　C++; . F9 [% z2 @. H, S1 F5 i
　　} + Q' H, G1 r$ T. ?+ d
　　else if(vl.y==w.y) 3 b( O5 U3 G6 g# E/ P
　　{ % ?# g! G' B7 ^6 W; k
　　k1 = (np+i-1)%np;
& [/ M# N( t0 i$ L9 \# I. w　　while(k1!=i && vl[k1].y==w.y)
4 o0 v# o: C. E4 w8 w/ p. g　　k1 = (np+k1-1)%np; 2 m7 r- X9 w1 x- ?* _) G
　　k2 = (i+1)%np; % J- u' r% _5 U1 F: y3 V" y/ ?
　　while(k2!=i && vl[k2].y==w.y) ! g& ]/ `6 w! V( p' ^( I" ]$ P) ?8 |- X
　　k2 = (k2+1)%np;
; V0 x3 p3 L: w( q- d　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) ! F" X4 U2 S  z- _
　　C++; ! j; H+ d  r$ U
　　if(k2 <= i) 9 ~5 k+ l4 l0 f
　　break; : ?* {( p3 Z4 I: C
　　i = k2; ' K. U# p: c& {0 k( k  w6 t3 o
　　}
3 ?# a2 y4 P. R4 Z' {0 J+ y　　}   e' Q$ _) B. i) _( I/ b+ j
　　return c%2;
; \% t7 f' a$ N6 O4 C3 T　　}
2 v+ M- `! u. ~
2 S9 ?8 Y9 i2 B1 o3 k
. Z: _2 u+ g" x0 v4 g1 ]% k　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。