src/util/glaresrc.c

/* Copyright (c) 1991 Regents of the University of California */

#ifndef lint
static char SCCSid[] = "$SunId$ LBL";
#endif

/*
 * Gather samples and compute glare sources.
 */

#include "glare.h"
#include "linregr.h"

#define vcont(vd,vu)    ((vu)-(vd)<=SEPS)
#define hcont(s1,s2)    ((s1)->r-(s2)->l>=-SEPS&&(s2)->r-(s1)->l>=-SEPS)

struct source   *curlist = NULL;        /* current source list */
struct source   *donelist = NULL;       /* finished sources */


struct srcspan *
newspan(l, r, v, sb)            /* allocate a new source span */
int     l, r, v;
float   *sb;
{
        register struct srcspan *ss;
        register int    i;

        ss = (struct srcspan *)malloc(sizeof(struct srcspan));
        if (ss == NULL)
                memerr("source spans");
        ss->l = l;
        ss->r = r;
        ss->v = v;
        ss->brsum = 0.0;
        for (i = l; i < r; i++)
                ss->brsum += sb[i+hsize];
        return(ss);
}


analyze()                       /* analyze our scene */
{
        int     h, v;
        int     left;
        float   *spanbr;

        spanbr = (float *)malloc((2*hsize+1)*sizeof(float));
        if (spanbr == NULL)
                memerr("view span brightness buffer");
        for (v = vsize; v >= -vsize; v--) {
                close_sources(v);
#ifndef DEBUG
                if (verbose)
                        fprintf(stderr, "%s: analyzing... %3ld%%\r",
                                progname, 100L*(vsize-v)/(2*vsize));
#endif
                getviewspan(v, spanbr);
                left = hsize + 1;
                for (h = -hsize; h <= hsize; h++) {
                        if (spanbr[h+hsize] < 0.0) {    /* off view */
                                if (left < h) {
                                        addsrcspan(newspan(left,h,v,spanbr));
                                        left = hsize + 1;
                                }
                                continue;
                        }
                        if (spanbr[h+hsize] > threshold) {      /* in source */
                                if (left > h)
                                        left = h;
                        } else {                        /* out of source */
                                if (left < h) {
                                        addsrcspan(newspan(left,h,v,spanbr));
                                        left = hsize + 1;
                                }
                                addindirect(h, v, spanbr[h+hsize]);
                        }
                }
                if (left < h)
                        addsrcspan(newspan(left,h,v,spanbr));
        }
        free((char *)spanbr);
        close_allsrcs();
}


addindirect(h, v, br)           /* add brightness to indirect illuminances */
int     h, v;
double  br;
{
        double  tanb, d;
        int     hl;
        register int    i;

        hl = hlim(v);
        if (h <= -hl) {                 /* left region */
                d = (double)(-h-hl)/sampdens;
                if (d >= 1.0-FTINY)
                        return;
                tanb = d/sqrt(1.0-d*d);
                for (i = 0; i < nglardirs; i++) {
                        d = indirect[i].lcos - tanb*indirect[i].lsin;
                        if (d > 0.0) {
                                indirect[i].sum += d * br;
                                indirect[i].n += d;
                        }
                }
                return;
        }
        if (h >= hl) {                  /* right region */
                d = (double)(-h+hl)/sampdens;
                if (d <= -1.0+FTINY)
                        return;
                tanb = d/sqrt(1.0-d*d);
                for (i = 0; i < nglardirs; i++) {
                        d = indirect[i].rcos - tanb*indirect[i].rsin;
                        if (d > 0.0) {
                                indirect[i].sum += d * br;
                                indirect[i].n += d;
                        }
                }
                return;
        }
                                        /* central region */
        for (i = 0; i < nglardirs; i++) {
                d = cos(h_theta(h,v) - indirect[i].theta);
                if (d > 0.0) {
                        indirect[i].sum += d * br;
                        indirect[i].n += d;
                }
        }
}


