捌玖网络工作室's Archiver

zw2004 发表于 2008-1-21 17:20

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

本文是采用射线法判断点是否在多边形内的C语言程序。多年前,我自己实现了这样一个算法。但是随着时间的推移,我决定重写这个代码。参考周培德的《计算几何》一书,结合我的实践和经验,我相信,在这个算法的实现上,这是你迄今为止遇到的最优的代码。k7UT+[v+Be
'b+\5|F.@/rm'pD:n)|K
  这是个C语言的小算法的实现程序,本来不想放到这里。可是,当我自己要实现这样一个算法的时候,想在网上找个现成的,考察下来竟然一个符合需要的也没有。我对自己大学读书时写的代码没有信心,所以,决定重新写一个,并把它放到这里,以飨读者。也增加一下BLOG的点击量。
,v+hr&kY:t9p
#fFTmR^x   首先定义点结构如下:
V~*\ S umST#}:C
w Z|V,T;G8Y ?)k,W4q 以下是引用片段:%m%R X p8L/_\
  /* Vertex structure */ .b!cp*B3|o"f1qg
  typedef struct
y2t\k,e   {
Y}IZnCI"]   double x, y;
Tj7w]_IL'} Ax   } vertex_t;
3~q:?"d o?5W a )Z"n6X#f.FA6D-c

?I6fOE6ck d   本算法里所指的多边形,是指由一系列点序列组成的封闭简单多边形。它的首尾点可以是或不是同一个点(不强制要求首尾点是同一个点)。这样的多边形可以是任意形状的,包括多条边在一条绝对直线上。因此,定义多边形结构如下: Q0B"?#Z R

e m\vb4rx 以下是引用片段:;f+E6|&e ir+c't
  /* Vertex list structure – polygon */ T7e7NPbr
  typedef struct
2q$['d5R-P:v(e i   {
"}&_ ZC-Y}   int num_vertices; /* Number of vertices in list */ B3x:e.X"Vy6i/wAC e~
  vertex_t *vertex; /* Vertex array pointer */ 0ZOhUI/bvs
  } vertexlist_t;
