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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
9 l6 I2 i4 o8 D' B) s6 N% e
, N% b; x: Q  v+ M　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
0 F2 h5 z' y( f% e3 `* f4 ^4 L
4 j+ @3 {8 C! v　　首先定义点结构如下：
+ f( {) m" A: d8 j/ z# j1 ?9 l  @- H5 @* Y7 X
以下是引用片段：
9 I" [" d3 H: ?　　/* Vertex structure */ 3 Q5 l( |( G: s
　　typedef struct 6 h  h+ f9 }" d# R+ v8 N
　　{ % R, m7 n7 q1 i( z
　　double x, y; - |4 n$ x# W" L" E) e
　　} vertex_t; ! \- W) W6 K; ?5 w3 t& t* Z

+ W. G0 m# Y* K1 j2 N( k
: ?  m$ [' b& O; G　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
( ?% M. T5 ~: K6 O% d
7 E$ R( l% n1 ?以下是引用片段：
8 w- g& \. D' h+ S* x3 q: J　　/* Vertex list structure – polygon */ & K# W3 R3 E( f1 R1 y6 L1 E
　　typedef struct
3 s0 s+ E  V/ Q% |5 h+ `1 B　　{
+ ?1 ]; g* L9 s' M( D# J6 k　　int num_vertices; /* Number of vertices in list */ 2 b: j0 [2 u( \: T
　　vertex_t *vertex; /* Vertex array pointer */
( m% _+ ~/ d6 r$ T* t　　} vertexlist_t;
! l3 m% V9 [+ h# e8 |  a/ ?6 x, O6 t2 o4 D2 V+ K3 i+ @8 P
% R* K% \/ `5 }; i( g
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：$ Q4 F1 E& `0 f) s" Y+ S6 I
+ p. o1 M3 U4 m' m1 u2 o; W" D
以下是引用片段：5 Z2 f; t$ s9 k) P9 h! r$ U1 l! M
　　/* bounding rectangle type */ # R4 y6 }, ^+ [/ ^9 v" L2 @" O# c
　　typedef struct
. I$ ^6 ^" }$ k( s* _( L6 q) Y+ h　　{ " P; F/ |% L* `3 x0 @7 G) d
　　double min_x, min_y, max_x, max_y; 9 f& i3 W0 }) f3 D5 p+ s/ T
　　} rect_t; 7 _1 R8 j) q. l. r: W
　　/* gets extent of vertices */   J1 n9 ^- \+ T$ ]! k' o; k8 h
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
$ v  P2 W( b( J　　rect_t* rc /* out extent*/ )
. n  F3 Q; ?! R7 N& g6 T+ _　　{ ' ^; h0 C2 M$ K' z
　　int i;
$ p8 i( G+ j% k6 v2 g　　if (np > 0){ ) p0 ^  V; ]. D, g3 N0 }/ \
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
) s2 n! z  g* A6 ?" j; b　　}else{
7 Z5 \. P4 c0 f6 K  h) y" D　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ : V, J4 n" M9 h. g5 p) F
　　} 8 _3 ?% a" R+ T0 F& F. f4 F0 [
　　for(i=1; i  
1 k' e- m, `$ L7 q4 c　　{ 7 U, _% B" p9 l5 E* v) E* h# W8 D! E
　　if(vl.x < rc->min_x) rc->min_x = vl.x; 5 K/ Y* a3 b* }+ I9 W" m. O
　　if(vl.y < rc->min_y) rc->min_y = vl.y;
# N! S" E4 F7 m: w  E3 G　　if(vl.x > rc->max_x) rc->max_x = vl.x; # a) y! g" p. O9 N  @; w6 ~( L
　　if(vl.y > rc->max_y) rc->max_y = vl.y;
8 Y' G7 y4 W  w0 Y2 y$ M　　}
3 \. K4 b4 X* |2 ?8 x, u　　} & p5 n$ f2 S# Q, D0 k6 G
, o$ B+ v! H' l9 E+ P& P7 J
* O  j, O3 n! q# R, v
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。* S5 D) \) D: x. x' t  e7 u
5 g% ?% w: d1 M
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
$ v" U. r  ^  w, Z4 D2 K, I4 `( Y
2 e0 E0 v( R) v% c, i$ S　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
! d, R! @, f* s* G( ~
4 B" o& m/ T  a; n* C2 ?1 r　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;, W7 z. A( }' ~6 I. _1 V% `) ]$ O
7 _/ g: @" p, u8 O# o
以下是引用片段：/ r$ q8 r2 }  Z3 ?' D8 y! s+ J
　　/* p, q is on the same of line l */ 9 s/ @8 U; a/ B1 \+ Q' C
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ * q( D/ ~. @% r. u$ s
　　const vertex_t* p, # v5 x5 V& L5 L" |
　　const vertex_t* q) 2 s$ e. F- k1 g1 o) K
　　{
# `: |1 x& p" K) v  {$ ]　　double dx = l_end->x - l_start->x; * O" S8 |" W/ R# h. J
　　double dy = l_end->y - l_start->y;
5 `, p% x3 A# d; X　　double dx1= p->x - l_start->x;
2 z* z! m3 k" k! h3 H! j　　double dy1= p->y - l_start->y; ' h$ r0 X6 t( B! T
　　double dx2= q->x - l_end->x;
1 J% d- y, n1 V  X. N3 s　　double dy2= q->y - l_end->y; ( i# P# L8 W1 C9 O* n
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); ( D( T- D/ N! N' x: A& M
　　}
1 V( X+ N; L$ j' T' Y　　/* 2 line segments (s1, s2) are intersect? */ 3 i$ W/ W1 R6 u
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, 0 ?+ b" F: u/ r5 H8 K8 f3 U
　　const vertex_t* s2_start, const vertex_t* s2_end) - f/ ~: D8 p+ F8 f6 c5 u. _
　　{ 9 f7 b  D2 L0 b0 w8 e
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && 4 k" j1 u: U$ B1 k& {2 s
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
" f- t: v( ?% Q$ B" l' L　　} + s/ ?+ Q5 j8 h  ^; o
4 s2 W: ]6 Q' \( R8 H7 ^

6 [& `/ H9 {. `8 p, j) W　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
5 }9 b4 k' C8 H7 u' R9 Z0 h8 c8 E4 _4 a, b. i1 v
以下是引用片段：
& @  t. @7 P; {　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
" ?/ _, k6 X6 U6 e　　const vertex_t* v) . D7 l) f+ p2 i$ P2 Z4 t9 A1 n
　　{
4 G$ T/ |7 ?3 ]% u- [5 k6 [# w　　int i, j, k1, k2, c; 9 y, g/ A* `1 _
　　rect_t rc;
3 G1 N* n3 d9 }1 |5 W. t　　vertex_t w;
9 V; a1 ~) d2 R/ H) V　　if (np < 3)
& {. o1 |* A$ N3 |! a: Q7 d　　return 0;
8 }# P& r# }! Q4 x# N' r  c　　vertices_get_extent(vl, np, &rc); * [) ^- w8 O* Q! C- U: j
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
3 g& T. N& f$ l' z- v6 _; ]　　return 0; % r! r! h2 P5 n3 L( ]
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
" B3 B; s4 j- O. D& `　　w.x = rc.max_x + DBL_EPSILON; 0 J! o" X+ u! ^$ i( k  Y8 {& \
　　w.y = v->y;
; ~/ m5 e# T( M( t5 X9 `　　c = 0; /* Intersection points counter */ 2 m* {# X+ v$ f: ~7 ^6 e1 n8 }( ?$ ?' A
　　for(i=0; i  
8 _+ ?* y1 s5 b2 ~& V6 Q- `3 D5 F　　{ . b, G, _3 ?: Y( d2 F
　　j = (i+1) % np;
! S7 H/ h9 D: U$ f　　if(is_intersect(vl+i, vl+j, v, &w))
8 ~  v0 Z+ W& A* i, n; f　　{ # B0 Q8 _3 \6 a" F: y7 O8 u1 w* b% }
　　C++; 8 ^8 l5 A; k  d) Y
　　} 2 U. ^5 f/ O  j) r
　　else if(vl.y==w.y)
0 p" b3 K( g- S0 t　　{ / J) p7 l5 I+ R% G1 C% m' N) O
　　k1 = (np+i-1)%np;
+ y. I' A) J/ d% W. |　　while(k1!=i && vl[k1].y==w.y) ) S7 P$ \, i1 P8 T/ x# Q
　　k1 = (np+k1-1)%np;
$ }5 g4 {1 O# \; D* t/ p　　k2 = (i+1)%np;
4 G/ s& A2 h7 c; @7 q" u- V　　while(k2!=i && vl[k2].y==w.y) % v7 B4 l% y! N! B' b
　　k2 = (k2+1)%np;
/ |( ]# C/ u5 V　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) 1 B2 s* V% v0 Q( \0 w6 C% x
　　C++;
; q: H& H6 {% {) L3 U　　if(k2 <= i)
/ d0 n9 Z: n5 L, V/ \) D. O% q5 O　　break;
( ^0 \1 L: L' M* _　　i = k2;
. G( q5 m# X( f7 i　　}
) [& D/ S/ o- K$ g, N　　}
. w+ S6 M. j1 y3 @9 @　　return c%2; $ W. X9 r4 d6 P! p( G3 J) F
　　} : N; P1 `4 B2 I5 o$ C4 {2 V  y3 v  d

; }/ }6 V( D' b: p6 P9 d: B+ B+ h5 N8 @) f* z6 i
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。