1# AO render benchmark
2# Original program (C) Syoyo Fujita in Javascript (and other languages)
3#      https://code.google.com/p/aobench/
4# Ruby(yarv2llvm) version by Hideki Miura
5# mruby version by Hideki Miura
6#
7
8IMAGE_WIDTH = 64
9IMAGE_HEIGHT = 64
10NSUBSAMPLES = 2
11NAO_SAMPLES = 8
12
13module Rand
14  # Use xorshift
15  @@x = 123456789
16  @@y = 362436069
17  @@z = 521288629
18  @@w = 88675123
19  BNUM = 1 << 29
20  BNUMF = BNUM.to_f
21  def self.rand
22    x = @@x
23    t = x ^ ((x & 0xfffff) << 11)
24    w = @@w
25    @@x, @@y, @@z = @@y, @@z, w
26    w = @@w = (w ^ (w >> 19) ^ (t ^ (t >> 8)))
27    (w % BNUM) / BNUMF
28  end
29end
30
31class Vec
32  def initialize(x, y, z)
33    @x = x
34    @y = y
35    @z = z
36  end
37
38  def x=(v); @x = v; end
39  def y=(v); @y = v; end
40  def z=(v); @z = v; end
41  def x; @x; end
42  def y; @y; end
43  def z; @z; end
44
45  def vadd(b)
46    Vec.new(@x + b.x, @y + b.y, @z + b.z)
47  end
48
49  def vsub(b)
50    Vec.new(@x - b.x, @y - b.y, @z - b.z)
51  end
52
53  def vcross(b)
54    Vec.new(@y * b.z - @z * b.y,
55            @z * b.x - @x * b.z,
56            @x * b.y - @y * b.x)
57  end
58
59  def vdot(b)
60    r = @x * b.x + @y * b.y + @z * b.z
61    r
62  end
63
64  def vlength
65    Math.sqrt(@x * @x + @y * @y + @z * @z)
66  end
67
68  def vnormalize
69    len = vlength
70    v = Vec.new(@x, @y, @z)
71    if len > 1.0e-17
72      v.x = v.x / len
73      v.y = v.y / len
74      v.z = v.z / len
75    end
76    v
77  end
78end
79
80
81class Sphere
82  def initialize(center, radius)
83    @center = center
84    @radius = radius
85  end
86
87  def center; @center; end
88  def radius; @radius; end
89
90  def intersect(ray, isect)
91    rs = ray.org.vsub(@center)
92    b = rs.vdot(ray.dir)
93    c = rs.vdot(rs) - (@radius * @radius)
94    d = b * b - c
95    if d > 0.0
96      t = - b - Math.sqrt(d)
97
98      if t > 0.0 and t < isect.t
99        isect.t = t
100        isect.hit = true
101        isect.pl = Vec.new(ray.org.x + ray.dir.x * t,
102                          ray.org.y + ray.dir.y * t,
103                          ray.org.z + ray.dir.z * t)
104        n = isect.pl.vsub(@center)
105        isect.n = n.vnormalize
106      end
107    end
108  end
109end
110
111class Plane
112  def initialize(p, n)
113    @p = p
114    @n = n
115  end
116
117  def intersect(ray, isect)
118    d = -@p.vdot(@n)
119    v = ray.dir.vdot(@n)
120    v0 = v
121    if v < 0.0
122      v0 = -v
123    end
124    if v0 < 1.0e-17
125      return
126    end
127
128    t = -(ray.org.vdot(@n) + d) / v
129
130    if t > 0.0 and t < isect.t
131      isect.hit = true
132      isect.t = t
133      isect.n = @n
134      isect.pl = Vec.new(ray.org.x + t * ray.dir.x,
135                        ray.org.y + t * ray.dir.y,
136                        ray.org.z + t * ray.dir.z)
137    end
138  end
139end
140
141class Ray
142  def initialize(org, dir)
143    @org = org
144    @dir = dir
145  end
146
147  def org; @org; end
148  def org=(v); @org = v; end
149  def dir; @dir; end
150  def dir=(v); @dir = v; end
151end
152
153class Isect
154  def initialize
155    @t = 10000000.0
156    @hit = false
157    @pl = Vec.new(0.0, 0.0, 0.0)
158    @n = Vec.new(0.0, 0.0, 0.0)
159  end
160
161  def t; @t; end
162  def t=(v); @t = v; end
163  def hit; @hit; end
164  def hit=(v); @hit = v; end
165  def pl; @pl; end
166  def pl=(v); @pl = v; end
167  def n; @n; end
168  def n=(v); @n = v; end
169end
170
171def clamp(f)
172  i = f * 255.5
173  if i > 255.0
174    i = 255.0
175  end
176  if i < 0.0
177    i = 0.0
178  end
179  i.to_i
180end
181
182def otherBasis(basis, n)
183  basis[2] = Vec.new(n.x, n.y, n.z)
184  basis[1] = Vec.new(0.0, 0.0, 0.0)
185
186  if n.x < 0.6 and n.x > -0.6
187    basis[1].x = 1.0
188  elsif n.y < 0.6 and n.y > -0.6
189    basis[1].y = 1.0
190  elsif n.z < 0.6 and n.z > -0.6
191    basis[1].z = 1.0
192  else
193    basis[1].x = 1.0
194  end
195
196  basis[0] = basis[1].vcross(basis[2])
197  basis[0] = basis[0].vnormalize
198
199  basis[1] = basis[2].vcross(basis[0])
200  basis[1] = basis[1].vnormalize
201end
202
203class Scene
204  def initialize
205    @spheres = Array.new
206    @spheres[0] = Sphere.new(Vec.new(-2.0, 0.0, -3.5), 0.5)
207    @spheres[1] = Sphere.new(Vec.new(-0.5, 0.0, -3.0), 0.5)
208    @spheres[2] = Sphere.new(Vec.new(1.0, 0.0, -2.2), 0.5)
209    @plane = Plane.new(Vec.new(0.0, -0.5, 0.0), Vec.new(0.0, 1.0, 0.0))
210  end
211
212  def ambient_occlusion(isect)
213    basis = Array.new(3)
214    otherBasis(basis, isect.n)
215
216    ntheta    = NAO_SAMPLES
217    nphi      = NAO_SAMPLES
218    eps       = 0.0001
219    occlusion = 0.0
220
221    p0 = Vec.new(isect.pl.x + eps * isect.n.x,
222                isect.pl.y + eps * isect.n.y,
223                isect.pl.z + eps * isect.n.z)
224    nphi.times do
225      ntheta.times do
226        r = Rand::rand
227        phi = 2.0 * 3.14159265 * Rand::rand
228        x = Math.cos(phi) * Math.sqrt(1.0 - r)
229        y = Math.sin(phi) * Math.sqrt(1.0 - r)
230        z = Math.sqrt(r)
231
232        rx = x * basis[0].x + y * basis[1].x + z * basis[2].x
233        ry = x * basis[0].y + y * basis[1].y + z * basis[2].y
234        rz = x * basis[0].z + y * basis[1].z + z * basis[2].z
235
236        raydir = Vec.new(rx, ry, rz)
237        ray = Ray.new(p0, raydir)
238
239        occisect = Isect.new
240        @spheres[0].intersect(ray, occisect)
241        @spheres[1].intersect(ray, occisect)
242        @spheres[2].intersect(ray, occisect)
243        @plane.intersect(ray, occisect)
244        if occisect.hit
245          occlusion = occlusion + 1.0
246        else
247          0.0
248        end
249      end
250    end
251
252    occlusion = (ntheta.to_f * nphi.to_f - occlusion) / (ntheta.to_f * nphi.to_f)
253    Vec.new(occlusion, occlusion, occlusion)
254  end
255
256  def render(w, h, nsubsamples)
257    cnt = 0
258    nsf = nsubsamples.to_f
259    h.times do |y|
260      w.times do |x|
261        rad = Vec.new(0.0, 0.0, 0.0)
262
263        # Subsmpling
264        nsubsamples.times do |v|
265          nsubsamples.times do |u|
266            cnt = cnt + 1
267            wf = w.to_f
268            hf = h.to_f
269            xf = x.to_f
270            yf = y.to_f
271            uf = u.to_f
272            vf = v.to_f
273
274            px = (xf + (uf / nsf) - (wf / 2.0)) / (wf / 2.0)
275            py = -(yf + (vf / nsf) - (hf / 2.0)) / (hf / 2.0)
276
277            eye = Vec.new(px, py, -1.0).vnormalize
278
279            ray = Ray.new(Vec.new(0.0, 0.0, 0.0), eye)
280
281            isect = Isect.new
282            @spheres[0].intersect(ray, isect)
283            @spheres[1].intersect(ray, isect)
284            @spheres[2].intersect(ray, isect)
285            @plane.intersect(ray, isect)
286            if isect.hit
287              col = ambient_occlusion(isect)
288              rad.x = rad.x + col.x
289              rad.y = rad.y + col.y
290              rad.z = rad.z + col.z
291            else
292              0.0
293            end
294          end
295        end
296
297        r = rad.x / (nsf * nsf)
298        g = rad.y / (nsf * nsf)
299        b = rad.z / (nsf * nsf)
300        printf("%c", clamp(r))
301        printf("%c", clamp(g))
302        printf("%c", clamp(b))
303      end
304    end
305  end
306end
307
308# File.open("ao.ppm", "w") do |fp|
309  printf("P6\n")
310  printf("%d %d\n", IMAGE_WIDTH, IMAGE_HEIGHT)
311  printf("255\n", IMAGE_WIDTH, IMAGE_HEIGHT)
312  Scene.new.render(IMAGE_WIDTH, IMAGE_HEIGHT, NSUBSAMPLES)
313#  Scene.new.render(256, 256, 2)
314# end
315