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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。3 y- K/ L) ~$ M7 A4 ]# j

) `& a9 h. [- }6 K& W　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。) }- B/ b- M  ]- }4 G/ Z
7 F* ?+ T: ?; ]' ]4 R! n' w
　　首先定义点结构如下：
5 l* U4 A9 S& n, Y- N2 a( w" w+ N$ @, v9 @4 u+ ]  i" R) J/ k
以下是引用片段：
7 @& t: S/ s( L" e9 F- V　　/* Vertex structure */
$ a1 h* W$ g2 p& m# Z　　typedef struct # Q& `- S7 y% v4 F! m' O) q
　　{ & Z7 y6 q& Q) e
　　double x, y;   q' x" N. C5 V4 t' A' e  W$ ]
　　} vertex_t;
7 {& n! D; I, ^3 Q. E: s7 S8 i/ T: `6 D( M* j' }
2 O' }/ e# D4 r6 y* y9 \+ I
　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
4 B4 U- o" ]. I! A, m, B
' N. B( S& R3 P" L" B以下是引用片段：# n8 g" r) z' D- ^; L2 S  l. `$ i
　　/* Vertex list structure – polygon */ & k3 R! ?6 t1 o9 O0 _/ Z; `3 _
　　typedef struct 0 i4 o5 T% ~- N& _5 v: v0 y
　　{
( q. q! T& `" w: J! R: P　　int num_vertices; /* Number of vertices in list */ ) n$ _7 N2 F2 V/ k
　　vertex_t *vertex; /* Vertex array pointer */
# x7 J0 C( {% V0 k/ ^5 ~6 y　　} vertexlist_t;
: ?8 A7 z/ _4 j- q* h9 B% v0 `/ u3 X" V, I! c2 L
- W: k/ _# C* {+ J, l0 j( W
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：: G# j% p4 Q; R; c4 f/ ~* o* q
) C, f4 e' B1 Q) Y: H# U! v( \2 n
以下是引用片段：3 @3 W8 I! t; e" d7 v; c# H
　　/* bounding rectangle type */
: Q1 M& s7 M2 d4 l　　typedef struct
  I2 J8 [$ u. I  `& q( _　　{ 1 @. u8 \+ Q8 W! z' `. L: \
　　double min_x, min_y, max_x, max_y; ; p' N1 ^5 ]6 w8 h
　　} rect_t;
8 i2 P7 ?1 y2 Q& G9 k6 F9 A- X# h　　/* gets extent of vertices */ 1 r& j  d: ^4 I$ N9 }9 }) D
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ 5 \! U2 g" O6 i% X9 \0 q9 K! L6 S( s1 g
　　rect_t* rc /* out extent*/ ) 5 M* B+ G. g. r! Q
　　{ # R4 C, \, |4 C( g1 @
　　int i; " N$ |& y5 Q/ Y( S
　　if (np > 0){ # m  m% P7 ^! f1 j: }
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; 2 Q& }- Z# M, j, A; @+ A- L1 [- ?5 c0 S
　　}else{
$ _1 A4 V& c2 t3 i/ d' p: b　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */   I- ?: D1 B/ W# h& Q: P, F6 Z* M* g
　　}
2 M. _! I. e' K1 ^& y0 ~5 G9 ^! t　　for(i=1; i  
+ D8 d; q" B& m" r  x0 k% q* P* p　　{   I( H+ J, Z) ^6 q5 l. o0 E% s
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
# C4 H! Q3 G# M6 N9 S　　if(vl.y < rc->min_y) rc->min_y = vl.y;
  l! }) G  Q* V3 K' r+ S% B) m　　if(vl.x > rc->max_x) rc->max_x = vl.x; 3 {. b8 M% \. ^- y% U
