Keskustelu: Koodit: VB.NET: Reaaliaikaista raytracingiä

Public Structure ColorFloat
    Dim r, g, b As Double
End Structure
Public Structure Surface
    Dim baseColor As ColorFloat
    Dim reflectivity As Double
End Structure
Public Structure LightSource
    Dim location As Vector.sVector
    Dim c As ColorFloat ' values from 0 to 1 are within normal range. above that makes light brighter
End Structure

Public Class frmRaytracer
    Inherits System.Windows.Forms.Form

    Dim backbuffer As Bitmap
    Dim a As Double

    Const numPrimitives As Integer = 15
    Dim primitives(numPrimitives) As Primitive.sPrimitive

    Const numLightSources As Integer = 4
    Dim lightSources(numLightSources) As LightSource

    Const numVertices As Integer = 4
    Dim vertices(numVertices) As Vector.sVector

    Dim directionTable() As Vector.sVector

    Dim cameraLocXDelta, cameraLocYDelta, cameraLocZDelta As Double
    Dim cameraRotX, cameraRoty, cameraRotz As Double

#Region " Windows Form Designer generated code "

    Public Sub New()
        MyBase.New()

        'This call is required by the Windows Form Designer.
        InitializeComponent()

        'Add any initialization after the InitializeComponent() call

    End Sub

    'Form overrides dispose to clean up the component list.
    Protected Overloads Overrides Sub Dispose(ByVal disposing As Boolean)
        If disposing Then
            If Not (components Is Nothing) Then
                components.Dispose()
            End If
        End If
        MyBase.Dispose(disposing)
    End Sub

    'Required by the Windows Form Designer
    Private components As System.ComponentModel.IContainer

    'NOTE: The following procedure is required by the Windows Form Designer
    'It can be modified using the Windows Form Designer.
    'Do not modify it using the code editor.
    Friend WithEvents picSurface As System.Windows.Forms.PictureBox
    <System.Diagnostics.DebuggerStepThrough()> Private Sub InitializeComponent()
        Me.picSurface = New System.Windows.Forms.PictureBox
        Me.SuspendLayout()
        '
        'picSurface
        '
        Me.picSurface.Location = New System.Drawing.Point(0, 0)
        Me.picSurface.Name = "picSurface"
        Me.picSurface.Size = New System.Drawing.Size(616, 520)
        Me.picSurface.TabIndex = 0
        Me.picSurface.TabStop = False
        '
        'frmRaytracer
        '
        Me.AutoScaleBaseSize = New System.Drawing.Size(5, 13)
        Me.ClientSize = New System.Drawing.Size(616, 518)
        Me.Controls.Add(Me.picSurface)
        Me.Name = "frmRaytracer"
        Me.Text = "Raytracer"
        Me.ResumeLayout(False)

    End Sub

