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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
) N4 c2 r9 Z: K5 O- j$ r$ g
" V" ]' Z  w$ X# m　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。( M$ L) r% Y  l' |6 I3 V

8 F8 m7 J6 b( W# ^# L　　首先定义点结构如下：% I, W- T6 S1 {" Z6 `2 M
% v8 A7 o; N: b. y" }/ c5 Y
以下是引用片段：
, N! ], N( K' [* @8 F　　/* Vertex structure */ 0 m  W' U  J" z4 K
　　typedef struct
% U  v; @/ J$ h- V1 R　　{
4 J5 f8 U" b/ ~　　double x, y;
+ M' b6 m0 X( l$ d" R$ c: w- F7 S　　} vertex_t;
; ^1 Q, I! A- ?6 O# Y; ^8 T& n2 P* N* n* O4 C1 y

" K( y& V9 s1 I4 h: G( b1 ]　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：( Q$ B% `5 h; Q

0 W, y4 ~8 C2 ~$ D2 r. {以下是引用片段：2 x3 L% i# A2 |; y
　　/* Vertex list structure – polygon */
% [) S/ E% _( p, A　　typedef struct
+ b% O9 l$ s' S8 C　　{ & ^4 p* f7 E- {4 A3 H9 n1 U! m
　　int num_vertices; /* Number of vertices in list */
3 L0 |. m5 s) }% v　　vertex_t *vertex; /* Vertex array pointer */ + J# ~0 U. r) o3 l. J' q
　　} vertexlist_t;
1 W& R' `0 u& C0 D
# Z3 S- C6 _6 t8 p: Y: o1 D! w8 ?  w. u" P  n. ?1 r" Z% E. T, w2 P2 n
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：, x1 M4 t) ^: i" p: P" x+ x; X

0 a! d- s/ o+ i, E$ N( G0 |以下是引用片段：  k. N) c+ w1 q4 W2 T; U
　　/* bounding rectangle type */ 2 d$ d  H* s! E% v
　　typedef struct
  ]- Y" [! v7 w　　{
