Rndr2024.julia

# myabs(x::Float64) = sign(x) * x

import Base.:-
import Base.:+

# convert from 0 based index to julia 1 based index
from0basedIdxToJulia(idx::Int) = idx + 1

pow(v, exponent) = v ^ exponent

struct Vec3
    x::Float32
    y::Float32
    z::Float32
end

-(v::Vec3) = Vec3(-v.x,-v.y,-v.z)
calcMag(v::Vec3) = sqrt(v.x*v.x + v.y*v.y + v.z*v.z)
dot(a::Vec3, b::Vec3) = a.x*b.x + a.y*b.y + a.z*b.z
scale(v::Vec3, s::Float32) = Vec3(v.x*s,v.y*s,v.z*s)
-(a::Vec3, b::Vec3) = Vec3(a.x-b.x,a.y-b.y,a.z-b.z)
+(a::Vec3, b::Vec3) = Vec3(a.x+b.x,a.y+b.y,a.z+b.z)
calcMagSquare(v::Vec3) = dot(v, v)

vecAbs(v::Vec3) = Vec3(abs(v.x), abs(v.y), abs(v.z))

function normalize(v::Vec3)::Vec3
    mag::Float32 = calcMag(v)
    return scale(v, 1.0f0/mag)
end

function cross(a::Vec3, b::Vec3)::Vec3
    x = a.y*b.z - a.z*b.y
    y = a.z*b.x - a.x*b.z
    z = a.x*b.y - a.y*b.x
    return Vec3(x, y, z)
end

# /param d direction
# /param n normal, must be normalized
function reflect(d::Vec3, n::Vec3)::Vec3
    # see https://math.stackexchange.com/questions/13261/how-to-get-a-reflection-vector
    return d - scale(n, 2.0f0*dot(d, n))
end

function maxComponentIndex(v::Vec3)::Int
    if v.x > v.y && v.x > v.z
        return 0
    elseif v.y > v.x && v.y > v.z
        return 1
    else
        return 2
    end
end

function vecPermutate(v::Vec3, x::Int, y::Int, z::Int)::Vec3
    xValue = v.x
    if x == 1
        xValue = v.y
    elseif x == 2
        xValue = v.z
    end

    yValue = v.x
    if y == 1
        yValue = v.y
    elseif y == 2
        yValue = v.z
    end

    zValue = v.x
    if z == 1
        zValue = v.y
    elseif z == 2
        zValue = v.z
    end

    return Vec3(xValue, yValue, zValue)
end

# return component by zero based index
function retComponentByIndex(v::Vec3, idx::Int)::Float32
    if idx == 0
        return v.x
    elseif idx == 1
        return v.y
    else
        return v.z
    end
end

struct Ray
    origin::Vec3
    direction::Vec3
    direction_inv::Vec3
    tMax::Float32
end

function ray__make(origin::Vec3, direction::Vec3, tMax::Float32)::Ray
    direction_inv::Vec3 = Vec3(1.0f0/direction.x, 1.0f0/direction.y, 1.0f0/direction.z)
    return Ray(origin, direction, direction_inv, tMax)
end

struct Sphere
    center::Vec3
    radius::Float32
end

function intersectRayVsSphere(ray::Ray, sphere::Sphere)
    a = dot(ray.direction, ray.direction)
    b = 2.0*dot(ray.direction, ray.origin - sphere.center)
    c = calcMagSquare(ray.origin - sphere.center) - sphere.radius*sphere.radius

    #println(a) # DBG
    #println(b) # DBG
    #println(c) # DBG

    inSqrt = b*b - 4.0f0*a*c
    if inSqrt < 0.0f0
        return nothing
    end

    sqrtRes = sqrt(inSqrt)
    ta = (-b - sqrtRes) / (2.0f0*a)
    tb = (-b + sqrtRes) / (2.0f0*a)

    return (ta, tb)
end






# manual test
if false
    rayA = ray__make(Vec3(0.0, 0.0, -5.0), Vec3(0.0, 0.0, 1.0), 1e20)

    sphereA = Sphere(Vec3(0.0, 0.0, 0.0), 1.0)

    res = intersectRayVsSphere(rayA, sphereA)

    println(res)

    if res === nothing
        println("res is nothing")
    else
        println("res is not nothing")
    end
end








struct Quaternion
    x::Float64
    y::Float64
    z::Float64
    w::Float64
end

function makeQuaternionId()::Quaternion
    return Quaternion(0.0, 0.0, 0.0, 1.0)
end

function makeQuaternionFromEuler(axis::Vec3, angle::Float64)
    # see https://en.wikipedia.org/wiki/Quaternions_and_spatial_rotation
    c::Float64 = cos(angle / 2.0)
    S::Float64 = sin(angle / 2.0)
    return Quaternion(axis.x*S, axis.y*S, axis.z*S, c)
end

# expects normalized Quaternion
function convQuaternionTo44Matrix(q::Quaternion)
    # see https://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm
    return [
        1.0 - 2.0*q.y*q.y - 2.0*q.z*q.z   2.0*q.x*q.y - 2.0*q.z*q.w         2.0*q.x*q.z + 2.0*q.y*q.w         0.0;
        2.0*q.x*q.y + 2.0*q.z*q.w  	      1.0 - 2.0*q.x*q.x - 2.0*q.z*q.z   2.0*q.y*q.z - 2.0*q.x*q.w         0.0;
        2.0*q.x*q.z - 2.0*q.y*q.w         2.0*q.y*q.z + 2.0*q.x*q.w         1.0 - 2.0*q.x*q.x - 2.0*q.y*q.y   0.0;
        0.0                               0.0                               0.0                               1.0
    ]
end

function matrixMul(m::Matrix{Float64}, v::Vec3)::Vec3
    resultM = m * [v.x; v.y; v.z; 1.0]
    return Vec3(resultM[1], resultM[2], resultM[3])
end

if false
    v = Vec3(1.0, 0.5, 0.5)

    q = makeQuaternionFromEuler(Vec3(1.0, 0.0, 0.0), 0.3)
    m = convQuaternionTo44Matrix(q)
    v = matrixMul(m, v)

    print(v)

    exit(0)
end






mutable struct Camera
    pos::Vec3
    dir::Vec3
    up::Vec3
    side::Vec3
end

# x and y between 0.0 and 1.0 inclusive
# result is not normalized
function cam__calcRayDir(cam::Camera, x::Float32, y::Float32)::Vec3
    rel::Vec3 = scale(cam.side, x*2.0f0-1.0f0) + scale(cam.up, y*2.0f0-1.0f0)
    return cam.dir + rel
end

function cam__calcAbsolutePosition(cam::Camera, dir::Vec3, t::Float32)::Vec3
    return cam.pos + scale(dir, t)
end

# manual test
if false
    cam = Camera(Vec3(0.0, 0.0, 0.0), Vec3(0.0, 0.0, 1.0), Vec3(0.0, 1.0, 0.0), Vec3(1.0, 0.0, 0.0))

    rayDirA::Vec3 = cam__calcRayDir(cam, 0.5f0, 0.5f0)

    println(rayDirA)
end



# manual test of simple rendering of "image" of sphere to test camera calculations
if false


    sizeX = 20
    sizeY = 20
    

    cam = Camera(Vec3(0.0, 0.0, 0.0), Vec3(0.0, 0.0, 1.0), Vec3(0.0, 1.0, 0.0), Vec3(1.0, 0.0, 0.0))

    sphereA = Sphere(Vec3(0.0, 0.0, 3.0), 1.0)



    arr = Float32[]

    
    for iy in 0:(sizeY-1)
        for ix in 0:(sizeX-1)

            rayDirA::Vec3 = cam__calcRayDir(cam, Float32(ix / (sizeX-1)), Float32(iy / (sizeY-1)))        

            rayA = ray__make(cam.pos, rayDirA, 1e20)

            res = intersectRayVsSphere(rayA, sphereA)


            if res === nothing
                push!(arr, -1.0f0)
            else
                push!(arr, res[0+1])
            end

            # DEBUG
            if res === nothing
                print(" ")
            else
                print("x")
            end
        end

        # DEBUG
        println("")
    end


end








# struct to store color
struct ColorRgb_64bit
    r::Float64
    g::Float64
    b::Float64
end

