Back to index

radiance  4R0+20100331
ambcomp.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: ambcomp.c,v 2.21 2007/12/31 18:19:42 greg Exp $";
00003 #endif
00004 /*
00005  * Routines to compute "ambient" values using Monte Carlo
00006  *
00007  *  Declarations of external symbols in ambient.h
00008  */
00009 
00010 #include "copyright.h"
00011 
00012 #include  "ray.h"
00013 
00014 #include  "ambient.h"
00015 
00016 #include  "random.h"
00017 
00018 
00019 void
00020 inithemi(                   /* initialize sampling hemisphere */
00021        register AMBHEMI  *hp,
00022        COLOR ac,
00023        RAY  *r,
00024        double  wt
00025 )
00026 {
00027        double d;
00028        register int  i;
00029                                    /* set number of divisions */
00030        if (ambacc <= FTINY &&
00031                      wt > (d = 0.8*intens(ac)*r->rweight/(ambdiv*minweight)))
00032               wt = d;                     /* avoid ray termination */
00033        hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
00034        i = ambacc > FTINY ? 3 : 1; /* minimum number of samples */
00035        if (hp->nt < i)
00036               hp->nt = i;
00037        hp->np = PI * hp->nt + 0.5;
00038                                    /* set number of super-samples */
00039        hp->ns = ambssamp * wt + 0.5;
00040                                    /* assign coefficient */
00041        copycolor(hp->acoef, ac);
00042        d = 1.0/(hp->nt*hp->np);
00043        scalecolor(hp->acoef, d);
00044                                    /* make axes */
00045        VCOPY(hp->uz, r->ron);
00046        hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
00047        for (i = 0; i < 3; i++)
00048               if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
00049                      break;
00050        if (i >= 3)
00051               error(CONSISTENCY, "bad ray direction in inithemi");
00052        hp->uy[i] = 1.0;
00053        fcross(hp->ux, hp->uy, hp->uz);
00054        normalize(hp->ux);
00055        fcross(hp->uy, hp->uz, hp->ux);
00056 }
00057 
00058 
00059 int
00060 divsample(                         /* sample a division */
00061        register AMBSAMP  *dp,
00062        AMBHEMI  *h,
00063        RAY  *r
00064 )
00065 {
00066        RAY  ar;
00067        int  hlist[3];
00068        double  spt[2];
00069        double  xd, yd, zd;
00070        double  b2;
00071        double  phi;
00072        register int  i;
00073                                    /* ambient coefficient for weight */
00074        if (ambacc > FTINY)
00075               setcolor(ar.rcoef, AVGREFL, AVGREFL, AVGREFL);
00076        else
00077               copycolor(ar.rcoef, h->acoef);
00078        if (rayorigin(&ar, AMBIENT, r, ar.rcoef) < 0)
00079               return(-1);
00080        if (ambacc > FTINY) {
00081               multcolor(ar.rcoef, h->acoef);
00082               scalecolor(ar.rcoef, 1./AVGREFL);
00083        }
00084        hlist[0] = r->rno;
00085        hlist[1] = dp->t;
00086        hlist[2] = dp->p;
00087        multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
00088        zd = sqrt((dp->t + spt[0])/h->nt);
00089        phi = 2.0*PI * (dp->p + spt[1])/h->np;
00090        xd = tcos(phi) * zd;
00091        yd = tsin(phi) * zd;
00092        zd = sqrt(1.0 - zd*zd);
00093        for (i = 0; i < 3; i++)
00094               ar.rdir[i] =  xd*h->ux[i] +
00095                             yd*h->uy[i] +
00096                             zd*h->uz[i];
00097        dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
00098        rayvalue(&ar);
00099        ndims--;
00100        multcolor(ar.rcol, ar.rcoef);      /* apply coefficient */
00101        addcolor(dp->v, ar.rcol);
00102                                    /* use rt to improve gradient calc */
00103        if (ar.rt > FTINY && ar.rt < FHUGE)
00104               dp->r += 1.0/ar.rt;
00105                                    /* (re)initialize error */
00106        if (dp->n++) {
00107               b2 = bright(dp->v)/dp->n - bright(ar.rcol);
00108               b2 = b2*b2 + dp->k*((dp->n-1)*(dp->n-1));
00109               dp->k = b2/(dp->n*dp->n);
00110        } else
00111               dp->k = 0.0;
00112        return(0);
00113 }
00114 
00115 
00116 static int
00117 ambcmp(                                   /* decreasing order */
00118        const void *p1,
00119        const void *p2
00120 )
00121 {
00122        const AMBSAMP *d1 = (const AMBSAMP *)p1;
00123        const AMBSAMP *d2 = (const AMBSAMP *)p2;
00124 
00125        if (d1->k < d2->k)
00126               return(1);
00127        if (d1->k > d2->k)
00128               return(-1);
00129        return(0);
00130 }
00131 
00132 
00133 static int
00134 ambnorm(                           /* standard order */
00135        const void *p1,
00136        const void *p2
00137 )
00138 {
00139        const AMBSAMP *d1 = (const AMBSAMP *)p1;
00140        const AMBSAMP *d2 = (const AMBSAMP *)p2;
00141        register int  c;
00142 
00143        if ( (c = d1->t - d2->t) )
00144               return(c);
00145        return(d1->p - d2->p);
00146 }
00147 
00148 
00149 double
00150 doambient(                         /* compute ambient component */
00151        COLOR  acol,
00152        RAY  *r,
00153        double  wt,
00154        FVECT  pg,
00155        FVECT  dg
00156 )
00157 {
00158        double  b, d;
00159        AMBHEMI  hemi;
00160        AMBSAMP  *div;
00161        AMBSAMP  dnew;
00162        register AMBSAMP  *dp;
00163        double  arad;
00164        int  divcnt;
00165        register int  i, j;
00166                                    /* initialize hemisphere */
00167        inithemi(&hemi, acol, r, wt);
00168        divcnt = hemi.nt * hemi.np;
00169                                    /* initialize */
00170        if (pg != NULL)
00171               pg[0] = pg[1] = pg[2] = 0.0;
00172        if (dg != NULL)
00173               dg[0] = dg[1] = dg[2] = 0.0;
00174        setcolor(acol, 0.0, 0.0, 0.0);
00175        if (divcnt == 0)
00176               return(0.0);
00177                                    /* allocate super-samples */
00178        if (hemi.ns > 0 || pg != NULL || dg != NULL) {
00179               div = (AMBSAMP *)malloc(divcnt*sizeof(AMBSAMP));
00180               if (div == NULL)
00181                      error(SYSTEM, "out of memory in doambient");
00182        } else
00183               div = NULL;
00184                                    /* sample the divisions */
00185        arad = 0.0;
00186        if ((dp = div) == NULL)
00187               dp = &dnew;
00188        divcnt = 0;
00189        for (i = 0; i < hemi.nt; i++)
00190               for (j = 0; j < hemi.np; j++) {
00191                      dp->t = i; dp->p = j;
00192                      setcolor(dp->v, 0.0, 0.0, 0.0);
00193                      dp->r = 0.0;
00194                      dp->n = 0;
00195                      if (divsample(dp, &hemi, r) < 0) {
00196                             if (div != NULL)
00197                                    dp++;
00198                             continue;
00199                      }
00200                      arad += dp->r;
00201                      divcnt++;
00202                      if (div != NULL)
00203                             dp++;
00204                      else
00205                             addcolor(acol, dp->v);
00206               }
00207        if (!divcnt) {
00208               if (div != NULL)
00209                      free((void *)div);
00210               return(0.0);         /* no samples taken */
00211        }
00212        if (divcnt < hemi.nt*hemi.np) {
00213               pg = dg = NULL;             /* incomplete sampling */
00214               hemi.ns = 0;
00215        } else if (arad > FTINY && divcnt/arad < minarad) {
00216               hemi.ns = 0;         /* close enough */
00217        } else if (hemi.ns > 0) {   /* else perform super-sampling? */
00218               comperrs(div, &hemi);                     /* compute errors */
00219               qsort(div, divcnt, sizeof(AMBSAMP), ambcmp);     /* sort divs */
00220                                           /* super-sample */
00221               for (i = hemi.ns; i > 0; i--) {
00222                      dnew = *div;
00223                      if (divsample(&dnew, &hemi, r) < 0) {
00224                             dp++;
00225                             continue;
00226                      }
00227                      dp = div;            /* reinsert */
00228                      j = divcnt < i ? divcnt : i;
00229                      while (--j > 0 && dnew.k < dp[1].k) {
00230                             *dp = *(dp+1);
00231                             dp++;
00232                      }
00233                      *dp = dnew;
00234               }
00235               if (pg != NULL || dg != NULL)      /* restore order */
00236                      qsort(div, divcnt, sizeof(AMBSAMP), ambnorm);
00237        }
00238                                    /* compute returned values */
00239        if (div != NULL) {
00240               arad = 0.0;          /* note: divcnt may be < nt*np */
00241               for (i = hemi.nt*hemi.np, dp = div; i-- > 0; dp++) {
00242                      arad += dp->r;
00243                      if (dp->n > 1) {
00244                             b = 1.0/dp->n;
00245                             scalecolor(dp->v, b);
00246                             dp->r *= b;
00247                             dp->n = 1;
00248                      }
00249                      addcolor(acol, dp->v);
00250               }
00251               b = bright(acol);
00252               if (b > FTINY) {
00253                      b = 1.0/b;    /* compute & normalize gradient(s) */
00254                      if (pg != NULL) {
00255                             posgradient(pg, div, &hemi);
00256                             for (i = 0; i < 3; i++)
00257                                    pg[i] *= b;
00258                      }
00259                      if (dg != NULL) {
00260                             dirgradient(dg, div, &hemi);
00261                             for (i = 0; i < 3; i++)
00262                                    dg[i] *= b;
00263                      }
00264               }
00265               free((void *)div);
00266        }
00267        if (arad <= FTINY)
00268               arad = maxarad;
00269        else
00270               arad = (divcnt+hemi.ns)/arad;
00271        if (pg != NULL) {           /* reduce radius if gradient large */
00272               d = DOT(pg,pg);
00273               if (d*arad*arad > 1.0)
00274                      arad = 1.0/sqrt(d);
00275        }
00276        if (arad < minarad) {
00277               arad = minarad;
00278               if (pg != NULL && d*arad*arad > 1.0) {    /* cap gradient */
00279                      d = 1.0/arad/sqrt(d);
00280                      for (i = 0; i < 3; i++)
00281                             pg[i] *= d;
00282               }
00283        }
00284        if ((arad /= sqrt(wt)) > maxarad)
00285               arad = maxarad;
00286        return(arad);
00287 }
00288 
00289 
00290 void
00291 comperrs(                   /* compute initial error estimates */
00292        AMBSAMP  *da, /* assumes standard ordering */
00293        register AMBHEMI  *hp
00294 )
00295 {
00296        double  b, b2;
00297        int  i, j;
00298        register AMBSAMP  *dp;
00299                             /* sum differences from neighbors */
00300        dp = da;
00301        for (i = 0; i < hp->nt; i++)
00302               for (j = 0; j < hp->np; j++) {
00303 #ifdef  DEBUG
00304                      if (dp->t != i || dp->p != j)
00305                             error(CONSISTENCY,
00306                                    "division order in comperrs");
00307 #endif
00308                      b = bright(dp[0].v);
00309                      if (i > 0) {         /* from above */
00310                             b2 = bright(dp[-hp->np].v) - b;
00311                             b2 *= b2 * 0.25;
00312                             dp[0].k += b2;
00313                             dp[-hp->np].k += b2;
00314                      }
00315                      if (j > 0) {         /* from behind */
00316                             b2 = bright(dp[-1].v) - b;
00317                             b2 *= b2 * 0.25;
00318                             dp[0].k += b2;
00319                             dp[-1].k += b2;
00320                      } else {             /* around */
00321                             b2 = bright(dp[hp->np-1].v) - b;
00322                             b2 *= b2 * 0.25;
00323                             dp[0].k += b2;
00324                             dp[hp->np-1].k += b2;
00325                      }
00326                      dp++;
00327               }
00328                             /* divide by number of neighbors */
00329        dp = da;
00330        for (j = 0; j < hp->np; j++)              /* top row */
00331               (dp++)->k *= 1.0/3.0;
00332        if (hp->nt < 2)
00333               return;
00334        for (i = 1; i < hp->nt-1; i++)            /* central region */
00335               for (j = 0; j < hp->np; j++)
00336                      (dp++)->k *= 0.25;
00337        for (j = 0; j < hp->np; j++)              /* bottom row */
00338               (dp++)->k *= 1.0/3.0;
00339 }
00340 
00341 
00342 void
00343 posgradient(                              /* compute position gradient */
00344        FVECT  gv,
00345        AMBSAMP  *da,               /* assumes standard ordering */
00346        register AMBHEMI  *hp
00347 )
00348 {
00349        register int  i, j;
00350        double  nextsine, lastsine, b, d;
00351        double  mag0, mag1;
00352        double  phi, cosp, sinp, xd, yd;
00353        register AMBSAMP  *dp;
00354 
00355        xd = yd = 0.0;
00356        for (j = 0; j < hp->np; j++) {
00357               dp = da + j;
00358               mag0 = mag1 = 0.0;
00359               lastsine = 0.0;
00360               for (i = 0; i < hp->nt; i++) {
00361 #ifdef  DEBUG
00362                      if (dp->t != i || dp->p != j)
00363                             error(CONSISTENCY,
00364                                    "division order in posgradient");
00365 #endif
00366                      b = bright(dp->v);
00367                      if (i > 0) {
00368                             d = dp[-hp->np].r;
00369                             if (dp[0].r > d) d = dp[0].r;
00370                                                  /* sin(t)*cos(t)^2 */
00371                             d *= lastsine * (1.0 - (double)i/hp->nt);
00372                             mag0 += d*(b - bright(dp[-hp->np].v));
00373                      }
00374                      nextsine = sqrt((double)(i+1)/hp->nt);
00375                      if (j > 0) {
00376                             d = dp[-1].r;
00377                             if (dp[0].r > d) d = dp[0].r;
00378                             mag1 += d * (nextsine - lastsine) *
00379                                           (b - bright(dp[-1].v));
00380                      } else {
00381                             d = dp[hp->np-1].r;
00382                             if (dp[0].r > d) d = dp[0].r;
00383                             mag1 += d * (nextsine - lastsine) *
00384                                           (b - bright(dp[hp->np-1].v));
00385                      }
00386                      dp += hp->np;
00387                      lastsine = nextsine;
00388               }
00389               mag0 *= 2.0*PI / hp->np;
00390               phi = 2.0*PI * (double)j/hp->np;
00391               cosp = tcos(phi); sinp = tsin(phi);
00392               xd += mag0*cosp - mag1*sinp;
00393               yd += mag0*sinp + mag1*cosp;
00394        }
00395        for (i = 0; i < 3; i++)
00396               gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])*(hp->nt*hp->np)/PI;
00397 }
00398 
00399 
00400 void
00401 dirgradient(                              /* compute direction gradient */
00402        FVECT  gv,
00403        AMBSAMP  *da,               /* assumes standard ordering */
00404        register AMBHEMI  *hp
00405 )
00406 {
00407        register int  i, j;
00408        double  mag;
00409        double  phi, xd, yd;
00410        register AMBSAMP  *dp;
00411 
00412        xd = yd = 0.0;
00413        for (j = 0; j < hp->np; j++) {
00414               dp = da + j;
00415               mag = 0.0;
00416               for (i = 0; i < hp->nt; i++) {
00417 #ifdef  DEBUG
00418                      if (dp->t != i || dp->p != j)
00419                             error(CONSISTENCY,
00420                                    "division order in dirgradient");
00421 #endif
00422                                                  /* tan(t) */
00423                      mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
00424                      dp += hp->np;
00425               }
00426               phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
00427               xd += mag * tcos(phi);
00428               yd += mag * tsin(phi);
00429        }
00430        for (i = 0; i < 3; i++)
00431               gv[i] = xd*hp->ux[i] + yd*hp->uy[i];
00432 }