#End Region

    Private Sub frmRaytracer_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
        Dim i As Integer

        ' allocate the ray direction lookup table
        directionTable = Tracer.GenerateRayDirectionTable
        SetupScene(primitives, numPrimitives, vertices, numVertices, lightSources, numLightSources)
        backbuffer = New Bitmap(picSurface.Width, picSurface.Height)
    End Sub

    Private Sub UpdateScene(ByVal primitives() As Primitive.sPrimitive, ByVal numPrimitives As Integer, ByVal vertices() As Vector.sVector, ByVal numVertices As Integer, ByVal lightSources() As LightSource, ByVal numLightSources As Integer, ByVal currTime As Double)
        Dim i As Integer
        For i = 0 To numPrimitives
            primitives(i).surface.baseColor.r = 128 * Math.Sin(currTime * i) + 128
            primitives(i).surface.baseColor.g = 128 * Math.Sin(currTime * i + 3 + Math.Sin(currTime)) + 128
            primitives(i).surface.baseColor.b = 128 * Math.Sin(currTime * i + 2 + 1.5 * Math.Sin(currTime)) + 128
            primitives(i).surface.reflectivity = 0.5 + 0.5 * Math.Sin(currTime + i)
        Next
        Dim normal As Vector.sVector
        Vector.VectorSetXYZ(normal, 0.14999999999999999, -1, 0)
        Vector.Rotate(normal, 0, currTime, 0)
        primitives(0).planeProperties.normal = normal

        Vector.VectorSetXYZ(vertices(0), 0, -100, 50)
        Vector.Rotate(vertices(0), currTime * 1.3999999999999999, currTime, currTime * 1.2352000000000001)
        Vector.VectorSetXYZ(vertices(1), 100, 100, 0)
        Vector.Rotate(vertices(1), currTime * 1.3999999999999999, currTime, currTime * 1.2352000000000001)
        Vector.VectorSetXYZ(vertices(2), 0, 100, 250)
        Vector.Rotate(vertices(2), currTime * 1.3999999999999999, currTime, currTime * 1.2352000000000001)
        Vector.VectorSetXYZ(vertices(3), -100, 100, 0)
        Vector.Rotate(vertices(3), currTime * 1.3999999999999999, currTime, currTime * 1.2352000000000001)

        primitives(5).cylinderProperties.radius = 60 + 40 * Math.Sin(currTime)

        For i = 6 To numPrimitives
            Vector.VectorSetXYZ(primitives(i).sphereProperties.center, 120 * Math.Sin(currTime + i * 2 + 0.56000000000000005), 120 * Math.Sin(currTime + i * 2 + 3 * Math.Sin(currTime * 0.20000000000000001)), 120 * Math.Sin(currTime + i * 5 + 1 + Math.Sin(currTime * 0.10000000000000001)))
            primitives(i).sphereProperties.radius = 30 * Math.Sin(currTime * 0.40000000000000002 + Math.Sin(currTime + i)) + 40
        Next

        For i = 0 To numLightSources

            lightSources(i).location.x = 250 * Math.Sin(currTime * 0.5 + i)
            lightSources(i).location.y = -100
            lightSources(i).location.z = 250 * Math.Cos(currTime * 0.5 + i)
            lightSources(i).c.r = 0.14999999999999999 * Math.Sin(currTime + i) + 0.29999999999999999
            lightSources(i).c.g = 0.14999999999999999 * Math.Sin(currTime + 3 + Math.Sin(currTime + i)) + 0.29999999999999999
            lightSources(i).c.b = 0.14999999999999999 * Math.Sin(currTime * i + 2 + 1.5 * Math.Sin(currTime)) + 0.29999999999999999
        Next
    End Sub
    Private Sub SetupScene(ByVal primitives() As Primitive.sPrimitive, ByVal numPrimitives As Integer, ByVal vertices() As Vector.sVector, ByVal numVertices As Integer, ByVal lightSources() As LightSource, ByVal numLightSources As Integer)
        Primitive.AssignPrimitiveType(primitives(0), Primitive.PrimitiveType.PLANE_TYPE)
        primitives(0).planeProperties.displacement = 150

        Primitive.AssignPrimitiveType(primitives(1), Primitive.PrimitiveType.TRIANGLE_TYPE)
        primitives(1).triangleProperties.v1 = vertices(0)
        primitives(1).triangleProperties.v2 = vertices(1)
        primitives(1).triangleProperties.v3 = vertices(2)

        Primitive.AssignPrimitiveType(primitives(2), Primitive.PrimitiveType.TRIANGLE_TYPE)
        primitives(2).triangleProperties.v1 = vertices(0)
        primitives(2).triangleProperties.v2 = vertices(3)
        primitives(2).triangleProperties.v3 = vertices(1)

        Primitive.AssignPrimitiveType(primitives(3), Primitive.PrimitiveType.TRIANGLE_TYPE)
        primitives(3).triangleProperties.v1 = vertices(0)
        primitives(3).triangleProperties.v2 = vertices(2)
        primitives(3).triangleProperties.v3 = vertices(3)

        Primitive.AssignPrimitiveType(primitives(4), Primitive.PrimitiveType.TRIANGLE_TYPE)
        primitives(4).triangleProperties.v1 = vertices(1)
        primitives(4).triangleProperties.v2 = vertices(3)
        primitives(4).triangleProperties.v3 = vertices(2)

        Primitive.AssignPrimitiveType(primitives(5), Primitive.PrimitiveType.CYLINDER_TYPE)
        primitives(5).cylinderProperties.axis = Cylinder.eAxis.Y_INFINITE
        Vector.VectorSetXYZ(primitives(5).cylinderProperties.center, -200, 0, 150)

        Dim i As Integer
        For i = 6 To numPrimitives
            Primitive.AssignPrimitiveType(primitives(i), Primitive.PrimitiveType.SPHERE_TYPE)
        Next
    End Sub

    Private Sub frmRaytracer_Activated(ByVal sender As Object, ByVal e As System.EventArgs) Handles MyBase.Activated
        Dim Cam As Camera.cCamera = New Camera.cCamera
        Camera.SetupCamera(Cam)
        Do While 1 = 1
            Cam.location.x += cameraLocXDelta
            Cam.location.y += cameraLocYDelta
            Cam.location.z += cameraLocZDelta

            Camera.CameraRotate(Cam, cameraRotX, cameraRotY, cameraRotZ)

            a += 0.10000000000000001
            If a >= 2 * Math.PI Then a = 0
            UpdateScene(primitives, numPrimitives, vertices, numVertices, lightSources, numLightSources, a * 0.5)
            Tracer.TraceScene(Cam, primitives, numPrimitives, lightSources, numLightSources, backbuffer, directionTable)
            picSurface.BackgroundImage = backbuffer
            picSurface.Refresh()
            Application.DoEvents()
        Loop
    End Sub

    Private Sub frmRaytracer_KeyDown(ByVal sender As Object, ByVal e As System.Windows.Forms.KeyEventArgs) Handles MyBase.KeyDown
        Select Case e.KeyCode
            Case Keys.Up
                cameraLocZDelta = 100
            Case Keys.Down
                cameraLocZDelta = -100
            Case Keys.Left
                cameraLocXDelta = -100
            Case Keys.Right
                cameraLocXDelta = 100
            Case Keys.NumPad4
                cameraRoty = -0.10000000000000001
            Case Keys.NumPad6
                cameraRoty = 0.10000000000000001
                'Case Keys.NumPad8
                '    cameraRotX = 0.10000000000000001
                'Case Keys.NumPad5
                '    cameraRotX = -0.10000000000000001
        End Select
    End Sub

    Private Sub frmRaytracer_KeyUp(ByVal sender As Object, ByVal e As System.Windows.Forms.KeyEventArgs) Handles MyBase.KeyUp
        Select Case e.KeyCode
            Case Keys.Up
                cameraLocZDelta = 0
            Case Keys.Down
                cameraLocZDelta = 0
            Case Keys.Left
                cameraLocXDelta = 0
            Case Keys.Right
                cameraLocXDelta = 0
            Case Keys.NumPad4
                cameraRoty = 0
            Case Keys.NumPad6
                cameraRoty = 0
                'Case Keys.NumPad8 ' bugi pallot katoo
                '    cameraRotX = 0
                'Case Keys.NumPad5
                '    cameraRotX = 0
        End Select
    End Sub
