src/rt/srcsupp.c

#ifndef lint
static const char       RCSid[] = "$Id: srcsupp.c,v 2.23 2016/04/21 00:40:35 greg Exp $";
#endif
/*
 *  Support routines for source objects and materials
 *
 *  External symbols declared in source.h
 */

#include "copyright.h"

#include  "ray.h"

#include  "otypes.h"

#include  "source.h"

#include  "cone.h"

#include  "face.h"

#define SRCINC          32              /* realloc increment for array */

SRCREC  *source = NULL;                 /* our list of sources */
int  nsources = 0;                      /* the number of sources */

SRCFUNC  sfun[NUMOTYPE];                /* source dispatch table */


void
initstypes(void)                        /* initialize source dispatch table */
{
        extern VSMATERIAL  mirror_vs, direct1_vs, direct2_vs;
        static SOBJECT  fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
        static SOBJECT  ssobj = {ssetsrc, nopart};
        static SOBJECT  sphsobj = {sphsetsrc, nopart};
        static SOBJECT  cylsobj = {cylsetsrc, cylpart};
        static SOBJECT  rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};

        sfun[MAT_MIRROR].mf = &mirror_vs;
        sfun[MAT_DIRECT1].mf = &direct1_vs;
        sfun[MAT_DIRECT2].mf = &direct2_vs;
        sfun[OBJ_FACE].of = &fsobj;
        sfun[OBJ_SOURCE].of = &ssobj;
        sfun[OBJ_SPHERE].of = &sphsobj;
        sfun[OBJ_CYLINDER].of = &cylsobj;
        sfun[OBJ_RING].of = &rsobj;
}


int
newsource(void)                 /* allocate new source in our array */
{
        if (nsources == 0)
                source = (SRCREC *)malloc(SRCINC*sizeof(SRCREC));
        else if (nsources%SRCINC == 0)
                source = (SRCREC *)realloc((void *)source,
                                (unsigned)(nsources+SRCINC)*sizeof(SRCREC));
        if (source == NULL)
                return(-1);
        source[nsources].sflags = 0;
        source[nsources].nhits = 1;
        source[nsources].ntests = 2;    /* initial hit probability = 50% */
#if SHADCACHE
        source[nsources].obscache = NULL;
#endif
        return(nsources++);
}


void
setflatss(                              /* set sampling for a flat source */
        SRCREC  *src
)
{
        double  mult;
        int  i;

        getperpendicular(src->ss[SU], src->snorm, rand_samp);
        mult = .5 * sqrt( src->ss2 );
        for (i = 0; i < 3; i++)
                src->ss[SU][i] *= mult;
        fcross(src->ss[SV], src->snorm, src->ss[SU]);
}


void
fsetsrc(                        /* set a face as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        FACE  *f;
        int  i, j;
        double  d;
        
        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
                                                /* get the face */
        f = getface(so);
        if (f->area == 0.)
                objerror(so, USER, "zero source area");
                                                /* find the center */
        for (j = 0; j < 3; j++) {
                src->sloc[j] = 0.0;
                for (i = 0; i < f->nv; i++)
                        src->sloc[j] += VERTEX(f,i)[j];
                src->sloc[j] /= (double)f->nv;
        }
        if (!inface(src->sloc, f))
                objerror(so, USER, "cannot hit source center");
        src->sflags |= SFLAT;
        VCOPY(src->snorm, f->norm);
        src->ss2 = f->area;
                                                /* find maximum radius */
        src->srad = 0.;
        for (i = 0; i < f->nv; i++) {
                d = dist2(VERTEX(f,i), src->sloc);
                if (d > src->srad)
                        src->srad = d;
        }
        src->srad = sqrt(src->srad);
                                                /* compute size vectors */
        if (f->nv == 4) {                               /* parallelogram case */
                for (j = 0; j < 3; j++) {
                        src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
                        src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
                }
        } else if (f->nv == 3) {                        /* triangle case */
                int     near0 = 2;
                double  dmin = dist2line(src->sloc, VERTEX(f,2), VERTEX(f,0));
                for (i = 0; i < 2; i++) {
                        double  d2 = dist2line(src->sloc, VERTEX(f,i), VERTEX(f,i+1));
                        if (d2 >= dmin)
                                continue;
                        near0 = i;
                        dmin = d2;                      /* radius = min distance */
                }
                if (dmin < .08*f->area)
                        objerror(so, WARNING, "triangular source with poor aspect");
                i = (near0 + 1) % 3;
                for (j = 0; j < 3; j++)
                        src->ss[SU][j] = VERTEX(f,i)[j] - VERTEX(f,near0)[j];
                normalize(src->ss[SU]);
                dmin = sqrt(dmin);
                for (j = 0; j < 3; j++)
                        src->ss[SU][j] *= dmin;
                fcross(src->ss[SV], f->norm, src->ss[SU]);
        } else
                setflatss(src);                         /* hope for convex! */
}


