等高线光滑与过滤.docx

等高线的光滑与过滤王海超郑先东文章摘要:在地形图绘制中我们需要对等高线进行光滑处理与过滤处理南方公司的Cass 软件有拟合与过滤功能但不能达到我们所要的效果，因此作者想到了另一种处理算法，下面 我们来阐述这种算法文章关键词:光滑，过滤等高线的光滑：光滑等高线有两个要求，第一是要保持等高线的原来位置不变，第二是不能 有尖角这是两个相互矛盾的要求，要满足这两个要求我们需要分两步走：第一步：例如图1所示，我们还从等高线的当前点出发，设当前点为n,第二个点n+1, 第三点n + 2,设n点坐标X0,Y0;设n + 1的坐标为X1,Y1;设n+2的坐标为X2,Y2;要把 n+1点处的尖角去掉，取a边和b边中距离短的一边的四分之一(Dis)，用这个距离作为参 数求出n+1点附近处的两个点n+3,n+4设n+3的坐标为X3,Y3;设n+4的坐标为 X4,Y4.求n+3,n+4点坐标的公式如下：X3=X0+(X1-X0)*dis/a Y3=Y0+(Y1-Y0)*dis/a X4=X2+(X1-X2)*dis/b Y4=Y2+(Y1-Y2)*dis/b图1 图2说明一下：图2为去除尖角后的图，上面的“1/4”是个参数，作者在几次试验后发现用这 个参数比较合理，既可以去掉尖角又不致超过限差的改变等高线的位置。

在使用过程中我们 可以多次使用此步，让等高线变得比较光滑第二步：使用3次B样拟合公式对等高线进行拟合处理，设等高线上的当前点为n，拟 合公式为：f0 = 1 / 6 * (-T A 3 + 3*TA2-3*T + 1)fl = 1 / 6*(3*TA3-6*TA2 + 4)f2 = 1 / 6 * (-3 * T A 3 + 3*T A 2 + 3*T + 1)f3 = 1 /6*T A3x = f0 * Px(n) + f1 * Px(n + 1) + f2 * Px(n + 2) + f3 * Px(n + 3)y = f0 * Py(n) + f1 * Py(n + 1) + f2 * Py(n + 2) + f3 * Py(n + 3)T为小于等于1的参数，计算公式为:i/m,其中”i”为从”0”到”m”的变量，”m”为要添加的节点 数目Px为横坐标点数组，Py为纵坐标点数组x,y为光滑后的等高线上的节点坐标上面的算法已用程序实现，以下为主要部分源代码：Sub CreateFittingLine()Dim Sset As AcadSelectionSetDim filterType(0) As IntegerDim filterData(0) As VariantfilterType(0) = 0filterData(0) = "*line”On Error Resume NextIf Not IsNull(ThisDrawing.SelectionSets.Item("Example")) ThenSet Sset = ThisDrawing.SelectionSets.Item("Example")Sset.DeleteEnd IfSet Sset = ThisDrawing.SelectionSets.Add("Example")Sset.SelectOnScreen filterType, filterDataDim Obj As AcadEntityDim k As IntegerDim a0, a1, a2, a3 As Integer '几个在闭合曲线中使用的参数For k = 0 To Sset.count - 1Set Obj = Sset.Item(k)'3次B样拟合算法Dim Px() As DoubleDim Py() As DoubleDim i, j, m, n As IntegerDim T As DoubleDim f0, f1, f2, f3 As DoubleDim x, y As DoubleDim gc As Long '等高线的高程m = GetVfertexCount(Obj) - 1Dim Zb1() As DoubleDim Zb2() As DoubleReDim Zb1(2 * m + 1) As DoubleFor i = 0 To mZb1(2 * i) = Obj.Coordinate(i)(0)Zb1(2 * i + 1) = Obj.Coordinate(i)(1)Next iCall changeOneRoleToFour(Zb1, Zb2)'去除距离太近的点'可以通过改变此函数中的距离参数来提高线的圆滑，可是有可能照成点线不一。

