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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。' N3 E8 g. H& \4 Q6 Q: D* q
/ Q9 `  w3 {& ~! O# l/ D
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。, e" h& J3 L1 R3 H3 i" |5 u3 `

# n; s3 U* s- `2 g, T) V　　首先定义点结构如下：
. X- h  u1 l7 C) Z( `: v% K
, _, h" W" x2 @3 B' l6 c+ n( @' B/ x以下是引用片段：# q, @7 g) \& C
　　/* Vertex structure */ 6 j& ~# N4 r) C' l2 w  L0 Z+ [
　　typedef struct " P. l6 l. t: Q& Q9 s, ]/ [
　　{
: Q, z* f! ^2 O' k/ }' ^: m3 i4 O　　double x, y;
! p  Q+ U; r/ o* O$ r$ Z/ ]/ w　　} vertex_t; & Z) B( f* D7 _

) {  t& D4 [; K; N& u/ s
! @; J& O/ P, E* o　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：0 u1 [5 I' Z/ h* Y

4 j* J6 U3 r+ |! ]以下是引用片段：8 h" o  W2 K& H# P
　　/* Vertex list structure – polygon */ $ b1 \0 ?; }$ E7 D% D3 ^7 ^
　　typedef struct
" r: b! \5 O. |　　{ & i; _2 [1 E+ }  L) V7 O- D
　　int num_vertices; /* Number of vertices in list */   A2 x" Q9 h# c" J
　　vertex_t *vertex; /* Vertex array pointer */ / B7 n" e" ~" c- T1 a
　　} vertexlist_t; $ r- ?4 ^, ?) Q. K
' N6 i5 t" \. R( j" W3 f. K

; @; u+ m1 K/ a, r4 ]" n" H　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：/ [: J# N; m9 a! w6 D, X: I
, D8 O, H& w/ v+ T- S+ X* ~
以下是引用片段：
& M" O5 v9 P% S' w/ \) q　　/* bounding rectangle type */ & k: M* b( r7 N, _/ [; y9 N& b
　　typedef struct 7 Q0 z. c4 F+ V+ a4 ]" I
　　{ ! U! f0 I2 p  g) @" G  c1 ]# G
　　double min_x, min_y, max_x, max_y;
. R1 _0 p0 ?# b+ E0 z  p, u* E　　} rect_t;
4 \' l4 f/ {! z3 M5 P" A: a　　/* gets extent of vertices */
; @( F; H; l( q- f( R. g- J　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ * M. O, n. `8 s1 @; S
　　rect_t* rc /* out extent*/ )
) p8 g6 Z, H3 [1 Y: H0 o2 U　　{
$ U# A! y" m. P% B% ^9 ~5 B  k1 d: `7 X  m　　int i; 9 e* F. C2 V2 ]- `0 E
　　if (np > 0){ ) R9 b3 i8 U1 q* W' e: w
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
) M* T5 y" E0 n) {4 t3 d0 b　　}else{
8 V+ K9 D- v! w　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
+ X8 t- K( z7 \$ N8 i　　} 1 ]  q" o5 q. j
　　for(i=1; i  
1 U! v# _( |. U! v/ B$ Z: V+ z　　{
2 I% \( j5 t: q　　if(vl.x < rc->min_x) rc->min_x = vl.x; ; q2 J0 f% w/ u1 X( L: ~
　　if(vl.y < rc->min_y) rc->min_y = vl.y; ! v: h) K: |2 @; X5 d
　　if(vl.x > rc->max_x) rc->max_x = vl.x; # H2 a8 w3 w" ?5 x4 v
　　if(vl.y > rc->max_y) rc->max_y = vl.y; 1 O$ C; q/ D9 |8 f; q
　　}
( u( v; k) O3 X　　} ( H7 H" `: F! @% p

/ ]4 g9 E0 P# B; j. r0 M' _
' M; S& S( u; M! M0 b3 D　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。! T8 N5 C+ y3 p- L! f, X
, V* m; i- P  M5 F" k# n/ A
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
" y! u, j8 `, g. S2 z; _( K) n8 s3 t5 Y5 \
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
2 V: N+ S  j! `4 T( l6 u  y! K, m$ c  J7 M/ {- F
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
3 |; A$ h/ L" u1 v5 a% _  U; o$ m9 `; z  y5 Q
以下是引用片段：
* D; @5 L6 {0 ?, X. H, `% ~　　/* p, q is on the same of line l */ # V8 I$ q3 W- p, j" {- ]3 X8 {
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
1 ^! T8 Y% {% B0 z- T! \　　const vertex_t* p,
/ K; a) U* E, K4 e* H, o, R' a　　const vertex_t* q)
4 A4 F  k% G; C3 {: [& S0 i/ d　　{
- M# C5 E: n, E- ^- a0 e　　double dx = l_end->x - l_start->x;
3 S1 D7 K9 k7 u, k+ \+ J3 U　　double dy = l_end->y - l_start->y; 1 V, s0 J' I1 j
　　double dx1= p->x - l_start->x; 2 Q7 S$ [7 K$ G( `. p2 \  c
　　double dy1= p->y - l_start->y;
  J7 |- S7 o8 Q/ n9 L6 g　　double dx2= q->x - l_end->x;
& B5 r" t3 Q, g9 z% T　　double dy2= q->y - l_end->y; + ?( b6 n8 s5 J) |) Z
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); ( d, J: ?! x7 ~  @: t9 ?" E3 K
　　}
- [) d) s" ?  d2 m　　/* 2 line segments (s1, s2) are intersect? */ 5 H  @  Z8 `) ~
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, 5 }/ E& F( r# d- D) w6 \
　　const vertex_t* s2_start, const vertex_t* s2_end) 6 w, S1 Z8 [6 e& t+ t. b
　　{
# a, T! X% O( x& r　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && : G6 z0 c5 J- r
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
7 ^: ]( I' a2 F8 H  o# _: U- `. m　　} : _8 u0 C# p' i( }4 [8 g
& `: H; F0 G. @- v- w5 H5 E  a
  ?1 _! \/ n# V+ F$ A9 N
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
) g: v6 R0 T: O: W& }" A  y/ I% L
以下是引用片段：. I8 c( v1 r* }; b, D
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ : B  D2 Y4 j8 V
　　const vertex_t* v)
. L7 c$ J6 s+ n6 w4 k8 r　　{ : P( I' r: ?. |( P) }: g& J! N% M* Z" @
　　int i, j, k1, k2, c; 3 x# }$ w0 V% G( m
　　rect_t rc;
( Q" N' k. d& S" m& U2 E　　vertex_t w;
9 z1 b0 m7 U6 \! y) `( H! U2 @6 s　　if (np < 3) : l; G3 J+ U$ }' J0 d8 D/ l2 ^0 }" {
　　return 0; 7 D- X$ ~4 m3 C! t
　　vertices_get_extent(vl, np, &rc);
1 V, d" n: U6 p: b7 a' {9 m8 e0 g* }& j8 c　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
$ ?/ F( J2 O; k) y* J) e　　return 0;
4 D& J5 o! A9 U3 G0 [, L  N6 i2 _! ?　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
8 |# O. E: ?8 i- P1 O8 V7 j. v　　w.x = rc.max_x + DBL_EPSILON;
' d, y6 \- U+ i* C; ^+ I1 H. a　　w.y = v->y;
+ s8 g6 b/ F" e+ _, _) Q# }) g　　c = 0; /* Intersection points counter */ # T. B% c5 A, b& l
　　for(i=0; i  4 Z0 f. t0 e$ o
　　{
3 E, s6 L; W/ U　　j = (i+1) % np; 1 t  Q1 [  Y: {3 Y
　　if(is_intersect(vl+i, vl+j, v, &w)) & S: ]8 ^2 H% G
　　{ 6 U  ~* C( [+ ?) t3 N0 d* g
　　C++;
; S0 Z) V' O+ r* K) l　　}
0 r) i2 v& p6 p' q8 }$ U6 z/ m　　else if(vl.y==w.y) / X9 C: ~. P! s7 K8 X
　　{
3 m. Z9 L# D) p* N" U　　k1 = (np+i-1)%np; 9 M% D- ?9 G1 l
　　while(k1!=i && vl[k1].y==w.y) % G9 y  }8 p- T$ q
　　k1 = (np+k1-1)%np;
/ K( J) t2 j$ O: c! Z  t; M　　k2 = (i+1)%np;
5 S8 g1 a5 T# n' \0 F& c　　while(k2!=i && vl[k2].y==w.y)
% @! u# g8 _+ `1 L5 ~　　k2 = (k2+1)%np;
6 u; n9 g+ a8 ~+ c( k" @　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) 6 }+ ]  J% q$ B% ^
　　C++;
7 T5 V. S* ?/ s; H/ O+ i$ X$ i' I　　if(k2 <= i) & Z' S' e- P' V8 e: _9 T" [
　　break;
8 k. P0 i9 e1 w　　i = k2;
( u: ]7 f; G# W( Q　　} 8 Q$ x$ f; t; o$ L
　　}
: C5 @; s- ^+ _: X" {7 \6 K0 `0 W9 q　　return c%2;
7 r  `5 w! s) R! ^　　} 3 r( |+ L1 J- H2 ?
- m, X2 e8 L; w+ U) d

4 B' Z5 L' P$ h9 X: x6 a6 J0 X　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。