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

zw2004

本文是采用射线法判断点是否在多边形内的C语言程序。多年前，我自己实现了这样一个算法。但是随着时间的推移，我决定重写这个代码。参考周培德的《计算几何》一书，结合我的实践和经验，我相信，在这个算法的实现上，这是你迄今为止遇到的最优的代码。
- p9 j# s0 S7 H6 i  ?/ `! l5 T9 O: \. Z5 {9 Q
　　这是个C语言的小算法的实现程序，本来不想放到这里。可是，当我自己要实现这样一个算法的时候，想在网上找个现成的，考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心，所以，决定重新写一个，并把它放到这里，以飨读者。也增加一下BLOG的点击量。# \& ~: Z' i9 A' W5 X

8 V+ S; r/ d+ C6 F6 B　　首先定义点结构如下：
, D8 U# z. z: e
7 P% V+ G9 @9 {- F以下是引用片段：
" b1 j, Y& a. N% r$ O! }4 z　　/* Vertex structure */
! A3 C9 g, b7 `, s  S5 J' u' w　　typedef struct
$ t5 q  N# P  H　　{
' X6 y5 y/ o  _* ~; F　　double x, y;
) }+ B; x% J! ~# _0 `. l& O　　} vertex_t;
, @5 G  M7 y% K" |
% R$ t  P: a# v5 k7 @, m% r& ]" m# f" _
4 e" ]5 ]; c$ A) P. ~　　本算法里所指的多边形，是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的，包括多条边在一条绝对直线上。因此，定义多边形结构如下：5 K) T# s( R. B' T" z- N9 h+ r