comp_thresh()                   /* compute glare threshold */
{
        int     h, v;
        int     nsamps;
        double  brsum, br;

        if (verbose)
                fprintf(stderr, "%s: computing glare threshold...\n",
                                progname);
        brsum = 0.0;
        nsamps = 0;
        for (v = vsize; v >= -vsize; v -= TSAMPSTEP) {
                for (h = -hsize; h <= hsize; h += TSAMPSTEP) {
                        if ((br = getviewpix(h, v)) < 0.0)
                                continue;
                        brsum += br;
                        nsamps++;
                }
        }
        if (nsamps == 0) {
                fprintf(stderr, "%s: no viewable scene!\n", progname);
                exit(1);
        }
        threshold = GLAREBR * brsum / nsamps;
        if (threshold <= FTINY) {
                fprintf(stderr, "%s: threshold zero!\n", progname);
                exit(1);
        }
        if (verbose) {
#ifdef DEBUG
                pict_stats();
#endif
                fprintf(stderr,
                        "%s: threshold set to %f cd/m2 from %d samples\n",
                                progname, threshold, nsamps);
        }
}


addsrcspan(nss)                 /* add new source span to our list */
struct srcspan  *nss;
{
        struct source   *last, *cs, *this;
        register struct srcspan *ss;

        cs = NULL;
        for (this = curlist; this != NULL; this = this->next) {
                for (ss = this->first; ss != NULL; ss = ss->next) {
                        if (!vcont(nss->v, ss->v))
                                break;
                        if (hcont(ss, nss)) {
                                if (cs == NULL)
                                        cs = this;
                                else {
                                        last->next = this->next;
                                        mergesource(cs, this);
                                        this = last;
                                }
                                break;
                        }
                }
                last = this;
        }
        if (cs == NULL) {
                cs = (struct source *)malloc(sizeof(struct source));
                if (cs == NULL)
                        memerr("source records");
                cs->dom = 0.0;
                cs->first = NULL;
                cs->next = curlist;
                curlist = cs;
        }
        nss->next = cs->first;
        cs->first = nss;
}


mergesource(sp, ap)             /* merge source ap into source sp */
struct source   *sp, *ap;
{
        struct srcspan  head;
        register struct srcspan *alp, *prev, *tp;

        head.next = sp->first;
        prev = &head;
        alp = ap->first;
        while (alp != NULL && prev->next != NULL) {
                if (prev->next->v > alp->v) {
                        tp = alp->next;
                        alp->next = prev->next;
                        prev->next = alp;
                        alp = tp;
                }
                prev = prev->next;
        }
        if (prev->next == NULL)
                prev->next = alp;
        sp->first = head.next;
        if (ap->dom > 0.0 && sp->dom > 0.0) {           /* sources are done */
                sp->dir[0] *= sp->dom;
                sp->dir[1] *= sp->dom;
                sp->dir[2] *= sp->dom;
                fvsum(sp->dir, sp->dir, ap->dir, ap->dom);
                normalize(sp->dir);
                sp->brt = (sp->brt*sp->dom + ap->brt*ap->dom)
                                / (sp->dom + ap->dom);
        }
        free((char *)ap);
}


close_sources(v)                /* close sources above v */
int     v;
{
        struct source   head;
        register struct source  *last, *this;

        head.next = curlist;
        last = &head;
        for (this = curlist; this != NULL; this = this->next)
                if (!vcont(v, this->first->v)) {
                        last->next = this->next;
                        donesource(this);
                        this = last;
                } else
                        last = this;
        curlist = head.next;
}


close_allsrcs()                 /* done with everything */
{
        register struct source  *this, *next;

        this = curlist;
        while (this != NULL) {
                next = this->next;
                donesource(this);
                this = next;
        }
        curlist = NULL;
}


struct srcspan *
splitspan(sso, h, v, m)         /* divide source span at point */
register struct srcspan *sso;
double  h, v, m;
{
        register struct srcspan *ssn;
        double  d;
        int     hs;

        d = h - m*(sso->v - v);
        hs = d < 0. ? d-.5 : d+.5;
        if (sso->l >= hs)
                return(NULL);
        if (sso->r <= hs)
                return(sso);
                                /* need to split it */
        ssn = (struct srcspan *)malloc(sizeof(struct srcspan));
        if (ssn == NULL)
                memerr("source spans in splitspan");
        ssn->brsum = (double)(hs - sso->l)/(sso->r - sso->l) * sso->brsum;
        sso->brsum -= ssn->brsum;
        ssn->v = sso->v;
        ssn->l = sso->l;
        ssn->r = sso->l = hs;
        return(ssn);
}


