src/util/glaresrc.c

#ifndef lint
static const char       RCSid[] = "$Id: glaresrc.c,v 2.6 2006/05/29 16:47:54 greg 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 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(struct srcspan        *sso, double    h, double       v, double       m);
static struct source * splitsource(struct source        *so);
static void donesource(struct source    *sp);
static struct source * findbuddy(struct source  *s, struct source       *l);
static void absorb(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
)
{
        struct srcspan  *ss;
        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);
}


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;
        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;
                }
        }
}


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;
        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;
        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;
        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 */
{
        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 */
        struct srcspan  *sso,
        double  h,
        double  v,
        double  m
)
{
        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;
        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 */
        struct source   *sp
)
{
        struct source   *newsrc;
        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*/
        struct source   *s,
        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);
}


void
absorb_specks(void)                     /* eliminate too-small sources */
{
        struct source   head, *buddy;
        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 */
        struct source   *s
)
{
        FVECT   dir;
        double  d;
        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
)
{
        struct srcspan  *ss;

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