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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。7 }8 t( O8 A8 e' `2 \. Q2 c
: A: o4 v! O+ X8 |& O2 n
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。
. f; M: A* Y. ?7 R4 [; c3 C$ I7 z6 a8 o7 U% P
　　首先定义点结构如下：( p3 H6 L# M/ s$ P1 ^7 D8 `  {
0 l: A: b' g& t% z
以下是引用片段：
. U: d( S  B& j7 T& s! h　　/* Vertex structure */
1 C3 ^$ L/ X$ `* U( F6 @; |, A　　typedef struct
3 J) d8 s& m) P6 T2 X　　{
" t7 [8 b+ m+ d0 q　　double x, y; " V8 p" X, U- T# q: N
　　} vertex_t;
2 A6 y( g3 R: P  F0 i* k' W- O  B4 P3 c& r* Z; \4 t
: n# z- i! c0 ~% M# r( F
　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：' h* m. C5 T7 c1 e' l1 c5 |+ s

& t3 o6 A9 X8 d以下是引用片段：7 f8 e, k+ V( S# u
　　/* Vertex list structure – polygon */
4 i3 p, Z5 r  O/ X) f4 o　　typedef struct
3 U/ O% |, c7 W4 R# K0 Y　　{ 2 F1 ^( }7 N+ B: D7 V8 n1 ~
　　int num_vertices; /* Number of vertices in list */
" ]4 `( b5 T1 L0 U2 C) G　　vertex_t *vertex; /* Vertex array pointer */ 6 z# ?( z5 L/ h1 E( f
　　} vertexlist_t;
8 m4 L8 J( t. e3 e3 S$ t, e/ A, v0 K, N7 e  t3 o- g4 _

1 y. R" a7 E+ Q4 B+ X　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：7 q, g. R3 ~) \8 i+ H8 ?$ K0 _
9 M, v+ J5 E* p- H4 l
以下是引用片段：
2 [5 H5 W4 K0 I$ [) Z; W+ @5 M+ M　　/* bounding rectangle type */
; [. ]7 Q0 l6 b4 X5 z5 `　　typedef struct 2 s# w9 J# n- B3 N2 m( w- N. `+ o
　　{
; q# {+ ]* j8 {# c　　double min_x, min_y, max_x, max_y; " n: c5 j; j/ W. j9 I+ Y: E6 H  q
　　} rect_t;
# Y( @0 k# D% E  V　　/* gets extent of vertices */ ( C9 j5 S5 f0 }2 S
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ : l) I; j& F  i& s. _+ j
　　rect_t* rc /* out extent*/ ) 9 g, G+ v, T1 ~' U) R
　　{ 9 Z$ p9 |  d" h2 u3 K7 \
　　int i; 5 l% I3 [) I( D- P/ Q- ^
　　if (np > 0){
6 g8 p  l8 ?1 Q　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
% r; w* ~! r4 F8 W　　}else{ & K: w3 I+ t- L" x, v5 O
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ : ^; b7 J3 y0 g
　　}
, D; }# z) B" D( j: h, A　　for(i=1; i  
$ m9 [/ P# S1 }+ q$ r8 x* X- H0 s  [　　{ ; ~* S8 {; {& y3 @8 a
　　if(vl.x < rc->min_x) rc->min_x = vl.x; # P, H: z5 K  y" J1 A
　　if(vl.y < rc->min_y) rc->min_y = vl.y; 2 `) [6 p! P) `  F+ f
　　if(vl.x > rc->max_x) rc->max_x = vl.x;
2 C$ J% ?' o$ O　　if(vl.y > rc->max_y) rc->max_y = vl.y;
& T2 c/ L! a2 o( ?　　} 5 q$ ?- z! K: F& D
　　} ) p( e3 e! f1 n; c1 a
8 ?1 F) K6 X3 V# c9 R2 ?7 }
5 ]) R1 b( a' a; J' h  ^* _( y: v
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。0 ~' ~9 o6 {% ?* b" r) \: @! ^( i

4 ^6 o2 A) W  `) [3 v7 v# C% y: h　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：4 _8 c( J$ d; N8 j6 V( c9 o
4 C+ z% n! D* j
　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;4 O9 E5 n7 F3 x1 r" N

! N0 I2 @+ w1 U3 T$ q4 N　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
2 B9 i9 x; Y! a: J+ y, C! M& T
; N0 |! C& \" z  e, }以下是引用片段：
# Q, `6 ^& C( M/ p　　/* p, q is on the same of line l */
0 M. t' G' Q% b$ e, g4 t　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ , \# z3 b) Q' }/ ?% G) c+ G% {' v1 M
　　const vertex_t* p,
& G7 ?. F% V8 A# l; m  q　　const vertex_t* q)
6 g# Y" F& J9 I( e' e" i　　{ 8 K1 b+ Y* j( s* c# C
　　double dx = l_end->x - l_start->x;
2 E! L/ E& |: b7 {# e% K　　double dy = l_end->y - l_start->y; ) n* x# D# t+ b; m/ Q. n
　　double dx1= p->x - l_start->x;
. z9 f  h6 W- E+ e　　double dy1= p->y - l_start->y; $ m5 n/ C" s5 j. g
　　double dx2= q->x - l_end->x; % S3 I7 Z" @7 G8 d" H
　　double dy2= q->y - l_end->y;
6 C/ Q3 i1 O" M: ^8 ]9 s# K; b　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); : S% p+ Z  x# e9 v
　　} + X' u+ @* S0 I6 `2 ?. i& H7 N; L
　　/* 2 line segments (s1, s2) are intersect? */
9 l. G+ o6 p& a　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
- n# |5 L6 ~  c; E+ v" ~4 @　　const vertex_t* s2_start, const vertex_t* s2_end) ' \1 \* U1 w8 N2 v- p8 m7 n8 z* ]
　　{
& v7 K. V0 [0 U, h! i% k　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && 1 L# l( Y0 }7 \: e! ~
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
4 q: O- v: {+ W* @: I% t　　}
& D7 ]* v; r$ c: Q; X; r
. J5 P" {; R7 t, T& l
# W4 V' E, d. r　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：
% f1 ~. K, S1 ~5 |' P8 n; v
, }; Q& U2 ?2 q4 M) E9 d: @以下是引用片段：" j1 f. w# Q  m0 m
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ : I& E9 [# L- }- X8 S4 i: Z6 u! `
　　const vertex_t* v)
' W8 a6 s/ w7 O  h7 ]& k　　{
, @9 z; `# f, ^　　int i, j, k1, k2, c; - t4 \- {1 d. c- L: ]3 x
　　rect_t rc; / U1 F7 @/ g. o! \' g
　　vertex_t w; & m7 f" G& ^& D
　　if (np < 3) ) j# m  |: r2 O2 [3 M
　　return 0;
1 E: \( }* G9 T/ n　　vertices_get_extent(vl, np, &rc);
2 t0 J) V! a. _3 m# e# t　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) , |1 D  U! h- O* B
　　return 0;
* U3 t" `. D" ?2 T$ L$ u　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ 1 r; N1 {4 B7 g: Q  V
　　w.x = rc.max_x + DBL_EPSILON; ! j5 ]7 c) @" @" Z# p
　　w.y = v->y; ) W- Z6 o& Q4 d, e
　　c = 0; /* Intersection points counter */ ; K+ u/ i. l, i5 l/ n
　　for(i=0; i  
& o  h4 o% @% c+ L  _" R6 c4 B　　{ ( O( I# l' ~, ]' q0 O& R& d8 l
　　j = (i+1) % np; 4 V0 l; B! {. \5 N1 r% ~: J+ F
　　if(is_intersect(vl+i, vl+j, v, &w))
  Y6 A; K% i/ f' |$ L　　{
+ i) A- G' ?* R8 y; u- j　　C++; 1 ?, a( M$ J6 F2 `: V! m% B, ?
　　} ) S; a  `- q! z& @, i: o
　　else if(vl.y==w.y) & H3 S( o. |/ R
　　{ 6 d4 W- B7 E+ B5 L: O* U1 j
　　k1 = (np+i-1)%np; " c* y) D4 f$ O
　　while(k1!=i && vl[k1].y==w.y) - Z% E) y# \* H
　　k1 = (np+k1-1)%np;
5 n  C( b! v+ x" E+ f　　k2 = (i+1)%np; 7 k9 \4 {! h0 g
　　while(k2!=i && vl[k2].y==w.y) ' e# }# `( t! L2 t- ^# v! Z9 Y6 ^5 O
　　k2 = (k2+1)%np; & P! J2 L8 @+ F/ q3 D9 Y4 }
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) - c9 ?1 I) B2 ?1 a* n
　　C++; 6 S9 d1 [* B( @8 p7 I' M$ a# {2 h
　　if(k2 <= i)
( F+ R% ?. M% q9 _2 [$ v　　break;
2 M! S, ^% V4 L/ l8 K: u# b8 e　　i = k2; 8 f- _3 B5 \* m+ w3 m
　　}
0 L2 K# H& U3 T, {0 J4 t& r; {　　} . G1 n8 v# u# c6 s, A9 ^! d; M/ I! r
　　return c%2; ! z  ^9 B6 Z+ b! ^! \* P" I
　　} ; o* R1 I! H0 B6 K/ O% ?/ w

) d! T8 _" r+ F5 L1 l# N/ l" k7 q
' o3 _' v5 ^7 w) m　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。