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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
" v2 ~; @7 D2 f0 V) c7 c$ i: i. G' P& s6 v
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。4 a. k6 h5 [6 M% V1 o& _9 F
% d, ^' u: e4 s) r, m7 G& M% ~
　　首先定义点结构如下：
- S( h9 p( `& @- k# u" w/ A2 B) G) U6 S
以下是引用片段：
! f1 W" d2 |0 ~) H- i4 L) n& _　　/* Vertex structure */ / O, {5 P0 u; ^3 @$ B
　　typedef struct * j' z- E* X; e5 y. M) B6 E4 B
　　{
* V1 R) J* |; W4 J  `  a; v　　double x, y;
; b: v: _6 s2 r8 y: b. s　　} vertex_t;
8 `! v, @: w7 a  n6 h% q; r' ?& m8 z% W# ]- F1 K  Y! T0 p8 I$ X

+ }4 c/ r" O: \+ z- ^　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：( J/ X5 N% Q: y! W& N. D5 B
1 v! Q/ o6 ]) `0 x
以下是引用片段：
  E) W, `" @& Y( S: v  i3 s　　/* Vertex list structure – polygon */
" Z0 G* j7 y7 X4 v. F　　typedef struct + T9 r" t) V. B% C2 \0 \
　　{
' H' N7 `7 O$ |　　int num_vertices; /* Number of vertices in list */
$ f  n! F0 E4 @0 [6 r; ^$ M- V　　vertex_t *vertex; /* Vertex array pointer */
0 a# I& g8 X. y7 a5 _4 o　　} vertexlist_t; 4 a+ _. h- |( A7 l& ?
5 t. K6 H. z7 X, d* _# i: i
1 ?; C2 Y- ?9 i, T  S
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：' b3 X  r& H2 ]7 G  t

' ?1 M( K& M3 x2 |以下是引用片段：
$ Z) {" F; J7 ?5 G, W　　/* bounding rectangle type */ $ A2 a) W. s$ J9 v% T) x
　　typedef struct ' {; R# L8 x) \2 T. k; z: P1 {
　　{
6 P( `! j, J" r+ [, Q1 f3 d! ?! C　　double min_x, min_y, max_x, max_y; 9 i/ |9 w+ ?! i7 ]7 o# ~
　　} rect_t;
8 M' k8 `3 f, A0 g0 w　　/* gets extent of vertices */ * Y/ H% {4 Y: t3 n- {' Y% D, F
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
) O6 a2 b3 k% U0 g　　rect_t* rc /* out extent*/ )
& b' s& t! O) v) u1 _+ @　　{ % M) u5 T0 e0 M
　　int i;
7 J" k9 w$ l" {$ {! Q" X　　if (np > 0){
& t! V2 V. o3 w　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
, x$ G, x1 R7 c* E9 V! h+ X/ s　　}else{
$ [* b4 ]4 x) ]( v' s% ~8 e0 V　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
- n* \/ E% ?+ N  a5 L& `/ x　　} " K4 [. m# m2 t: }$ v
　　for(i=1; i  
7 L8 o( G: U4 l& b　　{ : W! g6 o3 K" c1 z
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
2 w: Y. f, E6 q* I! b9 E　　if(vl.y < rc->min_y) rc->min_y = vl.y; , z+ {. Y  n- O$ a2 k) r
　　if(vl.x > rc->max_x) rc->max_x = vl.x;
/ j; N4 n9 R) a& _　　if(vl.y > rc->max_y) rc->max_y = vl.y; * _( o6 ~8 ~5 E
　　} ; Q, ~- W# x8 Q7 U
　　}
: T9 N/ V$ Z1 m
5 k: d" o1 A2 M9 P0 r
0 l  `* C- i3 e3 c1 m$ u5 x　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。- S7 E% H4 X  _5 D/ K5 J
0 X" T  Z: j. _# M, ?8 m! m" ?
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：0 v9 h* v( K- C  \6 F- ~

4 x- f1 V! w4 |3 P( B　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;. Z; w& T7 g/ i7 b
7 Q  b3 q9 z( E; \8 f* e
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;8 s  d" q$ p; D# e) A
! E& X4 B4 O" B' k5 e5 x
以下是引用片段：0 r. ]: u9 S) d+ U' n' T
　　/* p, q is on the same of line l */
  \) O% f  M; p  [* z! j  R　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ " `- G5 Y" n/ o1 _; `6 F8 ]3 @7 f
　　const vertex_t* p, / B% Z" J  W/ r  F$ |0 F
　　const vertex_t* q)
8 F( A9 Y- F, i& i; l( E/ o　　{
9 H* `# T! U* Z- P　　double dx = l_end->x - l_start->x; ) K. {* {% m  ^# v: _7 I( t9 i
　　double dy = l_end->y - l_start->y;
! ]1 K$ Q3 V: V8 {$ m　　double dx1= p->x - l_start->x;
1 }4 K7 B1 `4 ^' M% u( C3 m, p　　double dy1= p->y - l_start->y;
9 }( G- y# o: Q! v. d3 J2 N( j, r　　double dx2= q->x - l_end->x;
1 ~0 I6 F6 f- o% `( t4 x. ?　　double dy2= q->y - l_end->y; * }+ `; h* c& ]) i
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); ' z% b% |1 S) H4 J" A0 G  p
　　} ' `7 @" ]3 \' n& T+ a9 ^
　　/* 2 line segments (s1, s2) are intersect? */ $ _8 x/ d- ?5 N* \) S9 b
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
- ?) a' r6 t# s$ `1 \9 H' s" j　　const vertex_t* s2_start, const vertex_t* s2_end) * i4 F- \! S* q! ~
　　{ * [1 b" D/ T0 w8 ^
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && * G& ?) D: b$ P! [$ @# z
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
: f, y) h8 o6 s% `0 B1 f" o: i　　} ( I! \& T. S7 I. l/ C6 H* O; F

, H! r6 [! m$ Y$ }1 s- O9 R. {2 `7 o6 m/ J. ]' i
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
* C, L* m* h! X6 A( W  [( ^( K- X  m3 |9 U2 N8 ]8 _
以下是引用片段：
& W9 C1 c) j: i& a% [+ P　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
% A4 A5 J8 W$ r8 o% r　　const vertex_t* v) 6 d3 o  @8 Z+ p% ^
　　{
1 j1 a0 t" w' v4 D' H+ q& X8 ?　　int i, j, k1, k2, c; 1 o8 _8 S8 _. [5 x
　　rect_t rc;
4 x0 p( ~0 t) e7 ~" m; T) M6 |　　vertex_t w; + g) z6 _3 f( B4 k
　　if (np < 3) , I; ]4 C/ T# A* B
　　return 0;
! |* K( m8 R# C( ], @) k# y　　vertices_get_extent(vl, np, &rc); & Z- y3 q+ p( R; v# q" j
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) 2 \+ j8 X2 M) W1 M: h1 _9 a8 V0 }
　　return 0;
: Q! f2 R. {( s) `　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ + }' _& k# w! d( w
　　w.x = rc.max_x + DBL_EPSILON;
5 O: U1 V1 a1 D! x0 M　　w.y = v->y;
4 h. P9 {- R# U. s9 W　　c = 0; /* Intersection points counter */
, `/ Y$ ?: a) {5 z; u) ~　　for(i=0; i  . @$ ~0 j( O. ^& m0 [
　　{ ! B& u6 i& g8 e
　　j = (i+1) % np; + `! T4 L" B% L+ S' v4 q3 }
　　if(is_intersect(vl+i, vl+j, v, &w)) ' ~2 G, ?  T+ w- `' C
　　{ - G1 ]! q% f' N: `" ?
　　C++; 2 L3 ~/ x+ S* ~
　　}
, _& c: |# ^% g- M0 Q5 O3 k　　else if(vl.y==w.y) " {5 G' @& N1 }7 q( [  x' W9 n" M
　　{
& p7 t, }" }( E# a　　k1 = (np+i-1)%np;
% e1 Q7 g- x3 q! U, J$ Y　　while(k1!=i && vl[k1].y==w.y) 2 ~6 ^9 E4 R' _) `
　　k1 = (np+k1-1)%np; - k' q% ?# n, u
　　k2 = (i+1)%np; 8 d5 a( Q1 g/ w" }, Q1 ]7 o
　　while(k2!=i && vl[k2].y==w.y) & q9 r: F; D# C, I4 M
　　k2 = (k2+1)%np;
9 G/ a: h3 j2 T0 z* ^* ^5 S　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
" [2 z$ a" w+ k: ?　　C++; 1 U) X  D: d' g% r4 N
　　if(k2 <= i) 1 S2 i! {9 Q! I9 H( S4 Q& i) q2 B
　　break;
, H9 @( o% {9 K( |1 F+ W　　i = k2; . d/ D0 k, M" T& T
　　}
+ {: h$ @8 J/ C, m　　}
& W8 i5 l/ }+ F/ |8 {- G9 y( L　　return c%2; , p) D9 D2 m) g# S% C& g
　　} 3 {3 a. }: _) U( Q$ p1 t+ J: `! Z
/ m* \( q5 G' R& H* J" @! y

& M: l, I4 W2 Q5 k5 H* a　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。