+(a::ColorRgb_64bit, b::ColorRgb_64bit) = ColorRgb_64bit(a.r+b.r,a.g+b.g,a.b+b.b)
scale(arg::ColorRgb_64bit, s::Float64) = ColorRgb_64bit(arg.r*s, arg.g*s, arg.b*s)

function colorRgb__mulComponents(a::ColorRgb_64bit, b::ColorRgb_64bit)::ColorRgb_64bit
    return ColorRgb_64bit(a.r*b.r, a.g*b.g, a.b*b.b)
end




# shading
function shading__shadeLambertian(vIn::Vec3, vOut::Vec3)::Float32
    h0 = dot(vIn, vOut)
    h1 = max(h0, 0.0f0)
    return h1
end









mutable struct RayInfo
    closestHitT::Float64
    isAnyHit::Bool
    surfaceNormal::Vec3

    # TODO : pull from Material over surfaceMaterialIdx
    surfaceBaseColor::ColorRgb_64bit
    # TODO : pull from Material over surfaceMaterialIdx
    surfaceEmission::ColorRgb_64bit
    # TODO : pull from Material over surfaceMaterialIdx
    surfaceAlbedoColor::ColorRgb_64bit

    surfaceMaterialIdx::Int64 # 0-based

    hitInstanceIdx::Int64 # 0-based
end

function rayInfo__make():RayInfo
    return RayInfo(1e20, false, Vec3(1.0,0.0,0.0),ColorRgb_64bit(0.0,0.0,0.0),ColorRgb_64bit(0.0,0.0,0.0),ColorRgb_64bit(0.0,0.0,0.0),-1,-1)
end


# manual test of rendering to see if basic rendering loop works
if false


    sizeX = 20
    sizeY = 20
    

    cam = Camera(Vec3(0.0, 0.0, 0.0), Vec3(0.0, 0.0, 1.0), Vec3(0.0, 1.0, 0.0), Vec3(1.0, 0.0, 0.0))

    sphereA = Sphere(Vec3(0.0, 0.0, 3.0), 1.0)



    arr::Array{ColorRgb_64bit} = ColorRgb_64bit[]

    
    for iy in 0:(sizeY-1)
        for ix in 0:(sizeX-1)

            rayInfo::RayInfo = rayInfo__make()

            rayDirA::Vec3 = cam__calcRayDir(cam, Float32(ix / (sizeX-1)), Float32(iy / (sizeY-1)))        

            rayA = ray__make(cam.pos, rayDirA, 1e20)



            raySphereIntersectionRes = intersectRayVsSphere(rayA, sphereA)

            if !(raySphereIntersectionRes === nothing)
                t::Float32 = raySphereIntersectionRes[0+1] # t of first intersection of ray with sphere

                if t >= 0.0f0 && t < rayInfo.closestHitT && t < rayA.tMax
                    rayInfo.closestHitT = t
                    rayInfo.isAnyHit = true

                    posWorld::Vec3 = cam__calcAbsolutePosition(cam, rayDirA, t)
                    
                    normal::Vec3 = scale(posWorld-sphereA.center, 1.0f0/sphereA.radius)
                    rayInfo.surfaceNormal = normal
                    rayInfo.surfaceBaseColor = ColorRgb_64bit(1.0, 1.0, 1.0)
                end
            end
            


            resColor::ColorRgb_64bit = ColorRgb_64bit(0.0, 0.0, 0.0)

            if rayInfo.isAnyHit

                # TODO  :  compute direction to light
                dirToLight::Vec3 = Vec3(0.0f0, 0.0f0, -1.0f0)

                shade::Float32 = shading__shadeLambertian(dirToLight, rayInfo.surfaceNormal)

                lightColor::ColorRgb_64bit = ColorRgb_64bit(1.0, 1.0, 1.0)
                helperColorA::ColorRgb_64bit = colorRgb__mulComponents(rayInfo.surfaceBaseColor, lightColor)

                resColor = colorRgb__mulComponents(ColorRgb_64bit(shade, shade, shade), helperColorA)
            end


            push!(arr, resColor)

        end
    end


    # DEBUG IMAGE
    for iColor in arr
        println(iColor)
    end
end


# helper
function ray__calcAbsolutePosition(originPos::Vec3, dir::Vec3, t::Float32)::Vec3
    return originPos + scale(dir, t)
end









# helper
function FMA(a::Float32, b::Float32, c::Float32)::Float32
    return a * b + c
end

# helper
function differenceOfProducts(a::Float32, b::Float32, c::Float32, d::Float32)::Float32
    cd::Float32 = c * d;
    differenceOfProducts2::Float32 = FMA(a, b, -cd);
    error::Float32 = FMA(-c, d, cd);
    return differenceOfProducts2 + error;
end



struct TriangleA
    p0::Vec3
    p1::Vec3
    p2::Vec3
end

struct TriangleIntersectionA
    b0::Float32
    b1::Float32
    b2::Float32
    t::Float32
end

# ray vs triangle intersection

function intersectRayVsTriangleA(ray::Ray, triangle::TriangleA)::Union{Nothing, TriangleIntersectionA}

    # properties
    #  * watertight intersection-test

    # see PBRT book

    # Ü Return no intersection if triangle is degenerate
    if calcMagSquare(cross(triangle.p2 - triangle.p0, triangle.p1 - triangle.p0)) == 0.0f0
        return nothing
    end

    # Ü translate vertices
    p0t::Vec3 = triangle.p0 - ray.origin
    p1t::Vec3 = triangle.p1 - ray.origin
    p2t::Vec3 = triangle.p2 - ray.origin
    
    # Ü permutate components of triangle vertices and ray direction
    kz::Int = maxComponentIndex(vecAbs(ray.direction))
    kx::Int = kz + 1
    if kx == 3
        kx = 0
    end
    ky::Int = kx + 1
    if ky == 3
        ky = 0
    end
    d::Vec3 = vecPermutate(ray.direction, kx, ky, kz)
    p0t = vecPermutate(p0t, kx, ky, kz)
    p1t = vecPermutate(p1t, kx, ky, kz)
    p2t = vecPermutate(p2t, kx, ky, kz)

    # Ü apply shear transformation to translated vertex positions
    Sx::Float32 = -ray.direction.x / ray.direction.z
    Sy::Float32 = -ray.direction.y / ray.direction.z
    Sz::Float32 = 1.0f0 / ray.direction.z
    p0t = Vec3(p0t.x + Sx * p0t.z, p0t.y + Sy * p0t.z, p0t.z)
    p1t = Vec3(p1t.x + Sx * p1t.z, p1t.y + Sy * p1t.z, p1t.z)
    p2t = Vec3(p2t.x + Sx * p2t.z, p2t.y + Sy * p2t.z, p2t.z)
    

    # now we can do the actual intersection


    # Ü compute edge function coefficients e0, e1, e2
    e0::Float32 = differenceOfProducts(p1t.x, p2t.y, p1t.y, p2t.x)
    e1::Float32 = differenceOfProducts(p2t.x, p0t.y, p2t.y, p0t.x)
    e2::Float32 = differenceOfProducts(p0t.x, p1t.y, p0t.y, p1t.x)

    # TODO LOW  :  implement the same using double precision and check if we need it.
    #              


    # Ü perform triangle edge and determinant tests
    if (e0 < 0.0f0 || e1 < 0.0f0 || e2 < 0.0f0) && (e0 > 0.0f0 || e1 > 0.0f0 || e2 > 0.0f0)
        #println("exit D")
        return nothing
    end
    det::Float32 = e0 + e1 + e2
    if det == 0.0f0
        #println("exit C")
        return nothing
    end

    # Ü compute scaled hit distance to triangle and test against ray t range
    p0t = Vec3(p0t.x, p0t.y, p0t.z * Sz)
    p1t = Vec3(p1t.x, p1t.y, p1t.z * Sz)
    p2t = Vec3(p2t.x, p2t.y, p2t.z * Sz)
    tScaled::Float32 = e0 * p0t.z + e1 * p1t.z + e2 * p2t.z
    #println(det)
    #println(tScaled)
    if det < 0.0f0 && (tScaled >= 0.0f0 || tScaled < ray.tMax * det)
        #println(det < 0)
        #println(tScaled >= 0 || tScaled < ray.tMax * det)
        #println("exit E")
        return nothing
    elseif det > 0.0f0 && (tScaled <= 0.0f0 || tScaled > ray.tMax * det)
        #println("exit F")
        return nothing
    end

    # now we have a valid intersection

    # Ü compute barycentric coordinates and t value for triangle intersection
    invDet::Float32 = 1.0f0 / det
    b0::Float32 = e0 * invDet
    b1::Float32 = e1 * invDet
    b2::Float32 = e2 * invDet
    t::Float32 = tScaled * invDet

    # TODO LOW  :   <<Ensure that computed triangle is conservatively greater than zero>> 

    return TriangleIntersectionA(b0, b1, b2, t)
