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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
5 ]( q+ g; E$ b' y/ g, j3 k7 a* g- `, F! Q* N* t
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。$ m) |8 @$ C  \: b9 d
5 f7 A; F" b8 |! M% F
　　首先定义点结构如下：$ ?1 h4 t% L, j9 Y
- W6 X- ^: I, K
以下是引用片段：+ t( S6 R% C. e4 F3 n9 ^) e
　　/* Vertex structure */
5 P4 @+ z# d: B' d4 L　　typedef struct
9 {3 t9 }! G- J9 `/ l. d7 ~  x4 k2 h　　{
4 c1 a9 F6 N: n# q3 v7 G　　double x, y;
/ S0 S  ^2 O! Z, l　　} vertex_t; 3 Z4 q5 n- @3 w# P: S- l

# M" d5 x6 l$ F" x
$ z9 \% c) a" `9 i8 W　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
/ ?* M; L% C% h) d! p6 c) a
9 s) L( C3 }" D2 p7 t  s4 n' J以下是引用片段：' L# v, ]$ |8 [7 r: h: V4 z2 K. ~
　　/* Vertex list structure – polygon */ 3 r8 K4 ]" U# N& X6 S
　　typedef struct 7 }1 O5 r  k8 Q9 A: w* L0 S# h
　　{
8 ]1 ^9 E# m; K+ S, G　　int num_vertices; /* Number of vertices in list */ 6 U+ N$ k/ A8 P1 F9 O+ E! s1 }/ D' c
　　vertex_t *vertex; /* Vertex array pointer */
+ t2 m: y# U6 Q- z4 w3 ]; V7 z1 }　　} vertexlist_t;
7 B5 f5 Y) D: E7 L. w- L
& ^; t* j! L- L: O$ O# L: M5 E: m0 L$ k. v8 q- S5 `- p
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
- E9 f4 e$ ^" y4 ~' C% I; t) r- Z( `( p  ?$ R9 t2 t" ?
以下是引用片段：
. E6 m, j, ?4 u　　/* bounding rectangle type */
/ {) h1 G9 R- ^) d' W; E* a　　typedef struct - A9 k# K4 P6 \4 v. h, b8 o! l# E
　　{ # Q. L& K& l4 f! ]1 e. s- x- J0 o" ]
　　double min_x, min_y, max_x, max_y;
# s3 |1 W) `/ a　　} rect_t; , h  }5 F5 C0 c4 ~" n6 \# y# D
　　/* gets extent of vertices */
( P, ]. ]: I# ~0 D0 I  v　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
! b4 O% _. Y- {% K( h4 w　　rect_t* rc /* out extent*/ ) 3 L/ r. h+ j8 b0 x( Z+ l
　　{ / ]' X) g. Q( Z1 {) O9 ~# S
　　int i; 8 Y: o& z- U: j
　　if (np > 0){ - H" f; ], F) T% u6 }
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
% S1 y2 T# l- h. h: O$ e1 O　　}else{ 0 u* U% s% I$ }: V* G) B4 o; ~1 B
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
. }1 `& m, H' B3 P6 l0 ?  K( s9 q2 J　　}
5 _- p! V9 V" ]! ]) B( @8 b　　for(i=1; i  : Q& b) n' b8 B  k2 K6 x
　　{ + z! e6 y5 X+ l, ?
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
7 Z) g! `! S+ T  c! C- F　　if(vl.y < rc->min_y) rc->min_y = vl.y; 4 V( J+ z1 C* F
　　if(vl.x > rc->max_x) rc->max_x = vl.x; 5 U1 v7 X: w* w4 }  a0 Y
　　if(vl.y > rc->max_y) rc->max_y = vl.y;
0 ]( K- C: ^" B) G0 X4 A* q. A　　} 8 ]$ y- }+ c; b% U2 a7 P7 h5 M
　　}
! U6 K9 f& i9 a% v, s- D" t$ ~/ i: Y% L! P' C% h0 o* v

3 d! b9 k. L8 D　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。& t, x% t8 }* Q$ s7 Y9 f

