--- ray/src/gen/mkillum2.c 2007/09/21 05:53:21 2.21 +++ ray/src/gen/mkillum2.c 2012/10/13 20:15:43 2.38 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: mkillum2.c,v 2.21 2007/09/21 05:53:21 greg Exp $"; +static const char RCSid[] = "$Id: mkillum2.c,v 2.38 2012/10/13 20:15:43 greg Exp $"; #endif /* * Routines to do the actual calculation for mkillum @@ -10,8 +10,12 @@ static const char RCSid[] = "$Id: mkillum2.c,v 2.21 20 #include "mkillum.h" #include "face.h" #include "cone.h" -#include "random.h" +#include "source.h" +#include "paths.h" +#ifndef NBSDFSAMPS +#define NBSDFSAMPS 256 /* BSDF resampling count */ +#endif COLORV * distarr = NULL; /* distribution array */ int distsiz = 0; @@ -22,16 +26,16 @@ newdist( /* allocate & clear distribution array */ int siz ) { - if (siz == 0) { + if (siz <= 0) { if (distsiz > 0) - free((void *)distarr); + free(distarr); distarr = NULL; distsiz = 0; return; } if (distsiz < siz) { if (distsiz > 0) - free((void *)distarr); + free(distarr); distarr = (COLORV *)malloc(sizeof(COLOR)*siz); if (distarr == NULL) error(SYSTEM, "out of memory in newdist"); @@ -42,16 +46,20 @@ newdist( /* allocate & clear distribution array */ int -process_ray(RAY *r, int rv) +process_ray( /* process a ray result or report error */ + RAY *r, + int rv +) { COLORV *colp; - if (rv == 0) + if (rv == 0) /* no result ready */ return(0); if (rv < 0) error(USER, "ray tracing process died"); if (r->rno >= distsiz) error(INTERNAL, "bad returned index in process_ray"); + multcolor(r->rcol, r->rcoef); /* in case it's a source ray */ colp = &distarr[r->rno * 3]; addcolor(colp, r->rcol); return(1); @@ -59,7 +67,7 @@ process_ray(RAY *r, int rv) void -raysamp( /* queue a ray sample */ +raysamp( /* queue a ray sample */ int ndx, FVECT org, FVECT dir @@ -73,7 +81,7 @@ raysamp( /* queue a ray sample */ VCOPY(myRay.rorg, org); VCOPY(myRay.rdir, dir); myRay.rmax = .0; - rayorigin(&myRay, PRIMARY, NULL, NULL); + rayorigin(&myRay, PRIMARY|SPECULAR, NULL, NULL); myRay.rno = ndx; /* queue ray, check result */ process_ray(&myRay, ray_pqueue(&myRay)); @@ -81,6 +89,51 @@ raysamp( /* queue a ray sample */ void +srcsamps( /* sample sources from this surface position */ + struct illum_args *il, + FVECT org, + FVECT nrm, + MAT4 ixfm +) +{ + int nalt=1, nazi=1; + SRCINDEX si; + RAY sr; + FVECT v; + double d; + int i, j; + /* get sampling density */ + if (il->sampdens > 0) { + i = PI * il->sampdens; + nalt = sqrt(i/PI) + .5; + nazi = PI*nalt + .5; + } + initsrcindex(&si); /* loop over (sub)sources */ + for ( ; ; ) { + VCOPY(sr.rorg, org); /* pick side to shoot from */ + d = 5.*FTINY; + VSUM(sr.rorg, sr.rorg, nrm, d); + samplendx++; /* increment sample counter */ + if (!srcray(&sr, NULL, &si)) + break; /* end of sources */ + /* index direction */ + if (ixfm != NULL) + multv3(v, sr.rdir, ixfm); + else + VCOPY(v, sr.rdir); + if (v[2] >= -FTINY) + continue; /* only sample transmission */ + v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2]; + sr.rno = flatindex(v, nalt, nazi); + d = nalt*nazi*(1./PI) * v[2]; + d *= si.dom; /* solid angle correction */ + scalecolor(sr.rcoef, d); + process_ray(&sr, ray_pqueue(&sr)); + } +} + + +void rayclean() /* finish all pending rays */ { RAY myRay; @@ -127,9 +180,9 @@ rounddir( /* compute uniform spherical direction */ } -static void +void flatdir( /* compute uniform hemispherical direction */ - register FVECT dv, + FVECT dv, double alt, double azi ) @@ -145,6 +198,26 @@ flatdir( /* compute uniform hemispherical direction * int +flatindex( /* compute index for hemispherical direction */ + FVECT dv, + int nalt, + int nazi +) +{ + double d; + int i, j; + + d = 1.0 - dv[2]*dv[2]; + i = d*nalt; + d = atan2(dv[1], dv[0]) * (0.5/PI); + if (d < 0.0) d += 1.0; + j = d*nazi + 0.5; + if (j >= nazi) j = 0; + return(i*nazi + j); +} + + +int my_default( /* default illum action */ OBJREC *ob, struct illum_args *il, @@ -166,17 +239,18 @@ my_face( /* make an illum face */ char *nm ) { -#define MAXMISS (5*n*il->nsamps) int dim[2]; - int n, nalt, nazi, h, alti; + int n, nalt, nazi, alti; double sp[2], r1, r2; + int h; FVECT dn, org, dir; FVECT u, v; double ur[2], vr[2]; MAT4 xfm; - int nmisses; + char xfrot[64]; + int nallow; FACE *fa; - register int i, j; + int i, j; /* get/check arguments */ fa = getface(ob); if (fa->area == 0.0) { @@ -191,15 +265,7 @@ my_face( /* make an illum face */ nalt = sqrt(n/PI) + .5; nazi = PI*nalt + .5; } - if (il->sd != NULL) { - if (!getBSDF_xfm(xfm, fa->norm, il->udir)) { - objerror(ob, WARNING, "illegal up direction"); - freeface(ob); - return(my_default(ob, il, nm)); - } - n = il->sd->ninc; - } else - n = nazi*nalt; + n = nazi*nalt; newdist(n); /* take first edge >= sqrt(area) */ for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) { @@ -228,50 +294,76 @@ my_face( /* make an illum face */ } dim[0] = random(); /* sample polygon */ - nmisses = 0; + nallow = 5*n*il->nsamps; for (dim[1] = 0; dim[1] < n; dim[1]++) for (i = 0; i < il->nsamps; i++) { - /* random direction */ + /* randomize direction */ h = ilhash(dim, 2) + i; - if (il->sd != NULL) { - r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm); - } else { - multisamp(sp, 2, urand(h)); - alti = dim[1]/nazi; - r1 = (alti + sp[0])/nalt; - r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi; - flatdir(dn, r1, r2); - for (j = 0; j < 3; j++) + multisamp(sp, 2, urand(h)); + alti = dim[1]/nazi; + r1 = (alti + sp[0])/nalt; + r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi; + flatdir(dn, r1, r2); + for (j = 0; j < 3; j++) dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*fa->norm[j]; - } - /* random location */ + /* randomize location */ do { - multisamp(sp, 2, urand(h+4862+nmisses)); + multisamp(sp, 2, urand(h+4862+nallow)); r1 = ur[0] + (ur[1]-ur[0]) * sp[0]; r2 = vr[0] + (vr[1]-vr[0]) * sp[1]; for (j = 0; j < 3; j++) org[j] = r1*u[j] + r2*v[j] + fa->offset*fa->norm[j]; - } while (!inface(org, fa) && nmisses++ < MAXMISS); - if (nmisses > MAXMISS) { + } while (!inface(org, fa) && nallow-- > 0); + if (nallow < 0) { objerror(ob, WARNING, "bad aspect"); rayclean(); freeface(ob); return(my_default(ob, il, nm)); } - if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY) - r1 = -1.0001*il->thick - .0001; - else - r1 = .0001; + r1 = 5.*FTINY; for (j = 0; j < 3; j++) org[j] += r1*fa->norm[j]; /* send sample */ raysamp(dim[1], org, dir); } + /* add in direct component? */ + if (il->flags & IL_LIGHT) { + MAT4 ixfm; + for (i = 3; i--; ) { + ixfm[i][0] = u[i]; + ixfm[i][1] = v[i]; + ixfm[i][2] = fa->norm[i]; + ixfm[i][3] = 0.; + } + ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.; + ixfm[3][3] = 1.; + dim[0] = random(); + nallow = 10*il->nsamps; + for (i = 0; i < il->nsamps; i++) { + /* randomize location */ + h = dim[0] + samplendx++; + do { + multisamp(sp, 2, urand(h+nallow)); + r1 = ur[0] + (ur[1]-ur[0]) * sp[0]; + r2 = vr[0] + (vr[1]-vr[0]) * sp[1]; + for (j = 0; j < 3; j++) + org[j] = r1*u[j] + r2*v[j] + + fa->offset*fa->norm[j]; + } while (!inface(org, fa) && nallow-- > 0); + if (nallow < 0) { + objerror(ob, WARNING, "bad aspect"); + rayclean(); + freeface(ob); + return(my_default(ob, il, nm)); + } + /* sample source rays */ + srcsamps(il, org, fa->norm, ixfm); + } + } + /* wait for all rays to finish */ rayclean(); - if (il->sd != NULL) /* run distribution through BSDF */ - redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm); /* write out the face and its distribution */ if (average(il, distarr, n)) { if (il->sampdens > 0) @@ -282,7 +374,6 @@ my_face( /* make an illum face */ /* clean up */ freeface(ob); return(0); -#undef MAXMISS } @@ -310,8 +401,6 @@ my_sphere( /* make an illum sphere */ nalt = sqrt(2./PI*n) + .5; nazi = PI/2.*nalt + .5; } - if (il->sd != NULL) - objerror(ob, WARNING, "BSDF ignored"); n = nalt*nazi; newdist(n); dim[0] = random(); @@ -340,6 +429,7 @@ my_sphere( /* make an illum sphere */ /* send sample */ raysamp(dim[1]*nazi+dim[2], org, dir); } + /* wait for all rays to finish */ rayclean(); /* write out the sphere and its distribution */ if (average(il, distarr, n)) { @@ -370,26 +460,18 @@ my_ring( /* make an illum ring */ FVECT u, v; MAT4 xfm; CONE *co; - register int i, j; + int i, j; /* get/check arguments */ co = getcone(ob, 0); /* set up sampling */ - if (il->sampdens <= 0) - nalt = nazi = 1; - else { + if (il->sampdens <= 0) { + nalt = nazi = 1; /* diffuse assumption */ + } else { n = PI * il->sampdens; nalt = sqrt(n/PI) + .5; nazi = PI*nalt + .5; } - if (il->sd != NULL) { - if (!getBSDF_xfm(xfm, co->ad, il->udir)) { - objerror(ob, WARNING, "illegal up direction"); - freecone(ob); - return(my_default(ob, il, nm)); - } - n = il->sd->ninc; - } else - n = nazi*nalt; + n = nazi*nalt; newdist(n); mkaxes(u, v, co->ad); dim[0] = random(); @@ -398,38 +480,57 @@ my_ring( /* make an illum ring */ for (i = 0; i < il->nsamps; i++) { /* next sample point */ h = ilhash(dim,2) + i; - /* random direction */ - if (il->sd != NULL) { - r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm); - } else { - multisamp(sp, 2, urand(h)); - alti = dim[1]/nazi; - r1 = (alti + sp[0])/nalt; - r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi; - flatdir(dn, r1, r2); - for (j = 0; j < 3; j++) + /* randomize direction */ + multisamp(sp, 2, urand(h)); + alti = dim[1]/nazi; + r1 = (alti + sp[0])/nalt; + r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi; + flatdir(dn, r1, r2); + for (j = 0; j < 3; j++) dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j]; - } - /* random location */ + /* randomize location */ multisamp(sp, 2, urand(h+8371)); r3 = sqrt(CO_R0(co)*CO_R0(co) + sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co))); r2 = 2.*PI*sp[1]; r1 = r3*cos(r2); r2 = r3*sin(r2); - if (il->sd != NULL && DOT(dir, co->ad) < -FTINY) - r3 = -1.0001*il->thick - .0001; - else - r3 = .0001; + r3 = 5.*FTINY; for (j = 0; j < 3; j++) org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] + r3*co->ad[j]; /* send sample */ raysamp(dim[1], org, dir); } + /* add in direct component? */ + if (il->flags & IL_LIGHT) { + MAT4 ixfm; + for (i = 3; i--; ) { + ixfm[i][0] = u[i]; + ixfm[i][1] = v[i]; + ixfm[i][2] = co->ad[i]; + ixfm[i][3] = 0.; + } + ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.; + ixfm[3][3] = 1.; + dim[0] = random(); + for (i = 0; i < il->nsamps; i++) { + /* randomize location */ + h = dim[0] + samplendx++; + multisamp(sp, 2, urand(h)); + r3 = sqrt(CO_R0(co)*CO_R0(co) + + sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co))); + r2 = 2.*PI*sp[1]; + r1 = r3*cos(r2); + r2 = r3*sin(r2); + for (j = 0; j < 3; j++) + org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j]; + /* sample source rays */ + srcsamps(il, org, co->ad, ixfm); + } + } + /* wait for all rays to finish */ rayclean(); - if (il->sd != NULL) /* run distribution through BSDF */ - redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm); /* write out the ring and its distribution */ if (average(il, distarr, n)) { if (il->sampdens > 0)