struct source *
splitsource(so)                 /* divide source in two if it's big and long */
struct source   *so;
{
        LRSUM   lr;
        LRLIN   fit;
        register struct srcspan *ss, *ssn;
        struct srcspan  *ssl, *ssnl, head;
        int     h;
        double  mh, mv;
        struct source   *sn;

        lrclear(&lr);
        for (ss = so->first; ss != NULL; ss = ss->next)
                for (h = ss->l; h < ss->r; h++)
                        lrpoint(h, ss->v, &lr);
        if ((double)lr.n/(sampdens*sampdens) < SABIG)
                return(NULL);                   /* too small */
        if (lrfit(&fit, &lr) < 0)
                return(NULL);                   /* can't fit a line */
        if (fit.correlation < LCORR && fit.correlation > -LCORR)
                return(NULL);
        if (verbose)
                fprintf(stderr, "%s: splitting large source\n", progname);
        mh = lrxavg(&lr);
        mv = lryavg(&lr);
        sn = (struct source *)malloc(sizeof(struct source));
        if (sn == NULL)
                memerr("source records in splitsource");
        sn->dom = 0.0;
        sn->first = NULL;
        ssnl = NULL;
        head.next = so->first;
        ssl = &head;
        for (ss = so->first; ss != NULL; ssl = ss, ss = ss->next)
                if ((ssn = splitspan(ss, mh, mv, fit.slope)) != NULL) {
                        if (ssn == ss) {        /* remove from old */
                                ssl->next = ss->next;
                                ss = ssl;
                        }
                        if (ssnl == NULL)       /* add to new */
                                sn->first = ssn;
                        else
                                ssnl->next = ssn;
                        ssn->next = NULL;
                        ssnl = ssn;
                }
        so->first = head.next;
        return(sn);
}


donesource(sp)                  /* finished with this source */
register struct source  *sp;
{
        struct source   *newsrc;
        register struct srcspan *ss;
        int     h, n;
        double  hsum, vsum, d;

        while ((newsrc = splitsource(sp)) != NULL)      /* split it? */
                donesource(newsrc);
        sp->dom = 0.0;
        hsum = vsum = 0.0;
        sp->brt = 0.0;
        n = 0;
        for (ss = sp->first; ss != NULL; ss = ss->next) {
                sp->brt += ss->brsum;
                n += ss->r - ss->l;
                for (h = ss->l; h < ss->r; h++) {
                        d = pixsize(h, ss->v);
                        hsum += d*h;
                        vsum += d*ss->v;
                        sp->dom += d;
                }
        }
        sp->brt /= (double)n;
        compdir(sp->dir, (int)(hsum/sp->dom), (int)(vsum/sp->dom));
        freespans(sp);
        sp->next = donelist;
        donelist = sp;
        if (verbose)
                fprintf(stderr,
        "%s: source at [%.3f,%.3f,%.3f], dw %.5f, br %.1f (%d samps)\n",
                        progname, sp->dir[0], sp->dir[1], sp->dir[2], 
                        sp->dom, sp->brt, n);
}


struct source *
findbuddy(s, l)                 /* find close enough source to s in l*/
register struct source  *s, *l;
{
        struct source   *bestbuddy = NULL;
        double  d, r, mindist = MAXBUDDY;

        r = sqrt(s->dom/PI);
        for ( ; l != NULL; l = l->next) {
                d = sqrt(dist2(l->dir, s->dir)) - sqrt(l->dom/PI) - r;
                if (d < mindist) {
                        bestbuddy = l;
                        mindist = d;
                }
        }
        return(bestbuddy);
}


absorb_specks()                 /* eliminate too-small sources */
{
        struct source   head, *buddy;
        register struct source  *last, *this;

        if (verbose)
                fprintf(stderr, "%s: absorbing small sources...\n", progname);
        head.next = donelist;
        last = &head;
        for (this = head.next; this != NULL; this = this->next)
                if (TOOSMALL(this)) {
                        last->next = this->next;
                        buddy = findbuddy(this, head.next);
                        if (buddy != NULL)
                                mergesource(buddy, this);
                        else
                                absorb(this);
                        this = last;
                } else
                        last = this;
        donelist = head.next;
}


