# Ray Tracing in a week-end - Julia vs SIMD-optimized C++

**URL:** <https://discourse.julialang.org/t/ray-tracing-in-a-week-end-julia-vs-simd-optimized-c/72958>\
**Category:** Performance\
**Tags:** gpu, graphics, benchmark, pluto, staticarrays\
**Created:** [December 11, 2021, 11:58pm UTC](https://discourse.julialang.org/t/ray-tracing-in-a-week-end-julia-vs-simd-optimized-c/72958 "2021-12-11T23:58:13Z")\
**Posts on this page:** 1\
**Showing post:** 42

<div class="post-metadata">

**Author:** ![c42f](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/c42f/32/52842_2.png) [@c42f](https://discourse.julialang.org/u/c42f)\
**Post date:** [December 16, 2021, 3:19am UTC](https://discourse.julialang.org/t/ray-tracing-in-a-week-end-julia-vs-simd-optimized-c/72958/42 "2021-12-16T03:19:07Z")

</div>

Great thread! Fun fact: The design of `StaticArrays` was indirectly informed by the arrays-for-geometry-types I wish I’d had back when I used to write a lot of C++ rendering code (in the [aqsis](https://github.com/aqsis/aqsis) project and several toy renderers).

Some thoughts after looking over the project:

### StaticArrays usage

In the README that you’ve got a section “using Base.getproperty() for `vec.x` instead of `vec[1]`”. You could do this very easily using [StaticArrays.FieldVector](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#FieldVector-1):

```julia-auto
struct Vec3{T} <: FieldVector{3, T}
    x::T
    y::T
    z::T
end

```

On the other hand, there’s something quite nice about using plain `SVector`. For instance, you can use the short syntax:

- `SA[1,2,3]` is a shorter version of `@SVector[1,2,3]`
- `SA{T}[1,2,3]` is shorthand for `@SVector T[1,2,3]`
- There’s also aliases `SA_F64` and `SA_F32` for `SA{Float64}` and `SA{Float32}` if you’re writing concrete code.

Actually… here’s a nice idea — what if we added a `getproperty` feature to `StaticArrays` so that `v.x`, `v.y` and `v.z` “just worked” by default like everyone wants them to? Here, I grant you a [letter of marque 📜](https://en.wikipedia.org/wiki/Letter_of_marque) to use the following `getproperty` implementation:

```julia
function Base.getproperty(v::SVector{3}, name::Symbol)
    name === :x ? getfield(v, :data)[1] :
    name === :y ? getfield(v, :data)[2] :
    name === :z ? getfield(v, :data)[3] : getfield(v, name)
end

```

### Notes on the code

It’s been mentioned above, but to mention it again - `const _rng = Xoroshiro128Plus()` is no good with threads. (There’s a lot of globals though! I’m not sure it’s used from your threads section?)

A few thoughts about the loop here: [RayTracingWeekend.jl/src/proto.jl at 9d88dfa7e48e82101e87988ee73305f130121f66 · claforte/RayTracingWeekend.jl · GitHub](https://github.com/claforte/RayTracingWeekend.jl/blob/9d88dfa7e48e82101e87988ee73305f130121f66/src/proto.jl#L242)

- If column major vs row major in the output image matters for perf… you know it’s a toy renderer or toy scene as this is the outer loop 🙂
- Your u,v coordinates for your samples are “off by half a pixel” - you may want `u = (j-T(0.5))/nx` instead. If you’re intent on microoptimizing this outer loop, you could consider computing `dx = T(1)/nx` outside the loop doing `u = (j-T(0.5))*dx` to avoid the division.
- The dot product at [RayTracingWeekend.jl/src/proto.jl at 9d88dfa7e48e82101e87988ee73305f130121f66 · claforte/RayTracingWeekend.jl · GitHub](https://github.com/claforte/RayTracingWeekend.jl/blob/9d88dfa7e48e82101e87988ee73305f130121f66/src/proto.jl#L277) seems unnecessary if your ray directions are normalized?
- Is there a reason to have `T` in the abstract `Material` type? You can probably remove & simplify [RayTracingWeekend.jl/src/proto.jl at 9d88dfa7e48e82101e87988ee73305f130121f66 · claforte/RayTracingWeekend.jl · GitHub](https://github.com/claforte/RayTracingWeekend.jl/blob/9d88dfa7e48e82101e87988ee73305f130121f66/src/proto.jl#L338)
- The comment at [RayTracingWeekend.jl/src/proto.jl at 9d88dfa7e48e82101e87988ee73305f130121f66 · claforte/RayTracingWeekend.jl · GitHub](https://github.com/claforte/RayTracingWeekend.jl/blob/9d88dfa7e48e82101e87988ee73305f130121f66/src/proto.jl#L342) seems a bit misleading. It’s the distance, not the vector?
- There may be a performance benefit from having the abstract `Material{T}` at [RayTracingWeekend.jl/src/proto.jl at 9d88dfa7e48e82101e87988ee73305f130121f66 · claforte/RayTracingWeekend.jl · GitHub](https://github.com/claforte/RayTracingWeekend.jl/blob/9d88dfa7e48e82101e87988ee73305f130121f66/src/proto.jl#L352) due to the small `Union` optimization. But this optimization may go away and performance fall off a cliff if you increase the number of `Material` subtypes beyond a few.
- `HittableList = Vector{Hittable}` has abstractly typed elements, even though there’s only one concretely typed `Sphere <: Hittable`. If you add other hittables you’ll probably find a big performance reduction here. What happens if you replace `Vector{Hittable}` with `Vector{Sphere}`? (Yes this makes things less flexible but if there’s a huge performance boost we can recover the abstraction with an alternative method.)
- Why are all your `Material` subtypes mutable? This seems unnecessary and might have a performance impact.
- Does this `@simd` actually do anything [RayTracingWeekend.jl/src/proto.jl at 9d88dfa7e48e82101e87988ee73305f130121f66 · claforte/RayTracingWeekend.jl · GitHub](https://github.com/claforte/RayTracingWeekend.jl/blob/9d88dfa7e48e82101e87988ee73305f130121f66/src/proto.jl#L447) ? I’d guess not because `Vector{Hittable}` has an abstract element type?

A general thought about SIMD (on which I’m not an expert!) — data layout can be important and 3-vectors and arrays of structs aren’t a great data layout. The problem is that such structs are a great programming abstraction! Packages like [`StructOfArrays`](https://github.com/JuliaArrays/StructArrays.jl) can help to change the data layout into a SIMD-friendly struct-of-arrays format while preserving the array-of-structs abstraction you want to use in the code.

### Making the materials and hit system fully concrete with tight inner loops

If you’re willing to make larger changes to your `ray_color` function and replace the (intuitive!) recursive formulation with a ray batch system, you could have fully concrete inner loops and you’d have a much higher chance of getting really tight loops which SIMD properly.

I think there’s a few steps to this:

1. Replace the recursion in `ray_color` with a loop and a stack of hit events (make sure that the stack gets reused to avoid allocations)!

2. Walk batches of rays together. Each time there’s a hit, instead of recursing into it immediately, store the HitRecord in the batch.

3. Remove the abstract `Material` from `HitRecord`, and replace it with something like the following **concrete** `MaterialRef` type. This allows hits to be grouped into like materials, with all the hits for the same type of material processed together in a batch.

```julia-auto
struct MaterialRef
    UInt16 material_type_index
    UInt16 index
end

```

Now you have the following (untested!) code sketch to store (“intern”) a new material into some existing arrays:

```julia-auto
function store_material(material_types, all_materials, new_material)
    material_type_index = findfirst(material_types, typeof(new_material))
    storage = all_materials[]
    push!(storage, new_material)
    return MaterialRef(material_type_index, lastindex(storage))
end

# Storage arrays for materials can be created like such
material_types = subtypes(Material)
all_materials = collect(Vector{T}() for T in material_types)

```

The point of this is that you can act on all materials **of the same type** in batch mode in a tight loop **without** dynamic dispatch and **with** inlining:

```julia-auto
function process_hits(hits, material_storage, material_index)
    for hit in hits
        if hit.material_ref.material_type_index == material_index
            material = mat_storage[hit.material_ref.index]
            # Do whatever you want with the material and the hit
            # The compiler knows the concrete type of `material_storage`
            # so use of `material` will be fast and inlined.
        end
    end
end

# Very rough usage
mat_idx = 1
mat_storage = all_materials[mat_idx]
process_hits(hits, all_materials[mat_idx], mat_idx)

```

This kind of batch-based pattern is likely to be helpful (perhaps necessary?) on the GPU.

#### Sampling patterns

You’re using `rand()` without stratification inside `render()` (and elsewhere like in `random_vec2_in_disk`) but this is very noisy for a given number of samples per pixel and you can consider stratified sampling. Simple `NxN` jittered sampling would be a lot better and is easy to implement, at least within `render()`, though it restricts `n_samples = N*N` for some `N`.

Alternatively it would be interesting to try a low discrepancy sequence from a package like [GitHub - JuliaMath/Sobol.jl: generation of Sobol low-discrepancy sequence (LDS) for the Julia language](https://github.com/stevengj/Sobol.jl).

---

_[View the full topic](https://discourse.julialang.org/t/ray-tracing-in-a-week-end-julia-vs-simd-optimized-c/72958)._