& G2 u7 o! A" S9 G5 F　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：$ j: _2 r! `, V0 w
0 ^6 f9 B1 W6 P. {* q- g" o
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;" k9 m4 R; C' i# a2 y' n

# C2 ?) l/ ^* V. D　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
2 n. X, W& q0 c- C- f+ }, V1 N) o6 z  I/ I1 z; v- l
以下是引用片段：- @; U4 m7 `% A5 G! [' s
　　/* p, q is on the same of line l */
  b' Q3 M6 I0 Q  \　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
4 C# _3 H% P" _) J9 l9 E　　const vertex_t* p, & F2 H% p! ^3 q% ~- a
　　const vertex_t* q) ' _% e& J2 ?/ L) M" V
　　{
5 P5 r0 A! ?0 k　　double dx = l_end->x - l_start->x;
; r0 \4 _1 j; a. X- q　　double dy = l_end->y - l_start->y;
5 P$ x4 m4 R7 J0 ]8 ?2 ~　　double dx1= p->x - l_start->x; ! d4 V7 {, h& K. G& ]3 a+ m
　　double dy1= p->y - l_start->y;
" q* D- q4 O1 @" e7 L　　double dx2= q->x - l_end->x; 4 Q. b, p! L- L+ ^' U1 R
　　double dy2= q->y - l_end->y;
  a" U8 J, v( d8 f3 \) P. @; ~7 [5 u　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
9 Z5 P1 Q6 C, q/ E* {　　}
$ Z  ~% U' V# q4 E5 I　　/* 2 line segments (s1, s2) are intersect? */
+ A6 Q" g8 f* [2 p　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
) \6 J7 n( E+ g0 N1 L- x- v　　const vertex_t* s2_start, const vertex_t* s2_end)
( `  y# ~  Z! H# a+ {/ U6 w　　{
; t! \" U: D, W5 y$ _, h　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && / N6 f+ ^& g2 {
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; $ }) v0 [# Q4 K" ~
　　}
; M0 M  {4 ~3 \; d* Q; M9 f, Y3 D' L. P' R" s+ a/ N1 S
6 O% t4 k" T/ x
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
9 x7 ^& P# E6 R0 L0 k0 K
5 P4 z- U; R( O( Q8 ?8 O3 S) E- m以下是引用片段：$ _7 r5 ?5 N% N5 f- B
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ & o' }) h8 ^, D; W0 p
　　const vertex_t* v)
/ _+ }  L6 o' v' B9 h. _+ e　　{
$ R/ p$ {* K% R% ?# f8 Y　　int i, j, k1, k2, c;
* U/ ]' Z& `, B" M& F　　rect_t rc; 4 _3 w7 Z7 m" N! i( e" I
　　vertex_t w;
2 N7 j' Y! }, O* [7 Y　　if (np < 3) 8 r' ?3 R4 h, s3 h. N0 `3 f
　　return 0;
+ ~- b, g* p2 q& e　　vertices_get_extent(vl, np, &rc); : {, Z8 l. C( c
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) / _& F7 l, m3 h, v) A! V& E
　　return 0;
! f' \  p! b) _- }3 v) M9 e　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
; ^8 _8 E" m5 X# a: ^! _　　w.x = rc.max_x + DBL_EPSILON; ) ^; m6 j# s3 W, m3 @
　　w.y = v->y; ! g% y" M7 b) t8 F$ C: j& r: I
　　c = 0; /* Intersection points counter */ 1 h  V3 q% i  @. o" o. Y7 L
　　for(i=0; i  
' [8 Y; ?$ {0 T; j5 s6 H: V　　{
$ m; I( a/ l9 A6 `　　j = (i+1) % np;
& T/ _" R' b" T0 P' x( o2 w　　if(is_intersect(vl+i, vl+j, v, &w)) " Q- p! W: |) j
　　{ 5 e6 V7 P  r  o9 }6 v6 {0 r1 ]
　　C++;
+ V) A) T# j$ l7 v' f" p, a- }% t　　} 0 x: H7 c+ _8 a. f
　　else if(vl.y==w.y)
. l& D! Q. N! q0 P4 ^　　{
. |7 x/ L% C- t- S; j　　k1 = (np+i-1)%np;
4 o/ ]4 y" s4 Z0 U# D( k8 F　　while(k1!=i && vl[k1].y==w.y)
# ?0 P' `7 E5 L# K2 _　　k1 = (np+k1-1)%np; * ]1 f0 {1 Q/ H( T
　　k2 = (i+1)%np;
3 M0 Z7 r" g6 t# D" y　　while(k2!=i && vl[k2].y==w.y)
6 N% u) W. p( G8 H* ?3 x6 u4 A6 G　　k2 = (k2+1)%np; ( I9 i. q7 S1 X$ O6 ^9 Y$ u; ]
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
# _+ Q& W) G, D8 |' {/ o& H9 q　　C++;
# e# N+ Y3 c. u2 o+ J; Z3 s% _　　if(k2 <= i)
9 g. b9 U" a: S' j! M　　break;
5 w( ?5 V8 x, ?2 U0 s$ w0 H# {- e+ x　　i = k2; 3 X3 x3 k* L- R3 S4 O1 Y4 s
　　}
, f8 l0 z- C6 o, z　　} 9 Y" `6 Q  J4 S4 J; {
　　return c%2; - t! [3 P( H" N& i/ l, i
　　}
* w5 _5 ^5 n& L3 m7 b. v! s- A
8 ?$ [+ d, X! ^! r4 C& }: t0 t6 z, y" A4 G2 t$ n0 T
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。