标题:
C语言中显示 点在多边形内 算法
[打印本页]
作者:
zw2004
时间:
2008-1-21 17:20
标题:
C语言中显示 点在多边形内 算法
本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
7 q1 y0 e7 }# P4 D3 i8 x
( s9 ?( c& `% l* x: s5 B
这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
; j w+ I" K* s- _9 t
" @1 y7 \$ [7 f
首先定义点结构如下：
& p+ l$ |! o) |4 `, m
& Z% c: r* q* o; u
以下是引用片段：
0 _9 G: _' N% ^7 N' A% t8 L, E7 k
/* Vertex structure */
^3 d' o- C F8 ^) X
typedef struct
# x, r1 V/ T3 L) ]
{
1 Q$ K2 p, U4 p. L; w
double x, y;
9 A0 p; d, ^( E& J5 U
} vertex_t;
# [) m& z3 Q4 s2 r' w$ Y
, ]5 Z1 Y! h7 r( k6 K4 U, _; a
0 U- j. [) j2 a1 V, U
本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
3 K& R- l+ K2 g; O
3 J% z1 n. g1 G& K6 O" f
以下是引用片段：
: G1 c8 Q' |8 K9 z
/* Vertex list structure – polygon */
! N: I8 `$ B2 h. G" F
typedef struct
7 C0 r% a/ u! f0 L
{
! v7 N& K% \; l6 e# m
int num_vertices; /* Number of vertices in list */
( u* b2 p" x9 V& Y9 y
vertex_t *vertex; /* Vertex array pointer */
- a" p; w! u' v" N+ N) P" V( d8 _- _
} vertexlist_t;
% n5 t1 E% O' n( n
4 H% L$ `" t! d) h0 K
$ ], X, ?" U# i
为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
/ L" _. i7 A! b
1 K/ F! y3 ?8 e: |, u5 {: d
以下是引用片段：
/ G& h3 Y( V3 |6 ?6 e" K, {
/* bounding rectangle type */
4 A% O2 g3 O% S& ~5 @
typedef struct
_' C9 K( W2 }) `
{
: h7 \( m* j' g+ j1 @( W4 z
double min_x, min_y, max_x, max_y;
& S7 b5 _' o9 `( W3 ~! i D9 q
} rect_t;
( g% R/ d) d, |$ z
/* gets extent of vertices */
% a3 [3 j0 X' m, K
void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
5 ]$ @0 B; j; f2 M3 T8 S# Z
rect_t* rc /* out extent*/ )
9 M2 D% J2 s9 l+ Z2 z" M5 }5 f
{
O6 t: ~$ \, J5 I8 B' X0 u
int i;
# r: u2 H( O1 F/ Q# t
if (np > 0){
. Q6 p8 \( i9 ~5 Q- Y+ P* H
rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
" t1 z* `* v" V) J* U
}else{
$ i0 H1 L, d) s6 Y4 s
rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
" P' C$ Q# m5 t
}
0 H2 {4 n" a7 Z8 {
for(i=1; i
% ~" \' ~7 q2 G5 q9 e ~$ E
{
- e; N( E; a1 U* Y
if(vl
.x < rc->min_x) rc->min_x = vl
.x;
9 F+ d' o* y8 j5 b- S
if(vl
.y < rc->min_y) rc->min_y = vl
.y;
7 a; c3 M* h4 g, M3 n k
if(vl
.x > rc->max_x) rc->max_x = vl
.x;
, Y2 Q$ Y5 X7 {" \ P
if(vl
.y > rc->max_y) rc->max_y = vl
.y;
3 E, p- U/ _/ u) G+ C
}
- u$ o! d4 W6 k! n4 e8 H
}
$ m/ F; R# k! l+ L8 \2 a! L
9 t7 s& M% |, q% x6 p7 l# h* I/ l3 e
" g* P1 T9 S8 x6 k0 j- j, L' m
当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
2 s) l% P* m" i3 D' u4 \) G/ D% l* T
) ]0 r1 m( O9 V- B# E) U/ |3 I
具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
* @ b3 x3 B, ~
8 q$ y7 q- [7 b+ l! f/ P" T
(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
5 b% h9 K7 \2 D+ u s
5 ?; D- Y* M" u
(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
. s, S/ d/ _ P
3 o7 h( H- b+ P. P( K4 Y
以下是引用片段：
: N0 [% _4 f/ v9 N0 F" `8 I0 p5 {
/* p, q is on the same of line l */
; Z& t$ o8 C# D. c$ z0 J7 T9 W
static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
: g" g9 B+ f" A l
const vertex_t* p,
+ \' X$ R1 b2 j7 j5 y4 s; |" e
const vertex_t* q)
5 e2 B! O8 B( c1 L. e
{
$ L. S) F3 X& c( ]6 b
double dx = l_end->x - l_start->x;
6 J& u; P7 ~' B; ?
double dy = l_end->y - l_start->y;
, a' Y6 `3 d1 N; d* \5 C
double dx1= p->x - l_start->x;
7 d8 b) v3 j( e) } E7 K, v
double dy1= p->y - l_start->y;
0 E2 |7 L$ }" ?4 Z6 W. E6 J
double dx2= q->x - l_end->x;
& R( o) o3 b( E$ ^
double dy2= q->y - l_end->y;
) E4 e6 X: a( e0 u. G
return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
+ d! x0 ~ v' L8 y' M% o
}
R" W3 e5 _3 o9 o8 J$ s
/* 2 line segments (s1, s2) are intersect? */
3 I, B' b5 I, Q" B- q
static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
0 _" S7 l' K7 I9 ~$ ~- D' h4 z" t' w
const vertex_t* s2_start, const vertex_t* s2_end)
+ }/ _) K; {. c. N. L1 l
{
0 X3 d5 u- S% i4 D8 F7 l* K5 h
return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
0 Y& X& s# N. ?& L
is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
/ N1 {" V% C3 t- c/ K
}
' P E% o7 c6 m9 B9 `: l( r
: R0 G$ i# n7 a# |& a* g5 V8 W1 X- N
% G9 K4 t. I1 v8 R8 a; n
下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
6 h' P* ~8 h: a/ D: {9 F! b8 ]6 X
! F- S8 T2 Z/ Y" q% Q0 |
以下是引用片段：
. O: S" Z( f# T
int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
4 i$ G* D' W. D, Z
const vertex_t* v)
\" E( c$ G# K: r
{
) N9 t) I7 e3 H1 o: `4 Z
int i, j, k1, k2, c;
7 C( g, j# k- m3 k) f# c$ A9 x
rect_t rc;
4 M4 A2 M9 @2 ~) ~# \! V
vertex_t w;
/ u$ D' [: k* S2 ]. \0 h3 j
if (np < 3)
: P' v, H# C, J5 [5 t
return 0;
' z/ ?- T4 \& Y) q$ m
vertices_get_extent(vl, np, &rc);
/ X3 I7 d, D2 D+ R8 ^! Z
if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
8 ]# c1 H( g; L' X+ K- ]
return 0;
1 Y( l( `" {$ ]1 C) U; D
/* Set a horizontal beam l(*v, w) from v to the ultra right */
# W2 q4 b) f2 s6 s
w.x = rc.max_x + DBL_EPSILON;
, V5 t0 Y" B8 f; z6 w" ^# E
w.y = v->y;
1 h' `6 ]$ n+ ^: N2 D
c = 0; /* Intersection points counter */
4 B( s0 s# x% q) m
for(i=0; i
) u) M" E8 h6 }' a3 }+ d. _" d
{
/ E8 P7 d2 K9 v: f
j = (i+1) % np;
0 M6 J9 o+ d" \8 y
if(is_intersect(vl+i, vl+j, v, &w))
% q p) S' b q$ P) Z2 \* S
{
3 {$ d8 h/ y5 }. [1 l; q" V E( V! r
C++;
) l7 o- b# L, W- [, ?. W5 X
}
( h- v' F* [7 z3 y/ ^4 m' i
else if(vl
.y==w.y)
( {+ M: p6 C! G& ]& P3 ^
{
/ }. Z/ v) x3 Q4 ^6 y* R9 _ E0 |
k1 = (np+i-1)%np;
9 \ h, p; `1 t) L3 W/ j
while(k1!=i && vl[k1].y==w.y)
: b: g, ~0 w0 F' J
k1 = (np+k1-1)%np;
5 F; j, c% V+ P" D4 v& q
k2 = (i+1)%np;
8 K# w+ J& m3 _* I6 C$ F! V* K
while(k2!=i && vl[k2].y==w.y)
% H+ j# L. A* l3 ~, p/ S3 X, N
k2 = (k2+1)%np;
3 ^* t5 P4 M" p
if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
9 b- _2 B1 p* y+ z
C++;
- o3 W8 z- A4 O$ W1 [
if(k2 <= i)
7 T& f+ ~8 ]# N) z
break;
+ K% `$ S) Z6 R+ ^7 R; |8 |
i = k2;
3 D9 p# R, q; a$ Q. u, v8 S" i2 D1 ]
}
1 Q# d- x: ]" O& @0 a N$ P2 [
}
7 a/ A. d; O) b! d, X0 m1 I ^( J9 P
return c%2;
9 g0 q W# W8 P( b) r: w; q
}
# p" |( |1 m+ o5 Z! N/ T* X
/ N" ~9 I# c0 E3 @& Q; u
: b* c+ R, L. P7 `( b
本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。
欢迎光临 捌玖网络工作室 (http://89w.org/)
Powered by Discuz! 7.2
