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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。; v4 U. y' ~2 u4 w+ K6 H
/ g, L) ]6 J8 i0 r* ~( j
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
. ~! x1 g7 X3 G9 K) N
$ G& X/ Z4 e) G& M& @　　首先定义点结构如下：
# Y5 {# x6 x9 q5 Q% u. o/ f& g. q- u0 H& D
以下是引用片段：: t0 w! p& G% j+ T6 \
　　/* Vertex structure */
& I" g. N9 M3 }2 _# J3 e　　typedef struct % Y6 S! }3 q( @+ A9 z5 ?* m2 j
　　{ ) h, s& l( \8 ?; |
　　double x, y; ! Y% i+ u( q% o* Q3 T4 g. @: D, S
　　} vertex_t; & a  F/ e" s! J. e

0 ~  Z4 A( `0 G* d5 }/ j# Y. D! U6 @. ?1 p3 a' F2 l( |
　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
7 H# ?' N, U! P7 t2 M  g: u2 d8 h
6 x8 a& h& f$ k, f4 S; a; g5 }6 C3 J* X以下是引用片段：
" N# n6 L' e4 V# }( Z　　/* Vertex list structure – polygon */
# {, {# _2 _: u3 f　　typedef struct : o% ~; C0 l: t* |( Q
　　{ 9 u9 w* Z4 L$ j* R# ^
　　int num_vertices; /* Number of vertices in list */
+ e5 z& F, O5 J) _* `- H* F　　vertex_t *vertex; /* Vertex array pointer */ - L: l; u( v' n
　　} vertexlist_t;
# d8 F9 c2 ], x' y
. R6 m8 i! }8 d7 H; Z' i# M) V6 N5 r' M' }' d
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：. Q" ?& k5 `- a

0 z" ]# f1 t% ?& j. x以下是引用片段：
/ x$ z+ H+ p2 w: G　　/* bounding rectangle type */
* u3 D; h# [9 c+ \9 L　　typedef struct
7 Y; d! E8 A* D% ]8 G5 b　　{ ( r: Y6 l) P* ]5 m, j
　　double min_x, min_y, max_x, max_y; % F* [$ l& {, f5 b8 I* a& b* u
　　} rect_t; 4 q8 n0 w* N  T# m# u
　　/* gets extent of vertices */ ( I9 b% C* L6 M% S. G# e  s3 K8 a% C
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ - g9 i* k, I5 \7 x0 R8 i8 N
　　rect_t* rc /* out extent*/ ) 1 w2 \$ t) }1 E. q2 t, b9 E6 m
　　{ # V: s( Z9 D) c6 f  K1 J9 i
　　int i;
/ E1 T  \" {7 Q- [　　if (np > 0){
* t& r7 J0 A+ B1 f5 Q) t9 d　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; 5 q- a6 m* G' l& S& J- Z4 J! d, Y
　　}else{ 5 y% ]$ c0 H+ P; N( z$ n
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ : |. h- P  `" @( @3 J' s4 h
　　} 9 I7 L+ U' x9 N6 D* X) g* m8 q
　　for(i=1; i    P: @' w4 l" f8 G/ t( z" ?; ^% v
　　{ 8 v/ t' o' i) `  h
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
5 L$ _) u2 c8 d% j0 d$ ?　　if(vl.y < rc->min_y) rc->min_y = vl.y;
1 P( Q" J- j" n7 P  c　　if(vl.x > rc->max_x) rc->max_x = vl.x;
$ ?# Y$ Z7 N( C. J5 E( [2 ~1 s2 I　　if(vl.y > rc->max_y) rc->max_y = vl.y; 4 g. V! A0 m7 S/ J
　　}
' b# A7 ~& A  P9 x5 @　　} 7 F7 G+ c& M- ^6 a) Y8 R) V; Z* `

