matrix.bas

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

    ' the Y-Z plane.
    m3Identity R1
    D = Sqr(n1 * n1 + n2 * n2)
    R1(1, 1) = n2 / D
    R1(1, 2) = n1 / D
    R1(2, 1) = -R1(1, 2)
    R1(2, 2) = R1(1, 1)
    
    m3Identity R1i
    R1i(1, 1) = R1(1, 1)
    R1i(1, 2) = -R1(1, 2)
    R1i(2, 1) = -R1(2, 1)
    R1i(2, 2) = R1(2, 2)
    
    ' Rotate around the X-axis until the normal
    ' lies along the Y axis.
    m3Identity r2
    L = Sqr(n1 * n1 + n2 * n2 + n3 * n3)
    r2(2, 2) = D / L
    r2(2, 3) = -n3 / L
    r2(3, 2) = -r2(2, 3)
    r2(3, 3) = r2(2, 2)
    
    m3Identity R2i
    R2i(2, 2) = r2(2, 2)
    R2i(2, 3) = -r2(2, 3)
    R2i(3, 2) = -r2(3, 2)
    R2i(3, 3) = r2(3, 3)

    ' Reflect across the X-Z plane.
    m3Identity S
    S(2, 2) = -1

    ' Combine the matrices.
    m3MatMultiply M12, t, R1
    m3MatMultiply M34, r2, S
    m3MatMultiply M1234, M12, M34
    m3MatMultiply M56, R2i, R1i
    m3MatMultiply M567, M56, Ti
    m3MatMultiply M, M1234, M567
End Sub


