标题:
C语言中显示 点在多边形内 算法
[打印本页]
作者:
zw2004
时间:
2008-1-21 17:20
标题:
C语言中显示 点在多边形内 算法
本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
/ ?& x& H# e1 a3 L9 L( G1 ^* Y
2 X3 Z. |8 a; y* I* ~% u
这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
/ r- ?, H) B! i% l4 {& D
7 s8 ^$ w! J6 R; I- l8 p5 @
首先定义点结构如下：
% C' Z: {: k6 }; d6 Z
4 m. @, w, Q- @6 O; W- j
以下是引用片段：
2 z5 W7 H. c& ]7 z( O% U O- o
/* Vertex structure */
6 W0 e- G8 v; F, u1 {
typedef struct
8 z; W% A0 u4 U2 I9 T
{
) g; n, v5 I4 I1 \1 y& t5 J
double x, y;
~$ D: w: n5 E9 |" f
} vertex_t;
, S m1 q8 ^3 y! _. _; D: A) I
" J4 n' O7 a8 j* n4 ?
" u D+ O$ g: k/ | M8 \- X
本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
; t( N E' @. @5 q0 _4 j
( W( y: z8 N% n$ C% I
以下是引用片段：
; \5 L) {% f$ r6 c) Y
/* Vertex list structure – polygon */
% P; [/ |/ p; b9 k7 F% C
typedef struct
$ d. [4 Y8 d# E4 s2 u
{
+ K' L9 r: K4 K( k
int num_vertices; /* Number of vertices in list */
; c' s4 F7 u; p$ W6 J: c" F( m
vertex_t *vertex; /* Vertex array pointer */
& e0 }# \/ L. e; n2 e- N9 l1 ~
} vertexlist_t;
7 c9 A# p& X7 h6 _5 i$ Z
8 m! K: w" y# S# b/ j
# f3 n% z! Y/ f& H
为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
2 Y7 T- r! s4 W/ g8 A" ?( Q3 a. M
) \+ H* K' z- L9 D
以下是引用片段：
) ~$ k" d8 A, E# s/ e
/* bounding rectangle type */
2 L" U" C+ X; @" a4 K
typedef struct
+ a1 I1 _7 V% h7 S) N
{
# g( R6 Y7 H s6 i$ Z2 {
double min_x, min_y, max_x, max_y;
$ v1 l2 X7 O$ x1 |# q
} rect_t;
( N! E5 I# m# j2 G: \, }& A
/* gets extent of vertices */
9 N1 F' i- z- X- |- t
void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
* P6 j. H& y) `! f
rect_t* rc /* out extent*/ )
! e- P! x( u9 d/ y4 f
{
$ E' `3 x/ k4 ~0 O5 C w6 f* E
int i;
9 E V/ b* Y. A3 c3 g: w' x
if (np > 0){
Q7 h. ?4 h# ?, v( e
rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
1 V$ h2 f6 O) R3 ]& y
}else{
" b& V( `% [. y& }
rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
$ A6 f: J( }( x& j& s9 j3 a
}
. M' A8 b& T* T* m6 f" L( d! g
for(i=1; i
) p' E: u# l% ?. C6 ^% H5 X
{
9 r m7 ? f* x5 w* ]& G
if(vl
.x < rc->min_x) rc->min_x = vl
.x;
! `: T% x7 L+ f+ C y+ v o
if(vl
.y < rc->min_y) rc->min_y = vl
.y;
3 f* F" Q- A1 `* G
if(vl
.x > rc->max_x) rc->max_x = vl
.x;
" O8 p% m3 b$ ]! v- t S9 Z* O
if(vl
.y > rc->max_y) rc->max_y = vl
.y;
& H u9 q& h# `! P# C
}
: {' D2 K( w5 N5 K& b
}
- W3 c; ~- i" Q2 [9 P5 F. H1 \
* M: ?4 U) g H" ]
+ b9 |9 @* ` J9 z
当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
7 D# Y: M4 \8 m5 w8 _7 K% }" \ O7 {
% n5 V0 z+ ^( u! f, H- q# V y2 h
具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
7 j8 u. o' I$ _4 a
& i! i' n, q) U% Y2 O5 {5 A
(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
8 Q1 [1 G9 {3 `% ]
* M8 r7 u+ P# l5 g
(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
/ X( S6 k; E6 l' W. b& h
" V' [1 \* k* n' K
以下是引用片段：
$ L4 l6 Q9 E; N0 G$ ?/ m4 u
/* p, q is on the same of line l */
$ Z; z( X5 _/ Z+ I
static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
9 s0 q' S9 `; j: A3 r# _5 c
const vertex_t* p,
/ l* a. q0 E3 o8 @0 _- C
const vertex_t* q)
1 J2 o+ P% J6 e# l0 h+ r; t$ ~
{
- | v5 s/ E% U7 H) t0 r
double dx = l_end->x - l_start->x;
) A' [0 k4 @4 x- l m& r
double dy = l_end->y - l_start->y;
U3 C/ o& Z; l7 o
double dx1= p->x - l_start->x;
0 `( D& I; G5 ?1 K+ ~3 E/ c* {7 q
double dy1= p->y - l_start->y;
" _/ s4 y, D0 P. M
double dx2= q->x - l_end->x;
* z) N x1 P, Y/ @8 X8 V1 s( J( r
double dy2= q->y - l_end->y;
5 Y0 n% e& {4 m4 U* P/ g
return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
& o. x7 p9 [ k( O" }4 Y- w5 Y% A
}
# d5 w1 s0 T) j
/* 2 line segments (s1, s2) are intersect? */
s8 C4 f- v% w+ G, T% \
static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
# x! E5 N# a/ |, P* v5 D5 W: [
const vertex_t* s2_start, const vertex_t* s2_end)
* b; ]( i4 y5 n/ a: X- h
{
! o5 n" m: J1 @0 [ a0 u# e
return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
/ Q' H+ m( v$ Z
is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
9 k* m, h1 I ^ o4 g
}
, ~4 P5 \8 `& h b. B
% ?1 R, g$ U) T& J& t- d( t5 p- _
7 v7 M* d) @& k" S; y
下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
6 d4 Y( n" |8 K. M# ?; V: c: c$ a
& A: T* X- Q' J: f$ l$ M
以下是引用片段：
, t' [ ?+ b, k; ]4 o
int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
* D$ ?7 b; w* d0 O y
const vertex_t* v)
1 ?( m5 J# U0 J l8 R
{
# i& O: ~6 d7 q/ b$ T
int i, j, k1, k2, c;
7 P0 T! Z* m0 q+ O
rect_t rc;
. F0 r. b( o) K
vertex_t w;
) h" Z7 F8 r8 H0 M+ O- B9 p5 a
if (np < 3)
' j6 e7 p% m5 a/ {1 T
return 0;
# r# P) M% H$ b# `7 B7 S8 Z- N
vertices_get_extent(vl, np, &rc);
: h4 W0 E2 u9 U
if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
- A8 w3 r$ W1 r6 ^) Y% K9 O5 I3 ^- X
return 0;
0 W4 ]/ D9 w* x! \4 w0 `: ]. U
/* Set a horizontal beam l(*v, w) from v to the ultra right */
5 v, {+ L' S: x, j$ y
w.x = rc.max_x + DBL_EPSILON;
% U8 d/ Y) @& y1 u# ]! u
w.y = v->y;
$ _9 n' a) y2 q' @' Q0 W; ^
c = 0; /* Intersection points counter */
) J1 @2 r: R# p: z9 ^
for(i=0; i
2 e( [' M3 I. l' o1 K
{
: }8 c% Y- E* ?9 k
j = (i+1) % np;
9 a4 _# k5 w# R w
if(is_intersect(vl+i, vl+j, v, &w))
8 _0 K7 O7 y' H2 N- z) W: `
{
k4 h# h% M& W2 O% h
C++;
) k9 m% N& O2 l2 b8 o9 y
}
* [- b* [! b+ W; E* ^1 h
else if(vl
.y==w.y)
. q& ~' i% T9 l$ d, H. _
{
9 p/ d, Y; B$ z6 ]
k1 = (np+i-1)%np;
( E. G1 M/ ?( O' H' T
while(k1!=i && vl[k1].y==w.y)
0 t" b1 v6 d; h( g9 r9 K3 w0 U
k1 = (np+k1-1)%np;
3 A$ I1 E b! r: M3 P4 e: M3 [
k2 = (i+1)%np;
1 u4 c( O$ X2 F/ c
while(k2!=i && vl[k2].y==w.y)
5 q% h3 r4 R3 D B/ |
k2 = (k2+1)%np;
! F5 E( \% ^/ L* k3 F' m
if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
' b# Z% F y( V
C++;
: T5 F0 T9 P' U* t: x
if(k2 <= i)
; P0 L" ?# s4 u
break;
' J% I( L( n1 H, A
i = k2;
- p% Y! `1 }5 k6 O- ?! H
}
! _) v/ P) p4 y( D
}
; a: g% I! t( a
return c%2;
! B3 x: M( _0 }
}
7 n4 h, G2 O% G2 |
4 a7 K; e9 t. a: r. P0 y, H( l5 V. t
- ~ o' `2 {# z0 d
本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。
欢迎光临 捌玖网络工作室 (http://www.89w.org/)
Powered by Discuz! 7.2