End Class
Public Class Tracer
    Public Shared Function GenerateRayDirectionTable() As Vector.sVector()

        Dim direction(640 * 480) As Vector.sVector
        Dim x, y As Integer
        For y = 0 To 479
            For x = 0 To 639
                direction(x + y * 640).x = x - 320
                direction(x + y * 640).y = y - 240
                direction(x + y * 640).z = 255 '58 + 0.40000000000000002 * Math.Sqrt((x - 320) * (x - 320) + (y - 240) * (y - 240)) ' fisheye
                Vector.VectorNormalize(direction(x + y * 640))
            Next
        Next
        Return direction
    End Function

    Public Shared Sub TraceScene(ByRef cam As Camera.cCamera, ByVal prims() As Primitive.sPrimitive, ByVal numPrimitives As Integer, ByVal lightSources() As LightSource, ByVal numLightSources As Integer, ByVal buffer As Bitmap, ByVal directionTable() As Vector.sVector)
        ' setup view rays
        Dim primaryRay As Ray.sRay
        primaryRay.origin.x = 0
        primaryRay.origin.y = 0
        primaryRay.origin.z = -570

        Dim xDiff As Integer = 80 ' equal to 1/2 of the vertical screen size
        Dim yDiff As Integer = 60

        Dim x, y As Integer

        For y = 240 - yDiff To 240 + yDiff
            For x = 320 - xDiff To 319 + xDiff
                primaryRay.direction = directionTable(x + (y << 9) + (y << 7)) ' implimenting the direction table added some fps @ 240x180....
                primaryRay = Camera.CalculateCameraRay(cam, x, y, directionTable)
                Dim ccolor As ColorFloat = TraceRay(-1, primaryRay, prims, numPrimitives, lightSources, numLightSources, 0)
                If ccolor.r > 255 Then ccolor.r = 255
                If ccolor.g > 255 Then ccolor.g = 255
                If ccolor.b > 255 Then ccolor.b = 255
                buffer.SetPixel(x, y, Color.FromArgb(ccolor.r, ccolor.g, ccolor.b))
            Next
        Next

    End Sub

    Public Shared Function TraceRay(ByVal ignorenum As Integer, ByVal pray As Ray.sRay, ByVal prims() As Primitive.sPrimitive, ByVal numPrimitives As Integer, ByVal lightsources() As LightSource, ByVal numLightsources As Integer, ByVal depth As Integer) As ColorFloat
        Dim returnColor As ColorFloat
        If depth > 5 Then ' prevent infinite reflection
            returnColor.r = 0
            returnColor.g = 0
            returnColor.b = 0
            Return returnColor
        End If

        Dim closestIntersectionDistance As Double = 1000000 ' an impossibly large value
        Dim closestIntersectedPrimitiveNum As Integer = -1
        Dim currResult As Ray.TraceResult

        Dim currPrimitiveNum As Integer
        For currPrimitiveNum = 0 To numPrimitives ' cycle through all of the spheres to find the closest interesction
            If currPrimitiveNum <> ignorenum Then

                currResult = Primitive.IntersectPrimitive(prims(currPrimitiveNum), pray)

                If currResult.hit Then
                    If currResult.distance < closestIntersectionDistance Then
                        closestIntersectionDistance = currResult.distance
                        closestIntersectedPrimitiveNum = currPrimitiveNum
                    End If
                End If
            End If
        Next

        If closestIntersectedPrimitiveNum = -1 Then ' nothing was intersected

            returnColor.r = 0
            returnColor.g = 0
            returnColor.b = 0

        Else
            returnColor = Shader.Shade(prims(closestIntersectedPrimitiveNum), _
                            closestIntersectedPrimitiveNum, _
                            pray, _
                            closestIntersectionDistance, _
                            lightsources, numLightsources, _
                            prims, _
                            numPrimitives, _
                            depth)
        End If
        Return returnColor

    End Function

    Public Shared Function IsShadowed(ByVal currPrimitiveNum As Integer, _
    ByVal raytolight As Ray.sRay, _
    ByVal distanceToLight As Double, _
    ByVal prims() As Primitive.sPrimitive, _
    ByVal numPrimitives As Integer) As Boolean

        Dim i As Integer
        Dim tr As Ray.TraceResult
        tr.hit = False
        ' check every other sphere
        For i = 0 To numPrimitives
            If i <> currPrimitiveNum Then ' dont self-shadow
                Select Case prims(i).type

                    Case Primitive.PrimitiveType.SPHERE_TYPE
                        tr = Sphere.IntersectSphere(prims(i).sphereProperties, raytolight)
                    Case Primitive.PrimitiveType.PLANE_TYPE
                        tr = Plane.IntersectPlane(prims(i).planeProperties, raytolight)
                    Case Primitive.PrimitiveType.CYLINDER_TYPE
                        tr = Cylinder.IntersectCylinder(prims(i).cylinderProperties, raytolight)
                    Case Primitive.PrimitiveType.CYLINDER_TYPE
                        tr = Triangle.IntersectTriangle(prims(i).triangleProperties, raytolight)
                End Select

                If tr.hit And tr.distance < distanceToLight Then Return True
            End If
        Next
        Return False
    End Function
