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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
0 T$ a# ^" s; ?7 r; W- b1 l: T. X- E+ k0 `# t+ r2 A, I
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
; Z3 ~! T; R5 H# x; B+ \5 M+ p; r5 H
　　首先定义点结构如下：
5 S( X$ N: x3 }9 l& m% r
$ l& @' c+ _: Y* ~+ _$ T以下是引用片段：$ X( A# _" r1 t& F/ g
　　/* Vertex structure */ # H/ G4 b! f" W
　　typedef struct
- f3 _3 x+ {  H3 A- M$ l　　{ + Q  c$ P0 t6 w  [
　　double x, y;
+ P8 [+ E3 N5 f. i　　} vertex_t; % P1 m6 c1 w5 Q
  p5 ]: A' e! \: [& ?/ q0 [

' y6 j5 Q1 n, P& H　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：( J. ^, S  V$ z3 z
* Q- T8 `* }7 e% B! `1 ]
以下是引用片段：
4 U6 t1 r% t0 M3 h+ C　　/* Vertex list structure – polygon */ % C" L) Z; t+ Z4 a
　　typedef struct 4 C3 X) y  b( @- R; B; Q
　　{
$ p1 W- u5 g# E" B) M4 J　　int num_vertices; /* Number of vertices in list */ ; r2 o9 a) k& _& b
　　vertex_t *vertex; /* Vertex array pointer */ & s0 ?  M$ ]% a# @2 `  {/ a' W- S
　　} vertexlist_t;
+ T0 n0 t' r; n1 Z, @% v# i$ D
, T% q! Y4 r0 F. _7 T- K, P2 c5 H' l3 n
' t; {7 z' Y8 d6 U) O2 u. l　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：+ F+ D0 H! |4 ^0 r# w; g
# T+ {* Z$ y1 z/ a$ a- A" g) O
以下是引用片段：3 v) X( g' `6 O4 t1 u
　　/* bounding rectangle type */ , H8 p1 q# {% I
　　typedef struct 4 B( {: k8 m4 J) o+ M' V! x+ `
　　{
# j: y( G2 @/ Z- g  `　　double min_x, min_y, max_x, max_y; ) y; w8 k0 T1 S" S
　　} rect_t;
# g9 ]6 y, n  ^1 l; V6 A8 f/ c　　/* gets extent of vertices */ 9 i: y4 w  B. K$ c- ~2 N
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */
0 [* V# i/ H0 m* ]( z+ Z) j3 V$ ~　　rect_t* rc /* out extent*/ ) " m2 f' I! f  F( n7 I4 ?8 r1 i0 r
　　{ " t) S5 W+ d+ Q  }2 ^
　　int i;
2 U! \5 {# }7 ?& e& L3 }　　if (np > 0){ 5 x! m" Q! \  ]# A
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
* I1 d4 a& L5 E$ a! I$ X+ l　　}else{
; T& \. N1 h" g0 Q  a9 V; C. h　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
4 {, k! o% [3 {( j, P# \　　} , T4 u) C+ s2 b2 e6 `3 h/ T8 E5 ^
　　for(i=1; i  - B3 G- A5 I' K) M. i+ Q1 \
　　{ ( m4 @7 m1 V* @
　　if(vl.x < rc->min_x) rc->min_x = vl.x; : ]: F& O5 Z) y
　　if(vl.y < rc->min_y) rc->min_y = vl.y;
: A5 e, D4 X1 M6 |6 s$ _7 R　　if(vl.x > rc->max_x) rc->max_x = vl.x;
+ k* L/ P; C# X+ P, o　　if(vl.y > rc->max_y) rc->max_y = vl.y; + {$ {# ]% ^, W, l
　　}
  q8 c( n9 L/ ?) i# d# I4 D& K　　} 4 I+ P& L+ s# A' n
) r/ T" J" @5 [4 L) i

7 M, ?, ]+ @3 n+ Z) @　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。( \, p6 y7 p5 }* z' r5 W8 B
1 S+ m5 f9 A! b: C$ n
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
. G9 g; H; K" P! i4 M/ Q" _
) D7 J7 ~* k7 i8 ^% u　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;0 _- Z$ P) Z$ D# ~: C" p: s9 n& F- U
' K2 ]+ U( W3 {+ w3 n+ v9 P8 w7 x: E/ Y
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;$ ^  b. n2 d" C

' Q; P- X; E* u3 C5 @- r& i1 R以下是引用片段：' C5 y. O( `: M
　　/* p, q is on the same of line l */ & H- |" |0 z- W% p8 _2 `9 H
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ # o% h* |0 y# a( ?
　　const vertex_t* p,
, o0 q# I% O" O　　const vertex_t* q) : f8 k) M9 ?" @! b6 |
　　{ 4 ?3 ?! f4 t" L0 h1 O- U# e/ G
　　double dx = l_end->x - l_start->x;
& ?& U, [8 \8 Z2 [$ g4 z2 g　　double dy = l_end->y - l_start->y; * c( P) d6 c0 U0 q
　　double dx1= p->x - l_start->x; + w  P7 U$ a$ ~" b6 Y- U
　　double dy1= p->y - l_start->y;
- m* x  F. u+ Z$ ~& c　　double dx2= q->x - l_end->x; + s2 x6 J3 @- D+ h2 Q; l0 b9 \
　　double dy2= q->y - l_end->y;   Y* t) e/ ~8 X
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); 0 ]1 T* f: W$ p  L, I9 p" S, x: F4 N2 P
　　}
+ a7 j  I+ _; l5 `# W( H: m　　/* 2 line segments (s1, s2) are intersect? */
" q# {4 Q7 s+ {. N5 j/ R9 H　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
5 [+ T# l: t/ [, H, k( k( F　　const vertex_t* s2_start, const vertex_t* s2_end) * @. [4 }2 }: s! ~* k6 u% h
　　{
* q4 Q  ?: V* ]0 \, C' A% N% v6 j　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && ' j! l1 k5 {8 R
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; 0 m  Z, V4 ~$ W: g/ {. T3 v3 v
　　}
7 L4 b! J. @- B% \. a" a$ I
9 P9 P& Q) z4 o" u' n* V1 _. W
8 |9 i6 l% K$ T) h, r1 p3 i& \" K7 n/ v　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
5 N& |8 m( Y, D, @8 T' Z: ?! G% `2 c9 {5 [2 o
以下是引用片段：: d& T* C* L6 G$ l7 n% c- n
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ ' |/ s3 M$ [3 [+ c8 R. E: ?5 c
　　const vertex_t* v)
7 P8 ^$ W1 M& A3 Q  I& n) Z　　{
3 D1 b5 s' V' b6 j3 h　　int i, j, k1, k2, c; : _5 n* I  ^- p# w/ W5 e
　　rect_t rc;
8 A7 v) o: P0 S+ u+ t! x$ U. h　　vertex_t w; 0 t" R1 ?, s+ T6 N+ d+ P
　　if (np < 3) * t9 d7 J. g" H5 B9 C$ S* d5 Q
　　return 0;
6 e% R% k' ?+ R3 e  j3 F: w  ?7 s0 P　　vertices_get_extent(vl, np, &rc); % y& G7 i3 ]- f
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) % }0 x0 d1 f9 O, d9 V% z
　　return 0;
6 }/ K7 i5 X2 ^5 u0 Z# D8 v4 e　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
3 e+ Z  G. ]3 A, V. ?6 h/ {: F　　w.x = rc.max_x + DBL_EPSILON; " Q, G' c4 u) G0 }1 F
　　w.y = v->y; % u, m# J0 n' J* T
　　c = 0; /* Intersection points counter */ - J9 r9 w) v5 k; P0 u
　　for(i=0; i  
/ D- q2 ]+ \& x; t& ^0 U　　{
7 ]- I, _- D, O& i6 _% h　　j = (i+1) % np; 9 m: a' R  Y* W$ a6 |+ S4 Q8 H
　　if(is_intersect(vl+i, vl+j, v, &w))
( l3 E1 @( L8 Y　　{
2 C0 y+ q5 i$ m, a2 u　　C++; 8 ]; Q, [& }* m# M% T: g
　　} ! h4 V: o9 ]+ ~8 |8 O+ P. n& g4 m
　　else if(vl.y==w.y) 4 b% }- L% q! s' o3 F
　　{
' g" y) Y! t  y6 s& ^　　k1 = (np+i-1)%np;
3 z8 k/ Q. d# I$ I* Y　　while(k1!=i && vl[k1].y==w.y)
5 s) s! b6 R* p) B9 _+ a$ R. i　　k1 = (np+k1-1)%np;
- [' @/ n+ O) Y' s% `4 y. M　　k2 = (i+1)%np;   V, S$ s. g8 X6 L1 [/ B& b
　　while(k2!=i && vl[k2].y==w.y) 5 p( ~& a8 e( D5 U; I! u1 i
　　k2 = (k2+1)%np;
# Q; u. R1 y6 A6 ^4 V　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) - @! o5 t* R, B
　　C++;
6 L3 d! Z: a3 [　　if(k2 <= i) 3 j2 ?( z' {! ~0 a+ I
　　break;
; C1 Q  x. M) ?4 }4 P　　i = k2; . s( P4 R) ?5 L! _: W
　　} 6 s+ [8 Z6 x$ f
　　}
, _( m4 B& N- o! o( a　　return c%2; : Q  R$ N; [- l$ B- _9 r
　　} - K/ ~! _' P+ J0 G
' ^  u+ G2 z' \: ]4 D
/ m  U5 B& C8 e$ q5 c
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。