|
  
- UID
- 133
- 帖子
- 51
- 精华
- 1
- 积分
- 186
- 金币
- 55
- 威望
- 2
- 贡献
- 0
|
C语言中显示 点在多边形内 算法
本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。$ _. h* z+ k7 `' m9 p( P
9 l3 Z2 Q5 W0 v. a3 e8 W' Q8 b- K
这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
; ]$ v4 j' J7 s
) | [; X9 u9 T4 e- L4 Y 首先定义点结构如下：
9 d ] F( M4 |9 C4 n0 x
- W$ @0 ~; v2 M1 W) {3 a8 E以下是引用片段：
" l2 `- J" C2 h4 M0 x- } /* Vertex structure */
5 A& D: @: D, D! \ typedef struct
5 m1 H' P+ M8 D; A) m3 q) ~# @5 @ {
9 u2 ?; i, Z h double x, y;
; `5 Y+ H+ z- I0 i3 d/ I } vertex_t;
& J6 M+ s% k0 }% t) `5 Z, D4 u3 u. ]/ ` F1 ~- Z6 L4 F9 k. `5 X
$ ^* ?- B7 I" _
本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：, o1 x! R! o0 A0 B
8 K& l+ e% [, y9 ]; ?9 j
以下是引用片段：
* n& }2 B u0 b: z$ t9 r4 p /* Vertex list structure – polygon */
0 C& s I& E( c5 P8 ^9 v typedef struct 7 H; R _6 h' _
{
, P; O$ v& b0 @8 s# K$ j int num_vertices; /* Number of vertices in list */ . g h0 B9 J0 B e7 X. x* y: v3 C
vertex_t *vertex; /* Vertex array pointer */ ; m* Z) ]6 g, R7 R/ a
} vertexlist_t;
8 L2 ?: G# d$ x+ p: l) Z! r. p/ e$ t/ f, F i
( f, e( @9 y7 c f' s7 i F( E
为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：- z0 F- f% {9 B5 [- g
+ q! {8 h' E; |7 N以下是引用片段：2 l" `, g* ]$ U3 a4 {
/* bounding rectangle type */
3 b+ c; {6 T/ F( O: t5 M typedef struct
2 G, [( [$ d$ S+ ?8 m {
/ m& T% I, E, o0 f, c$ B1 \ double min_x, min_y, max_x, max_y;
* ]4 F# @0 ~- a. d } rect_t;
6 r6 S1 a! @( U /* gets extent of vertices */ 2 E. R1 f% r S [3 E
void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ # T5 Y! ?& X4 P6 L9 l! W+ w: u
rect_t* rc /* out extent*/ ) ( p, F3 G/ _0 O3 q* E0 d* X
{
, M8 c2 l0 {) |! h1 @ int i;
( q4 ~4 t: s7 [: _ [/ H5 L if (np > 0){ 0 S$ `( X& V$ x5 `6 B% W
rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; 0 F4 }4 [9 |( v. ]% @& t. K
}else{ ! r# E' o k, a% J9 I" j
rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
+ p' N# T d& l) B }
: J. e$ h: J$ D for(i=1; i + a% i# p" D" S4 q, y5 D
{ ' @; g- O# w3 Y8 [# a3 ]+ Y
if(vl.x < rc->min_x) rc->min_x = vl.x; ) e% g# z1 R" V/ h; @
if(vl.y < rc->min_y) rc->min_y = vl.y; 2 l+ ~* k% x" }8 Q6 x% m# Z/ R
if(vl.x > rc->max_x) rc->max_x = vl.x;
4 o. s5 A0 T# V6 F8 G, W! E5 h if(vl.y > rc->max_y) rc->max_y = vl.y; * Z. C3 y, ~% s! d
} 5 y, l$ {9 V6 e( X
}
5 v/ y7 x* Z* `- u5 I+ U/ K0 c* k# v+ L1 ^8 t+ P
1 S- ]( p! f3 ^ [ 当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。! M$ m n9 v+ Y* K" |0 K
( g3 U% h( l8 a% f 具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：( f7 |6 ?7 `5 y* s8 o
. c3 M3 {; e( s% R (1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
* `3 M, Z3 b& R# d; a6 `
" C8 d; r- \( Z (2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
1 |2 b {8 |; Y+ v' d7 s
, F% l( t" S+ E! K6 m6 ]以下是引用片段：
7 i/ W( _9 O$ y% Y6 d: N /* p, q is on the same of line l */
8 l& b4 v# z; J; N# O static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ 4 ~3 a e3 L: A6 \/ g$ v( i
const vertex_t* p,
' a- [* l4 R7 j3 s& D1 m: _ const vertex_t* q) + f% R |( U. g9 h# D2 A
{ : t6 z& U( c- \" L7 t n$ z
double dx = l_end->x - l_start->x; 8 b$ q" w, g7 j
double dy = l_end->y - l_start->y;
) ?- z! a( E+ |0 l) } double dx1= p->x - l_start->x;
' G, C$ V& i9 S. _6 R) G4 S double dy1= p->y - l_start->y; # F+ {7 r7 x1 C3 F7 H( g6 _* H j1 Q
double dx2= q->x - l_end->x; 5 c& c: w' D! a0 D9 W
double dy2= q->y - l_end->y; " C: B8 v7 }. x" g, ~+ ]7 r
return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
$ ~4 Y& J! ]# ~. @5 B; a+ a4 S }
, r- |# p- h% _& B% V /* 2 line segments (s1, s2) are intersect? */
, }( A# J# m7 ]1 @1 ~( Z; b/ V static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
& a/ }# y8 j+ p" O" s const vertex_t* s2_start, const vertex_t* s2_end)
) [4 a. i8 T! c7 h; \ {
1 _( f4 A% G, z* g; K& r return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
. l- t" v- o5 U( }& M! \9 J- } is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; - s1 J: K/ g5 J8 \& u
} ; F4 v' \1 M" y8 a
. X0 u# A; q9 V9 x/ m3 Q
, o$ ~- S; a! d0 ~, f* u
下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
% j( l1 g. t/ A* A
' ?, I$ i3 z5 Q- |+ @% K" u以下是引用片段：
+ Y3 r( X2 H2 n. M int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
2 R7 Q; X- a0 P. s9 m const vertex_t* v)
" P6 n3 V- E. S: U$ @8 q/ w' V { , v: D( x) }, |: L7 y$ T$ ~. k
int i, j, k1, k2, c; : y2 B9 Q4 ? T# K7 U; H" b
rect_t rc;
y c; h- {& X; S( w vertex_t w; . b; N( r! e$ G1 w G" v
if (np < 3)
: Q# Y# P$ q/ [, g5 F return 0; 9 R" n, C5 A8 w1 \3 ?
vertices_get_extent(vl, np, &rc); 0 v. p2 H% M S% Z4 E# S
if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) $ H7 i6 u8 j! C+ x5 X# C/ C
return 0;
; j2 w/ C& K, \; W1 ^ /* Set a horizontal beam l(*v, w) from v to the ultra right */
. M0 c% O8 k2 D0 W; w* Y# g. }4 Y- N w.x = rc.max_x + DBL_EPSILON; : T. w7 ~5 I' l/ Y! s# X/ L6 g$ Y
w.y = v->y;
& }3 a% c0 }5 c6 |& F [) h# I c = 0; /* Intersection points counter */ 5 |8 V: U8 R0 G0 k4 a: A
for(i=0; i
! Y1 R% C- k9 ?9 U5 H6 ^$ W6 R& A { 7 d3 `" _" G' Q
j = (i+1) % np; 3 X' `, n: F7 Y, k
if(is_intersect(vl+i, vl+j, v, &w)) ! c: _. k: H6 l# i
{ T u- G9 Z' X2 P6 D3 O
C++;
6 j, p6 J; x# b: o3 C# o } 8 e" Z" V1 T1 |7 l, x7 v
else if(vl.y==w.y) 1 s4 ^! X, K1 w1 e2 v5 M+ p. d
{ 1 k) _: ~# B8 ?2 _: f! D% M
k1 = (np+i-1)%np;
& K" X Q4 g& F; F: `$ H while(k1!=i && vl[k1].y==w.y)
' q: I( _4 f g& d1 i+ u8 H k1 = (np+k1-1)%np;
+ @' `% T. Y, N C \ k2 = (i+1)%np;
/ \! O( d$ N8 w* U; a while(k2!=i && vl[k2].y==w.y)
3 `: M" N" f0 M k2 = (k2+1)%np; & E8 s1 ^3 m8 H8 ~
if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
2 J& k$ ^$ c8 ]* Y4 y/ L C++;
) e3 y6 _ N: b& @, w" k0 ? if(k2 <= i)
, E& t% b0 |* Q* ~! ?5 W break; Z, R: F7 l# \5 O! Q( ^+ T
i = k2;
) v, a# p$ s8 b9 j6 n7 d& I }
2 ^7 v) [0 g* ~% \% @' b& J! J5 { } ' c* M# X6 }( V5 e* T) @: D
return c%2; * w+ N3 S4 w7 a2 E
}
$ T/ e0 e; m% b" }7 a- X0 Q7 d
1 @9 P7 b3 N" F/ }9 q1 v2 Y
: r2 {; \/ ~+ {7 \ 本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。 |
|
发表于 2008-1-21 17:20
|