end

# manual test to see if ray vs triangle intersection works
if false

    triangle0::TriangleA = TriangleA(Vec3(5.0f0, 5.0f0, 5.01f0), Vec3(5.0f0, -5.0f0, 5.02f0), Vec3(-5.0f0, 0.0f0, 5.03f0))

    #ray0::Ray = ray__make(Vec3(0.001f0, 0.000222f0, 0.000002232f0), Vec3(0.0f0, 0.000001f0, 0.9999999f0), 1e20) # works

    rayDir::Vec3 = normalize(Vec3(0.052631617f0, 0.552631617f0, 1.0f0))
    rayDir = normalize(Vec3(0.0f0, 0.000001f0, 0.9999999f0))
    #rayDir = Vec3(-0.57735026f0, -0.57735026f0, 0.57735026f0)
    
    if true # codeblock
        ray0::Ray = ray__make(Vec3(0.001f0, 0.000222f0, 0.000002232f0), rayDir, 1e20)

        rayVsTriangleRes = intersectRayVsTriangleA(ray0, triangle0)

        println(rayVsTriangleRes)
    end

    if true # codeblock
        ray0::Ray = ray__make(Vec3(0.001f0, 3.000222f0, 0.000002232f0), rayDir, 1e20)

        rayVsTriangleRes = intersectRayVsTriangleA(ray0, triangle0)

        println(rayVsTriangleRes)
    end

    #exit()
end















function calcBrdf__ABg(scatteredRay::Vec3, normal::Vec3, outgoingDir::Vec3, a::Float32, b::Float32, g::Float32)::Float32
    # see https://gist.github.com/PtrMan/18b46321f5cfb6ff3e26413c5d8a541d
    
    specularRay::Vec3 = reflect(normalize(outgoingDir), normal) # outgoing specular ray
    tangent::Vec3 = normalize(cross(normal, specularRay))
    ortogonal::Vec3 = normalize(cross(normal, tangent))

    specularProjected::Vec3 = Vec3(dot(specularRay, tangent), dot(specularRay, ortogonal), 0.0f0) # projected specular ray
    scatteredProjected::Vec3 = Vec3(dot(scatteredRay, tangent), dot(scatteredRay, ortogonal), 0.0f0) # projected specular ray

    #float a = 0.001f;
    #float b = 0.0025f;
    #float g = 1.5;

    projectedVecDistance::Float32 = calcMag(specularProjected-scatteredProjected)
    #projectedVecDistance = 0.005f



    res::Float32 = b / pow(projectedVecDistance, g)
    return res
end

# manual test for ABg-BRDF
if false
    normal::Vec3 = Vec3(1.0f0, 0.0000001f0, 0.0f0)
    outgoingDir::Vec3 = Vec3(0.999999f0, 0.0f0, 0.0f0)
    scatteredRay::Vec3 = Vec3(-1.0f0, 0.0f0, 0.0f0)

    a = 0.01f0
    b = 0.000000025f0
    g = 1.5f0
    brdfResult::Float32 = calcBrdf__ABg(scatteredRay, normal, outgoingDir, a, b, g)
    println(brdfResult)
end















struct MeshTriangle
    vertexIdx0::Int # zero based index
    vertexIdx1::Int # zero based index
    vertexIdx2::Int # zero based index

    #normalIdx0::Int # zero based index into normals array
    #normalIdx1::Int # zero based index into normals array
    #normalIdx2::Int # zero based index into normals array
end

mutable struct Mesh
    vertices::Array{Vec3}
    # normals::Array{Vec3}
    triangles::Array{MeshTriangle}
end


mutable struct Instance
    mesh::Mesh

    worldPosition::Vec3

    rotation::Quaternion

    # TODO : move into material
    emission::ColorRgb_64bit
    # TODO : move into material
    albedoColor::ColorRgb_64bit # color of the surface

    materialIdx::Int64 # index to material of the surface of this instance . 0-based
end

struct ExplicitLight
    pos::Vec3
    intensity::Float32
end

mutable struct MaterialA
    # index for the used BRDF of this material
    # 0 : lambertian BRDF
    # 1 : ABg BRDF
    brdfType::Int64
end

mutable struct Scene
    camera::Camera
    instances::Array{Instance}
    lights::Array{ExplicitLight}

    materials::Array{MaterialA}

    # for sampling the light sources which are instances with emissive material
    lightsourceInstanceIndices::Array{Int64}
end


#= commented because outdated and not used
function intersectRayVsScene(ray::Ray, rayInfo::RayInfo, scene::Scene)
    # inefficient: iterate over instances and polygons

    for iInstance in scene.instances
        for iMeshTriangle in iInstance.mesh.triangles

            vertexPos0::Vec3 = iInstance.worldPosition + iInstance.mesh.vertices[ iMeshTriangle.vertexIdx0+1 ]
            vertexPos1::Vec3 = iInstance.worldPosition + iInstance.mesh.vertices[ iMeshTriangle.vertexIdx1+1 ]
            vertexPos2::Vec3 = iInstance.worldPosition + iInstance.mesh.vertices[ iMeshTriangle.vertexIdx2+1 ]

            triangle0::TriangleA = TriangleA(vertexPos0, vertexPos1, vertexPos2)

            rayVsTriangleRes = intersectRayVsTriangleA(ray, triangle0)
            
            if !(rayVsTriangleRes === nothing)

                t::Float32 = rayVsTriangleRes.t

                if t >= 0.0f0 && t < rayInfo.closestHitT && t < ray.tMax
                    rayInfo.closestHitT = t
                    rayInfo.isAnyHit = true
                    
                    normal::Vec3 = cross(vertexPos1-vertexPos0, vertexPos2-vertexPos0)
                    rayInfo.surfaceNormal = normal
                    rayInfo.surfaceBaseColor = ColorRgb_64bit(1.0, 1.0, 1.0)
                    rayInfo.surfaceMaterialIdx = iInstance.materialIdx
                    rayInfo.hitInstanceIdx = -1 # TODO TODO TODO
                end
            end
        end
    end
end
=#









mutable struct Map2d
    arr::Array{ColorRgb_64bit}
    width::Int32
end

function map2d__make(width::Int32, height::Int32)::Map2d
    res::Map2d = Map2d([], width)
    for z in 1:width*height
        push!(res.arr, ColorRgb_64bit(0.0, 0.0, 0.0))
    end
    return res
end

function map2d__readAt(y::Int64, x::Int64, map2d::Map2d)::ColorRgb_64bit
    return map2d.arr[1 + (x + y*map2d.width)]
end

function map2d__retWidth(map2d::Map2d)::Int32
    return map2d.width
end

function map2d__retHeight(map2d::Map2d)::Int32
    return size(map2d.arr)[1] / map2d.width
end

function storePpm(img::Map2d)
    s::String="P3\n"
    s *= string(map2d__retWidth(img))*" "*string(map2d__retHeight(img))*"\n"
    s *= "255\n" # maximum value
    
    for iy in 0:map2d__retHeight(img)-1
        for ix in 0:map2d__retWidth(img)-1
            valRgb::ColorRgb_64bit = map2d__readAt(iy, ix, img)
            
            r::Float32 = min(valRgb.r, 1.0)
            g::Float32 = min(valRgb.g, 1.0)
            b::Float32 = min(valRgb.b, 1.0)

            s *= string(trunc(Int, r*255.0f0))*"\n"*string(trunc(Int, g*255.0f0))*"\n"*string(trunc(Int, b*255.0f0))*"\n"
        end
    end
    
    open("a0.ppm", "w") do f
        write(f, s)
    end
end















function retLength(arr)::Int64
    return size(arr)[1]
end