void
ssetsrc(                        /* set a source as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        double  theta;
        
        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
        if (so->oargs.nfargs != 4)
                objerror(so, USER, "bad arguments");
        src->sflags |= (SDISTANT|SCIR);
        VCOPY(src->sloc, so->oargs.farg);
        if (normalize(src->sloc) == 0.0)
                objerror(so, USER, "zero direction");
        theta = PI/180.0/2.0 * so->oargs.farg[3];
        if (theta <= FTINY)
                objerror(so, USER, "zero size");
        src->ss2 = 2.0*PI * (1.0 - cos(theta));
                                        /* the following is approximate */
        src->srad = sqrt(src->ss2/PI);
        VCOPY(src->snorm, src->sloc);
        setflatss(src);                 /* hey, whatever works */
}


void
sphsetsrc(                      /* set a sphere as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        int  i;

        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
        if (so->oargs.nfargs != 4)
                objerror(so, USER, "bad # arguments");
        if (so->oargs.farg[3] <= FTINY)
                objerror(so, USER, "illegal source radius");
        src->sflags |= SCIR;
        VCOPY(src->sloc, so->oargs.farg);
        src->srad = so->oargs.farg[3];
        src->ss2 = PI * src->srad * src->srad;
        memset(src->ss, 0, sizeof(src->ss));
        for (i = 0; i < 3; i++)
                src->ss[i][i] = 0.7236 * so->oargs.farg[3];
}


void
rsetsrc(                        /* set a ring (disk) as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        CONE  *co;
        
        src->sa.success = 2*AIMREQT-1;          /* bitch on second failure */
        src->so = so;
                                                /* get the ring */
        co = getcone(so, 0);
        if (co == NULL)
                objerror(so, USER, "illegal source");
        if (CO_R1(co) <= FTINY)
                objerror(so, USER, "illegal source radius");
        VCOPY(src->sloc, CO_P0(co));
        if (CO_R0(co) > 0.0)
                objerror(so, USER, "cannot hit source center");
        src->sflags |= (SFLAT|SCIR);
        VCOPY(src->snorm, co->ad);
        src->srad = CO_R1(co);
        src->ss2 = PI * src->srad * src->srad;
        setflatss(src);
}


void
cylsetsrc(                      /* set a cylinder as a source */
        SRCREC  *src,
        OBJREC  *so
)
{
        CONE  *co;
        int  i;
        
        src->sa.success = 4*AIMREQT-1;          /* bitch on fourth failure */
        src->so = so;
                                                /* get the cylinder */
        co = getcone(so, 0);
        if (co == NULL)
                objerror(so, USER, "illegal source");
        if (CO_R0(co) <= FTINY)
                objerror(so, USER, "illegal source radius");
        if (CO_R0(co) > .2*co->al)              /* heuristic constraint */
                objerror(so, WARNING, "source aspect too small");
        src->sflags |= SCYL;
        for (i = 0; i < 3; i++)
                src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
        src->srad = .5*co->al;
        src->ss2 = 2.*CO_R0(co)*co->al;
                                                /* set sampling vectors */
        for (i = 0; i < 3; i++)
                src->ss[SU][i] = .5 * co->al * co->ad[i];
        getperpendicular(src->ss[SW], co->ad, rand_samp);
        for (i = 0; i < 3; i++)
                src->ss[SW][i] *= .8559 * CO_R0(co);
        fcross(src->ss[SV], src->ss[SW], co->ad);
}


