-
Notifications
You must be signed in to change notification settings - Fork 2
/
aobench-mpl.sml
313 lines (289 loc) · 9.83 KB
/
aobench-mpl.sml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
(**
* Parallelized AObench with parallel-for
*
* MaPLe (an implementation of ParallelML)
* $ mpl aobench-mpl.mlb
*
* SML# 3.5.0
* $ make -f makefile-smlsharp
*)
structure AObench :
sig
val main : (string * string list) -> int
end
=
struct
type range = { begin : int, limit : int, step : int -> int }
fun for {begin,limit,step} f =
let
fun for' i = if i<limit
then
(f i;
for' (step i))
else
()
in
for' begin
end
fun succ x = x + 1
val (WIDTH, HEIGHT) = (256, 256)
val (NSUBSAMPLES, NAO_SAMPLES) = (2, 8)
val img = Word8Array.array (WIDTH * HEIGHT * 3, 0w0)
type vec = { x : real, y : real, z : real }
type sphere = { center : vec
, radius : real }
type plane = { p : vec, n : vec }
type 'a triple = 'a * 'a * 'a
type ray = { org : vec, dir : vec }
type isect = { t : real
, p : vec
, n : vec
, hit : bool
}
fun vec_to_string {x,y,z} =
let
val real = Real.fmt (StringCvt.FIX (SOME 3))
in
concat["vec={", real x, ",", real y, ",", real z, "}"]
end
fun ray_to_string {org, dir} =
concat ["ray=(", vec_to_string org, ",", vec_to_string dir, ")"]
fun isect_to_string {t,p,n,hit} =
let
val real = Real.fmt (StringCvt.FIX (SOME 3))
in
concat["isect{t=", real t, " p=", vec_to_string p, ", n=", vec_to_string n, " hit=", Bool.toString hit, "}"]
end
fun init_scene () : sphere triple * plane =
let
val ss = ({ center={ x= ~2.0, y= 0.0, z= ~3.5 }, radius= 0.5 }
,{ center={ x= ~0.5, y= 0.0, z= ~3.0 }, radius= 0.5 }
,{ center={ x= 1.0, y= 0.0, z= ~2.2 }, radius= 0.5 }
)
val plane = { p={ x= 0.0, y= ~0.5, z= 0.0 }
, n={ x= 0.0, y= 1.0, z= 0.0 }
}
in
(ss, plane)
end
infix 7 <*>
fun (v0:vec) <*> (v1:vec) : real =
(#x v0 * #x v1) + (#y v0 * #y v1) + (#z v0 * #z v1)
fun vcross (v0:vec, v1:vec) : vec =
{ x= #y v0 * #z v1 - #z v0 * #y v1
, y= #z v0 * #x v1 - #x v0 * #z v1
, z= #x v0 * #y v1 - #y v0 * #x v1
}
fun vnormalize (c as {x,y,z}:vec) =
let
val length = Math.sqrt (c <*> c)
in
if abs length > 1.0e~17
then
{ x= x / length
, y= y / length
, z= z / length
}
else
c
end
fun ray_sphere_intersect isect ({org,dir}:ray) ({center,radius}:sphere) =
let
val rs = { x= #x org - #x center
, y= #y org - #y center
, z= #z org - #z center
}
val B = rs <*> dir
val C = rs <*> rs - radius * radius
val D = B * B - C
in
if D > 0.0
then
let val t = ~B - Math.sqrt D
in
if t > 0.0 andalso t < #t isect
then
let
val p = { x= #x org + #x dir * t
, y= #y org + #y dir * t
, z= #z org + #z dir * t }
val n = { x= #x p - #x center
, y= #y p - #y center
, z= #z p - #z center }
in
{ t= t
, p= p
, n= vnormalize n
, hit= true
}
end
else
isect
end
else
isect
end
fun ray_plane_intersect isect ({dir,org}:ray) ({n,p}:plane) =
let
val d = Real.~ (p <*> n)
val v = dir <*> n
in
if abs v < 1.0e~17 then isect
else
let
val t = Real.~ ((org <*> n) + d) / v
in
if t > 0.0 andalso t < #t isect
then
{ t= t
, p= { x = #x org + #x dir * t
, y = #y org + #y dir * t
, z = #z org + #z dir * t }
, n= n
, hit= true
}
else
isect
end
end
fun orthoBasis (n as {x,y,z}) =
let
val basis2 = n
val basis1 = if x < 0.6 andalso x > ~0.6 then
{ x= 1.0, y= 0.0, z= 0.0 }
else if y < 0.6 andalso y > ~0.6 then
{ x= 0.0, y= 1.0, z= 0.0 }
else if z < 0.6 andalso z > ~0.6 then
{ x= 0.0, y= 0.0, z= 1.0 }
else
{ x= 1.0, y= 0.0, z= 0.0 }
val basis0 = vnormalize (vcross (basis1, basis2))
val basis1 = vnormalize (vcross (basis2, basis0))
in
(basis0, basis1, basis2)
end
fun mk_rand seed =
let
val rand = Random.rand (seed, Option.getOpt (Int.maxInt, 1073741823))
in
fn () => Random.randReal rand
end
fun ambient_occlusion drand48 ({p,n,...}:isect) (spheres:sphere triple) plane =
let
val ntheta = NAO_SAMPLES
val nphi = NAO_SAMPLES
val eps = 0.0001
val p = { x= #x p + eps * #x n
, y= #y p + eps * #y n
, z= #z p + eps * #z n
}
val (basis0, basis1, basis2) = orthoBasis n
val occlusion = ref 0.0
in
for {begin=0, limit=ntheta, step=succ} (fn j=>
for {begin=0, limit=nphi, step=succ} (fn i=>
let
val theta = Math.sqrt(drand48())
val phi = 2.0 * Math.pi * drand48()
val (x,y,z) = ( (Math.cos phi) * theta
, (Math.sin phi) * theta
, Math.sqrt (1.0 - theta * theta))
val rx = x * (#x basis0) + y * (#x basis1) + z * (#x basis2)
val ry = x * (#y basis0) + y * (#y basis1) + z * (#y basis2)
val rz = x * (#z basis0) + y * (#z basis1) + z * (#z basis2)
val zero = { x=0.0, y=0.0, z=0.0 }
val ray = { org= p, dir= { x= rx, y= ry, z= rz } }
val occIsect = { t= 1.0e17, p=zero, n=zero, hit= false }
val occIsect = ray_sphere_intersect occIsect ray (#1 spheres)
val occIsect = ray_sphere_intersect occIsect ray (#2 spheres)
val occIsect = ray_sphere_intersect occIsect ray (#3 spheres)
val occIsect = ray_plane_intersect occIsect ray plane
in
occlusion := !occlusion + (if #hit occIsect then 1.0 else 0.0)
end
));
let val occlusion = (real (ntheta * nphi) - !occlusion)
/ real (ntheta * nphi)
in
{ x=occlusion, y=occlusion, z=occlusion }
end
end
fun clamp f =
let val i = Real.floor (f * 255.5)
in
Word8.fromInt (Int.min (255, Int.max (i, 0)))
end
fun render grain img w h nsubsamples (spheres:sphere triple) plane =
let
val fimg = RealArray.array(w*h*3, 0.0)
in
ForkJoin.parfor grain (0, h) (fn y =>
let val drand48 = mk_rand (48271 + y) in
for {begin=0, limit=w, step=succ} (fn x=>
(
for {begin=0, limit=nsubsamples, step=succ} (fn v=>
for {begin=0, limit=nsubsamples, step=succ} (fn u=>
let
val zero = { x=0.0, y=0.0, z= 0.0 }
val px = (real x + (real u / real nsubsamples) - (real w/2.0)) / (real w/2.0)
val py = ~(real y + (real v / real nsubsamples) - (real h/2.0)) / (real h/2.0)
val ray = { org=zero
, dir=vnormalize { x=px, y=py, z= ~1.0 }
}
val isect = { t= 1.0e17, p=zero, n=zero, hit= false }
val isect = ray_sphere_intersect isect ray (#1 spheres)
val isect = ray_sphere_intersect isect ray (#2 spheres)
val isect = ray_sphere_intersect isect ray (#3 spheres)
val isect = ray_plane_intersect isect ray plane
in
if #hit isect
then
let
val col = ambient_occlusion drand48 isect spheres plane
val real = Real.fmt (StringCvt.FIX (SOME 3))
in
RealArray.update (fimg, 3*(y*w+x)+0, RealArray.sub(fimg, 3*(y*w+x)+0) + #x col);
RealArray.update (fimg, 3*(y*w+x)+1, RealArray.sub(fimg, 3*(y*w+x)+1) + #y col);
RealArray.update (fimg, 3*(y*w+x)+2, RealArray.sub(fimg, 3*(y*w+x)+2) + #z col)
end
else
()
end
));
RealArray.update (fimg, 3*(y*w+x)+0, RealArray.sub(fimg, 3*(y*w+x)+0) / real (nsubsamples * nsubsamples));
RealArray.update (fimg, 3*(y*w+x)+1, RealArray.sub(fimg, 3*(y*w+x)+1) / real (nsubsamples * nsubsamples));
RealArray.update (fimg, 3*(y*w+x)+2, RealArray.sub(fimg, 3*(y*w+x)+2) / real (nsubsamples * nsubsamples));
Word8Array.update (img, 3*(y*w+x)+0, clamp (RealArray.sub(fimg, 3*(y*w+x)+0)));
Word8Array.update (img, 3*(y*w+x)+1, clamp (RealArray.sub(fimg, 3*(y*w+x)+1)));
Word8Array.update (img, 3*(y*w+x)+2, clamp (RealArray.sub(fimg, 3*(y*w+x)+2)))
))
end
)
end
fun saveppm fname w h img =
let
val fp = BinIO.openOut fname
fun write_str fp ss = Substring.app (fn c=> BinIO.output1 (fp, Word8.fromInt (ord c))) (Substring.full ss)
in
write_str fp "P6\n";
write_str fp (Int.toString w^" "^Int.toString h^"\n");
write_str fp "255\n";
BinIO.output (fp, Word8Array.vector img);
BinIO.closeOut fp
end
fun main (name,args) =
let
val _ = print "init_scene...\n"
val (spheres,plane) = init_scene()
val grain = case args
of _::grain::args' => valOf (Int.fromString grain)
| _ => 1 (* fully parallel *)
in
print (concat ["rendering(grain size: ", Int.toString grain, ")...\n"]);
render grain img WIDTH HEIGHT NSUBSAMPLES spheres plane;
saveppm
(case args of file::_ => file | _=> "aosml.ppm")
WIDTH HEIGHT img;
0
end
end