absorb(s)                       /* absorb a source into indirect */
register struct source  *s;
{
        FVECT   dir;
        double  d;
        register int    i;

        for (i = 0; i < nglardirs; i++) {
                spinvector(dir, ourview.vdir, ourview.vup, indirect[i].theta);
                d = DOT(dir,s->dir)*s->dom*(sampdens*sampdens);
                if (d <= 0.0)
                        continue;
                indirect[i].sum += d * s->brt;
                indirect[i].n += d;
        }
        freespans(s);
        free((char *)s);
}


freespans(sp)                   /* free spans associated with source */
struct source   *sp;
{
        register struct srcspan *ss;

        for (ss = sp->first; ss != NULL; ss = ss->next)
                free((char *)ss);
        sp->first = NULL;
}
Revision:	1.15
Committed:	Tue Apr 30 16:39:39 1991 UTC (34 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.14:	+97 -4 lines
Log Message:	added routines for splitting long sources
#	User	Rev	Content
1	greg	1.1	/* Copyright (c) 1991 Regents of the University of California */
2
3			#ifndef lint
4			static char SCCSid[] = "$SunId$ LBL";
5			#endif
6
7			/*
8			* Gather samples and compute glare sources.
9			*/
10
11			#include "glare.h"
12	greg	1.15	#include "linregr.h"
13	greg	1.1
14			#define vcont(vd,vu) ((vu)-(vd)<=SEPS)
15			#define hcont(s1,s2) ((s1)->r-(s2)->l>=-SEPS&&(s2)->r-(s1)->l>=-SEPS)
16
17			struct source curlist = NULL; / current source list */
18			struct source donelist = NULL; / finished sources */
19
20
21			struct srcspan *
22			newspan(l, r, v, sb) /* allocate a new source span */
23			int l, r, v;
24			float *sb;
25			{
26			register struct srcspan *ss;
27			register int i;
28
29			ss = (struct srcspan *)malloc(sizeof(struct srcspan));
30			if (ss == NULL)
31			memerr("source spans");
32			ss->l = l;
33			ss->r = r;
34			ss->v = v;
35			ss->brsum = 0.0;
36			for (i = l; i < r; i++)
37			ss->brsum += sb[i+hsize];
38			return(ss);
39			}
40
41
42			analyze() /* analyze our scene */
43			{
44			int h, v;
45			int left;
46			float *spanbr;
47
48			spanbr = (float )malloc((2hsize+1)*sizeof(float));
49			if (spanbr == NULL)
50			memerr("view span brightness buffer");
51			for (v = vsize; v >= -vsize; v--) {
52	greg	1.5	close_sources(v);
53	greg	1.8	#ifndef DEBUG
54			if (verbose)
55			fprintf(stderr, "%s: analyzing... %3ld%%\r",
56			progname, 100L(vsize-v)/(2vsize));
57			#endif
58	greg	1.1	getviewspan(v, spanbr);
59			left = hsize + 1;
60	greg	1.2	for (h = -hsize; h <= hsize; h++) {
61	greg	1.1	if (spanbr[h+hsize] < 0.0) { /* off view */
62			if (left < h) {
63	greg	1.2	addsrcspan(newspan(left,h,v,spanbr));
64	greg	1.1	left = hsize + 1;
65			}
66			continue;
67			}
68			if (spanbr[h+hsize] > threshold) { /* in source */
69	greg	1.2	if (left > h)
70	greg	1.1	left = h;
71			} else { /* out of source */
72			if (left < h) {
73	greg	1.2	addsrcspan(newspan(left,h,v,spanbr));
74	greg	1.1	left = hsize + 1;
75			}
76	greg	1.13	addindirect(h, v, spanbr[h+hsize]);
77	greg	1.1	}
78	greg	1.2	}
79	greg	1.1	if (left < h)
80	greg	1.2	addsrcspan(newspan(left,h,v,spanbr));
81	greg	1.1	}
82	greg	1.7	free((char *)spanbr);
83	greg	1.1	close_allsrcs();
84			}
85
86
87	greg	1.13	addindirect(h, v, br) /* add brightness to indirect illuminances */
88			int h, v;
89	greg	1.1	double br;
90			{
91			double tanb, d;
92	greg	1.13	int hl;
93	greg	1.1	register int i;
94
95	greg	1.13	hl = hlim(v);
96			if (h <= -hl) { /* left region */
97			d = (double)(-h-hl)/sampdens;
98	greg	1.9	if (d >= 1.0-FTINY)
99	greg	1.1	return;
100			tanb = d/sqrt(1.0-d*d);
101			for (i = 0; i < nglardirs; i++) {
102			d = indirect[i].lcos - tanb*indirect[i].lsin;
103			if (d > 0.0) {
104			indirect[i].sum += d * br;
105	greg	1.10	indirect[i].n += d;
106	greg	1.1	}
107			}
108			return;
109			}
110	greg	1.13	if (h >= hl) { /* right region */
111			d = (double)(-h+hl)/sampdens;
112	greg	1.9	if (d <= -1.0+FTINY)
113	greg	1.1	return;
114			tanb = d/sqrt(1.0-d*d);
115			for (i = 0; i < nglardirs; i++) {
116			d = indirect[i].rcos - tanb*indirect[i].rsin;
117			if (d > 0.0) {
118			indirect[i].sum += d * br;
119	greg	1.10	indirect[i].n += d;
120	greg	1.1	}
121			}
122			return;
123			}
124			/* central region */
125			for (i = 0; i < nglardirs; i++) {
126	greg	1.13	d = cos(h_theta(h,v) - indirect[i].theta);
127	greg	1.1	if (d > 0.0) {
128			indirect[i].sum += d * br;
129	greg	1.10	indirect[i].n += d;
130	greg	1.1	}
131			}
132			}
133
134
135			comp_thresh() /* compute glare threshold */
136			{
137			int h, v;
138			int nsamps;
139			double brsum, br;
140
141			if (verbose)
142			fprintf(stderr, "%s: computing glare threshold...\n",
143			progname);
144			brsum = 0.0;
145			nsamps = 0;
146	greg	1.12	for (v = vsize; v >= -vsize; v -= TSAMPSTEP) {
147	greg	1.1	for (h = -hsize; h <= hsize; h += TSAMPSTEP) {
148	greg	1.13	if ((br = getviewpix(h, v)) < 0.0)
149	greg	1.1	continue;
150			brsum += br;
151			nsamps++;
152			}
153	greg	1.12	}
154	greg	1.1	if (nsamps == 0) {
155			fprintf(stderr, "%s: no viewable scene!\n", progname);
156			exit(1);
157			}
158			threshold = GLAREBR * brsum / nsamps;
159			if (threshold <= FTINY) {
160			fprintf(stderr, "%s: threshold zero!\n", progname);
161			exit(1);
162			}
163	greg	1.2	if (verbose) {
164	greg	1.8	#ifdef DEBUG
165	greg	1.2	pict_stats();
166	greg	1.8	#endif
167	greg	1.1	fprintf(stderr,
168			"%s: threshold set to %f cd/m2 from %d samples\n",
169			progname, threshold, nsamps);
170	greg	1.2	}
171	greg	1.1	}
172
173
174			addsrcspan(nss) /* add new source span to our list */
175			struct srcspan *nss;
176			{
177			struct source last, cs, *this;
178			register struct srcspan *ss;
179
180			cs = NULL;
181			for (this = curlist; this != NULL; this = this->next) {
182			for (ss = this->first; ss != NULL; ss = ss->next) {
183			if (!vcont(nss->v, ss->v))
184			break;
185			if (hcont(ss, nss)) {
186			if (cs == NULL)
187			cs = this;
188			else {
189	greg	1.7	last->next = this->next;
190	greg	1.1	mergesource(cs, this);
191	greg	1.5	this = last;
192	greg	1.1	}
193			break;
194			}
195			}
196			last = this;
197			}
198			if (cs == NULL) {
199			cs = (struct source *)malloc(sizeof(struct source));
200			if (cs == NULL)
201			memerr("source records");
202	greg	1.7	cs->dom = 0.0;
203	greg	1.1	cs->first = NULL;
204			cs->next = curlist;
205			curlist = cs;
206			}
207			nss->next = cs->first;
208			cs->first = nss;
209			}
210
211
212			mergesource(sp, ap) /* merge source ap into source sp */
213			struct source sp, ap;
214			{
215			struct srcspan head;
216			register struct srcspan alp, prev, *tp;
217
218			head.next = sp->first;
219			prev = &head;
220			alp = ap->first;
221			while (alp != NULL && prev->next != NULL) {
222			if (prev->next->v > alp->v) {
223			tp = alp->next;
224			alp->next = prev->next;
225			prev->next = alp;
226			alp = tp;
227			}
228			prev = prev->next;
229			}
230			if (prev->next == NULL)
231			prev->next = alp;
232			sp->first = head.next;
233	greg	1.7	if (ap->dom > 0.0 && sp->dom > 0.0) { /* sources are done */
234			sp->dir[0] *= sp->dom;
235			sp->dir[1] *= sp->dom;
236			sp->dir[2] *= sp->dom;
237			fvsum(sp->dir, sp->dir, ap->dir, ap->dom);
238			normalize(sp->dir);
239			sp->brt = (sp->brtsp->dom + ap->brtap->dom)
240			/ (sp->dom + ap->dom);
241			}
242			free((char *)ap);
243	greg	1.1	}
244
245
246			close_sources(v) /* close sources above v */
247			int v;
248			{
249			struct source head;
250			register struct source last, this;
251
252			head.next = curlist;
253			last = &head;
254			for (this = curlist; this != NULL; this = this->next)
255			if (!vcont(v, this->first->v)) {
256			last->next = this->next;
257			donesource(this);
258			this = last;
259			} else
260			last = this;
261			curlist = head.next;
262			}
263
264
265			close_allsrcs() /* done with everything */
266			{
267	greg	1.5	register struct source this, next;
268	greg	1.1
269	greg	1.5	this = curlist;
270			while (this != NULL) {
271			next = this->next;
272	greg	1.1	donesource(this);
273	greg	1.5	this = next;
274	greg	1.4	}
275	greg	1.1	curlist = NULL;
276			}
277
278
279	greg	1.15	struct srcspan *
280			splitspan(sso, h, v, m) /* divide source span at point */
281			register struct srcspan *sso;
282			double h, v, m;
283			{
284			register struct srcspan *ssn;
285			double d;
286			int hs;
287
288			d = h - m*(sso->v - v);
289			hs = d < 0. ? d-.5 : d+.5;
290			if (sso->l >= hs)
291			return(NULL);
292			if (sso->r <= hs)
293			return(sso);
294			/* need to split it */
295			ssn = (struct srcspan *)malloc(sizeof(struct srcspan));
296			if (ssn == NULL)
297			memerr("source spans in splitspan");
298			ssn->brsum = (double)(hs - sso->l)/(sso->r - sso->l) * sso->brsum;
299			sso->brsum -= ssn->brsum;
300			ssn->v = sso->v;
301			ssn->l = sso->l;
302			ssn->r = sso->l = hs;
303			return(ssn);
304			}
305
306
307			struct source *
308			splitsource(so) /* divide source in two if it's big and long */
309			struct source *so;
310			{
311			LRSUM lr;
312			LRLIN fit;
313			register struct srcspan ss, ssn;
314			struct srcspan ssl, ssnl, head;
315			int h;
316			double mh, mv;
317			struct source *sn;
318
319			lrclear(&lr);
320			for (ss = so->first; ss != NULL; ss = ss->next)
321			for (h = ss->l; h < ss->r; h++)
322			lrpoint(h, ss->v, &lr);
323			if ((double)lr.n/(sampdens*sampdens) < SABIG)
324			return(NULL); /* too small */
325			if (lrfit(&fit, &lr) < 0)
326			return(NULL); /* can't fit a line */
327			if (fit.correlation < LCORR && fit.correlation > -LCORR)
328			return(NULL);
329			if (verbose)
330			fprintf(stderr, "%s: splitting large source\n", progname);
331			mh = lrxavg(&lr);
332			mv = lryavg(&lr);
333			sn = (struct source *)malloc(sizeof(struct source));
334			if (sn == NULL)
335			memerr("source records in splitsource");
336			sn->dom = 0.0;
337			sn->first = NULL;
338			ssnl = NULL;
339			head.next = so->first;
340			ssl = &head;
341			for (ss = so->first; ss != NULL; ssl = ss, ss = ss->next)
342			if ((ssn = splitspan(ss, mh, mv, fit.slope)) != NULL) {
343			if (ssn == ss) { /* remove from old */
344			ssl->next = ss->next;
345			ss = ssl;
346			}
347			if (ssnl == NULL) /* add to new */
348			sn->first = ssn;
349			else
350			ssnl->next = ssn;
351			ssn->next = NULL;
352			ssnl = ssn;
353			}
354			so->first = head.next;
355			return(sn);
356			}
357
358
359	greg	1.1	donesource(sp) /* finished with this source */
360			register struct source *sp;
361			{
362	greg	1.15	struct source *newsrc;
363	greg	1.1	register struct srcspan *ss;
364			int h, n;
365	greg	1.13	double hsum, vsum, d;
366	greg	1.1
367	greg	1.15	while ((newsrc = splitsource(sp)) != NULL) /* split it? */
368			donesource(newsrc);
369	greg	1.1	sp->dom = 0.0;
370	greg	1.13	hsum = vsum = 0.0;
371	greg	1.1	sp->brt = 0.0;
372	greg	1.3	n = 0;
373	greg	1.1	for (ss = sp->first; ss != NULL; ss = ss->next) {
374			sp->brt += ss->brsum;
375			n += ss->r - ss->l;
376	greg	1.14	for (h = ss->l; h < ss->r; h++) {
377	greg	1.13	d = pixsize(h, ss->v);
378			hsum += d*h;
379			vsum += d*ss->v;
380			sp->dom += d;
381			}
382	greg	1.1	}
383			sp->brt /= (double)n;
384	greg	1.13	compdir(sp->dir, (int)(hsum/sp->dom), (int)(vsum/sp->dom));
385	greg	1.15	freespans(sp);
386	greg	1.1	sp->next = donelist;
387			donelist = sp;
388			if (verbose)
389			fprintf(stderr,
390	greg	1.12	"%s: source at [%.3f,%.3f,%.3f], dw %.5f, br %.1f (%d samps)\n",
391	greg	1.2	progname, sp->dir[0], sp->dir[1], sp->dir[2],
392	greg	1.5	sp->dom, sp->brt, n);
393	greg	1.7	}
394
395
396			struct source *
397			findbuddy(s, l) /* find close enough source to s in l*/
398			register struct source s, l;
399			{
400			struct source *bestbuddy = NULL;
401			double d, r, mindist = MAXBUDDY;
402
403			r = sqrt(s->dom/PI);
404			for ( ; l != NULL; l = l->next) {
405			d = sqrt(dist2(l->dir, s->dir)) - sqrt(l->dom/PI) - r;
406			if (d < mindist) {
407			bestbuddy = l;
408			mindist = d;
409			}
410			}
411			return(bestbuddy);
412			}
413
414
415			absorb_specks() /* eliminate too-small sources */
416			{
417			struct source head, *buddy;
418			register struct source last, this;
419
420			if (verbose)
421			fprintf(stderr, "%s: absorbing small sources...\n", progname);
422			head.next = donelist;
423			last = &head;
424	greg	1.11	for (this = head.next; this != NULL; this = this->next)
425	greg	1.7	if (TOOSMALL(this)) {
426			last->next = this->next;
427	greg	1.11	buddy = findbuddy(this, head.next);
428	greg	1.7	if (buddy != NULL)
429			mergesource(buddy, this);
430			else
431			absorb(this);
432			this = last;
433			} else
434			last = this;
435			donelist = head.next;
436			}
437
438
439			absorb(s) /* absorb a source into indirect */
440			register struct source *s;
441			{
442			FVECT dir;
443	greg	1.10	double d;
444			register int i;
445	greg	1.7
446			for (i = 0; i < nglardirs; i++) {
447			spinvector(dir, ourview.vdir, ourview.vup, indirect[i].theta);
448	greg	1.10	d = DOT(dir,s->dir)s->dom(sampdens*sampdens);
449			if (d <= 0.0)
450	greg	1.7	continue;
451	greg	1.10	indirect[i].sum += d * s->brt;
452			indirect[i].n += d;
453	greg	1.7	}
454	greg	1.15	freespans(s);
455	greg	1.7	free((char *)s);
456	greg	1.15	}
457
458
459			freespans(sp) /* free spans associated with source */
460			struct source *sp;
461			{
462			register struct srcspan *ss;
463
464			for (ss = sp->first; ss != NULL; ss = ss->next)
465			free((char *)ss);
466			sp->first = NULL;
467	greg	1.1	}