; B0 I8 g* w; {以下是引用片段：
' n6 |  R  R7 j: e0 x* V　　/* Vertex list structure – polygon */
+ F. D; W3 C/ W3 O2 \" B4 t　　typedef struct % ?9 r' ^4 F0 i* g8 I5 ~  `
　　{
( T2 Y- e8 p! n+ l* l　　int num_vertices; /* Number of vertices in list */ # S3 v9 K2 ]% {
　　vertex_t *vertex; /* Vertex array pointer */
; }# j+ k- Y% }- n3 |( ~' H　　} vertexlist_t; + d' D  o  \- j' H' z1 f- ?
& e8 f1 [6 p" z/ ?$ [. W

# m4 D% o2 H: W+ y　　为加快判别速度，首先计算多边形的外包矩形(rect_t)，判断点是否落在外包矩形内，只有满足落在外包矩形内的条件的点，才进入下一步的计算。为此，引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent，代码如下：8 n3 x4 [2 Z* p7 c  |: T
* r, f7 K4 A  S& I+ P
以下是引用片段：
+ \2 I' ?9 ?9 X0 S3 n1 u　　/* bounding rectangle type */
/ n9 l( K2 ]% }. b6 `　　typedef struct
& k( [+ a: n  `$ }$ Q' l& [! m　　{ 2 q/ l3 J/ Q! l8 z5 s
　　double min_x, min_y, max_x, max_y;
8 v2 M% j  ?' s2 `6 T/ v　　} rect_t; : P  C8 {9 h" |
　　/* gets extent of vertices */ ' M4 t' K5 j' H
　　void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ 1 F& d" D$ K7 G6 z
　　rect_t* rc /* out extent*/ ) " |  P2 G$ x9 w5 `. |7 A
　　{
, w1 q  u* v! u　　int i; 8 g% {" ~# D; H0 \$ t) z" X4 T' |
　　if (np > 0){
# B4 W: l( A/ l4 M) y9 ], Z　　rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y; 7 q9 l" Q; Q9 p2 }0 O* O
　　}else{ . F3 M# T6 C( z5 G) K* y8 l1 S# _
　　rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */
, a; l. [* f- W. G( k0 B　　} 8 |! E; g, x# g& d2 B5 A1 h- Z
　　for(i=1; i  ) m* L; ?% Z( T% ^7 z* E! I; I: ~
　　{ 5 U4 m8 Y* d  a! A" `& b; Y1 i
　　if(vl.x < rc->min_x) rc->min_x = vl.x; , E0 r. u" K* L- D" L, ~  t
　　if(vl.y < rc->min_y) rc->min_y = vl.y; : }6 ]0 a4 R7 \' k0 _
　　if(vl.x > rc->max_x) rc->max_x = vl.x; ) Z$ w4 Q, u& ]! o& b
　　if(vl.y > rc->max_y) rc->max_y = vl.y; 1 k/ w! Z6 u2 ~+ ?; K7 j1 `5 F0 v
　　} # B- y) k1 A5 I! [2 u
　　}
4 N! Z) x# l4 a( H; Q3 L% q8 A$ F& W+ l' d( h+ I  b* J

8 S, w! e* m' w　　当点满足落在多边形外包矩形内的条件，要进一步判断点(v)是否在多边形(vl：np)内。本程序采用射线法，由待测试点(v)水平引出一条射线B(v，w)，计算B与vl边线的交点数目，记为c，根据奇内偶外原则(c为奇数说明v在vl内，否则v不在vl内)判断点是否在多边形内。
$ h1 d1 L: x( L9 q6 t  ]- d$ N
+ H3 J! V  p1 I5 A5 ^3 p$ D2 W　　具体原理就不多说。为计算线段间是否存在交点，引入下面的函数：5 e# |) @( |+ u# @; `9 @5 F& B

7 R3 g% S2 D% Y8 Z6 P, J　　(1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start，l_end)的同侧;& G6 k: @) |' @; W% \. `5 Z

  e( z+ h5 B+ ~# b　　(2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;) K) m6 p8 w3 H4 N$ d

  _; w. L% b1 P以下是引用片段：5 \* g$ I* i" w, z3 C! L
　　/* p, q is on the same of line l */
, i7 O: D$ Z) e( a( Q* v　　static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
, J! g' p8 x) B8 Q' g　　const vertex_t* p, . k1 F6 {" S  D# V! }
　　const vertex_t* q)
" a' L6 w* {, r$ f; f　　{
5 e/ J; B# S' Z; r6 l( i$ L　　double dx = l_end->x - l_start->x; " O+ v: c7 u8 P) ?5 \
　　double dy = l_end->y - l_start->y; 0 ^+ e6 G- w2 G% s& {
　　double dx1= p->x - l_start->x; ' h) V7 T. ~' e  I, e, e
　　double dy1= p->y - l_start->y; 1 M* b$ G6 W1 v/ w, N
　　double dx2= q->x - l_end->x;
. H7 i. K' b( A* \) n　　double dy2= q->y - l_end->y;
, i- f  ~. R( e# S8 Y  l6 y9 u9 ~+ }4 C　　return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
5 j( e- _+ ^' a- ]3 k* ~! k　　}
) b( w  r6 v5 ^4 V" L5 H2 {　　/* 2 line segments (s1, s2) are intersect? */ & |, ]9 V9 g6 f6 G' c3 s4 D. W* w9 m
　　static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
! u& O) A3 }; R" M' Z/ H　　const vertex_t* s2_start, const vertex_t* s2_end) ) u% n3 x! O1 U" }( K
　　{ ! c  h* A% @7 z7 v) ?
　　return (is_same(s1_start, s1_end, s2_start, s2_end)==0 && 1 G/ r* ~4 t) {0 F- s- }; Y: p- @
　　is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0; 1 e. \8 D# Z8 T- ~$ Q
　　} % K. G$ F6 @) g. p( B
5 d* }2 D1 {9 d6 t- n/ X  J
# r+ z$ |1 a5 N, q% v
　　下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl：np)内的程序：! M1 M5 I4 W1 G. ~9 N

$ F4 W7 d$ \3 v- R以下是引用片段：6 e/ F' ~7 q( ~; C7 R; P
　　int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */
( Q! A) B1 c0 n0 g　　const vertex_t* v)
2 J6 E* y# f/ s; Y1 \　　{
! N3 V1 d7 }# A8 f" G　　int i, j, k1, k2, c;
2 d; E3 V' [5 o) U& W5 d% B0 w　　rect_t rc;
) d, b+ z; e9 i+ C　　vertex_t w; , E' U0 B9 I$ ]- q- \
　　if (np < 3) $ n% n% a2 D& R- b
　　return 0;
# J% X9 z  M: v6 H& S; k  i2 q! P　　vertices_get_extent(vl, np, &rc); / m% y% a* u9 h
　　if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
$ _, r& g* B2 ~1 U+ `9 l3 l　　return 0; - s& V+ I4 x# J9 u3 \0 G2 R/ B
　　/* Set a horizontal beam l(*v, w) from v to the ultra right */
4 h/ n2 {0 i- @+ O0 i　　w.x = rc.max_x + DBL_EPSILON; . o2 ]8 ]& {) c# r4 ]% S1 f
　　w.y = v->y;
  G: r( C$ u; f& A　　c = 0; /* Intersection points counter */ ! c- E; @- M% {0 D2 v5 f
　　for(i=0; i  ! s$ x5 K( t  R# E& `+ e5 V$ }- E
　　{ 0 _% A, S. A; A! y/ k
　　j = (i+1) % np; 8 c" H2 |" [6 q, W' h  T/ m
　　if(is_intersect(vl+i, vl+j, v, &w)) * Z, r7 r0 Y" @" G. d- E# v4 V) R
　　{
/ j% l% C: Z) U0 P# o　　C++; : F7 m7 r. Q! g9 L& }
　　}
3 @7 z4 ^; S# V9 s3 I/ @9 z7 _　　else if(vl.y==w.y)
9 b4 H( U0 `/ u$ r9 O; g. C　　{ : S1 O' X& K+ R( ^
　　k1 = (np+i-1)%np;
$ n: _, v( A, ]  ^　　while(k1!=i && vl[k1].y==w.y)
$ {5 X, e7 I0 ]' L, u( i　　k1 = (np+k1-1)%np; 3 t% s/ H0 s; |$ @- x
　　k2 = (i+1)%np;
+ }  H% k  \0 Z4 }7 x( O+ _　　while(k2!=i && vl[k2].y==w.y) 3 \7 R5 q& s4 x/ `  ?% U
　　k2 = (k2+1)%np; ( Q4 J2 z& o( `. n5 p9 d
　　if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) - j5 V! g6 Y2 H, c
　　C++; ( X, E' |- W! B4 m6 E* w- B
　　if(k2 <= i) ! D/ r' A  L7 e. s5 V5 P; m9 i
　　break;
5 k' g" F4 P# n3 n5 `- e　　i = k2; 3 X+ i3 z2 |& d" ^2 u' G
　　}
1 k- q# z. l4 P　　} " L" w$ `6 v2 A) |4 e" M
　　return c%2; ( ?; e2 F8 T$ k) j1 h, ~
　　} 3 w' F$ X' ?  q3 m9 ]

3 H0 v, D4 {5 u, \
' K2 b& @: s7 ~3 I/ c　　本想配些插图说明问题，但是，CSDN的文章里放图片我还没用过。以后再试吧!实践证明，本程序算法的适应性极强。但是，对于点正好落在多边形边上的极端情形，有可能得出2种不同的结果。