ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/gen/mkillum2.c
Revision: 2.38
Committed: Sat Oct 13 20:15:43 2012 UTC (11 years, 6 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad4R2, rad4R2P1
Changes since 2.37: +51 -327 lines
Log Message:
Corrected errors in XML interpreter and genBSDF and removed mkillum BSDF code

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: mkillum2.c,v 2.37 2011/08/16 18:09:53 greg Exp $";
3 #endif
4 /*
5 * Routines to do the actual calculation for mkillum
6 */
7
8 #include <string.h>
9
10 #include "mkillum.h"
11 #include "face.h"
12 #include "cone.h"
13 #include "source.h"
14 #include "paths.h"
15
16 #ifndef NBSDFSAMPS
17 #define NBSDFSAMPS 256 /* BSDF resampling count */
18 #endif
19
20 COLORV * distarr = NULL; /* distribution array */
21 int distsiz = 0;
22
23
24 void
25 newdist( /* allocate & clear distribution array */
26 int siz
27 )
28 {
29 if (siz <= 0) {
30 if (distsiz > 0)
31 free(distarr);
32 distarr = NULL;
33 distsiz = 0;
34 return;
35 }
36 if (distsiz < siz) {
37 if (distsiz > 0)
38 free(distarr);
39 distarr = (COLORV *)malloc(sizeof(COLOR)*siz);
40 if (distarr == NULL)
41 error(SYSTEM, "out of memory in newdist");
42 distsiz = siz;
43 }
44 memset(distarr, '\0', sizeof(COLOR)*siz);
45 }
46
47
48 int
49 process_ray( /* process a ray result or report error */
50 RAY *r,
51 int rv
52 )
53 {
54 COLORV *colp;
55
56 if (rv == 0) /* no result ready */
57 return(0);
58 if (rv < 0)
59 error(USER, "ray tracing process died");
60 if (r->rno >= distsiz)
61 error(INTERNAL, "bad returned index in process_ray");
62 multcolor(r->rcol, r->rcoef); /* in case it's a source ray */
63 colp = &distarr[r->rno * 3];
64 addcolor(colp, r->rcol);
65 return(1);
66 }
67
68
69 void
70 raysamp( /* queue a ray sample */
71 int ndx,
72 FVECT org,
73 FVECT dir
74 )
75 {
76 RAY myRay;
77 int rv;
78
79 if ((ndx < 0) | (ndx >= distsiz))
80 error(INTERNAL, "bad index in raysamp");
81 VCOPY(myRay.rorg, org);
82 VCOPY(myRay.rdir, dir);
83 myRay.rmax = .0;
84 rayorigin(&myRay, PRIMARY|SPECULAR, NULL, NULL);
85 myRay.rno = ndx;
86 /* queue ray, check result */
87 process_ray(&myRay, ray_pqueue(&myRay));
88 }
89
90
91 void
92 srcsamps( /* sample sources from this surface position */
93 struct illum_args *il,
94 FVECT org,
95 FVECT nrm,
96 MAT4 ixfm
97 )
98 {
99 int nalt=1, nazi=1;
100 SRCINDEX si;
101 RAY sr;
102 FVECT v;
103 double d;
104 int i, j;
105 /* get sampling density */
106 if (il->sampdens > 0) {
107 i = PI * il->sampdens;
108 nalt = sqrt(i/PI) + .5;
109 nazi = PI*nalt + .5;
110 }
111 initsrcindex(&si); /* loop over (sub)sources */
112 for ( ; ; ) {
113 VCOPY(sr.rorg, org); /* pick side to shoot from */
114 d = 5.*FTINY;
115 VSUM(sr.rorg, sr.rorg, nrm, d);
116 samplendx++; /* increment sample counter */
117 if (!srcray(&sr, NULL, &si))
118 break; /* end of sources */
119 /* index direction */
120 if (ixfm != NULL)
121 multv3(v, sr.rdir, ixfm);
122 else
123 VCOPY(v, sr.rdir);
124 if (v[2] >= -FTINY)
125 continue; /* only sample transmission */
126 v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2];
127 sr.rno = flatindex(v, nalt, nazi);
128 d = nalt*nazi*(1./PI) * v[2];
129 d *= si.dom; /* solid angle correction */
130 scalecolor(sr.rcoef, d);
131 process_ray(&sr, ray_pqueue(&sr));
132 }
133 }
134
135
136 void
137 rayclean() /* finish all pending rays */
138 {
139 RAY myRay;
140
141 while (process_ray(&myRay, ray_presult(&myRay, 0)))
142 ;
143 }
144
145
146 static void
147 mkaxes( /* compute u and v to go with n */
148 FVECT u,
149 FVECT v,
150 FVECT n
151 )
152 {
153 register int i;
154
155 v[0] = v[1] = v[2] = 0.0;
156 for (i = 0; i < 3; i++)
157 if (n[i] < 0.6 && n[i] > -0.6)
158 break;
159 v[i] = 1.0;
160 fcross(u, v, n);
161 normalize(u);
162 fcross(v, n, u);
163 }
164
165
166 static void
167 rounddir( /* compute uniform spherical direction */
168 register FVECT dv,
169 double alt,
170 double azi
171 )
172 {
173 double d1, d2;
174
175 dv[2] = 1. - 2.*alt;
176 d1 = sqrt(1. - dv[2]*dv[2]);
177 d2 = 2.*PI * azi;
178 dv[0] = d1*cos(d2);
179 dv[1] = d1*sin(d2);
180 }
181
182
183 void
184 flatdir( /* compute uniform hemispherical direction */
185 FVECT dv,
186 double alt,
187 double azi
188 )
189 {
190 double d1, d2;
191
192 d1 = sqrt(alt);
193 d2 = 2.*PI * azi;
194 dv[0] = d1*cos(d2);
195 dv[1] = d1*sin(d2);
196 dv[2] = sqrt(1. - alt);
197 }
198
199
200 int
201 flatindex( /* compute index for hemispherical direction */
202 FVECT dv,
203 int nalt,
204 int nazi
205 )
206 {
207 double d;
208 int i, j;
209
210 d = 1.0 - dv[2]*dv[2];
211 i = d*nalt;
212 d = atan2(dv[1], dv[0]) * (0.5/PI);
213 if (d < 0.0) d += 1.0;
214 j = d*nazi + 0.5;
215 if (j >= nazi) j = 0;
216 return(i*nazi + j);
217 }
218
219
220 int
221 my_default( /* default illum action */
222 OBJREC *ob,
223 struct illum_args *il,
224 char *nm
225 )
226 {
227 sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
228 nm, ofun[ob->otype].funame, ob->oname);
229 error(WARNING, errmsg);
230 printobj(il->altmat, ob);
231 return(1);
232 }
233
234
235 int
236 my_face( /* make an illum face */
237 OBJREC *ob,
238 struct illum_args *il,
239 char *nm
240 )
241 {
242 int dim[2];
243 int n, nalt, nazi, alti;
244 double sp[2], r1, r2;
245 int h;
246 FVECT dn, org, dir;
247 FVECT u, v;
248 double ur[2], vr[2];
249 MAT4 xfm;
250 char xfrot[64];
251 int nallow;
252 FACE *fa;
253 int i, j;
254 /* get/check arguments */
255 fa = getface(ob);
256 if (fa->area == 0.0) {
257 freeface(ob);
258 return(my_default(ob, il, nm));
259 }
260 /* set up sampling */
261 if (il->sampdens <= 0) {
262 nalt = nazi = 1; /* diffuse assumption */
263 } else {
264 n = PI * il->sampdens;
265 nalt = sqrt(n/PI) + .5;
266 nazi = PI*nalt + .5;
267 }
268 n = nazi*nalt;
269 newdist(n);
270 /* take first edge >= sqrt(area) */
271 for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
272 u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
273 u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
274 u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
275 if ((r1 = DOT(u,u)) >= fa->area-FTINY)
276 break;
277 }
278 if (i < fa->nv) { /* got one! -- let's align our axes */
279 r2 = 1.0/sqrt(r1);
280 u[0] *= r2; u[1] *= r2; u[2] *= r2;
281 fcross(v, fa->norm, u);
282 } else /* oh well, we'll just have to wing it */
283 mkaxes(u, v, fa->norm);
284 /* now, find limits in (u,v) coordinates */
285 ur[0] = vr[0] = FHUGE;
286 ur[1] = vr[1] = -FHUGE;
287 for (i = 0; i < fa->nv; i++) {
288 r1 = DOT(VERTEX(fa,i),u);
289 if (r1 < ur[0]) ur[0] = r1;
290 if (r1 > ur[1]) ur[1] = r1;
291 r2 = DOT(VERTEX(fa,i),v);
292 if (r2 < vr[0]) vr[0] = r2;
293 if (r2 > vr[1]) vr[1] = r2;
294 }
295 dim[0] = random();
296 /* sample polygon */
297 nallow = 5*n*il->nsamps;
298 for (dim[1] = 0; dim[1] < n; dim[1]++)
299 for (i = 0; i < il->nsamps; i++) {
300 /* randomize direction */
301 h = ilhash(dim, 2) + i;
302 multisamp(sp, 2, urand(h));
303 alti = dim[1]/nazi;
304 r1 = (alti + sp[0])/nalt;
305 r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
306 flatdir(dn, r1, r2);
307 for (j = 0; j < 3; j++)
308 dir[j] = -dn[0]*u[j] - dn[1]*v[j] -
309 dn[2]*fa->norm[j];
310 /* randomize location */
311 do {
312 multisamp(sp, 2, urand(h+4862+nallow));
313 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
314 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
315 for (j = 0; j < 3; j++)
316 org[j] = r1*u[j] + r2*v[j]
317 + fa->offset*fa->norm[j];
318 } while (!inface(org, fa) && nallow-- > 0);
319 if (nallow < 0) {
320 objerror(ob, WARNING, "bad aspect");
321 rayclean();
322 freeface(ob);
323 return(my_default(ob, il, nm));
324 }
325 r1 = 5.*FTINY;
326 for (j = 0; j < 3; j++)
327 org[j] += r1*fa->norm[j];
328 /* send sample */
329 raysamp(dim[1], org, dir);
330 }
331 /* add in direct component? */
332 if (il->flags & IL_LIGHT) {
333 MAT4 ixfm;
334 for (i = 3; i--; ) {
335 ixfm[i][0] = u[i];
336 ixfm[i][1] = v[i];
337 ixfm[i][2] = fa->norm[i];
338 ixfm[i][3] = 0.;
339 }
340 ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
341 ixfm[3][3] = 1.;
342 dim[0] = random();
343 nallow = 10*il->nsamps;
344 for (i = 0; i < il->nsamps; i++) {
345 /* randomize location */
346 h = dim[0] + samplendx++;
347 do {
348 multisamp(sp, 2, urand(h+nallow));
349 r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
350 r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
351 for (j = 0; j < 3; j++)
352 org[j] = r1*u[j] + r2*v[j]
353 + fa->offset*fa->norm[j];
354 } while (!inface(org, fa) && nallow-- > 0);
355 if (nallow < 0) {
356 objerror(ob, WARNING, "bad aspect");
357 rayclean();
358 freeface(ob);
359 return(my_default(ob, il, nm));
360 }
361 /* sample source rays */
362 srcsamps(il, org, fa->norm, ixfm);
363 }
364 }
365 /* wait for all rays to finish */
366 rayclean();
367 /* write out the face and its distribution */
368 if (average(il, distarr, n)) {
369 if (il->sampdens > 0)
370 flatout(il, distarr, nalt, nazi, u, v, fa->norm);
371 illumout(il, ob);
372 } else
373 printobj(il->altmat, ob);
374 /* clean up */
375 freeface(ob);
376 return(0);
377 }
378
379
380 int
381 my_sphere( /* make an illum sphere */
382 register OBJREC *ob,
383 struct illum_args *il,
384 char *nm
385 )
386 {
387 int dim[3];
388 int n, nalt, nazi;
389 double sp[4], r1, r2, r3;
390 FVECT org, dir;
391 FVECT u, v;
392 register int i, j;
393 /* check arguments */
394 if (ob->oargs.nfargs != 4)
395 objerror(ob, USER, "bad # of arguments");
396 /* set up sampling */
397 if (il->sampdens <= 0)
398 nalt = nazi = 1;
399 else {
400 n = 4.*PI * il->sampdens;
401 nalt = sqrt(2./PI*n) + .5;
402 nazi = PI/2.*nalt + .5;
403 }
404 n = nalt*nazi;
405 newdist(n);
406 dim[0] = random();
407 /* sample sphere */
408 for (dim[1] = 0; dim[1] < nalt; dim[1]++)
409 for (dim[2] = 0; dim[2] < nazi; dim[2]++)
410 for (i = 0; i < il->nsamps; i++) {
411 /* next sample point */
412 multisamp(sp, 4, urand(ilhash(dim,3)+i));
413 /* random direction */
414 r1 = (dim[1] + sp[0])/nalt;
415 r2 = (dim[2] + sp[1] - .5)/nazi;
416 rounddir(dir, r1, r2);
417 /* random location */
418 mkaxes(u, v, dir); /* yuck! */
419 r3 = sqrt(sp[2]);
420 r2 = 2.*PI*sp[3];
421 r1 = r3*ob->oargs.farg[3]*cos(r2);
422 r2 = r3*ob->oargs.farg[3]*sin(r2);
423 r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
424 for (j = 0; j < 3; j++) {
425 org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
426 r3*dir[j];
427 dir[j] = -dir[j];
428 }
429 /* send sample */
430 raysamp(dim[1]*nazi+dim[2], org, dir);
431 }
432 /* wait for all rays to finish */
433 rayclean();
434 /* write out the sphere and its distribution */
435 if (average(il, distarr, n)) {
436 if (il->sampdens > 0)
437 roundout(il, distarr, nalt, nazi);
438 else
439 objerror(ob, WARNING, "diffuse distribution");
440 illumout(il, ob);
441 } else
442 printobj(il->altmat, ob);
443 /* clean up */
444 return(1);
445 }
446
447
448 int
449 my_ring( /* make an illum ring */
450 OBJREC *ob,
451 struct illum_args *il,
452 char *nm
453 )
454 {
455 int dim[2];
456 int n, nalt, nazi, alti;
457 double sp[2], r1, r2, r3;
458 int h;
459 FVECT dn, org, dir;
460 FVECT u, v;
461 MAT4 xfm;
462 CONE *co;
463 int i, j;
464 /* get/check arguments */
465 co = getcone(ob, 0);
466 /* set up sampling */
467 if (il->sampdens <= 0) {
468 nalt = nazi = 1; /* diffuse assumption */
469 } else {
470 n = PI * il->sampdens;
471 nalt = sqrt(n/PI) + .5;
472 nazi = PI*nalt + .5;
473 }
474 n = nazi*nalt;
475 newdist(n);
476 mkaxes(u, v, co->ad);
477 dim[0] = random();
478 /* sample disk */
479 for (dim[1] = 0; dim[1] < n; dim[1]++)
480 for (i = 0; i < il->nsamps; i++) {
481 /* next sample point */
482 h = ilhash(dim,2) + i;
483 /* randomize direction */
484 multisamp(sp, 2, urand(h));
485 alti = dim[1]/nazi;
486 r1 = (alti + sp[0])/nalt;
487 r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
488 flatdir(dn, r1, r2);
489 for (j = 0; j < 3; j++)
490 dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
491 /* randomize location */
492 multisamp(sp, 2, urand(h+8371));
493 r3 = sqrt(CO_R0(co)*CO_R0(co) +
494 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
495 r2 = 2.*PI*sp[1];
496 r1 = r3*cos(r2);
497 r2 = r3*sin(r2);
498 r3 = 5.*FTINY;
499 for (j = 0; j < 3; j++)
500 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
501 r3*co->ad[j];
502 /* send sample */
503 raysamp(dim[1], org, dir);
504 }
505 /* add in direct component? */
506 if (il->flags & IL_LIGHT) {
507 MAT4 ixfm;
508 for (i = 3; i--; ) {
509 ixfm[i][0] = u[i];
510 ixfm[i][1] = v[i];
511 ixfm[i][2] = co->ad[i];
512 ixfm[i][3] = 0.;
513 }
514 ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
515 ixfm[3][3] = 1.;
516 dim[0] = random();
517 for (i = 0; i < il->nsamps; i++) {
518 /* randomize location */
519 h = dim[0] + samplendx++;
520 multisamp(sp, 2, urand(h));
521 r3 = sqrt(CO_R0(co)*CO_R0(co) +
522 sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
523 r2 = 2.*PI*sp[1];
524 r1 = r3*cos(r2);
525 r2 = r3*sin(r2);
526 for (j = 0; j < 3; j++)
527 org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j];
528 /* sample source rays */
529 srcsamps(il, org, co->ad, ixfm);
530 }
531 }
532 /* wait for all rays to finish */
533 rayclean();
534 /* write out the ring and its distribution */
535 if (average(il, distarr, n)) {
536 if (il->sampdens > 0)
537 flatout(il, distarr, nalt, nazi, u, v, co->ad);
538 illumout(il, ob);
539 } else
540 printobj(il->altmat, ob);
541 /* clean up */
542 freecone(ob);
543 return(1);
544 }