|
  
- UID
- 133
- 帖子
- 51
- 精华
- 1
- 积分
- 186
- 金币
- 55
- 威望
- 2
- 贡献
- 0
|
C语言中显示 点在多边形内 算法
本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
: e" b( F. q# J( [' V- G: f' O# g4 Y r2 C8 ^
这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
2 [+ K" M& [! {0 Z! C$ o0 X) E8 f1 J% s h; `
首先定义点结构如下：
0 m# W6 P! G# a
& R4 p# F+ d; y6 S* L以下是引用片段：6 Y6 _$ x8 r% c
/* Vertex structure */
. V r8 H2 B$ d( r) U! f typedef struct & J: m0 f8 l+ S* r3 U
{ / q; ]/ C3 X$ n- O% ^
double x, y;
% u4 g0 W, ]0 B# `2 y# R } vertex_t;
3 \4 _# p0 |5 ^& M/ ?* e$ O! S+ ?$ [3 X( J
* z0 m4 }, x$ w6 A" x9 X* d+ R8 o 本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
0 @$ H! _- g$ |% x
+ H! E0 W" v1 s! Y$ D+ S以下是引用片段：
5 ]: D9 Z3 ^9 ]! k" r% I' ^. m /* Vertex list structure – polygon */
$ T: p% h$ S+ l5 \( j, \: G typedef struct 6 a- ]0 L) i2 ?& f$ v! k
{ , U& S( q: _: k6 U: x0 q: s
int num_vertices; /* Number of vertices in list */
" ^- G: l M, G/ I, d vertex_t *vertex; /* Vertex array pointer */
1 b8 _, x1 O4 ?! s } vertexlist_t;
% e9 [1 l4 |7 @# Q" |% N$ r% y. P6 H' J0 H
0 \1 ~( `* }: ~5 |% m 为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
$ q3 i V8 Z9 @# c6 Y6 o
% e/ {: Y5 U1 K以下是引用片段：
! O9 @( o( B9 [1 w0 W8 E/ ]+ t /* bounding rectangle type */
: Y/ q- \0 G# ^% J& L typedef struct
$ s3 S$ n4 ^1 ]9 G) p- |6 |0 {+ B {
a0 Z! a5 ^5 `& m' `. i- o8 p double min_x, min_y, max_x, max_y; 6 M& _, s8 S) B, ^
} rect_t; 2 E6 }$ Q! Y- ]. U9 [ z
/* gets extent of vertices */ 1 v+ a2 |- }. m' `' K
void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ # B# E, f \( L# g$ @& e& K
rect_t* rc /* out extent*/ ) & o& S5 t, R1 L9 b P0 t5 z
{
- p: ]7 f5 l8 Z2 F4 B- J8 x int i; 2 h8 ^3 d- X9 y8 A+ W+ V: Y c
if (np > 0){
- G; Y1 [+ Z2 F$ w( L- o9 x8 a rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; ) {" e4 x4 O9 X; v; G
}else{
$ [9 R! V2 N2 s! h. l* V! ^/ [! x rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ ; l, g" O |2 ?
}
; r0 W# R; X' o) ^ for(i=1; i
2 T" u# E+ b( Q8 X: ~) f+ o { : s, h) i5 w( y ~5 J; \
if(vl.x < rc->min_x) rc->min_x = vl.x; ' `* c5 M, X( g4 h
if(vl.y < rc->min_y) rc->min_y = vl.y;
4 } Z/ e5 F/ U* y9 B if(vl.x > rc->max_x) rc->max_x = vl.x; $ m3 s4 u( b) u/ j: i7 j
if(vl.y > rc->max_y) rc->max_y = vl.y;
) A- f; b7 X& ?/ F( G+ e } / ~$ _: ]. e0 A& P
} " K# N: T, @1 v4 E. R# O
( K$ `& Y+ A0 e) H
K5 I! H0 G( C5 {/ ?, B 当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
: {( s, x N4 B% r' I+ X b; w4 S5 A, L5 ~) @
具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
0 ]2 C# a% ]$ O# S4 L) l4 @% ~5 o5 D. p9 c0 e* L
(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
5 h ]3 M+ T4 O: d
% x& b! U. k4 e (2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;2 E" b: U& w. Y& {
% S; A+ c* i/ F以下是引用片段：
0 ^7 ]/ O+ U6 g- M* i /* p, q is on the same of line l */ 7 x& n: y, Y$ U5 i3 L$ l; M
static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ 8 U; g* ?4 X1 r$ p0 Y$ e6 A8 D
const vertex_t* p, 6 Z. U. o, ~4 k3 T9 c- M
const vertex_t* q)
9 j& Y7 u, F, e. V C {
( Z) ]* X/ F E) ^$ R% r P( q0 `% Q double dx = l_end->x - l_start->x; * W' B A( \3 H' z
double dy = l_end->y - l_start->y; % S8 x" }+ }$ ?- l" H/ Y' f
double dx1= p->x - l_start->x; + h" M' g' O, x4 B* d; n0 C1 ?: H
double dy1= p->y - l_start->y; 5 ] P9 w+ v* e, R. Y
double dx2= q->x - l_end->x;
$ B7 k9 U7 r" u5 S8 f- ^' d7 k double dy2= q->y - l_end->y; # ^( c$ {0 r/ z8 L9 q
return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); 0 o3 u! X0 H/ }% o
}
9 D" K" G' X( @6 s /* 2 line segments (s1, s2) are intersect? */
$ i2 d5 K, Y! y" S1 L% B8 O/ n static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, 6 v$ K7 ~* s5 u6 R" K- [( X# ^1 C
const vertex_t* s2_start, const vertex_t* s2_end)
* M# c( Z$ C( ~; C {
, S/ n% z3 U1 Q' U' ^ return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && 3 E7 g& L: y; v4 U
is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; , ^2 L& K( U) } p/ U& |: B
}
k/ J6 \# c) G& P" E- D
' u% W4 |$ W* S s. V; Z6 p
, {9 L0 B, S# w1 Z m h& W6 \ 下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
4 u5 J3 {. J/ ~- }* H. O+ t8 b* I8 s. h7 ^
以下是引用片段：
& Z" p. |7 [& R% z int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ / F" K8 r, p6 V2 \4 {5 L
const vertex_t* v)
2 D) C6 {% @8 s5 Z' z { ( @6 ^" e- m* B* N$ L' i) D
int i, j, k1, k2, c; ) H# f! Y y& N' X _$ a
rect_t rc;
5 N5 f/ H0 c: c3 s* u4 D vertex_t w; 0 a/ H9 `) m% z
if (np < 3)
" B. n. p7 Z/ w7 T9 w return 0; ; f% e4 y1 R5 B4 _) T2 w3 v$ {- z% e
vertices_get_extent(vl, np, &rc);
) W( ~0 w. t8 `- U) d' W/ Q" t+ y: n if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
8 y# t9 ?) b- S. g! @! G1 | return 0;
: R4 i+ d2 D7 t2 ? /* Set a horizontal beam l(*v, w) from v to the ultra right */ _" a; U) D3 T4 o8 [
w.x = rc.max_x + DBL_EPSILON;
6 F6 ?: G1 y1 G$ g& n: H9 D w.y = v->y; 8 E* K0 t5 r1 e& [9 K5 G
c = 0; /* Intersection points counter */ 3 b, e" M5 \2 G1 O0 A
for(i=0; i 1 R( C2 w+ b0 |/ u2 d
{
1 z- W7 F! l0 ]+ S1 B9 h j = (i+1) % np;
, I" x9 U' Y# y I0 k* G if(is_intersect(vl+i, vl+j, v, &w)) * c9 s W9 a9 }2 }1 O6 X
{
; h. o1 ]3 X9 f. i% V, v C++; 6 b0 |+ Z. r& f+ {+ u
} " Z' j7 y9 U% {/ `" v6 M* c
else if(vl.y==w.y)
; y0 {8 w, p5 R! E {
- a% c: A. ?( }, W' g0 W k1 = (np+i-1)%np; 1 c: B* y' @; Z% ]/ O
while(k1!=i && vl[k1].y==w.y) 1 b+ @. ^4 L9 ~3 V- c4 v) Z
k1 = (np+k1-1)%np;
& Y( s9 L- \8 @5 R k2 = (i+1)%np;
& @% h- H3 w5 R$ X while(k2!=i && vl[k2].y==w.y) 6 h& S7 T$ `, L- v: v/ }
k2 = (k2+1)%np; ' d) C' k. V5 ]& Z$ l
if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) - z( G" c3 Y% J9 [ [& X' {
C++;
( _4 Z( ]; d4 Q# h" J; D5 ] if(k2 <= i)
3 t' {: Q) q/ Q break; " R/ R" w# Q0 d' s# b6 F
i = k2; % w% H% }. x1 D
} ' T) D, p2 _( w# i. F y5 H- R
}
, g: G3 x) ?; J( l7 A+ y9 O return c%2;
5 q, p8 w2 j1 Z6 a- G }
3 D+ e) K/ ]9 g+ _
( ^: e. [1 ?' {8 e/ Y! ^
3 x+ G4 b* |$ ^; _& w& n: q: ]: w, v 本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。 |
|
发表于 2008-1-21 17:20
|