End Class
Public Class Shader
    Public Shared Function CalculateIntersection(ByVal pray As Ray.sRay, ByVal distance As Double) As Vector.sVector
        Dim intersection As Vector.sVector
        ' calculate the location of the intersection between the primitive and the ray
        intersection.x = pray.origin.x + distance * pray.direction.x
        intersection.y = pray.origin.y + distance * pray.direction.y
        intersection.z = pray.origin.z + distance * pray.direction.z
        Return intersection
    End Function

    Public Shared Function CalculateLightRay(ByVal lightLoc As Vector.sVector, ByVal intersection As Vector.sVector, ByRef rayToLight As Ray.sRay) As Double ' puts the resulting vector in lightDir

        Dim lightDir As Vector.sVector = Vector.VectorSub(lightLoc, intersection)
        ' because we need the distance to the light for the shadow calculations, we will normalize lightDir manually
        Dim distanceToLight As Double = Math.Sqrt((lightDir.x * lightDir.x) + (lightDir.y * lightDir.y) + (lightDir.z * lightDir.z))
        Dim lightDirMagnitudeReciprocal As Double = 1 / distanceToLight
        lightDir.x *= lightDirMagnitudeReciprocal
        lightDir.y *= lightDirMagnitudeReciprocal
        lightDir.z *= lightDirMagnitudeReciprocal

        rayToLight.origin = intersection
        rayToLight.direction = lightDir

        Return distanceToLight
    End Function

    Public Shared Function CalculateLightingCoef(ByVal isShadowed As Boolean, ByVal directionToLight As Vector.sVector, ByVal normal As Vector.sVector) As Double

        If isShadowed Then
            ' no light reaches the intersection
            Return 0
        Else
            ' only calculate how much light reaches the intersection if it is not in shadow
            Dim lightCoef As Double = Vector.VectorDot(directionToLight, normal)
            If lightCoef < 0 Then lightCoef = 0
            Return lightCoef
        End If
    End Function

    Public Shared Function CalculateReflection(ByVal pray As Ray.sRay, ByVal intersection As Vector.sVector, _
    ByVal normal As Vector.sVector, ByVal currPrimitiveNum As Integer, ByVal prims() As Primitive.sPrimitive, _
    ByVal numPrimitives As Integer, ByVal lightsources() As LightSource, ByVal numLightSources As Integer, _
    ByVal depth As Integer) As ColorFloat

        ' R = I - 2(N.I)*N
        Dim reflectedRay As Ray.sRay
        Dim nDotI As Double = 2 * ((normal.x * pray.direction.x) + (normal.y * pray.direction.y) + (normal.z * pray.direction.z))
        reflectedRay.direction.x = pray.direction.x - (nDotI * normal.x)
        reflectedRay.direction.y = pray.direction.y - (nDotI * normal.y)
        reflectedRay.direction.z = pray.direction.z - (nDotI * normal.z)
        reflectedRay.origin = intersection

        Return Tracer.TraceRay(currPrimitiveNum, reflectedRay, prims, numPrimitives, lightsources, numLightSources, depth + 1)
    End Function

    Public Shared Function Shade(ByRef prim As Primitive.sPrimitive, _
    ByVal currPrimitiveNum As Integer, _
    ByVal pray As Ray.sRay, _
    ByVal distance As Double, _
    ByVal lightSources() As LightSource, ByVal numLightSources As Integer, _
    ByVal prims() As Primitive.sPrimitive, _
    ByVal numPrimitives As Integer, _
    ByVal depth As Integer) As ColorFloat

        Dim returnColor As ColorFloat
        returnColor.r = 0
        returnColor.g = 0
        returnColor.b = 0

        Dim intersection As Vector.sVector = CalculateIntersection(pray, distance)

        Dim normal As Vector.sVector = Primitive.CalculateNormal(prim, intersection)

        ' add specular components
        Dim i As Integer
        For i = 0 To numLightSources
            Dim rayToLight As Ray.sRay
            Dim distanceToLight As Double = CalculateLightRay(lightSources(i).location, intersection, rayToLight) ' sets rayToLight
            Dim lightCoef As Double = CalculateLightingCoef(Tracer.IsShadowed(currPrimitiveNum, rayToLight, distanceToLight, prims, numPrimitives), rayToLight.direction, normal)
            returnColor.r += prim.surface.baseColor.r * lightCoef * lightSources(i).c.r ' try checking first if lightCoef is 0... the check will probably be amortized over the
            returnColor.g += prim.surface.baseColor.g * lightCoef * lightSources(i).c.g ' cost of all these multiplications
            returnColor.b += prim.surface.baseColor.b * lightCoef * lightSources(i).c.b
        Next

        ' add reflective components
        If prim.surface.reflectivity <> 0 Then
            Dim reflectedColor As ColorFloat = CalculateReflection(pray, intersection, normal, currPrimitiveNum, prims, numPrimitives, lightSources, numLightSources, depth)
            returnColor.r += reflectedColor.r * prim.surface.reflectivity
            returnColor.g += reflectedColor.g * prim.surface.reflectivity
            returnColor.b += reflectedColor.b * prim.surface.reflectivity
        End If
        Return returnColor
    End Function
