home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 8 / CDASC08.ISO / NEWS / RADIANCE / SRC / GEN / MKILLUM2.C < prev    next >
C/C++ Source or Header  |  1993-10-07  |  9KB  |  365 lines
  1. /* Copyright (c) 1991 Regents of the University of California */
  2.  
  3. #ifndef lint
  4. static char SCCSid[] = "@(#)mkillum2.c 2.5 8/13/92 LBL";
  5. #endif
  6.  
  7. /*
  8.  * Routines to do the actual calculation for mkillum
  9.  */
  10.  
  11. #include  "mkillum.h"
  12.  
  13. #include  "face.h"
  14.  
  15. #include  "cone.h"
  16.  
  17. #include  "random.h"
  18.  
  19.  
  20. o_default(ob, il, rt, nm)    /* default illum action */
  21. OBJREC  *ob;
  22. struct illum_args  *il;
  23. struct rtproc  *rt;
  24. char  *nm;
  25. {
  26.     sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
  27.             nm, ofun[ob->otype].funame, ob->oname);
  28.     error(WARNING, errmsg);
  29.     printobj(il->altmat, ob);
  30. }
  31.  
  32.  
  33. o_face(ob, il, rt, nm)        /* make an illum face */
  34. OBJREC  *ob;
  35. struct illum_args  *il;
  36. struct rtproc  *rt;
  37. char  *nm;
  38. {
  39. #define MAXMISS        (5*n*il->nsamps)
  40.     int  dim[3];
  41.     int  n, nalt, nazi, h;
  42.     float  *distarr;
  43.     double  sp[2], r1, r2;
  44.     FVECT  dn, org, dir;
  45.     FVECT  u, v;
  46.     double  ur[2], vr[2];
  47.     int  nmisses;
  48.     register FACE  *fa;
  49.     register int  i, j;
  50.                 /* get/check arguments */
  51.     fa = getface(ob);
  52.     if (fa->area == 0.0) {
  53.         freeface(ob);
  54.         o_default(ob, il, rt, nm);
  55.         return;
  56.     }
  57.                 /* set up sampling */
  58.     if (il->sampdens <= 0)
  59.         nalt = nazi = 1;
  60.     else {
  61.         n = PI * il->sampdens;
  62.         nalt = sqrt(n/PI) + .5;
  63.         nazi = PI*nalt + .5;
  64.     }
  65.     n = nalt*nazi;
  66.     distarr = (float *)calloc(n, 3*sizeof(float));
  67.     if (distarr == NULL)
  68.         error(SYSTEM, "out of memory in o_face");
  69.                 /* take first edge longer than sqrt(area) */
  70.     for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
  71.         u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
  72.         u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
  73.         u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
  74.         if ((r1 = DOT(u,u)) >= fa->area-FTINY)
  75.             break;
  76.     }
  77.     if (i < fa->nv) {    /* got one! -- let's align our axes */
  78.         r2 = 1.0/sqrt(r1);
  79.         u[0] *= r2; u[1] *= r2; u[2] *= r2;
  80.         fcross(v, fa->norm, u);
  81.     } else            /* oh well, we'll just have to wing it */
  82.         mkaxes(u, v, fa->norm);
  83.                 /* now, find limits in (u,v) coordinates */
  84.     ur[0] = vr[0] = FHUGE;
  85.     ur[1] = vr[1] = -FHUGE;
  86.     for (i = 0; i < fa->nv; i++) {
  87.         r1 = DOT(VERTEX(fa,i),u);
  88.         if (r1 < ur[0]) ur[0] = r1;
  89.         if (r1 > ur[1]) ur[1] = r1;
  90.         r2 = DOT(VERTEX(fa,i),v);
  91.         if (r2 < vr[0]) vr[0] = r2;
  92.         if (r2 > vr[1]) vr[1] = r2;
  93.     }
  94.     dim[0] = random();
  95.                 /* sample polygon */
  96.     nmisses = 0;
  97.     for (dim[1] = 0; dim[1] < nalt; dim[1]++)
  98.         for (dim[2] = 0; dim[2] < nazi; dim[2]++)
  99.         for (i = 0; i < il->nsamps; i++) {
  100.                     /* random direction */
  101.             h = ilhash(dim, 3) + i;
  102.             multisamp(sp, 2, urand(h));
  103.             r1 = (dim[1] + sp[0])/nalt;
  104.             r2 = (dim[2] + sp[1] - .5)/nazi;
  105.             flatdir(dn, r1, r2);
  106.             for (j = 0; j < 3; j++)
  107.             dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*fa->norm[j];
  108.                     /* random location */
  109.             do {
  110.             multisamp(sp, 2, urand(h+4862+nmisses));
  111.             r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
  112.             r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
  113.             for (j = 0; j < 3; j++)
  114.                 org[j] = r1*u[j] + r2*v[j]
  115.                     + fa->offset*fa->norm[j];
  116.             } while (!inface(org, fa) && nmisses++ < MAXMISS);
  117.             if (nmisses > MAXMISS) {
  118.             objerror(ob, WARNING, "bad aspect");
  119.             rt->nrays = 0;
  120.             freeface(ob);
  121.             free((char *)distarr);
  122.             o_default(ob, il, rt, nm);
  123.             return;
  124.             }
  125.             for (j = 0; j < 3; j++)
  126.             org[j] += .001*fa->norm[j];
  127.                     /* send sample */
  128.             raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
  129.         }
  130.     rayflush(rt);
  131.                 /* write out the face and its distribution */
  132.     if (average(il, distarr, nalt*nazi)) {
  133.         if (il->sampdens > 0)
  134.             flatout(il, distarr, nalt, nazi, u, v, fa->norm);
  135.         illumout(il, ob);
  136.     } else
  137.         printobj(il->altmat, ob);
  138.                 /* clean up */
  139.     freeface(ob);
  140.     free((char *)distarr);
  141. #undef MAXMISS
  142. }
  143.  
  144.  
  145. o_sphere(ob, il, rt, nm)    /* make an illum sphere */
  146. register OBJREC  *ob;
  147. struct illum_args  *il;
  148. struct rtproc  *rt;
  149. char  *nm;
  150. {
  151.     int  dim[3];
  152.     int  n, nalt, nazi;
  153.     float  *distarr;
  154.     double  sp[4], r1, r2, r3;
  155.     FVECT  org, dir;
  156.     FVECT  u, v;
  157.     register int  i, j;
  158.                 /* check arguments */
  159.     if (ob->oargs.nfargs != 4)
  160.         objerror(ob, USER, "bad # of arguments");
  161.                 /* set up sampling */
  162.     if (il->sampdens <= 0)
  163.         nalt = nazi = 1;
  164.     else {
  165.         n = 4.*PI * il->sampdens;
  166.         nalt = sqrt(n/PI) + .5;
  167.         nazi = PI*nalt + .5;
  168.     }
  169.     n = nalt*nazi;
  170.     distarr = (float *)calloc(n, 3*sizeof(float));
  171.     if (distarr == NULL)
  172.         error(SYSTEM, "out of memory in o_sphere");
  173.     dim[0] = random();
  174.                 /* sample sphere */
  175.     for (dim[1] = 0; dim[1] < nalt; dim[1]++)
  176.         for (dim[2] = 0; dim[2] < nazi; dim[2]++)
  177.         for (i = 0; i < il->nsamps; i++) {
  178.                     /* next sample point */
  179.             multisamp(sp, 4, urand(ilhash(dim,3)+i));
  180.                     /* random direction */
  181.             r1 = (dim[1] + sp[0])/nalt;
  182.             r2 = (dim[2] + sp[1] - .5)/nazi;
  183.             rounddir(dir, r1, r2);
  184.                     /* random location */
  185.             mkaxes(u, v, dir);        /* yuck! */
  186.             r3 = sqrt(sp[2]);
  187.             r2 = 2.*PI*sp[3];
  188.             r1 = r3*ob->oargs.farg[3]*cos(r2);
  189.             r2 = r3*ob->oargs.farg[3]*sin(r2);
  190.             r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
  191.             for (j = 0; j < 3; j++) {
  192.             org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
  193.                     r3*dir[j];
  194.             dir[j] = -dir[j];
  195.             }
  196.                     /* send sample */
  197.             raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
  198.         }
  199.     rayflush(rt);
  200.                 /* write out the sphere and its distribution */
  201.     if (average(il, distarr, nalt*nazi)) {
  202.         if (il->sampdens > 0)
  203.             roundout(il, distarr, nalt, nazi);
  204.         else
  205.             objerror(ob, WARNING, "diffuse distribution");
  206.         illumout(il, ob);
  207.     } else
  208.         printobj(il->altmat, ob);
  209.                 /* clean up */
  210.     free((char *)distarr);
  211. }
  212.  
  213.  
  214. o_ring(ob, il, rt, nm)        /* make an illum ring */
  215. OBJREC  *ob;
  216. struct illum_args  *il;
  217. struct rtproc  *rt;
  218. char  *nm;
  219. {
  220.     int  dim[3];
  221.     int  n, nalt, nazi;
  222.     float  *distarr;
  223.     double  sp[4], r1, r2, r3;
  224.     FVECT  dn, org, dir;
  225.     FVECT  u, v;
  226.     register CONE  *co;
  227.     register int  i, j;
  228.                 /* get/check arguments */
  229.     co = getcone(ob, 0);
  230.                 /* set up sampling */
  231.     if (il->sampdens <= 0)
  232.         nalt = nazi = 1;
  233.     else {
  234.         n = PI * il->sampdens;
  235.         nalt = sqrt(n/PI) + .5;
  236.         nazi = PI*nalt + .5;
  237.     }
  238.     n = nalt*nazi;
  239.     distarr = (float *)calloc(n, 3*sizeof(float));
  240.     if (distarr == NULL)
  241.         error(SYSTEM, "out of memory in o_ring");
  242.     mkaxes(u, v, co->ad);
  243.     dim[0] = random();
  244.                 /* sample disk */
  245.     for (dim[1] = 0; dim[1] < nalt; dim[1]++)
  246.         for (dim[2] = 0; dim[2] < nazi; dim[2]++)
  247.         for (i = 0; i < il->nsamps; i++) {
  248.                     /* next sample point */
  249.             multisamp(sp, 4, urand(ilhash(dim,3)+i));
  250.                     /* random direction */
  251.             r1 = (dim[1] + sp[0])/nalt;
  252.             r2 = (dim[2] + sp[1] - .5)/nazi;
  253.             flatdir(dn, r1, r2);
  254.             for (j = 0; j < 3; j++)
  255.             dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
  256.                     /* random location */
  257.             r3 = sqrt(CO_R0(co)*CO_R0(co) +
  258.                 sp[2]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
  259.             r2 = 2.*PI*sp[3];
  260.             r1 = r3*cos(r2);
  261.             r2 = r3*sin(r2);
  262.             for (j = 0; j < 3; j++)
  263.             org[j] = CO_P0(co)[j] + r1*u[j] + r1*v[j] +
  264.                     .001*co->ad[j];
  265.  
  266.                     /* send sample */
  267.             raysamp(distarr+3*(dim[1]*nazi+dim[2]), org, dir, rt);
  268.         }
  269.     rayflush(rt);
  270.                 /* write out the ring and its distribution */
  271.     if (average(il, distarr, nalt*nazi)) {
  272.         if (il->sampdens > 0)
  273.             flatout(il, distarr, nalt, nazi, u, v, co->ad);
  274.         illumout(il, ob);
  275.     } else
  276.         printobj(il->altmat, ob);
  277.                 /* clean up */
  278.     freecone(ob);
  279.     free((char *)distarr);
  280. }
  281.  
  282.  
  283. raysamp(res, org, dir, rt)    /* compute a ray sample */
  284. float  res[3];
  285. FVECT  org, dir;
  286. register struct rtproc  *rt;
  287. {
  288.     register float  *fp;
  289.  
  290.     if (rt->nrays == rt->bsiz)
  291.         rayflush(rt);
  292.     rt->dest[rt->nrays] = res;
  293.     fp = rt->buf + 6*rt->nrays++;
  294.     *fp++ = org[0]; *fp++ = org[1]; *fp++ = org[2];
  295.     *fp++ = dir[0]; *fp++ = dir[1]; *fp = dir[2];
  296. }
  297.  
  298.  
  299. rayflush(rt)            /* flush buffered rays */
  300. register struct rtproc  *rt;
  301. {
  302.     register int  i;
  303.  
  304.     if (rt->nrays <= 0)
  305.         return;
  306.     bzero(rt->buf+6*rt->nrays, 6*sizeof(float));
  307.     errno = 0;
  308.     if ( process(rt->pd, (char *)rt->buf, (char *)rt->buf,
  309.             3*sizeof(float)*rt->nrays,
  310.             6*sizeof(float)*(rt->nrays+1)) <
  311.             3*sizeof(float)*rt->nrays )
  312.         error(SYSTEM, "error reading from rtrace process");
  313.     i = rt->nrays;
  314.     while (i--) {
  315.         rt->dest[i][0] += rt->buf[3*i];
  316.         rt->dest[i][1] += rt->buf[3*i+1];
  317.         rt->dest[i][2] += rt->buf[3*i+2];
  318.     }
  319.     rt->nrays = 0;
  320. }
  321.  
  322.  
  323. mkaxes(u, v, n)            /* compute u and v to go with n */
  324. FVECT  u, v, n;
  325. {
  326.     register int  i;
  327.  
  328.     v[0] = v[1] = v[2] = 0.0;
  329.     for (i = 0; i < 3; i++)
  330.         if (n[i] < 0.6 && n[i] > -0.6)
  331.             break;
  332.     v[i] = 1.0;
  333.     fcross(u, v, n);
  334.     normalize(u);
  335.     fcross(v, n, u);
  336. }
  337.  
  338.  
  339. rounddir(dv, alt, azi)        /* compute uniform spherical direction */
  340. register FVECT  dv;
  341. double  alt, azi;
  342. {
  343.     double  d1, d2;
  344.  
  345.     dv[2] = 1. - 2.*alt;
  346.     d1 = sqrt(1. - dv[2]*dv[2]);
  347.     d2 = 2.*PI * azi;
  348.     dv[0] = d1*cos(d2);
  349.     dv[1] = d1*sin(d2);
  350. }
  351.  
  352.  
  353. flatdir(dv, alt, azi)        /* compute uniform hemispherical direction */
  354. register FVECT  dv;
  355. double  alt, azi;
  356. {
  357.     double  d1, d2;
  358.  
  359.     d1 = sqrt(alt);
  360.     d2 = 2.*PI * azi;
  361.     dv[0] = d1*cos(d2);
  362.     dv[1] = d1*sin(d2);
  363.     dv[2] = sqrt(1. - alt);
  364. }
  365.