C语言中显示 点在多边形内 算法

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。: R( {' V9 M+ y0 ^# g5 k8 o
9 o6 h2 ~) D9 }
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。% m, x4 l+ ?$ S1 q7 D. e+ m
* g# \2 P4 r: i& N) }- A
　　首先定义点结构如下：% h( X" t4 {: x  d" B8 d" M' {
+ g! |! k. p7 {4 V8 k% r
以下是引用片段：" j# a# y2 a- v
　　/* Vertex structure */
6 @+ X  l% E- h: n7 {  {+ e　　typedef struct
' Z6 j, m7 E. k6 C　　{
4 \5 c/ ]) \  f$ n# C　　double x, y;
2 t4 q9 Y9 l6 ?- X/ |: r0 _　　} vertex_t; 4 t+ u' i- j# d" P& v

& z: Z8 d. i/ n/ \7 l5 j; C# O+ h$ ?+ f6 n" |+ m
　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
" E9 f+ s3 [+ Q' t
% L! n, a* V: [" [- V6 N以下是引用片段：
( X4 J; Q$ U, E. l　　/* Vertex list structure – polygon */ / c+ E4 \& p9 y9 _8 u$ P
　　typedef struct
- H9 u+ \% z8 ?. G# i# y! a　　{ & H; l9 N% x, F7 g7 R+ t
　　int num_vertices; /* Number of vertices in list */ ( @9 k* u3 @6 A  B& C
　　vertex_t *vertex; /* Vertex array pointer */ 7 M& T" `' d& Y0 B; O7 {1 ]+ u
　　} vertexlist_t; * o) a6 Q' C% R% o0 Q1 _$ J

4 h! Q! C7 n  W1 X  f/ E: ~  k- w! ~. X# Y8 ~3 _) w
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：3 a4 n( s  W: M( J+ K
: Q0 k; a3 E; I* @  W: e" E
以下是引用片段：$ C5 o) p! W; B1 r4 i3 T0 @
　　/* bounding rectangle type */
, I  e$ F; i6 T) m7 z4 s6 o' O$ w　　typedef struct * H' d7 s: ~  c
　　{ * C4 D' N" r  X! n3 T+ `
　　double min_x, min_y, max_x, max_y; % e4 A1 Y$ J- J6 s9 h
　　} rect_t;
) `1 d, U' b0 x' I　　/* gets extent of vertices */ . P1 u2 y0 V9 u( K5 `6 Q& Z
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ 9 w+ ], \0 T; k& j
　　rect_t* rc /* out extent*/ ) ) P6 f* H8 m8 N
　　{ # V1 I1 z$ j; |# v9 n  x  x, t
　　int i;
' J) }$ d- J8 h  N) r" ^: `　　if (np > 0){
$ ?% m2 _" F* {" S　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; ; l$ U2 F& B. ]2 j  {% L$ d- N
　　}else{
. k: o9 p+ G/ p  a$ X5 @4 s5 x　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
. w! C) D6 h2 q( o1 ~/ H　　} % ]$ d2 m( K* V, q) I6 }
　　for(i=1; i  
/ c5 a" \. B; J" X9 `　　{ " [; s# ?) V4 a. ~( p+ _: I" x  H8 V
　　if(vl.x < rc->min_x) rc->min_x = vl.x; % ^3 j0 d% S' g) O6 q: l
　　if(vl.y < rc->min_y) rc->min_y = vl.y;
2 u2 w% P; Y+ e* k7 ], g- |　　if(vl.x > rc->max_x) rc->max_x = vl.x;
2 F) A+ J4 W" ]% q9 q% k　　if(vl.y > rc->max_y) rc->max_y = vl.y; ' q: q) m5 }2 S
　　} # j# t/ H# o! J5 E
　　} ( `: Q" W) n+ G: [. H! N% ], s

8 o4 p6 R0 ^4 n9 [0 T/ `' Z3 v; @/ |1 G/ `+ K0 \* C9 O  H
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
7 j3 ?& j# s9 W! T, c# e( _
: t3 \$ F  w+ N　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：+ p5 W  @, O. t+ D3 w: |  D
: t; m* |: T7 r- [% T6 f% ^
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;* D; U+ ~! s$ ~+ |1 T1 E: D

( |$ r. W" \( \7 s/ v7 u/ W　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;  R4 t; s* @( p9 t
) [7 F3 I; P$ b, B/ M
以下是引用片段：5 s, C( L3 A# i0 o4 Q" n) o
　　/* p, q is on the same of line l */
# I4 t2 f& w: ?9 R: N1 `　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ 1 ]. Z; j8 B  m6 k+ \
　　const vertex_t* p,
+ z. J8 y$ [9 d2 ^3 i　　const vertex_t* q) 0 u8 X' F0 O% \0 O% ]
　　{ ! K2 M; ~" p, T2 l. j) b3 \
　　double dx = l_end->x - l_start->x; & m& n, G/ E# E& A1 j
　　double dy = l_end->y - l_start->y;
* H$ R, I- a3 @) j; `　　double dx1= p->x - l_start->x;
0 Y6 H/ @) D% @& F5 e( s5 K　　double dy1= p->y - l_start->y; " y# \& J1 N. X7 P  Q6 G
　　double dx2= q->x - l_end->x; 6 ]6 z- ~9 ^( W; @" b: ^
　　double dy2= q->y - l_end->y;
- v6 e0 i" Y* E6 T4 s' @' [; `　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
+ \. v. n; S+ [' s+ t3 O　　}
8 ?" J7 K! H; }) M2 y# y+ J) p　　/* 2 line segments (s1, s2) are intersect? */ 6 }6 j7 N& l6 Y1 j; v
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, ( P* `# C8 M) V+ w3 X
　　const vertex_t* s2_start, const vertex_t* s2_end)
2 e( ?! Y6 R0 j1 r: S　　{ 7 ]/ m% }2 q* A) p
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && $ Y: Z" U7 v, ~6 d* i# r" i# W3 T5 `
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
8 e8 ]' ]4 m" |, n* _　　} " F, w  g9 `4 D5 A. {# C. f& Q

0 [8 ]: Y2 l" O! e2 d) v, d
! Q7 b+ ~$ B9 R0 z　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
7 v7 p* c8 U- M. L5 ~+ A) G: ?
& E! W3 `5 H# {以下是引用片段：
2 s  h8 ^6 ?2 V! @3 ?; W　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
1 E  a" Q# w5 w, y　　const vertex_t* v) 1 o0 f3 A+ u- E1 T/ t8 Q
　　{
3 i: u: E6 H/ m　　int i, j, k1, k2, c;
" z3 Z4 t$ Z6 h# h9 j2 K# d4 L　　rect_t rc; 2 T7 R; j( h8 A* J" l1 i
　　vertex_t w;   F- N& J- ]8 |6 Y2 }$ M* c! j$ c1 _) l
　　if (np < 3) 4 p& S7 f- r: u- G" ~
　　return 0;
# c5 z- j: {: P+ y8 Q; f　　vertices_get_extent(vl, np, &rc);
$ e7 n' c/ q- Y$ g$ o7 S# @: z　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) + j' Q, a& ~) i& Q6 T9 J
　　return 0;
! d9 K! a0 G" p" e: z4 i( P7 \　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
6 {+ v3 ?# l* ?; J; G　　w.x = rc.max_x + DBL_EPSILON;
( i. z# s+ m. e; y- t　　w.y = v->y; # ~, E5 z# ^# m
　　c = 0; /* Intersection points counter */ - X! a5 V# J! x+ R
　　for(i=0; i  0 J% i$ t. t: e; @
　　{
* {& U& a% Y/ h0 G4 Z+ w  A6 c* s' y　　j = (i+1) % np;
9 z' u3 a. U5 W8 R* s+ R: e0 }* Q3 L　　if(is_intersect(vl+i, vl+j, v, &w))
" [6 c) Z& G6 A# v( d　　{
1 [/ Q3 n5 D! `　　C++; 1 A% s" y2 M  B
　　}
6 B- S8 s) J% S2 T　　else if(vl.y==w.y)
. f4 D" h) s/ Z1 ~1 C　　{ 1 Z! g% `6 x2 O& E  l8 p" L0 J
　　k1 = (np+i-1)%np; - ?% Z8 T9 g# L) X  }" c
　　while(k1!=i && vl[k1].y==w.y) + V2 B  Z! E) ]  o6 j
　　k1 = (np+k1-1)%np; 9 d8 S% ]" ?8 f& {( ^' J" \
　　k2 = (i+1)%np; % z, _# m, N9 y
　　while(k2!=i && vl[k2].y==w.y)
0 g$ Q7 p' Z' o　　k2 = (k2+1)%np; ' A5 v+ |, ?1 Z* Y" ^  Q  j
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) 7 O8 J4 g0 c" Y, c: H
　　C++;
" B' r" |( z: S1 m1 c　　if(k2 <= i)
0 c4 i) ?6 @! }　　break;
  D: j+ ^  D6 f/ H3 F0 g　　i = k2;
, x. r- x" ~6 x5 V　　} 8 w- K: c. P& P- j! c% H/ P' H
　　}
1 x9 y$ f" A+ |7 R( r4 x5 t　　return c%2;   z. z* M! m  }( M
　　}
5 K3 n. |" V% t1 j! O; B( C+ K/ f. @' \, {4 R3 o" V0 C

' b$ ]! E9 V" R, j　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。