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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
" @# a9 D% v( k* P) @* Y
1 m2 @- m8 e  `6 V9 d　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。4 a. R8 s, N5 P1 l+ j
  p7 o' W6 z$ Y8 L* l) Z3 W6 p; M
　　首先定义点结构如下：
$ B7 M- [: a  h) d+ Y5 e  W  A4 |& G* _
以下是引用片段：( T# S  g4 d+ O4 v1 _& A" U
　　/* Vertex structure */ ) E1 `. A7 z% n8 n5 t+ M  z3 ?. S0 X& N
　　typedef struct
! U# b+ l  {# D　　{
4 V# e4 N! w2 A6 ^　　double x, y; ; p) J8 z) ^1 g- T
　　} vertex_t;
" J! [! H0 Q  Y$ U6 K7 ^# z3 B. q( N% z; i; ^" m( N

$ d+ j  }% d0 J9 {) T% E6 n　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
( L" R$ ]6 H# N. Z
1 r6 p# ?+ M6 K6 n/ q( [# \以下是引用片段：
, n) h5 T8 n$ }" A, [! ^$ \　　/* Vertex list structure – polygon */ 1 e( {* D6 o+ I6 G, [
　　typedef struct ! k/ R& V7 n& C' M" @
　　{ % b' M' u+ I7 @; @
　　int num_vertices; /* Number of vertices in list */
( [0 S0 H% O' M　　vertex_t *vertex; /* Vertex array pointer */ / U- k' ?1 t' o5 H( N, o
　　} vertexlist_t; # V/ M* x6 {  ^; z7 Y3 x8 h, }6 r
( L! ?8 d, o; l7 m+ I
: }" z1 i& t; {- _
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
- p+ [( {' v1 h& y- ]5 ~7 C" _! d; ^( S$ M! a& j- Y! n
以下是引用片段：7 Z4 U! B' w9 ^8 o/ T# P+ o' H
　　/* bounding rectangle type */
3 ^) f( O- k. V9 B9 s　　typedef struct
# Z2 b$ Y) K' Z8 I/ m　　{
5 P2 l' Z' `/ t& |　　double min_x, min_y, max_x, max_y;
8 A  o5 \0 I6 U' x; t# \  i　　} rect_t;
& x) x& x% c/ ^9 Q　　/* gets extent of vertices */ $ ]" L* N% r3 S0 Y" U8 F
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ ( |+ D" o3 k. n1 L' d5 e
　　rect_t* rc /* out extent*/ )
. i& b% {: s6 s7 M　　{
" {( \" G: |. o4 r% v　　int i; 3 h7 V( i3 C3 Z' ~: {3 m: i
　　if (np > 0){ ' I5 E/ U1 }2 b. e) N
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
( O9 E/ ^3 l( @* s$ C. @　　}else{
1 `* p, n/ G, S$ u# a0 f( a5 u: ~　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ 9 T0 ~% {( z4 Z' ~/ }" S
　　}
! K; U$ i" F5 C　　for(i=1; i  
  c' K; \1 g6 V* ~" r& Y3 @, D9 p　　{ 0 H0 E6 r, X, H; |$ O3 W3 m4 H
　　if(vl.x < rc->min_x) rc->min_x = vl.x;   K- o0 x/ Z' g
　　if(vl.y < rc->min_y) rc->min_y = vl.y; $ Z) K1 ~+ [  y! V4 W5 h
　　if(vl.x > rc->max_x) rc->max_x = vl.x; 2 g; e/ c3 M4 B8 q
　　if(vl.y > rc->max_y) rc->max_y = vl.y;
' x& k" {8 X9 {: F3 P: k& `( A* x　　} 2 Y3 A' W5 p8 n; M6 e. R
　　} 1 `( [  d& f- y: i
. g. k( Z- q- Z% F

6 e9 v2 B, W* k1 o9 {* Q　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
/ O, i( {# X4 J) g! Q6 n
. ]( m1 H; }2 N% C! i5 ?- i　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：7 l3 }1 p# z* Y( S0 z

( m( x- l# l! g5 ^* [2 V7 o　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;. F3 v  }+ p9 I

5 @) p6 ^5 x& ]" m% n　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
% j! U5 C9 m2 ?9 t7 b" T2 _
, J& i4 H6 J7 ], X- ~, d以下是引用片段：
& F4 d& n' `, L　　/* p, q is on the same of line l */
# \7 H, s# I* K% T　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
7 [, E2 D+ w& K! b$ l5 [, A% U　　const vertex_t* p,
6 f, m$ K& t8 p" l1 @) p　　const vertex_t* q) % P0 p. ?+ d; t# Q; P+ l8 {' f3 |$ l, j
　　{
# O4 R+ Y2 n+ J3 E7 {8 ?　　double dx = l_end->x - l_start->x; * T1 y) Q  m9 N' P) t9 S+ @- b
　　double dy = l_end->y - l_start->y;
, T, y- V/ \4 @# u1 b0 o) |　　double dx1= p->x - l_start->x;
3 B$ b/ A' j# G- L/ b* z　　double dy1= p->y - l_start->y; ) w( R$ S5 Y9 v4 @/ F
　　double dx2= q->x - l_end->x; 7 S! v6 R" L: i& V  M
　　double dy2= q->y - l_end->y; 1 m  P; n6 G7 n) K4 z. Y
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
! ?0 O& C7 Z6 X3 F3 B( ^　　}
0 L8 Q  l& h' a% @　　/* 2 line segments (s1, s2) are intersect? */ + [8 y5 c$ i+ k" w; N
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, % t4 j- b6 G9 i9 v" f1 n
　　const vertex_t* s2_start, const vertex_t* s2_end)
2 F$ L/ v# f* m　　{ " r1 O9 S: I5 x  i; ?! x& q
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && + E" v0 y# A/ x2 J
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
6 u. [& O# r/ k　　} + q8 U' O/ i! f8 j6 I

1 B7 v) P8 W* I
3 z( N; u2 H: ]2 }7 @: `& m　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：+ U' v+ B! a; a# R

+ R! G" o9 r# y以下是引用片段：" e( b, w1 ]& M3 H. S" n
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
8 u2 h- p: f" R! J; R$ u+ ?　　const vertex_t* v)
# i, N$ G: \1 v# Y7 H; L3 z7 v0 w# h　　{
1 A) |  a6 M% ^8 s- m& Y4 u　　int i, j, k1, k2, c; ' J6 k1 V( z1 m8 m4 @* V
　　rect_t rc;
8 G* V3 V4 V' _* B) _+ A; m9 u　　vertex_t w; $ f7 s) T9 H( R1 P
　　if (np < 3) . P$ w: y5 q5 Y4 _
　　return 0; . S- n, J7 p% g- S0 j
　　vertices_get_extent(vl, np, &rc);
7 \2 }& l, W% v' m　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) 0 A  b# X! [* X# e8 }- L" w/ Z
　　return 0; 4 w+ h; I/ S; K" V7 j0 o
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ ' C% ^: I# {( |8 E
　　w.x = rc.max_x + DBL_EPSILON; ( w5 l( U& {( O+ T+ X
　　w.y = v->y; # E* U4 d$ y7 Y1 [# ^
　　c = 0; /* Intersection points counter */
$ S% m1 [1 Q7 y* |( D) w* F　　for(i=0; i    C  P- k2 j$ X, O
　　{
( m& y: l" q; G. C  v　　j = (i+1) % np; ' g. Q: ~8 {$ O* ~" U+ Y: h
　　if(is_intersect(vl+i, vl+j, v, &w))
, D3 @6 q  Z$ _" w9 d2 f  }# }　　{ 0 O7 M( z: p/ p& X& @
　　C++;
3 i7 N, s/ V8 ^( T  j6 F　　}
8 ~0 C- ]* r; X( a　　else if(vl.y==w.y)   w3 r; E- E4 F; _
　　{
6 t) @( [- i, b6 k8 o) h& R4 @　　k1 = (np+i-1)%np;
5 X% k! S/ o  S/ J/ P# U% z- c# u　　while(k1!=i && vl[k1].y==w.y) $ v/ ~& R* O6 S& _6 I' P
　　k1 = (np+k1-1)%np; ) Y. a/ n5 h9 V! F0 |& b" b& [
　　k2 = (i+1)%np; 0 z0 j7 S; q( J; r/ d0 o* y
　　while(k2!=i && vl[k2].y==w.y) " T  D& F/ Z. ^, d1 D8 r1 w/ w
　　k2 = (k2+1)%np; ( b2 Z4 U, J. G! {1 m4 [. t
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) ! K& @0 [# @3 k7 u
　　C++; ; R: ?3 K2 N' \. c7 _; |: \
　　if(k2 <= i)
' x2 }/ h  B, G% _  x　　break; & a. o# g3 i/ U1 k0 ~9 K! s
　　i = k2;
) a, m2 W; |' J) q. x* [　　} ( V! r4 {+ ^; q  b. @+ `0 ?
　　} 8 ]! T4 L1 h  D( _* k. W
　　return c%2; % {# j/ d- T# V$ R6 v+ m. @  C
　　}
7 C, {$ @9 c- j+ H/ a9 G* e7 e  \& n* ^* C
' m8 k( @, h2 Y- ~
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。