离散点生成三角网的VBA程序

minocean

下面是一个读入离散点文件，根据Delaunay算法绘制三角网格的VBA程序。但是程序运行结果一直有问题~~~~~比如应该绘制1000个网格，可是每次就只绘制了一小部分，不知道哪里出错了……
请大家帮忙看看，给点意见：）



Option Base 1
Public Sub TIN_DEM()
'读取数据
Dim a, b, n As Long, i As Long, fname As String
Dim Point(), ordata() As Double
' 输入生成TIN网的数据
fname = InputBox("请输入高程文件位置如：D:\1.dat")

'获得文本文件的行数，实际行数为n+1
Open fname For Binary As #1
a = StrConv(InputB(LOF(1), 1), vbUnicode)
Close #1
b = Split(a, vbCrLf)
n = UBound(b)
ReDim ordata(n + 1, 5)
ReDim Point(n + 2, 3)
'将动态数组定义为二维数组
'开始从文件中读取数据，赋值给二维数组ordata
i = 1
Open fname For Input As #1
Do While Not EOF(1)
Input #1, ordata(i, 1), ordata(i, 2), ordata(i, 3), ordata(i, 4), ordata(i, 5)
Point(i, 1) = ordata(i, 3)
Point(i, 2) = ordata(i, 4)
Point(i, 3) = ordata(i, 5)
i = i + 1
Loop
Close #1

'选取源点，并计算最近点 k1
ReDim Tri(1000, 4) As Long, cr_tri(4) As Long
Dim dist, distmin As Double, k1, k2, k3 As Long
distmin = 1E+17
For i = 2 To n + 1
dist = distance(Point(1, 1), Point(1, 2), 0, Point(i, 1), Point(i, 2), 0)
If dist < distmin Then
distmin = dist
k1 = i
End If
Next i
'根据最大角原则选取三角形另外一个点
Dim MAXangle As Double
  MAXangle = 1
For i = 2 To n + 1
If i = k1 Then
i = i + 1
End If
If Cangle(Point(i, 1), Point(i, 2), Point(k1, 1), Point(k1, 2), Point(1, 1), Point(1, 2)) < MAXangle Then
  MAXangle = Cangle(Point(i, 1), Point(i, 2), Point(k1, 1), Point(k1, 2), Point(1, 1), Point(1, 2))
  k2 = i
End If
Next i
'首个三角形包含的点的记录
Tri(1, 1) = MIN(1, k1, k2): Tri(1, 2) = Sort2(k1, k2, 1): Tri(1, 3) = MAX(k1, k2, 1)
cr_tri(4) = 1

'记录当前工作区三角形
Dim j As Long
j = 1
Do While cr_tri(4) >= j
cr_tri(1) = Tri(j, 1): cr_tri(2) = Tri(j, 2): cr_tri(3) = Tri(j, 3)
j = j + 1
'以当前工作区三角形三边为起点寻找异侧角度最大的点
Dim judge, sign As Long, angle1, Mangle As Double
Mangle = 1
For i = 1 To n + 1
If i = cr_tri(1) Then
i = i + 1
End If
If i = cr_tri(2) Then
i = i + 1
End If
  judge = Judgement(Point(cr_tri(1), 1), Point(cr_tri(1), 2), Point(cr_tri(2), 1), Point(cr_tri(2), 2), Point(cr_tri(3), 1), Point(cr_tri(3), 2), Point(i, 1), Point(i, 2))
If judge = 0 Then
  angle1 = Cangle(Point(i, 1), Point(i, 2), Point(cr_tri(1), 1), Point(cr_tri(1), 2), Point(cr_tri(2), 1), Point(cr_tri(2), 2))
If angle1 < Mangle Then
  Mangle = angle1
  k3 = i
End If
End If
Next i
'判断新三角形是否与原有三角形全等
If Mangle < 0.9848 Then
Tri(cr_tri(4) + 1, 1) = MIN(cr_tri(1), cr_tri(2), k3): Tri(cr_tri(4) + 1, 2) = Sort2(cr_tri(1), cr_tri(2), k3): Tri(cr_tri(4) + 1, 3) = MAX(cr_tri(1), cr_tri(2), k3)
cr_tri(4) = cr_tri(4) + 1
sign = cr_tri(4)
For i = 1 To cr_tri(4) - 1
  If Tri(sign, 1) = Tri(i, 1) Then
   If Tri(sign, 2) = Tri(i, 2) Then
    If Tri(sign, 3) = Tri(i, 3) Then
     cr_tri(4) = cr_tri(4) - 1
     End If
    End If
   End If