;] kL)g-B ra6^
aIE"[` AX@
.o'~V5T-C   为加快判别速度,首先计算多边形的外包矩形(rect_t),判断点是否落在外包矩形内,只有满足落在外包矩形内的条件的点,才进入下一步的计算。为此,引入外包矩形结构rect_t和求点集合的外包矩形内的方法vertices_get_extent,代码如下:
a [qH${$\a Q%Nl'ed3N4n:yX
以下是引用片段: Sew9jZ${f
  /* bounding rectangle type */
&Qf+?#]?(ko;S   typedef struct
:ksWO*{   { B+R#N:R j6h0ad
  double min_x, min_y, max_x, max_y;
"P];]0px;V6a ?/QNc gE_   } rect_t; ERe(q+t[1l2pA(r
  /* gets extent of vertices */ IB8ng0c,w)bf7h
  void vertices_get_extent (const vertex_t* vl, int np, /* in vertices */ D6Yg%a2])eM H(k
  rect_t* rc /* out extent*/ ) 3U4V c^K6x C\
  {
Z^~V]6K   int i; T*mzS5I
  if (np > 0){ 9Tk`:E\0V#|
  rc->min_x = rc->max_x = vl[0].x; rc->min_y = rc->max_y = vl[0].y;
.t7jEr4Y P#mW   }else{
| pX |q^f.n$A#c&o   rc->min_x = rc->min_y = rc->max_x = rc->max_y = 0; /* =0 ? no vertices at all */ WE&c9b#?&YS Mn!y
  } k0XD}S4?b.W#K
  for(i=1; i  M$eyp/JOY `
  {
Jp:B J'r[q;QHj   if(vl[i].x < rc->min_x) rc->min_x = vl[i].x;
hdt jV%V   if(vl[i].y < rc->min_y) rc->min_y = vl[i].y; f!gsVj\"Ccp
  if(vl[i].x > rc->max_x) rc->max_x = vl[i].x;
-?ES8~"RV   if(vl[i].y > rc->max_y) rc->max_y = vl[i].y; 6Q*@F~k
  } :jy q+U~jq
  } gDd7K7xnb,z x
#r B4Tr#?Hrk

4Q/r0a5ti-` ja`   当点满足落在多边形外包矩形内的条件,要进一步判断点(v)是否在多边形(vl:np)内。本程序采用射线法,由待测试点(v)水平引出一条射线B(v,w),计算B与vl边线的交点数目,记为c,根据奇内偶外原则(c为奇数说明v在vl内,否则v不在vl内)判断点是否在多边形内。d#c]9[|0m&N

)H%^XD fRimA   具体原理就不多说。为计算线段间是否存在交点,引入下面的函数:
"i5bU7r| u"f"}FN
y4?t1@`I   (1)is_same判断2(p、q)个点是(1)否(0)在直线l(l_start,l_end)的同侧;
R$WHc%h|G%s
.UBa$]S8| W\   (2)is_intersect用来判断2条线段(不是直线)s1、s2是(1)否(0)相交;
P7y0~&ge tUe+cB0R pi:U?D,D6Y3a
以下是引用片段:
U0X,c["Hh0Jt   /* p, q is on the same of line l */ (c:eX$L!R
  static int is_same(const vertex_t* l_start, const vertex_t* l_end, /* line l */
3f _;t I'x(@+l_R   const vertex_t* p, P_2Vrw#P*n
  const vertex_t* q) ;by)p1m8?f
  { qI%O l/GlN
  double dx = l_end->x - l_start->x; 'E"qXG_
  double dy = l_end->y - l_start->y; 9pxD1T[$cN\
  double dx1= p->x - l_start->x;
_{ ?9Xr3j }G   double dy1= p->y - l_start->y;
\+S/M8iU1Q)o   double dx2= q->x - l_end->x; %YN? S+]"g:oL
  double dy2= q->y - l_end->y;
}e#\(RQw1x   return ((dx*dy1-dy*dx1)*(dx*dy2-dy*dx2) > 0? 1 : 0);
b[+b/x?wE)M   } ^ `k)S9LJy
  /* 2 line segments (s1, s2) are intersect? */
OU q_uNO a9V   static int is_intersect(const vertex_t* s1_start, const vertex_t* s1_end,
:^!?R.iH(X$yH-@y   const vertex_t* s2_start, const vertex_t* s2_end)
U.E6El+zv   {
8R X#s+s{L/l   return (is_same(s1_start, s1_end, s2_start, s2_end)==0 &&
o_rxnJ#e   is_same(s2_start, s2_end, s1_start, s1_end)==0)? 1: 0;
j[$s[WU nz   } W&ki/\`r6B8x
XoY[ } [L\&w
P:ef&l&n)C v&p
  下面的函数pt_in_poly就是判断点(v)是(1)否(0)在多边形(vl:np)内的程序:m!F Go1i B}+BT l a7s

\j'z4L4Po'HK&A(d 以下是引用片段:
cy5bR:^ [   int pt_in_poly ( const vertex_t* vl, int np, /* polygon vl with np vertices */ 1rcA'YD
  const vertex_t* v)
*c)F"jR xk!v;E{2`vwK   { J#X%r"ly)Hh
  int i, j, k1, k2, c; b&}&{~ D mKl
  rect_t rc;
6Wk5S.D'J.k?   vertex_t w; /n l6F'U9p
  if (np < 3)
8]Q*r t|   return 0;
VOI2z/n{+g   vertices_get_extent(vl, np, &rc);
z+iw1| M m&Yk U;?   if (v->x < rc.min_x || v->x > rc.max_x || v->y < rc.min_y || v->y > rc.max_y)
o7{@c JCi{   return 0;
4d$nj'b j+b&ag   /* Set a horizontal beam l(*v, w) from v to the ultra right */
j%^,~2V:_ X   w.x = rc.max_x + DBL_EPSILON;
(T$G`"K2pU @G&F   w.y = v->y;
1av^*| G   c = 0; /* Intersection points counter */ (}f^ xJ Ux
  for(i=0; i  Gk(["B m+E;j(V
  { ,M(F p/_c$C
  j = (i+1) % np; 5KQrU3t#K+A+r
  if(is_intersect(vl+i, vl+j, v, &w))
)P/?Y[x G a   {
m8PQiU{E   C++; )] Z3S7F#o,f$zT
  } {w&p9} z S
  else if(vl[i].y==w.y)
Rh-qF5Z"D   {
"C7Ie7T)c   k1 = (np+i-1)%np;
m`O!{i ~   while(k1!=i && vl[k1].y==w.y)
8t;Fz^ i   k1 = (np+k1-1)%np;
g4IE)_.X,E-^ jq6NbQ   k2 = (i+1)%np;
(d-X1U?@/cr#r   while(k2!=i && vl[k2].y==w.y) -U/M"bBlY u
  k2 = (k2+1)%np;
)nRY-|o N   if(k1 != k2 && is_same(v, &w, vl+k1, vl+k2)==0) 6|$y;s%YAL G8T
  C++; ?/Q'Y+B"[att
  if(k2 <= i) /Y6K&q%of4Z
  break; z9A7te8nCwa
  i = k2;
.^~T8P)B^_?   } \^m#JzXrtZ{
  } Ju3l/sF_
  return c%2; d$c'QWoX
  } OZ`2i*} C w Y,s

$rt Hid&GE4I
$CaPm:n)o H8K2y   本想配些插图说明问题,但是,CSDN的文章里放图片我还没用过。以后再试吧!实践证明,本程序算法的适应性极强。但是,对于点正好落在多边形边上的极端情形,有可能得出2种不同的结果。

页: [1]
【捌玖网络】已经运行:


Powered by Discuz! Archiver 7.2  © 2001-2009 Comsenz Inc.