’把尖角变成4个角m = (UBound(Zb2) + 1) / 2 - 1 '现在有 m+1 个点ReDim Px(m) As DoubleReDim Py(m) As DoubleFor i = 0 To mPx(i) = Zb2(2 * i)Py(i) = Zb2(2 * i + 1)Next iDim Zb() As Double'MsgBox Px(m)n = -1If Obj.Closed = False Then '不闭合的曲线与闭合的曲线处理起来是不一样的，先处理不 闭合的曲线For j = 0 To m - 2If j = 0 Then '把起点带上a0 = jal = ja2 = ja3 = jElself j > m - 3 Then '把终点带上a0 = mal = ma2 = ma3 = mElsea0 = jal = j + 1a2 = j + 2a3 = j + 3End IfFor i = 0 To 3 '"3是个参数，参数越大等高线上节点越多”T = i / 3f0 = 1 / 6 * (-T A 3 + 3 * T A 2 - 3 * T + 1)fl = 1 / 6 * (3 * T a 3 - 6 * T a 2 + 4)f2 = 1 / 6 * (-3 * T a 3 + 3 * T a 2 + 3 * T + 1)f3 = 1 / 6 * T a 3x = f0 * Px(a0) + f1 * Px(a1) + f2 * Px(a2) + f3 * Px(a3)y = f0 * Py(a0) + f1 * Py(a1) + f2 * Py(a2) + f3 * Py(a3)n = n + 1ReDim Preserve Zb(2 * n + 1) As DoubleZb(2 * n) = xZb(2 * n + 1) = yNext iNext jElself Obj.Closed = True Then 处理闭合曲线For j = 0 To mSelect Case jCase Is < m - 2a0 = ja1 = j + 1a2 = j + 2a3 = j + 3Case m - 2a0 = m - 2al = m - 1a2 = ma3 = 0Case m - 1a0 = m - 1al = ma2 = 0a3 = 1Case ma0 = ma1 = 0a2 = 1a3 = 2End SelecdFor i = 0 To 3 '"3是个参数，参数越大等高线上节点越多”T = i / 3f0 = 1 / 6 * (-T A 3 + 3 * T A 2 - 3 * T + 1)f1 = 1 / 6 * (3 * T a 3 - 6 * T a 2 + 4)f2 = 1 / 6 * (-3 * T a 3 + 3 * T a 2 + 3 * T + 1)f3 = 1 / 6 * T a 3x = f0 * Px(a0) + f1 * Px(a1) + f2 * Px(a2) + f3 * Px(a3)y = f0 * Py(a0) + f1 * Py(a1) + f2 * Py(a2) + f3 * Py(a3)n = n + 1ReDim Preserve Zb(2 * n + 1) As DoubleZb(2 * n) = xZb(2 * n + 1) = yNextNext jEnd If'分闭合曲线与不闭合曲线'gc = obj.ElevationDim Pts() As DoubleCall ThreeFiltrates(Zb, Pts)'对生成的坐标过滤Call AddNewdgx(Pts, Obj)Obj.DeleteNext kEnd Sub图3是用上面的例子用此程序处理后的效果图。

等线的过滤：CASS中也有等高线过滤的功能，即“复合线滤波”，但此功能并不能很好的 完成既过滤掉等高线上多余的节点又保持等线的原来位置，例如当滤波值设小了过滤不完等 高线上多余的节点；大了又改变了等高线的位置，如图1，图2所示图4 图5图4为过滤前的等高线，图6是用“复合线滤波”过滤后的等高线于是作者想到了另一个算法即能过滤掉等高线上节点又能保持等高线位置如图6所示图6我们从等高线的起点开始，设当前点为”n”，第二个为” n+1”，第三个为” n+2”，我们把三个点组成 一个三角形，算出”n+1”点到底边的垂直距离”h”我们可以设置一个参数，用于决定是否保 留”n+1 ”点，当h大于这参数时去掉，小于时保留计算”h”的公式：h = 2 * area / c area = -i-s (s - a)(s - b)(s - c)其中”area”为上面三角形的面积,”s”为三角形的周长，”a”为点”n”到”n+1”的边，”b”为”n+1” 到”n+2”的边，”c”为点”n+2”到”n”的边此算法的好处在于可以保证等高线的位置不变又可以充分的过滤掉多余的节点，例如 我们把参数设为“0.2”，等高线部分位置挪动0.2m，对等高线来说没有影响，但是却可以过 滤掉绝大过多数多余的节点。

这个算法已经用代码实现，下面是部分源码：Public Sub FilterDgxCoordinates(ByVal Coords As Variant, ByRef Pts() As Double)Dim i, j, m, n As Integerm = (UBound(Coords) + 1) / 2 - 1 '等高线上节点个数减一，因为是从0开始循环n = 2 * m + 1Dim Zb() As DoubleReDim Zb(n) As DoubleFor i = 0 To n先把值传递过来Zb(i) = Coords(i)Next iIf m < 2 ThenExit SubEnd If'For i = 0 To n'Debug.Print Zb(i)'Next iDim a, b, c As Double '边Dim Area As Double 面积Dim P As Doubl e '周长Dim Cx As Double 垂 线Dim Pt(1) As DoubleDim Pt1(1) As DoubleDim Pt2(1) As Double'i = -1 '先初始化i'MsgBox "ok"'MsgBox m - 2For i = 0 To。