SPOT *
makespot(                       /* make a spotlight */
        OBJREC  *m
)
{
        SPOT  *ns;

        if ((ns = (SPOT *)m->os) != NULL)
                return(ns);
        if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
                return(NULL);
        if (m->oargs.farg[3] <= FTINY)
                objerror(m, USER, "zero angle");
        ns->siz = 2.0*PI * (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
        VCOPY(ns->aim, m->oargs.farg+4);
        if ((ns->flen = normalize(ns->aim)) == 0.0)
                objerror(m, USER, "zero focus vector");
        m->os = (char *)ns;
        return(ns);
}


int
spotout(                        /* check if we're outside spot region */
        RAY  *r,
        SPOT  *s
)
{
        double  d;
        FVECT  vd;
        
        if (s == NULL)
                return(0);
        if (s->flen < -FTINY) {         /* distant source */
                vd[0] = s->aim[0] - r->rorg[0];
                vd[1] = s->aim[1] - r->rorg[1];
                vd[2] = s->aim[2] - r->rorg[2];
                d = DOT(r->rdir,vd);
                /*                      wrong side?
                if (d <= FTINY)
                        return(1);      */
                d = DOT(vd,vd) - d*d;
                if (PI*d > s->siz)
                        return(1);      /* out */
                return(0);      /* OK */
        }
                                        /* local source */
        if (s->siz < 2.0*PI * (1.0 + DOT(s->aim,r->rdir)))
                return(1);      /* out */
        return(0);      /* OK */
}


double
fgetmaxdisk(            /* get center and squared radius of face */
        FVECT  ocent,
        OBJREC  *op
)
{
        double  maxrad2;
        double  d;
        int  i, j;
        FACE  *f;
        
        f = getface(op);
        if (f->area == 0.)
                return(0.);
        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;
        }
        d = DOT(ocent,f->norm);
        for (i = 0; i < 3; i++)
                ocent[i] += (f->offset - d)*f->norm[i];
        maxrad2 = 0.;
        for (j = 0; j < f->nv; j++) {
                d = dist2(VERTEX(f,j), ocent);
                if (d > maxrad2)
                        maxrad2 = d;
        }
        return(maxrad2);
}


double
rgetmaxdisk(            /* get center and squared radius of ring */
        FVECT  ocent,
        OBJREC  *op
)
{
        CONE  *co;
        
        co = getcone(op, 0);
        if (co == NULL)
                return(0.);
        VCOPY(ocent, CO_P0(co));
        return(CO_R1(co)*CO_R1(co));
}


double
fgetplaneq(                     /* get plane equation for face */
        FVECT  nvec,
        OBJREC  *op
)
{
        FACE  *fo;

        fo = getface(op);
        VCOPY(nvec, fo->norm);
        return(fo->offset);
}


double
rgetplaneq(                     /* get plane equation for ring */
        FVECT  nvec,
        OBJREC  *op
)
{
        CONE  *co;

        co = getcone(op, 0);
        if (co == NULL) {
                memset(nvec, 0, sizeof(FVECT));
                return(0.);
        }
        VCOPY(nvec, co->ad);
        return(DOT(nvec, CO_P0(co)));
}


int
commonspot(             /* set sp1 to intersection of sp1 and sp2 */
        SPOT  *sp1,
        SPOT  *sp2,
        FVECT  org
)
{
        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);
}


int
commonbeam(             /* set sp1 to intersection of sp1 and sp2 */
        SPOT  *sp1,
        SPOT  *sp2,
        FVECT  dir
)
{
        FVECT  cent, c1, c2;
        double  rad2, d;
                                        /* move centers to common plane */
        d = DOT(sp1->aim, dir);
        VSUM(c1, sp1->aim, dir, -d);
        d = DOT(sp2->aim, dir);
        VSUM(c2, sp2->aim, dir, -d);
                                        /* 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);
}


int
checkspot(                      /* check spotlight for behind source */
        SPOT  *sp,      /* spotlight */
        FVECT  nrm      /* source surface normal */
)
{
        double  d, d1;

        d = DOT(sp->aim, nrm);
        if (d > FTINY)                  /* center in front? */
                return(1);
                                        /* else check horizon */
        d1 = 1. - sp->siz/(2.*PI);
        return(1.-FTINY-d*d < d1*d1);
}


double
spotdisk(               /* intersect spot with object op */
        FVECT  oc,
        OBJREC  *op,
        SPOT  *sp,
        FVECT  pos
)
{
        FVECT  onorm;
        double  offs, d, dist;

        offs = getplaneq(onorm, op);
        d = -DOT(onorm, sp->aim);
        if (d >= -FTINY && d <= FTINY)
                return(0.);
        dist = (DOT(pos, onorm) - offs)/d;
        if (dist < 0.)
                return(0.);
        VSUM(oc, pos, sp->aim, dist);
        return(sp->siz*dist*dist/PI/(d*d));
}


double
beamdisk(               /* intersect beam with object op */
        FVECT  oc,
        OBJREC  *op,
        SPOT  *sp,
        FVECT  dir
)
{
        FVECT  onorm;
        double  offs, d, dist;

        offs = getplaneq(onorm, op);
        d = -DOT(onorm, dir);
        if (d >= -FTINY && d <= FTINY)
                return(0.);
        dist = (DOT(sp->aim, onorm) - offs)/d;
        VSUM(oc, sp->aim, dir, dist);
        return(sp->siz/PI/(d*d));
}


