等高線光滑與過濾

1、等高線的光滑與過濾 王海超 鄭先東文章摘要: 在地形圖繪制中我們需要對(duì)等高線進(jìn)行光滑處理與過濾處理。南方公司的Cass軟件有擬合與過濾功能但不能達(dá)到我們所要的效果，因此作者想到了另一種處理算法，下面我們來闡述這種算法。文章關(guān)鍵詞: 光滑，過濾。等高線的光滑：光滑等高線有兩個(gè)要求，第一是要保持等高線的原來位置不變，第二是不能有尖角。這是兩個(gè)相互矛盾的要求，要滿足這兩個(gè)要求我們需要分兩步走： 第一步：例如圖1所示，我們還從等高線的當(dāng)前點(diǎn)出發(fā),設(shè)當(dāng)前點(diǎn)為n,第二個(gè)點(diǎn)n+1,第三點(diǎn)n+2, 設(shè)n點(diǎn)坐標(biāo)X0,Y0;設(shè)n+1的坐標(biāo)為X1,Y1; 設(shè)n+2的坐標(biāo)為X2,Y2;要把n+1點(diǎn)處的尖角去掉，取a邊

2、和b邊中距離短的一邊的四分之一(Dis)，用這個(gè)距離作為參數(shù)求出n+1點(diǎn)附近處的兩個(gè)點(diǎn)n+3,n+4。設(shè)n+3的坐標(biāo)為X3,Y3;設(shè)n+4的坐標(biāo)為X4,Y4.求n+3,n+4點(diǎn)坐標(biāo)的公式如下：X3=X0+(X1-X0)*dis/aY3=Y0+(Y1-Y0)*dis/aX4=X2+(X1-X2)*dis/bY4=Y2+(Y1-Y2)*dis/b 圖1 圖2說明一下：圖2為去除尖角后的圖，上面的“1/4”是個(gè)參數(shù)，作者在幾次試驗(yàn)后發(fā)現(xiàn)用這個(gè)參數(shù)比較合理，既可以去掉尖角又不致超過限差的改變等高線的位置。在使用過程中我們可以多次使用此步，讓等高線變得比較光滑。第二步：使用3次B樣擬合公式對(duì)等高線進(jìn)行擬

3、合處理，設(shè)等高線上的當(dāng)前點(diǎn)為n，擬合公式為： T為小于等于1的參數(shù)，計(jì)算公式為:i/m,其中”i”為從”0”到”m”的變量，”m”為要添加的節(jié)點(diǎn)數(shù)目。Px為橫坐標(biāo)點(diǎn)數(shù)組，Py為縱坐標(biāo)點(diǎn)數(shù)組。x,y為光滑后的等高線上的節(jié)點(diǎn)坐標(biāo)。上面的算法已用程序?qū)崿F(xiàn)，以下為主要部分源代碼：Sub CreateFittingLine() Dim Sset As AcadSelectionSet Dim filterType(0) As Integer Dim filterData(0) As Variant filterType(0) = 0 filterData(0) = *line On Error Resu

