|
  
- UID
- 133
- 帖子
- 51
- 精华
- 1
- 积分
- 186
- 金币
- 55
- 威望
- 2
- 贡献
- 0

|
C语言中显示 点在多边形内 算法
本文是采用射线法判断点是否在多边形内的C语言程序。多年前,我自己实现了这样一个算法。但是随着时间的推移,我决定重写这个代码。参考周培德的《计算几何》一书,结合我的实践和经验,我相信,在这个算法的实现上,这是你迄今为止遇到的最优的代码。
# ^) i5 j) h# K) e( N& |% A- I- I
& \1 z2 q4 `/ M- u F' } 这是个C语言的小算法的实现程序,本来不想放到这里。可是,当我自己要实现这样一个算法的时候,想在网上找个现成的,考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心,所以,决定重新写一个,并把它放到这里,以飨读者。也增加一下BLOG的点击量。
. F4 h7 F+ _1 B# | c5 Z
# v6 W/ t& B8 z" \# H 首先定义点结构如下:
$ U* p2 H) w+ P% ~' G4 H, k: W- t) D/ _, Z. L) ^
以下是引用片段:5 n* W% P6 T4 i5 O
/* Vertex structure */
' a' U1 D3 z; {, | typedef struct
* Q% Y4 d; [5 r( v6 h1 X) C { [) E/ S1 s$ D
double x, y;
6 r. ]+ ]; J2 V6 r& i } vertex_t;
. Y; h: F9 N& p. f8 S( l. F8 V! K5 ]. N3 _" G# W
5 d& ]/ K. e, i; R; k e( L 本算法里所指的多边形,是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的,包括多条边在一条绝对直线上。因此,定义多边形结构如下:6 o4 q$ T* `( R% U* X+ B
" j9 Q/ }3 f7 s$ g6 c' P以下是引用片段:$ z. I' q- q$ q/ j! e6 z# m% _
/* Vertex list structure – polygon */
% u' {1 |# l. l9 j, \ typedef struct
# s5 a# l2 \4 W' a0 b' x, X( n4 n {
8 ?# G0 J' w6 P7 b0 z; v int num_vertices; /* Number of vertices in list */
, u3 D1 W8 M v8 X1 P, d! c vertex_t *vertex; /* Vertex array pointer */
2 y4 Z+ T& M; @ v } vertexlist_t; 6 X- E6 E- N& E- n
3 ]) ^7 V0 |9 H
9 o* }% `- e x- E5 w% N 为加快判别速度,首先计算多边形的外包矩形(rect_t),判断点是否落在外包矩形内,只有满足落在外包矩形内的条件的点,才进入下一步的计算。为此,引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent,代码如下:: t+ e/ l: _: }# d# ?
7 E% B9 {: s: F4 X, v. J C& }4 t以下是引用片段:
) }# {2 I4 U0 p /* bounding rectangle type */
4 O/ J2 _) w7 U7 \4 H4 D% { typedef struct
0 l- S- c* K8 y) f0 Q& d$ J% f { / H$ \' D0 M) X9 P
double min_x, min_y, max_x, max_y; $ z; E9 _ E, P6 H: K) d
} rect_t; 7 A1 R$ s" D7 g j1 L! v+ v: y% p
/* gets extent of vertices */ 0 X4 ~' I8 B8 c
void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ 5 u+ u& G4 P" _: R
rect_t* rc /* out extent*/ )
( ?$ L: ? z: Q( g* N5 ~ { 8 e% }( a7 n$ s! g
int i;
. f$ R8 g& m- A5 J' Q0 Z$ y if (np > 0){
/ \2 V5 L G* h rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
- g* s/ q- a0 g7 _0 I! e }else{
! m& @+ Y. G J# j: _- t6 Y rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
$ R! W0 R1 X8 q# e+ x% ?1 B1 V } - m# t9 U7 X" `7 S
for(i=1; i
8 q! @; K" [6 p6 [8 G, R/ B2 [ { # {4 A# L1 e5 K$ p0 `
if(vl.x < rc->min_x) rc->min_x = vl.x; - d6 N7 f9 o9 e1 C
if(vl.y < rc->min_y) rc->min_y = vl.y;
$ Z9 l# f) W( o0 v( b; }+ J! d- n if(vl.x > rc->max_x) rc->max_x = vl.x; ; g o+ ~7 u# d+ Q3 F+ @
if(vl.y > rc->max_y) rc->max_y = vl.y;
( K% Z( a5 J0 M7 R }
' H% [8 N0 @% V }
/ M4 k0 |6 ~6 g4 \) l. ]0 I% e, J( j w( \0 c" b* G0 ]+ w a) x
j2 h+ L k+ ~ v/ ^7 { 当点满足落在多边形外包矩形内的条件,要进一步判断点(v)是否在多边形(vl:np)内。本程序采用射线法,由待测试点(v)水平引出一条射线B(v,w),计算B与vl边线的交点数目,记为c,根据奇内偶外原则(c为奇数说明v在vl内,否则v不在vl内)判断点是否在多边形内。: m1 z* M: K' e
7 c/ @; {& |; n% a' J; S 具体原理就不多说。为计算线段间是否存在交点,引入下面的函数:$ n9 m0 w% I. r* y! z: k
0 n) B0 \8 W+ Z* I( a' ]7 R
(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start,l_end)的同侧;+ V+ w; k: n3 N2 R/ Q
' w y* a: K) ~8 h1 m (2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
* T; d9 Y5 \, P6 u; U5 ~9 G9 D7 ~, N' E) H3 }2 n4 `
以下是引用片段:( |, X, f) Q) Z+ r
/* p, q is on the same of line l */ " f4 @( G9 Y g# z& X7 U: n
static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ . ]& P# t5 L9 K/ \
const vertex_t* p, / V( [ W: q% ^' I) G3 G
const vertex_t* q) & l1 Q! \# d/ \2 u, n* c
{ + k+ z6 A. N$ d- V2 X: m9 a7 p
double dx = l_end->x - l_start->x; ( l7 b7 g& L( v8 q
double dy = l_end->y - l_start->y;
1 E5 C2 s, G! e1 ^ a- h double dx1= p->x - l_start->x; , C) r: ^7 H: g
double dy1= p->y - l_start->y;
@: E9 ^7 T b v1 H double dx2= q->x - l_end->x;
# o4 l; w# N# b' l8 j double dy2= q->y - l_end->y;
* f3 N2 r; \/ S( L& A# n% t return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
, }" E4 z0 D4 E- f- d } ! Q; R4 j. h+ v8 ~5 u. h9 o s3 m. o
/* 2 line segments (s1, s2) are intersect? */ 2 p' w+ D7 v8 U" f* P7 `
static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
' k F; A" T# g! W1 b const vertex_t* s2_start, const vertex_t* s2_end)
9 k1 k3 p+ V) Y& U {
3 y" v' V6 ^! ?" y* r$ E& I return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
; K' M7 @9 v! H, M' ` is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
$ T3 a- y' ~ | X: W8 Z' o }
! }9 x% f- f# r* ^1 L& W, p7 R. E/ z2 m9 x
% ~1 b" p+ B* n7 f( W2 X7 T+ ~
下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl:np)内的程序:+ M6 A8 Y: d3 U& g# @7 c
* u8 E. U" k3 Y* N6 e2 O2 O
以下是引用片段:
3 A g& o: S1 t& a1 K. u3 K int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ # ?+ }) H7 L# `1 l$ U- E7 X
const vertex_t* v)
- K$ F( S8 @7 B9 o" b4 w" X6 D* ] { . H& k3 u" u0 [" Z7 L p. _, F$ p
int i, j, k1, k2, c;
D7 u0 Z% K$ q; C0 _+ Q rect_t rc; / p6 A- I1 C, B' Q6 a- X! C
vertex_t w;
8 \ P6 H" S( r |( P if (np < 3) 7 M s1 E4 l2 m# ?$ s
return 0; ( J' S7 @4 i' n' |& T) S
vertices_get_extent(vl, np, &rc); & I) T w; a+ i9 i- x
if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) , n5 b2 Y/ d1 D, T( o6 Q9 |+ D& _
return 0;
/ S# q0 ?' N+ K /* Set a horizontal beam l(*v, w) from v to the ultra right */ * S+ _( R+ X p% P
w.x = rc.max_x + DBL_EPSILON;
?, W2 D. G0 C w.y = v->y;
# t9 Q$ m2 I3 N7 w7 p c = 0; /* Intersection points counter */
! X& _5 U5 J7 g for(i=0; i
# s: ^* O' q) ^: ^# S { 3 N. @! x/ c- g D
j = (i+1) % np; 6 c9 i) c Y; Y0 `( E0 w
if(is_intersect(vl+i, vl+j, v, &w))
1 J6 t v6 C1 z9 o3 s0 d4 [* P { $ ?2 @1 u! s) b4 v0 U
C++; 3 b: A6 q$ ]9 r2 }# \
}
) e0 z# f, m" a else if(vl.y==w.y)
6 C+ w' E" O! R( d2 t$ S# s5 q& a { " S& @' y- R6 _
k1 = (np+i-1)%np; 8 G; c1 _8 J5 M. g2 X! R7 P
while(k1!=i && vl[k1].y==w.y) 3 ?8 z- z0 _/ c- p* Y# {2 q2 b
k1 = (np+k1-1)%np;
2 @1 A& j1 B u$ G# c k2 = (i+1)%np;
; j1 w% ?2 E" C# v5 E while(k2!=i && vl[k2].y==w.y)
$ y3 l6 Z) ]3 y3 H k2 = (k2+1)%np;
9 c, M/ `) Z, R5 K6 Y% A if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) 6 d; T `1 x) n
C++;
/ ^& ^( u4 D5 D! W if(k2 <= i)
$ E8 |) H% m8 s' g5 x' N break;
5 Z0 W6 z$ b! n! t/ O i = k2;
2 P9 I4 i$ O! c9 p+ B" z4 w8 d } 0 Z+ m: ~0 V, {( C- B
} 0 H+ c/ @" k6 v
return c%2; N' y" L( `! |( U4 V/ K0 g
}
% p; s; [! C4 ?9 s/ t: ^+ R$ J5 i- Q* F% o& v& T6 \3 g
k& D2 Z( N2 A9 d l5 _. x7 Z
本想配些插图说明问题,但是,CSDN的文章里放图片我还没用过。以后再试吧!实践证明,本程序算法的适应性极强。但是,对于点正好落在多边形边上的极端情形,有可能得出2种不同的结果。 |
|