Next i
End If

Mangle = 1
For i = 1 To n + 1
If i = cr_tri(1) Then
i = i + 1
End If
If i = cr_tri(3) Then
i = i + 1
End If
  judge = Judgement(Point(cr_tri(1), 1), Point(cr_tri(1), 2), Point(cr_tri(3), 1), Point(cr_tri(3), 2), Point(cr_tri(2), 1), Point(cr_tri(2), 2), Point(i, 1), Point(i, 2))
If judge = 0 Then
  angle1 = Cangle(Point(i, 1), Point(i, 2), Point(cr_tri(1), 1), Point(cr_tri(1), 2), Point(cr_tri(3), 1), Point(cr_tri(3), 2))
If angle1 < Mangle Then
  Mangle = angle1
  k3 = i
End If
End If
Next i
'判断新三角形是否与原三角形全等
If Mangle < 0.9848 Then
Tri(cr_tri(4) + 1, 1) = MIN(cr_tri(1), cr_tri(3), k3): Tri(cr_tri(4) + 1, 2) = Sort2(cr_tri(1), cr_tri(3), k3): Tri(cr_tri(4) + 1, 3) = MAX(cr_tri(1), cr_tri(3), k3)
cr_tri(4) = cr_tri(4) + 1
sign = cr_tri(4)
For i = 1 To cr_tri(4) - 1
  If Tri(sign, 1) = Tri(i, 1) Then
   If Tri(sign, 2) = Tri(i, 2) Then
    If Tri(sign, 3) = Tri(i, 3) Then
     cr_tri(4) = cr_tri(4) - 1
     End If
    End If
   End If
Next i
End If

Mangle = 1
For i = 1 To n + 1
If i = cr_tri(2) Then
i = i + 1
End If
If i = cr_tri(3) Then
i = i + 1
End If
  judge = Judgement(Point(cr_tri(2), 1), Point(cr_tri(2), 2), Point(cr_tri(3), 1), Point(cr_tri(3), 2), Point(cr_tri(1), 1), Point(cr_tri(1), 2), Point(i, 1), Point(i, 2))
If judge = 0 Then
  angle1 = Cangle(Point(i, 1), Point(i, 2), Point(cr_tri(2), 1), Point(cr_tri(2), 2), Point(cr_tri(3), 1), Point(cr_tri(3), 2))
If angle1 < Mangle Then
  Mangle = angle1
  k3 = i
End If
End If
Next i
'判断新三角形是否与原三角形全等
If Mangle < 0.9848 Then
Tri(cr_tri(4) + 1, 1) = MIN(cr_tri(2), cr_tri(3), k3): Tri(cr_tri(4) + 1, 2) = Sort2(cr_tri(2), cr_tri(3), k3): Tri(cr_tri(4) + 1, 3) = MAX(cr_tri(2), cr_tri(3), k3)
cr_tri(4) = cr_tri(4) + 1
sign = cr_tri(4)
For i = 1 To cr_tri(4) - 1
  If Tri(sign, 1) = Tri(i, 1) Then
   If Tri(sign, 2) = Tri(i, 2) Then
    If Tri(sign, 3) = Tri(i, 3) Then
      cr_tri(4) = cr_tri(4) - 1
     End If
    End If
   End If
Next i
End If
Loop

'绘制TIN网
Dim line As AcadLine
Dim startpoint(3) As Double, endpoint(3) As Double
For i = 1 To cr_tri(4) - 1
startpoint(1) = Point(Tri(i, 1), 1): startpoint(2) = Point(Tri(i, 1), 2): startpoint(3) = Point(Tri(i, 1), 3)
endpoint(1) = Point(Tri(i, 2), 1): endpoint(2) = Point(Tri(i, 2), 2): endpoint(3) = Point(Tri(i, 2), 3)
Set line = ThisDrawing.ModelSpace.AddLine(startpoint, endpoint)
startpoint(1) = Point(Tri(i, 1), 1): startpoint(2) = Point(Tri(i, 1), 2): startpoint(3) = Point(Tri(i, 1), 3)
endpoint(1) = Point(Tri(i, 3), 1): endpoint(2) = Point(Tri(i, 3), 2): endpoint(3) = Point(Tri(i, 3), 3)
Set line = ThisDrawing.ModelSpace.AddLine(startpoint, endpoint)
startpoint(1) = Point(Tri(i, 2), 1): startpoint(2) = Point(Tri(i, 2), 2): startpoint(3) = Point(Tri(i, 2), 3)
endpoint(1) = Point(Tri(i, 3), 1): endpoint(2) = Point(Tri(i, 3), 2): endpoint(3) = Point(Tri(i, 3), 3)
Set line = ThisDrawing.ModelSpace.AddLine(startpoint, endpoint)
Next i

