home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / radsrc22 / src / rt / ambcomp.c next >
Encoding:
C/C++ Source or Header  |  1992-10-02  |  8.9 KB  |  398 lines
  1. /* Copyright (c) 1991 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)ambcomp.c 2.4 10/2/92 LBL";
  5. #endif
  6.  
  7. /*
  8.  * Routines to compute "ambient" values using Monte Carlo
  9.  */
  10.  
  11. #include  "ray.h"
  12.  
  13. #include  "ambient.h"
  14.  
  15. #include  "random.h"
  16.  
  17. typedef struct {
  18.     short  t, p;        /* theta, phi indices */
  19.     COLOR  v;        /* value sum */
  20.     float  r;        /* 1/distance sum */
  21.     float  k;        /* variance for this division */
  22.     int  n;            /* number of subsamples */
  23. }  AMBSAMP;        /* ambient sample division */
  24.  
  25. typedef struct {
  26.     FVECT  ux, uy, uz;    /* x, y and z axis directions */
  27.     short  nt, np;        /* number of theta and phi directions */
  28. }  AMBHEMI;        /* ambient sample hemisphere */
  29.  
  30.  
  31. static int
  32. ambcmp(d1, d2)                /* decreasing order */
  33. AMBSAMP  *d1, *d2;
  34. {
  35.     if (d1->k < d2->k)
  36.         return(1);
  37.     if (d1->k > d2->k)
  38.         return(-1);
  39.     return(0);
  40. }
  41.  
  42.  
  43. static int
  44. ambnorm(d1, d2)                /* standard order */
  45. AMBSAMP  *d1, *d2;
  46. {
  47.     register int  c;
  48.  
  49.     if (c = d1->t - d2->t)
  50.         return(c);
  51.     return(d1->p - d2->p);
  52. }
  53.  
  54.  
  55. divsample(dp, h, r)            /* sample a division */
  56. register AMBSAMP  *dp;
  57. AMBHEMI  *h;
  58. RAY  *r;
  59. {
  60.     RAY  ar;
  61.     int  hlist[3];
  62.     double  spt[2];
  63.     double  xd, yd, zd;
  64.     double  b2;
  65.     double  phi;
  66.     register int  i;
  67.  
  68.     if (rayorigin(&ar, r, AMBIENT, AVGREFL) < 0)
  69.         return(-1);
  70.     hlist[0] = r->rno;
  71.     hlist[1] = dp->t;
  72.     hlist[2] = dp->p;
  73.     multisamp(spt, 2, urand(ilhash(hlist,3)+dp->n));
  74.     zd = sqrt((dp->t + spt[0])/h->nt);
  75.     phi = 2.0*PI * (dp->p + spt[1])/h->np;
  76.     xd = cos(phi) * zd;
  77.     yd = sin(phi) * zd;
  78.     zd = sqrt(1.0 - zd*zd);
  79.     for (i = 0; i < 3; i++)
  80.         ar.rdir[i] =    xd*h->ux[i] +
  81.                 yd*h->uy[i] +
  82.                 zd*h->uz[i];
  83.     dimlist[ndims++] = dp->t*h->np + dp->p + 90171;
  84.     rayvalue(&ar);
  85.     ndims--;
  86.     addcolor(dp->v, ar.rcol);
  87.     if (ar.rt > FTINY && ar.rt < FHUGE)
  88.         dp->r += 1.0/ar.rt;
  89.                     /* (re)initialize error */
  90.     if (dp->n++) {
  91.         b2 = bright(dp->v)/dp->n - bright(ar.rcol);
  92.         b2 = b2*b2 + dp->k*((dp->n-1)*(dp->n-1));
  93.         dp->k = b2/(dp->n*dp->n);
  94.     } else
  95.         dp->k = 0.0;
  96.     return(0);
  97. }
  98.  
  99.  
  100. double
  101. doambient(acol, r, wt, pg, dg)        /* compute ambient component */
  102. COLOR  acol;
  103. RAY  *r;
  104. double  wt;
  105. FVECT  pg, dg;
  106. {
  107.     double  b, d;
  108.     AMBHEMI  hemi;
  109.     AMBSAMP  *div;
  110.     AMBSAMP  dnew;
  111.     register AMBSAMP  *dp;
  112.     double  arad;
  113.     int  ndivs, ns;
  114.     register int  i, j;
  115.                     /* initialize color */
  116.     setcolor(acol, 0.0, 0.0, 0.0);
  117.                     /* initialize hemisphere */
  118.     inithemi(&hemi, r, wt);
  119.     ndivs = hemi.nt * hemi.np;
  120.     if (ndivs == 0)
  121.         return(0.0);
  122.                     /* set number of super-samples */
  123.     ns = ambssamp * wt + 0.5;
  124.     if (ns > 0 || pg != NULL || dg != NULL) {
  125.         div = (AMBSAMP *)malloc(ndivs*sizeof(AMBSAMP));
  126.         if (div == NULL)
  127.             error(SYSTEM, "out of memory in doambient");
  128.     } else
  129.         div = NULL;
  130.                     /* sample the divisions */
  131.     arad = 0.0;
  132.     if ((dp = div) == NULL)
  133.         dp = &dnew;
  134.     for (i = 0; i < hemi.nt; i++)
  135.         for (j = 0; j < hemi.np; j++) {
  136.             dp->t = i; dp->p = j;
  137.             setcolor(dp->v, 0.0, 0.0, 0.0);
  138.             dp->r = 0.0;
  139.             dp->n = 0;
  140.             if (divsample(dp, &hemi, r) < 0)
  141.                 goto oopsy;
  142.             if (div != NULL)
  143.                 dp++;
  144.             else {
  145.                 addcolor(acol, dp->v);
  146.                 arad += dp->r;
  147.             }
  148.         }
  149.     if (ns > 0) {            /* perform super-sampling */
  150.         comperrs(div, &hemi);            /* compute errors */
  151.         qsort(div, ndivs, sizeof(AMBSAMP), ambcmp);    /* sort divs */
  152.                         /* super-sample */
  153.         for (i = ns; i > 0; i--) {
  154.             copystruct(&dnew, div);
  155.             if (divsample(&dnew, &hemi, r) < 0)
  156.                 goto oopsy;
  157.                             /* reinsert */
  158.             dp = div;
  159.             j = ndivs < i ? ndivs : i;
  160.             while (--j > 0 && dnew.k < dp[1].k) {
  161.                 copystruct(dp, dp+1);
  162.                 dp++;
  163.             }
  164.             copystruct(dp, &dnew);
  165.         }
  166.         if (pg != NULL || dg != NULL)    /* restore order */
  167.             qsort(div, ndivs, sizeof(AMBSAMP), ambnorm);
  168.     }
  169.                     /* compute returned values */
  170.     if (div != NULL) {
  171.         for (i = ndivs, dp = div; i-- > 0; dp++) {
  172.             arad += dp->r;
  173.             if (dp->n > 1) {
  174.                 b = 1.0/dp->n;
  175.                 scalecolor(dp->v, b);
  176.                 dp->r *= b;
  177.                 dp->n = 1;
  178.             }
  179.             addcolor(acol, dp->v);
  180.         }
  181.         b = bright(acol);
  182.         if (b > FTINY) {
  183.             b = ndivs/b;
  184.             if (pg != NULL) {
  185.                 posgradient(pg, div, &hemi);
  186.                 for (i = 0; i < 3; i++)
  187.                     pg[i] *= b;
  188.             }
  189.             if (dg != NULL) {
  190.                 dirgradient(dg, div, &hemi);
  191.                 for (i = 0; i < 3; i++)
  192.                     dg[i] *= b;
  193.             }
  194.         } else {
  195.             if (pg != NULL)
  196.                 for (i = 0; i < 3; i++)
  197.                     pg[i] = 0.0;
  198.             if (dg != NULL)
  199.                 for (i = 0; i < 3; i++)
  200.                     dg[i] = 0.0;
  201.         }
  202.         free((char *)div);
  203.     }
  204.     b = 1.0/ndivs;
  205.     scalecolor(acol, b);
  206.     if (arad <= FTINY)
  207.         arad = maxarad;
  208.     else
  209.         arad = (ndivs+ns)/arad;
  210.     if (pg != NULL) {        /* reduce radius if gradient large */
  211.         d = DOT(pg,pg);
  212.         if (d*arad*arad > 1.0)
  213.             arad = 1.0/sqrt(d);
  214.     }
  215.     if (arad < minarad) {
  216.         arad = minarad;
  217.         if (pg != NULL && d*arad*arad > 1.0) {    /* cap gradient */
  218.             d = 1.0/arad/sqrt(d);
  219.             for (i = 0; i < 3; i++)
  220.                 pg[i] *= d;
  221.         }
  222.     }
  223.     if ((arad /= sqrt(wt)) > maxarad)
  224.         arad = maxarad;
  225.     return(arad);
  226. oopsy:
  227.     if (div != NULL)
  228.         free((char *)div);
  229.     return(0.0);
  230. }
  231.  
  232.  
  233. inithemi(hp, r, wt)        /* initialize sampling hemisphere */
  234. register AMBHEMI  *hp;
  235. RAY  *r;
  236. double  wt;
  237. {
  238.     register int  i;
  239.                     /* set number of divisions */
  240.     if (wt < (.25*PI)/ambdiv+FTINY) {
  241.         hp->nt = hp->np = 0;
  242.         return;            /* zero samples */
  243.     }
  244.     hp->nt = sqrt(ambdiv * wt / PI) + 0.5;
  245.     hp->np = PI * hp->nt + 0.5;
  246.                     /* make axes */
  247.     VCOPY(hp->uz, r->ron);
  248.     hp->uy[0] = hp->uy[1] = hp->uy[2] = 0.0;
  249.     for (i = 0; i < 3; i++)
  250.         if (hp->uz[i] < 0.6 && hp->uz[i] > -0.6)
  251.             break;
  252.     if (i >= 3)
  253.         error(CONSISTENCY, "bad ray direction in inithemi");
  254.     hp->uy[i] = 1.0;
  255.     fcross(hp->ux, hp->uy, hp->uz);
  256.     normalize(hp->ux);
  257.     fcross(hp->uy, hp->uz, hp->ux);
  258. }
  259.  
  260.  
  261. comperrs(da, hp)        /* compute initial error estimates */
  262. AMBSAMP  *da;        /* assumes standard ordering */
  263. register AMBHEMI  *hp;
  264. {
  265.     double  b, b2;
  266.     int  i, j;
  267.     register AMBSAMP  *dp;
  268.                 /* sum differences from neighbors */
  269.     dp = da;
  270.     for (i = 0; i < hp->nt; i++)
  271.         for (j = 0; j < hp->np; j++) {
  272. #ifdef  DEBUG
  273.             if (dp->t != i || dp->p != j)
  274.                 error(CONSISTENCY,
  275.                     "division order in comperrs");
  276. #endif
  277.             b = bright(dp[0].v);
  278.             if (i > 0) {        /* from above */
  279.                 b2 = bright(dp[-hp->np].v) - b;
  280.                 b2 *= b2 * 0.25;
  281.                 dp[0].k += b2;
  282.                 dp[-hp->np].k += b2;
  283.             }
  284.             if (j > 0) {        /* from behind */
  285.                 b2 = bright(dp[-1].v) - b;
  286.                 b2 *= b2 * 0.25;
  287.                 dp[0].k += b2;
  288.                 dp[-1].k += b2;
  289.             } else {        /* around */
  290.                 b2 = bright(dp[hp->np-1].v) - b;
  291.                 b2 *= b2 * 0.25;
  292.                 dp[0].k += b2;
  293.                 dp[hp->np-1].k += b2;
  294.             }
  295.             dp++;
  296.         }
  297.                 /* divide by number of neighbors */
  298.     dp = da;
  299.     for (j = 0; j < hp->np; j++)        /* top row */
  300.         (dp++)->k *= 1.0/3.0;
  301.     if (hp->nt < 2)
  302.         return;
  303.     for (i = 1; i < hp->nt-1; i++)        /* central region */
  304.         for (j = 0; j < hp->np; j++)
  305.             (dp++)->k *= 0.25;
  306.     for (j = 0; j < hp->np; j++)        /* bottom row */
  307.         (dp++)->k *= 1.0/3.0;
  308. }
  309.  
  310.  
  311. posgradient(gv, da, hp)                /* compute position gradient */
  312. FVECT  gv;
  313. AMBSAMP  *da;            /* assumes standard ordering */
  314. register AMBHEMI  *hp;
  315. {
  316.     register int  i, j;
  317.     double  nextsine, lastsine, b, d;
  318.     double  mag0, mag1;
  319.     double  phi, cosp, sinp, xd, yd;
  320.     register AMBSAMP  *dp;
  321.  
  322.     xd = yd = 0.0;
  323.     for (j = 0; j < hp->np; j++) {
  324.         dp = da + j;
  325.         mag0 = mag1 = 0.0;
  326.         lastsine = 0.0;
  327.         for (i = 0; i < hp->nt; i++) {
  328. #ifdef  DEBUG
  329.             if (dp->t != i || dp->p != j)
  330.                 error(CONSISTENCY,
  331.                     "division order in posgradient");
  332. #endif
  333.             b = bright(dp->v);
  334.             if (i > 0) {
  335.                 d = dp[-hp->np].r;
  336.                 if (dp[0].r > d) d = dp[0].r;
  337.                             /* sin(t)*cos(t)^2 */
  338.                 d *= lastsine * (1.0 - (double)i/hp->nt);
  339.                 mag0 += d*(b - bright(dp[-hp->np].v));
  340.             }
  341.             nextsine = sqrt((double)(i+1)/hp->nt);
  342.             if (j > 0) {
  343.                 d = dp[-1].r;
  344.                 if (dp[0].r > d) d = dp[0].r;
  345.                 mag1 += d * (nextsine - lastsine) *
  346.                         (b - bright(dp[-1].v));
  347.             } else {
  348.                 d = dp[hp->np-1].r;
  349.                 if (dp[0].r > d) d = dp[0].r;
  350.                 mag1 += d * (nextsine - lastsine) *
  351.                         (b - bright(dp[hp->np-1].v));
  352.             }
  353.             dp += hp->np;
  354.             lastsine = nextsine;
  355.         }
  356.         mag0 *= 2.0*PI / hp->np;
  357.         phi = 2.0*PI * (double)j/hp->np;
  358.         cosp = cos(phi); sinp = sin(phi);
  359.         xd += mag0*cosp - mag1*sinp;
  360.         yd += mag0*sinp + mag1*cosp;
  361.     }
  362.     for (i = 0; i < 3; i++)
  363.         gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])/PI;
  364. }
  365.  
  366.  
  367. dirgradient(gv, da, hp)                /* compute direction gradient */
  368. FVECT  gv;
  369. AMBSAMP  *da;            /* assumes standard ordering */
  370. register AMBHEMI  *hp;
  371. {
  372.     register int  i, j;
  373.     double  mag;
  374.     double  phi, xd, yd;
  375.     register AMBSAMP  *dp;
  376.  
  377.     xd = yd = 0.0;
  378.     for (j = 0; j < hp->np; j++) {
  379.         dp = da + j;
  380.         mag = 0.0;
  381.         for (i = 0; i < hp->nt; i++) {
  382. #ifdef  DEBUG
  383.             if (dp->t != i || dp->p != j)
  384.                 error(CONSISTENCY,
  385.                     "division order in dirgradient");
  386. #endif
  387.                             /* tan(t) */
  388.             mag += bright(dp->v)/sqrt(hp->nt/(i+.5) - 1.0);
  389.             dp += hp->np;
  390.         }
  391.         phi = 2.0*PI * (j+.5)/hp->np + PI/2.0;
  392.         xd += mag * cos(phi);
  393.         yd += mag * sin(phi);
  394.     }
  395.     for (i = 0; i < 3; i++)
  396.         gv[i] = (xd*hp->ux[i] + yd*hp->uy[i])/(hp->nt*hp->np);
  397. }
  398.