matrix.bas

纯软件的光线追踪算法 可以满足你的要求

Attribute VB_Name = "Matrix"
Option Explicit

Public Const PI = 3.14159265358979
Public Const INFINITY = 2147483647

Public Const m3Parallel = 0
Public Const m3Perspective = 1

Public Type Point3D
    Coord(1 To 4) As Single
    Trans(1 To 4) As Single
End Type

Public Type Segment3D
    pt1 As Integer
    pt2 As Integer
End Type
' Convert the spherical coordinates into
' Cartesian coordinates.
Public Sub m3SphericalToCartesian(ByVal R As Single, ByVal theta As Single, ByVal phi As Single, ByRef X As Single, ByRef Y As Single, ByRef Z As Single)
Dim r2 As Single
    
    ' Create a line to the center of projection.
    Y = R * Sin(phi)
    r2 = R * Cos(phi)
    X = r2 * Cos(theta)
    Z = r2 * Sin(theta)
End Sub
' Create a transformation matrix for an oblique
' projection onto the X-Y plane.
Public Sub m3ObliqueXY(M() As Single, ByVal S As Single, ByVal theta As Single)
    m3Identity M
    M(3, 1) = -S * Cos(theta)
    M(3, 2) = -S * Sin(theta)
    M(3, 3) = 0
End Sub


' Create a transformation matrix for orthographic
' projection along the X axis.
Public Sub m3OrthoSide(M() As Single)
    m3Identity M
    M(1, 1) = 0
    M(3, 1) = -1
    M(3, 3) = 0
End Sub
' Create a transformation matrix for orthographic
' projection along the Y axis.
Public Sub m3OrthoTop(M() As Single)
    m3Identity M
    M(2, 2) = 0
    M(3, 2) = -1
    M(3, 3) = 0
End Sub

' Create a transformation matrix for orthographic
' projection along the Z axis.
Public Sub m3OrthoFront(M() As Single)
    m3Identity M
    M(3, 3) = 0
End Sub

' Create an identity matrix.
Public Sub m3Identity(M() As Single)
Dim i As Integer
Dim j As Integer

    For i = 1 To 4
        For j = 1 To 4
            If i = j Then
                M(i, j) = 1
            Else
                M(i, j) = 0
            End If
        Next j
    Next i
End Sub

' Normalize a 3-D point vector.
Public Sub m3NormalizeCoords(ByRef X As Single, ByRef Y As Single, ByRef Z As Single, ByRef S As Single)
    X = X / S
    Y = Y / S
    Z = Z / S
    S = 1
End Sub

' Normalize a 3-D point vector.
Public Sub m3NormalizePoint(P() As Single)
Dim i As Integer
Dim value As Single

    value = P(4)
    For i = 1 To 3
        P(i) = P(i) / value
    Next i
    P(4) = 1
End Sub


' Normalize a 3-D transformation matrix.
Public Sub m3NormalizeMatrix(M() As Single)
Dim i As Integer
Dim j As Integer
Dim value As Single

    value = M(4, 4)
    For i = 1 To 4
        For j = 1 To 4
            M(i, j) = M(i, j) / value
        Next j
    Next i
End Sub




' Create a 3-D transformation matrix for a
' perspective projection along the Z axis onto
' the X-Y plane with focus at the origin and the
' center of projection at distance (0, 0, D).
Public Sub m3PerspectiveXZ(M() As Single, ByVal D As Single)
    m3Identity M
    If D <> 0 Then M(3, 4) = -1 / D
End Sub