# prototyping area for OBJ importer
function convertObjTextToMesh(str::String)::Mesh
    
    # see https://en.wikipedia.org/wiki/Wavefront_.obj_file

    vertices::Array{Vec3} = []
    verticeNormals::Array{Vec3} = []
    faces = []

    lines::Array{String} = split(str, '\n')
    for iLine in lines

        parts = split(iLine, ' ') # Ü  split by spaces
        type_ = parts[1+0]

        if type_ == "v" || type_ == "vn"

            # Ü  Convert the vector by parsing float
            floatArr = parse.(Float32, parts[1+1:end])

            # Output the results
            #println("Type: $type_")
            #println("floatArr: $floatArr")

            if type_ == "v"
                push!(vertices, Vec3(floatArr[0+1], floatArr[1+1], floatArr[2+1]))
            elseif type_ == "vn"
                push!(verticeNormals, Vec3(floatArr[0+1], floatArr[1+1], floatArr[2+1]))
            end
        elseif type_ == "f" # face

            vertexIndices = []
            textureIndices = []
            normalIndices = []

            # loop is necessary to handle more complicated format of face with "//" encoding
            for iPart in parts[1+1:end]
                indicesText = split(iPart, '/')

                parsedAsIntArr = []
                for iIndexText in indicesText
                    if iIndexText == ""
                        push!(parsedAsIntArr, nothing)
                    else
                        parsedIntValue = parse.(Int32, iIndexText)
                        push!(parsedAsIntArr, parsedIntValue)
                    end
                end

                #println(parsedAsIntArr) # DBG

                vertexIndex = parsedAsIntArr[0+1]-1 # vertex index : convert to zero based
                textureIndex = nothing
                if retLength(parsedAsIntArr) >= 2 && !(parsedAsIntArr[1+1] === nothing) # is there a texture index ?
                    textureIndex = parsedAsIntArr[1+1]-1 # convert to zero based
                end
                normalIndex = nothing
                if retLength(parsedAsIntArr) >= 3 && !(parsedAsIntArr[2+1] === nothing) # is there a normal index ? 
                    normalIndex = parsedAsIntArr[2+1]-1 # convert to zero based
                end

                push!(vertexIndices, vertexIndex)
                push!(textureIndices, textureIndex)
                push!(normalIndices, normalIndex)
            end

            #println(vertexIndices) # DBG
            #println(textureIndices) # DBG
            #println(normalIndices) # DBG

            # TODO LOW  :  handle faces with 4 vertices!

            meshTriangle::MeshTriangle = MeshTriangle(
                vertexIndices[0+1], vertexIndices[1+1], vertexIndices[2+1]
                # , normalIndices[0+1], normalIndices[1+1], normalIndices[2+1]
            )

            push!(faces, meshTriangle)
        else
            # ignore
        end
    end

    # * generate actual mesh from it
    resultMesh::Mesh = Mesh(vertices, faces)
    return resultMesh
end

# manual test for OBJ parsing
if false
    # The input string
    str::String = ""
    str = "v 5.3455 -3.4345 -24.5664\nv 3.2 2.1 1.23\nf 1//1 2//2 3//3\n" # face
    
    mesh::mesh = convertObjTextToMesh(str)
end


# function to read obj file and convert the text to actual mesh
function readObj(path::String)::Mesh
    str::String = ""
    open(path) do f
        str = read(f, String)
    end
    return convertObjTextToMesh(str)
end

if false
    path = "C:\\Users\\rober\\fsRoot\\art_blender\\geoA\\boxA.obj"
    meshA::Mesh = readObj(path)
end


# TODO LOW  :  mesh: support array of vertex normals, so that we can load the normals from the OBJ file for real!



















struct Aabb
    min::Vec3
    max::Vec3
end

# ray vs AABB test
# necessary for BVH
function rayVsAabbIntersection(box::Aabb, ray::Ray)::Bool
    # see https://tavianator.com/2015/ray_box_nan.html
    
    t1 = (retComponentByIndex(box.min, 0) - retComponentByIndex(ray.origin, 0)) * retComponentByIndex(ray.direction_inv, 0)
    t2 = (retComponentByIndex(box.max, 0) - retComponentByIndex(ray.origin, 0)) * retComponentByIndex(ray.direction_inv, 0)

    tmin = min(t1, t2)
    tmax = max(t1, t2)

    for i in 1:2
        t1 = (retComponentByIndex(box.min, i) - retComponentByIndex(ray.origin, i)) * retComponentByIndex(ray.direction_inv, i)
        t2 = (retComponentByIndex(box.max, i) - retComponentByIndex(ray.origin, i)) * retComponentByIndex(ray.direction_inv, i)

        tmin = max(tmin, min(min(t1, t2), tmax))
        tmax = min(tmax, max(max(t1, t2), tmin))
    end

    return tmax > max(tmin, 0.0)
end


struct BvhLeafGeometry
    # in world coordinates
    vertexPos0::Vec3
    vertexPos1::Vec3
    vertexPos2::Vec3

    # TODO : vertex normals for all 3 vertices!

    emission::ColorRgb_64bit # emission of whole polygon
    albedoColor::ColorRgb_64bit # color of the surface

    materialIdx::Int64 # index to material of the surface of this triangle . 0-based
    instanceIdx::Int64 # index of instance to which the traingle belongs to . 0-based
end

mutable struct BvhNodeA
    aabb::Aabb

    # children which are container by this Bvh node
    children::Array{BvhNodeA}

    # leaf geometry which are raw triangles with all data on it
    leafGeometry::Array{BvhLeafGeometry}
end

function isLeaf(bvhNode::BvhNodeA)::Bool
    return size(bvhNode.children)[1] == 0 # is a leaf when there are no children
end


function rayVsBvh(ray::Ray, rayInfo::RayInfo,  bvhNode::BvhNodeA)
    isHitBoundingBox::Bool = rayVsAabbIntersection(bvhNode.aabb, ray)
    if !isHitBoundingBox
        return
    end

    if isLeaf(bvhNode)

        # trace ray against leaf geometry

        for iTriangleLeafGeometry in bvhNode.leafGeometry

            triangle0::TriangleA = TriangleA(iTriangleLeafGeometry.vertexPos0, iTriangleLeafGeometry.vertexPos1, iTriangleLeafGeometry.vertexPos2)

            rayVsTriangleRes = intersectRayVsTriangleA(ray, triangle0)
            
            if !(rayVsTriangleRes === nothing)

                t::Float32 = rayVsTriangleRes.t

                if t >= 0.0f0 && t < rayInfo.closestHitT && t < ray.tMax
                    rayInfo.closestHitT = t
                    rayInfo.isAnyHit = true
                    
                    normal::Vec3 = cross(iTriangleLeafGeometry.vertexPos1-iTriangleLeafGeometry.vertexPos0, iTriangleLeafGeometry.vertexPos2-iTriangleLeafGeometry.vertexPos0)
                    rayInfo.surfaceNormal = normal
                    rayInfo.surfaceBaseColor = ColorRgb_64bit(1.0, 1.0, 1.0)
                    rayInfo.surfaceEmission = iTriangleLeafGeometry.emission
                    rayInfo.surfaceAlbedoColor = iTriangleLeafGeometry.albedoColor
                    rayInfo.surfaceMaterialIdx = iTriangleLeafGeometry.materialIdx
                    rayInfo.hitInstanceIdx = iTriangleLeafGeometry.instanceIdx
                end
            end
        end
    else
        # Ü  call recursive
        for iChildren in bvhNode.children
            rayVsBvh(ray, rayInfo, iChildren)
        end
    end
end







# PROTOTYPING AERA for BVH building

struct BvhGeometryWithCenter
    geo::BvhLeafGeometry
    center::Vec3
end


