注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

简约男人

简约,不能简单

 
 
 

日志

 
 
关于我

一个过分渴望被理解的人其实就是一个软弱的人, 勇往直前的力量来自斩钉截铁的决心,不是来自别人的理解.

网易考拉推荐

多种方格网DEM插值算法分析与实现  

2009-10-03 10:28:17|  分类: 算法学习 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

 

说明:为了从抽样空间曲面数据(DEM、TIN)中导出更多信息,需要进行曲面的拟合。此回提交的代码用c++和GCC编译器实现了用于规则方格网的双线性多项式、双三次多项式曲面插值算法。算法的实现中优化掉了许多乘法,所以算法十分有效率。由于使用的是c99标准,用vc6.0编译时需要做一些适应性设置。

 DEM描述的是地面高程信息,它在测绘、水文、气象、地貌、地质、土壤、工程建设、通讯、气象、军事等国民经济和国防建设以及人文和自然科学领域有着广泛的应用。

如在工程建设上,可用于如土方量计算、通视分析等;在防洪减灾方面,是进行水文分析如汇水区分析、水系网络分析、降雨分析、蓄洪计算、淹没分析等的基础;DEM 在无线通讯上,可用于蜂窝电话的基站分析等等;DEM在军事上,可用于确定观察哨所应建立的位置、导航武器的瞄准等。

从DEM绘制透视立体图,能够直观的反映地形的形态。比用等高线表示的地形形态更加接近人的视觉感受,应用起来也更加方便。应用中经常关注的是DEM的某一局部,此时就需要对透视图进行局部放大。无论是规则格网DEM还是不规则三角网(TIN)都只是对空间曲面的一种抽样描述。为了保证放大区域的连续性(是曲面而不是离散的点)需要获得更加细致的高程信息。此时就需要用某种方式对放大区域进行曲面拟合。

此外,还有其它许多需要知道DEM中没有直接给出高程的应用情形。为了从抽样空间曲面数据(DEM)中导出更多信息,需要进行曲面的拟合。曲面的拟合有插值和逼近两种。此回提交的代码用c++和GCC编译器实现了用于规则方格网的双线性多项式、双三次多项式曲面插值算法。算法的实现中优化掉了许多乘法,所以算法十分有效率。由于使用的是C++的c99标准,最好不用vc6.0编译,编译出的结果可能不能正确的执行。

 

 一、方格网DEM的存储

 方格网DEM的是一种规则的高程模型。DEM中存储了地面上在X方向和Y方向等水平间隔的点的高程数据。

                                      +━━━━━+━━━━━+━━━━━+

                                      ┃                ┃                ┃                ┃

                                      ┃                ┃                ┃                ┃

                                      +━━━━━+━━━━━+━━━━━+

                                      ┃                ┃                ┃                ┃

                                      ┃                ┃                ┃                ┃  

                                       +━━━━━+━━━━━+━━━━━+

就像上图,+号表示给出了高程的点,水平和竖直方向上+号的间隔分别相等。如果要求DEM都表示矩形的地面范围,并且在DEM文件的头部指定了左上角点的地理坐标和水平与竖直方向上相邻给出坐标的点之间的间隔。那么就可以不必存储每一点的地面坐标,而只用一个二维数组来存储所有点的高程。需要地理位置数据时可以用某点在数组中的位置和DEM文件头中给出的左上角点的地理坐标和水平与竖直方向上相邻给出坐标的点之间的间隔计算得出。而某点在数组中的位置就是某点的文件坐标。由于文件坐标和地理坐标是一一对应的所以插值中只使用文件坐标,这样可以减少计算量加快计算效率。

如果把高程看作灰度信息,那么DEM文件不就和位图文件差不多了吗。仅有的区别就是文件头部指定了不同的信息,数据区存储的数值不是整数且范围不在0-255之间。ERDAS等许多RS、GIS和制图软件都采用了这种DEM数据存储方式。这些软件中的DEM可以显示为平面的图像和透视图。在显示中需要0-225之间的灰度值,是将高程数据减去某个数再除以某个数,然后取整使之分布在0-255之间。可就像图像处理中灰度拉伸的逆过程。

由于不太清楚各个软件的具体DEM数据格式,所以演示程序是将ERDAS中给出的一个示例DEM文件导出为BMP格式,然后对BMP图像存储的灰度值进行插值运算。

从上边的分析中可以知道DEM的关键数据存储为二维数组,这和灰度图像中的灰度数据存储方式是一样的。所以可以用于DEM插值的算法也可以用于位图放大的插值中。DEM的插值也可以用位图的灰度插值演示。

所以本算法处理的数据就是这种表示高程信息的二维数组,演示程序处理的是表示灰度的二维数组——BMP图像的像素灰度信息。 

 

二、DEM曲面插值算法分析 

1、双线性曲面插值算法

双线性曲面插值的模型为:

Z = f(x,y) = a0 + a1×x + a2×y + a3×x×y;

此模型的四个系数,可以通过方格网的4个顶点的值来确定。为了简化计算,为每一个格网指定一个局部坐标系,如下图:

                                   ∧                                       
                                          ┃                                       
                                          ┃ P1(1,0)                              
                                          ┣━━━━━━━━━━━━┓ P2(1,1)      
                                          ┃                                          ┃              
                                          ┃                                          ┃               
                                          ┃                                          ┃               
                                          ┃                   .P(xp,yp)         ┃               
                                          ┃                                          ┃       
                                          ┃P4(0,0)                            ┃ P3(0,1)      
                                          ╋━━━━━━━━━━━━┻━━━━━━━>

方格网的边长是一个定值,假设它是一个单位长度是可以的。将其带入插值模型解得:

