Profil de WXYZWXYZ的异人馆PhotosBlogListes Outils Aide

Blog


30 septembre

计算几何:点是否在平面多边形内

    又是周赛,今天状态还不错。
    今天遇到一题判断两个多边形是否有重叠部分。计算几何,这部分没学过,做题过程中临时上网找资料来看了。具体题目副本有PKU的3082题。
    判断两个多边形(假设为A和B)是否有重叠部分,可以先求A是否有某个顶点在B内,或B是否有某个顶点在A内,满足其中之一,就可说AB有重叠了。这样问题就转化为一个点是否在多边形内的问题。
    上网找找,判断点是否在多边形内有三个步骤:(转自csdn)
 
    第一步:判断这个点是不是就是多边形的端点;  
    第二步:判断这个点是不是落在多边形的边界上;  
    第三步:通过这个点横向作一平行射线,判断与多边形的交点数,如果交点是顶点,则交点数加一,结果如果是奇数,则该点落在多边形之内,如果是偶数,则反之。
 
    具体算法涉及向量叉积,具体这部分不详细说了,上网轻易查到,下面贴过主算法函数吧,参考很好用。
    还是来自于csdn的。
 
 
 
  const   double   INFINITY     =   1e10;    
  const   double   ESP   =   1e-5;    
  const   int   MAX_N           =   1000;    
   
   
  struct   Point   {    
        double   x,   y;    
  };    
   
  struct   LineSegment   {    
        Point   pt1,   pt2;    
  };
 
  inline   double   max(double   x,   double   y)    
  {    
        return   (x   >   y   ?   x   :   y);    
  }    
   
  inline   double   min(double   x,   double   y)    
  {    
        return   (x   <   y   ?   x   :   y);    
  }    
   
  //   计算叉乘   |P1P0|   ×   |P2P0|    
  double   Multiply(Point   p1,   Point   p2,   Point   p0)    
  {    
  return   (   (p1.x   -   p0.x)   *   (p2.y   -   p0.y)   -   (p2.x   -   p0.x)   *   (p1.y   -   p0.y)   );    
  }    
   
  //   判断线段是否包含点point    
  bool   IsOnline(Point   point,   LineSegment   line)    
  {    
  return(   (   fabs(Multiply(line.pt1,   line.pt2,   point))   <   ESP   )   &&    
  (   (   point.x   -   line.pt1.x   )   *   (   point.x   -   line.pt2.x   )   <=   0   )   &&    
  (   (   point.y   -   line.pt1.y   )   *   (   point.y   -   line.pt2.y   )   <=   0   )   );    
  }    
   
  //   判断线段相交    
  bool   Intersect(LineSegment   L1,   LineSegment   L2)    
  {    
  return(   (max(L1.pt1.x,   L1.pt2.x)   >=   min(L2.pt1.x,   L2.pt2.x))   &&    
                          (max(L2.pt1.x,   L2.pt2.x)   >=   min(L1.pt1.x,   L1.pt2.x))   &&    
                          (max(L1.pt1.y,   L1.pt2.y)   >=   min(L2.pt1.y,   L2.pt2.y))   &&    
                          (max(L2.pt1.y,   L2.pt2.y)   >=   min(L1.pt1.y,   L1.pt2.y))   &&    
                          (Multiply(L2.pt1,   L1.pt2,   L1.pt1)   *   Multiply(L1.pt2,   L2.pt2,   L1.pt1)   >=   0)   &&    
                          (Multiply(L1.pt1,   L2.pt2,   L2.pt1)   *   Multiply(L2.pt2,   L1.pt2,   L2.pt1)   >=   0)    
                      );    
  }    
   
   
  //   判断点在多边形内    
  bool   InPolygon(Point   polygon[],   int   n,   Point   point)    
  {    
          if   (n   ==   1)   {    
                return   (   (fabs(polygon[0].x   -   point.x)   <   ESP)   &&   (fabs(polygon[0].y   -   point.y)   <   ESP)   );    
  }   else   if   (n   ==   2)   {    
  LineSegment   side;    
  side.pt1   =   polygon[0];    
  side.pt2   =   polygon[1];    
  return   IsOnline(point,   side);    
  }    
   
  int   count   =   0;    
  LineSegment   line;    
  line.pt1   =   point;    
  line.pt2.y   =   point.y;    
  line.pt2.x   =   -   INFINITY;    
   
  for(   int   i   =   0;   i   <   n;   i++   )   {    
  //   得到多边形的一条边    
  LineSegment   side;    
  side.pt1   =   polygon[i];    
  side.pt2   =   polygon[(i   +   1)   %   n];    
   
  if(   IsOnline(point,   side)   )   {    
  return   true;    
  }    
   
  //   如果side平行x轴则不作考虑    
  if(   fabs(side.pt1.y   -   side.pt2.y)   <   ESP   )   {    
  continue;    
  }    
   
  if(   IsOnline(side.pt1,   line)   )   {    
  if(   side.pt1.y   >   side.pt2.y   )   count++;    
                  }   else   if(   IsOnline(side.pt2,   line)   )   {    
                          if(   side.pt2.y   >   side.pt1.y   )   count++;    
                  }   else   if(   Intersect(line,   side)   )   {    
                          count++;    
                  }    
          }    
   
          return   (   count   %   2   ==   1   );    
  }    

    最后再贴下某人在csdn上他的计算几何库目录,自己参考一下,有时间自己写个。
 
  ㈠   点的基本运算  
  1.   平面上两点之间距离 1  
  2.   判断两点是否重合     1  
  3.   矢量叉乘 1  
  4.   矢量点乘 2  
  5.   判断点是否在线段上 2  
  6.   求一点饶某点旋转后的坐标 2  
  7.   求矢量夹角 2  
   
  ㈡   线段及直线的基本运算  
  1.   点与线段的关系 3  
  2.   求点到线段所在直线垂线的垂足 4  
  3.   点到线段的最近点 4  
  4.   点到线段所在直线的距离 4  
  5.   点到折线集的最近距离 4  
  6.   判断圆是否在多边形内 5  
  7.   求矢量夹角余弦 5  
  8.   求线段之间的夹角 5  
  9.   判断线段是否相交 6  
  10.判断线段是否相交但不交在端点处 6  
  11.求线段所在直线的方程 6  
  12.求直线的斜率 7  
  13.求直线的倾斜角 7  
  14.求点关于某直线的对称点 7  
  15.判断两条直线是否相交及求直线交点 7  
  16.判断线段是否相交,如果相交返回交点 7  
   
  ㈢   多边形常用算法模块  
  1.   判断多边形是否简单多边形 8  
  2.   检查多边形顶点的凸凹性 9  
  3.   判断多边形是否凸多边形 9  
  4.   求多边形面积 9  
  5.   判断多边形顶点的排列方向,方法一 10  
  6.   判断多边形顶点的排列方向,方法二 10  
  7.   射线法判断点是否在多边形内 10  
  8.   判断点是否在凸多边形内 11  
  9.   寻找点集的graham算法 12  
  10.寻找点集凸包的卷包裹法 13  
  11.判断线段是否在多边形内 14  
  12.求简单多边形的重心 15  
  13.求凸多边形的重心 17  
  14.求肯定在给定多边形内的一个点 17  
  15.求从多边形外一点出发到该多边形的切线 18  
  16.判断多边形的核是否存在 19  
   
  ㈣   圆的基本运算  
  1   .点是否在圆内 20  
  2   .求不共线的三点所确定的圆 21  
   
  ㈤   矩形的基本运算  
  1.已知矩形三点坐标,求第4点坐标 22  
   
  ㈥   常用算法的描述 22  
   
  ㈦   补充  
  1.两圆关系: 24  
  2.判断圆是否在矩形内: 24  
  3.点到平面的距离: 25  
  4.点是否在直线同侧: 25  
  5.镜面反射线: 25  
  6.矩形包含: 26  
  7.两圆交点: 27  
  8.两圆公共面积: 28  
  9.   圆和直线关系: 29  
  10.   内切圆: 30  
  11.   求切点: 31  
  12.   线段的左右旋: 31  
  13.公式: 32  
