src/util/glaresrc.c

#ifndef lint
static const char       RCSid[] = "$Id: glaresrc.c,v 2.4 2003/02/22 02:07:30 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 */

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