--- ray/src/rt/ambcomp.c 2014/04/30 18:27:14 2.40 +++ 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.40 2014/04/30 18:27:14 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.40 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.40 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, 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; @@ -87,66 +158,247 @@ 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) - goto badsample; + 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--; - /* limit vertex distance */ + 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; - else if (ar.rt <= FTINY) /* should never happen! */ - goto badsample; VSUM(ap->p, ar.rorg, ar.rdir, ar.rt); - multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ copycolor(ap->v, ar.rcol); return(ap); -badsample: - setcolor(ap->v, 0., 0., 0.); - VCOPY(ap->p, hp->rp->rop); - return(NULL); } +/* 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, FVECT ap0, FVECT ap1, FVECT rop) +comp_fftri(FFTRI *ftp, AMBHEMI *hp, int i, int j, int dbit, const uby8 *vflags) { - double rdot_cp, dot_e, dot_er, rdot_r, rdot_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); + 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); @@ -158,8 +410,9 @@ comp_fftri(FFTRI *ftp, FVECT ap0, FVECT ap1, FVECT rop 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 (i = 3; i--; ) - ftp->rI2_eJ2[i] = ftp->I2*ftp->r_i[i] + J2*ftp->e_i[i]; + for (ii = 3; ii--; ) + ftp->rI2_eJ2[ii] = ftp->I2*ftp->r_i[ii] + J2*ftp->e_i[ii]; + ftp->valid++; } @@ -185,6 +438,11 @@ comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) 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); @@ -208,7 +466,7 @@ comp_hessian(FVECT hess[3], FFTRI *ftp, FVECT nrm) hess[i][j] = m1[i][j] + d1*( I3*m2[i][j] + K3*m3[i][j] + 2.0*J3*m4[i][j] ); hess[i][j] += d2*(i==j); - hess[i][j] *= 1.0/PI; + hess[i][j] *= -1.0/PI; } } @@ -230,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; @@ -248,10 +506,14 @@ comp_gradient(FVECT grad, FFTRI *ftp, FVECT nrm) double f1; int 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*ncp[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]) @@ -362,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; @@ -383,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) @@ -396,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) @@ -405,12 +642,10 @@ 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); @@ -422,17 +657,14 @@ ambHessian( /* anisotropic radii & pos. gradient */ 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); @@ -452,6 +684,7 @@ 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); @@ -466,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++) { @@ -487,6 +720,54 @@ ambdirgrad(AMBHEMI *hp, FVECT uv[2], float dg[2]) } +/* 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 */ @@ -498,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, K, 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); @@ -517,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) { @@ -528,17 +810,23 @@ 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. */ } - if (bright(acol) > FTINY) { /* normalize Y values */ - d = 0.99*cnt/bright(acol); + if ((d = bright(acol)) > FTINY) { /* normalize Y values */ + d = 0.99*(hp->ns*hp->ns)/d; K = 0.01; - } else { /* geometric Hessian fall-back */ - d = 0.0; + } else { /* or fall back on geometric Hessian */ K = 1.0; pg = NULL; dg = NULL; @@ -564,6 +852,7 @@ doambient( /* compute ambient component */ 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)