surfaces.bas

用浮动水平线算法绘剖面图！在环境工程领域或三维图形中可以看到！

Attribute VB_Name = "Module4"
Option Explicit
'视点方位角、俯视角
Dim AngXY As Single, AngZ As Single
'用于存放计算用户坐标系点到观察坐标系点的坐标值公式中的正余弦值
Dim AL1 As Single, AL2 As Single, AL3 As Single, AM1 As Single, AM2 As Single, AM3 As Single, AN1 As Single, AN2 As Single, AN3 As Single
'点的投影坐标
Dim XProj As Single, YProj As Single
'上、下浮动水平线一维数组，存放先前曲线上同一x对应的最大y值和最小y值
Dim Hmax(0 To 800) As Integer, Hmin(0 To 800) As Integer
'窗口坐标
Dim WindowMinX As Integer, WindowMaxX As Integer, WindowMinY As Integer, WindowMaxY As Integer
Dim WindowMinDX As Single, WindowMinDY As Single
Dim XprojMin As Single, XprojMax As Single, YprojMin As Single, YprojMax As Single
'观测物理量
Dim Xmin As Single, Xmax As Single, Ymin As Single, Ymax As Single
Dim DX As Single, DY As Single
Dim Zmin As Single, Zmax As Single, S As Single
'其它变量
Dim Xd As Integer, Yd As Integer
Dim Xg As Integer, Yg As Integer
Dim NX As Integer, NY As Integer
Dim IX As Integer, IY As Integer
Dim X As Single, Y As Single, Z As Single
'数据初始化
Private Sub InitData()
    Dim Aux As Single
    Dim I As Integer

    DX = (Xmax - Xmin) / NX 'X方向增量
    DY = (Ymax - Ymin) / NY 'Y方向增量

    '视口赋边界值
    XprojMin = 1E+37
    XprojMax = -1E+37
    YprojMin = 1E+37
    YprojMax = -1E+37
    Xg = -1
    Yg = -1
    Xd = -1
    Yd = -1

    '水平线数组赋初值
    For I = 0 To WindowMaxX
        Hmax(I) = 0
        Hmin(I) = WindowMaxY
    Next I

    If (AngXY > 270) Then
    ElseIf (AngXY > 180) Then
        Aux = Xmin
        Xmin = Xmax
        Xmax = Aux
        DX = -DX
    
        Aux = Ymin
        Ymin = Ymax
        Ymax = Aux
        DY = -DY
    ElseIf (AngXY > 90) Then
        Aux = Xmin
        Xmin = Xmax
        Xmax = Aux
        DX = -DX
        DY = -DY
    End If
End Sub
'计算新坐标系坐标轴的方向余弦
'AngXY  方位角
'AngZ   俯视角
Private Sub InitProject()
    Dim Pi As Single

    Pi = 4# * Atn(1#) / 180#
    AL1 = -Sin(AngZ * Pi) * Sin(AngXY * Pi)
    AL2 = -Sin(AngZ * Pi) * Cos(AngXY * Pi)
    AL3 = Cos(AngZ * Pi)
    AM1 = Cos(AngXY * Pi)
    AM2 = -Sin(AngXY * Pi)
    AM3 = 0#
    AN1 = -AL3 * AM2
    AN2 = AL3 * AM1
    AN3 = Sin(AngZ * Pi)
End Sub
'求点的透视投影和平行投影坐标值
Private Sub Project(X As Single, Y As Single, Z As Single)
    XProj = AM1 * X + AM2 * Y
    YProj = AN1 * X + AN2 * Y + AN3 * (Z - Zmin) * S
End Sub
'此函数对大于0的数返回1，小于0的数返回-1，等于0的数仍取0，它用来确定填充水平线数组的方向
Private Function Sign(X As Single) As Integer
    Dim Sn As Integer
    
    If (X > 0) Then
        Sn = 1
    ElseIf (X < 0) Then
        Sn = -1
    Else
        Sn = 0
    End If
    Sign = Sn