4、me Next If Not IsNull(ThisDrawing.SelectionSets.Item(Example) Then Set Sset = ThisDrawing.SelectionSets.Item(Example) Sset.Delete End If Set Sset = ThisDrawing.SelectionSets.Add(Example) Sset.SelectOnScreen filterType, filterData Dim Obj As AcadEntity Dim k As Integer Dim a0, a1, a2, a3 As Integer 幾

5、個(gè)在閉合曲線中使用的參數(shù) For k = 0 To Sset.count - 1 Set Obj = Sset.Item(k) 3次B樣擬合算法 Dim Px() As Double Dim Py() As Double Dim i, j, m, n As Integer Dim T As Double Dim f0, f1, f2, f3 As Double Dim x, y As Double Dim gc As Long 等高線的高程 m = GetVertexCount(Obj) - 1 Dim Zb1() As Double Dim Zb2() As Double ReDim Zb1

6、(2 * m + 1) As Double For i = 0 To m Zb1(2 * i) = Obj.Coordinate(i)(0) Zb1(2 * i + 1) = Obj.Coordinate(i)(1) Next i Call changeOneRoleToFour(Zb1, Zb2) 去除距離太近的點(diǎn)可以通過改變此函數(shù)中的距離參數(shù)來提高線的圓滑，可是有可能照成點(diǎn)線不一。把尖角變成4個(gè)角 m = (UBound(Zb2) + 1) / 2 - 1 現(xiàn)在有m+1個(gè)點(diǎn) ReDim Px(m) As Double ReDim Py(m) As Double For i = 0 To m

7、 Px(i) = Zb2(2 * i) Py(i) = Zb2(2 * i + 1) Next i Dim Zb() As Double MsgBox Px(m) n = -1 If Obj.Closed = False Then 不閉合的曲線與閉合的曲線處理起來是不一樣的，先處理不閉合的曲線。 For j = 0 To m - 2 If j = 0 Then 把起點(diǎn)帶上 a0 = j a1 = j a2 = j a3 = j ElseIf j m - 3 Then 把終點(diǎn)帶上 a0 = m a1 = m a2 = m a3 = m Else a0 = j a1 = j + 1 a2 = j

8、+ 2 a3 = j + 3 End If For i = 0 To 3 3是個(gè)參數(shù)，參數(shù)越大等高線上節(jié)點(diǎn)越多 T = i / 3 f0 = 1 / 6 * (-T 3 + 3 * T 2 - 3 * T + 1) f1 = 1 / 6 * (3 * T 3 - 6 * T 2 + 4) f2 = 1 / 6 * (-3 * T 3 + 3 * T 2 + 3 * T + 1) f3 = 1 / 6 * T 3 x = f0 * Px(a0) + f1 * Px(a1) + f2 * Px(a2) + f3 * Px(a3) y = f0 * Py(a0) + f1 * Py(a1) + f2

9、 * Py(a2) + f3 * Py(a3) n = n + 1 ReDim Preserve Zb(2 * n + 1) As Double Zb(2 * n) = x Zb(2 * n + 1) = y Next i Next j ElseIf Obj.Closed = True Then 處理閉合曲線 For j = 0 To m Select Case j Case Is m - 2 a0 = j a1 = j + 1 a2 = j + 2 a3 = j + 3 Case m - 2 a0 = m - 2 a1 = m - 1 a2 = m a3 = 0 Case m - 1 a0

10、= m - 1 a1 = m a2 = 0 a3 = 1 Case m a0 = m a1 = 0 a2 = 1 a3 = 2 End Selecd For i = 0 To 3 3是個(gè)參數(shù)，參數(shù)越大等高線上節(jié)點(diǎn)越多 T = i / 3 f0 = 1 / 6 * (-T 3 + 3 * T 2 - 3 * T + 1) f1 = 1 / 6 * (3 * T 3 - 6 * T 2 + 4) f2 = 1 / 6 * (-3 * T 3 + 3 * T 2 + 3 * T + 1) f3 = 1 / 6 * T 3 x = f0 * Px(a0) + f1 * Px(a1) + f2 * Px

11、(a2) + f3 * Px(a3) y = f0 * Py(a0) + f1 * Py(a1) + f2 * Py(a2) + f3 * Py(a3) n = n + 1 ReDim Preserve Zb(2 * n + 1) As Double Zb(2 * n) = x Zb(2 * n + 1) = y Next Next j End If 分閉合曲線與不閉合曲線 gc = obj.Elevation Dim Pts() As Double Call ThreeFiltrates(Zb, Pts) 對(duì)生成的坐標(biāo)過濾 Call AddNewdgx(Pts, Obj) Obj.Delet

12、e Next kEnd Sub 圖3圖3是用上面的例子用此程序處理后的效果圖。等線的過濾： CASS中也有等高線過濾的功能，即“復(fù)合線濾波”，但此功能并不能很好的完成既過濾掉等高線上多余的節(jié)點(diǎn)又保持等線的原來位置，例如當(dāng)濾波值設(shè)小了過濾不完等高線上多余的節(jié)點(diǎn)；大了又改變了等高線的位置，如圖1，圖2所示。 圖4 圖5圖4為過濾前的等高線，圖6是用“復(fù)合線濾波”過濾后的等高線。于是作者想到了另一個(gè)算法即能過濾掉等高線上節(jié)點(diǎn)又能保持等高線位置。如圖6所示。 圖6我們從等高線的起點(diǎn)開始，設(shè)當(dāng)前點(diǎn)為”n”,第二個(gè)為”n+1”,第三個(gè)為”n+2”,我們把三個(gè)點(diǎn)組成一個(gè)三角形，算出”n+1”點(diǎn)到底邊的垂直距

13、離”h”。我們可以設(shè)置一個(gè)參數(shù)，用于決定是否保留”n+1”點(diǎn)，當(dāng)h大于這參數(shù)時(shí)去掉，小于時(shí)保留。計(jì)算”h”的公式：其中”area”為上面三角形的面積,”s”為三角形的周長，”a”為點(diǎn)”n”到”n+1”的邊，”b”為”n+1”到”n+2”的邊，”c”為點(diǎn)”n+2”到”n”的邊。此算法的好處在于可以保證等高線的位置不變又可以充分的過濾掉多余的節(jié)點(diǎn)，例如我們把參數(shù)設(shè)為“0.2”，等高線部分位置挪動(dòng)0.2m，對(duì)等高線來說沒有影響，但是卻可以過濾掉絕大過多數(shù)多余的節(jié)點(diǎn)。這個(gè)算法已經(jīng)用代碼實(shí)現(xiàn)，下面是部分源碼：Public Sub FilterDgxCoordinates(ByVal Coords As

14、Variant, ByRef Pts() As Double) Dim i, j, m, n As Integer m = (UBound(Coords) + 1) / 2 - 1 等高線上節(jié)點(diǎn)個(gè)數(shù)減一,因?yàn)槭菑?開始循環(huán) n = 2 * m + 1 Dim Zb() As Double ReDim Zb(n) As Double For i = 0 To n 先把值傳遞過來 Zb(i) = Coords(i) Next i If m 0 Then Pt(0) = Zb(2 * i) Pt(1) = Zb(2 * i + 1) MsgBox Pt(0) ? Pt1(0) = Zb(2 * i

15、+ 2) Pt1(1) = Zb(2 * i + 3) MsgBox Pt1(0) ? Pt2(0) = Zb(2 * i + 4) Pt2(1) = Zb(2 * i + 5) MsgBox Pt2(0) ? a = Distance(Pt, Pt1) c = Distance(Pt, Pt2) b = Distance(Pt1, Pt2) P = (a + b + c) / 2 If P - c 0.001 Then 判斷p和c是否相等不能用等號(hào)，因?yàn)樗鼈兪恰癲ouble”型變量 Area = 0 Cx = 0 Else Area = Sqr(P * (P - a) * (P - b) * (P - c) Cx = 2 * Area / c End If If Cx 0 Then n = n + 1 ReDim Preserve Pts(2 * n + 1) As Double Pts(2 * n) = Zb(2 * i) Pts(2 * n + 1) = Zb(2 * i + 1) End If Next iEnd Sub 圖7圖7即是上面例子中的等高線用此方法過濾后的效果圖。總結(jié)：等高線的光滑中，第一步加點(diǎn)是為了保證使用擬合公式后等高線的位置不至移動(dòng)太大，可以保持在正確的位置上同時(shí)初步光滑；第二步使用擬合公式擬合是讓等高線完全光滑。大量的事實(shí)也