# build BVH without computing AABB
# /param splitDimIdx index of dimension which is used for splitting
function buildBvhRaw(toSegment::Array{BvhGeometryWithCenter}, splitDimIdx::Int64)::BvhNodeA
    

    if retLength(toSegment) <= 2
        # build leaf and return

        createdBvhNodeA::BvhNodeA = BvhNodeA(Aabb(Vec3(0.0f0, 0.0f0, 0.0f0), Vec3(0.0f0, 0.0f0, 0.0f0)), [], [])
        
        # Ü  add polygon geometry
        for iGeo in toSegment
            push!(createdBvhNodeA.leafGeometry, iGeo.geo)
        end

        return createdBvhNodeA
    end

    verticesAll::Array{Vec3} = []
    #push!(verticesAll, Vec3(0.0f0, 0.0f0, 0.0f0))
    #push!(verticesAll, Vec3(1.0f0, 1.0f0, 1.0f0))
    #push!(verticesAll, Vec3(2.0f0, 2.0f0, 2.0f0))
    #push!(verticesAll, Vec3(3.0f0, 3.0f0, 3.0f0))

    for iToSegment in toSegment
        # HACKY but quick
        push!(verticesAll, iToSegment.center)
    end


    #splitDimIdx::Int32 = 1

    sort!(verticesAll, by = iv -> retComponentByIndex(iv, splitDimIdx))

    # select middle vertex position. this is where we split the geometry
    middleVertex::Vec3 = verticesAll[1+trunc(Int, (retLength(verticesAll)/2))]

    #println(middleVertex) # DBG

    segmentLow::Array{BvhGeometryWithCenter} = []
    segmentHigh::Array{BvhGeometryWithCenter} = []


    for iToSegment in toSegment
        if retComponentByIndex(iToSegment.center, splitDimIdx) < retComponentByIndex(middleVertex, splitDimIdx)
            push!(segmentLow, iToSegment)
        else
            push!(segmentHigh, iToSegment)
        end
    end

    # HACKY : avoid unsplitted cases
    if retLength(segmentLow) == 0
        segmentLow2::Array{BvhGeometryWithCenter} = []
        segmentHigh2::Array{BvhGeometryWithCenter} = []

        push!(segmentLow2, segmentHigh[0+1])

        for idx0 in 1:retLength(segmentHigh)-1
            push!(segmentHigh2, segmentHigh[idx0+1])
        end
        
        segmentLow = segmentLow2
        segmentHigh = segmentHigh2
    elseif retLength(segmentHigh) == 0
        segmentLow2 = []
        segmentHigh2 = []

        push!(segmentHigh2, segmentLow[0+1])

        for idx0 in 1:retLength(segmentLow)-1
            push!(segmentLow2, segmentLow[idx0+1])
        end
        
        segmentLow = segmentLow2
        segmentHigh = segmentHigh2
    end


    # Ü  quick and dirty algorithm to decide the dimension of the split
    splitNextDimIdx::Int64 = mod(splitDimIdx+1, 3)

    nodeLow = buildBvhRaw(segmentLow, splitNextDimIdx)
    nodeHigh = buildBvhRaw(segmentHigh, splitNextDimIdx)

    # Ü build BVH node
    createdBvhNodeB::BvhNodeA = BvhNodeA(Aabb(Vec3(0.0f0, 0.0f0, 0.0f0), Vec3(0.0f0, 0.0f0, 0.0f0)), [nodeLow, nodeHigh], [])
    return createdBvhNodeB
end











# helper function to compute AABB of BvhLeafGeometry
function bvhLeafGeometry__calcAabb(geo::Array{BvhLeafGeometry})::Aabb
    min_::Vec3 = Vec3(Float32(1e20), Float32(1e20), Float32(1e20))
    max_::Vec3 = Vec3(Float32(-1e20), Float32(-1e20), Float32(-1e20))
    
    for iGeo in geo
        min_ = Vec3(min(iGeo.vertexPos0.x, min_.x), min(iGeo.vertexPos0.y, min_.y), min(iGeo.vertexPos0.z, min_.z))
        max_ = Vec3(max(iGeo.vertexPos0.x, max_.x), max(iGeo.vertexPos0.y, max_.y), max(iGeo.vertexPos0.z, max_.z))
        min_ = Vec3(min(iGeo.vertexPos1.x, min_.x), min(iGeo.vertexPos1.y, min_.y), min(iGeo.vertexPos1.z, min_.z))
        max_ = Vec3(max(iGeo.vertexPos1.x, max_.x), max(iGeo.vertexPos1.y, max_.y), max(iGeo.vertexPos1.z, max_.z))
        min_ = Vec3(min(iGeo.vertexPos2.x, min_.x), min(iGeo.vertexPos2.y, min_.y), min(iGeo.vertexPos2.z, min_.z))
        max_ = Vec3(max(iGeo.vertexPos2.x, max_.x), max(iGeo.vertexPos2.y, max_.y), max(iGeo.vertexPos2.z, max_.z))
    end

    return Aabb(min_, max_)
end


# merge two Aabb
function aabb__merge(a::Aabb, b::Aabb)::Aabb
    min_::Vec3 = Vec3(min(a.min.x, b.min.x), min(a.min.y, b.min.y), min(a.min.z, b.min.z))
    max_::Vec3 = Vec3(max(a.max.x, b.max.x), max(a.max.y, b.max.y), max(a.max.z, b.max.z))
    return Aabb(min_, max_)
end



# update AABB of BVH recursivly
function updateBvhAabbRecursive(bvhNode::BvhNodeA)
    if isLeaf(bvhNode)

        # * compute AABB of geometry inside leaf
        bvhNode.aabb = bvhLeafGeometry__calcAabb(bvhNode.leafGeometry)
    else
        # * compute AABB recursivly of children nodes, then compute AABB of THIS bvhNode by merging the AABB's of the childrens

        updateBvhAabbRecursive(bvhNode.children[0+1])
        updateBvhAabbRecursive(bvhNode.children[1+1])

        bvhNode.aabb = aabb__merge(bvhNode.children[0+1].aabb, bvhNode.children[1+1].aabb)
    end
end





function manualTestBuildBvh()

    toSegment::Array{BvhGeometryWithCenter} = []
    push!(toSegment, BvhGeometryWithCenter(BvhLeafGeometry(Vec3(0.0f0, 0.0f0, 0.0f0),Vec3(0.0f0, 0.0f0, 0.0f0),Vec3(0.0f0, 0.0f0, 0.0f0),ColorRgb_64bit(1.0,1.0,1.0)), Vec3(0.0f0, 0.0f0, 0.0f0)) )
    push!(toSegment, BvhGeometryWithCenter(BvhLeafGeometry(Vec3(1.0f0, 1.0f0, 1.0f0),Vec3(1.0f0, 1.0f0, 1.0f0),Vec3(1.0f0, 1.0f0, 1.0f0),ColorRgb_64bit(1.0,1.0,1.0)), Vec3(1.0f0, 1.0f0, 1.0f0)) )
    push!(toSegment, BvhGeometryWithCenter(BvhLeafGeometry(Vec3(2.0f0, 2.0f0, 2.0f0),Vec3(2.0f0, 2.0f0, 2.0f0),Vec3(2.0f0, 2.0f0, 2.0f0),ColorRgb_64bit(1.0,1.0,1.0)), Vec3(2.0f0, 2.0f0, 2.0f0)) )
    

    # * build actual BVH
    splitDimIdx::Int64 = 0
    rootBvhNode::BvhNodeA = buildBvhRaw(toSegment::Array{BvhGeometryWithCenter}, splitDimIdx)


    # * compute AABB's of all nodes
    updateBvhAabbRecursive(rootBvhNode)
end

if false
    manualTestBuildBvh()
end








function buildBvh(scene::Scene)::BvhNodeA
    toSegment::Array{BvhGeometryWithCenter} = []

    instanceIdx::Int64 = 0
    for iInstance in scene.instances

        rotationMatrix = convQuaternionTo44Matrix(scene.instances[ from0basedIdxToJulia(instanceIdx) ].rotation)

        for iMeshTriangle in iInstance.mesh.triangles

            vertexPos0::Vec3 = iInstance.mesh.vertices[ iMeshTriangle.vertexIdx0+1 ]
            vertexPos1::Vec3 = iInstance.mesh.vertices[ iMeshTriangle.vertexIdx1+1 ]
            vertexPos2::Vec3 = iInstance.mesh.vertices[ iMeshTriangle.vertexIdx2+1 ]

            vertexPos0 = matrixMul(rotationMatrix, vertexPos0)   + iInstance.worldPosition
            vertexPos1 = matrixMul(rotationMatrix, vertexPos1)   + iInstance.worldPosition
            vertexPos2 = matrixMul(rotationMatrix, vertexPos2)   + iInstance.worldPosition

            center::Vec3 = scale(vertexPos0 + vertexPos1 + vertexPos2, 1.0f0/3.0f0)

            geoWithCenter = BvhGeometryWithCenter(BvhLeafGeometry(vertexPos0, vertexPos1, vertexPos2,  iInstance.emission, iInstance.albedoColor, iInstance.materialIdx, instanceIdx), center)
            push!(toSegment, geoWithCenter)
        end
        
        instanceIdx = instanceIdx + 1
    end

    # * build actual BVH
    splitDimIdx::Int64 = 0
    rootBvhNode::BvhNodeA = buildBvhRaw(toSegment::Array{BvhGeometryWithCenter}, splitDimIdx)


    # * compute AABB's of all nodes
    updateBvhAabbRecursive(rootBvhNode)

    return rootBvhNode
