|
  
- UID
- 133
- 帖子
- 51
- 精华
- 1
- 积分
- 186
- 金币
- 55
- 威望
- 2
- 贡献
- 0
|
C语言中显示 点在多边形内 算法
本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。: W b0 x* E4 \0 c/ |% U! _
; ]# c- k2 f4 i( _) R7 U: p) E
这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。+ `8 }; o3 V- f! ~1 t% Z k
/ ]& y+ Z; l4 Z, N/ Z: Q
首先定义点结构如下：' Z7 E, n+ w* X) }0 M
0 t/ ~$ W. m+ D6 M) t ^以下是引用片段：
+ T% H) X4 m: a, ]% U /* Vertex structure */ 4 V4 B+ F# H! \: J$ |, ]( b
typedef struct 1 Y/ W7 r7 x5 C
{ ! O4 c; V3 B4 f' T% [. f
double x, y;
6 k( H/ m/ ~3 ~/ |+ P# w, k } vertex_t;
, t& m3 {+ }7 P9 f
9 W6 [- `6 Y+ z W3 ?+ B2 ] e" g! _, c4 L2 w* f4 u
本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：+ m- M5 K3 n0 Y+ l6 b# V. A3 b
- F/ P; K* `8 w) ]
以下是引用片段：
' _: M& w( S- C- y# d /* Vertex list structure – polygon */
M+ {1 }( D q8 T$ B8 M typedef struct
1 X. B. T6 }+ e- {7 ^4 m! { {
) t5 M* T* M/ e2 ^" ?8 E int num_vertices; /* Number of vertices in list */ " M+ D( Q; E8 ]. | U( B/ t
vertex_t *vertex; /* Vertex array pointer */
! b0 n0 Z4 w' u1 \; g; I! P5 n/ e* r } vertexlist_t; * q" c" _0 t; _; k
( ^7 X: S, V' O
3 y% ?+ I' ~( J- Z9 |* ]8 V' W 为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
1 R3 Y3 c# U( Z5 k2 u
. _ ?7 i& A' e$ z5 y3 m以下是引用片段：* l; X; t3 r) z7 J% D
/* bounding rectangle type */ 4 c# H% Q1 R/ @1 N- n
typedef struct
. N/ J! x: |+ @7 w% x {
2 N' V9 B/ S! Y4 o3 f" R double min_x, min_y, max_x, max_y;
1 K7 b" G8 B( ~5 ?9 |3 M! Y } rect_t; 9 Z. ], _, D( m9 J* _6 I: w
/* gets extent of vertices */ + M8 w+ K- x. q1 o% P
void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ o9 X% z7 Z- O3 l
rect_t* rc /* out extent*/ ) / d3 @( q5 H! ?! z; l- F
{
9 V& y5 Y" f* r. p2 { int i; - J$ J9 S; i$ h) T
if (np > 0){ / N, ^3 B1 Y$ @6 G9 k
rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; * r6 ~# C6 ]1 j1 e- n) W
}else{ 8 [& N6 W* _4 r" T& H4 T& _# B' }/ G
rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
$ H8 @: d+ h4 y8 a }
* s7 u- n. }8 }/ Z6 y for(i=1; i 4 [' m1 [/ I2 n
{ & ]( T- K3 E5 S7 [3 o$ ~
if(vl.x < rc->min_x) rc->min_x = vl.x; " Z' a) {, @$ Q/ p6 T/ p( m5 F
if(vl.y < rc->min_y) rc->min_y = vl.y; $ c# p/ l, x/ b0 F; B
if(vl.x > rc->max_x) rc->max_x = vl.x; 6 G, P& }' q$ P
if(vl.y > rc->max_y) rc->max_y = vl.y;
0 q# r2 d6 j1 D7 X! R# Q& C7 w }
3 f+ Y/ M" p4 z& g# ?( k } 2 G9 E& \- u+ ~, O7 A
2 L/ N; |$ V. m1 m. G" ?- Q8 N3 v4 F! x
当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。7 k+ s" L! w i: J
3 a& `8 U9 k% K* p. S- I& S* R4 Y1 v/ l
具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：0 b9 X3 @0 D% ^
6 ]1 Y3 c) u+ H4 n9 G. A1 p6 _ _9 H (1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
( H# @: _. C6 R1 Q1 ^) i$ I
1 i8 F7 M" }8 c ~# Z5 y" Q (2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;; ~* N( M1 D$ n8 z' x% W: }
, c4 E4 G6 |3 c
以下是引用片段：
' v' ^' }" ?* q8 r- A, f' t /* p, q is on the same of line l */
% a# P1 S" _! ]+ ]2 i+ G$ T7 T8 W static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ + ~! @5 x6 e% L8 r( Y# A
const vertex_t* p,
$ ]% K2 H8 |: Z2 T, G+ q1 V. X7 k' Y# ] const vertex_t* q)
: N& I5 H9 T% ^! E, m { 5 t4 J' V2 G) \6 i) F
double dx = l_end->x - l_start->x; & B; ~- Y5 R+ F, K+ U' B
double dy = l_end->y - l_start->y; 1 A# n. S" u4 b6 ?: W
double dx1= p->x - l_start->x;
3 t6 x) r+ _& w* d double dy1= p->y - l_start->y; # e4 x7 o3 W F0 u& ?
double dx2= q->x - l_end->x; ! A) U s' t- `# z `+ Z, _) v
double dy2= q->y - l_end->y; 6 s- B: a! p0 c7 C) L
return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
# n, v" ^1 Q+ B. ~# J0 W$ A } ( t% `+ V) T- q8 s! Z% I7 q
/* 2 line segments (s1, s2) are intersect? */ ' w( d# I5 Y. w7 d
static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
. x$ O$ d0 I+ |& \ const vertex_t* s2_start, const vertex_t* s2_end) + E9 Y) }2 W& P; w& X8 O
{ / R# g! b- i: j6 O/ ?" x
return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && * ]# G/ r0 R! l6 g8 _, m
is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
; u- a8 M- u& n8 E } " S# ^6 y4 o3 g) z- g
; A- S8 I5 T$ v+ e. S1 k
. d% t, w. M$ T3 \. m4 x( j( \ 下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
) {. G! C0 ?1 D/ k5 t \5 b
R; C* ?& G5 `% o& P& c以下是引用片段：& P8 F8 d! [6 D: c8 C
int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ ! A2 r8 Q4 t4 ~% m G
const vertex_t* v) ' Y6 m5 d- d; P6 i
{
# v) U. |( d' W1 v int i, j, k1, k2, c;
) T- D8 w5 S4 Z1 e0 y rect_t rc; 8 m% g3 J$ a: ~% o- P8 J* P3 I
vertex_t w;
+ {( g' a- t! M/ f if (np < 3)
! E! C4 O; S. Q, H return 0;
2 E+ R7 Z& o9 c9 ` vertices_get_extent(vl, np, &rc);
3 }) I+ \: t& ?; |' w; R if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
& d! d7 |# \& p I4 i: G; r" F return 0;
0 u# G% N: {; U1 N /* Set a horizontal beam l(*v, w) from v to the ultra right */
; z0 o- p# V3 v' y5 l w.x = rc.max_x + DBL_EPSILON;
T8 L( I3 _ Y3 I# U& N w.y = v->y; 2 I a+ G( X$ Y# ~0 Y
c = 0; /* Intersection points counter */ ' t# I4 t2 m: F
for(i=0; i . T9 Z5 B, b; B2 v
{
# S9 Z7 [) \$ D j = (i+1) % np;
9 ?* S m2 C; N2 c if(is_intersect(vl+i, vl+j, v, &w))
" j! z5 @) J& [2 q: d; L { / g# t( z+ `5 T/ m! C9 i2 Z
C++; 5 s) a4 R$ j+ G. q: O
}
/ e# d& X8 r4 R P3 S. z1 v. Z& h else if(vl.y==w.y)
' u; C, O8 I& E, o$ h6 x3 A {
% _5 P* y8 M5 D4 m+ @. l k1 = (np+i-1)%np;
3 `" X1 \0 H0 F2 J while(k1!=i && vl[k1].y==w.y) ( O; A0 q, X- Z: S7 q
k1 = (np+k1-1)%np; $ a8 e' w6 _- U8 @. ]
k2 = (i+1)%np; ) `" y8 H+ Z! W5 K( m% z- J7 Y
while(k2!=i && vl[k2].y==w.y) 4 S d. P, H( i+ Z4 d6 W# n
k2 = (k2+1)%np; / L+ J F1 T) f" [2 ]7 l4 @/ @: H
if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
( _" n* k+ N. ~; t$ z C++;
. e: s) c: ?. D- w# }) E if(k2 <= i) 3 {9 I$ P: O1 z: w' H* h
break; : V9 J$ M# e+ L$ z/ x# d
i = k2;
* X# y S) ^+ E6 q* }, V9 ^& E E } + u+ f9 ~2 G6 Y5 ?1 d
} m: M! u8 f. Q* I. x( @
return c%2;
" T4 @: M2 r4 u. c } 8 w9 I8 d6 g. R
% S$ S3 Y4 ^* Q3 \1 J7 }8 @+ h8 ^8 I$ h4 t/ w& o: O1 C# v; |2 X
本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。 |
|
发表于 2008-1-21 17:20
|