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

|
C语言中显示 点在多边形内 算法
本文是采用射线法判断点是否在多边形内的C语言程序。多年前,我自己实现了这样一个算法。但是随着时间的推移,我决定重写这个代码。参考周培德的《计算几何》一书,结合我的实践和经验,我相信,在这个算法的实现上,这是你迄今为止遇到的最优的代码。
& s/ \/ ?( E5 X2 J2 \4 D' e0 q" E' b1 t! g' l" [ i; U3 G
这是个C语言的小算法的实现程序,本来不想放到这里。可是,当我自己要实现这样一个算法的时候,想在网上找个现成的,考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心,所以,决定重新写一个,并把它放到这里,以飨读者。也增加一下BLOG的点击量。
0 K5 |7 E: B# ~2 p" V3 r0 F, M( {8 J |! z5 h* `1 I# a, W# ^$ N( o
首先定义点结构如下:
5 ]% p. T3 x$ i0 `3 f) w! S1 R: H
以下是引用片段:5 I" J5 ~' H" y' w5 G2 V9 S
/* Vertex structure */
5 t; W. _$ U2 c" _ typedef struct
6 H/ Y, }! I& u" q0 H* L { U/ A. ~% H* o. q$ m; i5 t+ w
double x, y;
- t# J. I2 _6 k6 K% n5 y3 B2 s } vertex_t;
# Q/ P2 F% S! ]! K( ]. b
. b+ n+ U; A& O) R. O& o
5 G& j$ W, h( ?& m! f( S8 K 本算法里所指的多边形,是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的,包括多条边在一条绝对直线上。因此,定义多边形结构如下:
/ s, R% e/ e" d
" ~& g# _' D8 {! W以下是引用片段:
+ F0 v6 w5 y6 X /* Vertex list structure – polygon */ 0 U% k& t( w: _% V* u3 M9 @
typedef struct
: {. f& W v% a9 i3 Y { 5 H( g0 J8 W( K1 ?1 @, c7 ?0 a- J2 B
int num_vertices; /* Number of vertices in list */
& B5 d6 v0 h/ P$ t0 g, Y0 @ vertex_t *vertex; /* Vertex array pointer */ " T, ^0 m& I2 ]
} vertexlist_t; " G& U! J7 e( N4 b) v& Z% }! I
- \* _. \# N( ~5 K! q3 j7 B5 u1 J: E q6 i7 ?5 J/ m. }+ }* F
为加快判别速度,首先计算多边形的外包矩形(rect_t),判断点是否落在外包矩形内,只有满足落在外包矩形内的条件的点,才进入下一步的计算。为此,引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent,代码如下:5 z: R! v& h# \, y8 n2 x
! Q w5 q, M! L6 H以下是引用片段:
9 L* Q1 }0 [6 U3 D1 ~ /* bounding rectangle type */ 5 a- D9 _; O0 v5 K3 c, U8 J5 `
typedef struct " Y3 G/ t& `% L- F7 s6 H
{ 3 }: F' `5 P6 F
double min_x, min_y, max_x, max_y;
4 z' \- f. w1 u6 G7 l } rect_t; . f, J `1 e* x8 o8 E1 _
/* gets extent of vertices */
- ]9 d* o) X) w2 Q& W void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
: ~: p) _0 g! x: f. P rect_t* rc /* out extent*/ )
. l; [3 \' S3 N {
7 E8 ^! b# a! Q5 t ? int i;
, B% V# ]& l! [. G3 y if (np > 0){ 0 F' g; a& L5 K: N! E
rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; ; Q5 r+ Z9 @5 L. T
}else{
1 J" E1 R; @% d rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
$ u7 b: v. i% y, k } 1 l* T! P- Z, D0 \3 s# N, O
for(i=1; i
1 t* ]8 A* Z) F9 o. O; h/ x$ S, @* @ { ) K l1 U& S- v5 J" e4 y5 E9 k0 L
if(vl.x < rc->min_x) rc->min_x = vl.x; 5 T z9 Q. J4 h' ?& a
if(vl.y < rc->min_y) rc->min_y = vl.y;
! v8 w1 ?: W# D6 N) B: t U if(vl.x > rc->max_x) rc->max_x = vl.x; 7 J) d+ G9 G# _ I( P' Y) Y6 b5 D
if(vl.y > rc->max_y) rc->max_y = vl.y;
) u* A6 w/ ?" I } 1 d* t7 A7 T+ l. ?1 Z( d
}
5 B8 B3 a: [ y3 x
7 A. o3 e! [$ G! S% x2 y9 d" k2 E6 M
当点满足落在多边形外包矩形内的条件,要进一步判断点(v)是否在多边形(vl:np)内。本程序采用射线法,由待测试点(v)水平引出一条射线B(v,w),计算B与vl边线的交点数目,记为c,根据奇内偶外原则(c为奇数说明v在vl内,否则v不在vl内)判断点是否在多边形内。8 A; b$ p* \; ~5 g" p. E& H
% q* B J, ]* P1 p9 e# { L
具体原理就不多说。为计算线段间是否存在交点,引入下面的函数:. M6 y s$ i: g" @+ L* F
; A1 |+ {$ e. T9 `% w" ^6 L
(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start,l_end)的同侧;0 m' B R6 g7 ]/ r" e1 E5 e4 m
* c! z! |* U# O3 ]) P/ d
(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
6 c% [6 F* V8 \* T( G0 D/ a& _6 ~2 N5 \6 E
以下是引用片段:4 F$ @" e2 ]% T, }9 Q
/* p, q is on the same of line l */
( Y+ E3 s: v1 A* r static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ - ^2 s# Q9 Z) ?' m0 H% _( N
const vertex_t* p, - E1 n1 F t! Z) ] k4 a. g; P
const vertex_t* q) 3 K7 B- s4 O8 [$ _2 U
{
2 ^2 a# `& A7 ?- f2 N double dx = l_end->x - l_start->x; " ?7 e( |0 ?3 E$ P
double dy = l_end->y - l_start->y;
' f5 I, h9 z$ i$ L( W8 { double dx1= p->x - l_start->x;
0 `4 J& ^$ m/ |6 g double dy1= p->y - l_start->y; : T5 O k& |8 z9 o, X2 a: l1 s4 t
double dx2= q->x - l_end->x; % p- V! F. p6 y# P' g3 ]
double dy2= q->y - l_end->y;
. Y4 n; Y# O- A# K return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
3 c8 @# `' i x6 t, t6 ` } : b" B% G! G. W3 `# v' L
/* 2 line segments (s1, s2) are intersect? */
6 ]- E. g8 M3 l% C5 E( _: o static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
; L/ Y/ p0 T6 i' L7 g/ ? const vertex_t* s2_start, const vertex_t* s2_end) & ?7 Y! Q% P& b E* z# r2 ~
{
3 Y" N; s" s* q+ J6 Z9 Y7 g return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && 4 g7 L+ Q1 h* n7 d+ ~! @
is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
W$ w& \' U, g0 A0 ?0 F }
7 A3 y* w6 x* m
% Q6 L% z0 O5 w: \1 u+ n
/ D% }- _/ X( m2 I* O- W8 G6 Z5 s 下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl:np)内的程序:
; T+ Z4 K8 O% I, w9 V
: p, }1 x# e/ ?# o+ M* P' d以下是引用片段:: \: B7 y8 S; x6 ~, ]
int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
" l! j! K" J% S2 d" W5 v! f3 G' @ const vertex_t* v)
8 I+ y! w% W7 b0 I" S/ ?9 @ {
/ Z* h: Q+ o8 F$ [ int i, j, k1, k2, c;
7 g9 A$ @4 X5 j1 Y rect_t rc;
) ^+ P4 |) ?' e$ d5 o9 H vertex_t w;
' \/ @* T8 a1 l8 M G if (np < 3)
L# L6 _) ?/ }: w return 0;
/ o+ q" [" @' c3 y; K! M: P" R* P vertices_get_extent(vl, np, &rc);
- T8 |. K( ]' [+ d- n% j8 m4 ~ if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
- G/ _9 m9 d) j/ Z2 \ ^, m# a2 Y return 0; 4 H8 m+ m; ^2 w9 Z+ N
/* Set a horizontal beam l(*v, w) from v to the ultra right */
# ]8 k) ^+ P. G w.x = rc.max_x + DBL_EPSILON; ( j) Y* U$ l9 J" j
w.y = v->y; % V2 k9 W2 r8 Q. z0 }
c = 0; /* Intersection points counter */
+ w1 |+ w' }# }8 `: x for(i=0; i
+ p9 x, l; ^3 s: I. H2 G3 b+ g6 L4 B { 3 T9 Q5 d% f& U" s
j = (i+1) % np;
* |& P" [% Q$ T1 w; v9 k& p if(is_intersect(vl+i, vl+j, v, &w)) : g! E6 K- j8 b) A1 C" U
{ / [, W0 s3 o2 l+ {8 d; F
C++;
9 B9 C) Z t3 b B- R% L5 a4 D1 z } * |0 Y, C$ }* e
else if(vl.y==w.y)
. H5 `1 `, x8 ? { 2 f6 s3 S$ |0 `( l
k1 = (np+i-1)%np;
P3 ]4 j" L+ S6 H% u3 Q/ [4 u4 u while(k1!=i && vl[k1].y==w.y)
, [6 ]9 a6 W5 I) Z+ s4 } k1 = (np+k1-1)%np; : {1 g7 a- a0 M# I& b' q6 N
k2 = (i+1)%np;
) w% o' X: [+ b while(k2!=i && vl[k2].y==w.y)
N5 F5 U+ [, f# |& N( q k2 = (k2+1)%np;
( |: Q: w; H7 k. |! x+ p if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) `* a4 d0 u! ]
C++;
: I, Q* \9 I" X1 m+ z if(k2 <= i)
7 [+ J: k, e E$ @* }9 U break;
" D7 G2 d! C& ~: U& M; ~) D) k1 c i = k2;
7 D( e1 M+ w+ m+ u, b# V } 5 B' L! w5 E, J; H
} ; @. j/ z$ ^& O1 V
return c%2; 1 g$ [! @. N8 M! s" u
} - `4 [' ^2 k$ G6 t/ M7 |
, T' D4 U3 B. J
9 N4 }2 {4 z3 }$ t9 a. K/ b 本想配些插图说明问题,但是,CSDN的文章里放图片我还没用过。以后再试吧!实践证明,本程序算法的适应性极强。但是,对于点正好落在多边形边上的极端情形,有可能得出2种不同的结果。 |
|