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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
; e+ Z; p0 p/ a3 A3 i3 A3 `4 q4 J5 j7 n1 M
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
0 ]* m, d( ?0 T' A/ w; _9 T6 Z
$ u/ j- q, O  A6 N) i' I+ T　　首先定义点结构如下：3 H* x/ N! s) }* Y1 W8 G

; P( m6 M' ^- J! u2 a以下是引用片段：. d9 h( f7 z% E5 [0 g# C
　　/* Vertex structure */
; [; {8 @7 W' C; f6 f$ o) w　　typedef struct
" J- i  T2 Q7 R- ^/ z　　{ ; O7 T- J9 ^: Y4 ?% C1 E  j
　　double x, y;
& y; l) m$ @* H1 V  M# Q　　} vertex_t; 2 L& v+ v1 S4 j, q. X8 f* F

3 X2 n5 M* J" P* Q5 m
" p. s! @7 ^" C' u5 _* _: U　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：2 E$ n6 V6 j1 R6 T1 w; X9 R
0 W+ e! V6 V# g4 s
以下是引用片段：1 H2 z  z" G' @$ J6 l
　　/* Vertex list structure – polygon */ ; {8 v( Y5 \% }9 r. V$ H" y2 P% d
　　typedef struct   A/ I1 {6 T- D: P! h( n* r4 w, E7 _
　　{ 5 @6 f, U/ E# N% ~$ k, A8 I
　　int num_vertices; /* Number of vertices in list */ " p# z1 [+ R2 n
　　vertex_t *vertex; /* Vertex array pointer */
/ W" U# s( O3 v8 ^/ B! N, @　　} vertexlist_t; ; E9 l  j) ^, l/ z2 c6 j) T
/ ^: v3 @& g: Z( j; O

. L, K, }& H( }# |2 S　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：3 R* N# s1 x+ A
3 d  m% A$ N+ V# z. l
以下是引用片段：% D0 G) h% w7 F7 M
　　/* bounding rectangle type */ : }- B2 c# ]! q2 ?
　　typedef struct - y, f( I9 F' s- \- x4 C2 e7 r2 Z: i
　　{ 6 L( K, N- ]1 S" G4 i2 [
　　double min_x, min_y, max_x, max_y;
2 w, l+ u  W6 O1 V- q( D　　} rect_t; ; ^/ i  W4 j1 I2 L7 f+ x
　　/* gets extent of vertices */ ! {* m% s1 T6 Z' F
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ ; Y( f" A( x& u* [6 W
　　rect_t* rc /* out extent*/ )
! T9 w; @, |. r　　{ / H8 Q9 J" n) E5 M- @  ]
　　int i; 2 s' n' c2 z4 g8 ^' G8 `
　　if (np > 0){ ; W: |9 s& ~) m" W" p+ h
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; 3 s2 [: y; [# R1 U9 I/ o
　　}else{ % t4 w+ F4 o. V# f! ]
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
$ X3 a: I, k2 v& z3 C& {) r6 y　　}
0 \# T4 d0 q3 v& a　　for(i=1; i  1 }) K' ]0 |* J+ W2 p: C
　　{ * M# Y2 c, W" M7 U& ^
　　if(vl.x < rc->min_x) rc->min_x = vl.x; - \4 V, |, M1 T7 T/ z; y
　　if(vl.y < rc->min_y) rc->min_y = vl.y;
5 P% W) G3 X- o7 q) x6 J　　if(vl.x > rc->max_x) rc->max_x = vl.x;
. U- M2 {2 W6 g0 p4 Q5 v1 H3 I　　if(vl.y > rc->max_y) rc->max_y = vl.y; ; Y* i7 a/ t; v; d. E* U
　　} # `2 k, F9 z2 H7 M4 l1 V+ ^$ b# B/ Q
　　} 9 C9 E- Y# W9 Z8 A& n& b5 e+ g

* Q5 X  V% E0 n" O- _0 v% U" f' |7 w+ c
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。2 G- P9 X7 c1 t8 `2 i
. ?4 F. x9 W) R8 `
　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：) M; Y* Z, I8 d$ n3 }2 A% ?