End Class
Public Class Ray
    Public Structure sRay
        Dim origin, direction As Vector.sVector
    End Structure

    Public Structure TraceResult
        Dim hit As Boolean
        Dim distance As Double
    End Structure
End Class
Public Class Camera
    Public Class cCamera
        Public location, right, up, forward As Vector.sVector ' forward, up and right define the direction of the camera
        Public rayRotationMatrix As Matrix.Matrix3x3 = New Matrix.Matrix3x3 ' the matrix to multiply a ray's direction to get its rotated direction
    End Class

    Public Shared Function CalculateCameraRay(ByVal cam As cCamera, ByVal x As Integer, ByVal y As Integer, ByVal directionTable() As Vector.sVector) As Ray.sRay
        Dim pray As Ray.sRay
        pray.origin = cam.location ' maybe optimize this statement
        pray.direction = directionTable(x + (y << 9) + (y << 7))

        pray.direction = Matrix.VectorMultMatrix(pray.direction, cam.rayRotationMatrix) ' rotate the ray
        Return pray
    End Function

    Public Shared Sub SetupCamera(ByRef cam As cCamera)
        cam.location.x = 0
        cam.location.y = 0
        cam.location.z = -2700
        Vector.VectorSetXYZ(cam.right, 1, 0, 0)
        Vector.VectorSetXYZ(cam.up, 0, -1, 0)
        Vector.VectorSetXYZ(cam.forward, 0, 0, 1)
        cam.rayRotationMatrix = Matrix.IdentityMatrix3x3() ' initally there is no rotation
    End Sub

    Public Shared Sub CameraRotate(ByRef cam As cCamera, ByVal thetaRight As Double, ByVal thetaUp As Double, ByVal thetaForward As Double)
        Dim rotationRight As Matrix.Matrix3x3 = New Matrix.Matrix3x3
        Dim rotationUp As Matrix.Matrix3x3 = New Matrix.Matrix3x3
        Dim rotationForward As Matrix.Matrix3x3 = New Matrix.Matrix3x3

        ' calculate rotation matrices about each axis that defines the camera's direction
        rotationRight = Matrix.CalculateArbitrayRotationMatrix(cam.right, thetaRight)
        rotationUp = Matrix.CalculateArbitrayRotationMatrix(cam.up, thetaUp)
        rotationForward = Matrix.CalculateArbitrayRotationMatrix(cam.forward, thetaForward)

        Dim combinedRotationMatrix As Matrix.Matrix3x3
        ' combinedeRotationMatrix = rotationRight * rotationUp * RotationForward
        'combinedRotationMatrix = MatrixMult3x3(MatrixMult3x3(rotationRight, rotationUp), rotationForward); gives a compiler error
        combinedRotationMatrix = Matrix.MatrixMult3x3(rotationRight, rotationUp)
        combinedRotationMatrix = Matrix.MatrixMult3x3(combinedRotationMatrix, rotationForward)

        ' rotate the defining axes by the combined rotation matrix
        cam.right = Matrix.VectorMultMatrix(cam.right, combinedRotationMatrix)
        cam.up = Matrix.VectorMultMatrix(cam.up, combinedRotationMatrix)
        cam.forward = Matrix.VectorMultMatrix(cam.forward, combinedRotationMatrix)

        ' combine the current ray rotation matrix with the new matrix
        cam.rayRotationMatrix = Matrix.MatrixMult3x3(cam.rayRotationMatrix, combinedRotationMatrix)
    End Sub

