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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。7 B; f5 X: A+ }& u) r0 {" `
3 r4 c1 l" @+ i7 ?+ k) C6 @
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。/ C6 |* S) a: ~0 E3 Z
+ `! X" S1 h5 O* B+ G
　　首先定义点结构如下：) ~- M3 `: F9 r* v
4 Q. @' i( Y8 N) ?& E% M* V& K" X
以下是引用片段：
# |1 @/ L  ~) j　　/* Vertex structure */ 4 t6 D! @( ]# v5 g  B
　　typedef struct # A& I1 F3 z( ^* h) I! h
　　{ 6 P) \- [  I1 a* _
　　double x, y; 8 h0 r0 T) Q  x9 x% h- E7 e" }* F
　　} vertex_t;
1 X0 L. X- V7 G! X3 s
3 l) B, e1 z* x) ^6 h3 W
' s/ C; q7 ?+ q0 V4 h) ~　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：1 |3 m- p7 C* i: N5 b
: b, B* b/ V! b( ]8 L9 t8 C. F
以下是引用片段：
# U' s2 K3 G' N; H　　/* Vertex list structure – polygon */
; x) W' C. b+ S0 Y　　typedef struct : k4 W5 s3 |3 b+ H8 h6 t$ O
　　{ 4 G% G* ]) G. J8 G/ n8 S/ `2 P
　　int num_vertices; /* Number of vertices in list */
& b" P4 r8 ?0 N5 ?- ^　　vertex_t *vertex; /* Vertex array pointer */
, x3 T4 e, F0 Q: T) n　　} vertexlist_t; . f/ [( t$ Z& l: s  t1 g4 A8 @9 a
9 M! e! P1 n! i3 Y' A* I/ X
/ Y, X* q9 Q$ p' ?& s6 U
　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：
) u4 U9 `9 D$ b+ }( q. {
0 z6 O# m1 w1 `( Y以下是引用片段：
+ c0 s6 o' E0 o　　/* bounding rectangle type */
# E3 J& H: m* a0 g　　typedef struct 3 ~. V# j  [/ {+ q4 m, L. O
　　{ / a2 B# F" f! a  \; ]
　　double min_x, min_y, max_x, max_y;
6 u( ~) ^4 c* W, E. I+ O　　} rect_t; 6 Y; }$ m( g6 S9 [
　　/* gets extent of vertices */
+ e, ^- k8 [" u% Q　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ ' h% \* z& ]2 @$ S: `+ ]* _- m3 @
　　rect_t* rc /* out extent*/ ) " m2 Z3 O3 s0 B- Y- d3 E
　　{
# t' @9 C7 M0 i) j* ]8 x* i5 W; _　　int i;
3 F- C7 k, o) E; N# G  v) F- a　　if (np > 0){ 2 h: R2 P& `. h. b# [& L
　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; * w9 @- n! X% G& B! b, s3 M! ?
　　}else{
( i- Z4 v. `6 P: m6 Y　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ 5 w" m9 W( C; G7 I" i. N
　　} " h& ?5 ?' k) i9 S$ r4 ~" J
　　for(i=1; i  7 C2 _8 @, T5 }4 w* _- a+ \: B4 B
　　{
5 Z; S# @) l/ h& ]　　if(vl.x < rc->min_x) rc->min_x = vl.x; , b4 f& C2 Z  n' r) c& Z: v; \9 i
　　if(vl.y < rc->min_y) rc->min_y = vl.y; 1 ^: C5 D& e- M4 n) \6 A1 Z
　　if(vl.x > rc->max_x) rc->max_x = vl.x;
3 R# z% \7 I) I( p　　if(vl.y > rc->max_y) rc->max_y = vl.y;
4 ^5 w+ w; d. ^; w& _$ p" r　　}
- T2 f- A. ?0 l) A$ j　　}
9 d7 ]" G: S0 e# l3 V
2 L! U' L9 {( h! d; G) m0 I4 i+ j' B) \4 h& A& z2 s
　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。8 U' f6 P0 Y7 ]- e- q- U

5 X# P4 `3 O) s/ D　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：
- B7 @, l2 ^+ P" B7 Z8 A9 p
2 ^3 W8 S' T( N. x) c; t3 R　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;2 ^$ |! j: `2 J: q6 z

9 J( q' e8 ?1 K$ @. W6 ~2 o& ]  }5 m　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
+ _" S9 X# P0 _+ f* L% P4 k0 ?9 e
* c- c/ z5 h2 g以下是引用片段：
' S8 U8 W- _  l3 R$ k　　/* p, q is on the same of line l */ . L& K- M8 J* d# [% D4 w1 o
　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */ % y  @# q% ]7 W" R
　　const vertex_t* p, - g+ t) {/ O9 J7 t
　　const vertex_t* q) ' T* }% j; y- j% j. w. j% J
　　{ 1 y0 D+ c2 F6 d- F' y
　　double dx = l_end->x - l_start->x;
! m7 H* @' U5 e$ t　　double dy = l_end->y - l_start->y;
* @3 C( @4 o% S8 R　　double dx1= p->x - l_start->x;   n5 q6 n1 |! D) q
　　double dy1= p->y - l_start->y;
2 l: [+ V) _' r" ^( m+ A" a  `　　double dx2= q->x - l_end->x; $ N+ E: h* N8 B. K# I
　　double dy2= q->y - l_end->y; 8 A5 V# A; b4 O
　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0); : f% _4 y6 _4 o" M4 Y
　　} * N6 X  G9 w4 Q4 p4 r
　　/* 2 line segments (s1, s2) are intersect? */
. C9 w% v1 W# Q9 F' }4 M　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
6 {# n" h& \1 s7 C　　const vertex_t* s2_start, const vertex_t* s2_end)
+ _, I0 `) T- J- L  q　　{ * ^# Q. G# N7 |/ x$ B. r* |  |1 H8 L
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
: j7 W. k0 n4 |% u% B. A　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
4 q% h/ Y# T1 ]9 @- B7 y　　}
) a& b# k5 y# h# a! P
& R5 [, V: K. g  J5 F
- r7 ^8 i- U  I# \  ?1 W　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：$ F, O$ o0 O% S" \) j: j
: {) J$ e8 q0 R" n$ g! _* w/ P1 E
以下是引用片段：' S+ s& `, t  J% Q. Q5 {( g
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ ; U$ e# ~% x) J1 k& {
　　const vertex_t* v)
: o2 D: b: ?+ [1 s; L3 J! w1 s　　{ ; t0 i( C% u. i: D( u9 p7 C
　　int i, j, k1, k2, c;
' a/ i/ m5 K# ]　　rect_t rc;
0 G  I' l' o; P8 y% n) t; T) _　　vertex_t w;
, B# V8 E0 A7 D3 q6 Q# e  J  F　　if (np < 3) " E$ _3 q# T9 C  e
　　return 0;
: ~2 |* p: o1 P, _: M# x　　vertices_get_extent(vl, np, &rc); - g! j% y0 i+ T  M
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y) 7 }- {3 A( |! H- b
　　return 0; # z# S! C9 l+ x; F
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */ 3 O. G! v( s( q/ ^2 _, ^5 s
　　w.x = rc.max_x + DBL_EPSILON;
, S; {6 `# {7 m2 A) u　　w.y = v->y; " o. k6 D( |6 Q6 ~+ ~5 m
　　c = 0; /* Intersection points counter */ 8 A; \3 `, p2 P' M8 D) c
　　for(i=0; i  / h$ a+ W$ X1 r- B* Y- @- o
　　{
2 R2 w6 |7 R/ s2 O% l+ H) N　　j = (i+1) % np;
) t: D6 m3 l& d% p' }$ u　　if(is_intersect(vl+i, vl+j, v, &w)) 6 J- {0 s. B+ |6 M6 l  _5 g
　　{ 1 b( d+ R- k4 h4 B# o. K6 M
　　C++;
1 c7 l2 p, t7 q! ^) _8 j" T8 H　　}
8 j7 o! W3 z) Y! o7 a　　else if(vl.y==w.y)
( R' m  v3 I! S1 l2 q1 f　　{
- ?5 a' x3 z7 A4 x3 j! a　　k1 = (np+i-1)%np; 0 Y$ x2 G: y7 R$ Y8 w: E/ \- `
　　while(k1!=i && vl[k1].y==w.y)
; y' i! K: u0 a: w3 W" f　　k1 = (np+k1-1)%np;
8 s0 {1 @: G# o+ f, e5 `# Q　　k2 = (i+1)%np;
' t5 L4 l' P' s/ ?$ ]% D1 s6 G　　while(k2!=i && vl[k2].y==w.y)
8 y8 V# M+ u; d, L! v2 y% f0 `' E　　k2 = (k2+1)%np; : ]! s- b( A" i0 n% r- \1 V. K( H- X
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0)
  T) x, B" q" @9 W　　C++;
. W; D. ?+ ^& a7 D; L& {; w( R　　if(k2 <= i) ) |9 T5 H6 o9 E, p& x3 H
　　break; % V- L+ ]: N7 v" p6 B$ p
　　i = k2; 6 k  R) A# \0 `9 w" |% q
　　} 3 y& `0 N2 z6 i( O* y& t+ T  k9 N
　　} : G! Z0 x- ^1 e4 i; J
　　return c%2;
5 r! c7 {8 i" x6 c/ S' N$ I* v　　}
; [) S( i( H  W; l) W
! j! h  \# V: |' ~' r! c: @5 z4 j5 L' y% l
　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。