Back to index

radiance  4R0+20100331
mksource.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char RCSid[] = "$Id: mksource.c,v 2.5 2007/08/24 04:41:36 greg Exp $";
00003 #endif
00004 /*
00005  * Generate distant sources corresponding to the given environment map
00006  */
00007 
00008 #include "ray.h"
00009 #include "random.h"
00010 #include "resolu.h"
00011 
00012 #define NTRUNKBR     4             /* number of branches at trunk */
00013 #define NTRUNKVERT   4             /* number of vertices at trunk */
00014 #define DEF_NSAMPS   262144L              /* default # sphere samples */
00015 #define DEF_MAXANG   15.           /* maximum source angle (deg.) */
00016 
00017 /* Data structure for geodesic samples */
00018 
00019 typedef struct tritree {
00020        int32         gdv[3];              /* spherical triangle vertex direc. */
00021        int32         sd;           /* sample direction if leaf */
00022        struct tritree       *kid;         /* 4 children if branch node */
00023        COLR          val;          /* sampled color value */
00024 } TRITREE;
00025 
00026 typedef struct lostlight {
00027        struct lostlight     *next; /* next in list */
00028        int32         sd;           /* lost source direction */
00029        COLOR         intens;              /* output times solid angle */
00030 } LOSTLIGHT;
00031 
00032 extern char   *progname;
00033 
00034 FVECT  scene_cent;          /* center of octree cube */
00035 RREAL  scene_rad;           /* radius to get outside cube from center */
00036 
00037 const COLR    blkclr = BLKCOLR;
00038 
00039 #define isleaf(node) ((node)->kid == NULL)
00040 
00041 /* Compute signum of signed volume for three vectors */
00042 int
00043 vol_sign(const FVECT v1, const FVECT v2, const FVECT v3)
00044 {
00045        double vol;
00046        
00047        vol  = v1[0]*(v2[1]*v3[2] - v2[2]*v3[1]);
00048        vol += v1[1]*(v2[2]*v3[0] - v2[0]*v3[2]);
00049        vol += v1[2]*(v2[0]*v3[1] - v2[1]*v3[0]);
00050        if (vol > .0)
00051               return(1);
00052        if (vol < .0)
00053               return(-1);
00054        return(0);
00055 }
00056 
00057 /* Is the given direction contained within the specified spherical triangle? */
00058 int
00059 intriv(const int32 trid[3], const FVECT sdir)
00060 {
00061        int    sv[3];
00062        FVECT  tri[3];
00063 
00064        decodedir(tri[0], trid[0]);
00065        decodedir(tri[1], trid[1]);
00066        decodedir(tri[2], trid[2]);
00067        sv[0] = vol_sign(sdir, tri[0], tri[1]);
00068        sv[1] = vol_sign(sdir, tri[1], tri[2]);
00069        sv[2] = vol_sign(sdir, tri[2], tri[0]);
00070        if ((sv[0] == sv[1]) & (sv[1] == sv[2]))
00071               return(1);
00072        return(!sv[0] | !sv[1] | !sv[2]);
00073 }
00074 
00075 /* Find leaf containing the given sample direction */
00076 TRITREE *
00077 findleaf(TRITREE *node, const FVECT sdir)
00078 {
00079        int    i;
00080 
00081        if (isleaf(node))
00082               return(intriv(node->gdv,sdir) ? node : (TRITREE *)NULL);
00083        for (i = 0; i < 4; i++) {
00084               TRITREE       *chknode = &node->kid[i];
00085               if (intriv(chknode->gdv, sdir))
00086                      return(isleaf(chknode) ? chknode :
00087                                    findleaf(chknode, sdir));
00088        }
00089        return(NULL);
00090 }
00091 
00092 /* Initialize leaf with random sample inside the given spherical triangle */
00093 void
00094 leafsample(TRITREE *leaf)
00095 {
00096        RAY    myray;
00097        RREAL  wt[3];
00098        FVECT  sdir;
00099        int    i, j;
00100                                           /* random point on triangle */
00101        i = random() % 3;
00102        wt[i] = frandom();
00103        j = random() & 1;
00104        wt[(i+2-j)%3] = 1. - wt[i] - 
00105                      (wt[(i+1+j)%3] = (1.-wt[i])*frandom());
00106        sdir[0] = sdir[1] = sdir[2] = .0;
00107        for (i = 0; i < 3; i++) {
00108               FVECT  vt;
00109               decodedir(vt, leaf->gdv[i]);
00110               VSUM(sdir, sdir, vt, wt[i]);
00111        }
00112        normalize(sdir);                   /* record sample direction */
00113        leaf->sd = encodedir(sdir);
00114                                           /* evaluate at inf. */
00115        VSUM(myray.rorg, scene_cent, sdir, scene_rad);
00116        VCOPY(myray.rdir, sdir);
00117        myray.rmax = 0.;
00118        ray_trace(&myray);
00119        setcolr(leaf->val, colval(myray.rcol,RED),
00120                      colval(myray.rcol,GRN),
00121                      colval(myray.rcol,BLU));
00122 }
00123 
00124 /* Initialize a branch node contained in the given spherical triangle */
00125 void
00126 subdivide(TRITREE *branch, const int32 dv[3])
00127 {
00128        FVECT  dvv[3], sdv[3];
00129        int32  sd[3];
00130        int    i;
00131        
00132        for (i = 0; i < 3; i++) {   /* copy spherical triangle */
00133               branch->gdv[i] = dv[i];
00134               decodedir(dvv[i], dv[i]);
00135        }
00136        for (i = 0; i < 3; i++) {   /* create new vertices */
00137               int    j = (i+1)%3;
00138               VADD(sdv[i], dvv[i], dvv[j]);
00139               normalize(sdv[i]);
00140               sd[i] = encodedir(sdv[i]);
00141        }
00142                                    /* allocate leaves */
00143        branch->kid = (TRITREE *)calloc(4, sizeof(TRITREE));
00144        if (branch->kid == NULL)
00145               error(SYSTEM, "out of memory in subdivide()");
00146                                    /* assign subtriangle directions */
00147        branch->kid[0].gdv[0] = dv[0];
00148        branch->kid[0].gdv[1] = sd[0];
00149        branch->kid[0].gdv[2] = sd[2];
00150        branch->kid[1].gdv[0] = sd[0];
00151        branch->kid[1].gdv[1] = dv[1];
00152        branch->kid[1].gdv[2] = sd[1];
00153        branch->kid[2].gdv[0] = sd[1];
00154        branch->kid[2].gdv[1] = dv[2];
00155        branch->kid[2].gdv[2] = sd[2];
00156        branch->kid[3].gdv[0] = sd[0];
00157        branch->kid[3].gdv[1] = sd[1];
00158        branch->kid[3].gdv[2] = sd[2];
00159 }
00160 
00161 /* Recursively subdivide the given node to the specified quadtree depth */
00162 void
00163 branchsample(TRITREE *node, int depth)
00164 {
00165        int    i;
00166 
00167        if (depth <= 0)
00168               return;
00169        if (isleaf(node)) {                /* subdivide leaf node */
00170               TRITREE       branch, *moved_leaf;
00171               FVECT  sdir;
00172               subdivide(&branch, node->gdv);
00173               decodedir(sdir, node->sd);
00174               moved_leaf = findleaf(&branch, sdir);
00175               if (moved_leaf != NULL) {   /* bequeath old sample */
00176                      moved_leaf->sd = node->sd;
00177                      copycolr(moved_leaf->val, node->val);
00178               }
00179               for (i = 0; i < 4; i++)            /* compute new samples */
00180                      if (&branch.kid[i] != moved_leaf)
00181                             leafsample(&branch.kid[i]);
00182               *node = branch;                    /* replace leaf with branch */
00183        }
00184        for (i = 0; i < 4; i++)                   /* subdivide children */
00185               branchsample(&node->kid[i], depth-1);
00186 }
00187 
00188 /* Sample sphere using triangular geodesic mesh */
00189 TRITREE       *
00190 geosample(int nsamps)
00191 {
00192        int    depth;
00193        TRITREE       *tree;
00194        FVECT  trunk[NTRUNKVERT];
00195        int    i, j;
00196                                    /* figure out depth */
00197        if ((nsamps -= 4) < 0)
00198               error(USER, "minimum number of samples is 4");
00199        nsamps = nsamps*3/NTRUNKBR; /* round up */
00200        for (depth = 0; nsamps > 1; depth++)
00201               nsamps >>= 2;
00202                                    /* make base tetrahedron */
00203        tree = (TRITREE *)malloc(sizeof(TRITREE));
00204        if (tree == NULL) goto memerr;
00205        trunk[0][0] = trunk[0][1] = 0; trunk[0][2] = 1;
00206        trunk[1][0] = 0;
00207        trunk[1][2] = cos(2.*asin(sqrt(2./3.)));
00208        trunk[1][1] = sqrt(1. - trunk[1][2]*trunk[1][2]);
00209        spinvector(trunk[2], trunk[1], trunk[0], 2.*PI/3.);
00210        spinvector(trunk[3], trunk[1], trunk[0], 4.*PI/3.);
00211        tree->gdv[0] = tree->gdv[1] = tree->gdv[2] = encodedir(trunk[0]);
00212        tree->kid = (TRITREE *)calloc(NTRUNKBR, sizeof(TRITREE));
00213        if (tree->kid == NULL) goto memerr;
00214                                    /* grow our tree from trunk */
00215        for (i = 0; i < NTRUNKBR; i++) {
00216               for (j = 0; j < 3; j++)            /* XXX works for tetra only */
00217                      tree->kid[i].gdv[j] = encodedir(trunk[(i+j)%NTRUNKVERT]);
00218               leafsample(&tree->kid[i]);
00219               branchsample(&tree->kid[i], depth);
00220        }
00221        return(tree);
00222 memerr:
00223        error(SYSTEM, "out of memory in geosample()");
00224        return NULL; /* dummy return */
00225 }
00226 
00227 /* Compute leaf exponent histogram */
00228 void
00229 get_ehisto(const TRITREE *node, long exphisto[256])
00230 {
00231        int    i;
00232 
00233        if (isleaf(node)) {
00234               ++exphisto[node->val[EXP]];
00235               return;
00236        }
00237        for (i = 0; i < 4; i++)
00238               get_ehisto(&node->kid[i], exphisto);
00239 }
00240 
00241 /* Get reasonable source threshold */
00242 double
00243 get_threshold(const TRITREE *tree)
00244 {
00245        long   exphisto[256];
00246        long   samptotal;
00247        int    i;
00248                                           /* compute sample histogram */
00249        memset((void *)exphisto, 0, sizeof(exphisto));
00250        for (i = 0; i < NTRUNKBR; i++)
00251               get_ehisto(&tree->kid[i], exphisto);
00252                                           /* use 98th percentile */
00253        for (i = 0; i < 256; i++)
00254               samptotal += exphisto[i];
00255        samptotal /= 50; 
00256        for (i = 256; (--i > 0) & (samptotal > 0); )
00257               samptotal -= exphisto[i];
00258        return(ldexp(.75, i-COLXS));
00259 }
00260 
00261 /* Find leaf containing the maximum exponent */
00262 TRITREE *
00263 findemax(TRITREE *node, int *expp)
00264 {
00265        if (!isleaf(node)) {
00266               TRITREE       *maxleaf;
00267               TRITREE *rleaf;
00268               maxleaf = findemax(&node->kid[0], expp);
00269               rleaf = findemax(&node->kid[1], expp);
00270               if (rleaf != NULL) maxleaf = rleaf;
00271               rleaf = findemax(&node->kid[2], expp);
00272               if (rleaf != NULL) maxleaf = rleaf;
00273               rleaf = findemax(&node->kid[3], expp);
00274               if (rleaf != NULL) maxleaf = rleaf;
00275               return(maxleaf);
00276        }
00277        if (node->val[EXP] <= *expp)
00278               return(NULL);
00279        *expp = node->val[EXP];
00280        return(node);
00281 }
00282 
00283 /* Compute solid angle of spherical triangle (approx.) */
00284 double
00285 tri_omegav(const int32 vd[3])
00286 {
00287        FVECT  v[3], e1, e2, vcross;
00288 
00289        decodedir(v[0], vd[0]);
00290        decodedir(v[1], vd[1]);
00291        decodedir(v[2], vd[2]);
00292        VSUB(e1, v[1], v[0]);
00293        VSUB(e2, v[2], v[1]);
00294        fcross(vcross, e1, e2);
00295        return(.5*VLEN(vcross));
00296 }
00297 
00298 /* Sum intensity times direction for above-threshold perimiter within radius */
00299 void
00300 vector_sum(FVECT vsum, TRITREE *node,
00301               FVECT cent, double maxr2, int ethresh)
00302 {
00303        if (isleaf(node)) {
00304               double intens;
00305               FVECT  sdir;
00306               if (node->val[EXP] < ethresh)
00307                      return;                            /* below threshold */
00308               if (fdir2diff(node->sd,cent) > maxr2)
00309                      return;                            /* too far away */
00310               intens = colrval(node->val,GRN) * tri_omegav(node->gdv);
00311               decodedir(sdir, node->sd);
00312               VSUM(vsum, vsum, sdir, intens);
00313               return;
00314        }
00315        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
00316                      fdir2diff(node->gdv[0],cent) < maxr2 &&
00317                      fdir2diff(node->gdv[1],cent) < maxr2 &&
00318                      fdir2diff(node->gdv[2],cent) < maxr2)
00319               return;                                   /* containing node */
00320        vector_sum(vsum, &node->kid[0], cent, maxr2, ethresh);
00321        vector_sum(vsum, &node->kid[1], cent, maxr2, ethresh);
00322        vector_sum(vsum, &node->kid[2], cent, maxr2, ethresh);
00323        vector_sum(vsum, &node->kid[3], cent, maxr2, ethresh);
00324 }
00325 
00326 /* Claim source contributions within the given solid angle */
00327 void
00328 claimlight(COLOR intens, TRITREE *node, FVECT cent, double maxr2)
00329 {
00330        int    remaining;
00331        int    i;
00332        if (isleaf(node)) {         /* claim contribution */
00333               COLOR  contrib;
00334               if (node->val[EXP] <= 0)
00335                      return;              /* already claimed */
00336               if (fdir2diff(node->sd,cent) > maxr2)
00337                      return;              /* too far away */
00338               colr_color(contrib, node->val);
00339               scalecolor(contrib, tri_omegav(node->gdv));
00340               addcolor(intens, contrib);
00341               copycolr(node->val, blkclr);
00342               return;
00343        }
00344        if (dir2diff(node->gdv[0],node->gdv[1]) > maxr2 &&
00345                      fdir2diff(node->gdv[0],cent) < maxr2 &&
00346                      fdir2diff(node->gdv[1],cent) < maxr2 &&
00347                      fdir2diff(node->gdv[2],cent) < maxr2)
00348               return;                     /* previously claimed node */
00349        remaining = 0;                     /* recurse on children */
00350        for (i = 0; i < 4; i++) {
00351               claimlight(intens, &node->kid[i], cent, maxr2);
00352               if (!isleaf(&node->kid[i]) || node->kid[i].val[EXP] != 0)
00353                      ++remaining;
00354        }
00355        if (remaining)
00356               return;
00357                                    /* consolidate empties */
00358        free((void *)node->kid); node->kid = NULL;
00359        copycolr(node->val, blkclr);
00360        node->sd = node->gdv[0];    /* doesn't really matter */
00361 }
00362 
00363 /* Add lost light contribution to the given list */
00364 void
00365 add2lost(LOSTLIGHT **llp, COLOR intens, FVECT cent)
00366 {
00367        LOSTLIGHT     *newll = (LOSTLIGHT *)malloc(sizeof(LOSTLIGHT));
00368 
00369        if (newll == NULL)
00370               return;
00371        copycolor(newll->intens, intens);
00372        newll->sd = encodedir(cent);
00373        newll->next = *llp;
00374        *llp = newll;
00375 }
00376 
00377 /* Check lost light list for contributions */
00378 void
00379 getlost(LOSTLIGHT **llp, COLOR intens, FVECT cent, double omega)
00380 {
00381        const double  maxr2 = omega/PI;
00382        LOSTLIGHT     lhead, *lastp, *thisp;
00383 
00384        lhead.next = *llp;
00385        lastp = &lhead;
00386        while ((thisp = lastp->next) != NULL)
00387               if (fdir2diff(thisp->sd,cent) <= maxr2) {
00388                      LOSTLIGHT     *mynext = thisp->next;
00389                      addcolor(intens, thisp->intens);
00390                      free((void *)thisp);
00391                      lastp->next = mynext;
00392               } else
00393                      lastp = thisp;
00394        *llp = lhead.next;
00395 }
00396 
00397 /* Create & print distant sources */
00398 void
00399 mksources(TRITREE *samptree, double thresh, double maxang)
00400 {
00401        const int     ethresh = (int)(log(thresh)/log(2.) + (COLXS+.5));
00402        const double  maxomega = 2.*PI*(1. - cos(PI/180./2.*maxang));
00403        const double  minintens = .05*thresh*maxomega;
00404        int           nsrcs = 0;
00405        LOSTLIGHT     *lostlightlist = NULL;
00406        int           emax;
00407        TRITREE              *startleaf;
00408        double        growstep;
00409        FVECT         curcent;
00410        double        currad;
00411        double        curomega;
00412        COLOR         curintens;
00413        double        thisthresh;
00414        int           thisethresh;
00415        int           i;
00416        /*
00417         * General algorithm:
00418         *     1) Start with brightest unclaimed pixel
00419         *     2) Grow source toward brightest unclaimed perimeter until:
00420         *            a) Source exceeds maximum size, or
00421         *            b) Perimeter values all below threshold, or
00422         *            c) Source average drops below threshold
00423         *     3) Loop until nothing over threshold
00424         *
00425         * Complexity added to absorb insignificant sources in larger ones.
00426         */
00427        if (thresh <= FTINY)
00428               return;
00429        for ( ; ; ) {
00430               emax = ethresh;             /* find brightest unclaimed */
00431               startleaf = NULL;
00432               for (i = 0; i < NTRUNKBR; i++) {
00433                      TRITREE       *bigger = findemax(&samptree->kid[i], &emax);
00434                      if (bigger != NULL)
00435                             startleaf = bigger;
00436               }
00437               if (startleaf == NULL)
00438                      break;
00439                                    /* claim it */
00440               decodedir(curcent, startleaf->sd);
00441               curomega = tri_omegav(startleaf->gdv);
00442               currad = sqrt(curomega/PI);
00443               growstep = 3.*currad;
00444               colr_color(curintens, startleaf->val);
00445               thisthresh = .15*bright(curintens);
00446               if (thisthresh < thresh) thisthresh = thresh;
00447               thisethresh = (int)(log(thisthresh)/log(2.) + (COLXS+.5));
00448               scalecolor(curintens, curomega);
00449               copycolr(startleaf->val, blkclr);
00450               do {                 /* grow source */
00451                      FVECT  vsum;
00452                      double movedist;
00453                      vsum[0] = vsum[1] = vsum[2] = .0;
00454                      for (i = 0; i < NTRUNKBR; i++)
00455                             vector_sum(vsum, &samptree->kid[i],
00456                                    curcent, 2.-2.*cos(currad+growstep),
00457                                           thisethresh);
00458                      if (normalize(vsum) == .0)
00459                             break;
00460                      movedist = acos(DOT(vsum,curcent));
00461                      if (movedist > growstep) {
00462                             VSUB(vsum, vsum, curcent);
00463                             movedist = growstep/VLEN(vsum); 
00464                             VSUM(curcent, curcent, vsum, movedist);
00465                             normalize(curcent);
00466                      } else
00467                             VCOPY(curcent, vsum);
00468                      currad += growstep;
00469                      curomega = 2.*PI*(1. - cos(currad));
00470                      for (i = 0; i < NTRUNKBR; i++)
00471                             claimlight(curintens, &samptree->kid[i],
00472                                           curcent, 2.-2.*cos(currad));
00473               } while (curomega < maxomega &&
00474                             bright(curintens)/curomega > thisthresh);
00475               if (bright(curintens) < minintens) {
00476                      add2lost(&lostlightlist, curintens, curcent);
00477                      continue;
00478               }
00479                                    /* absorb lost contributions */
00480               getlost(&lostlightlist, curintens, curcent, curomega);
00481               ++nsrcs;             /* output source */
00482               scalecolor(curintens, 1./curomega);
00483               printf("\nvoid illum IBLout\n");
00484               printf("0\n0\n3 %f %f %f\n",
00485                             colval(curintens,RED),
00486                             colval(curintens,GRN),
00487                             colval(curintens,BLU));
00488               printf("\nIBLout source IBLsrc%d\n", nsrcs);
00489               printf("0\n0\n4 %f %f %f %f\n",
00490                             curcent[0], curcent[1], curcent[2],
00491                             2.*180./PI*currad);
00492        }
00493 }
00494 
00495 int
00496 main(int argc, char *argv[])
00497 {
00498        long   nsamps = DEF_NSAMPS;
00499        double maxang = DEF_MAXANG;
00500        TRITREE       *samptree;
00501        double thresh = 0;
00502        int    i;
00503        
00504        progname = argv[0];
00505        for (i = 1; i < argc && argv[i][0] == '-'; i++)
00506               switch (argv[i][1]) {
00507               case 'd':            /* number of samples */
00508                      if (i >= argc-1) goto userr;
00509                      nsamps = atol(argv[++i]);
00510                      break;
00511               case 't':            /* manual threshold */
00512                      if (i >= argc-1) goto userr;
00513                      thresh = atof(argv[++i]);
00514                      break;
00515               case 'a':            /* maximum source angle */
00516                      if (i >= argc-1) goto userr;
00517                      maxang = atof(argv[++i]);
00518                      if (maxang <= FTINY)
00519                             goto userr;
00520                      if (maxang > 180.)
00521                             maxang = 180.;
00522                      break;
00523               default:
00524                      goto userr;
00525               }
00526        if (i < argc-1)
00527               goto userr;
00528                                    /* start our ray calculation */
00529        directvis = 0;
00530        ray_init(i == argc-1 ? argv[i] : (char *)NULL);
00531        VCOPY(scene_cent, thescene.cuorg);
00532        scene_cent[0] += 0.5*thescene.cusize;
00533        scene_cent[1] += 0.5*thescene.cusize;
00534        scene_cent[2] += 0.5*thescene.cusize;
00535        scene_rad = 0.86603*thescene.cusize;
00536                                    /* sample geodesic mesh */
00537        samptree = geosample(nsamps);
00538                                    /* get source threshold */
00539        if (thresh <= FTINY)
00540               thresh = get_threshold(samptree);
00541                                    /* done with ray samples */
00542        ray_done(1);
00543                                    /* print header */
00544        printf("# ");
00545        printargs(argc, argv, stdout);
00546                                    /* create & print sources */
00547        mksources(samptree, thresh, maxang);
00548                                    /* all done, no need to clean up */
00549        return(0);
00550 userr:
00551        fprintf(stderr, "Usage: %s [-d nsamps][-t thresh][-a maxang] [octree]\n",
00552                      argv[0]);
00553        exit(1);
00554 }
00555 
00556 void
00557 eputs(char  *s)
00558 {
00559        static int  midline = 0;
00560 
00561        if (!*s)
00562               return;
00563        if (!midline++) {
00564               fputs(progname, stderr);
00565               fputs(": ", stderr);
00566        }
00567        fputs(s, stderr);
00568        if (s[strlen(s)-1] == '\n') {
00569               fflush(stderr);
00570               midline = 0;
00571        }
00572 }
00573 
00574 void
00575 wputs(char *s)
00576 {
00577        /* no warnings */
00578 }