Back to index

radiance  4R0+20100331
mkillum2.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: mkillum2.c,v 2.34 2009/09/09 15:32:20 greg Exp $";
00003 #endif
00004 /*
00005  * Routines to do the actual calculation for mkillum
00006  */
00007 
00008 #include <string.h>
00009 
00010 #include  "mkillum.h"
00011 #include  "face.h"
00012 #include  "cone.h"
00013 #include  "source.h"
00014 
00015 #ifndef NBSDFSAMPS
00016 #define NBSDFSAMPS   256           /* BSDF resampling count */
00017 #endif
00018 
00019 COLORV *      distarr = NULL;             /* distribution array */
00020 int           distsiz = 0;
00021 COLORV *      direct_discount = NULL;     /* amount to take off direct */
00022 
00023 
00024 void
00025 newdist(                    /* allocate & clear distribution array */
00026        int siz
00027 )
00028 {
00029        if (siz <= 0) {
00030               if (distsiz > 0)
00031                      free((void *)distarr);
00032               distarr = NULL;
00033               distsiz = 0;
00034               return;
00035        }
00036        if (distsiz < siz) {
00037               if (distsiz > 0)
00038                      free((void *)distarr);
00039               distarr = (COLORV *)malloc(sizeof(COLOR)*siz);
00040               if (distarr == NULL)
00041                      error(SYSTEM, "out of memory in newdist");
00042               distsiz = siz;
00043        }
00044        memset(distarr, '\0', sizeof(COLOR)*siz);
00045 }
00046 
00047 
00048 static void
00049 new_discount()                     /* allocate space for direct contrib. record */
00050 {
00051        if (distsiz <= 0)
00052               return;
00053        direct_discount = (COLORV *)calloc(distsiz, sizeof(COLOR));
00054        if (direct_discount == NULL)
00055               error(SYSTEM, "out of memory in new_discount");
00056 }
00057 
00058 
00059 static void
00060 done_discount()                    /* clear off direct contrib. record */
00061 {
00062        if (direct_discount == NULL)
00063               return;
00064        free((void *)direct_discount);
00065        direct_discount = NULL;
00066 }
00067 
00068 
00069 int
00070 process_ray(                /* process a ray result or report error */
00071        RAY *r,
00072        int rv
00073 )
00074 {
00075        COLORV *colp;
00076 
00077        if (rv == 0)                /* no result ready */
00078               return(0);
00079        if (rv < 0)
00080               error(USER, "ray tracing process died");
00081        if (r->rno >= distsiz)
00082               error(INTERNAL, "bad returned index in process_ray");
00083        multcolor(r->rcol, r->rcoef);      /* in case it's a source ray */
00084        colp = &distarr[r->rno * 3];
00085        addcolor(colp, r->rcol);
00086        if (r->rsrc >= 0 &&         /* remember source contrib. */
00087                      direct_discount != NULL) {
00088               colp = &direct_discount[r->rno * 3];
00089               addcolor(colp, r->rcol);
00090        }
00091        return(1);
00092 }
00093 
00094 
00095 void
00096 raysamp(                    /* queue a ray sample */
00097        int  ndx,
00098        FVECT  org,
00099        FVECT  dir
00100 )
00101 {
00102        RAY    myRay;
00103        int    rv;
00104 
00105        if ((ndx < 0) | (ndx >= distsiz))
00106               error(INTERNAL, "bad index in raysamp");
00107        VCOPY(myRay.rorg, org);
00108        VCOPY(myRay.rdir, dir);
00109        myRay.rmax = .0;
00110        rayorigin(&myRay, PRIMARY, NULL, NULL);
00111        myRay.rno = ndx;
00112                                    /* queue ray, check result */
00113        process_ray(&myRay, ray_pqueue(&myRay));
00114 }
00115 
00116 
00117 void
00118 srcsamps(                   /* sample sources from this surface position */
00119        struct illum_args *il,
00120        FVECT org,
00121        FVECT nrm,
00122        MAT4 ixfm
00123 )
00124 {
00125        int  nalt, nazi;
00126        SRCINDEX  si;
00127        RAY  sr;
00128        FVECT  v;
00129        double  d;
00130        int  i, j;
00131                                           /* get sampling density */
00132        if (il->sampdens <= 0) {
00133               nalt = nazi = 1;
00134        } else {
00135               i = PI * il->sampdens;
00136               nalt = sqrt(i/PI) + .5;
00137               nazi = PI*nalt + .5;
00138        }
00139        initsrcindex(&si);                 /* loop over (sub)sources */
00140        for ( ; ; ) {
00141               VCOPY(sr.rorg, org);        /* pick side to shoot from */
00142               if (il->sd != NULL) {
00143                      int  sn = si.sn;
00144                      if (si.sp+1 >= si.np) ++sn;
00145                      if (sn >= nsources) break;
00146                      if (source[sn].sflags & SDISTANT)
00147                             d = DOT(source[sn].sloc, nrm);
00148                      else {
00149                             VSUB(v, source[sn].sloc, org);
00150                             d = DOT(v, nrm);
00151                      }
00152               } else
00153                      d = 1.0;             /* only transmission */
00154               if (d < 0.0)
00155                      d = -1.0001*il->thick - 5.*FTINY;
00156               else
00157                      d = 5.*FTINY;
00158               for (i = 3; i--; )
00159                      sr.rorg[i] += d*nrm[i];
00160               samplendx++;                /* increment sample counter */
00161               if (!srcray(&sr, NULL, &si))
00162                      break;               /* end of sources */
00163                                           /* index direction */
00164               if (ixfm != NULL)
00165                      multv3(v, sr.rdir, ixfm);
00166               else
00167                      VCOPY(v, sr.rdir);
00168               if (il->sd != NULL) {
00169                      i = getBSDF_incndx(il->sd, v);
00170                      if (i < 0)
00171                             continue;     /* must not be important */
00172                      sr.rno = i;
00173                      d = 1.0/getBSDF_incohm(il->sd, i);
00174               } else {
00175                      if (v[2] >= -FTINY)
00176                             continue;     /* only sample transmission */
00177                      v[0] = -v[0]; v[1] = -v[1]; v[2] = -v[2];
00178                      sr.rno = flatindex(v, nalt, nazi);
00179                      d = nalt*nazi*(1./PI) * v[2];
00180               }
00181               d *= si.dom;                /* solid angle correction */
00182               scalecolor(sr.rcoef, d);
00183               process_ray(&sr, ray_pqueue(&sr));
00184        }
00185 }
00186 
00187 
00188 void
00189 rayclean()                  /* finish all pending rays */
00190 {
00191        RAY    myRay;
00192 
00193        while (process_ray(&myRay, ray_presult(&myRay, 0)))
00194               ;
00195 }
00196 
00197 
00198 static void
00199 mkaxes(                     /* compute u and v to go with n */
00200        FVECT  u,
00201        FVECT  v,
00202        FVECT  n
00203 )
00204 {
00205        register int  i;
00206 
00207        v[0] = v[1] = v[2] = 0.0;
00208        for (i = 0; i < 3; i++)
00209               if (n[i] < 0.6 && n[i] > -0.6)
00210                      break;
00211        v[i] = 1.0;
00212        fcross(u, v, n);
00213        normalize(u);
00214        fcross(v, n, u);
00215 }
00216 
00217 
00218 static void
00219 rounddir(            /* compute uniform spherical direction */
00220        register FVECT  dv,
00221        double  alt,
00222        double  azi
00223 )
00224 {
00225        double  d1, d2;
00226 
00227        dv[2] = 1. - 2.*alt;
00228        d1 = sqrt(1. - dv[2]*dv[2]);
00229        d2 = 2.*PI * azi;
00230        dv[0] = d1*cos(d2);
00231        dv[1] = d1*sin(d2);
00232 }
00233 
00234 
00235 void
00236 flatdir(             /* compute uniform hemispherical direction */
00237        FVECT  dv,
00238        double  alt,
00239        double  azi
00240 )
00241 {
00242        double  d1, d2;
00243 
00244        d1 = sqrt(alt);
00245        d2 = 2.*PI * azi;
00246        dv[0] = d1*cos(d2);
00247        dv[1] = d1*sin(d2);
00248        dv[2] = sqrt(1. - alt);
00249 }
00250 
00251 int
00252 flatindex(           /* compute index for hemispherical direction */
00253        FVECT  dv,
00254        int    nalt,
00255        int    nazi
00256 )
00257 {
00258        double d;
00259        int    i, j;
00260        
00261        d = 1.0 - dv[2]*dv[2];
00262        i = d*nalt;
00263        d = atan2(dv[1], dv[0]) * (0.5/PI);
00264        if (d < 0.0) d += 1.0;
00265        j = d*nazi + 0.5;
00266        if (j >= nazi) j = 0;
00267        return(i*nazi + j);
00268 }
00269 
00270 
00271 int
00272 my_default(   /* default illum action */
00273        OBJREC  *ob,
00274        struct illum_args  *il,
00275        char  *nm
00276 )
00277 {
00278        sprintf(errmsg, "(%s): cannot make illum for %s \"%s\"",
00279                      nm, ofun[ob->otype].funame, ob->oname);
00280        error(WARNING, errmsg);
00281        printobj(il->altmat, ob);
00282        return(1);
00283 }
00284 
00285 
00286 int
00287 my_face(             /* make an illum face */
00288        OBJREC  *ob,
00289        struct illum_args  *il,
00290        char  *nm
00291 )
00292 {
00293        int  dim[2];
00294        int  n, nalt, nazi, alti;
00295        double  sp[2], r1, r2;
00296        int  h;
00297        FVECT  dn, org, dir;
00298        FVECT  u, v;
00299        double  ur[2], vr[2];
00300        MAT4  xfm;
00301        int  nallow;
00302        FACE  *fa;
00303        int  i, j;
00304                             /* get/check arguments */
00305        fa = getface(ob);
00306        if (fa->area == 0.0) {
00307               freeface(ob);
00308               return(my_default(ob, il, nm));
00309        }
00310                             /* set up sampling */
00311        if (il->sd != NULL) {
00312               if (!getBSDF_xfm(xfm, fa->norm, il->udir)) {
00313                      objerror(ob, WARNING, "illegal up direction");
00314                      freeface(ob);
00315                      return(my_default(ob, il, nm));
00316               }
00317               n = il->sd->ninc;
00318        } else {
00319               if (il->sampdens <= 0) {
00320                      nalt = nazi = 1;     /* diffuse assumption */
00321               } else {
00322                      n = PI * il->sampdens;
00323                      nalt = sqrt(n/PI) + .5;
00324                      nazi = PI*nalt + .5;
00325               }
00326               n = nazi*nalt;
00327        }
00328        newdist(n);
00329                             /* take first edge >= sqrt(area) */
00330        for (j = fa->nv-1, i = 0; i < fa->nv; j = i++) {
00331               u[0] = VERTEX(fa,i)[0] - VERTEX(fa,j)[0];
00332               u[1] = VERTEX(fa,i)[1] - VERTEX(fa,j)[1];
00333               u[2] = VERTEX(fa,i)[2] - VERTEX(fa,j)[2];
00334               if ((r1 = DOT(u,u)) >= fa->area-FTINY)
00335                      break;
00336        }
00337        if (i < fa->nv) {    /* got one! -- let's align our axes */
00338               r2 = 1.0/sqrt(r1);
00339               u[0] *= r2; u[1] *= r2; u[2] *= r2;
00340               fcross(v, fa->norm, u);
00341        } else               /* oh well, we'll just have to wing it */
00342               mkaxes(u, v, fa->norm);
00343                             /* now, find limits in (u,v) coordinates */
00344        ur[0] = vr[0] = FHUGE;
00345        ur[1] = vr[1] = -FHUGE;
00346        for (i = 0; i < fa->nv; i++) {
00347               r1 = DOT(VERTEX(fa,i),u);
00348               if (r1 < ur[0]) ur[0] = r1;
00349               if (r1 > ur[1]) ur[1] = r1;
00350               r2 = DOT(VERTEX(fa,i),v);
00351               if (r2 < vr[0]) vr[0] = r2;
00352               if (r2 > vr[1]) vr[1] = r2;
00353        }
00354        dim[0] = random();
00355                             /* sample polygon */
00356        nallow = 5*n*il->nsamps;
00357        for (dim[1] = 0; dim[1] < n; dim[1]++)
00358               for (i = 0; i < il->nsamps; i++) {
00359                                    /* randomize direction */
00360                   h = ilhash(dim, 2) + i;
00361                   if (il->sd != NULL) {
00362                      r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
00363                   } else {
00364                      multisamp(sp, 2, urand(h));
00365                      alti = dim[1]/nazi;
00366                      r1 = (alti + sp[0])/nalt;
00367                      r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
00368                      flatdir(dn, r1, r2);
00369                      for (j = 0; j < 3; j++)
00370                          dir[j] = -dn[0]*u[j] - dn[1]*v[j] -
00371                                           dn[2]*fa->norm[j];
00372                   }
00373                                    /* randomize location */
00374                   do {
00375                      multisamp(sp, 2, urand(h+4862+nallow));
00376                      r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
00377                      r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
00378                      for (j = 0; j < 3; j++)
00379                          org[j] = r1*u[j] + r2*v[j]
00380                                    + fa->offset*fa->norm[j];
00381                   } while (!inface(org, fa) && nallow-- > 0);
00382                   if (nallow < 0) {
00383                      objerror(ob, WARNING, "bad aspect");
00384                      rayclean();
00385                      freeface(ob);
00386                      return(my_default(ob, il, nm));
00387                   }
00388                   if (il->sd != NULL && DOT(dir, fa->norm) < -FTINY)
00389                      r1 = -1.0001*il->thick - 5.*FTINY;
00390                   else
00391                      r1 = 5.*FTINY;
00392                   for (j = 0; j < 3; j++)
00393                      org[j] += r1*fa->norm[j];
00394                                    /* send sample */
00395                   raysamp(dim[1], org, dir);
00396               }
00397                             /* add in direct component? */
00398        if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
00399               MAT4   ixfm;
00400               if (il->sd == NULL) {
00401                      for (i = 3; i--; ) {
00402                             ixfm[i][0] = u[i];
00403                             ixfm[i][1] = v[i];
00404                             ixfm[i][2] = fa->norm[i];
00405                             ixfm[i][3] = 0.;
00406                      }
00407                      ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
00408                      ixfm[3][3] = 1.;
00409               } else {
00410                      if (!invmat4(ixfm, xfm))
00411                             objerror(ob, INTERNAL,
00412                                    "cannot invert BSDF transform");
00413                      if (!(il->flags & IL_LIGHT))
00414                             new_discount();
00415               }
00416               dim[0] = random();
00417               nallow = 10*il->nsamps;
00418               for (i = 0; i < il->nsamps; i++) {
00419                                    /* randomize location */
00420                   h = dim[0] + samplendx++;
00421                   do {
00422                      multisamp(sp, 2, urand(h+nallow));
00423                      r1 = ur[0] + (ur[1]-ur[0]) * sp[0];
00424                      r2 = vr[0] + (vr[1]-vr[0]) * sp[1];
00425                      for (j = 0; j < 3; j++)
00426                          org[j] = r1*u[j] + r2*v[j]
00427                                    + fa->offset*fa->norm[j];
00428                   } while (!inface(org, fa) && nallow-- > 0);
00429                   if (nallow < 0) {
00430                      objerror(ob, WARNING, "bad aspect");
00431                      rayclean();
00432                      freeface(ob);
00433                      return(my_default(ob, il, nm));
00434                   }
00435                                    /* sample source rays */
00436                   srcsamps(il, org, fa->norm, ixfm);
00437               }
00438        }
00439                             /* wait for all rays to finish */
00440        rayclean();
00441        if (il->sd != NULL) {       /* run distribution through BSDF */
00442               nalt = sqrt(il->sd->nout/PI) + .5;
00443               nazi = PI*nalt + .5;
00444               redistribute(il->sd, nalt, nazi, u, v, fa->norm, xfm);
00445               done_discount();
00446               if (!il->sampdens)
00447                      il->sampdens = nalt*nazi/PI + .999;
00448        }
00449                             /* write out the face and its distribution */
00450        if (average(il, distarr, n)) {
00451               if (il->sampdens > 0)
00452                      flatout(il, distarr, nalt, nazi, u, v, fa->norm);
00453               illumout(il, ob);
00454        } else
00455               printobj(il->altmat, ob);
00456                             /* clean up */
00457        freeface(ob);
00458        return(0);
00459 }
00460 
00461 
00462 int
00463 my_sphere(    /* make an illum sphere */
00464        register OBJREC  *ob,
00465        struct illum_args  *il,
00466        char  *nm
00467 )
00468 {
00469        int  dim[3];
00470        int  n, nalt, nazi;
00471        double  sp[4], r1, r2, r3;
00472        FVECT  org, dir;
00473        FVECT  u, v;
00474        register int  i, j;
00475                             /* check arguments */
00476        if (ob->oargs.nfargs != 4)
00477               objerror(ob, USER, "bad # of arguments");
00478                             /* set up sampling */
00479        if (il->sampdens <= 0)
00480               nalt = nazi = 1;
00481        else {
00482               n = 4.*PI * il->sampdens;
00483               nalt = sqrt(2./PI*n) + .5;
00484               nazi = PI/2.*nalt + .5;
00485        }
00486        if (il->sd != NULL)
00487               objerror(ob, WARNING, "BSDF ignored");
00488        n = nalt*nazi;
00489        newdist(n);
00490        dim[0] = random();
00491                             /* sample sphere */
00492        for (dim[1] = 0; dim[1] < nalt; dim[1]++)
00493            for (dim[2] = 0; dim[2] < nazi; dim[2]++)
00494               for (i = 0; i < il->nsamps; i++) {
00495                                    /* next sample point */
00496                   multisamp(sp, 4, urand(ilhash(dim,3)+i));
00497                                    /* random direction */
00498                   r1 = (dim[1] + sp[0])/nalt;
00499                   r2 = (dim[2] + sp[1] - .5)/nazi;
00500                   rounddir(dir, r1, r2);
00501                                    /* random location */
00502                   mkaxes(u, v, dir);             /* yuck! */
00503                   r3 = sqrt(sp[2]);
00504                   r2 = 2.*PI*sp[3];
00505                   r1 = r3*ob->oargs.farg[3]*cos(r2);
00506                   r2 = r3*ob->oargs.farg[3]*sin(r2);
00507                   r3 = ob->oargs.farg[3]*sqrt(1.01-r3*r3);
00508                   for (j = 0; j < 3; j++) {
00509                      org[j] = ob->oargs.farg[j] + r1*u[j] + r2*v[j] +
00510                                    r3*dir[j];
00511                      dir[j] = -dir[j];
00512                   }
00513                                    /* send sample */
00514                   raysamp(dim[1]*nazi+dim[2], org, dir);
00515               }
00516                             /* wait for all rays to finish */
00517        rayclean();
00518                             /* write out the sphere and its distribution */
00519        if (average(il, distarr, n)) {
00520               if (il->sampdens > 0)
00521                      roundout(il, distarr, nalt, nazi);
00522               else
00523                      objerror(ob, WARNING, "diffuse distribution");
00524               illumout(il, ob);
00525        } else
00526               printobj(il->altmat, ob);
00527                             /* clean up */
00528        return(1);
00529 }
00530 
00531 
00532 int
00533 my_ring(             /* make an illum ring */
00534        OBJREC  *ob,
00535        struct illum_args  *il,
00536        char  *nm
00537 )
00538 {
00539        int  dim[2];
00540        int  n, nalt, nazi, alti;
00541        double  sp[2], r1, r2, r3;
00542        int  h;
00543        FVECT  dn, org, dir;
00544        FVECT  u, v;
00545        MAT4  xfm;
00546        CONE  *co;
00547        int  i, j;
00548                             /* get/check arguments */
00549        co = getcone(ob, 0);
00550                             /* set up sampling */
00551        if (il->sd != NULL) {
00552               if (!getBSDF_xfm(xfm, co->ad, il->udir)) {
00553                      objerror(ob, WARNING, "illegal up direction");
00554                      freecone(ob);
00555                      return(my_default(ob, il, nm));
00556               }
00557               n = il->sd->ninc;
00558        } else {
00559               if (il->sampdens <= 0) {
00560                      nalt = nazi = 1;     /* diffuse assumption */
00561               } else {
00562                      n = PI * il->sampdens;
00563                      nalt = sqrt(n/PI) + .5;
00564                      nazi = PI*nalt + .5;
00565               }
00566               n = nazi*nalt;
00567        }
00568        newdist(n);
00569        mkaxes(u, v, co->ad);
00570        dim[0] = random();
00571                             /* sample disk */
00572        for (dim[1] = 0; dim[1] < n; dim[1]++)
00573               for (i = 0; i < il->nsamps; i++) {
00574                                    /* next sample point */
00575                   h = ilhash(dim,2) + i;
00576                                    /* randomize direction */
00577                   if (il->sd != NULL) {
00578                      r_BSDF_incvec(dir, il->sd, dim[1], urand(h), xfm);
00579                   } else {
00580                      multisamp(sp, 2, urand(h));
00581                      alti = dim[1]/nazi;
00582                      r1 = (alti + sp[0])/nalt;
00583                      r2 = (dim[1] - alti*nazi + sp[1] - .5)/nazi;
00584                      flatdir(dn, r1, r2);
00585                      for (j = 0; j < 3; j++)
00586                             dir[j] = -dn[0]*u[j] - dn[1]*v[j] - dn[2]*co->ad[j];
00587                   }
00588                                    /* randomize location */
00589                   multisamp(sp, 2, urand(h+8371));
00590                   r3 = sqrt(CO_R0(co)*CO_R0(co) +
00591                          sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
00592                   r2 = 2.*PI*sp[1];
00593                   r1 = r3*cos(r2);
00594                   r2 = r3*sin(r2);
00595                   if (il->sd != NULL && DOT(dir, co->ad) < -FTINY)
00596                      r3 = -1.0001*il->thick - 5.*FTINY;
00597                   else
00598                      r3 = 5.*FTINY;
00599                   for (j = 0; j < 3; j++)
00600                      org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j] +
00601                                           r3*co->ad[j];
00602                                    /* send sample */
00603                   raysamp(dim[1], org, dir);
00604               }
00605                             /* add in direct component? */
00606        if (!directvis && (il->flags & IL_LIGHT || il->sd != NULL)) {
00607               MAT4   ixfm;
00608               if (il->sd == NULL) {
00609                      for (i = 3; i--; ) {
00610                             ixfm[i][0] = u[i];
00611                             ixfm[i][1] = v[i];
00612                             ixfm[i][2] = co->ad[i];
00613                             ixfm[i][3] = 0.;
00614                      }
00615                      ixfm[3][0] = ixfm[3][1] = ixfm[3][2] = 0.;
00616                      ixfm[3][3] = 1.;
00617               } else {
00618                      if (!invmat4(ixfm, xfm))
00619                             objerror(ob, INTERNAL,
00620                                    "cannot invert BSDF transform");
00621                      if (!(il->flags & IL_LIGHT))
00622                             new_discount();
00623               }
00624               dim[0] = random();
00625               for (i = 0; i < il->nsamps; i++) {
00626                                    /* randomize location */
00627                   h = dim[0] + samplendx++;
00628                   multisamp(sp, 2, urand(h));
00629                   r3 = sqrt(CO_R0(co)*CO_R0(co) +
00630                          sp[0]*(CO_R1(co)*CO_R1(co) - CO_R0(co)*CO_R0(co)));
00631                   r2 = 2.*PI*sp[1];
00632                   r1 = r3*cos(r2);
00633                   r2 = r3*sin(r2);
00634                   for (j = 0; j < 3; j++)
00635                      org[j] = CO_P0(co)[j] + r1*u[j] + r2*v[j];
00636                                    /* sample source rays */
00637                   srcsamps(il, org, co->ad, ixfm);
00638               }
00639        }
00640                             /* wait for all rays to finish */
00641        rayclean();
00642        if (il->sd != NULL) {       /* run distribution through BSDF */
00643               nalt = sqrt(il->sd->nout/PI) + .5;
00644               nazi = PI*nalt + .5;
00645               redistribute(il->sd, nalt, nazi, u, v, co->ad, xfm);
00646               done_discount();
00647               if (!il->sampdens)
00648                      il->sampdens = nalt*nazi/PI + .999;
00649        }
00650                             /* write out the ring and its distribution */
00651        if (average(il, distarr, n)) {
00652               if (il->sampdens > 0)
00653                      flatout(il, distarr, nalt, nazi, u, v, co->ad);
00654               illumout(il, ob);
00655        } else
00656               printobj(il->altmat, ob);
00657                             /* clean up */
00658        freecone(ob);
00659        return(1);
00660 }
00661 
00662 
00663 void
00664 redistribute(        /* pass distarr ray sums through BSDF */
00665        struct BSDF_data *b,
00666        int nalt,
00667        int nazi,
00668        FVECT u,
00669        FVECT v,
00670        FVECT w,
00671        MAT4 xm
00672 )
00673 {
00674        int    nout = 0;
00675        MAT4   mymat, inmat;
00676        COLORV *idist;
00677        COLORV *cp;
00678        FVECT  dv;
00679        double wt;
00680        int    i, j, k, c, o;
00681        COLOR  col, cinc;
00682                                    /* copy incoming distribution */
00683        if (b->ninc > distsiz)
00684               error(INTERNAL, "error 1 in redistribute");
00685        idist = (COLORV *)malloc(sizeof(COLOR)*b->ninc);
00686        if (idist == NULL)
00687               error(SYSTEM, "out of memory in redistribute");
00688        memcpy(idist, distarr, sizeof(COLOR)*b->ninc);
00689                                    /* compose direction transform */
00690        for (i = 3; i--; ) {
00691               mymat[i][0] = u[i];
00692               mymat[i][1] = v[i];
00693               mymat[i][2] = w[i];
00694               mymat[i][3] = 0.;
00695        }
00696        mymat[3][0] = mymat[3][1] = mymat[3][2] = 0.;
00697        mymat[3][3] = 1.;
00698        if (xm != NULL)
00699               multmat4(mymat, xm, mymat);
00700        for (i = 3; i--; ) {        /* make sure it's normalized */
00701               wt = 1./sqrt( mymat[0][i]*mymat[0][i] +
00702                             mymat[1][i]*mymat[1][i] +
00703                             mymat[2][i]*mymat[2][i]     );
00704               for (j = 3; j--; )
00705                      mymat[j][i] *= wt;
00706        }
00707        if (!invmat4(inmat, mymat)) /* need inverse as well */
00708               error(INTERNAL, "cannot invert BSDF transform");
00709        newdist(nalt*nazi);         /* resample distribution */
00710        for (i = b->ninc; i--; ) {
00711               int    direct_out = -1;
00712               COLOR  cdir;
00713               getBSDF_incvec(dv, b, i);   /* compute incident irrad. */
00714               multv3(dv, dv, mymat);
00715               if (dv[2] < 0.0) {
00716                      dv[0] = -dv[0]; dv[1] = -dv[1]; dv[2] = -dv[2];
00717                      direct_out += (direct_discount != NULL);
00718               }
00719               wt = getBSDF_incohm(b, i);
00720               wt *= dv[2];                /* solid_angle*cosine(theta) */
00721               cp = &idist[3*i];
00722               copycolor(cinc, cp);
00723               scalecolor(cinc, wt);
00724               if (!direct_out) {          /* discount direct contr. */
00725                      cp = &direct_discount[3*i];
00726                      copycolor(cdir, cp);
00727                      scalecolor(cdir, -wt);
00728                      if (b->nout != b->ninc)
00729                             direct_out = flatindex(dv, nalt, nazi);
00730                      else
00731                             direct_out = i;      /* assumes dist. mirroring */
00732               }
00733               for (k = nalt; k--; )              /* loop over distribution */
00734                 for (j = nazi; j--; ) {
00735                   int       rstart = random();
00736                   for (c = NBSDFSAMPS; c--; ) {
00737                      double  sp[2];
00738                      multisamp(sp, 2, urand(rstart+c));
00739                      flatdir(dv, (k + sp[0])/nalt,
00740                                    (j + .5 - sp[1])/nazi);
00741                      multv3(dv, dv, inmat);
00742                                           /* evaluate BSDF @ outgoing */
00743                      o = getBSDF_outndx(b, dv);
00744                      if (o < 0) {
00745                             nout++;
00746                             continue;
00747                      }
00748                      wt = BSDF_value(b, i, o) * (1./NBSDFSAMPS);
00749                      copycolor(col, cinc);
00750                      if (b->nout != b->ninc)
00751                             o = k*nazi + j;
00752                      if (o == direct_out)
00753                             addcolor(col, cdir); /* minus direct */
00754                      scalecolor(col, wt);
00755                      cp = &distarr[3*(k*nazi + j)];
00756                      addcolor(cp, col);   /* sum into distribution */
00757                   }
00758                 }
00759        }
00760        free(idist);                /* free temp space */
00761        if (nout) {
00762               sprintf(errmsg, "missing %.1f%% of BSDF directions",
00763                             100.*nout/(b->ninc*nalt*nazi*NBSDFSAMPS));
00764               error(WARNING, errmsg);
00765        }
00766 }