' Create a transformation atrix for rotating
' through angle theta around a line passing
' through (p1, p2, p3) in direction <d1, d2, d3>.
' Theta is measured counterclockwise as you look
' down the line opposite the line's direction.
Public Sub m3LineRotate(M() As Single, ByVal p1 As Single, ByVal p2 As Single, ByVal p3 As Single, ByVal d1 As Single, ByVal d2 As Single, ByVal d3 As Single, ByVal theta 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 Rot3(1 To 4, 1 To 4) As Single  ' Rotate.
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 line is in
    ' the Y-Z plane.
    m3Identity R1
    D = Sqr(d1 * d1 + d2 * d2)
    R1(1, 1) = d2 / D
    R1(1, 2) = d1 / D
    R1(2, 1) = -R1(1, 2)
    R1(2, 2) = R1(1, 1)
    
    m3Identity R1i
    R1i(1, 1) = R1(1, 1)
    R1i(1, 2) = -R1(1, 2)
    R1i(2, 1) = -R1(2, 1)
    R1i(2, 2) = R1(2, 2)
    
    ' Rotate around the X-axis until the line
    ' lies along the Y axis.
    m3Identity r2
    L = Sqr(d1 * d1 + d2 * d2 + d3 * d3)
    r2(2, 2) = D / L
    r2(2, 3) = -d3 / L
    r2(3, 2) = -r2(2, 3)
    r2(3, 3) = r2(2, 2)
    
    m3Identity R2i
    R2i(2, 2) = r2(2, 2)
    R2i(2, 3) = -r2(2, 3)
    R2i(3, 2) = -r2(3, 2)
    R2i(3, 3) = r2(3, 3)

    ' Rotate around the line (Y axis).
    m3YRotate Rot3, theta

    ' Combine the matrices.
    m3MatMultiply M12, t, R1
    m3MatMultiply M34, r2, Rot3
    m3MatMultiply M1234, M12, M34
    m3MatMultiply M56, R2i, R1i
    m3MatMultiply M567, M56, Ti
    m3MatMultiply M, M1234, M567
End Sub

' Create a 3-D transformation matrix for scaling
' by scale factors Sx, Sy, and Sz.
Public Sub m3Scale(M() As Single, ByVal Sx As Single, ByVal Sy As Single, ByVal Sz As Single)
    m3Identity M
    M(1, 1) = Sx
    M(2, 2) = Sy
    M(3, 3) = Sz
End Sub

' Create a 3-D transformation matrix for
' translation by Tx, Ty, and Tz.
Public Sub m3Translate(M() As Single, ByVal Tx As Single, ByVal Ty As Single, ByVal Tz As Single)
    m3Identity M
    M(4, 1) = Tx
    M(4, 2) = Ty
    M(4, 3) = Tz
End Sub

' Create a 3-D transformation matrix for rotation
' around the X axis (angle measured in radians).
Public Sub m3XRotate(M() As Single, ByVal theta As Single)
    m3Identity M
    M(2, 2) = Cos(theta)
    M(3, 3) = M(2, 2)
    M(2, 3) = Sin(theta)
    M(3, 2) = -M(2, 3)
End Sub

' Create a 3-D transformation matrix for rotation
' around the Y axis (angle measured in radians).
Public Sub m3YRotate(M() As Single, ByVal theta As Single)
    m3Identity M
    M(1, 1) = Cos(theta)
    M(3, 3) = M(1, 1)
    M(3, 1) = Sin(theta)
    M(1, 3) = -M(3, 1)
End Sub

' Create a 3-D transformation matrix for rotation
' around the Z axis (angle measured in radians).
Public Sub m3ZRotate(M() As Single, ByVal theta As Single)
    m3Identity M
    M(1, 1) = Cos(theta)
    M(2, 2) = M(1, 1)
    M(1, 2) = Sin(theta)
    M(2, 1) = -M(1, 2)
End Sub

' Create a matrix that rotates around the Y axis
' so the point (x, y, z) lies in the X-Z plane.
Public Sub m3YRotateIntoXZ(Result() As Single, ByVal X As Single, ByVal Y As Single, ByVal Z As Single)
Dim D As Single

    m3Identity Result
    D = Sqr(X * X + Y * Y)
    Result(1, 1) = X / D
    Result(1, 2) = -Y / D
    Result(2, 1) = -Result(1, 2)
    Result(2, 2) = Result(1, 1)
End Sub

' Set copy = orig.
Public Sub m3MatCopy(copy() As Single, orig() As Single)
Dim i As Integer
Dim j As Integer

    For i = 1 To 4
        For j = 1 To 4
            copy(i, j) = orig(i, j)
        Next j
    Next i
End Sub

' Apply a transformation matrix to a point where
' the transformation may not have 0, 0, 0, 1 in
' its final column. Normalize only the X and Y
' components of the result to preserve the Z
' information.
Public Sub m3ApplyFull(V() As Single, M() As Single, Result() As Single)
Dim i As Integer
Dim j As Integer
Dim value As Single

    For i = 1 To 4
        value = 0#
        For j = 1 To 4
            value = value + V(j) * M(j, i)
        Next j
        Result(i) = value
    Next i
    
    ' Renormalize the point.
    ' Note that value still holds Result(4).
    If value <> 0 Then
        Result(1) = Result(1) / value
        Result(2) = Result(2) / value
        ' Do not transform the Z component.
    Else
        ' Make the Z value greater than that of
        ' the center of projection so the point
        ' will be clipped.
        Result(3) = INFINITY
    End If
    Result(4) = 1#
End Sub




' Apply a transformation matrix to a point.
Public Sub m3Apply(V() As Single, M() As Single, Result() As Single)
    Result(1) = V(1) * M(1, 1) + _
                V(2) * M(2, 1) + _
                V(3) * M(3, 1) + M(4, 1)
    Result(2) = V(1) * M(1, 2) + _
                V(2) * M(2, 2) + _
                V(3) * M(3, 2) + M(4, 2)
    Result(3) = V(1) * M(1, 3) + _
                V(2) * M(2, 3) + _
                V(3) * M(3, 3) + M(4, 3)
    Result(4) = 1#
End Sub

' Multiply two matrices together. The matrices
' may not contain 0, 0, 0, 1 in their last
' columns.
Public Sub m3MatMultiplyFull(Result() As Single, A() As Single, B() As Single)
Dim i As Integer
Dim j As Integer
Dim k As Integer
Dim value As Single

    For i = 1 To 4
        For j = 1 To 4
            value = 0#
            For k = 1 To 4
                value = value + A(i, k) * B(k, j)
            Next k
            Result(i, j) = value
        Next j
    Next i
End Sub
' Multiply two matrices together.
Public Sub m3MatMultiply(Result() As Single, A() As Single, B() As Single)
    Result(1, 1) = A(1, 1) * B(1, 1) + A(1, 2) * B(2, 1) + A(1, 3) * B(3, 1)
    Result(1, 2) = A(1, 1) * B(1, 2) + A(1, 2) * B(2, 2) + A(1, 3) * B(3, 2)
    Result(1, 3) = A(1, 1) * B(1, 3) + A(1, 2) * B(2, 3) + A(1, 3) * B(3, 3)
    Result(1, 4) = 0#
    Result(2, 1) = A(2, 1) * B(1, 1) + A(2, 2) * B(2, 1) + A(2, 3) * B(3, 1)
    Result(2, 2) = A(2, 1) * B(1, 2) + A(2, 2) * B(2, 2) + A(2, 3) * B(3, 2)
    Result(2, 3) = A(2, 1) * B(1, 3) + A(2, 2) * B(2, 3) + A(2, 3) * B(3, 3)
    Result(2, 4) = 0#
    Result(3, 1) = A(3, 1) * B(1, 1) + A(3, 2) * B(2, 1) + A(3, 3) * B(3, 1)
    Result(3, 2) = A(3, 1) * B(1, 2) + A(3, 2) * B(2, 2) + A(3, 3) * B(3, 2)
    Result(3, 3) = A(3, 1) * B(1, 3) + A(3, 2) * B(2, 3) + A(3, 3) * B(3, 3)
    Result(3, 4) = 0#
    Result(4, 1) = A(4, 1) * B(1, 1) + A(4, 2) * B(2, 1) + A(4, 3) * B(3, 1) + B(4, 1)
    Result(4, 2) = A(4, 1) * B(1, 2) + A(4, 2) * B(2, 2) + A(4, 3) * B(3, 2) + B(4, 2)
    Result(4, 3) = A(4, 1) * B(1, 3) + A(4, 2) * B(2, 3) + A(4, 3) * B(3, 3) + B(4, 3)
    Result(4, 4) = 1#
End Sub

' Compute the cross product of two vectors.
' Set <X, Y, Z> = <x1, y1, z1> X <x2, y2, z2>.
Public Sub m3Cross(ByRef X As Single, ByRef Y As Single, ByRef Z As Single, ByVal x1 As Single, ByVal y1 As Single, ByVal z1 As Single, ByVal x2 As Single, ByVal y2 As Single, ByVal z2 As Single)
    X = y1 * z2 - z1 * y2
    Y = z1 * x2 - x1 * z2
    Z = x1 * y2 - y1 * x2
End Sub

' Give the vector the indicated length.
Public Sub m3SizeVector(ByVal L As Single, ByRef X As Single, ByRef Y As Single, ByRef Z As Single)
    L = L / Sqr(X * X + Y * Y + Z * Z)
    X = X * L
    Y = Y * L
    Z = Z * L
End Sub