src/util/glaresrc.c

#ifndef lint
static const char       RCSid[] = "$Id: glaresrc.c,v 2.5 2004/01/02 12:51:54 schorsch Exp $";
#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 */

void    pict_stats(void);

static struct srcspan * newspan(int     l, int  r, int  v, float        *sb);
static struct srcspan * newspan(int     l, int  r, int  v, float        *sb);
static void addindirect(int     h, int  v, double       br);
static void addsrcspan(struct srcspan   *nss);
static void mergesource(struct source   *sp, struct source      *ap);
static void close_sources(int   v);
static void close_allsrcs(void);
static struct srcspan * splitspan(register struct srcspan       *sso, double    h, double       v, double       m);
static struct source * splitsource(struct source        *so);
static void donesource(register struct source   *sp);
static struct source * findbuddy(register struct source *s, register struct source      *l);
static void absorb(register struct source       *s);
static void freespans(struct source     *sp);


static struct srcspan *
newspan(                /* allocate a new source span */
        int     l,
        int     r,
        int     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);
}


extern void
analyze(void)                   /* 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((void *)spanbr);
        close_allsrcs();
}


static void
addindirect(            /* add brightness to indirect illuminances */
        int     h,
        int     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;
                }
        }
}


extern void
comp_thresh(void)                       /* 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);
        }
}


static void
addsrcspan(                     /* 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;
}


static void
mergesource(            /* merge source ap into source sp */
        struct source   *sp,
        struct source   *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((void *)ap);
}


static void
close_sources(          /* 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;
}


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

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


static struct srcspan *
splitspan(              /* divide source span at point */
        register struct srcspan *sso,
        double  h,
        double  v,
        double  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);
}


static struct source *
splitsource(                    /* 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);
}


static void
donesource(                     /* 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((void *)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);
}


static struct source *
findbuddy(                      /* find close enough source to s in l*/
        register struct source  *s,
        register struct source  *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);
}


extern void
absorb_specks(void)                     /* 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;
}


static void
absorb(                 /* 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((void *)s);
}


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

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