16 septembre

每周一赛

    刚刚的晚上sicily每周一赛,十题做了八题。总的来说题目简单,两小时不到就已经做了八题,但竟然有某题,用到二分图最佳权匹配。死反,没好好学习啊,两周前不会到现在还是不会,只了努力靠后面三小时来狂攻它了。我才拿起我那本图论书,再在网上找一找,就研究起来了。可惜三小时过去,还是不行。结束了,还是八题,15名。唉,可惜啊。还有题没做出来的,也已经想象到了,有时间刷起它就行了。
    改天研究好这个二分图最佳权匹配再写个解题报告上来。
11 septembre

公选课

    昨天上了我有生以来第一堂公选课,叫做中国人移居海外历史与现状。
    一开始我报这门课程是因为我是台山人,想到移居海外,我想我会比较感兴趣,并且应该会经常提到台山,为了了解台山我就报这个了。
    第一堂课,理论我就不管它了,考试能过就OK了,更感兴趣的事实本身。听说我们是分地区说的,先是介绍东南亚。东南亚其实也是最多华人移居的地方,说是都掌握当地的民生经济命脉,这个也是真的。但局势动乱,排华事件多也是真的。
    这些还是不说多了,说说我们上课时放的资料片吧。首先我记得的是在越南西贡的邓老先生,访问他时他是用粤语回答的,如果没有字幕的话,在场众多不是广东的学生肯定就茫然了。50年代去到越南,辛苦劳动20年,终于在当地开起了个小有名气的金铺,但是,随后却是一场暴乱,所有人被赶出城市,而当他回来后,20年的心血已经被人占领的,并且根据战后法例,他是不能要回原来的房子。而在这场动乱中,他的亲人也全部去世了,采访他时,他只有一只小狗陪伴他,在乡下的别人弃置的破旧磨坊中度过晚年,看着他拿着他的全家福,说:“以前来的时候一大群,现在只剩我一个孤零零。如果他们都在的话,一定也已经成材了。心酸啊。
    然后是还有些在缅甸的,以前是国民党的军队,二战时到缅甸参战一起对抗日本,二战后,士兵们只想平静地生活,而不愿继续在中国参加内战,因此就在当地定居下来。当地政府却要让他们对抗反政府军,才能换来在当地的居住权,他们只好战斗。虽然后来是胜利了,但士兵牺牲却很多。其中有个老兵说:“在丛林里,人吃的东西不多,吃人的东西倒不少。”全班笑倒,但其中话语的悲惨,却让人不能保持笑容。又说到一个士兵当时被车从山坡上掉下来压断了腿。而山路却不好走,“一个人那山已经够爬了,何况还抬个人!”无奈得只能按照指挥员说的给他一颗子弹。那士兵也非常清楚他不能拖累其他人,也一直在叫同伴们扔下他。采访到这里,老兵也控制不住自己的泪水了。当时埋下那个同伴时,全部在场的人都哭了,指挥员说:“你们不要以为他很悲哀,明天我们死了,谁来埋葬我们?”看到这里,气氛是非常沉重的,坐在我前面那个女同学,好像还有点泪水,拿出纸来抹。课就在这高潮中结束了。
    下课了,下一课呢?
