--- ray/src/gen/mkillum2.c 2007/09/18 03:35:09 2.19 +++ ray/src/gen/mkillum2.c 2007/11/05 23:40:26 2.23 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: mkillum2.c,v 2.19 2007/09/18 03:35:09 greg Exp $"; +static const char RCSid[] = "$Id: mkillum2.c,v 2.23 2007/11/05 23:40:26 greg Exp $"; #endif /* * Routines to do the actual calculation for mkillum @@ -13,16 +13,11 @@ static const char RCSid[] = "$Id: mkillum2.c,v 2.19 20 #include "random.h" -static void mkaxes(FVECT u, FVECT v, FVECT n); -static void rounddir(FVECT dv, double alt, double azi); -static void flatdir(FVECT dv, double alt, double azi); +COLORV * distarr = NULL; /* distribution array */ +int distsiz = 0; -static COLORV * distarr = NULL; /* distribution array */ -static int distsiz = 0; - - -static void +void newdist( /* allocate & clear distribution array */ int siz ) @@ -35,17 +30,18 @@ newdist( /* allocate & clear distribution array */ return; } if (distsiz < siz) { - free((void *)distarr); - distarr = (COLORV *)malloc(sizeof(COLORV)*3*siz); + if (distsiz > 0) + free((void *)distarr); + distarr = (COLORV *)malloc(sizeof(COLOR)*siz); if (distarr == NULL) - error(SYSTEM, "Out of memory in newdist"); + error(SYSTEM, "out of memory in newdist"); distsiz = siz; } - memset(distarr, '\0', sizeof(COLORV)*3*siz); + memset(distarr, '\0', sizeof(COLOR)*siz); } -static int +int process_ray(RAY *r, int rv) { COLORV *colp; @@ -53,16 +49,16 @@ process_ray(RAY *r, int rv) if (rv == 0) return(0); if (rv < 0) - error(USER, "Ray tracing process died"); + error(USER, "ray tracing process died"); if (r->rno >= distsiz) - error(INTERNAL, "Bad returned index in process_ray"); + error(INTERNAL, "bad returned index in process_ray"); colp = &distarr[r->rno * 3]; addcolor(colp, r->rcol); return(1); } -static void +void raysamp( /* queue a ray sample */ int ndx, FVECT org, @@ -73,7 +69,7 @@ raysamp( /* queue a ray sample */ int rv; if ((ndx < 0) | (ndx >= distsiz)) - error(INTERNAL, "Bad index in raysamp"); + error(INTERNAL, "bad index in raysamp"); VCOPY(myRay.rorg, org); VCOPY(myRay.rdir, dir); myRay.rmax = .0; @@ -84,7 +80,7 @@ raysamp( /* queue a ray sample */ } -static void +void rayclean() /* finish all pending rays */ { RAY myRay; @@ -94,6 +90,60 @@ rayclean() /* finish all pending rays */ } +static void +mkaxes( /* compute u and v to go with n */ + FVECT u, + FVECT v, + FVECT n +) +{ + register int i; + + v[0] = v[1] = v[2] = 0.0; + for (i = 0; i < 3; i++) + if (n[i] < 0.6 && n[i] > -0.6) + break; + v[i] = 1.0; + fcross(u, v, n); + normalize(u); + fcross(v, n, u); +} + + +static void +rounddir( /* compute uniform spherical direction */ + register FVECT dv, + double alt, + double azi +) +{ + double d1, d2; + + dv[2] = 1. - 2.*alt; + d1 = sqrt(1. - dv[2]*dv[2]); + d2 = 2.*PI * azi; + dv[0] = d1*cos(d2); + dv[1] = d1*sin(d2); +} + + +static void +flatdir( /* compute uniform hemispherical direction */ + register FVECT dv, + double alt, + double azi +) +{ + double d1, d2; + + d1 = sqrt(alt); + d2 = 2.*PI * azi; + dv[0] = d1*cos(d2); + dv[1] = d1*sin(d2); + dv[2] = sqrt(1. - alt); +} + + int my_default( /* default illum action */ OBJREC *ob, @@ -117,14 +167,15 @@ my_face( /* make an illum face */ ) { #define MAXMISS (5*n*il->nsamps) - int dim[3]; - int n, nalt, nazi, h; + int dim[2]; + int n, nalt, nazi, h, alti; double sp[2], r1, r2; FVECT dn, org, dir; FVECT u, v; double ur[2], vr[2]; + MAT4 xfm; int nmisses; - register FACE *fa; + FACE *fa; register int i, j; /* get/check arguments */ fa = getface(ob); @@ -133,16 +184,25 @@ my_face( /* make an illum face */ return(my_default(ob, il, nm)); } /* set up sampling */ - if (il->sampdens <= 0) - nalt = nazi = 1; - else { - n = PI * il->sampdens; - 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 { + if (il->sampdens <= 0) { + nalt = nazi = 1; /* diffuse assumption */ + } else { + n = PI * il->sampdens; + nalt = sqrt(n/PI) + .5; + nazi = PI*nalt + .5; + } + n = nazi*nalt; } - n = nalt*nazi; newdist(n); - /* take first edge longer than sqrt(area) */ + /* take first edge >= sqrt(area) */ for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) { u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0]; u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1]; @@ -170,18 +230,23 @@ my_face( /* make an illum face */ dim[0] = random(); /* sample polygon */ nmisses = 0; - for (dim[1] = 0; dim[1] < nalt; dim[1]++) - for (dim[2] = 0; dim[2] < nazi; dim[2]++) + for (dim[1] = 0; dim[1] < n; dim[1]++) for (i = 0; i < il->nsamps; i++) { - /* random direction */ - h = ilhash(dim, 3) + i; - multisamp(sp, 2, urand(h)); - r1 = (dim[1] + sp[0])/nalt; - r2 = (dim[2] + 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 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++) + dir[j] = -dn[0]*u[j] - dn[1]*v[j] - + dn[2]*fa->norm[j]; + } + /* randomize location */ do { multisamp(sp, 2, urand(h+4862+nmisses)); r1 = ur[0] + (ur[1]-ur[0]) * sp[0]; @@ -194,17 +259,25 @@ my_face( /* make an illum face */ objerror(ob, WARNING, "bad aspect"); rayclean(); freeface(ob); - free((void *)distarr); return(my_default(ob, il, nm)); } + if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY) + r1 = -1.0001*il->thick - .0001; + else + r1 = .0001; for (j = 0; j < 3; j++) - org[j] += .001*fa->norm[j]; + org[j] += r1*fa->norm[j]; /* send sample */ - raysamp(dim[1]*nazi+dim[2], org, dir); + raysamp(dim[1], org, dir); } rayclean(); + if (il->sd != NULL) { /* run distribution through BSDF */ + nalt = sqrt(il->sd->nout/PI) + .5; + nazi = PI*nalt + .5; + redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm); + } /* write out the face and its distribution */ - if (average(il, distarr, nalt*nazi)) { + if (average(il, distarr, n)) { if (il->sampdens > 0) flatout(il, distarr, nalt, nazi, u, v, fa->norm); illumout(il, ob); @@ -241,6 +314,8 @@ 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(); @@ -271,7 +346,7 @@ my_sphere( /* make an illum sphere */ } rayclean(); /* write out the sphere and its distribution */ - if (average(il, distarr, nalt*nazi)) { + if (average(il, distarr, n)) { if (il->sampdens > 0) roundout(il, distarr, nalt, nazi); else @@ -291,55 +366,80 @@ my_ring( /* make an illum ring */ char *nm ) { - int dim[3]; - int n, nalt, nazi; - double sp[4], r1, r2, r3; + int dim[2]; + int n, nalt, nazi, alti; + double sp[2], r1, r2, r3; + int h; FVECT dn, org, dir; FVECT u, v; - register CONE *co; + MAT4 xfm; + CONE *co; register int i, j; /* get/check arguments */ co = getcone(ob, 0); /* set up sampling */ - if (il->sampdens <= 0) - nalt = nazi = 1; - 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 { + if (il->sampdens <= 0) { + nalt = nazi = 1; /* diffuse assumption */ + } else { + n = PI * il->sampdens; + nalt = sqrt(n/PI) + .5; + nazi = PI*nalt + .5; + } + n = nazi*nalt; } - n = nalt*nazi; newdist(n); mkaxes(u, v, co->ad); dim[0] = random(); /* sample disk */ - for (dim[1] = 0; dim[1] < nalt; dim[1]++) - for (dim[2] = 0; dim[2] < nazi; dim[2]++) + for (dim[1] = 0; dim[1] < n; dim[1]++) for (i = 0; i < il->nsamps; i++) { /* next sample point */ - multisamp(sp, 4, urand(ilhash(dim,3)+i)); - /* random direction */ - r1 = (dim[1] + sp[0])/nalt; - r2 = (dim[2] + sp[1] - .5)/nazi; - flatdir(dn, r1, r2); - for (j = 0; j < 3; j++) + h = ilhash(dim,2) + i; + /* randomize 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++) 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[2]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co))); - r2 = 2.*PI*sp[3]; + 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; for (j = 0; j < 3; j++) org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] + - .001*co->ad[j]; - + r3*co->ad[j]; /* send sample */ - raysamp(dim[1]*nazi+dim[2], org, dir); + raysamp(dim[1], org, dir); } rayclean(); + if (il->sd != NULL) { /* run distribution through BSDF */ + nalt = sqrt(il->sd->nout/PI) + .5; + nazi = PI*nalt + .5; + redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm); + } /* write out the ring and its distribution */ - if (average(il, distarr, nalt*nazi)) { + if (average(il, distarr, n)) { if (il->sampdens > 0) flatout(il, distarr, nalt, nazi, u, v, co->ad); illumout(il, ob); @@ -348,58 +448,4 @@ my_ring( /* make an illum ring */ /* clean up */ freecone(ob); return(1); -} - - -static void -mkaxes( /* compute u and v to go with n */ - FVECT u, - FVECT v, - FVECT n -) -{ - register int i; - - v[0] = v[1] = v[2] = 0.0; - for (i = 0; i < 3; i++) - if (n[i] < 0.6 && n[i] > -0.6) - break; - v[i] = 1.0; - fcross(u, v, n); - normalize(u); - fcross(v, n, u); -} - - -static void -rounddir( /* compute uniform spherical direction */ - register FVECT dv, - double alt, - double azi -) -{ - double d1, d2; - - dv[2] = 1. - 2.*alt; - d1 = sqrt(1. - dv[2]*dv[2]); - d2 = 2.*PI * azi; - dv[0] = d1*cos(d2); - dv[1] = d1*sin(d2); -} - - -static void -flatdir( /* compute uniform hemispherical direction */ - register FVECT dv, - double alt, - double azi -) -{ - double d1, d2; - - d1 = sqrt(alt); - d2 = 2.*PI * azi; - dv[0] = d1*cos(d2); - dv[1] = d1*sin(d2); - dv[2] = sqrt(1. - alt); }