/ w0 e0 r" W; J- B5 ]% B$ c) p1 }/ x$ i5 x; c
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。( m2 Z. e$ V. G/ c

6 s% g$ K: w  n6 u, h1 q3 U　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
9 o9 {. ?' l$ {  ?  ]9 L
, c4 y  e1 ?! u　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
3 e  z) G1 F! @, {, Q: N9 v' j" r& E. S7 [. ?/ W# c3 \
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
4 m- D$ l. m5 e/ _+ V- v' Z1 G+ P/ b, Q& p+ X
以下是引用片段：+ M, W. M3 M/ e' k& E/ ^" a! N. B
　　/* p, q is on the same of line l */ ' q1 X, s5 s& H( V9 B* F
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ 8 o9 g6 s& O$ T& [, r
　　const vertex_t* p,
2 p) L! ]  d$ ~5 L, P　　const vertex_t* q) 6 V. ]! h3 Y2 A9 T
　　{
) H. W1 A: g1 \6 e7 `4 [　　double dx = l_end->x - l_start->x; ! `" z3 W. o2 T
　　double dy = l_end->y - l_start->y;
+ V! [& O6 R/ q0 u% A3 R; x　　double dx1= p->x - l_start->x; 4 b8 n3 Q( k7 i# M4 ]  u" c( l( i
　　double dy1= p->y - l_start->y;
/ e4 s+ Y1 R2 O# c  w- {　　double dx2= q->x - l_end->x; : B, }- z2 m) G1 F, v
　　double dy2= q->y - l_end->y;
' f, e2 E# R% ]& Q　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); & {! f" W% g/ _* W2 L
　　}
  M9 O/ u8 k0 v1 D6 @! l　　/* 2 line segments (s1, s2) are intersect? */
" N( M/ g5 C/ `( z3 B3 v1 A　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, 0 q- K5 a- l- _/ V/ s( o
　　const vertex_t* s2_start, const vertex_t* s2_end) 4 Z( _$ Y! f- o) L, ?. F/ S
　　{
/ z6 w; u. K; O* o$ K  M4 @# S　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && " u( T" K& o1 q( b
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
! z% @7 p# U. J4 S5 n/ j+ a　　} : ?" e8 V; B# w5 i1 S

) p& g2 `8 @. ^6 t+ l' l
% B/ t  w! c  z5 R' `  `9 L  l. g, d　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：. {4 I5 p1 g. K, C8 {) U2 }
6 W- I4 W8 d: s" t
以下是引用片段：
3 E' c- Q9 d8 @! [& ^! N  g. I　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
" Y+ [5 N( O) d* K" f8 |: R3 Z　　const vertex_t* v) , L  n* V9 a  h9 t9 L8 v
　　{
# B, u6 j8 T9 ^8 B+ Z$ _7 e　　int i, j, k1, k2, c;
4 j+ c9 V7 l. {. y+ F　　rect_t rc; ! B2 c- \  O. Q
　　vertex_t w;
! O, S* P: n' s9 b2 B' X3 u　　if (np < 3) . ~7 {4 _& Z8 O" |
　　return 0; 3 i/ o2 _& N" D8 {2 `( x! S
　　vertices_get_extent(vl, np, &rc); , e3 W( }3 E2 l7 c
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
9 z, Q( Y3 g) u3 R　　return 0;
" ?2 f7 @% T. U: v　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ ! Y! Z  \, d% h! e1 X5 C9 o5 N2 D+ h
　　w.x = rc.max_x + DBL_EPSILON;
1 D: d2 Z3 w  ?5 _7 r" D　　w.y = v->y; ( R' \$ x6 ]5 N
　　c = 0; /* Intersection points counter */
7 p  g0 w7 D1 U1 b4 R　　for(i=0; i  
9 R. ]# q% l: _' z　　{ 8 K7 T# ^( K0 g( O- Z- n& {, w5 d
　　j = (i+1) % np; 3 l0 G5 w: ?# L  O) Y$ \
　　if(is_intersect(vl+i, vl+j, v, &w)) % o+ c8 g6 P+ d) H4 o+ ]  Q
　　{ : i# q$ X9 i( T$ J
　　C++; & T9 G. |$ z0 v: _4 W
　　} - x/ m4 u& q* b* |% E
　　else if(vl.y==w.y) 1 \4 T: k1 v" Y* X( b
　　{
! J. O) D# `' f4 Y5 D　　k1 = (np+i-1)%np; $ ]* H: F, }5 c9 {# n: r- n
　　while(k1!=i && vl[k1].y==w.y) 9 l) |4 f1 b; H
　　k1 = (np+k1-1)%np;
& M/ C' V& _5 Y# [) j3 p- e　　k2 = (i+1)%np;
5 b4 e4 J* @% O  [) Q) a) j　　while(k2!=i && vl[k2].y==w.y)
, O4 ]( U. W9 ]! a  {8 I4 z  x　　k2 = (k2+1)%np;
2 n; o2 f( n0 d& a% V) n; _" @　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
3 b" |. o/ J2 u" |5 _( L　　C++;
2 W$ v7 s" g% O, z+ c7 a/ |8 L　　if(k2 <= i)
5 U! W/ S3 t8 d9 s, F　　break;
6 d1 z' d- @+ q3 e0 M　　i = k2;
/ x6 ]5 C9 ]  k* B  Q0 x　　}
; |9 M% M. R0 ]4 N4 ?' w' I　　} 8 q3 \4 C# V' `( ?
　　return c%2;
6 Z7 y  o; M. B& a　　}
8 o. _$ n* e  Q5 N; ~9 V6 `: r1 y  ?0 E" |0 ~

! I+ ^0 n- d$ ?5 G& c　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。