散乱数据的可视化完整.ppt

,第五章 散乱数据的可视化,散乱数据指的是在二维平面上或三维空间中,无规则的、随机分布的数据。 散乱数据的可视化是对散乱数据进行插值或拟合,形成曲线或曲面并用图形或图象表示出来的技术。 散乱数据主要来源于3个方面: 一是物理量的测量数据; 二是科学实验所得数据; 三是科学计算或工程计算的结果数据 散乱数据的可视化有着广泛的应用领域。例如,地质勘探数据、测井数据、油藏数据、气象数据以及有限元计算结果中非结构化数据的显示等,散乱数据的分类按其复杂程度可分为单自变量、双自变量及多自变量。 其可视化的方法又可以分为散乱数据的插值及拟合。 设在二维平面上有 个点 ，并有 ， 插值问题就是要构造一个具有 连续的函数 ,使其在 点的函数值为 ，即 。 本章将主要讨论双自变量散乱数据的插值问题，首先介绍几种双自变量散乱数据的插值方法。然后再讨论大规模散乱数据的插值问题。,5.1 中、小规模散乱数据的插值,5.1.1 与距离成反比的加权法 这一方法首先是由气象学及地质学工作者提出来的,后来由于D.Shepard 的工作被称为Shepard方法。其基本思想是将插值函数F(x,y)定义为各数据点函数值fk的加权平均,即 式中 表示由(x,y)点到(xk,yk)点的距离。值一般取为2。 还可将(5.1)式重写为,式中 这是一种与距离成反比的加权方法，点 的值 对 的 影响与 至 的距离成反比。 （5.3） 式中的权函数有如下性质： （1） ，非负值。 （2） ，是 连续的。 （3） ，当 时， 否则 。 （4） ，具有加权性质。,图5.2,图5.3,图5.4,1 2 3 4 5,2 2 2 2 2,10 1 5 2 3,20 20 20 20 20,假设 点处的 值为 ，如果对（5.1）式求导，则有如下结论 1）如果 ，则 处不存在一阶偏导数。即在该点处形成角点或尖点。 2）如果 ，则 处的一阶偏导数为零。即在该点处的切平面平行于 平面，形成了“平台”效应。 图5.1给出了在单自变量时不同 值得插值结果。图5.2-图5.4则给出了双自变量时不同 值的插值结果。各个图中的 值如表5.1所示。,从以上讨论可以看出,Shepard方法的插值结果只能是C0连续。而且,当增加、删除或改变一个点时,权函数Wk(x,y)均需重新计算,因而该方法是一个全局插值算法。,为了克服Shepard方法的上述缺陷,Franke及Nielson提出了MQS(Modified Quadratic Shepard，改进的二次方程式Shepard)方法,它仍然是一个与距离成反比的加权方法。它的改进主要在以下两个方面: 第1,将(5.1)式中的权函数dk作适当修改,使其只能在局部范围内起作用,以改变Shepard方法的全局插值性质。这时的权函数定义为 式中,rw为一个常数。而,因此，当 点与某一点的距离大于rw时，权值就为零。 第2，用节点函数 代替(5.1)式中的 , 是一个插值于 点的二次多项式，即有 。而且 在点 附近与函数值 具有局部近似的性质。因此，如果认为距离 较远的点对 影响不大，则可以认为在 点附近， 就可近似地表示 了。 根据上述性质， 可由下式表示： 式中 是按最小二乘法由下式得出的优化解：,式中, ， 分别为 及 点的函数值，而 可按下式选取 其中rq为一常数，而 这说明，在（5.4）中，当 时，点 的值就不起作用了。于是 只与相邻点的值有关，因而是函数值得局部近似。 在求出 后，插值函数可表示为,由于MQS方法消除了Shepard方法中的一些缺陷，因而在散乱点插值中得到广泛的应用。但是，为了求得 ，需要多次求解线性方程组，计算量大，因此，一般只用于中、小规模散乱点的插值运算。 5.1.2 径向基函数插值法 径向基函数的来源：即基函数是由单个变量的函数构成的。一个点(x,y)的这种基函数的形式往往是hk(x,y)=h(dk)，这里的dk表示由点(x,y)至第k个数据点的距离。一般来说，这种方法不具有多项式精度，但只要稍加改进，即可获得具有多项式精度的插值公式 式中，qk(x,y)是一个多项式基，其阶次小于m。 上式中的系数ak,bk应满足下面的联立方程组,（5.8）式中的n个方程式满足了插值要求，而（5.9）式中的m个方程式则保证了多项式精度。二式中共有m+n个未知数，同时存在m+n个方程式，联立求解，即可得出待定系数。 下面介绍两种主要的径向基函数插值法。 1.Multiquadric(MQ)方法 Multiquadric 方法是由R.L.Hardy在1971年提出来的。它是最早提出并且应用最为成功的一种径向基函数插值法。它采用的插值函数为 以单自变量简单情况为例，如图5.5所示，设对应于x=xi的点，有值f(xi)。那么这种单自变量情况下的插值函数为,从理论上说，n可以与给定点的数目不等，ai也可以有任意值，但为了求解方便，令n等于给定点的数目，Xj则等于各给定点的x坐标值。此时，如将各Xi点的值fi带入上式，则可得一组联立方程 求解此组联立方程，即可得出 aj。图5.5中表示的是n=6时的结果。图中，不仅表示出F(x)，而且用虚线表示出F(x)的各分量。值得注意的是，构成各分量的绝对值函数均在各Xj处有斜率变化。 如果用双曲线函数代替绝对值函数，那么，同样的点也可以用连续可微曲线来进行插值，这时，插值函数为,式中 表示任意常数。图5.6表示出 时的插值结果。大量的计算表明，如果相对于数据点之间的距离， 的最值过大，则会导致病态的系数矩阵，影响了求解。在大多数情况下， 取值小些则结果较好。,以上的插值公式很容易推广到双自变量的情况，此时有 如果将n个点(xk,yk)的值fk代入上式，即可得到一组联立方程 求解此联立方程，即可求出系数aj(j=1,2,n)的值。 如令Qj表示任意的二次基函数，aj表示系数，fj表示给定点的值。则可将(5.15)式改写为 其矩阵形式为 其解为,2.薄板样条法 薄板样条法就是用此方法求出的散乱点的插值函数使下面这一泛函表达式具有最小值 I(F)表示受限于插值点的无限弹性薄板的弯曲能量。因此，这一方法的实质从力学观点看是使插值函数所代表的弹性薄板受限于插值点，并且具有最小的弯曲能量。这是一个泛函求极值的问题。这一变分问题的解即为我们所需要插值函数，具有式(5.8)和(5.9)的形式，其基函数具有hk(x,y)=h(dk)=dk2logdk 的形式。,那么，插值后，任意一点p的值为 （5.17）,5.1.3有限元方法 由于其基本原理与求解偏微分方程的有限元方法相同，因而被称为散乱点插值的有限元方法。在给出具有双自变量的散乱点和函数值，然后先求出二维平面上散乱点的凸包，并进行三角剖分，形成一系列的三角形，然后构成一系列的面片，使其插值于都有点和对应的函数值。 本节将主要讨论在生成二维点集的三角剖分之后，如何构造插值于各点函数的面片问题。一种最简单的方法是，构造出插值于各点函数的平面三角面片，这在各个三角面片之间只能是C0连续，不能满足要求。因而需要讨论如何利用高阶多项式进行插值，且在各面片具有C1连续的插值方法。 1.Clough-Tocher(克拉夫-托赫尔)插值法 如果用一个三次曲面片对一个三角形的顶点进行插值，该三次曲面片可表示为一个双自变量的三次多项式，即,该式中共有10个待定的系数。 如图5.8(a)所示，设A,B,C是三角形的3个顶点，将3个顶点的坐标值代入上式，可得出3个关系式。如果可以得出3个顶点处相对于x和y的偏导数，并将其代入上式，则又可得到6个关系式。为了使相邻的曲面片光滑连接，那么共享一条边的二个面片，在与该边正交的方向上的偏导数应该相同，这样又可以建立3个关系式。于是我们有了12个关系式，但是只有10个未知数，称为超定方程，无法求解。,2.Herron插值法 G.Herron提出了一种在理论上有创新的方法，他所定义的三角面片插值函数具有如下的性质：,3.偏导数的估计 上面两种插值方法都需要知道在各散乱点处得一阶偏导数。求偏导数的方法有很多，这里仅介绍一种求偏导数近似值得简单方法。,在此基础上 ，构造一个平均平面，该平均平面的各参数值是所有V0ViVj平面的参数的平均值。设平均平面的方程式为,则有,其中,平均平面的偏导数被认为是插值曲面在P0处的偏导数,有了插值曲面在V0处相对于x,y的偏导数之后，求方程导数就不难了。当然，用这个方法求得的结果只是偏导数的近似值，可以直接应用，也可以作为用迭代方法求偏导数精确值时的初始给定值。,4.最小模网格法（MNN法）,由G.M,Nielson提出的最小模网格法是一种比较成功的基于三角面片的散乱点插值方法，有较为广泛的应用。,5.1.4 实例与讨论,前面三节我们介绍了7个不同的算法。现在,用一个例子对其中的4个算法做些比较 。,图5.12(b)中给出了一个测试函数,并随机提取了33个采样点(如图5.12(a)所示)。在图5.12(c）,(d),(f)中,分别给出了用Shepard方法、MNN方法、MQS方法以及MQ方法的插值结果。,从图中可以看出,Shepard方法效果最差。因此,目前这一方法已不再使用。而MNN方法和MQ方法的插值效果最好。而且这两种方法所用的计算时间相差不多。MNN法比MQ法的计算速度稍快一些,但实现起来前者比后者更为复杂。至于MQS方法,其插值质量较Shepard方法虽有所改进,但仍不理想,实现起来相对简单。,前面介绍过的7个不同的算法，可以分为全局插值和局部插值两类。在全局插值算法中，插值函数决定于全部数据点的函数值，当增加、删除或改变一个点的位置或函数值时，插值函数必须重新计算。显然，Shepard方法、MQ方法、MNN方法及薄板样条法等属于这一类方法。在局部插值算法中，一个点的函数值只影响该点附近的一个子域。因而在增加、删除或改变一个点时，只需进行局部的重新计算。前面提到的MQS方法就属于这一类。但是这一优点往往很难实现，因为对于散乱分布的点而言，常常需要将一个点和其他所有的点相比较，以确定某一个点会影响到插值曲面的哪一部分。而且，由全局插值方法得到的插值曲面往往比用局部插值算法的结果更为美观。 在已经介绍过的算法中，插值曲面要通过求解联立方程组来得到，而方程式的数目至少等于散乱点的数目。为了求解的可能性，散乱点的数目必然受到限制。因此，上述算法均只能应用于几百个点的小规模散乱数据的插值。对于散乱点数目达到数千个甚至于数万个的插值问题，则必须另辟途径。,5.2 大规模散乱数据的插值,5.2.1 基于多层B样条的散乱数据插值方法 本节只要讨论基于多层B样条的散乱数据逼近及插值算法。 1.基于算法,图5.16是利用本算法所得的一个例子。图5.16(a)是原始数据。图5.16(b),(c),(d)分别表示m=n=16,m=n=8及m=n=32时的双三次B样条曲面。显然,控制点网格的密度影响着近似函数f的形状。当控制点网格密度较稀疏时,每个控制点的邻近点集覆盖较多的数据点,这就导致较多的点的影响融合在一起,从而使曲面形状比较光滑,但近似的精度却较低。反之,当控制点网格密度较高时,一个数据点的影响范围限制在较小的邻域内,这就使得曲面的近似程度较高,但是在邻近数据点的范围内会出现局部的峰值。由图5.16(c),(d)可以明显地看出这一特性。,2.多层B样条的散乱数据逼近 为了使得用来逼近散乱点的B样条函数既具有比较光滑的形状,又具有较高的近似精度,则需要采用多层B样条函数。假设在矩形域上定义了层次化的控制点网格0 1 .h,并设第一层控制点网格0的间距已经给定,而且后继各层网格的间距均为前一层的一半。那么,如果k是一个(m+3)(n+3)网格,则下一层k+1将有(2m+3)(2n+3)的控制点网格。 当用多层B样条来逼近散乱数据时,首先用上述基本算法在最稀疏的网格0上求得控制点,形成双三次B样条曲面来逼近散乱点集P。这一双三次B样条函数称为f0,它与P中的每一点(xc,yc,zc）存在着距离,设为1zc=zc-f0(xc,yc)。,然后,在较为精细的网格1上求取控制点,以形成双三次B样条函数f1,用来近似差值p1=(xc,yc,1zc)于是,f0+f1与P中的每一点就具有较小的差值,即2zc=zc-f0(xc,yc) -f1(xc,yc)。如此类推,在第k层时,我们可以在网格k上得到控制点,以近似表示pk=(xc,yc,kzc),其中kZc=k-1zc-fk-1(xc,yc)而最终的近似函数可表示为,3.多层B样条的散乱数据插值 前面介绍了用多层B样条逼近散乱数据的算法。那么,现在用多层B样条实现散乱数据的插值。,令d为散乱数据集pk中所有点对之间的最小距离。如果d4,那么,任何两个点在其44的邻域内均无共同的控制点。在这种情况下, k 中的每一个控制点在其邻近点集中最多只包含一个点,因而在应用上述基本算法构造B样条函数时,(5.37)式简化为(5.36)式,从而使B样条函数插值于pk中的每一个点。图5.18表示出对数据点p1,p2,进行插值和逼近的例子。值得注意的是,对数据点作插值或是逼近决定于两点间控制点网格的数目,并非决定于它们之间的实际距离。,5.2.2 自适应的层次B样条散乱点的插值方法 层次B样条散乱点的插值方法的缺陷:当在比较稀疏的网格上求出控制点并生成B样条曲面对给定的散乱点进行逼近时,如果曲面的某些部分误差较大,即便是一小部分的误差不能满足要求,也需要在全局范围内对网格加以细分,并重新计算控制点。从而导致相当一部分的计算是不必要的。如果能使网格的细分及控制点计算仅在误差超过容许值的子区域内进行,势必大大节约计算时间及空间,且不断提高逼近精度。 为了克服这一缺点, 提出了关于大规模散乱数据自适应层次B样条曲面的逼近和插值算法。该算法用文献LEE97中提出的方法计算出逼近散乱数据的B样条控制点网格,并计算出每一个散乱点处的误差。对于误差大于容许值的子区域,局部地构造出更精细的网格,按同一方法计算出控制点,得到误差值较小的逼近。对于误差小于容许值的区域,则不作处理。这一过程将递归地进行,直至满足所要求的精度或误差值为零。,1.几个关键问题 （1）不同层次B样条曲面的C2 连续性 网格的细化是在误差大于容许值的子区域中进行的,因此,必须保证由不同层次的网格所定义的B样条曲面之间具有C2连续性。假设在某一子区域中,某些点的逼近误差大于门限值t。我们将该子区域的范围沿4个方向向外扩展一个网格的距离。这里,所谓的网格距离是指父网格的宽度,因而是子网格间距离的两倍。如图5.19所示,子区域 扩展为 。子区域内控制点网格的密度是父区域的两倍.的边界与 的边界之间的区域称为约束区,这一区域内的控制点网格密度与子区域 内相同(如图5.20所示)。,（2）约束条件下B样条曲面的逼近及插值,于是 ，根据（5.40）式可以由父区域的控制网格矩阵Z,直