Doing Raytracer in Crystal

I came across this project today that compares different implementations of raytracer code.

I copied and ran the Ruby implementation in three different versions.

They need two rubygems to do the image creation. Here are times below.

# ruby-2.6.6
ruby raytracer.rb
Completed in 14373.036497 ms

# ruby 2.7.2
ruby raytracer.rb
Completed in 16478.567917 ms

# jruby-9.2.16.0
ruby raytracer.rb
raytracer.rb:348: warning: `-' after local variable or literal is interpreted as binary operator
raytracer.rb:348: warning: even though it seems like unary operator
Completed in 13724.8 ms

The Ruby code is only 351 loc and easy to read.
This looks like something fun, and straightforward, to do in Crystal, but needs the equivalent of the gems imageruby and imageruby-bmp for Crystal.

It creates a ruby-tracer.bmp file, which I converted to a *.jpeg to upload.

ruby-raytracer-266

Do such shards, or other projects, exist for Crystal?

6 Likes

There’s one in the samples directory of Crystal: crystal/raytracer.cr at master · crystal-lang/crystal · GitHub

It uses SDL, though. But combining it with stumpy_png should be straight-forward.

There’s also a text ray tracer there: crystal/text_raytracer.cr at master · crystal-lang/crystal · GitHub

5 Likes

Here’s that raytracer code combined with stumpy_png. It runs in 0.19s on my machine.

require "stumpy_png"

WIDTH     = 1280
HEIGHT    =  720
FOV       = 45.0
MAX_DEPTH =    6

struct Vec3
  getter :x
  getter :y
  getter :z

  def initialize
    @x, @y, @z = 0.0, 0.0, 0.0
  end

  def initialize(value)
    @x, @y, @z = value, value, value
  end

  def initialize(@x, @y, @z)
  end

  {% for op in %w(+ - * /) %}
    def {{op.id}}(other : Vec3)
      Vec3.new(@x {{op.id}} other.x, @y {{op.id}} other.y, @z {{op.id}} other.z)
    end

    def {{op.id}}(other : Float)
      Vec3.new(@x {{op.id}} other, @y {{op.id}} other, @z {{op.id}} other)
    end
  {% end %}

  def -
    Vec3.new(-@x, -@y, -@z)
  end

  def dot(other)
    @x * other.x + @y * other.y + @z * other.z
  end

  def magnitude
    Math.sqrt(dot(self))
  end

  def normalize
    m = magnitude
    Vec3.new(@x / m, @y / m, @z / m)
  end
end

record Ray, start : Vec3, dir : Vec3

class Sphere
  getter :color
  getter :reflection
  getter :transparency

  def initialize(@center : Vec3, @radius : Float64, @color : Vec3, @reflection = 0.0, @transparency = 0.0)
  end

  def intersects?(ray)
    vl = @center - ray.start
    a = vl.dot(ray.dir)
    return false if a < 0

    b2 = vl.dot(vl) - a * a
    r2 = @radius * @radius
    return false if b2 > r2

    true
  end

  def intersect(ray, distance)
    vl = @center - ray.start
    a = vl.dot(ray.dir)
    return nil if a < 0

    b2 = vl.dot(vl) - a * a
    r2 = @radius * @radius
    return nil if b2 > r2

    c = Math.sqrt(r2 - b2)
    near = a - c
    far = a + c
    near < 0 ? far : near
  end

  def normalize(v)
    (v - @center).normalize
  end
end

record Light, position : Vec3, color : Vec3
record Scene, objects : Array(Sphere), lights : Array(Light)

def trace(ray, scene, depth)
  nearest = 1e9
  obj = nil
  result = Vec3.new

  scene.objects.each do |o|
    distance = 1e9
    if (distance = o.intersect(ray, distance)) && distance < nearest
      nearest = distance
      obj = o
    end
  end

  if obj
    point_of_hit = ray.dir * nearest
    point_of_hit += ray.start
    normal = obj.normalize(point_of_hit)
    inside = false
    dot_normal_ray = normal.dot(ray.dir)
    if dot_normal_ray > 0
      inside = true
      normal = -normal
      dot_normal_ray = -dot_normal_ray
    end

    reflection_ratio = obj.reflection
    normE5 = normal * 1.0e-5

    scene.lights.each do |lgt|
      light_direction = (lgt.position - point_of_hit).normalize
      r = Ray.new(point_of_hit + normE5, light_direction)

      # go through the scene check whether we're blocked from the lights
      blocked = scene.objects.any? &.intersects? r

      unless blocked
        temp = lgt.color
        temp *= Math.max(0.0, normal.dot(light_direction))
        temp *= obj.color
        temp *= (1.0 - reflection_ratio)
        result += temp
      end
    end

    facing = Math.max(0.0, -dot_normal_ray)
    fresneleffect = reflection_ratio + (1.0 - reflection_ratio) * ((1.0 - facing) ** 5.0)

    # compute reflection
    if depth < MAX_DEPTH && reflection_ratio > 0
      reflection_direction = ray.dir - normal * 2.0 * dot_normal_ray
      reflection = trace(Ray.new(point_of_hit + normE5, reflection_direction), scene, depth + 1)
      result += reflection * fresneleffect
    end

    # compute refraction
    if depth < MAX_DEPTH && (obj.transparency > 0.0)
      ior = 1.5
      ce = ray.dir.dot(normal) * -1.0
      ior = inside ? 1.0 / ior : ior
      eta = 1.0 / ior
      gf = (ray.dir + normal * ce) * eta
      sin_t1_2 = 1.0 - ce * ce
      sin_t2_2 = sin_t1_2 * (eta * eta)
      if sin_t2_2 < 1.0
        gc = normal * Math.sqrt(1 - sin_t2_2)
        refraction_direction = gf - gc
        refraction = trace(Ray.new(point_of_hit - normal * 1.0e-4, refraction_direction),
          scene, depth + 1)
        result += refraction * (1.0 - fresneleffect) * obj.transparency
      end
    end
  end

  result
end

def render(scene, canvas)
  eye = Vec3.new
  h = Math.tan(FOV / 360.0 * 2.0 * Math::PI / 2.0) * 2.0
  ww = canvas.width.to_f64
  hh = canvas.height.to_f64
  w = h * ww / hh

  i = 0
  HEIGHT.times do |y|
    yy = y.to_f64
    WIDTH.times do |x|
      xx = x.to_f64
      dir = Vec3.new((xx - ww / 2.0) / ww * w,
        (hh / 2.0 - yy) / hh * h,
        -1.0).normalize
      pixel = trace(Ray.new(eye, dir), scene, 0.0)
      r = Math.min(255, (pixel.x * 255.0).round.to_i)
      g = Math.min(255, (pixel.y * 255.0).round.to_i)
      b = Math.min(255, (pixel.z * 255.0).round.to_i)
      canvas.set(x, y, StumpyCore::RGBA.from_rgb(r, g, b))
    end
  end
end

scene = Scene.new(
  [
    Sphere.new(Vec3.new(0.0, -10002.0, -20.0), 10000.0, Vec3.new(0.8, 0.8, 0.8)),
    Sphere.new(Vec3.new(0.0, 2.0, -20.0), 4.0, Vec3.new(0.8, 0.5, 0.5), 0.5),
    Sphere.new(Vec3.new(5.0, 0.0, -15.0), 2.0, Vec3.new(0.3, 0.8, 0.8), 0.2),
    Sphere.new(Vec3.new(-5.0, 0.0, -15.0), 2.0, Vec3.new(0.3, 0.5, 0.8), 0.2),
    Sphere.new(Vec3.new(-2.0, -1.0, -10.0), 1.0, Vec3.new(0.1, 0.1, 0.1), 0.1, 0.8),
  ],
  [
    Light.new(Vec3.new(-10.0, 20.0, 30.0), Vec3.new(2.0, 2.0, 2.0)),
  ]
)

canvas = StumpyCore::Canvas.new(WIDTH, HEIGHT)

render scene, canvas

StumpyPNG.write(canvas, "scene.png")

The output looks like this:

9 Likes

Now that I have some free time :smile:, I can try to do the raytracer for the comparison page.

1 Like

Hey @asterite I ran your raytracer on my system and it consistently does between 151ms ~ 155ms. Of course same picture.

EDIT:
I added this code at the end to time it.

t1 = Time.monotonic
canvas = StumpyCore::Canvas.new(WIDTH, HEIGHT)

render scene, canvas
t2 = (Time.monotonic - t1).total_seconds.round(6)
  
puts "Completed in #{t2*1000} ms"

StumpyPNG.write(canvas, "scene.png")
4 Likes

Cool! I think porting the original Ruby code with the cooler picture (in my opinion) might give a better idea of how it compares against Ruby. Maybe I’ll do it if I have some time…

2 Likes

A few years ago when I first discovered Crystal, I happened to be reading this book at the time:

The Ray Tracer Challenge

And I went through the main chapters and did the whole thing in Crystal - it was a BLAST and in the back of my mind I was eventually (and still might) write up a blog post / review of the book.

Anyway, this thread reminded me of it and I just dug up the code and posted it to GitHub. Take a peek!

Disclaimer: I just copied from my working folder and threw it up on a repository a few minutes ago. It might not work! But if it’s useful / interesting in any way, enjoy :slight_smile:

2 Likes

Here’s the original Ruby code ported to Crystal with StumpyPNG:

require "stumpy_png"

struct Vector
  getter x, y, z

  def initialize(@x : Float64, @y : Float64, @z : Float64)
  end

  def self.times(k, v)
    Vector.new(k * v.x, k * v.y, k * v.z)
  end

  def self.times(k, v)
    Vector.new(k * v.x, k * v.y, k * v.z)
  end

  def self.minus(v1, v2)
    Vector.new(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z)
  end

  def self.plus(v1, v2)
    Vector.new(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z)
  end

  def self.dot(v1, v2)
    v1.x * v2.x + v1.y * v2.y + v1.z * v2.z
  end

  def self.mag(v)
    Math.sqrt(v.x * v.x + v.y * v.y + v.z * v.z)
  end

  def self.norm(v)
    mag = Vector.mag(v)
    if (mag == 0)
      div = Float64::INFINITY
    else
      div = 1.0 / mag
    end
    Vector.times(div, v)
  end

  def self.cross(v1, v2)
    Vector.new(v1.y * v2.z - v1.z * v2.y, v1.z * v2.x - v1.x * v2.z, v1.x * v2.y - v1.y * v2.x)
  end
end

struct Color
  getter r, g, b

  def initialize(@r : Float64, @g : Float64, @b : Float64)
  end

  def self.scale(k, v)
    Color.new(k * v.r, k * v.g, k * v.b)
  end

  def self.plus(v1, v2)
    Color.new(v1.r + v2.r, v1.g + v2.g, v1.b + v2.b)
  end

  def self.times(v1, v2)
    Color.new(v1.r * v2.r, v1.g * v2.g, v1.b * v2.b)
  end

  def self.toDrawingColor(c)
    clamp = ->(d : Float64) do
      return 1 if d > 1
      return 0 if d < 0
      return d
    end
    r = (clamp.call(c.r)*255).floor
    g = (clamp.call(c.g)*255).floor
    b = (clamp.call(c.b)*255).floor
    {r, g, b}
  end
end

Color_white        = Color.new(1.0, 1.0, 1.0)
Color_grey         = Color.new(0.5, 0.5, 0.5)
Color_black        = Color.new(0.0, 0.0, 0.0)
Color_background   = Color_black
Color_defaultColor = Color_black

class Camera
  getter pos, forward, right, up

  def initialize(pos : Vector, lookAt)
    down = Vector.new(0.0, -1.0, 0.0)
    @pos = pos
    @forward = Vector.norm(Vector.minus(lookAt, @pos))
    @right = Vector.times(1.5, Vector.norm(Vector.cross(@forward, down)))
    @up = Vector.times(1.5, Vector.norm(Vector.cross(@forward, @right)))
  end
end

record Ray, start : Vector, dir : Vector

record Intersection, thing : Thing, ray : Ray, dist : Float64

record Light, pos : Vector, color : Color

abstract class Thing
end

class Sphere < Thing
  @radius2 : Float64

  def initialize(@center : Vector, radius : Float64, @_surface : Surface)
    @radius2 = radius*radius
  end

  def normal(pos)
    Vector.norm(Vector.minus(pos, @center))
  end

  def surface
    @_surface
  end

  def intersect(ray)
    eo = Vector.minus(@center, ray.start)
    v = Vector.dot(eo, ray.dir)
    dist = 0.0
    if (v >= 0)
      disc = @radius2 - (Vector.dot(eo, eo) - v * v)
      if (disc >= 0)
        dist = v - Math.sqrt(disc)
      end
    end
    if (dist == 0)
      return nil
    end
    Intersection.new(self, ray, dist)
  end
end

class Plane < Thing
  def initialize(@_norm : Vector, @offset : Float64, @_surface : Surface)
  end

  def normal(pos)
    return @_norm
  end

  def intersect(ray)
    denom = Vector.dot(@_norm, ray.dir)
    return nil if (denom > 0)
    dist = (Vector.dot(@_norm, ray.start) + @offset) / (-denom)
    return Intersection.new(self, ray, dist)
  end

  def surface
    @_surface
  end
end

abstract class Surface
end

class ShinySurface < Surface
  def diffuse(pos)
    return Color_white
  end

  def specular(pos)
    return Color_grey
  end

  def reflect(pos)
    return 0.7
  end

  def roughness
    return 250
  end
end

class CheckerboardSurface < Surface
  def diffuse(pos)
    return Color_white if ((pos.z).floor + (pos.x).floor) % 2 != 0
    return Color_black
  end

  def specular(pos)
    return Color_white
  end

  def reflect(pos)
    return 0.1 if ((pos.z).floor + (pos.x).floor) % 2 != 0
    return 0.7
  end

  def roughness
    return 250
  end
end

Surface_shiny        = ShinySurface.new
Surface_checkerboard = CheckerboardSurface.new

class RayTracer
  MaxDepth = 5

  def intersections(ray, scene)
    closest = Float64::INFINITY
    closestInter = nil
    scene.things.each do |item|
      inter = item.intersect(ray)
      if inter && inter.dist < closest
        closestInter = inter
        closest = inter.dist
      end
    end
    return closestInter
  end

  def testRay(ray, scene)
    isect = self.intersections(ray, scene)
    return isect.dist if isect
    return nil
  end

  def traceRay(ray, scene, depth)
    isect = self.intersections(ray, scene)
    return Color_background if isect.nil?
    return self.shade(isect, scene, depth)
  end

  def shade(isect : Intersection, scene, depth)
    d = isect.ray.dir
    pos = Vector.plus(Vector.times(isect.dist, d), isect.ray.start)
    normal = isect.thing.normal(pos)
    reflectDir = Vector.minus(d, Vector.times(2, Vector.times(Vector.dot(normal, d), normal)))
    naturalColor = Color.plus(Color_background, self.getNaturalColor(isect.thing, pos, normal, reflectDir, scene))

    if (depth >= MaxDepth)
      reflectedColor = Color_grey
    else
      reflectedColor = self.getReflectionColor(isect.thing, pos, normal, reflectDir, scene, depth)
    end
    return Color.plus(naturalColor, reflectedColor)
  end

  def getReflectionColor(thing, pos, normal, rd, scene, depth)
    return Color.scale(thing.surface.reflect(pos), self.traceRay(Ray.new(pos, rd), scene, depth + 1))
  end

  def getNaturalColor(thing, pos, norm, rd, scene)
    color = Color_defaultColor
    scene.lights.each do |light|
      color = self.addLight(color, light, pos, norm, scene, thing, rd)
    end
    return color
  end

  def addLight(col, light, pos, norm, scene, thing, rd)
    ldis = Vector.minus(light.pos, pos)
    livec = Vector.norm(ldis)
    neatIsect = self.testRay(Ray.new(pos, livec), scene)

    isInShadow = false
    isInShadow = neatIsect <= Vector.mag(ldis) if neatIsect

    return col if isInShadow

    illum = Vector.dot(livec, norm)

    lcolor = Color_defaultColor
    lcolor = Color.scale(illum, light.color) if illum > 0

    specular = Vector.dot(livec, Vector.norm(rd))
    scolor = Color_defaultColor
    scolor = Color.scale(specular ** thing.surface.roughness, light.color) if (specular > 0)

    return Color.plus(col, Color.plus(Color.times(thing.surface.diffuse(pos), lcolor), Color.times(thing.surface.specular(pos), scolor)))
  end

  def getPoint(x : Int32, y : Int32, camera, screenWidth, screenHeight)
    recenterX = ->(x : Int32) do
      (x - (screenWidth / 2.0)) / 2.0 / screenWidth
    end
    recenterY = ->(y : Int32) do
      -(y - (screenHeight / 2.0)) / 2.0 / screenHeight
    end
    return Vector.norm(Vector.plus(camera.forward, Vector.plus(Vector.times(recenterX.call(x), camera.right), Vector.times(recenterY.call(y), camera.up))))
  end

  def render(scene, image, screenWidth, screenHeight)
    (0..screenHeight - 1).each do |y|
      (0..screenWidth - 1).each do |x|
        color = self.traceRay(Ray.new(scene.camera.pos, self.getPoint(x, y, scene.camera, screenWidth, screenHeight)), scene, 0)
        r, g, b = Color.toDrawingColor(color)
        image.set(x, y, StumpyCore::RGBA.from_rgb(r, g, b))
      end
    end
  end
end

class DefaultScene
  getter :things, :lights, :camera

  def initialize
    @things = [
      Plane.new(Vector.new(0.0, 1.0, 0.0), 0.0, Surface_checkerboard),
      Sphere.new(Vector.new(0.0, 1.0, -0.25), 1.0, Surface_shiny),
      Sphere.new(Vector.new(-1.0, 0.5, 1.5), 0.5, Surface_shiny),
    ]
    @lights = [
      Light.new(Vector.new(-2.0, 2.5, 0.0), Color.new(0.49, 0.07, 0.07)),
      Light.new(Vector.new(1.5, 2.5, 1.5), Color.new(0.07, 0.07, 0.49)),
      Light.new(Vector.new(1.5, 2.5, -1.5), Color.new(0.07, 0.49, 0.071)),
      Light.new(Vector.new(0.0, 3.5, 0.0), Color.new(0.21, 0.21, 0.35)),
    ]
    @camera = Camera.new(Vector.new(3.0, 2.0, 4.0), Vector.new(-1.0, 0.5, 0.0))
  end
end

width = 500
height = 500
image = StumpyCore::Canvas.new(width, height)

t1 = Time.monotonic
rayTracer = RayTracer.new
scene = DefaultScene.new
rayTracer.render(scene, image, width, height)
t2 = Time.monotonic

puts "Completed in #{(t2 - t1).total_milliseconds}ms"

StumpyPNG.write(image, "ruby.png")

Compiled in release mode it runs in about 100ms, so more than 100x faster than the Ruby version. Here’s the output I get:

ruby

6 Likes

Cool! Are you planning on adding a PR to GitHub - edin/raytracer: Performance comparison of various compilers?

It’d be nice to see the comparative times added for Crystal! :slight_smile:

1 Like

I’d like to but I don’t have time for that. Also, the code is not idiomatic (like, there are returns when they are not needed.) It also outputs a png instead of a bmp.

But feel free to do so.

If you don’t mind, I can submit it (since I initiated this thread).
There was no requirement to output a *.bmp, as at least one other did a *.png too.

What I can do is first create a gist of the code, and clean it up a little (there were two identical implementations of the times method in the Ruby version, when only on was needed, which I removed before running, which can also be removed from the Crystal version). Everyone can then run it, and comment and suggest improvements before we submit it for inclusion in their list.

I’ve run about 5 of the other implementations from the list on my laptop to establish a baseline for the times. Now that the Crystal version is done I’ll post the list of times from my system to create a common point of reference for all the versions. And I can do it with the just released Crystal 1.0!. :smile:

6 Likes

I’ve run a number of the different languages, per their compilation *.bat files, with the results on my system, which creates a reference baseline going forward. You can compare the time differences shown on the website.

Times shown are fastest among at least 5 - 10 runs, depending on their length.

System: Linux 5.10.25, GCC 10.2, LLVM 10.0.1
Hardware: i7-6700HQ (2016), 3.5 GHz, 4C|8T

language (version)  | times (ms) | file |
--------------------|------------|------|
Nim   (1.44) [fast] |     42.48  |  bmp |
C     (g++ 10.2)    |     75.0   |  bmp |
C++   (g++ 10.2)    |     78.0   |  bmp |
Rust  (1.50)        |     84.92  |  bmp |
Crystal (1.0.0)     |    106.76  |  png |
Nim   (1.44) [orig] |    118.74  |  bmp |
D     (ldc2 1.25.1) |    160.0   |  bmp |
Julia (1.6.0)       |    179.31  |  bmp |
Go    (1.16)        |    244.98  |  png |
JRuby (9.2.16.0)    |  13059.19  |  bmp |
Ruby  (2.6.6)       |  13769.55  |  bmp |
Ruby  (2.7.2)       |  15816.52  |  bmp |
Python (3.6.5)      |  27000.00  |  png |

Notes:
The C version didn’t create a valid *.bmp file (no picture shows).
The Rust picture had different shading (darker) from the rest.
Someone submitted a highly optimized (cheating) version for Nim 6 days ago.
Luckily, I had the old (original) code, and present both time results.
The Crystal version is @asterite’s original posted here.
Now we can play with it to see if it can be made comparatively faster.

EDIT:
2021/3/25: Replaced Julia 1.5.4|233 ms, with new Julia 1.6.0|179 ms.

9 Likes

Performance-wise, Crystal is looking good! I see a mix of int’s and float’s in the ported code. I wonder if going all float (where currently mixed, thus cutting out auto-conversions) would help performance a bit more?

How were you running Ruby, with JIT on or off?

All the Ruby’s were run without jits.

I have 3.0, via rvm, but it has an install bug that causes it to not access gems, so it can’t get the 2 gems needed to run the program. Truffleruby also can’t use those gems so I can’t run it either.

Hey @asterite, in reviewing your code you used abstract classes.
I found this forum thread about them, but am not clear why you used them.
Were they necessary (is there another way to do it?), or was it just easier using them. I’m trying to understand your code before I try to see if I can make it faster.

We have two types of Surface and we need to store those in an instance variable. We can use a union or an abstract type. I chose an abstract type because supposedly there could be more Surface types.

Also, using unions or abstract types makes no difference in regards to performance.

One thing I would try to look in the snippet is avoiding creating procs. In a couple of places Ruby was creating procs, and at least one of them was a closure. That requires a heap allocation. You can try removing those procs and inlining things (they are really pointless in this case) and see if things speed up a bit.

Thanks, that’s a little bit clearer about the abstract issues, though I still need to read/study more about them to become really clear.

I also wondered about the use of lambdas (Ruby) and procs (Crystal). I think we can get rid of them too. I’ll try.

Quick (working) simplification so far.

Original

struct Color
  getter r, g, b

  def initialize(@r : Float64, @g : Float64, @b : Float64)
  end

  def self.scale(k, v)
    Color.new(k * v.r, k * v.g, k * v.b)
  end

  def self.plus(v1, v2)
    Color.new(v1.r + v2.r, v1.g + v2.g, v1.b + v2.b)
  end

  def self.times(v1, v2)
    Color.new(v1.r * v2.r, v1.g * v2.g, v1.b * v2.b)
  end

  def self.toDrawingColor(c)
    clamp = ->(d : Float64) do
      return 1 if d > 1
      return 0 if d < 0
      return d
    end
    r = (clamp.call(c.r)*255).floor
    g = (clamp.call(c.g)*255).floor
    b = (clamp.call(c.b)*255).floor
    {r, g, b}
  end
end

Simplification

struct Color
  getter r, g, b

  def initialize(@r : Float64, @g : Float64, @b : Float64) end

  def self.scale(k, v);   Color.new(k * v.r, k * v.g, k * v.b) end

  def self.plus(v1, v2);  Color.new(v1.r + v2.r, v1.g + v2.g, v1.b + v2.b) end

  def self.times(v1, v2); Color.new(v1.r * v2.r, v1.g * v2.g, v1.b * v2.b) end

  def self.toDrawingColor(c)
    clamp = ->(d : Float64) do
      return 255 if d > 1
      return   0 if d < 0
      (d * 255).floor
    end
    
    {clamp.call(c.r), clamp.call(c.g), clamp.call(c.b)}
  end
end