8 septembre

激情

    今天晚上豪爽了一下,上去玩CS了。
    由于久没碰,生疏了,打得很差,都不敢用WXYZ的名字上战场了。后来和我同宿舍的,以前和我一起拿到班CS冠军的人也上来打了,这我才发挥出来。虽然大不如前,还勉强看得过去。最后我们又改回我们以前一起用的名字了,74系列的,我的是74LS10。打得淋漓尽致,有他在后面护着,可以放心杀上去。那是一种信任。上学期这种激情回来了。
    近期做题多了,有点麻木,决定参加完星期日的POJ Monthly就暂时封笔了。需要到外面的世界寻找下激情。虽然说平静的生活才是美好的,但怎么说这也是大学,不想平静,大家都有躁动的心,都有无处发泄的热情。来个超新星爆发也好吧。
    现在心情如同放假一般,除了形式上上个课以外,行动上也是和放假一样。这不是我想要的,大学就要有大学的生活。没有激情就不是年轻人了。所以我要把老龄化的意识改变过来。 
4 septembre

[转]KMP算法详解

KMP算法详解

作者:matrix67 日期:2006-11-29

    如果机房马上要关门了,或者你急着要和MM约会,请直接跳到第六个自然段。

    我们这里说的KMP不是拿来放电影的(虽然我很喜欢这个软件),而是一种算法。KMP算法是拿来处理字符串匹配的。换句话说,给你两个字符串,你需要回答,B串是否是A串的子串(A串是否包含B串)。比如,字符串A="I'm matrix67",字符串B="matrix",我们就说B是A的子串。你可以委婉地问你的MM:“假如你要向你喜欢的人表白的话,我的名字是你的告白语中的子串吗?”
    解决这类问题,通常我们的方法是枚举从A串的什么位置起开始与B匹配,然后验证是否匹配。假如A串长度为n,B串长度为m,那么这种方法的复杂度是O (mn)的。虽然很多时候复杂度达不到mn(验证时只看头一两个字母就发现不匹配了),但我们有许多“最坏情况”,比如,A= "aaaaaaaaaaaaaaaaaaaaaaaaaab",B="aaaaaaaab"。我们将介绍的是一种最坏情况下O(n)的算法(这里假设 m<=n),即传说中的KMP算法。
    之所以叫做KMP,是因为这个算法是由Knuth、Morris、Pratt三个提出来的,取了这三个人的名字的头一个字母。这时,或许你突然明白了AVL 树为什么叫AVL,或者Bellman-Ford为什么中间是一杠不是一个点。有时一个东西有七八个人研究过,那怎么命名呢?通常这个东西干脆就不用人名字命名了,免得发生争议,比如“3x+1问题”。扯远了。
    个人认为KMP是最没有必要讲的东西,因为这个东西网上能找到很多资料。但网上的讲法基本上都涉及到“移动(shift)”、“Next函数”等概念,这非常容易产生误解(至少一年半前我看这些资料学习KMP时就没搞清楚)。在这里,我换一种方法来解释KMP算法。

    假如,A="abababaababacb",B="ababacb",我们来看看KMP是怎么工作的。我们用两个指针i和j分别表示,A[i-j+ 1..i]与B[1..j]完全相等。也就是说,i是不断增加的,随着i的增加j相应地变化,且j满足以A[i]结尾的长度为j的字符串正好匹配B串的前 j个字符(j当然越大越好),现在需要检验A[i+1]和B[j+1]的关系。当A[i+1]=B[j+1]时,i和j各加一;什么时候j=m了,我们就说B是A的子串(B串已经整完了),并且可以根据这时的i值算出匹配的位置。当A[i+1]<>B[j+1],KMP的策略是调整j的位置(减小j值)使得A[i-j+1..i]与B[1..j]保持匹配且新的B[j+1]恰好与A[i+1]匹配(从而使得i和j能继续增加)。我们看一看当 i=j=5时的情况。

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B = a b a b a c b
    j = 1 2 3 4 5 6 7


    此时,A[6]<>B[6]。这表明,此时j不能等于5了,我们要把j改成比它小的值j'。j'可能是多少呢?仔细想一下,我们发现,j'必须要使得B[1..j]中的头j'个字母和末j'个字母完全相等(这样j变成了j'后才能继续保持i和j的性质)。这个j'当然要越大越好。在这里,B [1..5]="ababa",头3个字母和末3个字母都是"aba"。而当新的j为3时,A[6]恰好和B[4]相等。于是,i变成了6,而j则变成了 4:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =     a b a b a c b
    j =     1 2 3 4 5 6 7


    从上面的这个例子,我们可以看到,新的j可以取多少与i无关,只与B串有关。我们完全可以预处理出这样一个数组P[j],表示当匹配到B数组的第j个字母而第j+1个字母不能匹配了时,新的j最大是多少。P[j]应该是所有满足B[1..P[j]]=B[j-P[j]+1..j]的最大值。
    再后来,A[7]=B[5],i和j又各增加1。这时,又出现了A[i+1]<>B[j+1]的情况:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =     a b a b a c b
    j =     1 2 3 4 5 6 7


    由于P[5]=3,因此新的j=3:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =         a b a b a c b
    j =         1 2 3 4 5 6 7


    这时,新的j=3仍然不能满足A[i+1]=B[j+1],此时我们再次减小j值,将j再次更新为P[3]:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =             a b a b a c b
    j =             1 2 3 4 5 6 7


    现在,i还是7,j已经变成1了。而此时A[8]居然仍然不等于B[j+1]。这样,j必须减小到P[1],即0:

    i = 1 2 3 4 5 6 7 8 9 ……
    A = a b a b a b a a b a b …
    B =               a b a b a c b
    j =             0 1 2 3 4 5 6 7


    终于,A[8]=B[1],i变为8,j为1。事实上,有可能j到了0仍然不能满足A[i+1]=B[j+1](比如A[8]="d"时)。因此,准确的说法是,当j=0了时,我们增加i值但忽略j直到出现A[i]=B[1]为止。
    这个过程的代码很短(真的很短),我们在这里给出:

