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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。" B0 i' m* ]- U$ H* O
! ]' O0 t; Z  t+ O: x, X
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。& t4 {5 X0 R- }' k/ u% S! p
! G8 D8 {; k# f6 a% _) Y
　　首先定义点结构如下：( a. F; v' o- M

: S0 t8 F+ ^( g) S6 M以下是引用片段：& z* b5 G% X* Z7 V! i" d1 q9 z6 G
　　/* Vertex structure */
6 d5 E+ w# f8 [: I; m　　typedef struct 9 k. y8 w% Y6 S. }  |
　　{ 0 \: h( ~; x$ B: I* M
　　double x, y;
/ V) o3 J: L- }! _" W　　} vertex_t;
" f8 f* Q& @- [9 {6 b7 z7 g4 F2 E6 c: Z7 g

9 _- M4 h9 O( m( K0 S' D　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：
. V9 y7 z# v+ _9 _( |) k8 e. A' y3 g( h# g0 N; E
以下是引用片段：
6 N5 D! ~$ }- |7 Q　　/* Vertex list structure – polygon */ 3 V( y! a, N" ~2 e0 K9 ^
　　typedef struct
- A. T5 ]# O+ d! L3 b　　{ / v. v8 J' R" C+ S6 k' ^; F
　　int num_vertices; /* Number of vertices in list */ 1 H7 H" H- {3 d: y- G
　　vertex_t *vertex; /* Vertex array pointer */
7 g. F6 m4 v! i8 S* G4 ^. m　　} vertexlist_t;
  K/ v: R( M; O& c2 J; m8 F
9 o% i# w9 r" t1 z; E0 ?- ^- Z3 x" S6 ~  c
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：% `! |$ w8 b1 k# x0 @! ]3 B& ]& l/ ~

0 s9 d) S+ `: k8 c& O5 V3 w以下是引用片段：% \/ W- V$ e; R; R9 T4 o
　　/* bounding rectangle type */
! W1 f3 r1 a! I8 z　　typedef struct + s" c$ N+ b( d& Y2 C- v" q2 I" H  j  Y
　　{ - h( x7 ~. n$ U2 E
　　double min_x, min_y, max_x, max_y; ' L2 B: M5 h+ z# d" Z
　　} rect_t;
- V- J3 m. B& o; w, A$ j　　/* gets extent of vertices */   w3 g1 z" \8 A& v% X  A9 b: Y6 p
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ # }2 I) r" M9 q9 _- O
　　rect_t* rc /* out extent*/ ) " p% |3 f/ {" R; k% b$ ~
　　{ : K! [1 \. K1 C; A: Q6 W  F: l( g" Z
　　int i;
1 \, b, v, j% Z0 T7 D" a$ R　　if (np > 0){
, T: N0 ~" S* w　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; ( O7 B% @/ r7 u6 f3 P* G
　　}else{
6 J7 N6 C, I$ q　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
  w5 w, s* I, w( \　　} ' o0 F. c. H$ h$ J
　　for(i=1; i  
8 j* u: c9 B" f+ c0 X　　{ # {- ^) M5 p) ^7 g5 Q
　　if(vl.x < rc->min_x) rc->min_x = vl.x;
0 [8 d* \0 e0 Y. y5 B　　if(vl.y < rc->min_y) rc->min_y = vl.y;   x# W5 Y' f" y- [! p4 b
　　if(vl.x > rc->max_x) rc->max_x = vl.x; / k4 s, |/ `+ s9 q$ m. c3 a1 @: @
　　if(vl.y > rc->max_y) rc->max_y = vl.y;
4 s9 m' B. K, U' C+ w; L　　}
8 R- `! N2 w6 F6 e8 I　　}
8 o, h0 t' d; `: V' f- W+ c6 q
; n2 q* W, L7 m; C1 b6 Q/ Q# T
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。5 w+ y( a2 l4 O$ r- J

  G# g; J& t% n5 F5 T9 I: [' U: p　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：! @, j: \' j( Z& S! I2 N/ N8 Z

& O4 U* G% W, M/ k  s& u　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
+ M, Z- |1 u( n6 k: d" S
' g$ T" j7 a+ ?  U& a5 @, g5 H0 {9 |　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
4 y' a, }# x2 p3 m+ X+ P
! r3 X# n; L. h/ C7 ]* Q以下是引用片段：3 m4 M9 E5 N* n2 s; }
　　/* p, q is on the same of line l */ 7 @+ l3 T/ V9 @% p
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
* o( g+ X, a* c% s; I  c　　const vertex_t* p,
# f3 @  {/ q/ C; f" @3 W　　const vertex_t* q)
) e* A  [. [* |, A. n& u8 p　　{ ! `( d0 }7 E$ G2 |
　　double dx = l_end->x - l_start->x; 7 o' e" u2 d9 \7 R6 n5 u
　　double dy = l_end->y - l_start->y; * ~& H! l' d) O# R
　　double dx1= p->x - l_start->x; 7 G. k: S  P. g9 S( K5 ~
　　double dy1= p->y - l_start->y; ! p) }. r9 q% D% O4 n2 ~
　　double dx2= q->x - l_end->x;
# N) r; ]$ W3 p* r/ A. ~! ]　　double dy2= q->y - l_end->y;
' G8 E3 D% _0 J" T8 ?　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
; x" c, V# n* W3 d6 h9 b　　} * M0 W$ I* p0 a7 @, x' f3 q
　　/* 2 line segments (s1, s2) are intersect? */
3 h' i' K0 Y) t* a( O& @1 a　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
. V7 b/ s1 ]0 t! k* O　　const vertex_t* s2_start, const vertex_t* s2_end) 0 |+ d$ c  O' e7 |: r* {
　　{ ) d" T1 g* g9 R% _9 g- i
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
0 ]% A; p( h& t. j+ m. G7 ^& U　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; * [# Y' O# D! Y
　　}
, J( w/ K9 ?, x+ A$ [) S! M9 p8 u, a  H
1 r9 h; O. ~9 A9 T
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：- `: \, S' o6 t

7 S* U+ i3 i* K. U4 p以下是引用片段：$ x% a! h0 e) B1 V; K9 q, U+ {' y
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ $ o% x$ o# \7 }! k  G3 y
　　const vertex_t* v) 0 s9 O" D  g) u2 h
　　{
4 v# E& D$ U, d. E. a+ O　　int i, j, k1, k2, c; " p% v  R7 U. e( _& R6 h6 j
　　rect_t rc;
6 x8 i) i* O4 B% D+ p* F  u# A　　vertex_t w;
* g0 X; J* E$ l) N" i& d' |　　if (np < 3)
9 O7 u( E9 A. ?/ ^　　return 0; " k9 z& O2 ?! N% |5 h/ @
　　vertices_get_extent(vl, np, &rc); / ?' y0 e4 D* k: i8 K
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)   `- c& S7 m9 k6 n0 b6 a- ^
　　return 0; * A' ^' E: V. o9 L
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
$ ?' T' H8 m: S  x# p　　w.x = rc.max_x + DBL_EPSILON; 3 e% L8 U3 U7 _. Q5 z7 E
　　w.y = v->y;
; f$ V9 j- R* a/ P　　c = 0; /* Intersection points counter */ 0 [7 t" Y1 p+ C% W: q5 p' T
　　for(i=0; i  1 W8 e# j) y8 c0 v: c7 ]' O
　　{
8 a" S" b0 C% j8 `　　j = (i+1) % np;
; b! }8 Z- d6 ^! e/ [　　if(is_intersect(vl+i, vl+j, v, &w))
* V9 o. {! _; R( ~! c0 o　　{
3 l  n" i1 ?* t3 K( n- ?　　C++; ; q. E3 s9 i$ q6 Q& ~. R
　　} ! i$ I4 u  R, q& ^! P  m
　　else if(vl.y==w.y) 6 }5 b* ~  s' d1 ]9 E; G: G6 F
　　{
  r9 m* k3 ?& L; u　　k1 = (np+i-1)%np;
1 {, A# U% C; y/ X9 K6 d　　while(k1!=i && vl[k1].y==w.y)
1 U) k4 q8 T5 v  R+ E2 x) x: D　　k1 = (np+k1-1)%np; " v" @  k0 N% B( T) }: L7 I3 w
　　k2 = (i+1)%np; 1 h$ q- n6 L# d2 G. L. |
　　while(k2!=i && vl[k2].y==w.y)
8 e& p5 V4 L3 s! S% ^　　k2 = (k2+1)%np;
" |' R( h7 w/ H& j  N  m　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) ' c3 W/ A- E7 d, u3 Q
　　C++;
' S& ~3 O0 I- J　　if(k2 <= i)
6 X! g1 j  ?6 {  z　　break;
* D0 Y. Y4 N% K& a　　i = k2; & y! M( ~+ Z. a- I5 }6 q
　　} 2 T# m7 Y9 d0 g$ b: r% p
　　}
7 _( n+ b3 U/ X- Z　　return c%2; : D% E/ k, t  e- e
　　}
- d$ E/ }! W3 ]* |7 y' u; T( n
- ~3 `- D; f2 B
, C6 C; g! v3 n- H; G: P+ [　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。