End Function
'求出投影后的窗口
Private Sub FindWindow(GridZ() As Single)
    For IY = 1 To NY
        Y = Ymin + (IY - 1#) * DY
        For IX = NX To 1 Step -1
            X = Xmin + (IX - 1) * DX
        
            Z = GridZ(IY, IX)       '取观测值
            Call Project(X, Y, Z)   '求点的投影

            '求投影后的窗口边界值
            If (XProj < XprojMin) Then XprojMin = XProj
            If (XProj > XprojMax) Then XprojMax = XProj
            If (YProj < YprojMin) Then YprojMin = YProj
            If (YProj > YprojMax) Then YprojMax = YProj
        Next IX
    Next IY
End Sub
'此函数用于设置屏幕x方向和y方向的比例因子为其较小比例因子
'（比例因子就是窗口内x方向或y方向增加一个单位，屏幕视口内x方向或y方向的增加量）
Private Sub CalculateScale()
WindowMinDX = (WindowMaxX - WindowMinX) / (XprojMax - XprojMin) 'x方向的比例因子
WindowMinDY = (WindowMaxY - WindowMinY) / (YprojMax - YprojMin) 'y方向的比例因子
'取小的比例因子
If (WindowMinDX < WindowMinDY) Then
    WindowMinDY = WindowMinDX
Else
    WindowMinDX = WindowMinDY
End If
End Sub
'求两个数中的较小者
Private Function Min(X1 As Integer, X2 As Integer) As Integer
If (X1 < X2) Then
    Min = X1
Else
    Min = X2
End If
End Function
'求两个数中的较大者
Private Function Max(X1 As Integer, X2 As Integer) As Integer
If (X1 > X2) Then
    Max = X1
Else
    Max = X2
End If
End Function


'给上、下水平线数组元素赋值，在可见点与可见点或可见点与交点之间，用线值填充数组元素值
Private Sub Horizon(X1 As Integer, Y1 As Integer, X2 As Integer, Y2 As Integer)
    Dim Slope As Single, Delta As Single
    Dim XX As Integer, YY As Integer
    Dim DeltaT As Integer

    Delta = X2 - X1
    DeltaT = Sign(Delta)
    If (DeltaT = 0) Then '无须插值
        '直接填充数组元素
        Hmax(X2 + 1) = Max(Hmax(X2), Y2)
        Hmin(X2 + 1) = Min(Hmin(X2), Y2)
    Else
        Slope = (Y2 - Y1) / (X2 - X1)
        For XX = X2 + 1 To X1
            YY = Slope * (XX - X1) + Y1 '线性插值法求值
            '填充数组元素
            Hmax(XX) = Max(Hmax(XX), YY)
            Hmin(XX) = Min(Hmin(XX), YY)
        Next XX
    End If
End Sub
'检测可见性，由上水平线数组得知可见，用1表示；由下水平线数组得知可见，用-1表示;不可见用0表示
Private Function Visibility(X As Integer, Y As Integer) As Integer
    Dim Vis As Integer
    
    If (Y < Hmax(X) And Y > Hmin(X)) Then '不可见
        Vis = 0
    ElseIf (Y >= Hmax(X)) Then '上方可见
        Vis = 1
    Else '下方可见
        Vis = -1
    End If
    Visibility = Vis
End Function
'求当前曲线上可见性发生变化的两点之间线段与先前曲线的交点
Private Sub Inter(X1 As Integer, Y1 As Integer, X2 As Integer, Y2 As Integer, TabAux() As Integer, Xi As Integer, Yi As Integer)
    Dim Y2t As Single, Y1t As Single, Den As Single
    Dim Xii As Single, Yii As Single

    If (X2 = X1) Then '两点为同一点
        Xii = X2 '任取一为交点X坐标值
        Yii = TabAux(X2) '取相应Y坐标
    Else
        Den = Y2 - Y1 - TabAux(X2) + TabAux(X1) '交点公式中的分母
        If (Den <> 0) Then
            '交点的X、Y坐标
            Y1t = Y1
            Y2t = Y2
            Xii = (CDbl(X1) * (Y2 - TabAux(X2)) + CDbl(X2) * (TabAux(X1) - Y1)) / Den
            Yii = (Y2t * TabAux(X1) - Y1t * TabAux(X2)) / Den
        Else
            Xii = X2
            Yii = Y2
        End If
    End If
    Xi = Xii + 0.5
    Yi = Yii + 0.5
End Sub
'绘出函数式曲面的图形
Private Sub Drawing(GridZ() As Single)
Dim VisCur As Integer, VisPre As Integer
Dim XPre As Integer, YPre As Integer, XCur As Integer, YCur As Integer
Dim Xi As Integer, Yi As Integer

'对NY方向网格循环
For IY = 1 To NY
    Y = Ymin + (IY - 1) * DY        '取-Y值
    X = Xmax                        '取一X值
    Z = GridZ(IY, NX)
    Call Project(X, Y, Z)           '对点进行投影
    XPre = 0.5 + (XProj - XprojMin) * WindowMinDX + WindowMinX
    YPre = 0.5 + (YProj - YprojMin) * WindowMinDY + WindowMinY
    Call FillEdge(XPre, YPre, Xd, Yd)
    VisPre = Visibility(XPre, YPre)
    '对NX方向网格循环
    For IX = NX To 1 Step -1
        X = Xmin + (IX - 1) * DX    '取一X值
        Z = GridZ(IY, IX)           '取观测值
        Call Project(X, Y, Z)       '对点进行投影
        
        '求当前点的X、Y坐标
        XCur = 0.5 + (XProj - XprojMin) * WindowMinDX + WindowMinX
        YCur = 0.5 + (YProj - YprojMin) * WindowMinDY + WindowMinY
        
        VisCur = Visibility(XCur, YCur)
        
        If (Hmax(XCur) = 0 Or Hmin(XCur) = WindowMaxY) Then VisCur = VisPre
        If (VisCur = VisPre) Then
            If (VisCur = 1 Or VisCur = -1) Then '可见
                If (XCur >= 0) Then
                    Call PlotLine(XPre, YPre, XCur, YCur)
                ElseIf (YCur >= 0) Then
                    Call PlotLine(XPre, YPre, XPre, YCur)
                Else
                    Call PlotLine(XPre, YPre, XPre, YPre)
                End If
                '填充水平线数组
                Call Horizon(XPre, YPre, XCur, YCur)
            End If
        Else
            If (VisCur = 0) Then '当前点不可见
                If (VisPre = 1) Then '先前点可见
                    '------------相交情形之一------------
                    Call Inter(XPre, YPre, XCur, YCur, Hmax, Xi, Yi)
                Else
                    '------------相交情形之二------------
                    Call Inter(XPre, YPre, XCur, YCur, Hmin, Xi, Yi)
                End If
                If (Xi >= 0) Then '绘出可见线
                    Call PlotLine(XPre, YPre, Xi, Yi)
                ElseIf (Yi >= 0) Then '绘出可见线
                    Call PlotLine(XPre, YPre, XPre, Yi)
                Else
                    Call PlotLine(XPre, YPre, XPre, YPre)
                End If
                Call Horizon(XPre, YPre, Xi, Yi)
            Else
                If (VisCur = 1) Then '当前点可见
                    If (VisPre = 0) Then '先前点不可见
                        '------------相交情形之三------------
                        Call Inter(XPre, YPre, XCur, YCur, Hmax, Xi, Yi)
                        If (Xi >= 0) Then '绘出可见线
                            Call PlotLine(Xi, Yi, XCur, YCur)
                        ElseIf (Yi >= 0) Then '绘出可见线
                            Call PlotLine(XCur, Yi, XCur, YCur)
                        Else
                            Call PlotLine(XCur, YCur, XCur, YCur)
                        End If
                        Call Horizon(Xi, Yi, XCur, YCur)
                    Else
                        '------------相交情形之四------------
                        Call Inter(XPre, YPre, XCur, YCur, Hmin, Xi, Yi)
                        If (Xi >= 0) Then '绘出可见线
                            Call PlotLine(XPre, YPre, Xi, Yi)
                        ElseIf (Yi >= 0) Then '绘出可见线
                            Call PlotLine(XPre, YPre, XPre, Yi)
                        Else
                            Call PlotLine(XPre, YPre, XPre, YPre)
                        End If
                        Call Horizon(XPre, YPre, Xi, Yi)
                        
                        '------------相交情形之五------------
                        Call Inter(XPre, YPre, XCur, YCur, Hmax, Xi, Yi)
                        If (Xi >= 0) Then '绘出可见线
                            Call PlotLine(Xi, Yi, XCur, YCur)
                        ElseIf (Yi >= 0) Then '绘出可见线
                            Call PlotLine(XCur, Yi, XCur, YCur)
                        Else
                            Call PlotLine(XCur, YCur, XCur, YCur)
                        End If
                        Call Horizon(Xi, Yi, XCur, YCur)
                    End If
                Else
                    If (VisPre = 0) Then '先前点不可见
                        '------------相交情形之六------------
                        Call Inter(XPre, YPre, XCur, YCur, Hmin, Xi, Yi)
                        If (Xi >= 0) Then '绘出可见线
                            Call PlotLine(Xi, Yi, XCur, YCur)
                        ElseIf (Yi >= 0) Then '绘出可见线
                            Call PlotLine(XCur, Yi, XCur, YCur)
                        Else
                            Call PlotLine(XCur, YCur, XCur, YCur)
                        End If
                        Call Horizon(Xi, Yi, XCur, YCur)
                    Else
                        '------------相交情形之七------------
                        Call Inter(XPre, YPre, XCur, YCur, Hmax, Xi, Yi)
                        If (Xi >= 0) Then '绘出可见线
                            Call PlotLine(XPre, YPre, Xi, Yi)
                        ElseIf (Yi >= 0) Then '绘出可见线
                            Call PlotLine(XPre, YPre, XPre, Yi)
                        Else
                            Call PlotLine(XPre, YPre, XPre, YPre)
                        End If
                        Call Horizon(XPre, YPre, Xi, Yi)
                            
                        '------------相交情形之八------------
                        Call Inter(XPre, YPre, XCur, YCur, Hmin, Xi, Yi)
                        If (Xi >= 0) Then '绘出可见线
                            Call PlotLine(Xi, Yi, XCur, YCur)
                        ElseIf (Yi >= 0) Then '绘出可见线
                            Call PlotLine(XCur, Yi, XCur, YCur)
                        Else
                            Call PlotLine(XCur, YCur, XCur, YCur)
                        End If
                        Call Horizon(Xi, Yi, XCur, YCur)
                    End If
                End If
            End If
        End If
        VisPre = VisCur
        XPre = XCur
        YPre = YCur
    Next IX
    Call FillEdge(XCur, YCur, Xg, Yg) '填充边界
Next IY
End Sub
Private Sub PlotLine(X1 As Integer, Y1 As Integer, X2 As Integer, Y2 As Integer)
fSurface.Picture1.Line (X1, WindowMaxY - Y1)-(X2, WindowMaxY - Y2)
End Sub


'用于填充边界值
Private Sub FillEdge(X As Integer, Y As Integer, XLateral As Integer, YLateral As Integer)
If (XLateral <> -1) Then Call Horizon(XLateral, YLateral, X, Y)
XLateral = X
YLateral = Y
End Sub

'网格点值               GridZ
'X区间端点值            Xmin,Xmax
'Y区间端点值            Ymin,Ymax
'Y方向网格数            NY
'X方向网格数            NX
'方位角                 AngXY
'俯视角                 AngZ
'放大系数               S
Public Sub Surface(GridZ() As Single, XminT As Single, XmaxT As Single, YminT As Single, YmaxT As Single, ZminT As Single, ZmaxT As Single, NYt As Integer, NXt As Integer, AngXYt As Single, AngZt As Single, St As Single)
    Dim I As Integer, X As Single
    
    Xmin = XminT
    Xmax = XmaxT
    Ymin = YminT
    Ymax = YmaxT
    Zmin = ZminT
    Zmax = ZmaxT

    WindowMinX = 0
    WindowMinY = 0
    WindowMaxX = 500
    WindowMaxY = 500

    NY = NYt
    NX = NXt
        
    If (Xmax - Xmin > Ymax - Ymin) Then
        S = (Ymax - Ymin) / (Zmax - Zmin) * St
    Else
        S = (Xmax - Xmin) / (Zmax - Zmin) * St
    End If
    
    AngXY = AngXYt
    AngZ = AngZt
    
    Call InitData
    Call InitProject
    Call FindWindow(GridZ)
    Call CalculateScale
    Call Drawing(GridZ)
End Sub