--- ray/src/rt/ambcomp.c 2014/04/23 17:30:10 2.30 +++ ray/src/rt/ambcomp.c 2014/05/04 01:02:13 2.48 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: ambcomp.c,v 2.30 2014/04/23 17:30:10 greg Exp $"; +static const char RCSid[] = "$Id: ambcomp.c,v 2.48 2014/05/04 01:02:13 greg Exp $"; #endif /* * Routines to compute "ambient" values using Monte Carlo @@ -8,6 +8,10 @@ static const char RCSid[] = "$Id: ambcomp.c,v 2.30 201 * for Irradiance Caching" by Schwarzhaupt, Wann Jensen, & Jarosz * from ACM SIGGRAPH Asia 2012 conference proceedings. * + * Added book-keeping optimization to avoid calculations that would + * cancel due to traversal both directions on edges that are adjacent + * to same-valued triangles. This cuts about half of Hessian math. + * * Declarations of external symbols in ambient.h */ @@ -19,27 +23,95 @@ static const char RCSid[] = "$Id: ambcomp.c,v 2.30 201 #ifdef NEWAMB +/* #define AHEM_MARG 1.2 /* hem margin */ + extern void SDsquare2disk(double ds[2], double seedx, double seedy); + /* vertex direction bit positions */ +#define VDB_xy 0 +#define VDB_y 01 +#define VDB_x 02 +#define VDB_Xy 03 +#define VDB_xY 04 +#define VDB_X 05 +#define VDB_Y 06 +#define VDB_XY 07 + /* get opposite vertex direction bit */ +#define VDB_OPP(f) (~(f) & 07) + /* adjacent triangle vertex flags */ +static const int adjacent_trifl[8] = { + 0, /* forbidden diagonal */ + 1<sa[(i)*(h)->ns + (j)] +#define ambndx(h,i,j) ((i)*(h)->ns + (j)) +#define ambsam(h,i,j) (h)->sa[ambndx(h,i,j)] typedef struct { - FVECT r_i, r_i1, e_i, rI2_eJ2; - double nf, I1, I2; + FVECT r_i, r_i1, e_i, rcp, rI2_eJ2; + double I1, I2; + int valid; } FFTRI; /* vectors and coefficients for Hessian calculation */ +/* Get index for adjacent vertex */ +static int +adjacent_verti(AMBHEMI *hp, int i, int j, int dbit) +{ + int i0 = i*hp->ns + j; + + switch (dbit) { + case VDB_y: return(i0 - hp->ns); + case VDB_x: return(i0 - 1); + case VDB_Xy: return(i0 - hp->ns + 1); + case VDB_xY: return(i0 + hp->ns - 1); + case VDB_X: return(i0 + 1); + case VDB_Y: return(i0 + hp->ns); + /* the following should never occur */ + case VDB_xy: return(i0 - hp->ns - 1); + case VDB_XY: return(i0 + hp->ns + 1); + } + return(-1); +} + + +/* Get vertex direction bit for the opposite edge to complete triangle */ +static int +vdb_edge(int db1, int db2) +{ + switch (db1) { + case VDB_x: return(db2==VDB_y ? VDB_Xy : VDB_Y); + case VDB_y: return(db2==VDB_x ? VDB_xY : VDB_X); + case VDB_X: return(db2==VDB_Xy ? VDB_y : VDB_xY); + case VDB_Y: return(db2==VDB_xY ? VDB_x : VDB_Xy); + case VDB_xY: return(db2==VDB_x ? VDB_y : VDB_X); + case VDB_Xy: return(db2==VDB_y ? VDB_x : VDB_Y); + } + error(INTERNAL, "forbidden diagonal in vdb_edge()"); + return(-1); +} + + static AMBHEMI * inithemi( /* initialize sampling hemisphere */ COLOR ac, @@ -59,8 +131,7 @@ inithemi( /* initialize sampling hemisphere */ if (n < i) n = i; /* allocate sampling array */ - hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + - sizeof(struct s_ambsamp)*(n*n - 1)); + hp = (AMBHEMI *)malloc(sizeof(AMBHEMI) + sizeof(AMBSAMP)*(n*n - 1)); if (hp == NULL) return(NULL); hp->rp = r; @@ -70,14 +141,14 @@ inithemi( /* initialize sampling hemisphere */ d = 1.0/(n*n); scalecolor(hp->acoef, d); /* make tangent plane axes */ - hp->uy[0] = 0.1 - 0.2*frandom(); - hp->uy[1] = 0.1 - 0.2*frandom(); - hp->uy[2] = 0.1 - 0.2*frandom(); - for (i = 0; i < 3; i++) - if (r->ron[i] < 0.6 && r->ron[i] > -0.6) + hp->uy[0] = 0.5 - frandom(); + hp->uy[1] = 0.5 - frandom(); + hp->uy[2] = 0.5 - frandom(); + for (i = 3; i--; ) + if ((-0.6 < r->ron[i]) & (r->ron[i] < 0.6)) break; - if (i >= 3) - error(CONSISTENCY, "bad ray direction in inithemi()"); + if (i < 0) + error(CONSISTENCY, "bad ray direction in inithemi"); hp->uy[i] = 1.0; VCROSS(hp->ux, hp->uy, r->ron); normalize(hp->ux); @@ -87,77 +158,261 @@ inithemi( /* initialize sampling hemisphere */ } -static struct s_ambsamp * -ambsample( /* sample an ambient direction */ - AMBHEMI *hp, - int i, - int j -) +/* Sample ambient division and apply weighting coefficient */ +static int +getambsamp(RAY *arp, AMBHEMI *hp, int i, int j, int n) { - struct s_ambsamp *ap = &ambsamp(hp,i,j); - RAY ar; - double spt[2], zd; - int ii; + int hlist[3], ii; + double spt[2], zd; /* ambient coefficient for weight */ if (ambacc > FTINY) - setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); + setcolor(arp->rcoef, AVGREFL, AVGREFL, AVGREFL); else - copycolor(ar.rcoef, hp->acoef); - if (rayorigin(&ar, AMBIENT, hp->rp, ar.rcoef) < 0) { - setcolor(ap->v, 0., 0., 0.); - VCOPY(ap->p, hp->rp->rop); - return(NULL); /* no sample taken */ - } + copycolor(arp->rcoef, hp->acoef); + if (rayorigin(arp, AMBIENT, hp->rp, arp->rcoef) < 0) + return(0); if (ambacc > FTINY) { - multcolor(ar.rcoef, hp->acoef); - scalecolor(ar.rcoef, 1./AVGREFL); + multcolor(arp->rcoef, hp->acoef); + scalecolor(arp->rcoef, 1./AVGREFL); } - /* generate hemispherical sample */ - SDsquare2disk(spt, (i+.1+.8*frandom())/hp->ns, - (j+.1+.8*frandom())/hp->ns ); + hlist[0] = hp->rp->rno; + hlist[1] = j; + hlist[2] = i; + multisamp(spt, 2, urand(ilhash(hlist,3)+n)); + if (!n) { /* avoid border samples for n==0 */ + if ((spt[0] < 0.1) | (spt[0] >= 0.9)) + spt[0] = 0.1 + 0.8*frandom(); + if ((spt[1] < 0.1) | (spt[1] >= 0.9)) + spt[1] = 0.1 + 0.8*frandom(); + } + SDsquare2disk(spt, (j+spt[1])/hp->ns, (i+spt[0])/hp->ns); zd = sqrt(1. - spt[0]*spt[0] - spt[1]*spt[1]); for (ii = 3; ii--; ) - ar.rdir[ii] = spt[0]*hp->ux[ii] + + arp->rdir[ii] = spt[0]*hp->ux[ii] + spt[1]*hp->uy[ii] + zd*hp->rp->ron[ii]; - checknorm(ar.rdir); - dimlist[ndims++] = i*hp->ns + j + 90171; - rayvalue(&ar); /* evaluate ray */ - ndims--; - multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ + checknorm(arp->rdir); + dimlist[ndims++] = ambndx(hp,i,j) + 90171; + rayvalue(arp); /* evaluate ray */ + ndims--; /* apply coefficient */ + multcolor(arp->rcol, arp->rcoef); + return(1); +} + + +static AMBSAMP * +ambsample( /* initial ambient division sample */ + AMBHEMI *hp, + int i, + int j +) +{ + AMBSAMP *ap = &ambsam(hp,i,j); + RAY ar; + /* generate hemispherical sample */ + if (!getambsamp(&ar, hp, i, j, 0) || ar.rt <= FTINY) { + memset(ap, 0, sizeof(AMBSAMP)); + return(NULL); + } + ap->d = 1.0/ar.rt; /* limit vertex distance */ + if (ar.rt > 10.0*thescene.cusize) + ar.rt = 10.0*thescene.cusize; + VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); copycolor(ap->v, ar.rcol); - if (ar.rt > 20.0*maxarad) /* limit vertex distance */ - VSUM(ap->p, ar.rorg, ar.rdir, 20.0*maxarad); - else - VCOPY(ap->p, ar.rop); return(ap); } +/* Estimate errors based on ambient division differences */ +static float * +getambdiffs(AMBHEMI *hp) +{ + float *earr = (float *)calloc(hp->ns*hp->ns, sizeof(float)); + float *ep; + AMBSAMP *ap; + double b, d2; + int i, j; + + if (earr == NULL) /* out of memory? */ + return(NULL); + /* compute squared neighbor diffs */ + for (ap = hp->sa, ep = earr, i = 0; i < hp->ns; i++) + for (j = 0; j < hp->ns; j++, ap++, ep++) { + b = bright(ap[0].v); + if (i) { /* from above */ + d2 = b - bright(ap[-hp->ns].v); + d2 *= d2; + ep[0] += d2; + ep[-hp->ns] += d2; + } + if (j) { /* from behind */ + d2 = b - bright(ap[-1].v); + d2 *= d2; + ep[0] += d2; + ep[-1] += d2; + } + } + /* correct for number of neighbors */ + earr[0] *= 2.f; + earr[hp->ns-1] *= 2.f; + earr[(hp->ns-1)*hp->ns] *= 2.f; + earr[(hp->ns-1)*hp->ns + hp->ns-1] *= 2.f; + for (i = 1; i < hp->ns-1; i++) { + earr[i*hp->ns] *= 4./3.; + earr[i*hp->ns + hp->ns-1] *= 4./3.; + } + for (j = 1; j < hp->ns-1; j++) { + earr[j] *= 4./3.; + earr[(hp->ns-1)*hp->ns + j] *= 4./3.; + } + return(earr); +} + + +/* Perform super-sampling on hemisphere (introduces bias) */ +static void +ambsupersamp(double acol[3], AMBHEMI *hp, int cnt) +{ + float *earr = getambdiffs(hp); + double e2sum = 0; + AMBSAMP *ap; + RAY ar; + COLOR asum; + float *ep; + int i, j, n; + + if (earr == NULL) /* just skip calc. if no memory */ + return; + /* add up estimated variances */ + for (ep = earr + hp->ns*hp->ns; ep-- > earr; ) + e2sum += *ep; + ep = earr; /* perform super-sampling */ + for (ap = hp->sa, i = 0; i < hp->ns; i++) + for (j = 0; j < hp->ns; j++, ap++) { + int nss = *ep/e2sum*cnt + frandom(); + setcolor(asum, 0., 0., 0.); + for (n = 1; n <= nss; n++) { + if (!getambsamp(&ar, hp, i, j, n)) { + nss = n-1; + break; + } + addcolor(asum, ar.rcol); + } + if (nss) { /* update returned ambient value */ + const double ssf = 1./(nss + 1); + for (n = 3; n--; ) + acol[n] += ssf*colval(asum,n) + + (ssf - 1.)*colval(ap->v,n); + } + e2sum -= *ep++; /* update remainders */ + cnt -= nss; + } + free(earr); +} + + +/* Compute vertex flags, indicating farthest in each direction */ +static uby8 * +vertex_flags(AMBHEMI *hp) +{ + uby8 *vflags = (uby8 *)calloc(hp->ns*hp->ns, sizeof(uby8)); + uby8 *vf; + AMBSAMP *ap; + int i, j; + + if (vflags == NULL) + error(SYSTEM, "out of memory in vertex_flags()"); + vf = vflags; + ap = hp->sa; /* compute farthest along first row */ + for (j = 0; j < hp->ns-1; j++, vf++, ap++) + if (ap[0].d <= ap[1].d) + vf[0] |= 1<ns; i++) { + for (j = 0; j < hp->ns-1; j++, vf++, ap++) { + if (ap[0].d <= ap[-hp->ns].d) /* row before */ + vf[0] |= 1<ns] |= 1<ns].d) /* diagonal we care about */ + vf[0] |= 1<ns] |= 1<ns].d) /* final column edge */ + vf[0] |= 1<ns] |= 1<sa[v0].v,CIEY)); + v1 = adjacent_verti(hp, i, j, dbit1); + if (vflags[v0] & 1<v2 */ + return(colval(hp->sa[v1].v,CIEY)); + v2 = adjacent_verti(hp, i, j, dbit2); + if (vflags[v0] & 1<v1 */ + return(colval(hp->sa[v2].v,CIEY)); + /* else check if v1>v2 */ + if (vflags[v1] & 1<sa[v1].v,CIEY)); + return(colval(hp->sa[v2].v,CIEY)); +} + + /* Compute vectors and coefficients for Hessian/gradient calcs */ static void -comp_fftri(FFTRI *ftp, float ap0[3], float ap1[3], FVECT rop) +comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags) { - FVECT vcp; - double dot_e, dot_er, dot_r, dot_r1, J2; - int i; + const int i0 = ambndx(hp,i,j); + double rdot_cp, dot_e, dot_er, rdot_r, rdot_r1, J2; + int i1, ii; - VSUB(ftp->r_i, ap0, rop); - VSUB(ftp->r_i1, ap1, rop); - VSUB(ftp->e_i, ap1, ap0); - VCROSS(vcp, ftp->e_i, ftp->r_i); - ftp->nf = 1.0/DOT(vcp,vcp); + ftp->valid = 0; /* check if we can skip this edge */ + ii = adjacent_trifl[dbit]; + if ((vflags[i0] & ii) == ii) /* cancels if vertex used as value */ + return; + i1 = adjacent_verti(hp, i, j, dbit); + ii = adjacent_trifl[VDB_OPP(dbit)]; + if ((vflags[i1] & ii) == ii) /* on either end (for both triangles) */ + return; + /* else go ahead with calculation */ + VSUB(ftp->r_i, hp->sa[i0].p, hp->rp->rop); + VSUB(ftp->r_i1, hp->sa[i1].p, hp->rp->rop); + VSUB(ftp->e_i, hp->sa[i1].p, hp->sa[i0].p); + VCROSS(ftp->rcp, ftp->r_i, ftp->r_i1); + rdot_cp = 1.0/DOT(ftp->rcp,ftp->rcp); dot_e = DOT(ftp->e_i,ftp->e_i); dot_er = DOT(ftp->e_i, ftp->r_i); - dot_r = DOT(ftp->r_i,ftp->r_i); - dot_r1 = DOT(ftp->r_i1,ftp->r_i1); - ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) / sqrt(dot_r*dot_r1) ) * - sqrt( ftp->nf ); - ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)/dot_r1 - dot_er/dot_r + - dot_e*ftp->I1 )*0.5*ftp->nf; - J2 = 0.5/dot_e*( 1.0/dot_r - 1.0/dot_r1 ) - dot_er/dot_e*ftp->I2; - for (i = 3; i--; ) - ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; + rdot_r = 1.0/DOT(ftp->r_i,ftp->r_i); + rdot_r1 = 1.0/DOT(ftp->r_i1,ftp->r_i1); + ftp->I1 = acos( DOT(ftp->r_i, ftp->r_i1) * sqrt(rdot_r*rdot_r1) ) * + sqrt( rdot_cp ); + ftp->I2 = ( DOT(ftp->e_i, ftp->r_i1)*rdot_r1 - dot_er*rdot_r + + dot_e*ftp->I1 )*0.5*rdot_cp; + J2 = ( 0.5*(rdot_r - rdot_r1) - dot_er*ftp->I2 ) / dot_e; + for (ii = 3; ii--; ) + ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii]; + ftp->valid++; } @@ -178,28 +433,32 @@ compose_matrix(FVECT mat[3], FVECT va, FVECT vb) static void comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) { - FVECT vcp; + FVECT ncp; FVECT m1[3], m2[3], m3[3], m4[3]; double d1, d2, d3, d4; double I3, J3, K3; int i, j; + + if (!ftp->valid) { /* preemptive test */ + memset(hess, 0, sizeof(FVECT)*3); + return; + } /* compute intermediate coefficients */ d1 = 1.0/DOT(ftp->r_i,ftp->r_i); d2 = 1.0/DOT(ftp->r_i1,ftp->r_i1); d3 = 1.0/DOT(ftp->e_i,ftp->e_i); d4 = DOT(ftp->e_i, ftp->r_i); - I3 = 0.25*ftp->nf*( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + - 3.0/d3*ftp->I2 ); + I3 = ( DOT(ftp->e_i, ftp->r_i1)*d2*d2 - d4*d1*d1 + 3.0/d3*ftp->I2 ) + / ( 4.0*DOT(ftp->rcp,ftp->rcp) ); J3 = 0.25*d3*(d1*d1 - d2*d2) - d4*d3*I3; K3 = d3*(ftp->I2 - I3/d1 - 2.0*d4*J3); /* intermediate matrices */ - VCROSS(vcp, nrm, ftp->e_i); - compose_matrix(m1, vcp, ftp->rI2_eJ2); + VCROSS(ncp, nrm, ftp->e_i); + compose_matrix(m1, ncp, ftp->rI2_eJ2); compose_matrix(m2, ftp->r_i, ftp->r_i); compose_matrix(m3, ftp->e_i, ftp->e_i); compose_matrix(m4, ftp->r_i, ftp->e_i); - VCROSS(vcp, ftp->r_i, ftp->e_i); - d1 = DOT(nrm, vcp); + d1 = DOT(nrm, ftp->rcp); d2 = -d1*ftp->I2; d1 *= 2.0; for (i = 3; i--; ) /* final matrix sum */ @@ -229,7 +488,7 @@ rev_hessian(FVECT hess[3]) /* Add to radiometric Hessian from the given triangle */ static void add2hessian(FVECT hess[3], FVECT ehess1[3], - FVECT ehess2[3], FVECT ehess3[3], COLORV v) + FVECT ehess2[3], FVECT ehess3[3], double v) { int i, j; @@ -243,15 +502,18 @@ add2hessian(FVECT hess[3], FVECT ehess1[3], static void comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm) { - FVECT vcp; + FVECT ncp; double f1; int i; - VCROSS(vcp, ftp->r_i, ftp->r_i1); - f1 = 2.0*DOT(nrm, vcp); - VCROSS(vcp, nrm, ftp->e_i); + if (!ftp->valid) { /* preemptive test */ + memset(grad, 0, sizeof(FVECT)); + return; + } + f1 = 2.0*DOT(nrm, ftp->rcp); + VCROSS(ncp, nrm, ftp->e_i); for (i = 3; i--; ) - grad[i] = (0.5/PI)*( ftp->I1*vcp[i] + f1*ftp->rI2_eJ2[i] ); + grad[i] = (0.5/PI)*( ftp->I1*ncp[i] + f1*ftp->rI2_eJ2[i] ); } @@ -267,7 +529,7 @@ rev_gradient(FVECT grad) /* Add to displacement gradient from the given triangle */ static void -add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, COLORV v) +add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, FVECT egrad3, double v) { int i; @@ -276,32 +538,6 @@ add2gradient(FVECT grad, FVECT egrad1, FVECT egrad2, F } -/* Return brightness of furthest ambient sample */ -static COLORV -back_ambval(struct s_ambsamp *ap1, struct s_ambsamp *ap2, - struct s_ambsamp *ap3, FVECT orig) -{ - COLORV vback; - FVECT vec; - double d2, d2best; - - VSUB(vec, ap1->p, orig); - d2best = DOT(vec,vec); - vback = colval(ap1->v,CIEY); - VSUB(vec, ap2->p, orig); - d2 = DOT(vec,vec); - if (d2 > d2best) { - d2best = d2; - vback = colval(ap2->v,CIEY); - } - VSUB(vec, ap3->p, orig); - d2 = DOT(vec,vec); - if (d2 > d2best) - return(colval(ap3->v,CIEY)); - return(vback); -} - - /* Compute anisotropic radii and eigenvector directions */ static int eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3]) @@ -319,11 +555,13 @@ eigenvectors(FVECT uv[2], float ra[2], FVECT hessian[3 hess2[0][1] = DOT(uv[0], b); hess2[1][0] = DOT(uv[1], a); hess2[1][1] = DOT(uv[1], b); - /* compute eigenvalues */ - if ( quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], - hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]) != 2 || - (evalue[0] = fabs(evalue[0])) <= FTINY*FTINY || - (evalue[1] = fabs(evalue[1])) <= FTINY*FTINY ) + /* compute eigenvalue(s) */ + i = quadratic(evalue, 1.0, -hess2[0][0]-hess2[1][1], + hess2[0][0]*hess2[1][1]-hess2[0][1]*hess2[1][0]); + if (i == 1) /* double-root (circle) */ + evalue[1] = evalue[0]; + if (!i || ((evalue[0] = fabs(evalue[0])) <= FTINY*FTINY) | + ((evalue[1] = fabs(evalue[1])) <= FTINY*FTINY) ) error(INTERNAL, "bad eigenvalue calculation"); if (evalue[0] > evalue[1]) { @@ -360,6 +598,7 @@ ambHessian( /* anisotropic radii & pos. gradient */ static char memerrmsg[] = "out of memory in ambHessian()"; FVECT (*hessrow)[3] = NULL; FVECT *gradrow = NULL; + uby8 *vflags; FVECT hessian[3]; FVECT gradient; FFTRI fftr; @@ -381,10 +620,11 @@ ambHessian( /* anisotropic radii & pos. gradient */ error(SYSTEM, memerrmsg); memset(gradient, 0, sizeof(gradient)); } + /* get vertex position flags */ + vflags = vertex_flags(hp); /* compute first row of edges */ for (j = 0; j < hp->ns-1; j++) { - comp_fftri(&fftr, ambsamp(hp,0,j).p, - ambsamp(hp,0,j+1).p, hp->rp->rop); + comp_fftri(&fftr, hp, 0, j, VDB_X, vflags); if (hessrow != NULL) comp_hessian(hessrow[j], &fftr, hp->rp->ron); if (gradrow != NULL) @@ -394,8 +634,7 @@ ambHessian( /* anisotropic radii & pos. gradient */ for (i = 0; i < hp->ns-1; i++) { FVECT hesscol[3]; /* compute first vertical edge */ FVECT gradcol; - comp_fftri(&fftr, ambsamp(hp,i,0).p, - ambsamp(hp,i+1,0).p, hp->rp->rop); + comp_fftri(&fftr, hp, i, 0, VDB_Y, vflags); if (hessrow != NULL) comp_hessian(hesscol, &fftr, hp->rp->ron); if (gradrow != NULL) @@ -403,34 +642,29 @@ ambHessian( /* anisotropic radii & pos. gradient */ for (j = 0; j < hp->ns-1; j++) { FVECT hessdia[3]; /* compute triangle contributions */ FVECT graddia; - COLORV backg; - backg = back_ambval(&ambsamp(hp,i,j), &ambsamp(hp,i,j+1), - &ambsamp(hp,i+1,j), hp->rp->rop); + double backg; + backg = back_ambval(hp, i, j, VDB_X, VDB_Y, vflags); /* diagonal (inner) edge */ - comp_fftri(&fftr, ambsamp(hp,i,j+1).p, - ambsamp(hp,i+1,j).p, hp->rp->rop); + comp_fftri(&fftr, hp, i, j+1, VDB_xY, vflags); if (hessrow != NULL) { comp_hessian(hessdia, &fftr, hp->rp->ron); rev_hessian(hesscol); add2hessian(hessian, hessrow[j], hessdia, hesscol, backg); } - if (gradient != NULL) { + if (gradrow != NULL) { comp_gradient(graddia, &fftr, hp->rp->ron); rev_gradient(gradcol); add2gradient(gradient, gradrow[j], graddia, gradcol, backg); } /* initialize edge in next row */ - comp_fftri(&fftr, ambsamp(hp,i+1,j+1).p, - ambsamp(hp,i+1,j).p, hp->rp->rop); + comp_fftri(&fftr, hp, i+1, j+1, VDB_x, vflags); if (hessrow != NULL) comp_hessian(hessrow[j], &fftr, hp->rp->ron); if (gradrow != NULL) comp_gradient(gradrow[j], &fftr, hp->rp->ron); /* new column edge & paired triangle */ - backg = back_ambval(&ambsamp(hp,i,j+1), &ambsamp(hp,i+1,j+1), - &ambsamp(hp,i+1,j), hp->rp->rop); - comp_fftri(&fftr, ambsamp(hp,i,j+1).p, ambsamp(hp,i+1,j+1).p, - hp->rp->rop); + backg = back_ambval(hp, i+1, j+1, VDB_x, VDB_y, vflags); + comp_fftri(&fftr, hp, i, j+1, VDB_Y, vflags); if (hessrow != NULL) { comp_hessian(hesscol, &fftr, hp->rp->ron); rev_hessian(hessdia); @@ -450,12 +684,13 @@ ambHessian( /* anisotropic radii & pos. gradient */ /* release row buffers */ if (hessrow != NULL) free(hessrow); if (gradrow != NULL) free(gradrow); + free(vflags); if (ra != NULL) /* extract eigenvectors & radii */ eigenvectors(uv, ra, hessian); - if (pg != NULL) { /* tangential position gradient/PI */ - pg[0] = DOT(gradient, uv[0]) / PI; - pg[1] = DOT(gradient, uv[1]) / PI; + if (pg != NULL) { /* tangential position gradient */ + pg[0] = DOT(gradient, uv[0]); + pg[1] = DOT(gradient, uv[1]); } } @@ -464,11 +699,11 @@ ambHessian( /* anisotropic radii & pos. gradient */ static void ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) { - struct s_ambsamp *ap; - double dgsum[2]; - int n; - FVECT vd; - double gfact; + AMBSAMP *ap; + double dgsum[2]; + int n; + FVECT vd; + double gfact; dgsum[0] = dgsum[1] = 0.0; /* sum values times -tan(theta) */ for (ap = hp->sa, n = hp->ns*hp->ns; n--; ap++) { @@ -476,15 +711,63 @@ ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) VSUB(vd, ap->p, hp->rp->rop); /* brightness over cosine factor */ gfact = colval(ap->v,CIEY) / DOT(hp->rp->ron, vd); - /* -sine = -proj_radius/vd_length */ - dgsum[0] += DOT(uv[1], vd) * gfact; - dgsum[1] -= DOT(uv[0], vd) * gfact; + /* sine = proj_radius/vd_length */ + dgsum[0] -= DOT(uv[1], vd) * gfact; + dgsum[1] += DOT(uv[0], vd) * gfact; } dg[0] = dgsum[0] / (hp->ns*hp->ns); dg[1] = dgsum[1] / (hp->ns*hp->ns); } +/* Make sure radii don't extend beyond what we see in our periphery */ +static int +hem_radii(AMBHEMI *hp, FVECT uv[2], float ra[2]) +{ +#ifdef AHEM_MARG +#define MAXDACCUM 47 + const double hemarg = AHEM_MARG*ambacc; /* hem margin */ + float radivisor2[MAXDACCUM+1]; + int i, j, k = hp->ns/10 + 1; /* around 5%ile */ + const int n2accum = (k < MAXDACCUM) ? k : MAXDACCUM ; + int na = 0; + double d; + /* circle around perimeter */ + for (i = 0; i < hp->ns; i++) + for (j = 0; j < hp->ns; j += !i|(i==hp->ns-1) ? 1 : hp->ns-1) { + AMBSAMP *ap = &ambsam(hp,i,j); + double radiv2 = 0; + FVECT vec; + if (ap->d <= FTINY) + continue; + VSUB(vec, ap->p, hp->rp->rop); + for (k = 2; k--; ) { + d = ap->d * DOT(vec, uv[k]) * ra[k]; + radiv2 += d*d; + } + radiv2 *= hemarg*hemarg * ap->d * ap->d; + if (radiv2 <= 1.0) + continue; + /* insert in percentile list */ + for (k = na; k && radiv2 > radivisor2[k-1]; k--) + radivisor2[k] = radivisor2[k-1]; + radivisor2[k] = radiv2; + na += (na < n2accum); + } + if (na < n2accum) /* current radii are OK? */ + return(0); + /* else apply divisor */ + d = 1.0/sqrt(radivisor2[na-1]); + ra[0] *= d; + ra[1] *= d; + return(1); +#undef MAXDACCUM +#else + return(0); +#endif +} + + int doambient( /* compute ambient component */ COLOR rcol, /* input/output color */ @@ -496,12 +779,12 @@ doambient( /* compute ambient component */ float dg[2] /* returned (optional) */ ) { - AMBHEMI *hp = inithemi(rcol, r, wt); - int cnt = 0; - FVECT my_uv[2]; - double d, acol[3]; - struct s_ambsamp *ap; - int i, j; + AMBHEMI *hp = inithemi(rcol, r, wt); + int cnt; + FVECT my_uv[2]; + double d, K, acol[3]; + AMBSAMP *ap; + int i, j; /* check/initialize */ if (hp == NULL) return(0); @@ -515,6 +798,7 @@ doambient( /* compute ambient component */ dg[0] = dg[1] = 0.0; /* sample the hemisphere */ acol[0] = acol[1] = acol[2] = 0.0; + cnt = 0; for (i = hp->ns; i--; ) for (j = hp->ns; j--; ) if ((ap = ambsample(hp, i, j)) != NULL) { @@ -526,20 +810,30 @@ doambient( /* compute ambient component */ free(hp); return(0); /* no valid samples */ } + if (cnt < hp->ns*hp->ns) { /* incomplete sampling? */ + copycolor(rcol, acol); + free(hp); + return(-1); /* return value w/o Hessian */ + } + cnt = ambssamp*wt + 0.5; /* perform super-sampling? */ + if (cnt > 0) + ambsupersamp(acol, hp, cnt); copycolor(rcol, acol); /* final indirect irradiance/PI */ - if (cnt < hp->ns*hp->ns || /* incomplete sampling? */ - (ra == NULL) & (pg == NULL) & (dg == NULL)) { + if ((ra == NULL) & (pg == NULL) & (dg == NULL)) { free(hp); return(-1); /* no radius or gradient calc. */ } - multcolor(acol, hp->acoef); /* normalize Y values */ - if ((d = bright(acol)) > FTINY) - d = 1.0/d; - else - d = 0.0; + if ((d = bright(acol)) > FTINY) { /* normalize Y values */ + d = 0.99*(hp->ns*hp->ns)/d; + K = 0.01; + } else { /* or fall back on geometric Hessian */ + K = 1.0; + pg = NULL; + dg = NULL; + } ap = hp->sa; /* relative Y channel from here on... */ for (i = hp->ns*hp->ns; i--; ap++) - colval(ap->v,CIEY) = bright(ap->v)*d + 0.0314; + colval(ap->v,CIEY) = bright(ap->v)*d + K; if (uv == NULL) /* make sure we have axis pointers */ uv = my_uv; @@ -550,6 +844,15 @@ doambient( /* compute ambient component */ ambdirgrad(hp, uv, dg); if (ra != NULL) { /* scale/clamp radii */ + if (pg != NULL) { + if (ra[0]*(d = fabs(pg[0])) > 1.0) + ra[0] = 1.0/d; + if (ra[1]*(d = fabs(pg[1])) > 1.0) + ra[1] = 1.0/d; + if (ra[0] > ra[1]) + ra[0] = ra[1]; + } + hem_radii(hp, uv, ra); if (ra[0] < minarad) { ra[0] = minarad; if (ra[1] < minarad) @@ -562,6 +865,14 @@ doambient( /* compute ambient component */ ra[1] = maxarad; if (ra[0] > maxarad) ra[0] = maxarad; + } + if (pg != NULL) { /* cap gradient if necessary */ + d = pg[0]*pg[0]*ra[0]*ra[0] + pg[1]*pg[1]*ra[1]*ra[1]; + if (d > 1.0) { + d = 1.0/sqrt(d); + pg[0] *= d; + pg[1] *= d; + } } } free(hp); /* clean up and return */