　　if(vl.y > rc->max_y) rc->max_y = vl.y; 0 S* x( f0 ?$ j9 m) ]6 |( H  }' V
　　}
( z0 z9 J6 \  U5 F3 `　　} 0 k& j1 X+ j. i1 F/ h; s0 T* ~. o
5 M, K( V! T+ P1 Y' ^
1 }' q) R: M/ n4 E! i
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
" `& O- d- |6 n' E% v4 j2 D+ q4 L) _( K( a
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
( E4 B/ A8 W9 ^& |2 B8 X" D' p# V- J2 R/ ^
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;  \; M$ b7 M3 }) @  m" \
5 [6 E# @, N2 @7 u/ ?
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;0 ^3 {1 T5 D8 ?2 \$ ^5 e+ }% X
4 N: r4 m8 i8 h; B1 k7 h; n
以下是引用片段：
) R; b' d0 ?) i( J3 R　　/* p, q is on the same of line l */ . G9 L( `$ O4 V1 h, i
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ # R( o' I. U: r/ i6 U" E8 ]
　　const vertex_t* p,
4 a& j. e; a$ h- R' W# d2 P' j　　const vertex_t* q) 4 l! B% `" X/ |0 G* a+ Z
　　{ & g+ \6 A2 [6 s; m# X, i. C: ?
　　double dx = l_end->x - l_start->x;
( s: H3 N( B6 p+ H! K" d4 `+ d　　double dy = l_end->y - l_start->y; : B! S' Z/ M8 q  X
　　double dx1= p->x - l_start->x; ) _' N& M" i; w% O
　　double dy1= p->y - l_start->y;
1 E0 `+ C! m4 ^' h$ c" a　　double dx2= q->x - l_end->x; ( K' _8 y% c2 F. n
　　double dy2= q->y - l_end->y;
( I/ E! _( J3 l# ]" H　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
( D7 ?, P6 H. h! z( L　　} # i$ a: L$ C! g  k
　　/* 2 line segments (s1, s2) are intersect? */
6 \9 R/ g: D5 r2 R! ]3 C　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end, - a* k7 G1 ]: [7 ?/ w
　　const vertex_t* s2_start, const vertex_t* s2_end) $ {+ e1 o% Q8 R# x. `
　　{ $ p4 e5 U. E; [' ~) L  g- J  T
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && % Z) ~" b$ P8 y% U+ H
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
* O6 m) T% o. e9 \* ?; R% ]7 i　　} 8 y; A( M8 ~- w' ~. K" r1 P

- @9 Y6 s3 R, U( R2 O% N
, ]  ~' Q: G8 S; Y  B3 s" P　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：: [, b: D; M* e4 J+ H
, [' @9 }2 d% m. _
以下是引用片段：
; U1 \- F6 h9 p1 h1 d9 z3 ^& j　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ 8 E0 t2 z! N4 W' o$ j9 m
　　const vertex_t* v) * X  B  n- H; y4 A$ q3 |3 I4 m$ |
　　{ " @1 p* @/ }4 G3 B: {: u
　　int i, j, k1, k2, c; 2 [7 X7 G# ~4 @+ p  H  F' l4 w
　　rect_t rc;
/ k, Q/ Q- H( C/ q　　vertex_t w;
. h. k2 x$ \8 ^' m8 w2 m　　if (np < 3) 4 S5 G" U* ?" O+ U
　　return 0;
/ C9 N  c: v" b- j　　vertices_get_extent(vl, np, &rc);
) j1 K( t+ M# d( L5 e5 I1 ]　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) ( A1 u' a- [* f' F# L/ U0 S( \2 h
　　return 0;
1 t4 @7 h2 v5 u6 f　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
  g3 S: S; \$ S$ K4 a  ^" M- J　　w.x = rc.max_x + DBL_EPSILON; ; |, D7 K& I) x/ a) K0 Y
　　w.y = v->y;
/ F1 ?5 V/ M; s/ x+ G　　c = 0; /* Intersection points counter */
% [0 t+ [3 Y9 c3 w　　for(i=0; i  ) \! {$ Z5 Y5 s$ r& K: w
　　{
- ~2 @0 f; W0 A* P0 m; `4 n0 O$ n　　j = (i+1) % np; 4 o) G3 Y/ i& V1 H3 b1 G
　　if(is_intersect(vl+i, vl+j, v, &w))
2 G3 `7 m+ I# D( s- t5 _, m　　{
# `1 i* p* M' @1 A　　C++;
, Q0 s& D. w& l) Z- l3 a　　} : X/ @; e3 r. {3 j0 k3 P0 b9 A
　　else if(vl.y==w.y) * l" _% B! a1 b& D* i- M" u
　　{ 4 K( n- E4 _1 g2 q. c7 R" R7 I5 n
　　k1 = (np+i-1)%np;
( t4 a* U( a1 W2 k3 S- o　　while(k1!=i && vl[k1].y==w.y)
; |( ], k' J; w, b+ v9 x, ]　　k1 = (np+k1-1)%np; ( p( g9 d; s, M2 _8 H$ ]3 Q3 w
　　k2 = (i+1)%np;
$ f( M4 A: H3 o) ~4 y5 Y　　while(k2!=i && vl[k2].y==w.y)
$ k, H$ Z" E  p: e3 Y! S5 a, x　　k2 = (k2+1)%np;
2 U, J) R$ \' e4 p8 g( Z　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
! Z. ?( Q9 S/ G7 t# w& R* E  E　　C++; 3 f/ x' h' N2 F
　　if(k2 <= i)
, I: r8 x3 Y9 i7 S+ D2 c- m& y7 j　　break; 4 ]. J" f1 P3 I( Y- _# D
　　i = k2; 2 {& b& S& F) \; @8 w1 Z: @# t! ^
　　}
9 p$ A6 x) A3 m  o' D　　} , E- F' B9 w* m- q1 Q1 B. H
　　return c%2; . K+ N" V- b6 t
　　} $ B5 ~5 ]) ?; O
* o# H' A% ]9 G/ [2 H7 w  Z5 J( `
2 t) J+ C) A7 l+ C& b
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。