end







# scene which is represented as a single BVH
mutable struct BvhScene
    camera::Camera
    bvh::BvhNodeA
    lights::Array{ExplicitLight}
end


# use BVH for ray traversal
function intersectRayVsScene(ray::Ray, rayInfo::RayInfo, scene::BvhScene)
    rayVsBvh(ray, rayInfo, scene.bvh)
end









###### LIGHTING HELPERS

###
#### compute PDF of the light source to be able to sample from the light source
#### /param areaOfLightsource area of the light
###function calcPdfOfAreaLight(distance::Float32, cosTheta::Float32, areaOfLightsource::Float32)::Float32
###    # see https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html#samplinglightsdirectly
###    return distance*distance / (cosTheta*areaOfLightsource)
###end









function vec3_rng():Vec3
    x::Float32 = rand(-1.0:eps():1.0)
    y::Float32 = rand(-1.0:eps():1.0)
    z::Float32 = rand(-1.0:eps():1.0)
    return Vec3(x,y,z)
end

function randomUnitVec3():Vec3
    while true
        p::Vec3 = vec3_rng()
        lenSquared::Float32 = calcMagSquare(p)
        if 1e-160 < lenSquared && lenSquared < 1.0f0
            return scale(p, 1.0f0 / sqrt(lenSquared))
        end
    end
end

function randomVecOnHemisphere(normal::Vec3):Vec3
    # see https://raytracing.github.io/books/RayTracingInOneWeekend.html#diffusematerials/asimplediffusematerial

    onUnitSphere::Vec3 = randomUnitVec3()
    if dot(onUnitSphere, normal) > 0.0f0
        return onUnitSphere
    else
        return -onUnitSphere
    end
end

# manual test
if false
    println(randomVecOnHemisphere(Vec3(1.0f0, 0.0f0, 0.0f0)))
end













function samplePointOnMesh(mesh::Mesh)
    # code to compute random point on mesh
    # HACKY OPTIMIZATION: (we dont take care of sampling by equal surface area here!)

    u::Float32 = rand(0.0:eps():1.0)
    v::Float32 = rand(0.0:eps():1.0)

    cntPolygons::Int64 = retLength(mesh.triangles)
    idxPolygon::Int64 = rand(0:cntPolygons-1) # rand() inclusive . 0 based index

    # UV interpolation over vertices
    # BUG LOW : this is probably buggy, but doesnt matter much

    w::Float32 = (2.0f0 - (u+v)) * 0.5f0

    uTick::Float32 = 1.0f0 - u
    vTick::Float32 = 1.0f0 - v
    wTick::Float32 = 1.0f0 - w

    p0::Vec3 = mesh.vertices[ mesh.triangles[1+idxPolygon].vertexIdx0+1 ]
    p1::Vec3 = mesh.vertices[ mesh.triangles[1+idxPolygon].vertexIdx1+1 ]
    p2::Vec3 = mesh.vertices[ mesh.triangles[1+idxPolygon].vertexIdx2+1 ]

    pointOnTriangle::Vec3 = scale(p0, uTick) + scale(p1, vTick) + scale(p2, wTick)

    return pointOnTriangle
end







# /param surfaceNormal normal (normalized)
function calcScatteringBrdf(dirToLight::Vec3, outgoingDir::Vec3, surfaceNormal::Vec3, material::MaterialA)
    if material.brdfType == 0 # lambertial BRDF model
        return 1.0
    elseif material.brdfType == 1 # ABg BRDF model
        a = 0.01f0
        b = 0.025f0
        g = 1.5f0

        # DEBUG sanity check
        #if dot(dirToLight, surfaceNormal) < 0.0
        #    println("err L")
        #end
        #if dot(outgoingDir, surfaceNormal) < 0.0
        #    println("err O")
        #end


        result = calcBrdf__ABg(dirToLight, surfaceNormal, outgoingDir, a, b, g)
        #print(result) # DBG
        return result
    end

end


