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));
                        fflush(stderr);
                }
#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;
                }
        }
        freespans(sp);
        if (sp->dom <= FTINY) {         /* must be right at edge of image */
                free((char *)sp);
                return;
        }
        sp->brt /= (double)n;
        compdir(sp->dir, (int)(hsum/sp->dom), (int)(vsum/sp->dom));
        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;

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