--- ray/src/rt/virtuals.c 1991/06/20 09:37:38 1.2 +++ ray/src/rt/virtuals.c 1991/07/11 16:43:38 1.16 @@ -11,21 +11,20 @@ static char SCCSid[] = "$SunId$ LBL"; #include "ray.h" -#include "source.h" +#include "octree.h" #include "otypes.h" -#include "cone.h" +#include "source.h" -#include "face.h" +#include "random.h" -extern int directrelay; /* maximum number of source relays */ +#define MINSAMPLES 5 /* minimum number of pretest samples */ +#define STESTMAX 30 /* maximum seeks per sample */ -double getplaneq(); -double getmaxdisk(); -double intercircle(); -SRCREC *makevsrc(); +double getdisk(); + static OBJECT *vobject; /* virtual source objects */ static int nvobjects = 0; /* number of virtual source objects */ @@ -40,10 +39,13 @@ markvirtuals() /* find and mark virtual sources */ /* find virtual source objects */ for (i = 0; i < nobjects; i++) { o = objptr(i); - if (o->omod == OVOID) + if (!issurface(o->otype) || o->omod == OVOID) continue; if (!isvlight(objptr(o->omod)->otype)) continue; + if (sfun[o->otype].of == NULL || + sfun[o->otype].of->getpleq == NULL) + objerror(o, USER, "illegal material"); if (nvobjects == 0) vobject = (OBJECT *)malloc(sizeof(OBJECT)); else @@ -55,380 +57,330 @@ markvirtuals() /* find and mark virtual sources */ } if (nvobjects == 0) return; +#ifdef DEBUG + fprintf(stderr, "found %d virtual source objects\n", nvobjects); +#endif /* append virtual sources */ for (i = nsources; i-- > 0; ) - if (!(source[i].sflags & SSKIP)) - addvirtuals(&source[i], directrelay); + addvirtuals(i, directrelay); /* done with our object list */ free((char *)vobject); nvobjects = 0; } -addvirtuals(sr, nr) /* add virtual sources associated with sr */ -SRCREC *sr; +addvirtuals(sn, nr) /* add virtuals associated with source */ +int sn; int nr; { register int i; /* check relay limit first */ if (nr <= 0) return; + if (source[sn].sflags & SSKIP) + return; /* check each virtual object for projection */ for (i = 0; i < nvobjects; i++) - vproject(objptr(i), sr, nr-1); /* calls us recursively */ + /* vproject() calls us recursively */ + vproject(objptr(vobject[i]), sn, nr-1); } -SRCREC * -makevsrc(op, sp, pm) /* make virtual source if reasonable */ +vproject(o, sn, n) /* create projected source(s) if they exist */ +OBJREC *o; +int sn; +int n; +{ + register int i; + register VSMATERIAL *vsmat; + MAT4 proj; + int ns; + + if (o == source[sn].so) /* objects cannot project themselves */ + return; + /* get virtual source material */ + vsmat = sfun[objptr(o->omod)->otype].mf; + /* project virtual sources */ + for (i = 0; i < vsmat->nproj; i++) + if ((*vsmat->vproj)(proj, o, &source[sn], i)) + if ((ns = makevsrc(o, sn, proj)) >= 0) { +#ifdef DEBUG + virtverb(ns, stderr); +#endif + addvirtuals(ns, n); + } +} + + +int +makevsrc(op, sn, pm) /* make virtual source if reasonable */ OBJREC *op; -register SRCREC *sp; +register int sn; MAT4 pm; { - register SRCREC *newsrc; - FVECT nsloc, ocent, nsnorm; - double maxrad2; - double d1, d2; + FVECT nsloc, nsnorm, ocent, v; + double maxrad2, d; + int nsflags; SPOT theirspot, ourspot; register int i; + + nsflags = source[sn].sflags | (SVIRTUAL|SSPOT|SFOLLOW); /* get object center and max. radius */ - maxrad2 = getmaxdisk(ocent, op); + maxrad2 = getdisk(ocent, op, sn); if (maxrad2 <= FTINY) /* too small? */ - return(NULL); + return(-1); /* get location and spot */ - if (sp->sflags & SDISTANT) { /* distant source */ - if (sp->sflags & SPROX) - return(NULL); /* should never get here! */ - multv3(nsloc, sp->sloc, pm); + if (source[sn].sflags & SDISTANT) { /* distant source */ + if (source[sn].sflags & SPROX) + return(-1); /* should never get here! */ + multv3(nsloc, source[sn].sloc, pm); VCOPY(ourspot.aim, ocent); ourspot.siz = PI*maxrad2; ourspot.flen = 0.; - if (sp->sflags & SSPOT) { - copystruct(&theirspot, sp->sl.s); - multp3(theirspot.aim, sp->sl.s->aim, pm); + if (source[sn].sflags & SSPOT) { + copystruct(&theirspot, source[sn].sl.s); + multp3(theirspot.aim, source[sn].sl.s->aim, pm); + d = ourspot.siz; if (!commonbeam(&ourspot, &theirspot, nsloc)) - return(NULL); /* no overlap */ + return(-1); /* no overlap */ + if (ourspot.siz < d-FTINY) { /* it shrunk */ + d = beamdisk(v, op, &ourspot, nsloc); + if (d <= FTINY) + return(-1); + if (d < maxrad2) { + maxrad2 = d; + VCOPY(ocent, v); + } + } } } else { /* local source */ - multp3(nsloc, sp->sloc, pm); + multp3(nsloc, source[sn].sloc, pm); for (i = 0; i < 3; i++) ourspot.aim[i] = ocent[i] - nsloc[i]; - if ((d1 = normalize(ourspot.aim)) == 0.) - return(NULL); /* at source!! */ - if (sp->sflags & SPROX && d1 > sp->sl.prox) - return(NULL); /* too far away */ - ourspot.siz = 2.*PI*(1. - d1/sqrt(d1*d1+maxrad2)); + if ((d = normalize(ourspot.aim)) == 0.) + return(-1); /* at source!! */ + if (source[sn].sflags & SPROX && d > source[sn].sl.prox) + return(-1); /* too far away */ ourspot.flen = 0.; - if (sp->sflags & SSPOT) { - copystruct(&theirspot, sp->sl.s); - multv3(theirspot.aim, sp->sl.s->aim, pm); - if (!commonspot(&ourspot, &theirspot, nsloc)) - return(NULL); /* no overlap */ - ourspot.flen = theirspot.flen; + if (d*d > maxrad2) + ourspot.siz = 2.*PI*(1. - sqrt(1.-maxrad2/(d*d))); + else + nsflags &= ~SSPOT; + if (source[sn].sflags & SSPOT) { + copystruct(&theirspot, source[sn].sl.s); + multv3(theirspot.aim, source[sn].sl.s->aim, pm); + if (nsflags & SSPOT) { + ourspot.flen = theirspot.flen; + d = ourspot.siz; + if (!commonspot(&ourspot, &theirspot, nsloc)) + return(-1); /* no overlap */ + } else { + nsflags |= SSPOT; + copystruct(&ourspot, &theirspot); + d = 2.*ourspot.siz; + } + if (ourspot.siz < d-FTINY) { /* it shrunk */ + d = spotdisk(v, op, &ourspot, nsloc); + if (d <= FTINY) + return(-1); + if (d < maxrad2) { + maxrad2 = d; + VCOPY(ocent, v); + } + } } - if (sp->sflags & SFLAT) { /* check for behind source */ - multv3(nsnorm, sp->snorm, pm); - if (checkspot(&ourspot, nsnorm) < 0) - return(NULL); + if (source[sn].sflags & SFLAT) { /* behind source? */ + multv3(nsnorm, source[sn].snorm, pm); + if (!checkspot(&ourspot, nsnorm)) + return(-1); } } - /* everything is OK, make source */ - if ((newsrc = newsource()) == NULL) + /* pretest visibility */ + nsflags = vstestvis(nsflags, op, ocent, maxrad2, sn); + if (nsflags & SSKIP) + return(-1); /* obstructed */ + /* it all checks out, so make it */ + if ((i = newsource()) < 0) goto memerr; - newsrc->sflags = sp->sflags | (SVIRTUAL|SSPOT|SFOLLOW); - VCOPY(newsrc->sloc, nsloc); - if (newsrc->sflags & SFLAT) - VCOPY(newsrc->snorm, nsnorm); - newsrc->ss = sp->ss; newsrc->ss2 = sp->ss2; - if ((newsrc->sl.s = (SPOT *)malloc(sizeof(SPOT))) == NULL) - goto memerr; - copystruct(newsrc->sl.s, &ourspot); - if (newsrc->sflags & SPROX) - newsrc->sl.prox = sp->sl.prox; - newsrc->sa.svnext = sp - source; - return(newsrc); + source[i].sflags = nsflags; + VCOPY(source[i].sloc, nsloc); + if (nsflags & SFLAT) + VCOPY(source[i].snorm, nsnorm); + source[i].ss = source[sn].ss; source[i].ss2 = source[sn].ss2; + if (nsflags & SSPOT) { + if ((source[i].sl.s = (SPOT *)malloc(sizeof(SPOT))) == NULL) + goto memerr; + copystruct(source[i].sl.s, &ourspot); + } + if (nsflags & SPROX) + source[i].sl.prox = source[sn].sl.prox; + source[i].sa.svnext = sn; + source[i].so = op; + return(i); memerr: error(SYSTEM, "out of memory in makevsrc"); } -commonspot(sp1, sp2, org) /* set sp1 to intersection of sp1 and sp2 */ -register SPOT *sp1, *sp2; -FVECT org; +double +getdisk(oc, op, sn) /* get visible object disk */ +FVECT oc; +OBJREC *op; +register int sn; { - FVECT cent; - double rad2, cos1, cos2; - - cos1 = 1. - sp1->siz/(2.*PI); - cos2 = 1. - sp2->siz/(2.*PI); - if (sp2->siz >= 2.*PI-FTINY) /* BIG, just check overlap */ - return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 - - sqrt((1.-cos1*cos1)*(1.-cos2*cos2))); - /* compute and check disks */ - rad2 = intercircle(cent, sp1->aim, sp2->aim, - 1./(cos1*cos1) - 1., 1./(cos2*cos2) - 1.); - if (rad2 <= FTINY || normalize(cent) == 0.) - return(0); - VCOPY(sp1->aim, cent); - sp1->siz = 2.*PI*(1. - 1./sqrt(1.+rad2)); - return(1); + double rad2, roffs, offs, d, rd, rdoto; + FVECT rnrm, nrm; + /* first, use object getdisk function */ + rad2 = getmaxdisk(oc, op); + if (!(source[sn].sflags & SVIRTUAL)) + return(rad2); /* all done for normal source */ + /* check for correct side of relay surface */ + roffs = getplaneq(rnrm, source[sn].so); + rd = DOT(rnrm, source[sn].sloc); /* source projection */ + if (!(source[sn].sflags & SDISTANT)) + rd -= roffs; + d = DOT(rnrm, oc) - roffs; /* disk distance to relay plane */ + if ((d > 0.) ^ (rd > 0.)) + return(rad2); /* OK if opposite sides */ + if (d*d >= rad2) + return(0.); /* no relay is possible */ + /* we need a closer look */ + offs = getplaneq(nrm, op); + rdoto = DOT(rnrm, nrm); + if (d*d >= rad2*(1.-rdoto*rdoto)) + return(0.); /* disk entirely on projection side */ + /* should shrink disk but I'm lazy */ + return(rad2); } -commonbeam(sp1, sp2, dir) /* set sp1 to intersection of sp1 and sp2 */ -register SPOT *sp1, *sp2; -FVECT dir; +int +vstestvis(f, o, oc, or2, sn) /* pretest source visibility */ +int f; /* virtual source flags */ +OBJREC *o; /* relay object */ +FVECT oc; /* relay object center */ +double or2; /* relay object radius squared */ +register int sn; /* target source number */ { - FVECT cent, c1, c2; - double rad2, d; - register int i; - /* move centers to common plane */ - d = DOT(sp1->aim, dir); - for (i = 0; i < 3; i++) - c1[i] = sp1->aim[i] - d*dir[i]; - d = DOT(sp2->aim, dir); - for (i = 0; i < 3; i++) - c2[i] = sp2->aim[i] - d*dir[i]; - /* compute overlap */ - rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI); - if (rad2 <= FTINY) - return(0); - VCOPY(sp1->aim, cent); - sp1->siz = PI*rad2; - return(1); -} - - -checkspot(sp, nrm) /* check spotlight for behind source */ -register SPOT *sp; -FVECT nrm; -{ - double d, d1; - - d = DOT(sp->aim, nrm); - if (d > FTINY) /* center in front? */ - return(0); - /* else check horizon */ - d1 = 1. - sp->siz/(2.*PI); - return(1.-FTINY-d*d > d1*d1); -} - - -mirrorproj(m, nv, offs) /* get mirror projection for surface */ -register MAT4 m; -FVECT nv; -double offs; -{ - register int i, j; - /* assign matrix */ - setident4(m); - for (i = 0; i < 3; i++) - for (j = 0; j < 3; j++) - m[i][j] -= 2.*nv[i]*nv[j]; - for (j = 0; j < 3; j++) - m[3][j] = 2.*offs*nv[j]; -} - - -double -intercircle(cc, c1, c2, r1s, r2s) /* intersect two circles */ -FVECT cc; /* midpoint (return value) */ -FVECT c1, c2; /* circle centers */ -double r1s, r2s; /* radii squared */ -{ - double a2, d2, l; - FVECT disp; - register int i; - - for (i = 0; i < 3; i++) - disp[i] = c2[i] - c1[i]; - d2 = DOT(disp,disp); - /* circle within overlap? */ - if (r1s < r2s) { - if (r2s >= r1s + d2) { - VCOPY(cc, c1); - return(r1s); - } + RAY sr; + FVECT onorm; + FVECT offsdir; + double or, d; + int infront; + int stestlim, ssn; + int nhit, nok; + register int i, n; + /* return if pretesting disabled */ + if (vspretest <= 0) + return(f); + /* get surface normal */ + getplaneq(onorm, o); + /* set number of rays to sample */ + if (source[sn].sflags & SDISTANT) { + n = (2./3.*PI*PI)*or2/(thescene.cusize*thescene.cusize)* + vspretest + .5; + infront = DOT(onorm, source[sn].sloc) > 0.; } else { - if (r1s >= r2s + d2) { - VCOPY(cc, c2); - return(r2s); - } + for (i = 0; i < 3; i++) + offsdir[i] = source[sn].sloc[i] - oc[i]; + n = or2/DOT(offsdir,offsdir)*vspretest + .5; + infront = DOT(onorm, offsdir) > 0.; } - a2 = .25*(2.*(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2); - /* no overlap? */ - if (a2 <= 0.) - return(0.); - /* overlap, compute center */ - l = sqrt((r1s - a2)/d2); - for (i = 0; i < 3; i++) - cc[i] = c1[i] + l*disp[i]; - return(a2); -} - - -/* - * The following routines depend on the supported OBJECTS: - */ - - -double -getmaxdisk(ocent, op) /* get object center and squared radius */ -FVECT ocent; -register OBJREC *op; -{ - double maxrad2; - - switch (op->otype) { - case OBJ_FACE: - { - double d2; - register int i, j; - register FACE *f = getface(op); - - for (i = 0; i < 3; i++) { - ocent[i] = 0.; - for (j = 0; j < f->nv; j++) - ocent[i] += VERTEX(f,j)[i]; - ocent[i] /= (double)f->nv; + if (n < MINSAMPLES) n = MINSAMPLES; +#ifdef DEBUG + fprintf(stderr, "pretesting source %d in object %s with %d rays\n", + sn, o->oname, n); +#endif + /* sample */ + or = sqrt(or2); + stestlim = n*STESTMAX; + ssn = 0; + nhit = nok = 0; + while (n-- > 0) { + /* get sample point */ + do { + if (ssn >= stestlim) { +#ifdef DEBUG + fprintf(stderr, "\ttoo hard to hit\n"); +#endif + return(f); /* too small a target! */ } - maxrad2 = 0.; - for (j = 0; j < f->nv; j++) { - d2 = dist2(VERTEX(f,j), ocent); - if (d2 > maxrad2) - maxrad2 = d2; - } + for (i = 0; i < 3; i++) + offsdir[i] = or*(1. - + 2.*urand(urind(931*i+5827,ssn))); + ssn++; + for (i = 0; i < 3; i++) + sr.rorg[i] = oc[i] + offsdir[i]; + d = DOT(offsdir,onorm); + if (infront) + for (i = 0; i < 3; i++) { + sr.rorg[i] -= (d-.0001)*onorm[i]; + sr.rdir[i] = -onorm[i]; + } + else + for (i = 0; i < 3; i++) { + sr.rorg[i] -= (d+.0001)*onorm[i]; + sr.rdir[i] = onorm[i]; + } + rayorigin(&sr, NULL, PRIMARY, 1.0); + } while (!(*ofun[o->otype].funp)(o, &sr)); + /* check against source */ + samplendx++; + if (srcray(&sr, NULL, sn) == 0.) + continue; + sr.revf = srcvalue; + rayvalue(&sr); + if (bright(sr.rcol) <= FTINY) + continue; + nok++; + /* check against obstructions */ + srcray(&sr, NULL, sn); + rayvalue(&sr); + if (bright(sr.rcol) > FTINY) + nhit++; + if (nhit > 0 && nhit < nok) { +#ifdef DEBUG + fprintf(stderr, "\tpartially occluded\n"); +#endif + return(f); /* need to shadow test */ } - return(maxrad2); - case OBJ_RING: - { - register CONE *co = getcone(op, 0); - - VCOPY(ocent, CO_P0(co)); - maxrad2 = CO_R1(co); - maxrad2 *= maxrad2; - } - return(maxrad2); } - objerror(op, USER, "illegal material"); -} - - -double -getplaneq(nvec, op) /* get plane equation for object */ -FVECT nvec; -OBJREC *op; -{ - register FACE *fo; - register CONE *co; - - switch (op->otype) { - case OBJ_FACE: - fo = getface(op); - VCOPY(nvec, fo->norm); - return(fo->offset); - case OBJ_RING: - co = getcone(op, 0); - VCOPY(nvec, co->ad); - return(DOT(nvec, CO_P0(co))); - } - objerror(op, USER, "illegal material"); -} - - -/* - * The following routines depend on the supported MATERIALS: - */ - - -vproject(o, s, n) /* create projected source(s) if they exist */ -OBJREC *o; -SRCREC *s; -int n; -{ - SRCREC *ns; - FVECT norm; - double offset; - MAT4 proj; - /* get surface normal and offset */ - offset = getplaneq(norm, o); - switch (objptr(o->omod)->otype) { - case MAT_MIRROR: /* mirror source */ - if (DOT(s->sloc, norm) <= (s->sflags & SDISTANT ? - FTINY : offset+FTINY)) - return; /* behind mirror */ - mirrorproj(proj, norm, offset); - if ((ns = makevsrc(o, s, proj)) != NULL) - addvirtuals(ns, n); - break; - } -} - - -vsrcrelay(rn, rv) /* relay virtual source ray */ -register RAY *rn, *rv; -{ - int snext; - register int i; - /* source we're aiming for here */ - snext = source[rv->rsrc].sa.svnext; - /* compute relayed ray direction */ - switch (objptr(rv->ro->omod)->otype) { - case MAT_MIRROR: /* mirror: singular reflection */ - rayorigin(rn, rv, REFLECTED, 1.); - /* ignore textures */ - for (i = 0; i < 3; i++) - rn->rdir[i] = rv->rdir[i] + 2.*rv->rod*rv->ron[i]; - break; + if (nhit == 0) { #ifdef DEBUG - default: - error(CONSISTENCY, "inappropriate material in vsrcrelay"); + fprintf(stderr, "\t0%% hit rate\n"); #endif + return(f | SSKIP); /* 0% hit rate: totally occluded */ } - rn->rsrc = snext; +#ifdef DEBUG + fprintf(stderr, "\t100%% hit rate\n"); +#endif + return(f & ~SFOLLOW); /* 100% hit rate: no occlusion */ } + - -m_mirror(m, r) /* shade mirrored ray */ -register OBJREC *m; -register RAY *r; +#ifdef DEBUG +virtverb(sn, fp) /* print verbose description of virtual source */ +register int sn; +FILE *fp; { - COLOR mcolor; - RAY nr; register int i; - if (m->oargs.nfargs != 3 || m->oargs.nsargs > 1) - objerror(m, USER, "bad number of arguments"); - if (r->rsrc >= 0) { /* aiming for somebody */ - if (source[r->rsrc].so != r->ro) - return; /* but not us */ - } else if (m->oargs.nsargs > 0) { /* else call substitute? */ - rayshade(r, modifier(m->oargs.sarg[0])); + fprintf(fp, "%s virtual source %d in %s %s\n", + source[sn].sflags & SDISTANT ? "distant" : "local", + sn, ofun[source[sn].so->otype].funame, + source[sn].so->oname); + fprintf(fp, "\tat (%f,%f,%f)\n", + source[sn].sloc[0], source[sn].sloc[1], source[sn].sloc[2]); + fprintf(fp, "\tlinked to source %d (%s)\n", + source[sn].sa.svnext, source[source[sn].sa.svnext].so->oname); + if (source[sn].sflags & SFOLLOW) + fprintf(fp, "\talways followed\n"); + else + fprintf(fp, "\tnever followed\n"); + if (!(source[sn].sflags & SSPOT)) return; - } - if (r->rod < 0.) /* back is black */ - return; - /* get modifiers */ - raytexture(r, m->omod); - /* assign material color */ - setcolor(mcolor, m->oargs.farg[0], - m->oargs.farg[1], - m->oargs.farg[2]); - multcolor(mcolor, r->pcol); - /* compute reflected ray */ - if (r->rsrc >= 0) /* relayed light source */ - vsrcrelay(&nr, r); - else { /* ordinary reflection */ - FVECT pnorm; - double pdot; - - if (rayorigin(&nr, r, REFLECTED, bright(mcolor)) < 0) - return; - pdot = raynormal(pnorm, r); /* use textures */ - for (i = 0; i < 3; i++) - nr.rdir[i] = r->rdir[i] + 2.*pdot*pnorm[i]; - } - rayvalue(&nr); - multcolor(nr.rcol, mcolor); - addcolor(r->rcol, nr.rcol); + fprintf(fp, "\twith spot aim (%f,%f,%f) and size %f\n", + source[sn].sl.s->aim[0], source[sn].sl.s->aim[1], + source[sn].sl.s->aim[2], source[sn].sl.s->siz); } +#endif