# recursive function to compute incomming energy of ray computed via path-tracing
#
# returns color-energy transmitted by the ray
function recursiveRay(incommingRay::Ray, currentDepth::Int64, maxDepth::Int64,  scene::BvhScene, unacceleratedScene::Scene):ColorRgb_64bit
    if currentDepth == 1
        #println("ray with depth = 1 happened!") # DBG
    end
    
    if currentDepth >= maxDepth
        #println("exit")
        # max depth exceeded. we don't receive any energy from this ray
        return ColorRgb_64bit(0.0,0.0,0.0)
    end

    # intersect with scene
    incommingRayInfo::RayInfo = rayInfo__make()
    rayVsSceneIntersectionRes = intersectRayVsScene(incommingRay, incommingRayInfo, scene)

    if incommingRayInfo.isAnyHit

        materialOfSurface::MaterialA = unacceleratedScene.materials[ from0basedIdxToJulia(incommingRayInfo.surfaceMaterialIdx) ]


        #println("any hit") # DBG

        surfaceNormal::Vec3 = normalize(incommingRayInfo.surfaceNormal) # Ü  we need to normalize here because it can be unnormalized

        # HACK for debugging
        #return ColorRgb_64bit(1.0, 0.0, 0.0)
        # HACK for debugging
        #return ColorRgb_64bit(surfaceNormal.x, surfaceNormal.y, surfaceNormal.z)


        emissionSum::Float64 = incommingRayInfo.surfaceEmission.r + incommingRayInfo.surfaceEmission.g + incommingRayInfo.surfaceEmission.b
        if emissionSum > 1.0e-5
            #println("emission hack") # DBG
            
            # HACK : we terminate the ray if it has a emission. We do this to safe a bit of compute.
            return incommingRayInfo.surfaceEmission
        end

        # compute outgoing direction
        outDir::Vec3 = randomVecOnHemisphere(surfaceNormal)
        
        outgoingRayOrigin::Vec3 = ray__calcAbsolutePosition(incommingRay.origin, incommingRay.direction, Float32(incommingRayInfo.closestHitT)) + scale(surfaceNormal, 0.0001f0)
        

        # multiple importance sampling : weights
        multipleImpSamplingWeightStrategyA::Float32 = 1.0f0
        multipleImpSamplingWeightStrategyB::Float32 = 0.0f0
        if currentDepth == 0
            multipleImpSamplingWeightStrategyA = 0.5f0
            multipleImpSamplingWeightStrategyB = 0.5f0
        end

        # always do MIS
        if true
            multipleImpSamplingWeightStrategyA = 0.5f0
            multipleImpSamplingWeightStrategyB = 0.5f0
        end

        


        

        incomingEnergy::ColorRgb_64bit = ColorRgb_64bit(0.0,0.0,0.0) # declaration
        outgoingEnergy::ColorRgb_64bit = ColorRgb_64bit(0.0,0.0,0.0) # declaration
        cosTheta::Float32 = 0.0f0 # declaration

        outgoingEnergySum::ColorRgb_64bit = ColorRgb_64bit(0.0,0.0,0.0)


        # sampling strategy : sample material
        outgoingRay::Ray = ray__make(outgoingRayOrigin, outDir, Float32(1e20))
        incomingEnergy = recursiveRay(outgoingRay, currentDepth+1, maxDepth,  scene, unacceleratedScene)
        
        cosTheta = dot(surfaceNormal, outDir)
        if cosTheta < 0.0f0
            #println("info: cosTheta is smaller than  0.0")
        end
        
        lightAttenuation::Float32 = 1.0
        pdfValue::Float32 = 1.0 / cosTheta

        scatteringPdf::Float32 = 1.0 # lambertian
        scatteringPdf = calcScatteringBrdf(outDir, -incommingRay.direction, surfaceNormal, materialOfSurface)

        energyFromScatter::ColorRgb_64bit = scale(incomingEnergy, Float64(lightAttenuation * (scatteringPdf / pdfValue)))

        outgoingEnergy = energyFromScatter
        outgoingEnergy = colorRgb__mulComponents(outgoingEnergy, incommingRayInfo.surfaceAlbedoColor)
        outgoingEnergySum = outgoingEnergySum + scale(outgoingEnergy, Float64(multipleImpSamplingWeightStrategyA))


        ###outgoingEnergy = colorRgb__mulComponents(incomingEnergy, incommingRayInfo.surfaceAlbedoColor)
        ###outgoingEnergy = colorRgb__mulComponents(outgoingEnergy, ColorRgb_64bit(brdfResult, brdfResult, brdfResult))
        ###outgoingEnergySum = outgoingEnergySum + scale(outgoingEnergy, Float64(multipleImpSamplingWeightStrategyA))




        if multipleImpSamplingWeightStrategyB > 1.0e-6
            # sampling strategy : sample light

            pointOnLight::Vec3 = Vec3(0.0f0, 0.0f0, 0.0f0) # sampled point on light

            lightsourceInstanceIdx::Int64 = -1
            # TODO MID : randomly select the index from the array
            lightsourceInstanceIdx = unacceleratedScene.lightsourceInstanceIndices[ from0basedIdxToJulia(0) ]

            # sample point on light
            lightsourceMesh::Mesh = unacceleratedScene.instances[ from0basedIdxToJulia(lightsourceInstanceIdx) ].mesh # mesh of the lightsource
            pointOnLight = samplePointOnMesh(lightsourceMesh)

            rotationMatrix = convQuaternionTo44Matrix(unacceleratedScene.instances[ from0basedIdxToJulia(lightsourceInstanceIdx) ].rotation)
            pointOnLight = matrixMul(rotationMatrix, pointOnLight)   + unacceleratedScene.instances[ from0basedIdxToJulia(lightsourceInstanceIdx) ].worldPosition


            outDiff::Vec3 = pointOnLight-outgoingRayOrigin
            outDir = normalize(outDiff)

            #println(outDir)
            #println(surfaceNormal)

            cosTheta = dot(surfaceNormal, outDir)
            if cosTheta < 0.0f0
                #println("info: cosTheta B is < 0.0") # DBG
            end
            ###cosTheta = max(cosTheta, 0.000001f0)
            
            #println(cosTheta)
            #exit(0)

            
            # how to sample the light
            # see https://raytracing.github.io/books/RayTracingTheRestOfYourLife.html#samplinglightsdirectly

            distanceToLight::Float32 = calcMag(outDiff)


            
            # compute visibilty term by shooting ray toward ligthsource and checking if it did hit exactly the emitter            
            visibilityTerm::Float32 = 1.0

            lightsourceRay::Ray = ray__make(outgoingRayOrigin, outDir, Float32(1e20))
            lightsourceRayInfo::RayInfo = rayInfo__make()
            rayVsSceneIntersectionRes = intersectRayVsScene(lightsourceRay, lightsourceRayInfo, scene)

            # now we check if the ray did hit the specific light source which is the instance
            isNoIntersectionWithLightsource::Bool = false
            if lightsourceRayInfo.isAnyHit && lightsourceRayInfo.hitInstanceIdx != lightsourceInstanceIdx
                visibilityTerm = 0.0
                isNoIntersectionWithLightsource = true
            end



            lightSourceArea::Float32 = 1.0f0 # TODO TODO TODO : use actual area of light

            ### DEPRECATED
            ###pdfLightsampling::Float32 = calcPdfOfAreaLight(distanceToLight, cosTheta, lightSourceArea)




            ## BEGIN calc PDF of light
            absCosThetaOut::Float32 = abs(cosTheta)

            # we need to compute the PDF with respect to the solid angle
            pdfWithRespectToSolidAngle::Float32 = absCosThetaOut / (distanceToLight*distanceToLight)

            #println(absCosThetaOut)
            #println(distanceToLight)
            #println(pdfWithRespectToSolidAngle) # DBG
            #exit(0)

            # compute PDF of sampling lightsource
            pdfLightsampling = 0.0
            if isNoIntersectionWithLightsource 
                pdfLightsampling = (1.0f0 / pdfWithRespectToSolidAngle) * (1.0f0 / lightSourceArea)
            end

            ## END calc PDF of light

            #println(pdfLightsampling) # DBG
            #exit(0)




            ###lightAttenuation = 1.0 / (distanceToLight*distanceToLight) # actual attenuation to lightsource

            scatteringPdf = 1.0 # lambertian
            scatteringPdf = calcScatteringBrdf(outDir, -incommingRay.direction, surfaceNormal, materialOfSurface)

            # color and strength of the light
            incomingEnergy = ColorRgb_64bit(20.0,20.0,20.0)
            ###incomingEnergy = scale(incomingEnergy, Float64(visibilityTerm))

            #energyFromScatter = scale(incomingEnergy, Float64(lightAttenuation * (scatteringPdf / pdfLightsampling)))
            if pdfLightsampling < 1e-7
                energyFromScatter = ColorRgb_64bit(0.0, 0.0, 0.0)
            else
                temp0 = (1.0f0 / pdfLightsampling)
                #println(temp0)
                energyFromScatter = scale(incomingEnergy, Float64(temp0))

                #energyFromScatter = ColorRgb_64bit(1.0, 1.0, 1.0) # HACK HACK HACK
            end

            outgoingEnergy = energyFromScatter
            outgoingEnergy = colorRgb__mulComponents(outgoingEnergy, incommingRayInfo.surfaceAlbedoColor)
            outgoingEnergySum = outgoingEnergySum + scale(outgoingEnergy, Float64(multipleImpSamplingWeightStrategyB))
        end

        return outgoingEnergySum

    else
        # ray didn't hit anything, we don't receive anything from the ray
        return ColorRgb_64bit(0.0,0.0,0.0)
    end

end





function linearToGammaCorrection(x::AbstractFloat)::AbstractFloat
    # see https://raytracing.github.io/books/RayTracingInOneWeekend.html
    if x < 0.0
       return 0.0
    end
    return sqrt(x)
 end