' Create a 3-D transformation matrix for a
' projection with:
'       center of projection    (cx, cy, cz)
'       focus                   (fx, fy, fx)
'       UP vector               <ux, yx, uz>
' ptype should be m3Perspective or m3Parallel.
Public Sub m3Project(M() As Single, ByVal ptype As Integer, ByVal Cx As Single, ByVal Cy As Single, ByVal Cz As Single, ByVal Fx As Single, ByVal Fy As Single, ByVal Fz As Single, ByVal ux As Single, ByVal uy As Single, ByVal uz As Single)
Static M1(1 To 4, 1 To 4) As Single
Static M2(1 To 4, 1 To 4) As Single
Static M3(1 To 4, 1 To 4) As Single
Static M4(1 To 4, 1 To 4) As Single
Static M5(1 To 4, 1 To 4) As Single
Static M12(1 To 4, 1 To 4) As Single
Static M34(1 To 4, 1 To 4) As Single
Static M1234(1 To 4, 1 To 4) As Single
Dim sin1 As Single
Dim cos1 As Single
Dim sin2 As Single
Dim cos2 As Single
Dim sin3 As Single
Dim cos3 As Single
Dim A As Single
Dim B As Single
Dim C As Single
Dim d1 As Single
Dim d2 As Single
Dim d3 As Single
Dim up1(1 To 4) As Single
Dim up2(1 To 4) As Single

    ' Translate the focus to the origin.
    m3Translate M1, -Fx, -Fy, -Fz

    A = Cx - Fx
    B = Cy - Fy
    C = Cz - Fz
    d1 = Sqr(A * A + C * C)
    If d1 <> 0 Then
        sin1 = -A / d1
        cos1 = C / d1
    End If
    d2 = Sqr(A * A + B * B + C * C)
    If d2 <> 0 Then
        sin2 = B / d2
        cos2 = d1 / d2
    End If
    
    ' Rotate around the Y axis to place the
    ' center of projection in the Y-Z plane.
    m3Identity M2
    
    ' If d1 = 0 then the center of projection
    ' already lies in the Y axis and thus the Y-Z plane.
    If d1 <> 0 Then
        M2(1, 1) = cos1
        M2(1, 3) = -sin1
        M2(3, 1) = sin1
        M2(3, 3) = cos1
    End If
    
    ' Rotate around the X axis to place the
    ' center of projection in the Z axis.
    m3Identity M3
    
    ' If d2 = 0 then the center of projection
    ' lies at the origin. This makes projection
    ' impossible.
    If d2 <> 0 Then
        M3(2, 2) = cos2
        M3(2, 3) = sin2
        M3(3, 2) = -sin2
        M3(3, 3) = cos2
    End If
    
    ' Apply the rotations to the UP vector.
    up1(1) = ux
    up1(2) = uy
    up1(3) = uz
    up1(4) = 1
    m3Apply up1, M2, up2
    m3Apply up2, M3, up1

    ' Rotate around the Z axis to put the UP
    ' vector in the Y-Z plane.
    d3 = Sqr(up1(1) * up1(1) + up1(2) * up1(2))
    m3Identity M4
    
    ' If d3 = 0 then the UP vector is a zero
    ' vector so do nothing.
    If d3 <> 0 Then
        sin3 = up1(1) / d3
        cos3 = up1(2) / d3
        M4(1, 1) = cos3
        M4(1, 2) = sin3
        M4(2, 1) = -sin3
        M4(2, 2) = cos3
    End If
    
    ' Project.
    If ptype = m3Perspective And d2 <> 0 Then
        m3PerspectiveXZ M5, d2
    Else
        m3Identity M5
    End If

    ' Combine the transformations.
    m3MatMultiply M12, M1, M2
    m3MatMultiply M34, M3, M4
    m3MatMultiply M1234, M12, M34
    If ptype = m3Perspective Then
        m3MatMultiplyFull M, M1234, M5
    Else
        m3MatMultiply M, M1234, M5
    End If
End Sub



' Create a 3-D transformation matrix for a
' perspective projection with:
'       center of projection    (r, phi, theta)
'       focus                   (fx, fy, fx)
'       up vector               <ux, uy, uz>
' ptype should be m3Perspective or m3Parallel.
Public Sub m3PProject(M() As Single, ByVal ptype As Integer, ByVal R As Single, ByVal phi As Single, ByVal theta As Single, ByVal Fx As Single, ByVal Fy As Single, ByVal Fz As Single, ByVal ux As Single, ByVal uy As Single, ByVal uz As Single)
Dim Cx As Single
Dim Cy As Single
Dim Cz As Single
Dim r2 As Single

    ' Convert to Cartesian coordinates.
    Cy = R * Sin(phi)
    r2 = R * Cos(phi)
    Cx = r2 * Cos(theta)
    Cz = r2 * Sin(theta)
    m3Project M, ptype, Cx, Cy, Cz, Fx, Fy, Fz, ux, uy, uz
End Sub

' Create a transformation matrix for reflecting
' across the plane passing through (p1, p2, p3)
' with normal vector <n1, n2, n3>.
Public Sub m3Reflect(M() As Single, ByVal p1 As Single, ByVal p2 As Single, ByVal p3 As Single, ByVal n1 As Single, ByVal n2 As Single, ByVal n3 As Single)
Dim t(1 To 4, 1 To 4) As Single     ' Translate.
Dim R1(1 To 4, 1 To 4) As Single    ' Rotate 1.
Dim r2(1 To 4, 1 To 4) As Single    ' Rotate 2.
Dim S(1 To 4, 1 To 4) As Single     ' Reflect.
Dim R2i(1 To 4, 1 To 4) As Single   ' Unrotate 2.
Dim R1i(1 To 4, 1 To 4) As Single   ' Unrotate 1.
Dim Ti(1 To 4, 1 To 4) As Single    ' Untranslate.
Dim D As Single
Dim L As Single
Dim M12(1 To 4, 1 To 4) As Single
Dim M34(1 To 4, 1 To 4) As Single
Dim M1234(1 To 4, 1 To 4) As Single
Dim M56(1 To 4, 1 To 4) As Single
Dim M567(1 To 4, 1 To 4) As Single

    ' Translate the plane to the origin.
    m3Translate t, -p1, -p2, -p3
    m3Translate Ti, p1, p2, p3

    ' Rotate around Z-axis until the normal is in