Back to index

radiance  4R0+20100331
virtuals.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: virtuals.c,v 2.17 2006/09/07 05:20:54 greg Exp $";
00003 #endif
00004 /*
00005  * Routines for simulating virtual light sources
00006  *     Thus far, we only support planar mirrors.
00007  *
00008  *  External symbols declared in source.h
00009  */
00010 
00011 #include "copyright.h"
00012 
00013 #include  "ray.h"
00014 
00015 #include  "otypes.h"
00016 
00017 #include  "source.h"
00018 
00019 #include  "random.h"
00020 
00021 #define  MINSAMPLES  16            /* minimum number of pretest samples */
00022 #define  STESTMAX    32            /* maximum seeks per sample */
00023 
00024 
00025 static OBJECT  *vobject;           /* virtual source objects */
00026 static int  nvobjects = 0;         /* number of virtual source objects */
00027 
00028 
00029 extern void
00030 markvirtuals(void)                 /* find and mark virtual sources */
00031 {
00032        register OBJREC  *o;
00033        register int  i;
00034                                    /* check number of direct relays */
00035        if (directrelay <= 0)
00036               return;
00037                                    /* find virtual source objects */
00038        for (i = 0; i < nsceneobjs; i++) {
00039               o = objptr(i);
00040               if (!issurface(o->otype) || o->omod == OVOID)
00041                      continue;
00042               if (!isvlight(vsmaterial(o)->otype))
00043                      continue;
00044               if (sfun[o->otype].of == NULL ||
00045                             sfun[o->otype].of->getpleq == NULL) {
00046                      objerror(o,WARNING,"secondary sources not supported");
00047                      continue;
00048               }
00049               if (nvobjects == 0)
00050                      vobject = (OBJECT *)malloc(sizeof(OBJECT));
00051               else
00052                      vobject = (OBJECT *)realloc((void *)vobject,
00053                             (unsigned)(nvobjects+1)*sizeof(OBJECT));
00054               if (vobject == NULL)
00055                      error(SYSTEM, "out of memory in addvirtuals");
00056               vobject[nvobjects++] = i;
00057        }
00058        if (nvobjects == 0)
00059               return;
00060 #ifdef DEBUG
00061        fprintf(stderr, "found %d virtual source objects\n", nvobjects);
00062 #endif
00063                                    /* append virtual sources */
00064        for (i = nsources; i-- > 0; )
00065               addvirtuals(i, directrelay);
00066                                    /* done with our object list */
00067        free((void *)vobject);
00068        nvobjects = 0;
00069 }
00070 
00071 
00072 extern void
00073 addvirtuals(         /* add virtuals associated with source */
00074        int  sn,
00075        int  nr
00076 )
00077 {
00078        register int  i;
00079                             /* check relay limit first */
00080        if (nr <= 0)
00081               return;
00082        if (source[sn].sflags & SSKIP)
00083               return;
00084                             /* check each virtual object for projection */
00085        for (i = 0; i < nvobjects; i++)
00086                                    /* vproject() calls us recursively */
00087               vproject(objptr(vobject[i]), sn, nr-1);
00088 }
00089 
00090 
00091 extern void
00092 vproject(            /* create projected source(s) if they exist */
00093        OBJREC  *o,
00094        int  sn,
00095        int  n
00096 )
00097 {
00098        register int  i;
00099        register VSMATERIAL  *vsmat;
00100        MAT4  proj;
00101        int  ns;
00102 
00103        if (o == source[sn].so)     /* objects cannot project themselves */
00104               return;
00105                             /* get virtual source material */
00106        vsmat = sfun[vsmaterial(o)->otype].mf;
00107                             /* project virtual sources */
00108        for (i = 0; i < vsmat->nproj; i++)
00109               if ((*vsmat->vproj)(proj, o, &source[sn], i))
00110                      if ((ns = makevsrc(o, sn, proj)) >= 0) {
00111                             source[ns].sa.sv.pn = i;
00112 #ifdef DEBUG
00113                             virtverb(ns, stderr);
00114 #endif
00115                             addvirtuals(ns, n);
00116                      }
00117 }
00118 
00119 
00120 extern OBJREC *
00121 vsmaterial(                 /* get virtual source material pointer */
00122        OBJREC  *o
00123 )
00124 {
00125        register int  i;
00126        register OBJREC  *m;
00127 
00128        i = o->omod;
00129        m = findmaterial(objptr(i));
00130        if (m == NULL)
00131               return(objptr(i));
00132        if (m->otype != MAT_ILLUM || m->oargs.nsargs < 1 ||
00133                      !strcmp(m->oargs.sarg[0], VOIDID) ||
00134                      (i = lastmod(objndx(m), m->oargs.sarg[0])) == OVOID)
00135               return(m);           /* direct modifier */
00136        return(objptr(i));          /* illum alternate */
00137 }
00138 
00139 
00140 extern int
00141 makevsrc(            /* make virtual source if reasonable */
00142        OBJREC  *op,
00143        register int  sn,
00144        MAT4  pm
00145 )
00146 {
00147        FVECT  nsloc, nsnorm, ocent, v;
00148        double  maxrad2, d;
00149        int  nsflags;
00150        SPOT  theirspot, ourspot;
00151        register int  i;
00152 
00153        nsflags = source[sn].sflags | (SVIRTUAL|SSPOT|SFOLLOW);
00154                                    /* get object center and max. radius */
00155        maxrad2 = getdisk(ocent, op, sn);
00156        if (maxrad2 <= FTINY)                     /* too small? */
00157               return(-1);
00158                                    /* get location and spot */
00159        if (source[sn].sflags & SDISTANT) {              /* distant source */
00160               if (source[sn].sflags & SPROX)
00161                      return(-1);          /* should never get here! */
00162               multv3(nsloc, source[sn].sloc, pm);
00163               normalize(nsloc);
00164               VCOPY(ourspot.aim, ocent);
00165               ourspot.siz = PI*maxrad2;
00166               ourspot.flen = -1.;
00167               if (source[sn].sflags & SSPOT) {
00168                      multp3(theirspot.aim, source[sn].sl.s->aim, pm);
00169                                           /* adjust for source size */
00170                      d = sqrt(dist2(ourspot.aim, theirspot.aim));
00171                      d = sqrt(source[sn].sl.s->siz/PI) + d*source[sn].srad;
00172                      theirspot.siz = PI*d*d;
00173                      ourspot.flen = theirspot.flen = source[sn].sl.s->flen;
00174                      d = ourspot.siz;
00175                      if (!commonbeam(&ourspot, &theirspot, nsloc))
00176                             return(-1);   /* no overlap */
00177                      if (ourspot.siz < d-FTINY) {       /* it shrunk */
00178                             d = beamdisk(v, op, &ourspot, nsloc);
00179                             if (d <= FTINY)
00180                                    return(-1);
00181                             if (d < maxrad2) {
00182                                    maxrad2 = d;
00183                                    VCOPY(ocent, v);
00184                             }
00185                      }
00186               }
00187        } else {                           /* local source */
00188               multp3(nsloc, source[sn].sloc, pm);
00189               for (i = 0; i < 3; i++)
00190                      ourspot.aim[i] = ocent[i] - nsloc[i];
00191               if ((d = normalize(ourspot.aim)) == 0.)
00192                      return(-1);          /* at source!! */
00193               if (source[sn].sflags & SPROX && d > source[sn].sl.prox)
00194                      return(-1);          /* too far away */
00195               ourspot.flen = 0.;
00196                                           /* adjust for source size */
00197               d = (sqrt(maxrad2) + source[sn].srad) / d;
00198               if (d < 1.-FTINY)
00199                      ourspot.siz = 2.*PI*(1. - sqrt(1.-d*d));
00200               else
00201                      nsflags &= ~SSPOT;
00202               if (source[sn].sflags & SSPOT) {
00203                      theirspot = *(source[sn].sl.s);
00204                      multv3(theirspot.aim, source[sn].sl.s->aim, pm);
00205                      normalize(theirspot.aim);
00206                      if (nsflags & SSPOT) {
00207                             ourspot.flen = theirspot.flen;
00208                             d = ourspot.siz;
00209                             if (!commonspot(&ourspot, &theirspot, nsloc))
00210                                    return(-1);   /* no overlap */
00211                      } else {
00212                             nsflags |= SSPOT;
00213                             ourspot = theirspot;
00214                             d = 2.*ourspot.siz;
00215                      }
00216                      if (ourspot.siz < d-FTINY) {       /* it shrunk */
00217                             d = spotdisk(v, op, &ourspot, nsloc);
00218                             if (d <= FTINY)
00219                                    return(-1);
00220                             if (d < maxrad2) {
00221                                    maxrad2 = d;
00222                                    VCOPY(ocent, v);
00223                             }
00224                      }
00225               }
00226               if (source[sn].sflags & SFLAT) {   /* behind source? */
00227                      multv3(nsnorm, source[sn].snorm, pm);
00228                      normalize(nsnorm);
00229                      if (nsflags & SSPOT && !checkspot(&ourspot, nsnorm))
00230                             return(-1);
00231               }
00232        }
00233                                    /* pretest visibility */
00234        nsflags = vstestvis(nsflags, op, ocent, maxrad2, sn);
00235        if (nsflags & SSKIP)
00236               return(-1);   /* obstructed */
00237                                    /* it all checks out, so make it */
00238        if ((i = newsource()) < 0)
00239               goto memerr;
00240        source[i].sflags = nsflags;
00241        VCOPY(source[i].sloc, nsloc);
00242        multv3(source[i].ss[SU], source[sn].ss[SU], pm);
00243        multv3(source[i].ss[SV], source[sn].ss[SV], pm);
00244        if (nsflags & SFLAT)
00245               VCOPY(source[i].snorm, nsnorm);
00246        else
00247               multv3(source[i].ss[SW], source[sn].ss[SW], pm);
00248        source[i].srad = source[sn].srad;
00249        source[i].ss2 = source[sn].ss2;
00250        if (nsflags & SSPOT) {
00251               if ((source[i].sl.s = (SPOT *)malloc(sizeof(SPOT))) == NULL)
00252                      goto memerr;
00253               *(source[i].sl.s) = ourspot;
00254        }
00255        if (nsflags & SPROX)
00256               source[i].sl.prox = source[sn].sl.prox;
00257        source[i].sa.sv.sn = sn;
00258        source[i].so = op;
00259        return(i);
00260 memerr:
00261        error(SYSTEM, "out of memory in makevsrc");
00262        return -1; /* pro forma return */
00263 }
00264 
00265 
00266 extern double
00267 getdisk(             /* get visible object disk */
00268        FVECT  oc,
00269        OBJREC  *op,
00270        register int  sn
00271 )
00272 {
00273        double  rad2, roffs, offs, d, rd, rdoto;
00274        FVECT  rnrm, nrm;
00275                             /* first, use object getdisk function */
00276        rad2 = getmaxdisk(oc, op);
00277        if (!(source[sn].sflags & SVIRTUAL))
00278               return(rad2);        /* all done for normal source */
00279                             /* check for correct side of relay surface */
00280        roffs = getplaneq(rnrm, source[sn].so);
00281        rd = DOT(rnrm, source[sn].sloc);   /* source projection */
00282        if (!(source[sn].sflags & SDISTANT))
00283               rd -= roffs;
00284        d = DOT(rnrm, oc) - roffs;  /* disk distance to relay plane */
00285        if ((d > 0.) ^ (rd > 0.))
00286               return(rad2);        /* OK if opposite sides */
00287        if (d*d >= rad2)
00288               return(0.);          /* no relay is possible */
00289                             /* we need a closer look */
00290        offs = getplaneq(nrm, op);
00291        rdoto = DOT(rnrm, nrm);
00292        if (d*d >= rad2*(1.-rdoto*rdoto))
00293               return(0.);          /* disk entirely on projection side */
00294                             /* should shrink disk but I'm lazy */
00295        return(rad2);
00296 }
00297 
00298 
00299 extern int
00300 vstestvis(           /* pretest source visibility */
00301        int  f,                     /* virtual source flags */
00302        OBJREC  *o,          /* relay object */
00303        FVECT  oc,           /* relay object center */
00304        double  or2,         /* relay object radius squared */
00305        register int  sn     /* target source number */
00306 )
00307 {
00308        RAY  sr;
00309        FVECT  onorm;
00310        FVECT  offsdir;
00311        SRCINDEX  si;
00312        double  or, d, d1;
00313        int  stestlim, ssn;
00314        int  nhit, nok;
00315        register int  i, n;
00316                             /* return if pretesting disabled */
00317        if (vspretest <= 0)
00318               return(f);
00319                             /* get surface normal */
00320        getplaneq(onorm, o);
00321                             /* set number of rays to sample */
00322        if (source[sn].sflags & SDISTANT) {
00323                                    /* 32. == heuristic constant */
00324               n = 32.*or2/(thescene.cusize*thescene.cusize)*vspretest + .5;
00325        } else {
00326               for (i = 0; i < 3; i++)
00327                      offsdir[i] = source[sn].sloc[i] - oc[i];
00328               d = DOT(offsdir,offsdir);
00329               if (d <= FTINY)
00330                      n = 2.*PI * vspretest + .5;
00331               else
00332                      n = 2.*PI * (1.-sqrt(1./(1.+or2/d)))*vspretest + .5;
00333        }
00334        if (n < MINSAMPLES) n = MINSAMPLES;
00335 #ifdef DEBUG
00336        fprintf(stderr, "pretesting source %d in object %s with %d rays\n",
00337                      sn, o->oname, n);
00338 #endif
00339                             /* sample */
00340        or = sqrt(or2);
00341        stestlim = n*STESTMAX;
00342        ssn = 0;
00343        nhit = nok = 0;
00344        initsrcindex(&si);
00345        while (n-- > 0) {
00346                                    /* get sample point */
00347               do {
00348                      if (ssn >= stestlim) {
00349 #ifdef DEBUG
00350                             fprintf(stderr, "\ttoo hard to hit\n");
00351 #endif
00352                             return(f);    /* too small a target! */
00353                      }
00354                      multisamp(offsdir, 3, urand(sn*931+5827+ssn));
00355                      for (i = 0; i < 3; i++)
00356                             offsdir[i] = or*(1. - 2.*offsdir[i]);
00357                      ssn++;
00358                      d = 1. - DOT(offsdir, onorm);
00359                      for (i = 0; i < 3; i++) {
00360                             sr.rorg[i] = oc[i] + offsdir[i] + d*onorm[i];
00361                             sr.rdir[i] = -onorm[i];
00362                      }
00363                      sr.rmax = 0.0;
00364                      rayorigin(&sr, PRIMARY, NULL, NULL);
00365               } while (!(*ofun[o->otype].funp)(o, &sr));
00366                                    /* check against source */
00367               VCOPY(sr.rorg, sr.rop);     /* starting from intersection */
00368               samplendx++;
00369               if (si.sp >= si.np-1 ||
00370                             !srcray(&sr, NULL, &si) || sr.rsrc != sn) {
00371                      si.sn = sn-1;        /* reset index to our source */
00372                      si.np = 0;
00373                      if (!srcray(&sr, NULL, &si) || sr.rsrc != sn)
00374                             continue;     /* can't get there from here */
00375               }
00376               sr.revf = srcvalue;
00377               rayvalue(&sr);                     /* check sample validity */
00378               if ((d = bright(sr.rcol)) <= FTINY)
00379                      continue;
00380               nok++;               /* got sample; check obstructions */
00381               rayclear(&sr);
00382               sr.revf = raytrace;
00383               rayvalue(&sr);
00384               if ((d1 = bright(sr.rcol)) > FTINY) {
00385                      if (d - d1 > FTINY) {
00386 #ifdef DEBUG
00387                             fprintf(stderr, "\tpartially shadowed\n");
00388 #endif
00389                             return(f);    /* intervening transmitter */
00390                      }
00391                      nhit++;
00392               }
00393               if (nhit > 0 && nhit < nok) {
00394 #ifdef DEBUG
00395                      fprintf(stderr, "\tpartially occluded\n");
00396 #endif
00397                      return(f);           /* need to shadow test */
00398               }
00399        }
00400        if (nhit == 0) {
00401 #ifdef DEBUG
00402               fprintf(stderr, "\t0%% hit rate\n");
00403 #endif
00404               return(f | SSKIP);   /* 0% hit rate:  totally occluded */
00405        }
00406 #ifdef DEBUG
00407        fprintf(stderr, "\t100%% hit rate\n");
00408 #endif
00409        return(f & ~SFOLLOW);              /* 100% hit rate:  no occlusion */
00410 }
00411        
00412 
00413 #ifdef DEBUG
00414 extern void
00415 virtverb(     /* print verbose description of virtual source */
00416        register int  sn,
00417        FILE  *fp
00418 )
00419 {
00420        register int  i;
00421 
00422        fprintf(fp, "%s virtual source %d in %s %s\n",
00423                      source[sn].sflags & SDISTANT ? "distant" : "local",
00424                      sn, ofun[source[sn].so->otype].funame,
00425                      source[sn].so->oname);
00426        fprintf(fp, "\tat (%f,%f,%f)\n",
00427               source[sn].sloc[0], source[sn].sloc[1], source[sn].sloc[2]);
00428        fprintf(fp, "\tlinked to source %d (%s)\n",
00429               source[sn].sa.sv.sn, source[source[sn].sa.sv.sn].so->oname);
00430        if (source[sn].sflags & SFOLLOW)
00431               fprintf(fp, "\talways followed\n");
00432        else
00433               fprintf(fp, "\tnever followed\n");
00434        if (!(source[sn].sflags & SSPOT))
00435               return;
00436        fprintf(fp, "\twith spot aim (%f,%f,%f) and size %f\n",
00437                      source[sn].sl.s->aim[0], source[sn].sl.s->aim[1],
00438                      source[sn].sl.s->aim[2], source[sn].sl.s->siz);
00439 }
00440 #endif