function testEntry_A()


    sizeX::Int32 = 320 # 40
    sizeY::Int32 = 320 # 40


    

    cam = Camera(Vec3(1.911f0, 0.01222f0, 0.01232f0), Vec3(0.0, 0.000001, 0.99999), Vec3(0.0, 0.5, 0.0), Vec3(0.5, 0.0, 0.0))

    #sphereA = Sphere(Vec3(0.0, 0.0, 3.0), 1.0)

    instancesA::Array{Instance} = Instance[]

    boxB_path = "C:\\Users\\rober\\fsRoot\\art_blender\\geoA\\boxB.obj"
    meshBoxB::Mesh = readObj(boxB_path)


    testMeshC_path = "C:\\Users\\rober\\fsRoot\\art_blender\\geoA\\test_object_intersectedA.obj"
    testMeshC::Mesh = readObj(testMeshC_path)

    if true # codeblock
        #meshVertices::Array{Vec3} = Vec3[]
        meshVertices = Vec3[]
        push!(meshVertices, Vec3(0.04f0, 2.02f0, 0.01f0))
        push!(meshVertices, Vec3(2.05f0, 0.01f0, 0.02f0))
        push!(meshVertices, Vec3(-0.02f0, 0.01f0, 0.03f0))

        #meshTriangles::Array{MeshTriangle} = MeshTriangle[]
        meshTriangles = MeshTriangle[]
        push!(meshTriangles, MeshTriangle(0, 1, 2))

        mesh = Mesh(meshVertices, meshTriangles)

        instanceWorldPosition::Vec3 = Vec3(0.0f0, 0.0f0, -5.0f0)
        instanceRotation = makeQuaternionFromEuler(Vec3(1.0, 0.0, 0.0), 0.0000001)
        albedoColor::ColorRgb_64bit = ColorRgb_64bit(1.0, 1.0, 1.0)
        emission::ColorRgb_64bit = ColorRgb_64bit(150.0, 150.0, 150.0)
        materialIdx::Int64 = 0
        createdInstance = Instance(mesh, instanceWorldPosition, instanceRotation,  emission, albedoColor, materialIdx)


        push!(instancesA, createdInstance)
    end



    #mesh::Mesh
    if true # codeblock
        meshVertices::Array{Vec3} = Vec3[]
        #push!(meshVertices, Vec3(0.04f0, 1.02f0, 0.01f0))
        #push!(meshVertices, Vec3(1.05f0, 0.01f0, 0.02f0))
        #push!(meshVertices, Vec3(-0.02f0, 0.01f0, 0.03f0))

        # polgygon from mesh of box (I believe it is swizzled)

        push!(meshVertices, Vec3(-0.999127f0, 0.999825f0, 1.001047f0))
        push!(meshVertices, Vec3(0.999825f0, -1.001221f0, 0.998952f0))
        push!(meshVertices, Vec3(-1.000175f0, -1.000175f0, 0.999651f0))
        
        
        

        meshTriangles::Array{MeshTriangle} = MeshTriangle[]
        push!(meshTriangles, MeshTriangle(0, 1, 2))
        
        #mesh::Mesh = Mesh(meshVertices, meshTriangles)
        
        mesh::Mesh = meshBoxB

        instanceWorldPosition = Vec3(0.0f0, 0.0f0, 5.0f0)
        instanceRotation = makeQuaternionFromEuler(Vec3(1.0, 0.0, 0.0), 0.0000001)
        albedoColor = ColorRgb_64bit(1.0, 1.0, 1.0)
        emission = ColorRgb_64bit(0.0, 0.0, 0.0)
        materialIdx = 1
        createdInstance = Instance(mesh, instanceWorldPosition, instanceRotation,  emission, albedoColor, materialIdx)


        push!(instancesA, createdInstance)
    end

    if true # codeblock
        mesh = testMeshC

        instanceWorldPosition = Vec3(1.0f0, 0.0f0, 7.0f0)
        instanceRotation = makeQuaternionFromEuler(Vec3(1.0, 0.0, 0.0), 1e-4)
        albedoColor = ColorRgb_64bit(1.0, 0.0, 0.0)
        emission = ColorRgb_64bit(0.0, 0.0, 0.0)
        materialIdx = 2 # metallic ABg material
        createdInstance = Instance(mesh, instanceWorldPosition, instanceRotation,  emission, albedoColor, materialIdx)

        push!(instancesA, createdInstance)
    end

    lights::Array{ExplicitLight} = ExplicitLight[]
    push!(lights, ExplicitLight(Vec3(0.0f0, 0.0f0, 0.0f0), 2.0f0) )

    materials::Array{MaterialA} = MaterialA[]


    # material for light
    push!(materials, MaterialA(0))

    push!(materials, MaterialA(0))
    push!(materials, MaterialA(1)) # metallic ABg material

    lightsourceInstanceIndices::Array{Int64} = Int64[]

    # TODO MID: iterate over all instances of the scene and filter for instances which have emissive materials
    # HACKY
    push!(lightsourceInstanceIndices, 2)

    unacceleratedScene::Scene = Scene(cam, instancesA, lights, materials, lightsourceInstanceIndices)
    bvh::BvhNodeA = buildBvh(unacceleratedScene)
    scene::BvhScene = BvhScene(cam, bvh, lights)



    #arr::Array{ColorRgb_64bit} = ColorRgb_64bit[]
    outputImg::Map2d = map2d__make(sizeX, sizeY)

    
    for iy in 0:(sizeY-1)
        for ix in 0:(sizeX-1)


            rayInfo::RayInfo = rayInfo__make()

            #println(string(Float32(ix) / Float32(sizeX-1)) * " " * string(Float32(iy) / Float32(sizeY-1)) )
            #println()
            rayDirA::Vec3 = cam__calcRayDir(scene.camera, Float32(ix) / Float32(sizeX-1), Float32(iy) / Float32(sizeY-1))
            rayDirA = normalize(rayDirA) # Ü  we need to normalize, because the ray-triangle intersection in software expect this to be normalized
            #println(rayDirA)

            rayA = ray__make(scene.camera.pos, rayDirA, Float32(1e20))

            # orthogonal camera projection
            #rayA = ray__make(Vec3((Float32(ix) / 40.0f0) * 5.0f0, (Float32(iy) / 40.0f0) * 5.0f0, 0.0f0), Vec3(0.0f0, 0.0f0, 1.0f0), 1e20)

            maxRayDepth::Int64 = 3
            

            primaryIntegrationRayCount::Int64 = 750 # 50 is good but noisy   # 10

            # simple integration
            incommingEnergyAccumulator::ColorRgb_64bit = ColorRgb_64bit(0.0,0.0,0.0)
            for integrationIt in 1:primaryIntegrationRayCount
                incommingEnergy::ColorRgb_64bit = recursiveRay(rayA, 0, maxRayDepth,  scene, unacceleratedScene)
                incommingEnergy = scale(incommingEnergy, 1.0/primaryIntegrationRayCount)
                
                incommingEnergyAccumulator = incommingEnergyAccumulator + incommingEnergy
            end

            resColor::ColorRgb_64bit = incommingEnergyAccumulator


            #rayVsSceneIntersectionRes = intersectRayVsScene(rayA, rayInfo, scene)
            #
            #resColor::ColorRgb_64bit = ColorRgb_64bit(0.0, 0.0, 0.0)
            #
            #if rayInfo.isAnyHit
            #
            #    intersectionWorldPosition::Vec3 = cam__calcAbsolutePosition(scene.camera, rayDirA, convert(Float32, rayInfo.closestHitT))
            #
            #    # compute direction to light
            #    diffFromSurfaceToLight::Vec3 = scene.lights[0+1].pos - intersectionWorldPosition
            #    dirToLight::Vec3 = normalize(diffFromSurfaceToLight)
            #
            #    surfaceNormal::Vec3 = normalize(rayInfo.surfaceNormal) # Ü  we need to normalize here because it can be unnormalized
            #
            #    distanceToLight::Float32 = calcMag(diffFromSurfaceToLight)
            #
            #    # * shadow ray
            #    rayForLight = ray__make(intersectionWorldPosition + scale(surfaceNormal, Float32(1.0e-4)), dirToLight, distanceToLight)
            #    rayForLightRayInfo::RayInfo = rayInfo__make()
            #    rayVsSceneIntersectionForLightrayRes = intersectRayVsScene(rayForLight, rayForLightRayInfo, scene)
            #
            #    incommingEnergy::ColorRgb_64bit = ColorRgb_64bit(scene.lights[0+1].intensity, scene.lights[0+1].intensity, scene.lights[0+1].intensity)
            #    if false && rayForLightRayInfo.isAnyHit # HACK  :  disabled shadow
            #        incommingEnergy = ColorRgb_64bit(0.0, 0.0, 0.0)
            #    end
            #
            #    shade::Float32 = 0.0
            #
            #    # lambertian BRDF
            #    #shade = shading__shadeLambertian(dirToLight, surfaceNormal)
            #
            #    # ABg BRDF
            #    outgoingDir::Vec3 = -rayDirA
            #    a = 0.01f0
            #    b = 0.0025f0
            #    g = 1.5f0
            #    shade = calcBrdf__ABg(dirToLight, surfaceNormal, outgoingDir, a, b, g)
            #
            #
            #    helperColorA::ColorRgb_64bit = colorRgb__mulComponents(rayInfo.surfaceBaseColor, incommingEnergy)
            #
            #    resColor = colorRgb__mulComponents(ColorRgb_64bit(shade, shade, shade), helperColorA)
            #end


            #push!(arr, resColor)

            #if rayInfo.isAnyHit
            #    print("x")
            #else
            #    print(".")
            #end

            # gamma-correction
            # see
            if true
                resColor = ColorRgb_64bit(linearToGammaCorrection(resColor.r), linearToGammaCorrection(resColor.g), linearToGammaCorrection(resColor.b))
            end

            outputImg.arr[1 + (ix + iy*sizeX)] = resColor

        end

        #println("")
    end


    # DEBUG IMAGE
    #for iColor in arr
    #    println(iColor)
    #end

    storePpm(outputImg)

end

if true
    testEntry_A()
end










# TODO  :  implement tone-mapping








# DONE :  implement ray-traced shadows!
# DONE  :  implement simple path tracing algorithm using recursive function






# BASICALLY DONE  :  sampling light: shoot ray and check if it did hit the specific light source or not
#          DONE              : check if light got intersected by looking at emissive of hit surface
#          TODO low priority : compute area of emissive instance
#          TODO low priority : compute interpolation on triangle correctly!