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