Back to index

radiance  4R0+20100331
srcsupp.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: srcsupp.c,v 2.17 2008/12/06 01:08:53 greg Exp $";
00003 #endif
00004 /*
00005  *  Support routines for source objects and materials
00006  *
00007  *  External symbols declared in source.h
00008  */
00009 
00010 #include "copyright.h"
00011 
00012 #include  "ray.h"
00013 
00014 #include  "otypes.h"
00015 
00016 #include  "source.h"
00017 
00018 #include  "cone.h"
00019 
00020 #include  "face.h"
00021 
00022 #define SRCINC              8             /* realloc increment for array */
00023 
00024 SRCREC  *source = NULL;                   /* our list of sources */
00025 int  nsources = 0;                 /* the number of sources */
00026 
00027 SRCFUNC  sfun[NUMOTYPE];           /* source dispatch table */
00028 
00029 
00030 void
00031 initstypes()                /* initialize source dispatch table */
00032 {
00033        extern VSMATERIAL  mirror_vs, direct1_vs, direct2_vs;
00034        static SOBJECT  fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
00035        static SOBJECT  ssobj = {ssetsrc, nopart};
00036        static SOBJECT  sphsobj = {sphsetsrc, nopart};
00037        static SOBJECT  cylsobj = {cylsetsrc, cylpart};
00038        static SOBJECT  rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};
00039 
00040        sfun[MAT_MIRROR].mf = &mirror_vs;
00041        sfun[MAT_DIRECT1].mf = &direct1_vs;
00042        sfun[MAT_DIRECT2].mf = &direct2_vs;
00043        sfun[OBJ_FACE].of = &fsobj;
00044        sfun[OBJ_SOURCE].of = &ssobj;
00045        sfun[OBJ_SPHERE].of = &sphsobj;
00046        sfun[OBJ_CYLINDER].of = &cylsobj;
00047        sfun[OBJ_RING].of = &rsobj;
00048 }
00049 
00050 
00051 int
00052 newsource()                 /* allocate new source in our array */
00053 {
00054        if (nsources == 0)
00055               source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
00056        else if (nsources%SRCINC == 0)
00057               source = (SRCREC *)realloc((void *)source,
00058                             (unsigned)(nsources+SRCINC)*sizeof(SRCREC));
00059        if (source == NULL)
00060               return(-1);
00061        source[nsources].sflags = 0;
00062        source[nsources].nhits = 1;
00063        source[nsources].ntests = 2;       /* initial hit probability = 50% */
00064 #if SHADCACHE
00065        source[nsources].obscache = NULL;
00066 #endif
00067        return(nsources++);
00068 }
00069 
00070 
00071 void
00072 setflatss(src)                            /* set sampling for a flat source */
00073 register SRCREC  *src;
00074 {
00075        double  mult;
00076        register int  i;
00077 
00078        src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
00079        for (i = 0; i < 3; i++)
00080               if (src->snorm[i] < 0.6 && src->snorm[i] > -0.6)
00081                      break;
00082        src->ss[SV][i] = 1.0;
00083        fcross(src->ss[SU], src->ss[SV], src->snorm);
00084        mult = .5 * sqrt( src->ss2 / DOT(src->ss[SU],src->ss[SU]) );
00085        for (i = 0; i < 3; i++)
00086               src->ss[SU][i] *= mult;
00087        fcross(src->ss[SV], src->snorm, src->ss[SU]);
00088 }
00089 
00090 
00091 void
00092 fsetsrc(src, so)                   /* set a face as a source */
00093 register SRCREC  *src;
00094 OBJREC  *so;
00095 {
00096        register FACE  *f;
00097        register int  i, j;
00098        double  d;
00099        
00100        src->sa.success = 2*AIMREQT-1;            /* bitch on second failure */
00101        src->so = so;
00102                                           /* get the face */
00103        f = getface(so);
00104        if (f->area == 0.0)
00105               objerror(so, USER, "zero source area");
00106                                           /* find the center */
00107        for (j = 0; j < 3; j++) {
00108               src->sloc[j] = 0.0;
00109               for (i = 0; i < f->nv; i++)
00110                      src->sloc[j] += VERTEX(f,i)[j];
00111               src->sloc[j] /= (double)f->nv;
00112        }
00113        if (!inface(src->sloc, f))
00114               objerror(so, USER, "cannot hit source center");
00115        src->sflags |= SFLAT;
00116        VCOPY(src->snorm, f->norm);
00117        src->ss2 = f->area;
00118                                           /* find maximum radius */
00119        src->srad = 0.;
00120        for (i = 0; i < f->nv; i++) {
00121               d = dist2(VERTEX(f,i), src->sloc);
00122               if (d > src->srad)
00123                      src->srad = d;
00124        }
00125        src->srad = sqrt(src->srad);
00126                                           /* compute size vectors */
00127        if (f->nv == 4)                           /* parallelogram case */
00128               for (j = 0; j < 3; j++) {
00129                      src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
00130                      src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
00131               }
00132        else
00133               setflatss(src);
00134 }
00135 
00136 
00137 void
00138 ssetsrc(src, so)                   /* set a source as a source */
00139 register SRCREC  *src;
00140 register OBJREC  *so;
00141 {
00142        double  theta;
00143        
00144        src->sa.success = 2*AIMREQT-1;            /* bitch on second failure */
00145        src->so = so;
00146        if (so->oargs.nfargs != 4)
00147               objerror(so, USER, "bad arguments");
00148        src->sflags |= (SDISTANT|SCIR);
00149        VCOPY(src->sloc, so->oargs.farg);
00150        if (normalize(src->sloc) == 0.0)
00151               objerror(so, USER, "zero direction");
00152        theta = PI/180.0/2.0 * so->oargs.farg[3];
00153        if (theta <= FTINY)
00154               objerror(so, USER, "zero size");
00155        src->ss2 = 2.0*PI * (1.0 - cos(theta));
00156                                    /* the following is approximate */
00157        src->srad = sqrt(src->ss2/PI);
00158        VCOPY(src->snorm, src->sloc);
00159        setflatss(src);                    /* hey, whatever works */
00160        src->ss[SW][0] = src->ss[SW][1] = src->ss[SW][2] = 0.0;
00161 }
00162 
00163 
00164 void
00165 sphsetsrc(src, so)                 /* set a sphere as a source */
00166 register SRCREC  *src;
00167 register OBJREC  *so;
00168 {
00169        register int  i;
00170 
00171        src->sa.success = 2*AIMREQT-1;            /* bitch on second failure */
00172        src->so = so;
00173        if (so->oargs.nfargs != 4)
00174               objerror(so, USER, "bad # arguments");
00175        if (so->oargs.farg[3] <= FTINY)
00176               objerror(so, USER, "illegal source radius");
00177        src->sflags |= SCIR;
00178        VCOPY(src->sloc, so->oargs.farg);
00179        src->srad = so->oargs.farg[3];
00180        src->ss2 = PI * src->srad * src->srad;
00181        for (i = 0; i < 3; i++)
00182               src->ss[SU][i] = src->ss[SV][i] = src->ss[SW][i] = 0.0;
00183        for (i = 0; i < 3; i++)
00184               src->ss[i][i] = 0.7236 * so->oargs.farg[3];
00185 }
00186 
00187 
00188 void
00189 rsetsrc(src, so)                   /* set a ring (disk) as a source */
00190 register SRCREC  *src;
00191 OBJREC  *so;
00192 {
00193        register CONE  *co;
00194        
00195        src->sa.success = 2*AIMREQT-1;            /* bitch on second failure */
00196        src->so = so;
00197                                           /* get the ring */
00198        co = getcone(so, 0);
00199        if (CO_R1(co) <= FTINY)
00200               objerror(so, USER, "illegal source radius");
00201        VCOPY(src->sloc, CO_P0(co));
00202        if (CO_R0(co) > 0.0)
00203               objerror(so, USER, "cannot hit source center");
00204        src->sflags |= (SFLAT|SCIR);
00205        VCOPY(src->snorm, co->ad);
00206        src->srad = CO_R1(co);
00207        src->ss2 = PI * src->srad * src->srad;
00208        setflatss(src);
00209 }
00210 
00211 
00212 void
00213 cylsetsrc(src, so)                 /* set a cylinder as a source */
00214 register SRCREC  *src;
00215 OBJREC  *so;
00216 {
00217        register CONE  *co;
00218        register int  i;
00219        
00220        src->sa.success = 4*AIMREQT-1;            /* bitch on fourth failure */
00221        src->so = so;
00222                                           /* get the cylinder */
00223        co = getcone(so, 0);
00224        if (CO_R0(co) <= FTINY)
00225               objerror(so, USER, "illegal source radius");
00226        if (CO_R0(co) > .2*co->al)         /* heuristic constraint */
00227               objerror(so, WARNING, "source aspect too small");
00228        src->sflags |= SCYL;
00229        for (i = 0; i < 3; i++)
00230               src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
00231        src->srad = .5*co->al;
00232        src->ss2 = 2.*CO_R0(co)*co->al;
00233                                           /* set sampling vectors */
00234        for (i = 0; i < 3; i++)
00235               src->ss[SU][i] = .5 * co->al * co->ad[i];
00236        src->ss[SV][0] = src->ss[SV][1] = src->ss[SV][2] = 0.0;
00237        for (i = 0; i < 3; i++)
00238               if (co->ad[i] < 0.6 && co->ad[i] > -0.6)
00239                      break;
00240        src->ss[SV][i] = 1.0;
00241        fcross(src->ss[SW], src->ss[SV], co->ad);
00242        normalize(src->ss[SW]);
00243        for (i = 0; i < 3; i++)
00244               src->ss[SW][i] *= .8559 * CO_R0(co);
00245        fcross(src->ss[SV], src->ss[SW], co->ad);
00246 }
00247 
00248 
00249 SPOT *
00250 makespot(m)                 /* make a spotlight */
00251 register OBJREC  *m;
00252 {
00253        register SPOT  *ns;
00254 
00255        if ((ns = (SPOT *)m->os) != NULL)
00256               return(ns);
00257        if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
00258               return(NULL);
00259        if (m->oargs.farg[3] <= FTINY)
00260               objerror(m, USER, "zero angle");
00261        ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
00262        VCOPY(ns->aim, m->oargs.farg+4);
00263        if ((ns->flen = normalize(ns->aim)) == 0.0)
00264               objerror(m, USER, "zero focus vector");
00265        m->os = (char *)ns;
00266        return(ns);
00267 }
00268 
00269 
00270 int
00271 spotout(r, s)               /* check if we're outside spot region */
00272 register RAY  *r;
00273 register SPOT  *s;
00274 {
00275        double  d;
00276        FVECT  vd;
00277        
00278        if (s == NULL)
00279               return(0);
00280        if (s->flen < -FTINY) {            /* distant source */
00281               vd[0] = s->aim[0] - r->rorg[0];
00282               vd[1] = s->aim[1] - r->rorg[1];
00283               vd[2] = s->aim[2] - r->rorg[2];
00284               d = DOT(r->rdir,vd);
00285               /*                   wrong side?
00286               if (d <= FTINY)
00287                      return(1);    */
00288               d = DOT(vd,vd) - d*d;
00289               if (PI*d > s->siz)
00290                      return(1);    /* out */
00291               return(0);    /* OK */
00292        }
00293                                    /* local source */
00294        if (s->siz < 2.0*PI * (1.0 + DOT(s->aim,r->rdir)))
00295               return(1);    /* out */
00296        return(0);    /* OK */
00297 }
00298 
00299 
00300 double
00301 fgetmaxdisk(ocent, op)             /* get center and squared radius of face */
00302 FVECT  ocent;
00303 OBJREC  *op;
00304 {
00305        double  maxrad2;
00306        double  d;
00307        register int  i, j;
00308        register FACE  *f;
00309        
00310        f = getface(op);
00311        if (f->area == 0.)
00312               return(0.);
00313        for (i = 0; i < 3; i++) {
00314               ocent[i] = 0.;
00315               for (j = 0; j < f->nv; j++)
00316                      ocent[i] += VERTEX(f,j)[i];
00317               ocent[i] /= (double)f->nv;
00318        }
00319        d = DOT(ocent,f->norm);
00320        for (i = 0; i < 3; i++)
00321               ocent[i] += (f->offset - d)*f->norm[i];
00322        maxrad2 = 0.;
00323        for (j = 0; j < f->nv; j++) {
00324               d = dist2(VERTEX(f,j), ocent);
00325               if (d > maxrad2)
00326                      maxrad2 = d;
00327        }
00328        return(maxrad2);
00329 }
00330 
00331 
00332 double
00333 rgetmaxdisk(ocent, op)             /* get center and squared radius of ring */
00334 FVECT  ocent;
00335 OBJREC  *op;
00336 {
00337        register CONE  *co;
00338        
00339        co = getcone(op, 0);
00340        VCOPY(ocent, CO_P0(co));
00341        return(CO_R1(co)*CO_R1(co));
00342 }
00343 
00344 
00345 double
00346 fgetplaneq(nvec, op)               /* get plane equation for face */
00347 FVECT  nvec;
00348 OBJREC  *op;
00349 {
00350        register FACE  *fo;
00351 
00352        fo = getface(op);
00353        VCOPY(nvec, fo->norm);
00354        return(fo->offset);
00355 }
00356 
00357 
00358 double
00359 rgetplaneq(nvec, op)               /* get plane equation for ring */
00360 FVECT  nvec;
00361 OBJREC  *op;
00362 {
00363        register CONE  *co;
00364 
00365        co = getcone(op, 0);
00366        VCOPY(nvec, co->ad);
00367        return(DOT(nvec, CO_P0(co)));
00368 }
00369 
00370 
00371 int
00372 commonspot(sp1, sp2, org)   /* set sp1 to intersection of sp1 and sp2 */
00373 register SPOT  *sp1, *sp2;
00374 FVECT  org;
00375 {
00376        FVECT  cent;
00377        double  rad2, cos1, cos2;
00378 
00379        cos1 = 1. - sp1->siz/(2.*PI);
00380        cos2 = 1. - sp2->siz/(2.*PI);
00381        if (sp2->siz >= 2.*PI-FTINY)              /* BIG, just check overlap */
00382               return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
00383                                    sqrt((1.-cos1*cos1)*(1.-cos2*cos2)));
00384                             /* compute and check disks */
00385        rad2 = intercircle(cent, sp1->aim, sp2->aim,
00386                      1./(cos1*cos1) - 1.,  1./(cos2*cos2) - 1.);
00387        if (rad2 <= FTINY || normalize(cent) == 0.)
00388               return(0);
00389        VCOPY(sp1->aim, cent);
00390        sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2));
00391        return(1);
00392 }
00393 
00394 
00395 int
00396 commonbeam(sp1, sp2, dir)   /* set sp1 to intersection of sp1 and sp2 */
00397 register SPOT  *sp1, *sp2;
00398 FVECT  dir;
00399 {
00400        FVECT  cent, c1, c2;
00401        double  rad2, d;
00402        register int  i;
00403                                    /* move centers to common plane */
00404        d = DOT(sp1->aim, dir);
00405        for (i = 0; i < 3; i++)
00406               c1[i] = sp1->aim[i] - d*dir[i];
00407        d = DOT(sp2->aim, dir);
00408        for (i = 0; i < 3; i++)
00409               c2[i] = sp2->aim[i] - d*dir[i];
00410                                    /* compute overlap */
00411        rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
00412        if (rad2 <= FTINY)
00413               return(0);
00414        VCOPY(sp1->aim, cent);
00415        sp1->siz = PI*rad2;
00416        return(1);
00417 }
00418 
00419 
00420 int
00421 checkspot(sp, nrm)          /* check spotlight for behind source */
00422 register SPOT  *sp;  /* spotlight */
00423 FVECT  nrm;          /* source surface normal */
00424 {
00425        double  d, d1;
00426 
00427        d = DOT(sp->aim, nrm);
00428        if (d > FTINY)                     /* center in front? */
00429               return(1);
00430                                    /* else check horizon */
00431        d1 = 1. - sp->siz/(2.*PI);
00432        return(1.-FTINY-d*d < d1*d1);
00433 }
00434 
00435 
00436 double
00437 spotdisk(oc, op, sp, pos)   /* intersect spot with object op */
00438 FVECT  oc;
00439 OBJREC  *op;
00440 register SPOT  *sp;
00441 FVECT  pos;
00442 {
00443        FVECT  onorm;
00444        double  offs, d, dist;
00445        register int  i;
00446 
00447        offs = getplaneq(onorm, op);
00448        d = -DOT(onorm, sp->aim);
00449        if (d >= -FTINY && d <= FTINY)
00450               return(0.);
00451        dist = (DOT(pos, onorm) - offs)/d;
00452        if (dist < 0.)
00453               return(0.);
00454        for (i = 0; i < 3; i++)
00455               oc[i] = pos[i] + dist*sp->aim[i];
00456        return(sp->siz*dist*dist/PI/(d*d));
00457 }
00458 
00459 
00460 double
00461 beamdisk(oc, op, sp, dir)   /* intersect beam with object op */
00462 FVECT  oc;
00463 OBJREC  *op;
00464 register SPOT  *sp;
00465 FVECT  dir;
00466 {
00467        FVECT  onorm;
00468        double  offs, d, dist;
00469        register int  i;
00470 
00471        offs = getplaneq(onorm, op);
00472        d = -DOT(onorm, dir);
00473        if (d >= -FTINY && d <= FTINY)
00474               return(0.);
00475        dist = (DOT(sp->aim, onorm) - offs)/d;
00476        for (i = 0; i < 3; i++)
00477               oc[i] = sp->aim[i] + dist*dir[i];
00478        return(sp->siz/PI/(d*d));
00479 }
00480 
00481 
00482 double
00483 intercircle(cc, c1, c2, r1s, r2s)  /* intersect two circles */
00484 FVECT  cc;                  /* midpoint (return value) */
00485 FVECT  c1, c2;                     /* circle centers */
00486 double  r1s, r2s;           /* radii squared */
00487 {
00488        double  a2, d2, l;
00489        FVECT  disp;
00490        register int  i;
00491 
00492        for (i = 0; i < 3; i++)
00493               disp[i] = c2[i] - c1[i];
00494        d2 = DOT(disp,disp);
00495                                    /* circle within overlap? */
00496        if (r1s < r2s) {
00497               if (r2s >= r1s + d2) {
00498                      VCOPY(cc, c1);
00499                      return(r1s);
00500               }
00501        } else {
00502               if (r1s >= r2s + d2) {
00503                      VCOPY(cc, c2);
00504                      return(r2s);
00505               }
00506        }
00507        a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
00508                                    /* no overlap? */
00509        if (a2 <= 0.)
00510               return(0.);
00511                                    /* overlap, compute center */
00512        l = sqrt((r1s - a2)/d2);
00513        for (i = 0; i < 3; i++)
00514               cc[i] = c1[i] + l*disp[i];
00515        return(a2);
00516 }