Day2计算几何与三维凸包.pptx

专题讨论：计算几何与凸包算法AC_Aerolight2013.4.30 PREFACE这一节的讨论集中于几个比较实用但关注度不高的凸包算法。穿插一些计算几何基础题目今天的目的依然是科普，因此内容会比较简单。已经对相关内容有所了解的同学可以跳过这些，但不要打扰其他同学。Lets begin. 计算几何基础直线表示标准式 Ax+By+C = 0斜率式 y = kx+b两点式 (x-x0)/(x1-x0) = (y-y0)/(y1-y0) x/A+y/B = 1 Balabala 计算几何基础点积二维情形：ab = axbx+ayby = |a|*|b|*cos =arccos(/|a|/|b|) Dis(a,b) = |B-A| = 叉积二维情形：ab = axby-aybx |ab| = |a|*|b|*sin 表示以a,b为边的平行四边形的有向面积。右手定则：用面积判定方向可以用于计算多边形面积 凸包：定义什么是凸多边形 As precise as you can.对于形内两点X和Y，线段XY一定完全位于形内对于多边形的每条边，整个多边形都在其一侧 凸包对于点集P，包含P的最小凸多边形 给定点集P=(Xi,Yi)，求点集P的凸包。 增量法按输入顺序将点加入点集维护当前点集的凸包(1) 新点在凸包内 O(n)判定(2) 新点在凸包外删除从新点可见的凸包边用链表维护凸包 O(n)删除和判断O(n2)在线算法 增量法改进时间复杂度用数据结构分别维护上下凸壳 O(nlogn)三点共线的情况可见性判定 JAVRIS步进法（卷包裹法）假设有一根无限长的绳子，一开始绳子靠在点集的一个外围点上。随后绳子沿某个方向（如逆时针）绕动，直到碰到一个点之后停止绕动以新点为中心继续沿指定方向绕动绳子当回到起点时，得到凸包选取一个一定在凸包上的点X，沿着点集逆时针走一圈，当走回X时，得到整个凸包 (1)X怎么选？ (2)如何找最右手边的射线 JAVRIS步进法（卷包裹法）复杂度分析 O(n*h) H是凸包上的节点数 Output-Sensitive GRAHAMS CONVEX HULL目前使用最为广泛的算法。确定极角序选取最左下的点作为参考点，按夹角对其他点排序按夹角排序？叉积定方向偏序关系 GRAHAMS CONVEX HULL维护凸壳节点按极角序入栈保持栈中节点凸壳性质弹出栈顶元素直到新边和原栈顶边成左转关系复杂度分析一个节点进栈一次出栈一次，O(n) 排序O(nlogn)总复杂度O(nlogn)常用技巧：维护水平序节点 DIVIDE-AND-CONQUER考虑对点集进行分治Solve(S)返回S的凸包。Solve(S)将S分为两个点集S1和S2，保证S1在S2左侧合并Solve(S1)和Solve(S2)返回的凸包 合并凸包的关键：找到切线上下切线可以分别考虑 DIVIDE-AND-CONQUERU-,U,V和U,V,V+都成逆时针顺序排列双指针扫描 从L的最右端和R的最左端开始维护U上可见点的最远点，直到一个点都看不见 QUICKHULL回忆Quicksort选定一个标准mid，将mid的部分分别排序定义过程QuickHull(a,b)保证A-B是凸包上的两个顶点。定义点集S，使S中的点除了a,b以外都在a-b的右侧。 A-B这条边被称为凸包的弦(chord)QuickHull(a,b)应当返回S的凸包。 QUICKHULLQuickHull(a,b)找到离ab最远的点c，那么c一定在凸包上 QuickHull(a,c) QuickHull(c,b)将点集分成三个部分 (1)三角形内点 (2)AC右侧点 (3)CB右侧点 对(2)(3)分治(1)的点全部抛弃(2)和(3)会有交集么？(1)(2) (3) AB C QUICKHULLQuickHull的效率三角形ABC内的点一定不在凸包上在点集随机的情况下，复杂度十分优秀能够被圆周撒点卡到O(nlogn)没有点在三角形内初始调用指定点集的最左/最右点x,y调用QuickHull(x,y)和QuickHull(y,x) 去递归 O(NLOGN) BOUND证明：凸包的时间复杂度不会低于O(nlogn)考虑点集Pi=(ai,ai2)对Pi求凸包，实际上就是对Ai的排序排序的复杂度不会低于O(nlogn)凸包的复杂度不会低于O(nlogn)Not Output-Sensitive CHANS ALGORITHM假设已知凸包上有M个点。Step 1：将点均分为N/M个部分对每个点集进行Graham-Scan，得到N/M个凸包复杂度统计 O(n/m)*(mlogm) = O(nlogm)Step 2：合并N/M个凸包。 CHANS ALGORITHM CHANS ALGORITHMStep 2：合并N/M个凸包卷包裹选取初始点x0，每次选择最右手边的一个点前进可能的点都在n/m个凸包上在每个凸包内部二分查询：O(logm)时间复杂度选取最右点：O(n/m*logm)凸包上有M个点：O(m*n/m*logm)=O(nlogm) Chans Algorithm求凸包的复杂度是O(nlogm)M是啥? CHANS ALGORITHM将求凸包过程记为Solve(N,M)假设最终凸包上的点数为H。从小到大枚举M？ O(nlog1)+O(nlog2)+O(nlogh) =O(n)*(log1+log2+logh) =O(n)*(log(h!)=O(n)*O(hlogh)=O(nhlogh) 只要M=H，算法就能正确出解 CHANS ALGORITHM倍增枚举M？M=2,4,8,16,2t, O(nlog2)+O(nlog4)+O(nlogh) =O(n)*(1+2+3+logh) =O(n)*O(log2h) =O(nlog2h)比指数更快的增长率M=2,4,16,256,2(2t). O(nlog2)+O(nlog4)+O(nlogh) =O(n)*(1+2+4+.+logh) =O(n)*O(logh)=O(nlogh) MINKOWSKI和给定两个点集P,Q，定义他们的Minkowski和P+Q = X1+X2,Y1+Y2 | (x1,y1)P,(x2,y2)Q给定点集P,Q,求他们Minkowski和的凸包面积|P|,|Q| = 105 MINKOWSKI和等价于P和Q凸包的Minkowski和凸包的边构成(Px+Py,P(x+1)+Py) or (Px+Py,Px+P(y+1)相当于用P和Q凸包的所有边做一个新凸包双指针扫描解决 PART II回顾：凸包增量法 O(n2) Javris步进法 O(nh) Graham扫描 O(nlogn) QuickHull O(n2)O(n) Chan分块法 O(nlogh) 将凸包扩展到三维情况。 计算几何基础：表达定义直角坐标系右手坐标系三维空间的点坐标: (x,y,z)点和向量有什么区别Ax+By+Cz+D = 0描述一个面 (A,B,C)是平面的法向量右手法则 描述一条直线基础：两个平面的交点 (A1x+B1y+C1z+D1=0,A2x+B2y+C2z+D2=0) 计算几何基础：表达直线：参数方程 L=X0+tv tR X0是L上的一个点 V是一个向量，表示L的方向直线：两点式 (x-x0)/(x-x1)=(y-y0)/(y-y1)=(z-z0)/(z-z1)令v=P1-P0，则等价于上面的式子平面：点法式/标准式 Ax+By+Cz+D=0 A(x-x0)+B(y-y0)+C(z-z0)=0判定点和平面关系 Ax0+By0+Cz0+D 0 -点与法向量同侧 计算几何基础：运算点积 = AxBx+AyBy+AzBz+= = |a|*|b|*cos = 0 a b也被称为内积(Inner Product) 计算几何基础：运算P1：求一个点Q到直线(Po,V)的距离 P0是直线上一点，v是方向向量直线上的点Q满足Q=P0+v，vRP2：求一个点Q到平面(Po,N)的距离 P0是平面上一点，N是法向量第一个问题：如何判断点X在不在平面上 = 0 或者 = 计算几何基础：运算Solution to P1：这个方法也适用于2D。 计算几何基础：运算Solution to P2：假设q是Q在上的投影。 计算几何基础：运算P3：给定两条异面直线，求它们之间的距离假设两条直线是L1=P1+us,L2=P2+vt (s,tR)垂线交点是Q1=P1+us,Q2=P2+vt那么有 = 0和 = 0 计算几何基础：运算 ARDENIA给定三维空间上的两条线段，求它们的最近距离。分情况讨论(1)线段所在直线是异面/相交直线计算异面直线距离当且仅当两个垂点都在线段上时取到(2)直线平行，或(1)的条件不满足分别计算四个端点到另一条线段的最短距离 视为平面问题 STAR WAR给定三维空间上的两个四面体，求它们的最近距离。保证四面体相离分情况讨论点-面距离线-面距离？面-面距离？ 计算几何基础：运算二维叉积 ab=AxBy-AyBx ab=|a|b|sin三维叉积：定义 c=ab=|a|b|sin*n n是和ab所在平面垂直的单位向量 ab = -ba也被称为外积(Outer Product) 确定n的方向右手定则/右手系右手四指从X扫到Y，大拇指方向为N 计算几何基础：运算平面上三个不共线点确定唯一平面。给定三个点M1,M2,M3，求平面解析式。N = (M2-M1)(M3-M1)平面上的点P应当满足 = 0 计算几何基础：运算计算P1(x1,y1,z1) P2(x2,y2,z2)的值假设i=(1,0,0),j=(0,1,0),k=(0,0,1)P1 P2 = (x1i+y1j+z1k) (x2i+y2j+z2k)=x1x2(i i)+x1y2(i j)+x1z2(i k)+y1x2(ji)+y1y2(j j)+y1z2(j k)+z1x2(k i)+z1y2(k j)+z1z2(k k)=(y1z2-y2z1)i+(z1x2-z2x1)j+(x1y2-x2y1)k 每一项看起来都像一个二维叉积。回忆：三阶行列式按第1行展开 计算几何基础：运算叉积的另一种写法：三阶行列式对上式进行求值，可以得到相同的结果。 222 11121 ZYX ZYX kjiPP 计算几何基础：运算2D：用ab表示以a和b为边的平行四边形面积3D：|ab|表示以a和b为边的平行四边形面积平行四边形加一维就是平行六面体给定向量a,b,c，求以a,b,c为三边的平行六面体体积 计算几何基础：运算求平行六面体的面积S(Base)=|bc|H=|a|cos也是bc和a的夹角V=SH=|a|*|bc|*cos=a (bc) a(bc)被称为向量的混合积，记作a,b,c 计算几何基础：运算为什么叫a,b,c实质是一个三阶行列式a,b,c = |det(a,b,c)|混合积满足轮换性，a,b,c=b,c,a=c,a,ba,b,c = -c,b,a什么时候a,b,c 0? A,B,C成右手系如何求以a,b,c为边的三棱锥体积？ V(三棱锥)=1/6V(六面体)三维空间中的三棱锥也被称为单纯形 A LETTER TO PROGRAMMERS给定一个向量，模拟以下操作(1)平移增量(x,y,z)(2)将向量放大到倍(3)将向量沿指定的旋转轴(x,y,z)逆时针旋转d度(4)将之前的x个操作重复执行y次序列长=1000 涉及数字除角度外均为整数且(d,e,f)，方向(u,v,w)是单位向量将点(x,y,z)带入，可以得到第二式 凸包：定义和实现在三维空间上定义凸多面体(convex polyhedron)对于形内两点X和Y，线段XY一定完全位于形内对于多边形的每个面，整个多边形都在其一侧 如何存储凸多面体？关键：存储凸多边形的面(facet)拆一条边为两条半边见右面的环绕标志 所有平面法线向外 AERODYNAMICS给定凸多面体S上的n个顶点，假设所有顶点的z坐标最小为zmin最大为zmax，求S在z=zmin到z=zmax中每个整数平面上截取的面积。N=100,x,y,z=100且均为整数. AERODYNAMICS需要计算三维凸包么对于zmin,zmax中的每个整数，计算所有边在z=z上的投影，计算凸包面积即可 SHADE OF HALLELUJAH MOUNTAIN给定凸多边形S，一个投影平面P0和光源L求S在P0上造成阴影的体积大小简单起见，保证平面和凸多边形不相交。N=100.平面旋转后求凸包即可不用旋转知识怎么做 卷包裹法(GIFT-WRAPPING)回忆2D下的卷包裹法Expand to 3D一张纸紧贴着凸包的一条边AB沿AB的右手边收拢，直到遇到顶点C为止那么ABC是凸包上的一个面下一个方向？递归：沿ac和cb方向搜索 去递归拓展性：能不能拓展到高维度上 卷包裹法(GIFT-WRAPPING)初始边的选择（见图(a)） 将点集投影到平面上求凸包，凸包的任一条边均可时间复杂度O(nh) 卷包裹法(GIFT-WRAPPING) void wrap(int i, int j) int k = i, l, m; for(m = 0; m h; m+) if (Im = i for(l = 0; l EPS) k = l; if(k = j) return; Ih = i; Jh = j; Kh = k; h+; wrap(k,j); wrap(i,k); 增量法(INCREMENTAL)考虑在线版本的问题给定点集P，每次加入一个点，动态维护P的凸包。假设当前凸包未退化（体积0）(1)点P在凸包内或凸包面上(2)点P在凸包外删除不必要的凸包面将新点加入凸包判定面(a,b,c)对点P可见 0 N=(b-a)(c-a) 0 - p-a,b-a,c-a 0 增量法(INCREMENTAL)白色的面从Pr可见，灰色的面从Pr不可见黑色的边被称为地平线 地平线的两端，一定是一个可见面一个不可见面 增量法(INCREMENTAL)删除可见面用新点连接地平线的每一条边 初始化：得到一个非退化的三棱锥初始化失败的情况 增量法(INCREMENTAL)证明：凸多面体的边数和面数是O(n)级别 (1)由欧拉公式，F(面数)+V(点数)-E(边数)=2 (2)一条边恰归属于两个面，一个面有=3条边，2E=3F联立(1)(2)可得E=3V-6,F 0 QUICKHULLQuickhull(P0,P1,P2)找到离平面(P0,P1,P2)最远的点K删除在三棱锥K-P0-P1-P2内的点和面 Quickhull(K,P0,P1) Quickhull(K,P1,P2) Quickhull(K,P2,P0)红色箭头为各平面法线 这个算法是有效的么？ QUICKHULL问题在哪里？观察点T的位置T同时在平面K-P0-P1和K-P2-P0外侧一个点会被递归多次 为什么2d下没有这个问题？递归的各个集合有交集处理方式：随意分配到一个集合内T QUICKHULL证明：随意分配可以保证算法正确性如果一个点是凸包上的点，那么永远不会被其他三棱锥包含而剔除因此随意分配不会导致漏算等情况Quickhull(P0,P1,P2,S)找到离平面(P0,P1,P2)最远的点K删除在三棱锥K-P0-P1-P2内的点和面将S中剩下的点分到三个外集中 Quickhull(K,P0,P1,S1) Quickhull(K,P1,P2,S2) Quickhull(K,P2,P0,S3) QUICKHULL IN KDQuickhull拓展性很好，可以解决任意维凸包。(D-1)维平面被称为凸包的facet，是凸包的组成部分。相当于二维下的边和三维下的面(D-2)维平面被称为凸包的ridge，是两个facet的交集。相当于二维下的点和三维下的边Simplified Beneath-Beyond Theorem若H是空间Rd上的一个凸包，P是H外一点，一个facet F是H+P的凸包当且仅当： (1) FH，P在F下方 (2) F H，且F由P和一个ridge组成，这个ridge关联的一个facet在p下方，其他在p上方 QUICKHULL IN KD 初始化 寻找最远点 查询地平线 创建facet 分配顶点 删除可见面复杂度：O(nlogn)(d=4) CHANS ALGORITHM回忆二维的Chans AlgorithmStep 1：将点均分为N/M个部分对每个点集进行Graham-Scan，得到N/M个凸包Step 2：合并N/M个凸包卷包裹选取初始点x0，每次选择最右手边的一个点前进在凸包上二分 Chans Algorithm能否扩张到三维？ CHANS ALGORITHMStep 1Graham-Scan不能扩展到三维三维中没有非常良好的序用随机增量/Quickhull代替，依然是O(nlogn)Step 2卷包裹算法不能优化到对每个凸包log m三维中没有良好的点序给定一个多面体，在O(logn)时间内得到卷包裹的下一个目标 给定一个凸点集，在O(logn)时间内，找到点P使得面ABP和面ABC二面角最大 DOBKIN-KIRKPATRICK HIERARCHY对点集进行D-K分层假设点集序列P0,P1,PkP0=P，是原凸多边形Pi+1由Pi删掉一个点集的独立集得到删掉的点中没有连边的Pk是一个单纯形（三棱锥或四面体） DOBKIN-KIRKPATRICK HIERARCHY 以上是一次找删除点集的操作可以证明这么做保证K最大为O(logn)级别 DOBKIN-KIRKPATRICK HIERARCHY查询点集 从Pk开始查询，Pk的查询是O(4)的每次退回Pi-1查询，只查询Pi-1中当前结果的邻居可以证明查询的复杂度也是O(logn) CHANS ALGORITHM解决了点集查询问题，剩下的步骤同以往时间复杂度是O(nlogh) CHANS D-C ALGORITHM基于2维分治法的3维凸包和上一个名字一样，加上词语以区别两个算法的作者是同一个人Timothy M.Chan回忆2维分治法 分别考虑上下凸壳 CHANS D-C ALGORITHM将三维问题转为二维问题考虑凸包的一个facet，所在的平面是z=sx+ty+b将t看成时间参数，定义t时刻的平面点集Pt=(x,z-ty)(1)点P在三维凸包的下凸壳上，当且仅当在某个时间t时，Pt落在直线y=sx+b上，且其他点在该直线上方(2)边P1P2在三维凸包的下凸壳上，当且仅当某个时间t时，P1P2是二维凸壳的一条边 求t=-inf+inf之间凸壳的动画 CHANS D-C ALGORITHM虽然时间t是无穷的，但凸包点集的变化次数有限 (1)点Pj进入凸包，删除边PiPk，加入PiPj/PjPk (2)点Pj离开凸包，删除PiPj/PjPk，加入PiPk一个点进入凸包一次离开凸包一次将点集按x排序，进行分治Solve(l,r) Solve(l,mid) Solve(mid+1,r)合并两段点集的动画 CHANS D-C ALGORITHM关键操作：合并点集L和R的凸壳动画为A初始动画：合并两个t=-inf时的下凸包二维的凸包合并维护可能发生的增删点事件 (1)A和B的增删点操作 (2)切线u-v的改动(1)的维护如果L发生了增删点，当且仅当涉及点在U左侧时，A需要增删对应的点 如果R发生了增删点，当且仅当涉及点在V右侧时，A需要增删对应的点以上两个事件对应的时刻由分治过程得到。 CHANS D-C ALGORITHM(2)的维护U-V是切线的充要条件是(1)u-uv是逆时针，(2)uu+v是顺时针，(3)uvv+是逆时针，(4)uv-v是顺时针。 (1)不满足：新的切线是u-v，从u-和v之间删除u (2)不满足：新的切线是u+v，在u和v之间插入u+ (3)不满足：新的切线是uv+，在u和v+之间删除v (4)不满足：新的切线是uv-，在u和v之间插入v- CHANS D-C ALGORITHM(2)的维护计算事件发生的时间：三点共线复杂度O(nlogn) Implementation CHANS D-C ALGORITHMHull(a,b)：a用于存储事件列表，b用于辅助 Turn(a,b,c)返回三点是顺时针还是逆时针假设此时last和next里存储的是两个t=-inf时的凸包 CHANS D-C ALGORITHM 合并凸壳的事件 CHANS D-C ALGORITHM时光倒流在prev和next中存下点j进入凸包时的邻接点，以及t=-inf时的凸包 用于返回上一层，以及输出答案用 CHANS D-C ALGORITHM输出答案对返回的事件序列进行模拟即可当一个点j被插入凸包时，假设两侧节点是i和k存在时刻t使i,j,k共线，且其他顶点在直线上方 I,j,k都在某个平面z=sx+ty+b上，且其他点在平面上方另一侧凸包可以类似处理 CHANS D-C ALGORITHM给定平面点集P，在线查询离输入点(a,b)最近的点。点距离用直线距离。Dist2=(a-x)2+(b-y)2 = (x2+y2)-2ax-2by+(a2+b2)A2+B2是常数，可以忽略在三维空间上定义点(X,Y)-(2X,2Y,X2+Y2) 存储分治中所有时刻的凸包点集，考虑时刻b的情况 CHANS D-C ALGORITHM在时刻B，(x,z-by) = (2X,X2+Y2-2BY)而此时求的dist等于x2+y2-2ax-2by恰好等于新坐标系中的-ax+y在新凸包上进行二分即可 PART II CONCLUSIONO(n2)卷包裹法 输出敏感，实现简单增量法 思路直观O(nlogn)随机增量 思路直观，效果不错，但较难实现 D-C（分治法） 做法直观，效率也很好，但较难实现 QuickHull 扩展性好，直接扩展到高维 Chans Algorithm 输出敏感，时间复杂度最优，但实现非常困难，需要辅助的定位结构 Chans D-C 思路精巧，实现灵活，但效率较低 重心P1：给定平面上的三角形，求其重心坐标。几何定义：顶点与对边中点的连线交点为重心 (Gx,Gy) = (Ax+Bx+Cx)/3,(Ay+By+Cy)/3) G=(A+B+C)/3P2：给定空间内的三棱锥，求其重心坐标。几何定义：顶点和对面重心的连线交点为重心。 G=(A+B+C+D)/4可以推广到N维单纯形上 P3：给定平面内的凸包，求凸包的重心坐标。P4：给定空间内的凸包，求凸包的重心坐标。 重心Solution to P3对凸包进行三角剖分对于每个三角形，求出其面积Si，重心Gi点上带权重的多边形重心 G=(SiGi)/(Si)以原点进行三角剖分 Gx=(Xi+Xi+1)*(PiPj)/6/SSolution to P4假设无四点共面，凸包每个面都是三角形 重心以原点进行三角剖分G=(ViGi)/(Vi)假设各个面是i-j-k Gi=(Pi+Pj+Pk)/4，Vi=Pi,Pj,Pk/6 G=(Pi+Pj+Pk)*Pi,Pj,Pk)/24V一个面上有多个点的情况对每个面进行三角剖分再带入上式 假设K个面，每个面有Si个节点 G=(Pi,1+Pi,j+Pi,j+1)*Pi,1,Pi,j,Pi,j+1)/24V顺便得到了多面体的体积。 WF2010 PAPERWEIGHT你的公司生产一种压纸器。压纸器本身由两个共面的四面体组成（5个顶点，6个面和9条边），压纸器的某个地方镶嵌了一个RFID芯片。一般用户会将压纸器放在计算机顶部，而RFID芯片需要足够接近计算机才能够正常工作。你的任务是计算在压纸器稳定放置于计算机顶部的前提下，芯片离计算机顶部的最远和最近距离。 称压纸器为稳定的，当且仅当压纸器的重心向任何方向移动0.2单位时，压纸器都处于平衡状态。 WF2010 PAPERWEIGHT Fig10：一个压纸器的示意图Fig11：两种放置方法 WF2010 PAPERWEIGHTC(5,3)枚举底面(1)两个顶点不在底面一侧(2)两个顶点在底面同一侧稳定性测试重心离底面三边的距离不小于0.2Tricky Case 四点共面底面扩大为四边形稳定性测试：重心离底面四边的距离不小于0.2 ASTEROIDS给定空间上的两个凸多面体，允许随意旋转移动，求它们重心间的最短距离Dist(G1,CH1)+Dist(G2,CH2)把两个垂点贴在一起就行了 INTERSECTION OF TWO PRISMS有两个无穷大的棱锥P1和P2，P1的所有侧棱和z轴平行，P2的所有侧棱和y轴平行。根据两个棱锥的性质，可以用P1和xOy的截面来定义P1，用P2和xOz的截面来定义P2。已知P1和P2的截面都是凸多边形。求P1和P2的交的体积。 N,M=100，坐标=100且均为整数。 INTERSECTION OF TWO PRISMS两个棱锥的截面图和实际样式 INTERSECTION OF TWO PRISMS INTERSECTION OF TWO PRISMS求出三维凸包的样式？寻找突破口Q：目标多面体在任意垂直X轴的平面截面是什么A：一定是长方形为什么？ 由于坐标是整数，可以分段求解两个相邻的截面中间是一个四棱台 INTERSECTION OF TWO PRISMS求四棱台的体积S=h*1/3(S1+S2+sqrt(S1S2) S1是下底面面积，S2是上底面面积避免根号用f(t)表示在x=t处截得的面积大小 F(t)=width(S1,t)*width(S2,t)S=1/6(F(l)+F(h)+4F(mid) 易知f(t)是二次函数，S就是该函数的积分在计算近似面积时，也可以利用这种法则