--- ray/src/rt/ambcomp.c 2005/05/28 22:27:54 2.15 +++ ray/src/rt/ambcomp.c 2005/10/28 16:16:33 2.20 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: ambcomp.c,v 2.15 2005/05/28 22:27:54 greg Exp $"; +static const char RCSid[] = "$Id: ambcomp.c,v 2.20 2005/10/28 16:16:33 greg Exp $"; #endif /* * Routines to compute "ambient" values using Monte Carlo @@ -19,26 +19,28 @@ static const char RCSid[] = "$Id: ambcomp.c,v 2.15 200 void inithemi( /* initialize sampling hemisphere */ register AMBHEMI *hp, + COLOR ac, RAY *r, - COLOR ac, double wt ) { + double d; register int i; /* set number of divisions */ - hp->nt = sqrt(ambdiv * wt * (1./PI/AVGREFL)) + 0.5; + if (ambacc <= FTINY && + wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight))) + wt = d; /* avoid ray termination */ + hp->nt = sqrt(ambdiv * wt / PI) + 0.5; i = ambacc > FTINY ? 3 : 1; /* minimum number of samples */ if (hp->nt < i) hp->nt = i; hp->np = PI * hp->nt + 0.5; /* set number of super-samples */ hp->ns = ambssamp * wt + 0.5; - /* assign coefficients */ + /* assign coefficient */ copycolor(hp->acoef, ac); - if (wt >= r->rweight) - hp->drc = 1.; - else - hp->drc = wt / r->rweight; + d = 1.0/(hp->nt*hp->np); + scalecolor(hp->acoef, d); /* make axes */ VCOPY(hp->uz, r->ron); hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0; @@ -69,12 +71,16 @@ divsample( /* sample a division */ double phi; register int i; /* ambient coefficient for weight */ - setcolor(ar.rcoef, h->drc, h->drc, h->drc); + if (ambacc > FTINY) + setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL); + else + copycolor(ar.rcoef, h->acoef); if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0) return(-1); - copycolor(ar.rcoef, h->acoef); /* correct coefficient for trace */ - b2 = 1.0/(h->nt*h->np + h->ns); /* XXX not uniform if ns > 0 */ - scalecolor(ar.rcoef, b2); + if (ambacc > FTINY) { + multcolor(ar.rcoef, h->acoef); + scalecolor(ar.rcoef, 1./AVGREFL); + } hlist[0] = r->rno; hlist[1] = dp->t; hlist[2] = dp->p; @@ -91,6 +97,7 @@ divsample( /* sample a division */ dimlist[ndims++] = dp->t*h->np + dp->p + 90171; rayvalue(&ar); ndims--; + multcolor(ar.rcol, ar.rcoef); /* apply coefficient */ addcolor(dp->v, ar.rcol); /* use rt to improve gradient calc */ if (ar.rt > FTINY && ar.rt < FHUGE) @@ -143,7 +150,6 @@ double doambient( /* compute ambient component */ COLOR acol, RAY *r, - COLOR ac, double wt, FVECT pg, FVECT dg @@ -155,18 +161,22 @@ doambient( /* compute ambient component */ AMBSAMP dnew; register AMBSAMP *dp; double arad; - int ndivs; + int divcnt; register int i, j; - /* initialize color */ - setcolor(acol, 0.0, 0.0, 0.0); /* initialize hemisphere */ - inithemi(&hemi, r, ac, wt); - ndivs = hemi.nt * hemi.np; - if (ndivs == 0) + inithemi(&hemi, acol, r, wt); + divcnt = hemi.nt * hemi.np; + /* initialize */ + if (pg != NULL) + pg[0] = pg[1] = pg[2] = 0.0; + if (dg != NULL) + dg[0] = dg[1] = dg[2] = 0.0; + setcolor(acol, 0.0, 0.0, 0.0); + if (divcnt == 0) return(0.0); /* allocate super-samples */ if (hemi.ns > 0 || pg != NULL || dg != NULL) { - div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP)); + div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP)); if (div == NULL) error(SYSTEM, "out of memory in doambient"); } else @@ -175,33 +185,44 @@ doambient( /* compute ambient component */ arad = 0.0; if ((dp = div) == NULL) dp = &dnew; + divcnt = 0; for (i = 0; i < hemi.nt; i++) for (j = 0; j < hemi.np; j++) { dp->t = i; dp->p = j; setcolor(dp->v, 0.0, 0.0, 0.0); dp->r = 0.0; dp->n = 0; - if (divsample(dp, &hemi, r) < 0) - goto oopsy; + if (divsample(dp, &hemi, r) < 0) { + if (div != NULL) + dp++; + continue; + } arad += dp->r; + divcnt++; if (div != NULL) dp++; else addcolor(acol, dp->v); } - if (hemi.ns > 0 && arad > FTINY && ndivs/arad < minarad) + if (!divcnt) + return(0.0); /* no samples taken */ + if (divcnt < hemi.nt*hemi.np) { + pg = dg = NULL; /* incomplete sampling */ + hemi.ns = 0; + } else if (arad > FTINY && divcnt/arad < minarad) { hemi.ns = 0; /* close enough */ - else if (hemi.ns > 0) { /* else perform super-sampling */ + } else if (hemi.ns > 0) { /* else perform super-sampling? */ comperrs(div, &hemi); /* compute errors */ - qsort(div, ndivs, sizeof(AMBSAMP), ambcmp); /* sort divs */ + qsort(div, divcnt, sizeof(AMBSAMP), ambcmp); /* sort divs */ /* super-sample */ for (i = hemi.ns; i > 0; i--) { dnew = *div; - if (divsample(&dnew, &hemi, r) < 0) - goto oopsy; - /* reinsert */ - dp = div; - j = ndivs < i ? ndivs : i; + if (divsample(&dnew, &hemi, r) < 0) { + dp++; + continue; + } + dp = div; /* reinsert */ + j = divcnt < i ? divcnt : i; while (--j > 0 && dnew.k < dp[1].k) { *dp = *(dp+1); dp++; @@ -209,12 +230,12 @@ doambient( /* compute ambient component */ *dp = dnew; } if (pg != NULL || dg != NULL) /* restore order */ - qsort(div, ndivs, sizeof(AMBSAMP), ambnorm); + qsort(div, divcnt, sizeof(AMBSAMP), ambnorm); } /* compute returned values */ if (div != NULL) { - arad = 0.0; - for (i = ndivs, dp = div; i-- > 0; dp++) { + arad = 0.0; /* note: divcnt may be < nt*np */ + for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) { arad += dp->r; if (dp->n > 1) { b = 1.0/dp->n; @@ -226,7 +247,7 @@ doambient( /* compute ambient component */ } b = bright(acol); if (b > FTINY) { - b = ndivs/b; + b = 1.0/b; /* compute & normalize gradient(s) */ if (pg != NULL) { posgradient(pg, div, &hemi); for (i = 0; i < 3; i++) @@ -237,22 +258,13 @@ doambient( /* compute ambient component */ for (i = 0; i < 3; i++) dg[i] *= b; } - } else { - if (pg != NULL) - for (i = 0; i < 3; i++) - pg[i] = 0.0; - if (dg != NULL) - for (i = 0; i < 3; i++) - dg[i] = 0.0; } free((void *)div); } - b = 1.0/ndivs; - scalecolor(acol, b); if (arad <= FTINY) arad = maxarad; else - arad = (ndivs+hemi.ns)/arad; + arad = (divcnt+hemi.ns)/arad; if (pg != NULL) { /* reduce radius if gradient large */ d = DOT(pg,pg); if (d*arad*arad > 1.0) @@ -269,10 +281,6 @@ doambient( /* compute ambient component */ if ((arad /= sqrt(wt)) > maxarad) arad = maxarad; return(arad); -oopsy: - if (div != NULL) - free((void *)div); - return(0.0); } @@ -382,7 +390,7 @@ posgradient( /* compute position gradient */ yd += mag0*sinp + mag1*cosp; } for (i = 0; i < 3; i++) - gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])/PI; + gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI; } @@ -417,5 +425,5 @@ dirgradient( /* compute direction gradient */ yd += mag * tsin(phi); } for (i = 0; i < 3; i++) - gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])/(hp->nt*hp->np); + gv[i] = xd*hp->ux[i] + yd*hp->uy[i]; }