End Sub
'计算距离过程
Function distance(ByVal x1 As Double, ByVal y1 As Double, ByVal z1 As Double, ByVal x2 As Double, ByVal y2 As Double, ByVal z2 As Double) As Double
distance = Sqr((x1 - x2) ^ 2 + (y1 - y2) ^ 2 + (z1 - z2) ^ 2)
End Function
'计算角度过程
Function Cangle(ByVal a1 As Double, ByVal a2 As Double, ByVal b1 As Double, ByVal b2 As Double, ByVal c1 As Double, ByVal c2 As Double) As Double
Dim AA As Double, BB As Double, CC As Double
  AA = distance(b1, b2, 0, c1, c2, 0)
  BB = distance(a1, a2, 0, c1, c2, 0)
  CC = distance(a1, a2, 0, b1, b2, 0)
  Cangle = (BB ^ 2 + CC ^ 2 - AA ^ 2) / (2 * BB * CC)
End Function
'判断异侧过程
Function Judgement(ByVal x1 As Double, ByVal y1 As Double, ByVal x2 As Double, ByVal y2 As Double, ByVal x3 As Double, ByVal y3 As Double, ByVal d1 As Double, ByVal d2 As Double) As Long
  Dim F(2), a, b, c As Double
  If x2 <> x1 Then
a = (y2 - y1) / (x2 - x1)
b = (y1 * x2 - y2 * x1) / (x2 - x1)
F(1) = y3 - a * x3 - b
F(2) = d2 - a * d1 - b
c = F(1) * F(2)
  End If
If c >= 0 Then
Judgement = 1
ElseIf c < 0 Then
Judgement = 0
End If
End Function
'排序求三个点中居中值
Function Sort2(ByVal a1 As Long, ByVal a2 As Long, ByVal a3 As Long) As Long
Dim t As Long
If a1 < a2 Then
t = a1: a1 = a2: a2 = t
End If
If a1 < a3 Then
t = a1: a1 = a3: a3 = t
End If
If a2 < a3 Then
t = a2: a2 = a3: a3 = t
End If
Sort2 = a2
End Function
'排序求三个点中最小值
Function MIN(ByVal a1 As Long, ByVal a2 As Long, ByVal a3 As Long) As Long
Dim t As Long
If a1 < a2 Then
t = a1: a1 = a2: a2 = t
End If
If a1 < a3 Then
t = a1: a1 = a3: a3 = t
End If
If a2 < a3 Then
t = a2: a2 = a3: a3 = t
End If
MIN = a3
End Function
'排序求三个点中最大值
Function MAX(ByVal a1 As Long, ByVal a2 As Long, ByVal a3 As Long) As Long
Dim t As Long
If a1 < a2 Then
t = a1: a1 = a2: a2 = t
End If
If a1 < a3 Then
t = a1: a1 = a3: a3 = t
End If
If a2 < a3 Then
t = a2: a2 = a3: a3 = t
End If
MAX = a1
End Function

minocean

自己顶一个吧  

quanguang

你的高程点数据文件格式是什么，发个例子上来，我调试下

107627ZJ

quanguang 发表于 2011-11-27 09:48
你的高程点数据文件格式是什么，发个例子上来，我调试下

调试好了能给我一个吗？QQ252113063   非常需要

momokill

调试好了分享一下可以吗？邮箱：ffg8435@163.com

3xxx

高人啊。

江湖远人

用CIVIL 3D做这些只是小菜一碟！

3xxx

这个代码有完善版的么？

jinrongch

完善版的有吗？求分享jinrongch@163.com谢谢！！！

15633587506

这个代码现在能正确运行了吗？能不能给我一份，谢谢！！！575870955@qq.com