程序代码 程序代码
j:=0;for i:=1 to n dobegin   while (j>0) and (B[j+1]<>A[i]) do j:=P[j];   if B[j+1]=A[i] then j:=j+1;   if j=m then   begin      writeln('Pattern occurs with shift ',i-m);      j:=P[j];   end;end;


    最后的j:=P[j]是为了让程序继续做下去,因为我们有可能找到多处匹配。
    这个程序或许比想像中的要简单,因为对于i值的不断增加,代码用的是for循环。因此,这个代码可以这样形象地理解:扫描字符串A,并更新可以匹配到B的什么位置。

    现在,我们还遗留了两个重要的问题:一,为什么这个程序是线性的;二,如何快速预处理P数组。
    为什么这个程序是O(n)的?其实,主要的争议在于,while循环使得执行次数出现了不确定因素。我们将用到时间复杂度的摊还分析中的主要策略,简单地说就是通过观察某一个变量或函数值的变化来对零散的、杂乱的、不规则的执行次数进行累计。KMP的时间复杂度分析可谓摊还分析的典型。我们从上述程序的j 值入手。每一次执行while循环都会使j减小(但不能减成负的),而另外的改变j值的地方只有第五行。每次执行了这一行,j都只能加1;因此,整个过程中j最多加了n个1。于是,j最多只有n次减小的机会(j值减小的次数当然不能超过n,因为j永远是非负整数)。这告诉我们,while循环总共最多执行了n次。按照摊还分析的说法,平摊到每次for循环中后,一次for循环的复杂度为O(1)。整个过程显然是O(n)的。这样的分析对于后面P数组预处理的过程同样有效,同样可以得到预处理过程的复杂度为O(m)。
    预处理不需要按照P的定义写成O(m^2)甚至O(m^3)的。我们可以通过P[1],P[2],...,P[j-1]的值来获得P[j]的值。对于刚才的B="ababacb",假如我们已经求出了P[1],P[2],P[3]和P[4],看看我们应该怎么求出P[5]和P[6]。P[4]=2,那么P [5]显然等于P[4]+1,因为由P[4]可以知道,B[1,2]已经和B[3,4]相等了,现在又有B[3]=B[5],所以P[5]可以由P[4] 后面加一个字符得到。P[6]也等于P[5]+1吗?显然不是,因为B[ P[5]+1 ]<>B[6]。那么,我们要考虑“退一步”了。我们考虑P[6]是否有可能由P[5]的情况所包含的子串得到,即是否P[6]=P[ P[5] ]+1。这里想不通的话可以仔细看一下:

        1 2 3 4 5 6 7
    B = a b a b a c b
    P = 0 0 1 2 3 ?


    P[5]=3是因为B[1..3]和B[3..5]都是"aba";而P[3]=1则告诉我们,B[1]和B[5]都是"a"。既然P[6]不能由P [5]得到,或许可以由P[3]得到(如果B[2]恰好和B[6]相等的话,P[6]就等于P[3]+1了)。显然,P[6]也不能通过P[3]得到,因为B[2]<>B[6]。事实上,这样一直推到P[1]也不行,最后,我们得到,P[6]=0。
    怎么这个预处理过程跟前面的KMP主程序这么像呢?其实,KMP的预处理本身就是一个B串“自我匹配”的过程。它的代码和上面的代码神似:

程序代码 程序代码
P[1]:=0;j:=0;for i:=2 to m dobegin   while (j>0) and (B[j+1]<>B[i]) do j:=P[j];   if B[j+1]=B[i] then j:=j+1;   P[i]:=j;end;


    最后补充一点:由于KMP算法只预处理B串,因此这种算法很适合这样的问题:给定一个B串和一群不同的A串,问B是哪些A串的子串。

    串匹配是一个很有研究价值的问题。事实上,我们还有后缀树,自动机等很多方法,这些算法都巧妙地运用了预处理,从而可以在线性的时间里解决字符串的匹配。我们以后来说。

    昨天发现一个特别晕的事,知道怎么去掉BitComet的广告吗?把界面语言设成英文就行了。
    还有,金山词霸和Dr.eye都可以去自杀了,Babylon素王道。

Matrix67原创
转贴请注明出处 
 


Trackback: http://tb.blog.csdn.net/TrackBack.aspx?PostId=1635175