End Class
Public Class Primitive
    Enum PrimitiveType
        SPHERE_TYPE
        PLANE_TYPE
        CYLINDER_TYPE
        TRIANGLE_TYPE
    End Enum
    Public Structure sPrimitive
        Dim surface As surface
        Dim type As PrimitiveType
        Dim planeProperties As Plane.PlaneProperties
        Dim sphereProperties As Sphere.SphereProperties
        Dim cylinderProperties As Cylinder.CylinderProperties
        Dim triangleProperties As Triangle.TriangleProperties
    End Structure

    Public Shared Sub AssignPrimitiveType(ByRef prim As sPrimitive, ByVal type As PrimitiveType)
        Select Case type
            Case PrimitiveType.SPHERE_TYPE
                prim.type = PrimitiveType.SPHERE_TYPE
                prim.sphereProperties = New Sphere.SphereProperties
            Case PrimitiveType.PLANE_TYPE
                prim.type = PrimitiveType.PLANE_TYPE
                prim.planeProperties = New Plane.PlaneProperties
            Case PrimitiveType.CYLINDER_TYPE
                prim.type = PrimitiveType.CYLINDER_TYPE
                prim.cylinderProperties = New Cylinder.CylinderProperties
            Case PrimitiveType.TRIANGLE_TYPE
                prim.type = PrimitiveType.TRIANGLE_TYPE
                prim.triangleProperties = New Triangle.TriangleProperties
        End Select
    End Sub
    Public Shared Sub AssignPrimitiveType(ByRef prim As sPrimitive)
        Select Case prim.type
            Case PrimitiveType.SPHERE_TYPE
                prim.sphereProperties = Nothing
            Case PrimitiveType.PLANE_TYPE
                prim.planeProperties = Nothing
            Case PrimitiveType.CYLINDER_TYPE
                prim.cylinderProperties = Nothing
            Case PrimitiveType.TRIANGLE_TYPE
                prim.triangleProperties = Nothing
        End Select
    End Sub

    Public Shared Function IntersectPrimitive(ByRef prim As sPrimitive, ByRef rray As Ray.sRay) As Ray.TraceResult
        Select Case prim.type
            Case PrimitiveType.SPHERE_TYPE
                Return Sphere.IntersectSphere(prim.sphereProperties, rray)
            Case PrimitiveType.PLANE_TYPE
                Return Plane.IntersectPlane(prim.planeProperties, rray)
            Case PrimitiveType.CYLINDER_TYPE
                Return Cylinder.IntersectCylinder(prim.cylinderProperties, rray)
            Case PrimitiveType.TRIANGLE_TYPE
                Return Triangle.IntersectTriangle(prim.triangleProperties, rray)
        End Select
    End Function
    Public Shared Function CalculateNormal(ByRef prim As sPrimitive, ByRef intersection As Vector.sVector) As Vector.sVector
        Select Case prim.type
            Case PrimitiveType.SPHERE_TYPE
                Return Sphere.SphereNormal(prim.sphereProperties, intersection)
            Case PrimitiveType.PLANE_TYPE
                Return Plane.PlaneNormal(prim.planeProperties)
            Case PrimitiveType.CYLINDER_TYPE
                Return Cylinder.CylinderNormal(prim.cylinderProperties, intersection)
            Case PrimitiveType.TRIANGLE_TYPE
                Return Triangle.TriangleNormal(prim.triangleProperties)
        End Select
    End Function
End Class
Public Class Vector
    Public Structure sVector
        Dim x, y, z As Double
    End Structure

    Public Shared Sub VectorSetXYZ(ByRef v As sVector, ByVal x As Double, ByVal y As Double, ByVal z As Double)
        v.x = x
        v.y = y
        v.z = z
    End Sub

    Public Shared Function VectorAdd(ByVal a As sVector, ByVal b As sVector) As sVector ' result = a + b
        Dim result As sVector
        result.x = a.x + b.x
        result.y = a.y + b.y
        result.z = a.z + b.z
        Return result
    End Function

    Public Shared Function VectorSub(ByVal a As sVector, ByVal b As sVector) As sVector ' result = a - b
        Dim result As sVector
        result.x = a.x - b.x
        result.y = a.y - b.y
        result.z = a.z - b.z
        Return result
    End Function

    Public Shared Function VectorDot(ByVal a As sVector, ByVal b As sVector) As Double
        Return (a.x * b.x) + (a.y * b.y) + (a.z * b.z)
    End Function

    Public Shared Function VectorCross(ByVal a As sVector, ByVal b As sVector) As sVector
        Dim c As sVector
        c.x = a.y * b.z - a.z * b.y
        c.y = a.z * b.x - a.x * b.z
        c.z = a.x * b.y - a.y * b.x
        Return c
    End Function

    Public Shared Sub VectorNormalize(ByRef v As sVector)
        Dim scaleFactor As Double = 1 / Math.Sqrt((v.x * v.x) + (v.y * v.y) + (v.z * v.z))
        v.x *= scaleFactor
        v.y *= scaleFactor
        v.z *= scaleFactor
    End Sub

    Public Shared Sub Rotate(ByRef v As sVector, ByVal ax As Double, ByVal ay As Double, ByVal az As Double)
        Dim temp As sVector

        temp.y = v.y
        v.y = v.y * Math.Cos(ax) - v.z * Math.Sin(ax)
        v.z = v.z * Math.Cos(ax) + temp.y * Math.Sin(ax)

        temp.z = v.z
        v.z = v.z * Math.Cos(ay) - v.x * Math.Sin(ay)
        v.x = v.x * Math.Cos(ay) + temp.z * Math.Sin(ay)

        temp.x = v.x
        v.x = v.x * Math.Cos(az) - v.y * Math.Sin(az)
        v.y = v.y * Math.Cos(az) + temp.x * Math.Sin(az)
    End Sub

