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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
' ~9 t5 r6 l" {, f
% |* a4 u8 S# F' ^, I　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
4 M# [$ u: Y  M; h5 W1 ~: g
& i8 q: w  F" u8 p" H2 |" o　　首先定义点结构如下：$ D6 V1 {, B$ ^( j0 d

- @. w2 n9 x' w. _$ f: ^以下是引用片段：
& f# [6 \: H/ X% n1 ?# S　　/* Vertex structure */
/ y6 n2 R+ q" C# Z3 p7 D% p　　typedef struct $ @9 B& {: m% O6 X
　　{
7 @8 M3 A9 ~: ~! _' F0 H+ L. w　　double x, y;
) ?3 x5 j+ r: Q4 x　　} vertex_t;
( q0 m# J- x, D  {* j2 }" V; a( O) b; n& Q) v( u

  H  a/ [2 q' n2 C- V5 D) o　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：& ^4 H6 f1 }* E) C# Y8 A: o

9 o4 A6 n3 R/ R2 G# M  l& s( _以下是引用片段：4 l$ L$ ], b2 P! |
　　/* Vertex list structure – polygon */
0 _2 L0 y1 z7 R2 A  e- f　　typedef struct
& ]6 H& s( v4 ~　　{ 6 z% F- l  y) ~1 r' |. p
　　int num_vertices; /* Number of vertices in list */
6 G3 D( I% \. [7 N  V& ^* b) _( Z+ o　　vertex_t *vertex; /* Vertex array pointer */ $ R) _& _: T' d! Y- ^1 u
　　} vertexlist_t; * C7 c, ]3 y9 |+ L- l

; k& ?  ~6 p3 \2 z. b# a0 q1 P2 Z7 e# j; @3 T1 i, _/ m4 p
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：; ?- _( l- ]0 \+ n4 D$ J) C

