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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
' N4 `/ ?+ h& _4 I! _5 n! E! r1 j+ J& q8 U0 |! c
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。4 g3 w  F% ]2 C0 q6 Z  u% }) D
- S1 e* r8 T1 }- o( k/ u
　　首先定义点结构如下：
6 H  ^$ o* a3 l7 X" x
3 t; n: q) U6 B, q以下是引用片段：/ G. U: ?7 f# i' D. H
　　/* Vertex structure */ 3 I, J1 n5 H; s2 ^
　　typedef struct
# ^: ]2 m( M) F　　{
" d4 j( ]4 x( s, Z$ a% {% c0 b　　double x, y;
) b$ I8 ?8 w! x* `/ V/ A　　} vertex_t; - E3 h& }4 E8 }) T  r# M( \
: ?) y7 K& l8 S6 M& U6 P

4 {- ]# u3 t* Z! e) I  [　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
" w5 Y7 G" q" }+ N( J& _- ^: I  t5 f* w" l0 L3 g
以下是引用片段：
4 @& j8 ~/ b# B* O　　/* Vertex list structure – polygon */
+ S. f$ r& [/ w! h5 e　　typedef struct 8 k- N8 M" |3 A1 ^7 a. V: v
　　{ 8 F3 [2 G. s% c* z5 A
　　int num_vertices; /* Number of vertices in list */
  T2 j' x! h' ~! P3 Q* s! A6 Q　　vertex_t *vertex; /* Vertex array pointer */
1 [6 q: d0 [5 I3 A% b) y3 {　　} vertexlist_t; , Z* ]) a8 p7 x) _* E2 M5 j/ j
8 C: F' }1 J) c, w4 I

- R- r& j: o( A, h! |+ E　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
% f/ G/ q5 f2 v* P' {
: M& C3 @) c6 j. @- B) E以下是引用片段：
  M' l9 M1 X0 y& k: p6 h5 y　　/* bounding rectangle type */
8 J7 n5 B! Y/ V3 I　　typedef struct
9 b" V  z$ ?1 O/ E! x　　{ ! C' |6 z) M9 c3 X
　　double min_x, min_y, max_x, max_y; + S1 t7 E) }) J  H* c
　　} rect_t; ) V2 {) _8 v" o1 l3 s% S
　　/* gets extent of vertices */ 4 l+ _3 r7 `) l) l
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
5 p: I7 V9 o( `  o/ p　　rect_t* rc /* out extent*/ ) # C, O5 P+ }4 L# O9 V, q6 v
　　{ 0 H. w9 l6 A, w, o+ ^1 Q; P6 X
　　int i;
' z% i  z2 R" y+ H3 H　　if (np > 0){ / t9 O2 w5 p5 Z+ a
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; ) m' ^7 @" T6 E  ~( @* W8 Z7 ]9 e
　　}else{
* Y) V3 d  x" E/ \: S+ m/ X, o　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ 0 H  K5 t& I$ Q+ k( k
　　}
8 k! J0 q: \+ y6 e　　for(i=1; i  * ?8 C8 o8 R8 J
　　{ 7 A6 i7 _0 B% r5 }
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
3 f* ]& _9 c/ f! E$ s　　if(vl.y < rc->min_y) rc->min_y = vl.y;
  M# O6 a+ C4 h' k4 a　　if(vl.x > rc->max_x) rc->max_x = vl.x; 6 g- J/ f2 u! ~7 H: L, C
　　if(vl.y > rc->max_y) rc->max_y = vl.y; " w7 x  ~% E" E" H
　　} $ E) A" Y* p/ C
　　} 9 c7 `/ L- z- |  J4 h+ J

& D: B9 L0 i6 w( q, Z" J  l5 |. m7 k2 ?) Y+ a) `
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。4 K7 s' ]6 S- \( i, x: i4 K# H
' o1 N; S' O" ~6 E1 l
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：; o6 z  w' r* c; N5 M1 Y$ Q
( ^; d( R* B5 M# Z/ ~8 X' f- X
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;& @: R3 _& a- v, y0 u- u
$ ]/ c* m) B& h% p2 G$ y3 x
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
  J& V0 R; T  u0 D) T) j! x