a0 = Z4; a1 = Z1 – Z4; a2 = Z3 – Z4; a3 = Z2 + Z4 – Z1 – Z3;

代入原模型,合并同类项得:

Zp = Z4(1-xp)(1-yp) + Z3(1-yp)(xp) + Z2(xp)(yp) + Z1(1-xp)(yp);

这不就是加权平均值吗。还可以看做卷积运算,如果把四个定点的高程组成一个2×2的矩阵,把上式中的系数看做权矩阵的话。(进一步思考,这难道不和图像几何变换中的重采样的模型一样吗,完全一样,所以这个算法对图像插值的结果就是图像的放大。)

用这种方式看待问题可以简化问题,进一步就是简化算法的复杂度。下面将要说的双三次曲面插值可以用同样的方式去思考,不过,看做加权平均的话,参与加权平均的值有16个;看做卷积的话,矩阵大小是4×4的。

对每一个需要插值的方格建立一个相同的局部坐标系,并进行归一化,即把方格的边长看作一个单位长度。假设,确定了要在每一个格网中在水平和竖直方向上要均匀插入点的数目X,Y;那么每一个格网都有相同的X×Y个4×4的权阵。既然相同就可以把权矩阵的计算提到循环外部。算法的效率又提高了好几倍。

 

2、双三次曲面插值算法

双三次曲面插值的模型:

多种方格网DEM插值算法分析与实现 - tr0217 - 简约男人

该模型需要知道16个系数,将需插值点周围4×4个点带入可以解出各个系数。为了求得系数需要解16个线性方程组成的方程组,每个方程组有16项。这个方程组可以用一些算法像高斯消元法、追赶法、迭代法等算法求解。这些算法不在此论述。

同样我们将为每一个需要插值的格网建立一个局部坐标并进行归一化。带入双三次的模型可以解得一个可以看作加权平均的式子。为了提高算法的效率,下面给出归一化后的近似求权值的方程。

P(x) = 1 – 2×x^2 + |x|^3               (0 <= |x| <= 1)

P(x) = 4 – 8×|x| + 5×x^2 - |x|^3      (1 <= |x| <= 2)

如下图:设p点的坐标为(xp,yp)

                                         ∧                                 
                                                ┃                                 
                                            y4┣━━━━┳━━━━┳━━━━┓              
                                                ┃              ┃              ┃              ┃
                                                ┃              ┃              ┃              ┃              
                                            y3┣━━━━╋━━━━╋━━━━┫              
                                                ┃              ┃      .p     ┃              ┃
                                                ┃              ┃              ┃              ┃              
                                            y2┣━━━━╋━━━━╋━━━━┫               
                                                ┃              ┃              ┃              ┃
                                                ┃              ┃              ┃              ┃             
                                            y1╋━━━━┻━━━━┻━━━━┻━━━━━━> 

                                         x1           x2         x3          x4

上述求权值的方程(函数)的原点定义于待插值的点上。归一化的格网的边长看作1,那么x的变化范围就在0到2之间。上述求权值的方程只是一维的,很显然不合要求。观察双三次插值模型,显然f(x,y)具有很好的对称性,那么在另一维上的近似求权方程也是这个。从而p点的高程值就是:(设P点周围的格网点高程组成的矩阵为I)

Zp = [P(xp - x1)  P(xp – x2)  P(x3 - xp) P(x4 - xp)]*I*[P(yp - y1)  P(yp – y2)  P(y3 - yp) P(y4 - yp)]^T(表示转置)

同双线性插值一样,在对每一个格网进行插值的循环外边计算出各个权阵。

另外一个问题。双三次曲面插值需要周围16个点的高程值,那么最外边一圈方格就不能用双三次曲面插值的方法。在下面的代码中最外边一圈的插值是使用双线性插值来完成的。然而双线性插值的结果和双三次插值的结果区别太大。实际做全局插值时应该用牛顿插值或者拉格朗日插值。不过在DEM的具体应用中出现这个问题的概率较小,一般这个DEM文件的边界格网是另一个DEM文件的内部格网,要使用边界数据时,可以换到另一个DEM文件中去。

 

三、DEM曲面插值算法源码

 

具体源码不在此列出。请到文件Interpolate.h中查阅。此文及算法源码,C语言bmp操作接口的下载地址:http://download.csdn.net/source/1712666

 

四、算法处理结果

下面的两幅图片,分别是对下图所示的DEM模型进行双线性(图一)和双三次(图二)插值的结果。(每个格网边又插入了8个值)

                                                        多种方格网DEM插值算法分析与实现 - tr0217 - 简约男人

从图中可以明显看出,双线性插值的结果要粗糙,有马赛克效果。双线性插值在理论上相邻插值单元的边界值是连续的,但是实际插值的过程中,不论插值的结果采样间隔多么精细,总是离散的。又在每一个插值单元的边界只是连续、不光滑;在边界两边的变化率相差较大,所以就出现了这中情况。

双三次插值的结果明显要细腻的多,没有马赛克效果。双三次插值在理论上相邻插值单元的边界值不仅是连续的,还是光滑的。虽然实际插值的结果总是离散的,但是离散的情况也能够反应连续时的性质。所以双三次插值的结果比双线性插值的结果的效果好。

但是双三次插值的时间复杂度理论上是双线性插值的四倍。但是对算法作出足够的优化,再考虑到现在的硬件执行效率的优势。双三次插值是能够使用在很多以前从时间方面考虑而采用了双线性插值的应用情形。

多种方格网DEM插值算法分析与实现 - tr0217 - 简约男人
                                                                  图(一)

多种方格网DEM插值算法分析与实现 - tr0217 - 简约男人

                                                        图(二)

  评论这张
 
阅读(2419)| 评论(2)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017