7 R' I/ b7 k. b! L　　double min_x, min_y, max_x, max_y;
5 m5 @4 q7 N3 ]2 f5 I　　} rect_t;
, L& X( B8 m" X3 W　　/* gets extent of vertices */
; @- z6 f* y2 g　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ 6 g; E5 n# g5 ~5 k7 _. \: U  A4 k
　　rect_t* rc /* out extent*/ ) 0 s8 u+ G. `5 U% n& n+ c5 C  l
　　{
8 w+ a+ L( f/ Y0 g　　int i;
% Y1 D- x) D# ?9 d9 Y( H　　if (np > 0){ " o, D5 g, d* k( [* E
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
9 m, M4 i9 Z7 G8 J+ c! M( N　　}else{ ; H( C# m9 h8 n2 q
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */   k- V# b( F8 W0 m2 }. u
　　} # r0 N# e8 q) Y# {$ z' R, o
　　for(i=1; i  # W' x4 ?  T4 G' p
　　{
1 M4 n% i5 X& ]2 M3 T　　if(vl.x < rc->min_x) rc->min_x = vl.x;
2 c9 S8 S) d5 N7 W! i　　if(vl.y < rc->min_y) rc->min_y = vl.y;
- C) c* n0 ]( h2 `+ u; z% d　　if(vl.x > rc->max_x) rc->max_x = vl.x; - I1 p) L1 j8 \& h% R0 E
　　if(vl.y > rc->max_y) rc->max_y = vl.y;
" l* y" K( y) d" ~. F) O　　} 0 w) V. l! J$ ^1 L+ H6 P  R
　　}
; g8 O7 n7 r. S+ T% U; h. F2 w, N. ?3 J* S, o  F9 t
& _/ {' @( D) p
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
( ^$ y3 d, C4 S1 ~0 [5 ?/ h5 Z
% e) q9 e! S6 ^+ D3 p. s　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：# E  a: o/ O: |& D: x6 q
9 I7 a7 L9 g$ O9 K0 y5 f
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;6 g* C( ?- R( y4 ~* g2 e

8 r/ n" k1 `0 e+ q0 i# w5 S  z　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
% t1 @) O5 y+ G: q' ~; _* B( s8 y9 w5 _! q4 j( V3 i
以下是引用片段：
9 I6 O6 Z! H% {  \4 M( q7 I! n( k　　/* p, q is on the same of line l */
6 \& t- m5 L- e9 y　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ ) d, `& D; W9 J( \2 P
　　const vertex_t* p, ' ]/ J! g$ Q( S1 B$ d: g/ S; R
　　const vertex_t* q) 7 U/ F  P' q, @! B1 i& V  v
　　{
( Q' g) |: L6 H! p* E　　double dx = l_end->x - l_start->x; 8 S0 m: H$ T1 U- q7 Z
　　double dy = l_end->y - l_start->y;
* O" p/ L, w; f- h2 v　　double dx1= p->x - l_start->x;
- J7 Y* y/ [3 {- @0 L7 V　　double dy1= p->y - l_start->y;
5 T5 }# E: T, G& X　　double dx2= q->x - l_end->x;
7 {# P. o* P/ `1 K+ _1 v　　double dy2= q->y - l_end->y;
; x* h2 w- [6 k) T2 a2 `　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
4 Q0 A% }- p' \4 c+ o　　}
8 Q3 V% ]+ N4 {　　/* 2 line segments (s1, s2) are intersect? */ , ^9 q! c' r* t1 i! V, p+ }1 }$ t# d
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
3 h: f4 L3 m' }6 w$ B% ]% i　　const vertex_t* s2_start, const vertex_t* s2_end)
- w3 P4 [+ k4 B' o/ {( o　　{ ; l4 d8 o+ A( H* k3 O( v
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
( e! C2 m  h; E　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
) P/ Z/ X( u1 Q8 e9 Y" f　　} ! N1 s4 _3 ~  u  L) M' O( _
* Y+ k* @$ S  j9 \4 e
5 ^' R0 z  {0 @( x$ w6 O
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：/ J1 D' |3 K) W& Z
6 D: {! v7 R" O" P  |
以下是引用片段：
- v. s; J+ @; d* G" u+ j; ^2 M　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
' b' `* w- q4 i! [3 _) ]　　const vertex_t* v)
8 _' X* h2 }1 f; d+ m# T　　{
* B- U% [9 L9 L　　int i, j, k1, k2, c; ' g% N( n2 `7 Y% q
　　rect_t rc;
# {2 x2 ^  Y6 G  ?3 W) ]　　vertex_t w; / m4 h0 j! d+ V) _3 B, i& }& d
　　if (np < 3)
" _  q9 T: p6 v7 r( ?9 Y/ o　　return 0;
4 i8 B6 s: q5 j* A　　vertices_get_extent(vl, np, &rc);
. T# g5 L  V0 W% [　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
& X" g" Y* |* X) i　　return 0; " k  l' l9 ~& p- f/ h
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
+ `0 I2 c6 |- Z* w4 T7 p3 u$ p6 t7 Q　　w.x = rc.max_x + DBL_EPSILON; 0 T" s. }- B) O  E* G& }
　　w.y = v->y;
3 u& V$ `9 D: ]　　c = 0; /* Intersection points counter */ : J4 K9 z0 l* R) A
　　for(i=0; i  
1 \" J6 K7 X7 r9 Z+ j7 @7 D# c　　{   l' v4 F1 B  b1 q- {/ ~
　　j = (i+1) % np; 2 A& x  w! h" p( M  T! h! T! Y
　　if(is_intersect(vl+i, vl+j, v, &w)) 6 ?: y( F2 S$ a' n6 ?* M# P; \1 O
　　{ ) x* x; q! P; o  E
　　C++; . S! ^0 p! e; J
　　} 2 P% J7 ?9 i, Z4 E- ~: J
　　else if(vl.y==w.y) ( i! K, m9 ~7 z, {  Y2 e: x4 x
　　{ . f9 k% L3 g9 T, P: x
　　k1 = (np+i-1)%np; + `% }* @7 L7 t
　　while(k1!=i && vl[k1].y==w.y) ( F& z3 V) T, c& P" N' Q
　　k1 = (np+k1-1)%np; " [/ p0 S) Q2 h. K: e9 G( R" b
　　k2 = (i+1)%np;
& t5 j; q0 w5 g. E  u, g; m　　while(k2!=i && vl[k2].y==w.y) 3 w& ?, h  J% {& k; Z5 x9 y
　　k2 = (k2+1)%np;
: R4 j& Q" w, B) S　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) + d; ~3 Z1 X% |- {4 P
　　C++;
  t) z" N7 r' ?8 K9 v0 f" o　　if(k2 <= i)
9 F4 g  X: V' ^7 m3 P* F9 g1 P6 ?　　break; $ y& o0 U3 X$ s& F: u5 \
　　i = k2;
- `* ?) L4 {. [& _# H" o, ^　　}
( a# C1 y2 v* J1 E3 x3 n+ T　　}
1 F( u: x. I8 W1 r, h. w( V: H　　return c%2;
+ d. G, q- y5 V, [+ g　　}
1 X7 }& M4 z- A' G6 U8 o; ?1 |; v7 @8 |- q. @& {( @

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