' U$ ^) \" x6 h6 j0 @" n+ S2 L　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;
. u2 I; N; c2 }1 b5 I4 K# Q3 W3 H4 B5 E8 W- P
　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
% @6 D" V0 o2 S# x/ ]9 ?; W; T. X! W( ~7 d# V3 V& x6 H
以下是引用片段：
% g2 @8 H& ?5 e( u$ X. w　　/* p, q is on the same of line l */ & ~3 ^3 ?4 R$ X: G9 T( A
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ ( a% T7 C6 R( H  ]$ k' H1 k
　　const vertex_t* p, & B3 J- I+ [+ g% d' v$ i3 I: v
　　const vertex_t* q) & z3 _! B0 O0 H/ I. P) K! @
　　{ $ M8 I6 Z' y) \) j) b
　　double dx = l_end->x - l_start->x;
1 w( U4 p. p  y: ]* h" t　　double dy = l_end->y - l_start->y; & t& B5 B( L  Z; y2 R% t! C
　　double dx1= p->x - l_start->x;
' O: w5 Z1 }9 l' |7 A$ g　　double dy1= p->y - l_start->y;
* m9 o& \. L6 Q+ Q　　double dx2= q->x - l_end->x;
! c4 H0 v* i) M. g8 q/ C　　double dy2= q->y - l_end->y;
- }# S0 M+ D1 M. q' o　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
/ G; u/ |& v# k　　}
1 y/ i3 D- W! |" Y8 i# s2 ]+ F　　/* 2 line segments (s1, s2) are intersect? */
7 |* m& \0 F; g* ^5 c8 K0 y　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
( K# V. P: V* Y7 D) _　　const vertex_t* s2_start, const vertex_t* s2_end) & K( n0 Z6 X% U2 I- P0 e
　　{ - H0 H3 N, x4 {- E0 }: ^* \
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
# z6 n% f5 n+ i; |( [# _( F/ P　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
- p3 H# m* r/ N1 f( a　　} - Y2 C2 q8 {2 L  j6 n1 D' r( g$ y
. V! R( o" r5 @! g
, h6 m  ]. v2 U
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
/ R0 ~" p. b! @/ Q* r% I" ?- `8 l0 P2 A: i2 _
以下是引用片段：1 z( ^1 Z/ \, q" d  r' N
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ & N* x& c8 d& v$ ~" s3 O
　　const vertex_t* v) 2 Y5 }0 W7 w: u) y# x7 L) Q$ J
　　{ $ r2 r" \3 x, Q. Q
　　int i, j, k1, k2, c; 8 k( X" B9 q4 i4 d8 l1 S& Y
　　rect_t rc; 2 Q" t; ?( g4 w- q; s" l
　　vertex_t w; , o2 V, [/ n' W9 x* o- @& w1 d8 F
　　if (np < 3) & K5 h8 V) f, E) ]! l
　　return 0;
, \/ k7 I* A) R7 ~　　vertices_get_extent(vl, np, &rc);
4 X7 m8 q7 ^4 {　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) 4 P" W3 g' ?( \+ G: O8 A9 I
　　return 0; 3 x. [% y1 z, K3 `8 C3 y
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ ( r$ r6 H4 k5 |8 n" k( K
　　w.x = rc.max_x + DBL_EPSILON;
, s. ~7 P; d4 S, n/ R　　w.y = v->y;
+ j! ^# e% j" ^( F+ y　　c = 0; /* Intersection points counter */ # v$ F( d5 U# i, m
　　for(i=0; i  
3 S$ ^) l* E5 B  t' W' j% i$ @* {　　{
! s# I6 M8 X9 V- j$ P' V　　j = (i+1) % np;   Q* j! [& p% r6 A, o
　　if(is_intersect(vl+i, vl+j, v, &w))
8 y( H! Y+ G1 B) O' W　　{
3 d% J+ o4 w2 d  m　　C++;
  |2 Y: q/ l3 |1 X9 Y- ]　　} ( Q8 y+ N! F9 v$ a7 J, o
　　else if(vl.y==w.y) 3 N, e* Y  i) B3 O: r# ]! R7 _  F
　　{ + o5 X- v! Z1 ?5 ~
　　k1 = (np+i-1)%np;
/ z. n: C" P) d# E1 `; ^　　while(k1!=i && vl[k1].y==w.y)
) H9 Q$ E: u% G2 z% S$ F2 Y　　k1 = (np+k1-1)%np; 9 L. v( q+ a: C9 t" T
　　k2 = (i+1)%np; ' |' B8 ?* K( j$ [4 }4 S$ [
　　while(k2!=i && vl[k2].y==w.y) 7 W9 `( Q7 T1 L7 o  S
　　k2 = (k2+1)%np; , V! t; Y3 ?4 N3 r2 W6 z
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)   ~: w* s% n# `  |7 q# q: ]3 S
　　C++;
9 w: V3 _; j# c　　if(k2 <= i)
9 N, O, E: O- S4 ]　　break;
, j6 L& p2 y0 W* c  m& E　　i = k2;
: D0 L5 I3 Z. J2 \/ [! Z4 x　　}
3 X$ @4 i( K  S% p  N# s; Q　　} ( g3 \( a* ~$ `, W9 U, ?
　　return c%2; # s7 |" f7 Z, P7 D
　　}
- d* G* C! b  d. E( a( Z& q7 T: \8 c
& T2 M* ]  {! w& R6 T4 V' r' s
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。