+ G7 a. W/ d0 g* h4 ]( w( {6 a# s以下是引用片段：
7 y' c/ c; g4 C4 Y; @. }4 j　　/* bounding rectangle type */ : [( {" ]% q4 v9 @" Y( s5 D
　　typedef struct
4 u  W' b, F2 r& f2 a. f+ I　　{ 5 N$ @9 t; M  g4 u) E4 p  L
　　double min_x, min_y, max_x, max_y;
% U8 h1 u  r3 v4 p8 l　　} rect_t; & q8 _" p$ H3 D0 b4 J
　　/* gets extent of vertices */
3 s2 w/ M' o( E! C1 F# S! I0 ]　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
9 z1 }4 N' X4 T4 G, k　　rect_t* rc /* out extent*/ )
3 d! F2 q# ^" R6 \1 H+ S( U　　{ * q! a& l0 m" r$ V* [9 }7 m/ G
　　int i; 8 p: ^: C: {& u
　　if (np > 0){ 9 x( ?- E+ ~5 L) i  I& [* j
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; ; W5 h7 [% F) @% [0 W$ y2 O
　　}else{ ' @$ j( i( M6 E# j# b" ~7 d
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
2 @3 _' N8 \9 i- {; V4 `' e) e/ X　　}
$ L8 N2 I9 o3 c" ], r8 {. B$ T　　for(i=1; i  
5 V1 Y/ ^, x9 F8 G7 ^  b6 [0 |4 c8 B　　{
+ O; }5 v0 L: ~4 p& ?　　if(vl.x < rc->min_x) rc->min_x = vl.x;
9 E/ I$ Y$ b. k& d　　if(vl.y < rc->min_y) rc->min_y = vl.y;
. c! n7 D. E0 k4 k% u7 A! }1 m3 G) y　　if(vl.x > rc->max_x) rc->max_x = vl.x;
; ]3 [, y, g7 G: _+ |　　if(vl.y > rc->max_y) rc->max_y = vl.y;
& E! W8 W- }3 w; w4 s- j& @　　} 1 G2 U$ M  h5 _7 \  u
　　}
8 j2 v6 m. a) i# Y$ J: E. ~  o0 z% ~& X
3 ]7 b3 @* N2 U9 [% W% K- ~! \, {& b; K" }
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
: f& d; o, w4 D- j8 [! {, H/ p
6 J0 H. k4 c' c( E  Z# X5 j　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：1 Y: j4 b2 w6 \4 ^

7 X$ c3 R/ s' @5 [( l　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;. Z( ^4 Y: a% j* I

; T$ _- I. z7 l! O- u　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
0 }7 X* c, H% H5 U5 y: m3 l
& b5 ~2 y# |4 X以下是引用片段：
. W: f& f& u$ W5 m/ K　　/* p, q is on the same of line l */ 3 I, v8 P7 h# B& C& G! d
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ & m* z& ]3 m1 @# k# \& j
　　const vertex_t* p, & c# u$ ^' U! Z1 W, N6 t
　　const vertex_t* q)
. H3 i  Q2 v2 F2 |( d　　{
% d/ x! h  J6 p5 Y# ^, I　　double dx = l_end->x - l_start->x; 4 p# T8 w) F4 T
　　double dy = l_end->y - l_start->y; 1 O4 i) X2 i! ?. D4 _& o: C
　　double dx1= p->x - l_start->x;
& h% T& \2 u) r8 S3 Q　　double dy1= p->y - l_start->y;
% K, _. [, J! ~　　double dx2= q->x - l_end->x;
/ l! T% y+ x8 y* l　　double dy2= q->y - l_end->y; 9 a3 I. ~. K/ T- ]
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); 9 d% v, v& q  s" l& Q& f8 L
　　}
7 C! x5 W- B6 p1 s, u, C0 W5 g1 u　　/* 2 line segments (s1, s2) are intersect? */
' ~  W. z: {1 y* ^3 H　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, - W# z3 p. y- }9 q- ~
　　const vertex_t* s2_start, const vertex_t* s2_end)
$ J4 P5 F' C2 [; P　　{
) e8 E0 |, N) S' ?* L! X　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && + M8 ~! w( D2 }/ p3 ?7 N
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; 6 [. m9 v+ {- S0 H
　　}
, b, ?, D/ U2 H: c- x3 {
- l/ \- t& C) s
/ R; q- z* |- U3 C- t3 S8 q　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
) z* m& L% c2 y" {1 a
/ ]! f, T1 F$ d+ F: I" n& W以下是引用片段：, X, M0 Z& i+ p, I& F+ e
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
* ~5 O( i! F+ |　　const vertex_t* v) ' }$ A. u2 d: B6 X1 j
　　{ ' h5 I! p9 v! b8 U  h9 Y3 t
　　int i, j, k1, k2, c; # y9 u- X9 A/ t. [9 ]* ^
　　rect_t rc;
6 p2 Y( C5 ^4 ^9 k* @4 F　　vertex_t w; 8 ~6 V" o; P/ G/ U6 }
　　if (np < 3)
- y& D) o9 G# ]3 P/ F4 g　　return 0; 1 o3 V+ k$ I' ]; |: U
　　vertices_get_extent(vl, np, &rc);
0 C$ j( K% c. k* ~$ T　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
3 D4 e. Y. ?; M- n& s: M! x4 Z　　return 0;
. K3 x3 L" R/ [7 U% ~　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ ; @1 W! Q8 J1 H3 `. C' W
　　w.x = rc.max_x + DBL_EPSILON; ! o& G0 Q4 w! ^
　　w.y = v->y;
5 U( N) v: u  C3 s) x9 ]　　c = 0; /* Intersection points counter */
* g+ \7 K4 B3 u# f9 E0 |8 p( F  }　　for(i=0; i  ) \, q5 k5 E  V5 @( ]
　　{ & C2 ]+ N( l: k
　　j = (i+1) % np;
- ^6 R  G  M0 K$ c　　if(is_intersect(vl+i, vl+j, v, &w)) 2 ~. v# i( s) W1 {* b( U
　　{ 7 v( P# b9 Z5 l3 _$ S& p# U
　　C++; + A) ]1 S7 b7 u+ t
　　} 3 k4 s4 Y7 e7 _8 g  S+ [
　　else if(vl.y==w.y)
5 ]; E5 b7 F9 O9 Y) ^, E; W; W. ?　　{   `& n! |) ]; A! @5 n
　　k1 = (np+i-1)%np;
; c; y4 v4 M. I  J6 g) r　　while(k1!=i && vl[k1].y==w.y)
/ s# c% i  J# _% i  j6 F# q　　k1 = (np+k1-1)%np; ( U1 J/ k/ W4 F4 R$ e
　　k2 = (i+1)%np;
4 k$ ?2 e1 ]7 p* \/ ^) h  Y/ ^　　while(k2!=i && vl[k2].y==w.y)
2 j- M7 \; V& h- j; K　　k2 = (k2+1)%np;
, @1 {  a/ S& L$ W; p2 n! ?# D　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) . o2 D1 n* b8 H, m
　　C++; 0 E7 q) g# h9 B
　　if(k2 <= i)
% g" }/ {+ h- Q% v　　break; 3 _7 V" G9 e  Z& }. f! Y
　　i = k2; " k: @( x% C+ b% c3 X/ V$ h+ x
　　}
" |, l. s4 H1 m6 B; U3 R) ?　　}
2 y7 F* C6 @4 ^　　return c%2; 0 z7 I% a6 A' ^$ b' l; m5 r
　　}
5 o1 ?% q" h5 \, I  a
5 h- @$ P0 B& K- _, I6 W0 S$ k/ W) u1 |. a6 A
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。