?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
move the smallpt code under the smp folder
[AROS.git] / test / exec / smp / smp-smallpt / smallpt.cc
bloba99a1579bf5d7f070007b317e90a5845e6633224
1 #define DEBUG 1
2 #include <math.h> // smallpt, a Path Tracer by Kevin Beason, 2008
3 #include <stdlib.h>
4 #include <exec/types.h>
5 #include <aros/debug.h>
6 #include <exec/tasks.h>
7 #include <proto/exec.h>
8
9 #include "renderer.h"
10
11 struct Vec { // Usage: time ./smallpt 5000 && xv image.ppm
12 double x;
13 double y;
14 double z; // position, also color (r,g,b)
15
16 Vec(double x_ = 0, double y_ = 0, double z_ = 0)
17 {
18 x = x_;
19 y = y_;
20 z = z_;
21 }
22
23 Vec operator+(const Vec &b) const { return Vec(x + b.x, y + b.y, z + b.z); }
24 Vec operator-(const Vec &b) const { return Vec(x - b.x, y - b.y, z - b.z); }
25 Vec operator*(double b) const { return Vec(x * b, y * b, z * b); }
26 Vec mult(const Vec &b) const { return Vec(x * b.x, y * b.y, z * b.z); }
27 Vec& norm() { return *this = *this * (1/sqrt(x * x + y * y + z * z)); }
28 double dot(const Vec &b) const { return x * b.x + y * b.y + z * b.z; } // cross:
29 Vec operator%(Vec &b){ return Vec(y * b.z - z * b.y, z * b.x - x * b.z, x * b.y - y * b.x);}
30 };
31
32 struct Ray {
33 Vec o, d;
34 Ray(Vec o_, Vec d_) : o(o_), d(d_) {}
35 };
36
37 enum Refl_t
38 {
39 DIFF,
40 SPEC,
41 REFR
42 }; // material types, used in radiance()
43
44 struct Sphere
45 {
46 double rad; // radius
47 Vec p, e, c; // position, emission, color
48 Refl_t refl; // reflection type (DIFFuse, SPECular, REFRactive)
49 Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_) : rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
50 double intersect(const Ray &r) const
51 { // returns distance, 0 if nohit
52 Vec op = p - r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
53 double t, eps = 1e-4, b = op.dot(r.d), det = b * b - op.dot(op) + rad * rad;
54 if (det < 0)
55 return 0;
56 else
57 det = sqrt(det);
58 return (t = b - det) > eps ? t : ((t = b + det) > eps ? t : 0);
59 }
60 };
61
62 Sphere spheres[] = {//Scene: radius, position, emission, color, material
63 Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12), Vec(), DIFF), //Lite
64 Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left
65 Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght
66 Sphere(1e5, Vec(50,40.8, 1e5), Vec(),Vec(.55,.55,.55),DIFF),//Back
67 Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(), DIFF),//Frnt
68 Sphere(1e5, Vec(50, 1e5, 81.6), Vec(),Vec(.75,.75,.75),DIFF),//Botm
69 Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top
70 Sphere(16.5,Vec(27,16.5,47), Vec(),Vec(.4,.4,.3)*.999, SPEC),//Mirr
71 Sphere(16.5,Vec(73,16.5,78), Vec(),Vec(.8,.7,.95)*.999, REFR),//Glas
72 Sphere(10.5,Vec(23,10.5,98), Vec(),Vec(.6,1,0.7)*.999, REFR),//Glas
73 Sphere(8.,Vec(50,8.,108), Vec(),Vec(1,0.6,0.7)*.999, REFR),//Glas
74 Sphere(6.5, Vec(53,6.5,48), Vec(),Vec(0.3,.4,.4)*.999, SPEC),//Mirr
75 };
76
77 const int numSpheres = sizeof(spheres)/sizeof(Sphere);
78
79 inline double clamp(double x)
80 {
81 return x<0 ? 0 : x>1 ? 1 : x;
82 }
83
84 inline int toInt(double x)
85 {
86 return int(pow(clamp(x),1/2.2)*255+.5);
87 }
88
89 inline bool intersect(const Ray &r, double &t, int &id)
90 {
91 double n=numSpheres, d, inf=t=1e20;
92 for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;}
93 return t<inf;
94 }
95
96 int maximal_ray_depth = 1000;
97
98 // ca. 650bytes per ray depth, 650KB stack required for max ray depth of 1000
99
100 Vec radiance_expl(const Ray &r, int depth, unsigned short *Xi,int E=1){
101 double t; // distance to intersection
102 int id=0; // id of intersected object
103 if (!intersect(r, t, id)) return Vec(); // if miss, return black
104 const Sphere &obj = spheres[id]; // the hit object
105 Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
106 double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
107 depth++;
108 // If depth larger than maximal_ray_depth do not use Russian roulette, give up unconditionally
109 // because AROS does not have automatic stack expansion
110 if (depth > maximal_ray_depth)
111 return obj.e*E;
112 else if (depth>10||!p) { // From depth 10 start Russian roulette
113 if (erand48(Xi)<p) f=f*(1/p); else return obj.e*E;
114 }
115 if (obj.refl == DIFF){ // Ideal DIFFUSE reflection
116 double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
117 Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
118 Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
119
120 // Loop over any lights
121 Vec e;
122 for (int i=0; i<numSpheres; i++){
123 const Sphere &s = spheres[i];
124 if (s.e.x<=0 && s.e.y<=0 && s.e.z<=0) continue; // skip non-lights
125
126 Vec sw=s.p-x, su=((fabs(sw.x)>.1?Vec(0,1):Vec(1))%sw).norm(), sv=sw%su;
127 double cos_a_max = sqrt(1-s.rad*s.rad/(x-s.p).dot(x-s.p));
128 double eps1 = erand48(Xi), eps2 = erand48(Xi);
129 double cos_a = 1-eps1+eps1*cos_a_max;
130 double sin_a = sqrt(1-cos_a*cos_a);
131 double phi = 2*M_PI*eps2;
132 Vec l = su*cos(phi)*sin_a + sv*sin(phi)*sin_a + sw*cos_a;
133 l.norm();
134 if (intersect(Ray(x,l), t, id) && id==i){ // shadow ray
135 double omega = 2*M_PI*(1-cos_a_max);
136 e = e + f.mult(s.e*l.dot(nl)*omega)*M_1_PI; // 1/pi for brdf
137 }
138 }
139
140 return obj.e*E+e+f.mult(radiance_expl(Ray(x,d),depth,Xi,0));
141 } else if (obj.refl == SPEC) // Ideal SPECULAR reflection
142 return obj.e + f.mult(radiance_expl(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
143 Ray reflRay(x, r.d-n*2*n.dot(r.d)); // Ideal dielectric REFRACTION
144 bool into = n.dot(nl)>0; // Ray from outside going in?
145 double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
146 if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0) // Total internal reflection
147 return obj.e + f.mult(radiance_expl(reflRay,depth,Xi));
148 Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
149 double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
150 double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
151 return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ? // Russian roulette
152 radiance_expl(reflRay,depth,Xi)*RP:radiance_expl(Ray(x,tdir),depth,Xi)*TP) :
153 radiance_expl(reflRay,depth,Xi)*Re+radiance_expl(Ray(x,tdir),depth,Xi)*Tr);
154 }
155
156 Vec radiance(const Ray &r, int depth, unsigned short *Xi)
157 {
158 double t; // distance to intersection
159 int id=0; // id of intersected object
160
161 if (!intersect(r, t, id))
162 return Vec(); // if miss, return black
163
164 const Sphere &obj = spheres[id]; // the hit object
165
166 Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
167
168 depth++;
169 // Above maximal_ray_depth break recursive loop unconditionally
170 if (depth > maximal_ray_depth)
171 return obj.e;
172 else if (depth>5) // From depth of 5 start Russian roulette
173 {
174 double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
175
176 if (erand48(Xi)<p)
177 f=f*(1/p);
178 else return obj.e; //R.R.
179 }
180 if (obj.refl == DIFF){ // Ideal DIFFUSE reflection
181 double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
182 Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
183 Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
184 return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
185 } else if (obj.refl == SPEC) // Ideal SPECULAR reflection
186 return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
187 Ray reflRay(x, r.d-n*2*n.dot(r.d)); // Ideal dielectric REFRACTION
188 bool into = n.dot(nl)>0; // Ray from outside going in?
189 double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
190 if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0) // Total internal reflection
191 return obj.e + f.mult(radiance(reflRay,depth,Xi));
192 Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
193 double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
194 double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
195 return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ? // Russian roulette
196 radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) :
197 radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
198 }
199
200 static inline struct MyMessage *AllocMsg(struct MinList *msgPool)
201 {
202 if (IsMinListEmpty(msgPool))
203 bug("!!! smallpt.cc - run out of free messages!\n");
204
205 return (struct MyMessage *)REMHEAD(msgPool);
206 }
207 void __prepare()
208 {
209 ULONG *ptr = (ULONG *)FindTask(NULL)->tc_SPLower;
210 ULONG *rsp = 0;
211
212 asm volatile("mov %%rsp,%0":"=r"(rsp));
213
214 while (ptr < rsp-1)
215 *ptr++ = 0xdeadbeef;
216 }
217
218 void __test()
219 {
220 ULONG *ptr = (ULONG *)FindTask(NULL)->tc_SPLower;
221 IPTR top = (IPTR)FindTask(NULL)->tc_SPUpper;
222
223 while(*ptr++ == 0xdeadbeef);
224
225 bug("--> USED STACK: %d\n", top - (IPTR)ptr);
226 }
227
228 extern "C" void RenderTile(struct ExecBase *SysBase, struct MsgPort *masterPort, struct MsgPort **myPort)
229 {
230 Vec *c;
231 struct MyMessage *msg;
232 struct MyMessage *m;
233 struct MsgPort *port = CreateMsgPort();
234 struct MsgPort *syncPort = CreateMsgPort();
235 struct MinList msgPool;
236
237 c = (Vec *)AllocMem(sizeof(Vec) * TILE_SIZE * TILE_SIZE, MEMF_ANY | MEMF_CLEAR);
238
239 *myPort = port;
240
241 FreeSignal(syncPort->mp_SigBit);
242 syncPort->mp_SigBit = -1;
243 syncPort->mp_Flags = PA_IGNORE;
244
245 NEWLIST(&msgPool);
246
247 D(bug("[SMP-SmallPT-Task] hello, msgport=%p\n", port));
248
249 msg = (struct MyMessage *)AllocMem(sizeof(struct MyMessage) * 20, MEMF_PUBLIC | MEMF_CLEAR);
250 for (int i=0; i < 20; i++)
251 ADDHEAD(&msgPool, &msg[i]);
252
253 if (port)
254 {
255 ULONG signals;
256 BOOL doWork = TRUE;
257 BOOL redraw = TRUE;
258
259 m = AllocMsg(&msgPool);
260
261 /* Tell renderer that we are bored and want to do some work */
262 m->mm_Message.mn_Length = sizeof(msg);
263 m->mm_Message.mn_ReplyPort = port;
264 m->mm_Type = MSG_HUNGRY;
265 PutMsg(masterPort, &m->mm_Message);
266
267 D(bug("[SMP-SmallPT-Task] Just told renderer I'm hungry\n"));
268
269 do {
270 signals = Wait(SIGBREAKF_CTRL_C | (1 << port->mp_SigBit));
271
272 if (signals & (1 << port->mp_SigBit))
273 {
274 while ((m = (struct MyMessage *)GetMsg(port)))
275 {
276 if (m->mm_Message.mn_Node.ln_Type == NT_REPLYMSG)
277 {
278 ADDHEAD(&msgPool, m);
279 continue;
280 }
281 else
282 {
283 if (m->mm_Type == MSG_DIE)
284 {
285 doWork = FALSE;
286 ReplyMsg(&m->mm_Message);
287 }
288 else if (m->mm_Type == MSG_RENDERTILE)
289 {
290 struct tileWork *tile = m->mm_Body.RenderTile.tile;
291 int w = m->mm_Body.RenderTile.width;
292 int h = m->mm_Body.RenderTile.height;
293 int samps = m->mm_Body.RenderTile.numberOfSamples;
294 ULONG *buffer = m->mm_Body.RenderTile.buffer;
295 int tile_x = m->mm_Body.RenderTile.tile->x;
296 int tile_y = m->mm_Body.RenderTile.tile->y;
297 struct MsgPort *guiPort = m->mm_Body.RenderTile.guiPort;
298 int explicit_mode = m->mm_Body.RenderTile.explicitMode;
299
300 if (explicit_mode)
301 spheres[0] = Sphere(5, Vec(50,81.6-16.5,81.6),Vec(4,4,4)*20, Vec(), DIFF);
302 else
303 spheres[0] = Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12), Vec(), DIFF);
304
305 ReplyMsg(&m->mm_Message);
306
307 __prepare();
308
309 Ray cam(Vec(50, 52, 295.6), Vec(0, -0.042612, -1).norm()); // cam pos, dir
310 Vec cx = Vec(w * .5135 / h), cy = (cx % cam.d).norm() * .5135, r;
311
312 for (int i=0; i < TILE_SIZE * TILE_SIZE; i++)
313 c[i] = Vec();
314
315 for (int _y=tile_y * 32; _y < (tile_y + 1) * 32; _y++)
316 {
317 int y = h - _y - 1;
318 for (unsigned short _x=tile_x * 32, Xi[3]={0,0,(UWORD)(y*y*y)}; _x < (tile_x + 1) * 32; _x++) // Loop cols
319 {
320 int x = _x; // w - _x - 1;
321 for (int sy=0, i=(_y-tile_y*32)*32+_x-tile_x*32; sy<2; sy++) // 2x2 subpixel rows
322 {
323 for (int sx=0; sx<2; sx++, r=Vec())
324 { // 2x2 subpixel cols
325 for (int s=0; s<samps; s++)
326 {
327 double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1);
328 double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2);
329 Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) +
330 cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d;
331 if (explicit_mode)
332 r = r + radiance_expl(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
333 else
334 r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps);
335 } // Camera rays are pushed ^^^^^ forward to start in interior
336 c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
337 }
338 }
339 }
340 int start_ptr = tile_y*32*w + tile_x*32;
341 for (int yy=0; yy < 32; yy++)
342 {
343 for (int xx=0; xx < 32; xx++)
344 {
345 buffer[start_ptr+xx] = ((toInt(c[(xx+32*yy)].z) & 0xff) << 24) |
346 ((toInt(c[(xx+32*yy)].y) & 0xff) << 16) | ((toInt(c[(xx + 32*yy)].x) & 0xff) << 8) | 0xff;
347 }
348 start_ptr += w;
349 }
350
351 #if 1
352 if (redraw)
353 {
354 m = AllocMsg(&msgPool);
355 m->mm_Message.mn_Length = sizeof(struct MyMessage);
356 m->mm_Message.mn_ReplyPort = syncPort;
357 m->mm_Type = MSG_REDRAWTILE;
358 m->mm_Body.RedrawTile.TileX = tile_x;
359 m->mm_Body.RedrawTile.TileY = tile_y;
360 PutMsg(guiPort, &m->mm_Message);
361 redraw = FALSE;
362 }
363 else if ((m = (struct MyMessage *)GetMsg(syncPort)))
364 {
365 ADDHEAD(&msgPool, m);
366 redraw = TRUE;
367 }
368 #else
369 (void)syncPort;
370 Signal((struct Task *)guiPort->mp_SigTask, SIGBREAKF_CTRL_D);
371 #endif
372 }
373 __test();
374 Signal((struct Task *)guiPort->mp_SigTask, SIGBREAKF_CTRL_D);
375
376 m = AllocMsg(&msgPool);
377 m->mm_Message.mn_Length = sizeof(struct MyMessage);
378 m->mm_Message.mn_ReplyPort = port;
379 m->mm_Type = MSG_RENDERREADY;
380 m->mm_Body.RenderTile.tile = tile;
381 PutMsg(masterPort, &m->mm_Message);
382
383 m = AllocMsg(&msgPool);
384 m->mm_Message.mn_Length = sizeof(struct MyMessage);
385 m->mm_Message.mn_ReplyPort = port;
386 m->mm_Type = MSG_HUNGRY;
387 PutMsg(masterPort, &m->mm_Message);
388 }
389 }
390 }
391 }
392
393 if (signals & SIGBREAKB_CTRL_C)
394 doWork = FALSE;
395
396 } while(doWork);
397 }
398 D(bug("[SMP-SmallPT-Task] cleaning up stuff\n"));
399 FreeMem(msg, sizeof(struct MyMessage) * 20);
400 DeleteMsgPort(port);
401 DeleteMsgPort(syncPort);
402
403 FreeMem(c, sizeof(Vec) * TILE_SIZE * TILE_SIZE);
404 }
Mirror of AROS' svn.