double
intercircle(                            /* intersect two circles */
        FVECT  cc,              /* midpoint (return value) */
        FVECT  c1,              /* circle centers */
        FVECT  c2,
        double  r1s,            /* radii squared */
        double  r2s
)
{
        double  a2, d2, l;
        FVECT  disp;

        VSUB(disp, c2, c1);
        d2 = DOT(disp,disp);
                                        /* circle within overlap? */
        if (r1s < r2s) {
                if (r2s >= r1s + d2) {
                        VCOPY(cc, c1);
                        return(r1s);
                }
        } else {
                if (r1s >= r2s + d2) {
                        VCOPY(cc, c2);
                        return(r2s);
                }
        }
        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);
        VSUM(cc, c1, disp, l);
        return(a2);
}
Revision:	2.24
Committed:	Thu Aug 4 22:43:46 2022 UTC (21 months, 1 week ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, HEAD
Changes since 2.23:	+25 -7 lines
Log Message:	fix: reduced aiming failures for triangular light sources
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	2.24	static const char RCSid[] = "$Id: srcsupp.c,v 2.23 2016/04/21 00:40:35 greg Exp $";
3	greg	1.1	#endif
4			/*
5			* Support routines for source objects and materials
6	greg	2.9	*
7			* External symbols declared in source.h
8			*/
9
10	greg	2.10	#include "copyright.h"
11	greg	1.1
12			#include "ray.h"
13
14			#include "otypes.h"
15
16			#include "source.h"
17
18			#include "cone.h"
19
20			#include "face.h"
21
22	greg	2.19	#define SRCINC 32 /* realloc increment for array */
23	greg	1.1
24			SRCREC source = NULL; / our list of sources */
25			int nsources = 0; /* the number of sources */
26
27			SRCFUNC sfun[NUMOTYPE]; /* source dispatch table */
28
29
30	greg	2.9	void
31	greg	2.20	initstypes(void) /* initialize source dispatch table */
32	greg	1.1	{
33	greg	1.9	extern VSMATERIAL mirror_vs, direct1_vs, direct2_vs;
34	greg	1.14	static SOBJECT fsobj = {fsetsrc, flatpart, fgetplaneq, fgetmaxdisk};
35			static SOBJECT ssobj = {ssetsrc, nopart};
36			static SOBJECT sphsobj = {sphsetsrc, nopart};
37			static SOBJECT cylsobj = {cylsetsrc, cylpart};
38			static SOBJECT rsobj = {rsetsrc, flatpart, rgetplaneq, rgetmaxdisk};
39	greg	1.1
40			sfun[MAT_MIRROR].mf = &mirror_vs;
41	greg	1.9	sfun[MAT_DIRECT1].mf = &direct1_vs;
42			sfun[MAT_DIRECT2].mf = &direct2_vs;
43	greg	1.1	sfun[OBJ_FACE].of = &fsobj;
44			sfun[OBJ_SOURCE].of = &ssobj;
45			sfun[OBJ_SPHERE].of = &sphsobj;
46	greg	1.14	sfun[OBJ_CYLINDER].of = &cylsobj;
47	greg	1.1	sfun[OBJ_RING].of = &rsobj;
48			}
49
50
51	greg	1.2	int
52	greg	2.20	newsource(void) /* allocate new source in our array */
53	greg	1.1	{
54			if (nsources == 0)
55	greg	1.13	source = (SRCREC )malloc(SRCINCsizeof(SRCREC));
56			else if (nsources%SRCINC == 0)
57	greg	2.11	source = (SRCREC )realloc((void )source,
58	greg	1.13	(unsigned)(nsources+SRCINC)*sizeof(SRCREC));
59	greg	1.1	if (source == NULL)
60	greg	1.2	return(-1);
61	greg	1.1	source[nsources].sflags = 0;
62	greg	2.14	source[nsources].nhits = 1;
63			source[nsources].ntests = 2; /* initial hit probability = 50% */
64	greg	2.12	#if SHADCACHE
65			source[nsources].obscache = NULL;
66			#endif
67	greg	1.2	return(nsources++);
68	greg	1.1	}
69
70
71	greg	2.9	void
72	greg	2.20	setflatss( /* set sampling for a flat source */
73			SRCREC *src
74			)
75	greg	1.14	{
76			double mult;
77	greg	2.20	int i;
78	greg	1.14
79	greg	2.22	getperpendicular(src->ss[SU], src->snorm, rand_samp);
80	greg	2.21	mult = .5 * sqrt( src->ss2 );
81	greg	1.14	for (i = 0; i < 3; i++)
82			src->ss[SU][i] *= mult;
83			fcross(src->ss[SV], src->snorm, src->ss[SU]);
84			}
85
86
87	greg	2.9	void
88	greg	2.20	fsetsrc( /* set a face as a source */
89			SRCREC *src,
90			OBJREC *so
91			)
92	greg	1.1	{
93	greg	2.20	FACE *f;
94			int i, j;
95	greg	1.14	double d;
96	greg	1.1
97			src->sa.success = 2AIMREQT-1; / bitch on second failure */
98			src->so = so;
99			/* get the face */
100			f = getface(so);
101	greg	2.23	if (f->area == 0.)
102	greg	2.16	objerror(so, USER, "zero source area");
103	greg	1.1	/* find the center */
104			for (j = 0; j < 3; j++) {
105			src->sloc[j] = 0.0;
106			for (i = 0; i < f->nv; i++)
107			src->sloc[j] += VERTEX(f,i)[j];
108			src->sloc[j] /= (double)f->nv;
109			}
110			if (!inface(src->sloc, f))
111	greg	2.16	objerror(so, USER, "cannot hit source center");
112	greg	1.1	src->sflags \|= SFLAT;
113			VCOPY(src->snorm, f->norm);
114			src->ss2 = f->area;
115	greg	1.14	/* find maximum radius */
116			src->srad = 0.;
117			for (i = 0; i < f->nv; i++) {
118			d = dist2(VERTEX(f,i), src->sloc);
119			if (d > src->srad)
120			src->srad = d;
121			}
122			src->srad = sqrt(src->srad);
123			/* compute size vectors */
124	greg	2.24	if (f->nv == 4) { /* parallelogram case */
125	greg	1.14	for (j = 0; j < 3; j++) {
126			src->ss[SU][j] = .5*(VERTEX(f,1)[j]-VERTEX(f,0)[j]);
127			src->ss[SV][j] = .5*(VERTEX(f,3)[j]-VERTEX(f,0)[j]);
128			}
129	greg	2.24	} else if (f->nv == 3) { /* triangle case */
130			int near0 = 2;
131			double dmin = dist2line(src->sloc, VERTEX(f,2), VERTEX(f,0));
132			for (i = 0; i < 2; i++) {
133			double d2 = dist2line(src->sloc, VERTEX(f,i), VERTEX(f,i+1));
134			if (d2 >= dmin)
135			continue;
136			near0 = i;
137			dmin = d2; /* radius = min distance */
138			}
139			if (dmin < .08*f->area)
140			objerror(so, WARNING, "triangular source with poor aspect");
141			i = (near0 + 1) % 3;
142			for (j = 0; j < 3; j++)
143			src->ss[SU][j] = VERTEX(f,i)[j] - VERTEX(f,near0)[j];
144			normalize(src->ss[SU]);
145			dmin = sqrt(dmin);
146			for (j = 0; j < 3; j++)
147			src->ss[SU][j] *= dmin;
148			fcross(src->ss[SV], f->norm, src->ss[SU]);
149			} else
150			setflatss(src); /* hope for convex! */
151	greg	1.1	}
152
153
154	greg	2.9	void
155	greg	2.20	ssetsrc( /* set a source as a source */
156			SRCREC *src,
157			OBJREC *so
158			)
159	greg	1.1	{
160			double theta;
161
162			src->sa.success = 2AIMREQT-1; / bitch on second failure */
163			src->so = so;
164			if (so->oargs.nfargs != 4)
165			objerror(so, USER, "bad arguments");
166	greg	2.17	src->sflags \|= (SDISTANT\|SCIR);
167	greg	1.1	VCOPY(src->sloc, so->oargs.farg);
168			if (normalize(src->sloc) == 0.0)
169			objerror(so, USER, "zero direction");
170			theta = PI/180.0/2.0 * so->oargs.farg[3];
171			if (theta <= FTINY)
172			objerror(so, USER, "zero size");
173			src->ss2 = 2.0PI (1.0 - cos(theta));
174	greg	1.14	/* the following is approximate */
175			src->srad = sqrt(src->ss2/PI);
176			VCOPY(src->snorm, src->sloc);
177			setflatss(src); /* hey, whatever works */
178	greg	1.1	}
179
180
181	greg	2.9	void
182	greg	2.20	sphsetsrc( /* set a sphere as a source */
183			SRCREC *src,
184			OBJREC *so
185			)
186	greg	1.1	{
187	greg	2.20	int i;
188	greg	1.14
189	greg	1.1	src->sa.success = 2AIMREQT-1; / bitch on second failure */
190			src->so = so;
191			if (so->oargs.nfargs != 4)
192			objerror(so, USER, "bad # arguments");
193			if (so->oargs.farg[3] <= FTINY)
194	greg	2.16	objerror(so, USER, "illegal source radius");
195	greg	2.17	src->sflags \|= SCIR;
196	greg	1.1	VCOPY(src->sloc, so->oargs.farg);
197	greg	1.14	src->srad = so->oargs.farg[3];
198			src->ss2 = PI * src->srad * src->srad;
199	greg	2.24	memset(src->ss, 0, sizeof(src->ss));
200	greg	1.14	for (i = 0; i < 3; i++)
201	greg	2.17	src->ss[i][i] = 0.7236 * so->oargs.farg[3];
202	greg	1.1	}
203
204
205	greg	2.9	void
206	greg	2.20	rsetsrc( /* set a ring (disk) as a source */
207			SRCREC *src,
208			OBJREC *so
209			)
210	greg	1.1	{
211	greg	2.20	CONE *co;
212	greg	1.1
213			src->sa.success = 2AIMREQT-1; / bitch on second failure */
214			src->so = so;
215			/* get the ring */
216			co = getcone(so, 0);
217	greg	2.23	if (co == NULL)
218			objerror(so, USER, "illegal source");
219	greg	2.16	if (CO_R1(co) <= FTINY)
220			objerror(so, USER, "illegal source radius");
221	greg	1.1	VCOPY(src->sloc, CO_P0(co));
222			if (CO_R0(co) > 0.0)
223	greg	2.16	objerror(so, USER, "cannot hit source center");
224	greg	2.17	src->sflags \|= (SFLAT\|SCIR);
225	greg	1.1	VCOPY(src->snorm, co->ad);
226	greg	1.14	src->srad = CO_R1(co);
227			src->ss2 = PI * src->srad * src->srad;
228			setflatss(src);
229			}
230
231
232	greg	2.9	void
233	greg	2.20	cylsetsrc( /* set a cylinder as a source */
234			SRCREC *src,
235			OBJREC *so
236			)
237	greg	1.14	{
238	greg	2.20	CONE *co;
239			int i;
240	greg	1.14
241			src->sa.success = 4AIMREQT-1; / bitch on fourth failure */
242			src->so = so;
243			/* get the cylinder */
244			co = getcone(so, 0);
245	greg	2.23	if (co == NULL)
246			objerror(so, USER, "illegal source");
247	greg	2.16	if (CO_R0(co) <= FTINY)
248			objerror(so, USER, "illegal source radius");
249	greg	1.14	if (CO_R0(co) > .2co->al) / heuristic constraint */
250			objerror(so, WARNING, "source aspect too small");
251	greg	1.15	src->sflags \|= SCYL;
252	greg	1.14	for (i = 0; i < 3; i++)
253			src->sloc[i] = .5 * (CO_P1(co)[i] + CO_P0(co)[i]);
254	greg	1.15	src->srad = .5*co->al;
255	greg	1.14	src->ss2 = 2.CO_R0(co)co->al;
256			/* set sampling vectors */
257			for (i = 0; i < 3; i++)
258			src->ss[SU][i] = .5 * co->al * co->ad[i];
259	greg	2.22	getperpendicular(src->ss[SW], co->ad, rand_samp);
260	greg	1.14	for (i = 0; i < 3; i++)
261	greg	1.15	src->ss[SW][i] = .8559 CO_R0(co);
262	greg	1.14	fcross(src->ss[SV], src->ss[SW], co->ad);
263	greg	1.1	}
264
265
266			SPOT *
267	greg	2.20	makespot( /* make a spotlight */
268			OBJREC *m
269			)
270	greg	1.1	{
271	greg	2.20	SPOT *ns;
272	greg	1.1
273	greg	2.5	if ((ns = (SPOT *)m->os) != NULL)
274			return(ns);
275	greg	1.1	if ((ns = (SPOT *)malloc(sizeof(SPOT))) == NULL)
276			return(NULL);
277	greg	2.16	if (m->oargs.farg[3] <= FTINY)
278			objerror(m, USER, "zero angle");
279	greg	1.1	ns->siz = 2.0PI (1.0 - cos(PI/180.0/2.0 * m->oargs.farg[3]));
280			VCOPY(ns->aim, m->oargs.farg+4);
281			if ((ns->flen = normalize(ns->aim)) == 0.0)
282			objerror(m, USER, "zero focus vector");
283	greg	2.5	m->os = (char *)ns;
284	greg	1.1	return(ns);
285			}
286
287
288	greg	2.9	int
289	greg	2.20	spotout( /* check if we're outside spot region */
290			RAY *r,
291			SPOT *s
292			)
293	greg	2.5	{
294			double d;
295			FVECT vd;
296
297			if (s == NULL)
298			return(0);
299	greg	2.8	if (s->flen < -FTINY) { /* distant source */
300	greg	2.5	vd[0] = s->aim[0] - r->rorg[0];
301			vd[1] = s->aim[1] - r->rorg[1];
302			vd[2] = s->aim[2] - r->rorg[2];
303			d = DOT(r->rdir,vd);
304			/* wrong side?
305			if (d <= FTINY)
306			return(1); */
307			d = DOT(vd,vd) - d*d;
308			if (PI*d > s->siz)
309			return(1); /* out */
310			return(0); /* OK */
311			}
312			/* local source */
313			if (s->siz < 2.0PI (1.0 + DOT(s->aim,r->rdir)))
314			return(1); /* out */
315			return(0); /* OK */
316			}
317
318
319	greg	1.1	double
320	greg	2.20	fgetmaxdisk( /* get center and squared radius of face */
321			FVECT ocent,
322			OBJREC *op
323			)
324	greg	1.1	{
325			double maxrad2;
326	greg	1.5	double d;
327	greg	2.20	int i, j;
328			FACE *f;
329	greg	1.1
330			f = getface(op);
331	greg	1.5	if (f->area == 0.)
332			return(0.);
333	greg	1.1	for (i = 0; i < 3; i++) {
334			ocent[i] = 0.;
335			for (j = 0; j < f->nv; j++)
336			ocent[i] += VERTEX(f,j)[i];
337			ocent[i] /= (double)f->nv;
338			}
339	greg	1.5	d = DOT(ocent,f->norm);
340			for (i = 0; i < 3; i++)
341			ocent[i] += (f->offset - d)*f->norm[i];
342	greg	1.1	maxrad2 = 0.;
343			for (j = 0; j < f->nv; j++) {
344	greg	1.5	d = dist2(VERTEX(f,j), ocent);
345			if (d > maxrad2)
346			maxrad2 = d;
347	greg	1.1	}
348			return(maxrad2);
349			}
350
351
352			double
353	greg	2.20	rgetmaxdisk( /* get center and squared radius of ring */
354			FVECT ocent,
355			OBJREC *op
356			)
357	greg	1.1	{
358	greg	2.20	CONE *co;
359	greg	1.1
360			co = getcone(op, 0);
361	greg	2.23	if (co == NULL)
362			return(0.);
363	greg	1.1	VCOPY(ocent, CO_P0(co));
364			return(CO_R1(co)*CO_R1(co));
365			}
366
367
368			double
369	greg	2.20	fgetplaneq( /* get plane equation for face */
370			FVECT nvec,
371			OBJREC *op
372			)
373	greg	1.1	{
374	greg	2.20	FACE *fo;
375	greg	1.1
376			fo = getface(op);
377			VCOPY(nvec, fo->norm);
378			return(fo->offset);
379			}
380
381
382			double
383	greg	2.20	rgetplaneq( /* get plane equation for ring */
384			FVECT nvec,
385			OBJREC *op
386			)
387	greg	1.1	{
388	greg	2.20	CONE *co;
389	greg	1.1
390			co = getcone(op, 0);
391	greg	2.23	if (co == NULL) {
392			memset(nvec, 0, sizeof(FVECT));
393			return(0.);
394			}
395	greg	1.1	VCOPY(nvec, co->ad);
396			return(DOT(nvec, CO_P0(co)));
397	greg	1.4	}
398
399
400	greg	2.9	int
401	greg	2.20	commonspot( /* set sp1 to intersection of sp1 and sp2 */
402			SPOT *sp1,
403			SPOT *sp2,
404			FVECT org
405			)
406	greg	1.4	{
407			FVECT cent;
408			double rad2, cos1, cos2;
409
410			cos1 = 1. - sp1->siz/(2.*PI);
411			cos2 = 1. - sp2->siz/(2.*PI);
412			if (sp2->siz >= 2.PI-FTINY) / BIG, just check overlap */
413			return(DOT(sp1->aim,sp2->aim) >= cos1*cos2 -
414			sqrt((1.-cos1cos1)(1.-cos2*cos2)));
415			/* compute and check disks */
416			rad2 = intercircle(cent, sp1->aim, sp2->aim,
417			1./(cos1cos1) - 1., 1./(cos2cos2) - 1.);
418			if (rad2 <= FTINY \|\| normalize(cent) == 0.)
419			return(0);
420			VCOPY(sp1->aim, cent);
421			sp1->siz = 2.PI(1. - 1./sqrt(1.+rad2));
422			return(1);
423			}
424
425
426	greg	2.9	int
427	greg	2.20	commonbeam( /* set sp1 to intersection of sp1 and sp2 */
428			SPOT *sp1,
429			SPOT *sp2,
430			FVECT dir
431			)
432	greg	1.4	{
433			FVECT cent, c1, c2;
434			double rad2, d;
435			/* move centers to common plane */
436			d = DOT(sp1->aim, dir);
437	greg	2.18	VSUM(c1, sp1->aim, dir, -d);
438	greg	1.4	d = DOT(sp2->aim, dir);
439	greg	2.18	VSUM(c2, sp2->aim, dir, -d);
440	greg	1.4	/* compute overlap */
441			rad2 = intercircle(cent, c1, c2, sp1->siz/PI, sp2->siz/PI);
442			if (rad2 <= FTINY)
443			return(0);
444			VCOPY(sp1->aim, cent);
445			sp1->siz = PI*rad2;
446			return(1);
447			}
448
449
450	greg	2.9	int
451	greg	2.20	checkspot( /* check spotlight for behind source */
452			SPOT sp, / spotlight */
453			FVECT nrm /* source surface normal */
454			)
455	greg	1.4	{
456			double d, d1;
457
458			d = DOT(sp->aim, nrm);
459			if (d > FTINY) /* center in front? */
460	greg	1.8	return(1);
461	greg	1.4	/* else check horizon */
462			d1 = 1. - sp->siz/(2.*PI);
463	greg	1.8	return(1.-FTINY-dd < d1d1);
464	greg	1.4	}
465
466
467			double
468	greg	2.20	spotdisk( /* intersect spot with object op */
469			FVECT oc,
470			OBJREC *op,
471			SPOT *sp,
472			FVECT pos
473			)
474	greg	1.6	{
475			FVECT onorm;
476			double offs, d, dist;
477
478			offs = getplaneq(onorm, op);
479			d = -DOT(onorm, sp->aim);
480			if (d >= -FTINY && d <= FTINY)
481			return(0.);
482			dist = (DOT(pos, onorm) - offs)/d;
483			if (dist < 0.)
484			return(0.);
485	greg	2.18	VSUM(oc, pos, sp->aim, dist);
486	greg	1.6	return(sp->sizdistdist/PI/(d*d));
487			}
488
489
490			double
491	greg	2.20	beamdisk( /* intersect beam with object op */
492			FVECT oc,
493			OBJREC *op,
494			SPOT *sp,
495			FVECT dir
496			)
497	greg	1.6	{
498			FVECT onorm;
499			double offs, d, dist;
500
501			offs = getplaneq(onorm, op);
502			d = -DOT(onorm, dir);
503			if (d >= -FTINY && d <= FTINY)
504			return(0.);
505			dist = (DOT(sp->aim, onorm) - offs)/d;
506	greg	2.18	VSUM(oc, sp->aim, dir, dist);
507	greg	1.6	return(sp->siz/PI/(d*d));
508			}
509
510
511			double
512	greg	2.20	intercircle( /* intersect two circles */
513			FVECT cc, /* midpoint (return value) */
514			FVECT c1, /* circle centers */
515			FVECT c2,
516			double r1s, /* radii squared */
517			double r2s
518			)
519	greg	1.4	{
520			double a2, d2, l;
521			FVECT disp;
522
523	greg	2.18	VSUB(disp, c2, c1);
524	greg	1.4	d2 = DOT(disp,disp);
525			/* circle within overlap? */
526			if (r1s < r2s) {
527			if (r2s >= r1s + d2) {
528			VCOPY(cc, c1);
529			return(r1s);
530			}
531			} else {
532			if (r1s >= r2s + d2) {
533			VCOPY(cc, c2);
534			return(r2s);
535			}
536			}
537			a2 = .25(2.(r1s+r2s) - d2 - (r2s-r1s)*(r2s-r1s)/d2);
538			/* no overlap? */
539			if (a2 <= 0.)
540			return(0.);
541			/* overlap, compute center */
542			l = sqrt((r1s - a2)/d2);
543	greg	2.18	VSUM(cc, c1, disp, l);
544	greg	1.4	return(a2);
545	greg	1.1	}