18638人阅读
算法系列（37）
&&&&&&& 椭圆和直线、圆一样，是图形学领域中的一种常见图元，椭圆的生成算法（光栅转换算法）也是图形学软件中最常见的生成算法之一。在平面解析几何中，椭圆的方程可以描述为，其中是圆心坐标，和是椭圆的长短轴，特别的，当就是坐标中心点时，椭圆方程可以简化为。在计算机图形学中，椭圆图形也存在在点阵输出设备上显示或输出的问题，因此也需要一套光栅扫描转换算法。为了简化，我们先考虑圆心在原点的椭圆的生成，对于中心不是原点的椭圆，可以通过坐标的平移变换获得相应位置的椭圆。&&&&&&& 在进行扫描转换之前，需要了解一下椭圆的对称性，如图（）所示：图（）椭圆的对称性中心在原点。焦点在坐标轴上的标准椭圆具有轴对称、轴对称和原点对称特性，已知椭圆上第一象限的点坐标是，则椭圆在另外三个象限的对称点分别是、和。因此，只要画出第一象限的四分之一椭圆，就可以利用这三个对称性得到整个椭圆。在光栅设备上输出椭圆有很多种方法，可以根据直角平面坐标方程直接求解点坐标，利用极坐标方程求解，但是因为涉及到浮点数取整，效果都不好，一般都不使用直接求解的方式。本文就介绍几种计算机图形学中两种比较常用的椭圆生成方法：中点画椭圆算法和椭圆生成算法。中点画椭圆法&&&&&&& 中点在坐标原点，焦点在坐标轴上（轴对齐）的椭圆的平面集合方程是：，也可以转化为如下非参数化方程形式：（方程）无论是中点画线算法、中点画圆算法还是本节要介绍的中点画椭圆算法，对选择方向像素Δ增量还是方向像素Δ增量都是很敏感的。举个例子，如果某段圆弧上，方向上增量个像素时，方向上的增量如果，则比较适合用中点算法，如果方向上的增量，就会产生一些跳跃的点，最后生成的光栅位图圆弧会有一些突变的点，看起来好像不在圆弧上。因此，对于中点画圆弧算法，要区分出椭圆弧上哪段Δ增量变化显著，哪段Δ增量变化显著，然后区别对待。由于椭圆的对称性，我们只考虑第一象限的椭圆圆弧，如图（）所示：图（）第一象限椭圆弧示意图定义椭圆弧上某点的切线法向量如下：对方程分别求偏导和偏导，最后得到椭圆弧上点处的法向量是（）。的点是椭圆弧上的分界点。此点之上的部分（橙褐色部分）椭圆弧法向量的分量比较大，即：；此点之下的部分（蓝紫色部分）椭圆弧法向量的分量比较大，即：。&&&&&&& 对于图（）中橙褐色标识的上部区域，方向每变化个单位，方向变化大于一个单位，因此中点算法需要沿着方向步进画点，每次增量加，求的值。同理，对于图（）中蓝紫色标识的下部区域，中点算法沿着方向反向步进，每次减，求的值。先来讨论上部区域椭圆弧的生成，如图（）所示：图（）中点画椭圆算法对上部区域处理示意图假设当前位置是，则下一个可能的点就是点右边的点或右下方的点，取舍的方法取决于判别式，的定义如下：若，表示像素点和的中点在椭圆内，这时可取为下一个像素点。此时，，代入判别式得到：计算出的增量是。同理，若，表示像素点和的中点在椭圆外，这时应当取为下一个像素点。此时，，代入判别式得到：计算出的增量是。计算的增量的目的是减少计算量，提高算法效率，每次判断一个点时，不必完整的计算判别式，只需在上一次计算出的判别式上增加一个增量即可。&&&&&&& 接下来继续讨论下部区域椭圆弧的生成，如图（）所示：图（）中点画椭圆算法对下部区域处理示意图假设当前位置是，则下一个可能的点就是点左下方的点或下方的点，取舍的方法同样取决于判别式，的定义如下：若，表示像素点和的中点在椭圆内，这时可取为下一个像素点。此时，，代入判别式得到：计算出的增量是。同理，若，表示像素点和的中点在椭圆外，这时应当取为下一个像素点。此时，，代入判别式得到：计算出的增量是。&&&&&&& 中点画椭圆算法从点开始，第一个中点是，判别式的初始值是：上部区域生成算法的循环终止条件是：，下部区域的循环终止条件是，至此，中点画椭圆算法就可以完整给出了：20&void MP_Ellipse(int xc , int yc , int a, int b)21&{& 22&&&& double sqa = a * a;23&&&& double sqb = b * b;24&25&&&& double d = sqb + sqa * (-b + 0.25);26&&&& int x = 0; 27&&&& int y = b;28&&&& EllipsePlot(xc, yc, x, y);29&&&& while( sqb * (x + 1) & sqa * (y - 0.5))30&&&& {&&&& 31&&&& &&& if (d & 0)32&&&& &&& {33&&&& &&& &&& d += sqb * (2 * x + 3); 34&&&& &&& }35&&&& &&& else& 36&&&& &&& {& 37&&&& &&& &&& d += (sqb * (2 * x + 3) + sqa * (-2 * y + 2));38&&&& &&& &&& y--;&&& 39&&&& &&& } 40&&&& &&& x++;& 41&&&& &&& EllipsePlot(xc, yc, x, y);42&&&& }43&&&& d = (b * (x + 0.5)) * 2 + (a * (y - 1)) * 2 - (a * b) * 2;44&&&& while(y & 0)45&&&& {& 46&&&& &&& if (d & 0) 47&&&& &&& { 48&&&& &&& &&& d += sqb * (2 * x + 2) + sqa * (-2 * y + 3);49&&&& &&& &&& x++;& 50&&&& &&& }51&&&& &&& else& 52&&&& &&& {53&&&& &&& &&& d += sqa * (-2 * y + 3);& 54&&&& &&& }55&&&& &&& y--; 56&&&& &&& EllipsePlot(xc, yc, x, y);57&&&& }58&}EllipsePlot()函数利用椭圆的三个对称性，一次完成四个对称点的绘制，因为简单，此处就不再列出代码。算法&&&&&&& 中点画椭圆法中，计算判别式使用了浮点运算，影响了椭圆的生成效率。如果能将判别式规约到整数运算，则可以简化计算，提高效率。于是人们针对中点画椭圆法进行了多种改进，提出了很多种中点生成椭圆的整数型算法，椭圆生成算法就是其中之一。&&&&&&& 在生成椭圆上部区域时，以轴为步进方向，如图（）所示：图（）椭圆生成算法判别式每步进一个单位，就需要在判断保持不变还是也步进减，算法定义判别式为：如果，则取为下一个点，否则，取为下一个点。采用判别式，避免了中点算法因而引入的浮点运算，使得判别式规约为全整数运算，算法效率得到了很大的提升。根据椭圆方程，可以计算出和分别是：以作为椭圆上部区域的起点，将其代入判别式可以得到如下递推关系：&在生成椭圆下部区域时，以轴为步进方向，如图（）所示，每步进减一个单位，就需要在判断保持不变还是步进加一个单位，对于下部区域，计算出和分别是：以作为椭圆下部区域的起点，将其代入判别式可以得到如下递推关系：根据以上分析，椭圆生成算法的实现就比较简单了：&61&void Bresenham_Ellipse(int xc , int yc , int a, int b)&62&{ &63&&&& int sqa = a * a;&64&&&& int sqb = b * b;&65&&66&&&& int x = 0;&67&&&& int y = b;&68&&&& int d = 2 * sqb - 2 * b * sqa + sqa;&69&&&& EllipsePlot(xc, yc, x, y);&70&&&& int P_x = ROUND_INT( (double)sqa/sqrt((double)(sqa+sqb)) );&71&&&& while(x &= P_x)&72&&&& {&73&&&& &&& if(d & 0)&74&&&& &&& {&75&&&& &&& &&& d += 2 * sqb * (2 * x + 3);&76&&&& &&& }&77&&&& &&& else&78&&&& &&& {&79&&&& &&& &&& d += 2 * sqb * (2 * x + 3) - 4 * sqa * (y - 1);&80&&&& &&& &&& y--;&81&&&& &&& }&82&&&& &&& x++;&83&&&& &&& EllipsePlot(xc, yc, x, y);&84&&&& }&85&&86&&&& d = sqb * (x * x + x) + sqa * (y * y - y) - sqa * sqb;&87&&&& while(y &= 0)&88&&&& { &89&&&& &&& EllipsePlot(xc, yc, x, y);&90&&&& &&& y--;&91&&&& &&& if(d & 0)&92&&&& &&& { &93&&&& &&& &&& x++; &94&&&& &&& &&& d = d - 2 * sqa * y - sqa + 2 * sqb * x + 2 * sqb; &95&&&& &&& }&96&&&& &&& else&97&&&& &&& { &98&&&& &&& &&& d = d - 2 * sqa * y - sqa; &99&&&& &&& } 100&&&& } 101&}总结一下，本文介绍了两种计算机图形学中常见的椭圆生成算法，实际上还有很多种改进算法，包括提高效率，引入反走样技术等等，但那已经不是本文的重点，能把算法的基本原理讲清楚才是本文的目的。参考资料：【】计算几何：算法设计与分析周培德清华大学出版社年【】计算几何：算法与应用德贝尔赫（邓俊辉译）清华大学出版社年【】计算机图形学孙家广、杨常贵清华大学出版社年
参考知识库
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：1163258次
积分：12600
积分：12600
排名：第698名
原创：104篇
译文：12篇
评论：1120条
阅读：58185
文章：31篇
阅读：622343
(2)(1)(1)(2)(1)(2)(1)(1)(1)(1)(1)(2)(1)(1)(2)(1)(1)(5)(2)(4)(1)(1)(5)(1)(1)(2)(3)(1)(1)(1)(1)(1)(3)(1)(1)(5)(1)(1)(1)(1)(3)(1)(1)(2)(1)(1)(2)(3)(1)(3)(2)(3)(8)(8)(1)(2)(2)(1)(5)18796人阅读
算法系列（37）
3.6 用矢量的叉积判断直线段是否有交&&&&&&& &&&&&&& 矢量叉积计算的另一个常用用途是直线段求交。求交算法是计算机图形学的核心算法，也是体现速度和稳定性的重要标志，高效并且稳定的求交算法是任何一个CAD软件都必需要重点关注的。求交包含两层概念，一个是判断是否相交，另一个是求出交点。直线（段）的求交算法相对来说是比较简单的，首先来看看如何判断两直线段是否相交。&&&&&&& 常规的代数计算通常分三步，首先根据线段还原出两条线段所在直线的方程，然后联立方程组求出交点，最后再判断交点是否在线段区间上。常规的代数方法非常繁琐，每次都要解方程组求交点，特别是交点不在线段区间的情况，计算交点就是做无用功。计算几何方法判断直线段是否有交点通常分两个步骤完成，这两个步骤分别是快速排斥试验和跨立试验。假设要判断线段P1P2和线段Q1Q2是否有交点，则：&（1）&&&&&& 快速排斥试验&&& 设以线段 P1P2 为对角线的矩形为R1， 设以线段 Q1Q2 为对角线的矩形为R2，如果R1和R2不相交，则两线段不会有交点；&（2）&&&&&& 跨立试验。如果两线段相交，则两线段必然相互跨立对方，所谓跨立，指的是一条线段的两端点分别位于另一线段所在直线的两边。判断是否跨立，还是要用到矢量叉积的几何意义。以图3为例，若P1P2跨立Q1Q2 ，则矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的两侧，即：&( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) & 0&上式可改写成：&( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) & 0&当 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 时，说明线段P1P2和Q1Q2共线（但是不一定有交点）。同理判断Q1Q2跨立P1P2的依据是：&( Q1 - P1 ) × ( P2 - P1 ) * ( Q2 - P1 ) × ( P2 - P1 ) & 0&具体情况如下图所示：&图3 直线段跨立试验示意图&&&&&&&& 根据矢量叉积的几何意义，跨立试验只能证明线段的两端点位于另一个线段所在直线的两边，但是不能保证是在另一直线段的两端，因此，跨立试验只是证明两条线段有交点的必要条件，必需和快速排斥试验一起才能组成直线段相交的充分必要条件。根据以上分析，两条线段有交点的完整判断依据就是：1)以两条线段为对角线的两个矩形有交集；2)两条线段相互跨立。&&&&&&& 判断直线段跨立用计算叉积算法的CrossProduct()函数即可，还需要一个判断两个矩形是否有交的算法。矩形求交也是最简单的求交算法之一，原理就是根据两个矩形的最大、最小坐标判断。图4展示了两个矩形没有交集的各种情况：&图4 矩形没有交集的几种情况&图5展示了两个矩形有交集的各种情况：图5 矩形有交集的几种情况&从图4和图5可以分析出来两个矩形有交集的几何坐标规律，就是在x坐标方向和y坐标方向分别满足最大值最小值法则，简单解释这个法则就是每个矩形在每个方向上的坐标最大值都要大于另一个矩形在这个坐标方向上的坐标最小值，否则在这个方向上就不能保证一定有位置重叠。由以上分析，判断两个矩形是否相交的算法就可以如下实现：186&bool IsRectIntersect(const Rect& rc1, const Rect& rc2)187&{188&&&& return ( (std::max(rc1.p1.x, rc1.p2.x) &= std::min(rc2.p1.x, rc2.p2.x))189&&&& &&& &&&& && (std::max(rc2.p1.x, rc2.p2.x) &= std::min(rc1.p1.x, rc1.p2.x))190&&&& &&& &&&& && (std::max(rc1.p1.y, rc1.p2.y) &= std::min(rc2.p1.y, rc2.p2.y))191&&&& &&& &&&& && (std::max(rc2.p1.y, rc2.p2.y) &= std::min(rc1.p1.y, rc1.p2.y)) );192&}&&&&&&&& 完成了排斥试验和跨立试验的算法，最后判断直线段是否有交点的算法就水到渠成了：204&bool IsLineSegmentIntersect(const LineSeg& ls1, const LineSeg& ls2)205&{206&&&& if(IsLineSegmentExclusive(ls1, ls2)) //排斥实验207&&&& {208&&&& &&& return false;209&&&& }210&&&& //( P1 - Q1 ) ×'a1?( Q2 - Q1 )211&&&& double p1xq = CrossProduct(ls1.ps.x - ls2.ps.x, ls1.ps.y - ls2.ps.y, 212&&&& &&& &&& &&& &&& &&& &&& && ls2.pe.x - ls2.ps.x, ls2.pe.y - ls2.ps.y); 213&&&& //( P2 - Q1 ) ×'a1?( Q2 - Q1 )214&&&& double p2xq = CrossProduct(ls1.pe.x - ls2.ps.x, ls1.pe.y - ls2.ps.y, 215&&&& &&& &&& &&& &&& &&& &&& && ls2.pe.x - ls2.ps.x, ls2.pe.y - ls2.ps.y); 216&217&&&& //( Q1 - P1 ) ×'a1?( P2 - P1 )218&&&& double q1xp = CrossProduct(ls2.ps.x - ls1.ps.x, ls2.ps.y - ls1.ps.y, 219&&&& &&& &&& &&& &&& &&& &&& && ls1.pe.x - ls1.ps.x, ls1.pe.y - ls1.ps.y); 220&&&& //( Q2 - P1 ) ×'a1?( P2 - P1 )221&&&& double q2xp = CrossProduct(ls2.pe.x - ls1.ps.x, ls2.pe.y - ls1.ps.y, 222&&&& &&& &&& &&& &&& &&& &&& && ls1.pe.x - ls1.ps.x, ls1.pe.y - ls1.ps.y); 223&224&&&& //跨立实验225&&&& return ( (p1xq * p2xq &= 0.0) && (q1xp * q2xp &= 0.0) ); 226&}&IsLineSegmentExclusive()函数就是调用IsRectIntersect()函数根据结果做排斥判断，此处不再列出代码。&3.7 点和多边形关系的算法实现&&&&&&&& 好了，现在我们已经了解了矢量叉积的意义，以及判断直线段是否有交点的算法，现在回过头看看文章开始部分的讨论的问题：如何判断一个点是否在多边形内部？根据射线法的描述，其核心是求解从P点发出的射线与多边形的边是否有交点。注意，这里说的是射线，而我们前面讨论的都是线段，好像不适用吧？没错，确实是不适用，但是我要介绍一种用计算机解决问题时常用的建模思想，应用了这种思想之后，我们前面讨论的方法就适用了。什么思想呢？就是根据问题域的规模和性质抽象和简化模型的思想，这可不是故弄玄虚，说说具体的思路吧。&&&&&&& 计算机是不能表示无穷大和无穷小，计算机处理的每一个数都有确定的值，而且必须有确定的值。我们面临的问题域是整个实数空间的坐标系，在每个维度上都是从负无穷到正无穷，比如射线，就是从坐标系中一个明确的点到无穷远处的连线。这就有点为难计算机了，为此我们需要简化问题的规模。假设问题中多边形的每个点的坐标都不会超过(-10000.0, +10000.0)区间（比如我们常见的图形输出设备都有大小的限制），我们就可以将问题域简化为(-10000.0, +10000.0)区间内的一小块区域，对于这块区域来说，&= 10000.0就意味着无穷远。你肯定已经明白了，数学模型经过简化后，算法中提到的射线就可以理解为从模型边界到内部点P之间的线段，前面讨论的关于线段的算法就可以使用了。&&&&&&& 射线法的基本原理是判断由P点发出的射线与多边形的交点个数，交点个数是奇数表示P点在多边形内（在多边形的边上也视为在多边形内部的特殊情况），正常情况下经过点P的射线应该如图6（a）所示那样，但是也可能碰到多种非正常情况，比如刚好经过多边形一个定点的情况，如图6 （b），这会被误认为和两条边都有交点，还可能与某一条边共线如图6 （c）和（d），共线就有无穷多的交点，导致判断规则失效。还要考虑凹多边形的情况，如图6（e）。&&图6 射线法可能遇到的各种交点情况&&&&&&&&& 针对这些特殊情况，在对多边形的每条边进行判断时，要考虑以下这些特殊情况，假设当前处理的边是P1P2，则有以下原则：&1）如果点P在边P1P2上，则直接判定点P在多边形内；2）如果从P发出的射线正好穿过P1或者P2，那么这个交点会被算作2次（因为在处理以P1或P2为端点的其它边时可能已经计算过这个点了），对这种情况的处理原则是：如果P的y坐标与P1、P2中较小的y坐标相同，则忽略这个交点；3）如果从P发出的射线与P1P2平行，则忽略这条边；&&&&&&&& 对于第三个原则，需要判断两条直线是否平行，通常的方法是计算两条直线的斜率，但是本算法因为只涉及到直线段（射线也被模型简化为长线段了），就简化了很多，判断直线是否水平，只要比较一下线段起始点的y坐标是否相等就行了，而判断直线是否垂直，也只要比较一下线段起始点的x坐标是否相等就行了。&&&&&&& 应用以上原则后，扫描线法判断点是否在多边形内的算法流程就完整了，图7就是算法的流程图：&&&&&&&& 最终扫描线法判断点是否在多边形内的算法实现如下：228&bool IsPointInPolygon(const Polygon& py, const Point& pt)229&{230&&&& assert(py.IsValid()); /*只考虑正常的多边形，边数&=3*/231&232&&&& int count = 0;233&&&& LineSeg ll = LineSeg(pt, Point(-INFINITE, pt.y)); /*射线L*/234&&&& for(int i = 0; i & py.GetPolyCount(); i++)235&&&& {236&&&& &&& /*当前点和下一个点组成线段P1P2*/237&&&& &&& LineSeg pp = LineSeg(py.pts[i], py.pts[(i + 1) % py.GetPolyCount()]); 238&&&& &&& if(IsPointOnLineSegment(pp, pt))239&&&& &&& {240&&&& &&& &&& return true;241&&&& &&& }242&243&&&& &&& if(!pp.IsHorizontal())244&&&& &&& {245&&&& &&& &&& if((IsSameFloatValue(pp.ps.y, pt.y)) && (pp.ps.y & pp.pe.y))246&&&& &&& &&& {247&&&& &&& &&& &&& count++;248&&&& &&& &&& }249&&&& &&& &&& else if((IsSameFloatValue(pp.pe.y, pt.y)) && (pp.pe.y & pp.ps.y))250&&&& &&& &&& {251&&&& &&& &&& &&& count++;252&&&& &&& &&& }253&&&& &&& &&& else254&&&& &&& &&& {255&&&& &&& &&& &&& if(IsLineSegmentIntersect(pp, ll))256&&&& &&& &&& &&& {257&&&& &&& &&& &&& &&& count++;258&&&& &&& &&& &&& }259&&&& &&& &&& }260&&&& &&& }261&&&& }262&263&&&& return ((count % 2) == 1);264&}&&&&&&&& 在图形学领域实施的真正工程代码，通常还会增加一个多边形的外包矩形快速判断，对点根本就不在多边形周围的情况做快速排除，提高算法效率。这又涉及到求多边形外包矩形的算法，这个算法也很简单，就是遍历多边形的所有节点，找出各个坐标方向上的最大最小值。以下就是求多边形外包矩形的算法：266&void GetPolygonEnvelopRect(const Polygon& py, Rect& rc)267&{268&&&& assert(py.IsValid()); /*只考虑正常的多边形，边数&=3*/269&270&&&& double minx = py.pts[0].x;271&&&& double maxx = py.pts[0].x;272&&&& double miny = py.pts[0].y;273&&&& double maxy = py.pts[0].y;274&&&& for(int i = 1; i & py.GetPolyCount(); i++)275&&&& {276&&&& &&& if(py.pts[i].x & minx)277&&&& &&& &&& minx = py.pts[i].x;278&&&& &&& if(py.pts[i].x & maxx)279&&&& &&& &&& maxx = py.pts[i].x;280&&&& &&& if(py.pts[i].y & miny)281&&&& &&& &&& miny = py.pts[i].y;282&&&& &&& if(py.pts[i].y & maxy)283&&&& &&& &&& maxy = py.pts[i].y;284&&&& }285&286&&&& rc = Rect(minx, miny, maxx, maxy);287&}&&&&&&&&& 除了扫描线法，还可以通过多边形边的法矢量方向、多边形面积以及角度和等方法判断点与多边形的关系。但是这些算法要么只支持凸多边形，要么需要复杂的三角函数运算（多边形边数小于44时，可采用近似公式计算夹角和，避免三角函数运算），使用的范围有限，只有扫描线法被广泛应用。&&&&&&& 至此，本文的内容已经完结，以上通过对点与矩形、点与圆、点与直线以及点与多边形位置关系判断算法的讲解，向大家展示了几种常见的计算几何算法实现，用简短而易懂的代码剖析了这些算法的实质。下一篇将介绍计算机图形学中最基本的直线生成算法。&&&&参考资料：&【1】计算几何：算法设计与分析 周培德& 清华大学出版社 2005年【2】计算几何：算法与应用 德贝尔赫（邓俊辉译）& 清华大学出版社 2005年【3】算法导论 Thomas H.Cormen等（潘金贵等译） 机械工业出版社 2006年【4】计算机图形学 孙家广、杨常贵 清华大学出版社 1995年&&&
参考知识库
* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
访问：1163260次
积分：12600
积分：12600
排名：第698名
原创：104篇
译文：12篇
评论：1120条
阅读：58185
文章：31篇
阅读：622343
(2)(1)(1)(2)(1)(2)(1)(1)(1)(1)(1)(2)(1)(1)(2)(1)(1)(5)(2)(4)(1)(1)(5)(1)(1)(2)(3)(1)(1)(1)(1)(1)(3)(1)(1)(5)(1)(1)(1)(1)(3)(1)(1)(2)(1)(1)(2)(3)(1)(3)(2)(3)(8)(8)(1)(2)(2)(1)(5)posts - 90,&
comments - 33,&
trackbacks - 0
计算机的出现使得很多原本十分繁琐的工作得以大幅度简化，但是也有一些在人们直观看来很容易的问题却需要拿出一套并不简单的通用解决方案，比如几何问题。作为计算机科学的一个分支，计算几何主要研究解决几何问题的算法。在现代工程和数学领域，计算几何在图形学、机器人技术、超大规模集成电路设计和统计等诸多领域有着十分重要的应用。在本文中，我们将对计算几何常用的基本算法做一个全面的介绍，希望对您了解并应用计算几何的知识解决问题起到帮助。
(本文整理的计算几何基本概念和常用算法包括如下内容：)
1. 矢量的概念
2. 矢量加减法
3. 矢量叉积
4. 折线段的拐向判断
5. 判断点是否在线段上
6. 判断两线段是否相交
7. 判断线段和直线是否相交
8. 判断矩形是否包含点
9. 判断线段、折线、多边形是否在矩形中
10. 判断矩形是否在矩形中
11. 判断圆是否在矩形中
12. 判断点是否在多边形中
13. 判断线段是否在多边形内
14. 判断折线是否在多边形内
15. 判断多边形是否在多边形内
16. 判断矩形是否在多边形内
17. 判断圆是否在多边形内
18. 判断点是否在圆内
19. 判断线段、折线、矩形、多边形是否在圆内
20. 判断圆是否在圆内
21. 计算点到线段的最近点
22. 计算点到折线、矩形、多边形的最近点
23. 计算点到圆的最近距离及交点坐标
24. 计算两条共线的线段的交点
25. 计算线段或直线与线段的交点
26. 求线段或直线与折线、矩形、多边形的交点
27. 求线段或直线与圆的交点
28. 凸包的概念
29. 凸包的求法
三、算法介绍
1. 矢量的概念：
如果一条线段的端点是有次序之分的，我们把这种线段成为有向线段(directed segment)。如果有向线段p1p2的起点p1在坐标原点，我们可以把它称为矢量(vector)p2。
2. 矢量加减法：
设二维矢量P = ( x1, y1 )，Q = ( x2 , y2 )，则矢量加法定义为： P + Q = ( x1 + x2 , y1 + y2 )，同样的，矢量减法定义为： P - Q = ( x1 - x2 , y1 - y2 )。显然有性质 P + Q = Q + P，P - Q = - ( Q - P )。
3. 矢量叉积：
计算矢量叉积是与直线和线段相关算法的核心部分。设矢量P = ( x1, y1 )，Q = ( x2, y2 )，则矢量叉积定义为由(0,0)、p1、p2和p1+p2所组成的平行四边形的带符号的面积，即：P × Q = x1*y2 - x2*y1，其结果是一个标量。显然有性质 P × Q = - ( Q × P ) 和 P × ( - Q ) = - ( P × Q )。一般在不加说明的情况下，本文下述算法中所有的点都看作矢量，两点的加减法就是矢量相加减，而点的乘法则看作矢量叉积。
叉积的一个非常重要性质是可以通过它的符号判断两矢量相互之间的顺逆时针关系：
若 P × Q & 0 , 则P在Q的顺时针方向。
若 P × Q & 0 , 则P在Q的逆时针方向。
若 P × Q = 0 , 则P与Q共线，但可能同向也可能反向。
4. 折线段的拐向判断：
折线段的拐向判断方法可以直接由矢量叉积的性质推出。对于有公共端点的线段p0p1和p1p2，通过计算(p2 - p0) × (p1 - p0)的符号便可以确定折线段的拐向：
若(p2 - p0) × (p1 - p0) & 0,则p0p1在p1点拐向右侧后得到p1p2。
若(p2 - p0) × (p1 - p0) & 0,则p0p1在p1点拐向左侧后得到p1p2。
若(p2 - p0) × (p1 - p0) = 0,则p0、p1、p2三点共线。
具体情况可参考下图：
5. 判断点是否在线段上：
设点为Q，线段为P1P2 ，判断点Q在该线段上的依据是：( Q - P1 ) × ( P2 - P1 ) = 0 且 Q 在以 P1，P2为对角顶点的矩形内。前者保证Q点在直线P1P2上，后者是保证Q点不在线段P1P2的延长线或反向延长线上，对于这一步骤的判断可以用以下过程实现：
ON-SEGMENT(pi,pj,pk)
if min(xi,xj) &= xk &= max(xi,xj) and min(yi,yj) &= yk &= max(yi,yj)
特别要注意的是，由于需要考虑水平线段和垂直线段两种特殊情况，min(xi,xj)&=xk&=max(xi,xj)和min(yi,yj)&=yk&=max(yi,yj)两个条件必须同时满足才能返回真值。
6. 判断两线段是否相交：
我们分两步确定两条线段是否相交：
(1)快速排斥试验
设以线段 P1P2 为对角线的矩形为R， 设以线段 Q1Q2 为对角线的矩形为T，如果R和T不相交，显然两线段不会相交。
(2)跨立试验
如果两线段相交，则两线段必然相互跨立对方。若P1P2跨立Q1Q2 ，则矢量 ( P1 - Q1 ) 和( P2 - Q1 )位于矢量( Q2 - Q1 ) 的两侧，即( P1 - Q1 ) × ( Q2 - Q1 ) * ( P2 - Q1 ) × ( Q2 - Q1 ) & 0。上式可改写成( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) & 0。当 ( P1 - Q1 ) × ( Q2 - Q1 ) = 0 时，说明 ( P1 - Q1 ) 和 ( Q2 - Q1 )共线，但是因为已经通过快速排斥试验，所以 P1 一定在线段 Q1Q2上；同理，( Q2 - Q1 ) ×(P2 - Q1 ) = 0 说明 P2 一定在线段 Q1Q2上。所以判断P1P2跨立Q1Q2的依据是：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) &= 0。同理判断Q1Q2跨立P1P2的依据是：( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) &= 0。具体情况如下图所示：
在相同的原理下，对此算法的具体的实现细节可能会与此有所不同，除了这种过程外，大家也可以参考《算法导论》上的实现。
7. 判断线段和直线是否相交：
有了上面的基础，这个算法就很容易了。如果线段P1P2和直线Q1Q2相交，则P1P2跨立Q1Q2，即：( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) &= 0。
8. 判断矩形是否包含点：
只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。
9. 判断线段、折线、多边形是否在矩形中：
因为矩形是个凸集，所以只要判断所有端点是否都在矩形中就可以了。
10. 判断矩形是否在矩形中：
只要比较左右边界和上下边界就可以了。
11. 判断圆是否在矩形中：
很容易证明，圆在矩形中的充要条件是：圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。
12. 判断点是否在多边形中：
判断点P是否在多边形中是计算几何中一个非常基本但是十分重要的算法。以点P为端点，向左方作射线L，由于多边形是有界的，所以射线L的左端一定在多边形外，考虑沿着L从无穷远处开始自左向右移动，遇到和多边形的第一个交点的时候，进入到了多边形的内部，遇到第二个交点的时候，离开了多边形，……所以很容易看出当L和多边形的交点数目C是奇数的时候，P在多边形内，是偶数的话P在多边形外。
但是有些特殊情况要加以考虑。如图下图 (a)(b)(c)(d)所示。在图(a)中，L和多边形的顶点相交，这时候交点只能计算一个；在图(b)中，L和多边形顶点的交点不应被计算；在图 (c)和(d) 中，L和多边形的一条边重合，这条边应该被忽略不计。如果L和多边形的一条边重合，这条边应该被忽略不计。
为了统一起见，我们在计算射线L和多边形的交点的时候，1。对于多边形的水平边不作考虑；2。对于多边形的顶点和L相交的情况，如果该顶点是其所属的边上纵坐标较大的顶点，则计数，否则忽略；3。对于P在多边形边上的情形，直接可判断P属于多边行。由此得出算法的伪代码如下：
count ← 0;
以P为端点，作从右向左的射线L;
for 多边形的每条边s
do if P在边s上
if s不是水平的
then if s的一个端点在L上
if 该端点是s两端点中纵坐标较大的端点
then count ← count+1
else if s和L相交
then count ← count+1;
if count mod 2 = 1
其中做射线L的方法是：设P'的纵坐标和P相同，横坐标为正无穷大（很大的一个正数），则P和P'就确定了射线L。
判断点是否在多边形中的这个算法的时间复杂度为O(n)。
另外还有一种算法是用带符号的三角形面积之和与多边形面积进行比较，这种算法由于使用浮点数运算所以会带来一定误差，不推荐大家使用。
13. 判断线段是否在多边形内：
线段在多边形内的一个必要条件是线段的两个端点都在多边形内，但由于多边形可能为凹，所以这不能成为判断的充分条件。如果线段和多边形的某条边内交（两线段内交是指两线段相交且交点不在两线段的端点），因为多边形的边的左右两侧分属多边形内外不同部分，所以线段一定会有一部分在多边形外(见图a)。于是我们得到线段在多边形内的第二个必要条件：线段和多边形的所有边都不内交。
线段和多边形交于线段的两端点并不会影响线段是否在多边形内；但是如果多边形的某个顶点和线段相交，还必须判断两相邻交点之间的线段是否包含于多边形内部（反例见图b)。
因此我们可以先求出所有和线段相交的多边形的顶点，然后按照X-Y坐标排序(X坐标小的排在前面，对于X坐标相同的点，Y坐标小的排在前面，这种排序准则也是为了保证水平和垂直情况的判断正确)，这样相邻的两个点就是在线段上相邻的两交点，如果任意相邻两点的中点也在多边形内，则该线段一定在多边形内。
证明如下：
命题1：
如果线段和多边形的两相邻交点P1 ，P2的中点P' 也在多边形内，则P1, P2之间的所有点都在多边形内。
假设P1,P2之间含有不在多边形内的点，不妨设该点为Q，在P1, P'之间，因为多边形是闭合曲线，所以其内外部之间有界，而P1属于多边行内部，Q属于多边性外部，P'属于多边性内部，P1-Q-P'完全连续，所以 P1Q和QP'一定跨越多边形的边界，因此在P1,P'之间至少还有两个该线段和多边形的交点，这和P1P2是相邻两交点矛盾，故命题成立。证毕。
由命题1直接可得出推论：
设多边形和线段PQ的交点依次为P1,P2,……Pn，其中Pi和Pi+1是相邻两交点，线段PQ在多边形内的充要条件是：P，Q在多边形内且对于i =1, 2,……, n-1，Pi ,Pi+1的中点也在多边形内。
在实际编程中，没有必要计算所有的交点，首先应判断线段和多边形的边是否内交，倘若线段和多边形的某条边内交则线段一定在多边形外；如果线段和多边形的每一条边都不内交，则线段和多边形的交点一定是线段的端点或者多边形的顶点，只要判断点是否在线段上就可以了。
至此我们得出算法如下：
if 线端PQ的端点不都在多边形内
点集pointSet初始化为空;
for 多边形的每条边s
do if 线段的某个端点在s上
then 将该端点加入pointS
else if s的某个端点在线段PQ上
then 将该端点加入pointS
else if s和线段PQ相交 // 这时候已经可以肯定是内交了
将pointSet中的点按照X-Y坐标排序;
for pointSet中每两个相邻点 pointSet[i] , pointSet[ i+1]
do if pointSet[i] , pointSet[ i+1] 的中点不在多边形中
这个过程中的排序因为交点数目肯定远小于多边形的顶点数目n，所以最多是常数级的复杂度，几乎可以忽略不计。因此算法的时间复杂度也是O(n)。
14. 判断折线是否在多边形内：
只要判断折线的每条线段是否都在多边形内即可。设折线有m条线段，多边形有n个顶点，则该算法的时间复杂度为O(m*n)。
15. 判断多边形是否在多边形内：
只要判断多边形的每条边是否都在多边形内即可。判断一个有m个顶点的多边形是否在一个有n个顶点的多边形内复杂度为O(m*n)。
16. 判断矩形是否在多边形内：
将矩形转化为多边形，然后再判断是否在多边形内。
17 判断圆是否在多边形内：
只要计算圆心到多边形的每条边的最短距离，如果该距离大于等于圆半径则该圆在多边形内。计算圆心到多边形每条边最短距离的算法在后文阐述。
18. 判断点是否在圆内：
计算圆心到该点的距离，如果小于等于半径则该点在圆内。
19. 判断线段、折线、矩形、多边形是否在圆内：
因为圆是凸集，所以只要判断是否每个顶点都在圆内即可。
20. 判断圆是否在圆内：
设两圆为O1,O2，半径分别为r1, r2，要判断O2是否在O1内。先比较r1，r2的大小，如果r1&r2则O2不可能在O1内；否则如果两圆心的距离大于r1 - r2 ，则O2不在O1内；否则O2在O1内。
21. 计算点到线段的最近点：
如果该线段平行于X轴（Y轴），则过点point作该线段所在直线的垂线，垂足很容易求得，然后计算出垂足，如果垂足在线段上则返回垂足，否则返回离垂足近的端点；如果该线段不平行于X轴也不平行于Y轴，则斜率存在且不为0。设线段的两端点为pt1和pt2，斜率为：k = ( pt2.y - pt1. y ) / (pt2.x - pt1.x );该直线方程为：y = k* ( x - pt1.x) + pt1.y。其垂线的斜率为 - 1 / k，垂线方程为：y = (-1/k) * (x - point.x) + point.y 。
联立两直线方程解得：x = ( k^2 * pt1.x + k * (point.y - pt1.y ) + point.x ) / ( k^2 + 1) ，y = k * ( x - pt1.x) + pt1.y;然后再判断垂足是否在线段上，如果在线段上则返回垂足；如果不在则计算两端点到垂足的距离，选择距离垂足较近的端点返回。
22. 计算点到折线、矩形、多边形的最近点：
只要分别计算点到每条线段的最近点，记录最近距离，取其中最近距离最小的点即可。
23. 计算点到圆的最近距离及交点坐标：
如果该点在圆心，因为圆心到圆周任一点的距离相等，返回UNDEFINED。
连接点P和圆心O，如果PO平行于X轴，则根据P在O的左边还是右边计算出最近点的横坐标为centerPoint.x - radius 或 centerPoint.x + radius。如果PO平行于Y轴，则根据P在O的上边还是下边计算出最近点的纵坐标为 centerPoint.y -+radius或 centerPoint.y - radius。如果PO不平行于X轴和Y轴，则PO的斜率存在且不为0，这时直线PO斜率为k = （ P.y - O.y ）/ ( P.x - O.x )。直线PO的方程为：y = k * ( x - P.x) + P.y。设圆方程为:(x - O.x ) ^2 + ( y - O.y ) ^2 = r ^2，联立两方程组可以解出直线PO和圆的交点，取其中离P点较近的交点即可。
24. 计算两条共线的线段的交点：
对于两条共线的线段，它们之间的位置关系有下图所示的几种情况。图(a)中两条线段没有交点；图 (b) 和 (d) 中两条线段有无穷焦点；图 (c) 中两条线段有一个交点。设line1是两条线段中较长的一条，line2是较短的一条，如果line1包含了line2的两个端点，则是图(d)的情况，两线段有无穷交点；如果line1只包含line2的一个端点，那么如果line1的某个端点等于被line1包含的line2的那个端点，则是图(c) 的情况，这时两线段只有一个交点，否则就是图(b)的情况，两线段也是有无穷的交点；如果line1不包含line2的任何端点，则是图(a)的情况，这时两线段没有交点。
25. 计算线段或直线与线段的交点:
设一条线段为L0 = P1P2，另一条线段或直线为L1 = Q1Q2 ，要计算的就是L0和L1的交点。
1． 首先判断L0和L1是否相交（方法已在前文讨论过），如果不相交则没有交点，否则说明L0和L1一定有交点，下面就将L0和L1都看作直线来考虑。
2． 如果P1和P2横坐标相同，即L0平行于Y轴
a) 若L1也平行于Y轴，
i. 若P1的纵坐标和Q1的纵坐标相同，说明L0和L1共线，假如L1是直线的话他们有无穷的交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）；
ii. 否则说明L0和L1平行，他们没有交点；
b) 若L1不平行于Y轴，则交点横坐标为P1的横坐标，代入到L1的直线方程中可以计算出交点纵坐标；
3． 如果P1和P2横坐标不同，但是Q1和Q2横坐标相同，即L1平行于Y轴，则交点横坐标为Q1的横坐标，代入到L0的直线方程中可以计算出交点纵坐标；
4． 如果P1和P2纵坐标相同，即L0平行于X轴
a) 若L1也平行于X轴，
i. 若P1的横坐标和Q1的横坐标相同，说明L0和L1共线，假如L1是直线的话他们有无穷的交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）；
ii. 否则说明L0和L1平行，他们没有交点；
b) 若L1不平行于X轴，则交点纵坐标为P1的纵坐标，代入到L1的直线方程中可以计算出交点横坐标；
5． 如果P1和P2纵坐标不同，但是Q1和Q2纵坐标相同，即L1平行于X轴，则交点纵坐标为Q1的纵坐标，代入到L0的直线方程中可以计算出交点横坐标；
6． 剩下的情况就是L1和L0的斜率均存在且不为0的情况
a) 计算出L0的斜率K0，L1的斜率K1 ；
b) 如果K1 = K2
i. 如果Q1在L0上，则说明L0和L1共线，假如L1是直线的话有无穷交点，假如L1是线段的话可用"计算两条共线线段的交点"的算法求他们的交点（该方法在前文已讨论过）；
ii. 如果Q1不在L0上，则说明L0和L1平行，他们没有交点。
c) 联立两直线的方程组可以解出交点来
这个算法并不复杂，但是要分情况讨论清楚，尤其是当两条线段共线的情况需要单独考虑，所以在前文将求两条共线线段的算法单独写出来。另外，一开始就先利用矢量叉乘判断线段与线段（或直线）是否相交，如果结果是相交，那么在后面就可以将线段全部看作直线来考虑。需要注意的是，我们可以将直线或线段方程改写为ax+by+c=0的形式，这样一来上述过程的部分步骤可以合并，缩短了代码长度，但是由于先要求出参数，这种算法将花费更多的时间。
26. 求线段或直线与折线、矩形、多边形的交点：
分别求与每条边的交点即可。
27. 求线段或直线与圆的交点:
设圆心为O，圆半径为r，直线（或线段）L上的两点为P1,P2。
1. 如果L是线段且P1，P2都包含在圆O内，则没有交点；否则进行下一步。
2. 如果L平行于Y轴，
a) 计算圆心到L的距离dis；
b) 如果dis & r 则L和圆没有交点；
c) 利用勾股定理，可以求出两交点坐标，但要注意考虑L和圆的相切情况。
3. 如果L平行于X轴，做法与L平行于Y轴的情况类似；
4. 如果L既不平行X轴也不平行Y轴，可以求出L的斜率K，然后列出L的点斜式方程，和圆方程联立即可求解出L和圆的两个交点；
5. 如果L是线段，对于2，3，4中求出的交点还要分别判断是否属于该线段的范围内。
28. 凸包的概念：
点集Q的凸包(convex hull)是指一个最小凸多边形，满足Q中的点或者在多边形边上或者在其内。下图中由红色线段表示的多边形就是点集Q={p0,p1,...p12}的凸包。
29. 凸包的求法：
现在已经证明了凸包算法的时间复杂度下界是O(n*logn),但是当凸包的顶点数h也被考虑进去的话，Krikpatrick和Seidel的剪枝搜索算法可以达到O(n*logh)，在渐进意义下达到最优。最常用的凸包算法是Graham扫描法和Jarvis步进法。本文只简单介绍一下Graham 扫描法，其正确性的证明和Jarvis步进法的过程大家可以参考《算法导论》。
对于一个有三个或以上点的点集Q，Graham扫描法的过程如下：
令p0为Q中Y-X坐标排序下最小的点
设&p1,p2,...pm&为对其余点按以p0为中心的极角逆时针排序所得的点集（如果有多个点有相同的极角，除了距p0最远的点外全部移除
for i ← 3 to m
do while 由S的栈顶元素的下一个元素、S的栈顶元素以及pi构成的折线段不拐向左侧
此过程执行后，栈S由底至顶的元素就是Q的凸包顶点按逆时针排列的点序列。需要注意的是，我们对点按极角逆时针排序时，并不需要真正求出极角，只需要求出任意两点的次序就可以了。而这个步骤可以用前述的矢量叉积性质实现。
㈠ 点的基本运算
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
/* 需要包含的头文件 */
#include&&cmath &
/* 常用的常量定义 */
const&double&INF = 1E200
const&double&EP = 1E-10
const&int&MAXV = 300
const&double&PI = 3.
/* 基本几何结构 */
struct&POINT
&&&&double&x;
&&&&double&y; POINT(double&a=0,&double&b=0) { x=a; y=b;}&//constructor
struct&LINESEG
&&& POINT LINESEG(POINT a, POINT b) { s=a; e=b;}
&&& LINESEG() { }
struct&LINE&// 直线的解析方程 a*x+b*y+c=0 为统一表示，约定 a &= 0
&&&&double&a;
&&&&double&b;
&&&&double&c; LINE(double&d1=1,&double&d2=-1,&double&d3=0) {a=d1; b=d2; c=d3;}
/********************/
* 点的基本运算 *
/********************/
double&dist(POINT p1,POINT p2)&// 返回两点之间欧氏距离
&&&&return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) );
bool&equal_point(POINT p1,POINT p2)&// 判断两个点是否重合
&&&&return&( (abs(p1.x-p2.x)&EP)&&(abs(p1.y-p2.y)&EP) );
/******************************************************************************
r=multiply(sp,ep,op),得到(sp-op)*(ep-op)的叉积
r&0:ep在矢量opsp的逆时针方向；
r=0：opspep三点共线；
r&0:ep在矢量opsp的顺时针方向
*******************************************************************************/
double&multiply(POINT sp,POINT ep,POINT op)
&&&&return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y));
/*******************************************************************************
r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积，如果两个矢量都非零矢量
r&0:两矢量夹角为锐角；r=0：两矢量夹角为直角；r&0:两矢量夹角为钝角
*******************************************************************************/
double&dotmultiply(POINT p1,POINT p2,POINT p0)
&&&&return&((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y));
/* 判断点p是否在线段l上，条件：(p在线段l所在的直线上)&& (点p在以线段l为对角线的矩形内) */
bool&online(LINESEG l,POINT p)
&&&&return((multiply(l.e,p,l.s)==0)
&&&&&&& &&( ( (p.x-l.s.x)*(p.x-l.e.x)&=0 )&&( (p.y-l.s.y)*(p.y-l.e.y)&=0 ) ) );
// 返回点p以点o为圆心逆时针旋转alpha(单位：弧度)后所在的位置
POINT rotate(POINT o,double&alpha,POINT p)
&&& p.x-=o.x;
&&& p.y-=o.y;
&&& tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x;
&&& tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y;
&&&&return&
/* 返回顶角在o点，起始边为os，终止边为oe的夹角(单位：弧度)
角度小于pi，返回正值
角度大于pi，返回负值
可以用于求线段之间的夹角
double&angle(POINT o,POINT s,POINT e)
&&&&double&cosfi,fi,
&&&&double&dsx = s.x - o.x;
&&&&double&dsy = s.y - o.y;
&&&&double&dex = e.x - o.x;
&&&&double&dey = e.y - o.y;
&&& cosfi=dsx*dex+dsy*
&&& norm=(dsx*dsx+dey*dey)*(dex*dex+dey*dey);
&&& cosfi /= sqrt( norm );
&&&&if&(cosfi &= 1.0 )&return&0;
&&&&if&(cosfi &= -1.0 )&return&-3.1415926;
&&& fi=acos(cosfi);
&&&&if&(dsx*dey-dsy*dex&0)&return&&// 说明矢量os 在矢量 oe的顺时针方向
&&&&return&-
/*****************************/
* 线段及直线的基本运算 *
/*****************************/
/* 判断点与线段的关系,用途很广泛
本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足
& AC dot AB
& r = ---------
& ||AB||^2
& (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
& = -------------------------------
&&& r has the following meaning:
&&&&& r=0 P = A
&&&&& r=1 P = B
&&&&& r&0 P is on the backward extension of AB
&&&&& r&1 P is on the forward extension of AB
&&&&& 0&r&1 P is interior to AB
double&relation(POINT p,LINESEG l)
&&& LINESEG
&&& tl.s=l.s;
&&& tl.e=p;
&&&&return&dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e));
// 求点C到线段AB所在直线的垂足 P
POINT perpendicular(POINT p,LINESEG l)
&&&&double&r=relation(p,l);
&&& tp.x=l.s.x+r*(l.e.x-l.s.x);
&&& tp.y=l.s.y+r*(l.e.y-l.s.y);
&&&&return&
/* 求点p到线段l的最短距离,并返回线段上距该点最近的点np
注意：np是线段l上到点p最近的点，不一定是垂足 */
double&ptolinesegdist(POINT p,LINESEG l,POINT &np)
&&&&double&r=relation(p,l);
&&&&if(r&0)
&&&&&&& np=l.s;
&&&&&&&&return&dist(p,l.s);
&&&&if(r&1)
&&&&&&& np=l.e;
&&&&&&&&return&dist(p,l.e);
&&& np=perpendicular(p,l);
&&&&return&dist(p,np);
// 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别
double&ptoldist(POINT p,LINESEG l)
&&&&return&abs(multiply(p,l.e,l.s))/dist(l.s,l.e);
/* 计算点到折线集的最近距离,并返回最近点.
注意：调用的是ptolineseg()函数 */
double&ptopointset(int&vcount,POINT pointset[],POINT p,POINT &q)
&&&&int&i;
&&&&double&cd=double(INF),
&&& LINESEG
&&& POINT tq,
&&&&for(i=0;i&vcount-1;i++)
&&&&&&& l.s=pointset[i];
&&&&&&& l.e=pointset[i+1];
&&&&&&& td=ptolinesegdist(p,l,tq);
&&&&&&&&if(td&cd)
&&&&&&&&&&& cd=
&&&&&&&&&&& cq=
&&&&return&
/* 判断圆是否在多边形内.ptolineseg()函数的应用2 */
bool&CircleInsidePolygon(int&vcount,POINT center,double&radius,POINT polygon[])
&&&&double&d;
&&& q.x=0;
&&& q.y=0;
&&& d=ptopointset(vcount,polygon,center,q);
&&&&if(d&radius||fabs(d-radius)&EP)
&&&&&&&&return&true;
&&&&&&&&return&false;
/* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意：如果想从余弦求夹角的话，注意反余弦函数的定义域是从 0到pi */
double&cosine(LINESEG l1,LINESEG l2)
&&&&return&(((l1.e.x-l1.s.x)*(l2.e.x-l2.s.x) +
&&&&&&& (l1.e.y-l1.s.y)*(l2.e.y-l2.s.y))/(dist(l1.e,l1.s)*dist(l2.e,l2.s))) );
// 返回线段l1与l2之间的夹角 单位：弧度 范围(-pi，pi)
double&lsangle(LINESEG l1,LINESEG l2)
&&& POINT o,s,e;
&&& o.x=o.y=0;
&&& s.x=l1.e.x-l1.s.x;
&&& s.y=l1.e.y-l1.s.y;
&&& e.x=l2.e.x-l2.s.x;
&&& e.y=l2.e.y-l2.s.y;
&&&&return&angle(o,s,e);
// 如果线段u和v相交(包括相交在端点处)时，返回true
bool&intersect(LINESEG u,LINESEG v)
&&&&return( (max(u.s.x,u.e.x)&=min(v.s.x,v.e.x))&&&//排斥实验
&&&&&&& (max(v.s.x,v.e.x)&=min(u.s.x,u.e.x))&&
&&&&&&& (max(u.s.y,u.e.y)&=min(v.s.y,v.e.y))&&
&&&&&&& (max(v.s.y,v.e.y)&=min(u.s.y,u.e.y))&&
&&&&&&& (multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)&=0)&&&//跨立实验
&&&&&&& (multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)&=0));
// (线段u和v相交)&&(交点不是双方的端点) 时返回true
bool&intersect_A(LINESEG u,LINESEG v)
&&&&return((intersect(u,v))&&
&&&&&&& (!online(u,v.s))&&
&&&&&&& (!online(u,v.e))&&
&&&&&&& (!online(v,u.e))&&
&&&&&&& (!online(v,u.s)));
// 线段v所在直线与线段u相交时返回true；方法：判断线段u是否跨立线段v
bool&intersect_l(LINESEG u,LINESEG v)
&&&&return&multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)&=0;
// 根据已知两点坐标，求过这两点的直线解析方程： a*x+b*y+c = 0 (a &= 0)
LINE makeline(POINT p1,POINT p2)
&&&&int&sign = 1;
&&& tl.a=p2.y-p1.y;
&&&&if(tl.a&0)
&&&&&&& sign = -1;
&&&&&&& tl.a=sign*tl.a;
&&& tl.b=sign*(p1.x-p2.x);
&&& tl.c=sign*(p1.y*p2.x-p1.x*p2.y);
&&&&return&
// 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200
double&slope(LINE l)
&&&&if(abs(l.a) & 1e-20)return&0;
&&&&if(abs(l.b) & 1e-20)return&INF;
&&&&return&-(l.a/l.b);
// 返回直线的倾斜角alpha ( 0 - pi)
double&alpha(LINE l)
&&&&if(abs(l.a)& EP)return&0;
&&&&if(abs(l.b)& EP)return&PI/2;
&&&&double&k=slope(l);
&&&&if(k&0)
&&&&&&&&return&atan(k);
&&&&&&&&return&PI+atan(k);
// 求点p关于直线l的对称点
POINT symmetry(LINE l,POINT p)
&&& tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b);
&&& tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b);
&&&&return&
// 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交，返回true，且返回交点p
bool&lineintersect(LINE l1,LINE l2,POINT &p)&// 是 L1，L2
&&&&double&d=l1.a*l2.b-l2.a*l1.b;
&&&&if(abs(d)&EP)&// 不相交
&&&&&&&&return&false;
&&& p.x = (l2.c*l1.b-l1.c*l2.b)/d;
&&& p.y = (l2.a*l1.c-l1.a*l2.c)/d;
&&&&return&true;
// 如果线段l1和l2相交，返回true且交点由(inter)返回，否则返回false
bool&intersection(LINESEG l1,LINESEG l2,POINT &inter)
&&& LINE ll1,ll2;
&&& ll1=makeline(l1.s,l1.e);
&&& ll2=makeline(l2.s,l2.e);
&&&&if(lineintersect(ll1,ll2,inter))
&&&&&&&&return&online(l1,inter);
&&&&&&&&return&false;
阅读(...) 评论()