! l9 ]- v* R) S以下是引用片段：
) r# ^/ j, B! k- S2 k; l　　/* p, q is on the same of line l */
1 f3 f: D+ b$ O: H, r　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
" t( |& r3 G% ]! P0 a# `$ e　　const vertex_t* p, 0 A. u1 e' i7 y' d/ K
　　const vertex_t* q)
! n6 V; A# z* X" `2 j　　{
- }8 I& e; N0 ~. C; p! \' u　　double dx = l_end->x - l_start->x;
$ _, a4 }' i, {! N, A/ h　　double dy = l_end->y - l_start->y; 5 A% d4 Y% `7 {
　　double dx1= p->x - l_start->x; 3 P% m. n; K, o
　　double dy1= p->y - l_start->y;
( g- i7 D: k4 _2 o' p/ g- E　　double dx2= q->x - l_end->x; $ N$ y, }. h/ `
　　double dy2= q->y - l_end->y; 9 G* ~$ J2 o7 V& e" j; e1 d
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
3 Z! G) z* z6 _: c; Y2 S9 \1 l　　}
% q$ C  T2 X# d3 A3 p& a9 e　　/* 2 line segments (s1, s2) are intersect? */ ) D- m2 P& k- V" |+ \: p1 _8 V
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
6 r: d. }7 _4 L/ O, C　　const vertex_t* s2_start, const vertex_t* s2_end) , Z; F; Q3 d) N6 q  [: a+ }
　　{
' u' y! \0 v$ W, Q0 X9 j　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && * g# \% u7 B1 c: o# C7 j
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; 2 `% }6 {( U" `! {  P3 n
　　}
5 y9 N; D- z, ~2 u+ @  T4 B# @" _% E- O* m8 f& {0 Q" y

; _; `) T( V; i$ @- v: B　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：  ^& b9 h- l1 I+ A: p1 Q
8 w1 x7 X% k4 l) A' g. r. K) [
以下是引用片段：
+ g# o* d2 _  r) r9 v) W! L) t　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ 8 [) }" ?7 E5 {! r" c% V* H
　　const vertex_t* v) 3 U" |  F' {4 S& Y  @8 }
　　{
& N* y& Q/ n! S. ^; B) q　　int i, j, k1, k2, c; ! V$ t1 o2 _5 A+ z0 J
　　rect_t rc;
' g8 I; s3 g# s/ |. W" d　　vertex_t w;
- Y2 C0 R8 y0 T2 ?! K3 j) M　　if (np < 3)   z2 \" |& y# e/ ~9 ?
　　return 0;
8 p8 s' A0 o7 A5 e) L5 [* I　　vertices_get_extent(vl, np, &rc);
4 g# X& G+ A( b  \- I　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) : z2 i! n4 y* W* \+ q  v
　　return 0;
6 h1 Z/ y9 U& l8 C$ z) @　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
$ p& T# _" c% _6 ~5 q　　w.x = rc.max_x + DBL_EPSILON;
: |# O8 O5 Y9 o4 U　　w.y = v->y; ! Q- `& x& A3 F, \3 M
　　c = 0; /* Intersection points counter */ ; l; x( j* f; e
　　for(i=0; i  2 N# u' \6 m) M! K% E. c+ P2 ?- x
　　{
9 ^+ q1 B# L' g　　j = (i+1) % np; 7 Y# f/ i9 x* B
　　if(is_intersect(vl+i, vl+j, v, &w)) : |5 M7 D7 m- T( E# @% a& u
　　{ " T" b3 Q' J! F, a: F& ^+ v
　　C++;
: D+ e& G' A0 ]. a8 z　　} ' m6 {/ s/ m2 I% p; j( S: t
　　else if(vl.y==w.y) / A# S2 Y* i( M( o8 N
　　{
" u, P& j8 `% A3 G　　k1 = (np+i-1)%np; ' X5 F* b- B# x, r1 \6 s
　　while(k1!=i && vl[k1].y==w.y)
8 g6 x6 V- v- K; A+ S　　k1 = (np+k1-1)%np; 0 }1 k" r$ \; c- q" _$ N
　　k2 = (i+1)%np;
  r( F; {* t: `; e! R% C　　while(k2!=i && vl[k2].y==w.y)
' _% ]$ F& Q2 w6 ~2 r1 `0 z, b　　k2 = (k2+1)%np; 6 F- v3 {$ _) D, B2 R4 S/ K- U0 f
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
9 \7 H5 q# \& t　　C++; ( ], c0 j" t) s6 r
　　if(k2 <= i) ! E' X6 e% s5 S  Z# @
　　break; ! ^. q6 P/ j, W) Y3 d; }7 ]9 X
　　i = k2;
, \6 |6 ^3 R3 \8 N# ?" r1 H　　} 1 J/ m- m0 [* @2 C& M; x
　　} 4 Y' \: p! L6 [4 t$ w
　　return c%2;
6 Q# u0 d' v$ ?5 o. s　　} 6 W$ ]& y6 d$ B/ x# U, z; a0 r* X) J
6 ?, n- n7 \% N

3 c7 C$ S1 m" w8 h. m! Q) q1 w7 [　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。