# Monte Carlo ray tracing flow chart algorithm pseudocode: A @goto-implementation and a recursion based implementation, any ideas for a more modern implementation?
**URL:** https://discourse.julialang.org/t/monte-carlo-ray-tracing-flow-chart-algorithm-pseudocode-a-goto-implementation-and-a-recursion-based-implementation-any-ideas-for-a-more-modern-implementation/108089
**Category:** Modelling & Simulations
**Tags:** monte-carlo
**Created:** [December 27, 2023, 8:45am UTC](https://discourse.julialang.org/t/monte-carlo-ray-tracing-flow-chart-algorithm-pseudocode-a-goto-implementation-and-a-recursion-based-implementation-any-ideas-for-a-more-modern-implementation/108089 "2023-12-27T08:45:04Z")
**Posts on this page:** 4
**Page:** 1
### Author:  [@NikoBiele](https://discourse.julialang.org/u/NikoBiele)
#### Post date: [December 27, 2023, 8:45am UTC](https://discourse.julialang.org/t/monte-carlo-ray-tracing-flow-chart-algorithm-pseudocode-a-goto-implementation-and-a-recursion-based-implementation-any-ideas-for-a-more-modern-implementation/108089/1 "2023-12-27T08:45:04Z")
Hi Julia-people 🙂
I am working with heat transfer, more specifically simulating thermal radiation heat transfer in participating media through Monte Carlo ray tracing. I am currently working with a 2D algorithm. The reason I am looking for inspiration and inputs is that I would like to port my code to be GPU code, and before I do that I wanted to ensure that my approach is optimal. I designed the following flow chart algorithm which should be general enough to work for both 1D, 2D and 3D.
By now I have designed two code implementations of it. The code is quite involved, so I will post only pseudo code here.
My first and original take on the code was a massive use of the @goto macro. I got this implementation to be very fast, ray tracing about 1 billion rays in five minutes on 16 threads in parallel. The code has been validated against analytical results. The pseudo code is as follows:
```julia
@label emitRay
emit ray() # emit ray
enclosure = whichEnclosure() # find out which enclosure we're in
walls = localWalls() # find the local walls
dist = distToSurface() # find the distance to the wall along direction
@label doesRayHitWall # label for use with @goto
if S < dist # if the distance to the wall is larger than the sampled distance
@scatter_or_absorb # then the ray is scattered or absorbed
point = updatePoint() # update the current point
if R_omega > omega # if a random number is larger than the scattering albedo
# then the ray is absorbed
enclosure = whichEnclosure()
absorb(enclosure) += 1
@goto emitRay
else # ray is scattered
scattering = sampledirection() # sample scattering direction
end
enclosure = whichEnclosure() # find out which enclosure we're in
walls = localWalls() # find the local walls
dist, wall_index = distToSurface() # find the distance to the wall along direction
if dist > S
@goto scatter_or_absorb # ray is scattered or absorbed
else
@goto rayHitsWall # go to label
end
else
@label rayHitsWall # ray hits a wall
solidwalls = solidwalls()
if wall_index == solidwalls # ray hits a solid wall
if R_alpha > alphaWall # then the ray is reflected by the wall
newDirection = sampleDirection()
else # the ray is absorbed by the wall
absorb(wall) += 1
end
@goto emitRay
else # ray hits an imaginary boundary
point = updatePoint()
enclosure = whichEnclosure()
walls = localWalls()
dist, index = distToSurface()
@goto doesRayHitWall
end
end
```
Even though I got the above implementation to execute very efficiently I recently learned that using @goto statements can be confusing and is an ‘oldschool’ way of coding. Therefore I also implemented the same flow chart in a recursion type code as described below:
```julia
function rayHitWall()
enclosure = whichEnclosure()
walls = localWalls()
dist, index = distToSurface()
if dist < S # if the remaining distance is shorted than the distance to the wall
solidwalls = solidWalls()
while index != solidwall && dist < S
point = updatePoint()
enclosure = whichEnclosure()
walls = localWalls()
dist, index = distToSurface()
solidwalls = solidWalls()
end
if dist > S # ray is absorbed or scattered
if R_omega > omega # ray is absorbed
enclosure = enclosure()
absorb(enclosure) += 1
return
else # ray is scattered
newdirection = sampleDirection()
rayHitWall() # call the function itself recursively
end
elseif R_alpha < alpha # ray is absorbed by the wall
enclosure = enclosure()
walls = localWalls()
dist, index = distToSurface()
solidwalls = solidWalls()
absorb(enclosure) += 1
return
else # ray is reflected
newdirection = sampleDirection()
rayHitWall() # call the function itself recursively
end
else
if R_omega > omega # ray is absorbed
enclosure = enclosure()
absorb(enclosure) += 1
return
else # ray is scattered
newdirection = sampleDirection()
rayHitWall() # call the function itself recursively
end
end
end
for i = 1:N_rays
emission = sampleEmission()
rayHitWall()
end
```
The above recursive code also works, however, it is significantly slower than the implementation based on @goto statements. Furthermore, there is a risk of getting ‘stackoverflow’ with this code.
But I was curious if someone else could come up with a more modern approach to implement the flowchart in (pseudo)code which would also be efficient, and possibly also applicable to GPU programming?
Any help would be appreciated. 🙂
---
### Author:  [@mihalybaci](https://discourse.julialang.org/u/mihalybaci)
#### Post date: [December 27, 2023, 2:31pm UTC](https://discourse.julialang.org/t/monte-carlo-ray-tracing-flow-chart-algorithm-pseudocode-a-goto-implementation-and-a-recursion-based-implementation-any-ideas-for-a-more-modern-implementation/108089/2 "2023-12-27T14:31:08Z")
I can’t help with code, but there have been a few other ray tracing codes presented here. [This thread](https://discourse.julialang.org/t/ray-tracing-in-a-week-end-julia-vs-simd-optimized-c/72958) has a long discussion, and if you check out the links right after the original post, there are other threads to follow. Maybe one of those will give you some new ideas.
---
### Author:  [@Zentrik](https://discourse.julialang.org/u/Zentrik)
#### Post date: [December 27, 2023, 10:56pm UTC](https://discourse.julialang.org/t/monte-carlo-ray-tracing-flow-chart-algorithm-pseudocode-a-goto-implementation-and-a-recursion-based-implementation-any-ideas-for-a-more-modern-implementation/108089/3 "2023-12-27T22:56:08Z")
This looks fairly similar to graphics ray tracing, my code in Julia and c++ for GPU implementation is here [GitHub - Zentrik/Renderer at GPU-WaveFront](https://github.com/Zentrik/Renderer/tree/GPU-WaveFront). It uses a technique called wavefront rendering which is more suitable for the GPU, [https://www.google.com/url?q=https://research.nvidia.com/sites/default/files/pubs/2013-07\_Megakernels-Considered-Harmful/laine2013hpg\_paper.pdf&sa=U&ved=2ahUKEwir-Ov\_rK-DAxWhYEEAHTEiDVAQFnoECAEQAg&usg=AOvVaw0grbKUz0v\_4BivlR7Y97po](https://www.google.com/url?q=https://research.nvidia.com/sites/default/files/pubs/2013-07_Megakernels-Considered-Harmful/laine2013hpg_paper.pdf&sa=U&ved=2ahUKEwir-Ov_rK-DAxWhYEEAHTEiDVAQFnoECAEQAg&usg=AOvVaw0grbKUz0v_4BivlR7Y97po).
Though I would suggest finding a existing approach even if it’s in another language and using that.
EDIT: My Julia GPU code has a simpler version here, [Renderer/src/simple.jl at GPU-WaveFront · Zentrik/Renderer · GitHub](https://github.com/Zentrik/Renderer/blob/GPU-WaveFront/src/simple.jl). This should have the same performance as the other file but is simpler to understand.
---
### Author:  [@NikoBiele](https://discourse.julialang.org/u/NikoBiele)
#### Post date: [December 29, 2023, 11:59pm UTC](https://discourse.julialang.org/t/monte-carlo-ray-tracing-flow-chart-algorithm-pseudocode-a-goto-implementation-and-a-recursion-based-implementation-any-ideas-for-a-more-modern-implementation/108089/4 "2023-12-29T23:59:58Z")
Thanks a lot for both your comments! 🙂