End Class
Public Class Matrix
    Public Class Matrix3x3
        Public elements(,) As Double
        Sub New()
            ReDim elements(2, 2)
        End Sub
    End Class


    ' multiples a vector by a 3x3 matrix to return a vector
    ' the vectors are treated as a 3x1 (meaning it has 3 columns and 1 row) matrix
    '                           [b00 b10 b20]
    ' [rx ry rz] = [ax ay az] * [b01 b11 b21]
    '                          [b02 b12 b22]
    Public Shared Function VectorMultMatrix(ByVal a As Vector.sVector, ByVal b As Matrix3x3) As Vector.sVector 'returns a*b. vectors are treated as matrices([x y z])
        Dim result As Vector.sVector
        result.x = a.x * b.elements(0, 0) + a.y * b.elements(0, 1) + a.z * b.elements(0, 2)
        result.y = a.x * b.elements(1, 0) + a.y * b.elements(1, 1) + a.z * b.elements(1, 2)
        result.z = a.x * b.elements(2, 0) + a.y * b.elements(2, 1) + a.z * b.elements(2, 2)
        Return result
    End Function

    Public Shared Function MatrixMult3x3(ByVal a As Matrix3x3, ByVal b As Matrix3x3) As Matrix3x3 'returns ab
        Dim result As Matrix3x3 = New Matrix3x3
        Dim i, j As Integer
        For i = 0 To 2
            For j = 0 To 2
                result.elements(i, j) = a.elements(0, j) * b.elements(i, 0) + a.elements(i, j) * b.elements(i, 1) + a.elements(2, j) * b.elements(i, 2)
            Next
        Next
        Return result
    End Function

    Public Shared Function CalculateArbitrayRotationMatrix(ByVal axis As Vector.sVector, ByVal theta As Double) As Matrix3x3
        Dim r As Matrix3x3 = New Matrix3x3
        Dim c As Double = Math.Cos(theta), s = Math.Sin(theta), t = 1 - Math.Cos(theta)
        r.elements(0, 0) = t * axis.x * axis.x + c
        r.elements(1, 0) = t * axis.x * axis.y - s * axis.z
        r.elements(2, 0) = t * axis.x * axis.z + s * axis.y
        r.elements(0, 1) = t * axis.x * axis.y + s * axis.z
        r.elements(1, 1) = t * axis.y * axis.y + c
        r.elements(2, 1) = t * axis.y * axis.z - s * axis.x
        r.elements(0, 2) = t * axis.x * axis.z - s * axis.y
        r.elements(1, 2) = t * axis.y * axis.z + s * axis.x
        r.elements(2, 2) = t * axis.z * axis.z + c
        Return r
    End Function

    Public Shared Function IdentityMatrix3x3() As Matrix3x3
        Dim i As Matrix3x3 = New Matrix3x3
        i.elements(0, 0) = 1
        i.elements(1, 0) = 0
        i.elements(2, 0) = 0
        i.elements(0, 1) = 0
        i.elements(1, 1) = 1
        i.elements(2, 1) = 0
        i.elements(0, 2) = 0
        i.elements(1, 2) = 0
        i.elements(2, 2) = 1
        Return i
    End Function
End Class
Public Class Triangle
    Public Structure TriangleProperties ' to optimize, precalculate normal
        Dim v1, v2, v3 As Vector.sVector ' vertices
    End Structure

    Public Shared Function IntersectTriangle(ByVal tri As TriangleProperties, ByVal pray As Ray.sRay) As Ray.TraceResult

        Dim tr As Ray.TraceResult
        Dim u, v As Double
        Dim edge1, edge2, tvec, pvec, qvec As Vector.sVector
        Dim det, invDet As Double

        edge1 = Vector.VectorSub(tri.v2, tri.v1)
        edge2 = Vector.VectorSub(tri.v3, tri.v1)
        pvec = Vector.VectorCross(pray.direction, edge2)

        det = Vector.VectorDot(edge1, pvec)

        If det > -0.00000099999999999999995 And det < 0.00000099999999999999995 Then
            tr.hit = False
            Return tr
        End If
        invDet = 1 / det

        tvec = Vector.VectorSub(pray.origin, tri.v1)

        u = Vector.VectorDot(tvec, pvec) * invDet
        If u < 0 Or u > 1 Then
            tr.hit = False
            Return tr
        End If

        qvec = Vector.VectorCross(tvec, edge1)

        v = Vector.VectorDot(pray.direction, qvec) * invDet
        If v < 0 Or (u + v) > 1 Then
            tr.hit = False
            Return tr
        End If

        tr.distance = Vector.VectorDot(edge2, qvec) * invDet
        If tr.distance < 0 Then
            tr.hit = False
            Return tr
        End If

        tr.hit = True
        Return tr
    End Function

    Public Shared Function TriangleNormal(ByVal tr As TriangleProperties) As Vector.sVector
        Dim edge1 As Vector.sVector = Vector.VectorSub(tr.v2, tr.v1), edge2 = Vector.VectorSub(tr.v3, tr.v1)
        Dim normal As Vector.sVector = Vector.VectorCross(edge1, edge2)
        Vector.VectorNormalize(normal)
        Return normal
    End Function
End Class
Public Class Plane
    Public Structure PlaneProperties
        Dim normal As Vector.sVector
        Dim displacement As Double
    End Structure

    Public Shared Function PlaneNormal(ByVal pl As PlaneProperties) As Vector.sVector
        Return pl.normal
    End Function

    Public Shared Function IntersectPlane(ByVal pl As PlaneProperties, ByVal pray As Ray.sRay) As Ray.TraceResult
        Dim tr As Ray.TraceResult
        Dim t As Double = -(pl.normal.x * pray.origin.x + pl.normal.y * pray.origin.y + pl.normal.z * pray.origin.z + pl.displacement) / _
            (pl.normal.x * pray.direction.x + pl.normal.y * pray.direction.y + pl.normal.z * pray.direction.z)
        If t < 0 Then
            tr.hit = False
            Return tr
        End If
        tr.hit = True
        tr.distance = t
        Return tr
    End Function
End Class
Public Class Sphere
    Public Structure SphereProperties
        Dim center As Vector.sVector
        Dim radius As Double
    End Structure

    Public Shared Function IntersectSphere(ByVal sp As SphereProperties, ByVal pray As Ray.sRay) As Ray.TraceResult
        Dim tr As Ray.TraceResult

        Dim rayToSphereCenter As Vector.sVector = Vector.VectorSub(sp.center, pray.origin)
        Dim lengthRTSC2 As Double = Vector.VectorDot(rayToSphereCenter, rayToSphereCenter) ' lengthRTSC2 = length of the ray from the ray's origin to the sphere's center squared

        Dim closestApproach As Double = Vector.VectorDot(rayToSphereCenter, pray.direction)
        If closestApproach < 0 Then ' "intersection" on säteen takana
            tr.hit = False
            Return tr
        End If

        ' halfCord2 = the distance squared from the closest approach of the ray to a perpendicular to the ray through the center of the sphere to the place where the ray actually intersects the sphere
        Dim halfCord2 As Double = (sp.radius * sp.radius) - lengthRTSC2 + (closestApproach * closestApproach)
        ' sphere.radius * sphere.radius could be precalced, but it might take longer to load it
        ' than to calculate it
        If halfCord2 < 0 Then ' the ray misses the sphere
            tr.hit = False
            Return tr
        End If

        tr.hit = True
        tr.distance = closestApproach - Math.Sqrt(halfCord2)
        Return tr
    End Function

    Public Shared Function SphereNormal(ByVal sp As SphereProperties, ByVal intersection As Vector.sVector) As Vector.sVector
        Dim normal As Vector.sVector
        'Laske normaali "intersectionin" osoittaamaan pisteeseen
        Dim oneOverRadius As Double = 1 / sp.radius
        normal.x = (intersection.x - sp.center.x) * oneOverRadius   ' sama kuin ( intersection.x - sphere.center.x ) / sphere.radius
        normal.y = (intersection.y - sp.center.y) * oneOverRadius
        normal.z = (intersection.z - sp.center.z) * oneOverRadius
        Return normal
    End Function

End Class
Public Class Cylinder
    Public Enum eAxis
        X_INFINITE
        Y_INFINITE
        Z_INFINITE
    End Enum
    Public Structure CylinderProperties
        Dim center As Vector.sVector
        Dim radius As Double
        Dim axis As Integer
    End Structure

    Public Shared Function IntersectCylinder(ByVal cy As CylinderProperties, ByVal pray As Ray.sRay) As Ray.TraceResult

        Dim tr As Ray.TraceResult
        Dim a, b, c As Double
        Select Case cy.axis
            Case eAxis.X_INFINITE
                a = pray.direction.y * pray.direction.y + pray.direction.z * pray.direction.z

                b = 2 * (pray.direction.y * (pray.origin.y - cy.center.y) + _
                        pray.direction.z * (pray.origin.z - cy.center.z))

                c = (pray.origin.y - cy.center.y) * (pray.origin.y - cy.center.y) + _
                        (pray.origin.z - cy.center.z) * (pray.origin.z - cy.center.z) - cy.radius * cy.radius
            Case eAxis.Y_INFINITE
                a = pray.direction.x * pray.direction.x + _
                        pray.direction.z * pray.direction.z

                b = 2 * (pray.direction.x * (pray.origin.x - cy.center.x) + _
                     pray.direction.z * (pray.origin.z - cy.center.z))

                c = (pray.origin.x - cy.center.x) * (pray.origin.x - cy.center.x) + _
                (pray.origin.z - cy.center.z) * (pray.origin.z - cy.center.z) - cy.radius _
                * cy.radius
            Case eAxis.Z_INFINITE
                a = pray.direction.x * pray.direction.x + _
                        pray.direction.y * pray.direction.y

                b = 2 * (pray.direction.x * (pray.origin.x - cy.center.x) + _
                        pray.direction.y * (pray.origin.y - cy.center.y))

                c = (pray.origin.x - cy.center.x) * (pray.origin.x - cy.center.x) + _
                    (pray.origin.y - cy.center.y) * (pray.origin.y - cy.center.y) - cy.radius * cy.radius
        End Select

        Dim discriminant As Double = b * b - 4 * a * c
        If discriminant < 0 Then
            tr.hit = False
            Return tr
        End If

        tr.distance = (-b - Math.Sqrt(discriminant)) / (2 * a)
        If tr.distance < 0 Then
            tr.hit = False
            Return tr
        End If
        tr.hit = True
        Return tr
    End Function
    Public Shared Function CylinderNormal(ByVal cy As CylinderProperties, ByVal intersection As Vector.sVector) As Vector.sVector
        Dim normal As Vector.sVector
        Dim oneOverRadius As Double = 1 / cy.radius
        Select Case cy.axis
            Case eAxis.X_INFINITE
                normal.x = 0
                normal.y = (intersection.y - cy.center.y) * oneOverRadius
                normal.z = (intersection.z - cy.center.z) * oneOverRadius
            Case eAxis.Y_INFINITE
                normal.x = (intersection.x - cy.center.x) * oneOverRadius
                normal.y = 0
                normal.z = (intersection.z - cy.center.z) * oneOverRadius
            Case eAxis.Z_INFINITE
                normal.x = (intersection.x - cy.center.x) * oneOverRadius
                normal.y = (intersection.y - cy.center.y) * oneOverRadius
                normal.z = 0
        End Select
        Return normal
    End Function
End Class