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"

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

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


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

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


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

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


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

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


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

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


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

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


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

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


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

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


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

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


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

        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;
                }
                free((char *)ss);
        }
        sp->first = NULL;
        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;
        }
        for ( ; s->first != NULL; s->first = s->first->next)
                free((char *)s->first);
        free((char *)s);
}
Revision:	1.14
Committed:	Mon Apr 22 08:20:58 1991 UTC (33 years ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 1.13:	+1 -2 lines
Log Message:	fixed minor bugs in last change
#	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
13			#define vcont(vd,vu) ((vu)-(vd)<=SEPS)
14			#define hcont(s1,s2) ((s1)->r-(s2)->l>=-SEPS&&(s2)->r-(s1)->l>=-SEPS)
15
16			struct source curlist = NULL; / current source list */
17			struct source donelist = NULL; / finished sources */
18
19
20			struct srcspan *
21			newspan(l, r, v, sb) /* allocate a new source span */
22			int l, r, v;
23			float *sb;
24			{
25			register struct srcspan *ss;
26			register int i;
27
28			ss = (struct srcspan *)malloc(sizeof(struct srcspan));
29			if (ss == NULL)
30			memerr("source spans");
31			ss->l = l;
32			ss->r = r;
33			ss->v = v;
34			ss->brsum = 0.0;
35			for (i = l; i < r; i++)
36			ss->brsum += sb[i+hsize];
37			return(ss);
38			}
39
40
41			analyze() /* analyze our scene */
42			{
43			int h, v;
44			int left;
45			float *spanbr;
46
47			spanbr = (float )malloc((2hsize+1)*sizeof(float));
48			if (spanbr == NULL)
49			memerr("view span brightness buffer");
50			for (v = vsize; v >= -vsize; v--) {
51	greg	1.5	close_sources(v);
52	greg	1.8	#ifndef DEBUG
53			if (verbose)
54			fprintf(stderr, "%s: analyzing... %3ld%%\r",
55			progname, 100L(vsize-v)/(2vsize));
56			#endif
57	greg	1.1	getviewspan(v, spanbr);
58			left = hsize + 1;
59	greg	1.2	for (h = -hsize; h <= hsize; h++) {
60	greg	1.1	if (spanbr[h+hsize] < 0.0) { /* off view */
61			if (left < h) {
62	greg	1.2	addsrcspan(newspan(left,h,v,spanbr));
63	greg	1.1	left = hsize + 1;
64			}
65			continue;
66			}
67			if (spanbr[h+hsize] > threshold) { /* in source */
68	greg	1.2	if (left > h)
69	greg	1.1	left = h;
70			} else { /* out of source */
71			if (left < h) {
72	greg	1.2	addsrcspan(newspan(left,h,v,spanbr));
73	greg	1.1	left = hsize + 1;
74			}
75	greg	1.13	addindirect(h, v, spanbr[h+hsize]);
76	greg	1.1	}
77	greg	1.2	}
78	greg	1.1	if (left < h)
79	greg	1.2	addsrcspan(newspan(left,h,v,spanbr));
80	greg	1.1	}
81	greg	1.7	free((char *)spanbr);
82	greg	1.1	close_allsrcs();
83			}
84
85
86	greg	1.13	addindirect(h, v, br) /* add brightness to indirect illuminances */
87			int h, v;
88	greg	1.1	double br;
89			{
90			double tanb, d;
91	greg	1.13	int hl;
92	greg	1.1	register int i;
93
94	greg	1.13	hl = hlim(v);
95			if (h <= -hl) { /* left region */
96			d = (double)(-h-hl)/sampdens;
97	greg	1.9	if (d >= 1.0-FTINY)
98	greg	1.1	return;
99			tanb = d/sqrt(1.0-d*d);
100			for (i = 0; i < nglardirs; i++) {
101			d = indirect[i].lcos - tanb*indirect[i].lsin;
102			if (d > 0.0) {
103			indirect[i].sum += d * br;
104	greg	1.10	indirect[i].n += d;
105	greg	1.1	}
106			}
107			return;
108			}
109	greg	1.13	if (h >= hl) { /* right region */
110			d = (double)(-h+hl)/sampdens;
111	greg	1.9	if (d <= -1.0+FTINY)
112	greg	1.1	return;
113			tanb = d/sqrt(1.0-d*d);
114			for (i = 0; i < nglardirs; i++) {
115			d = indirect[i].rcos - tanb*indirect[i].rsin;
116			if (d > 0.0) {
117			indirect[i].sum += d * br;
118	greg	1.10	indirect[i].n += d;
119	greg	1.1	}
120			}
121			return;
122			}
123			/* central region */
124			for (i = 0; i < nglardirs; i++) {
125	greg	1.13	d = cos(h_theta(h,v) - indirect[i].theta);
126	greg	1.1	if (d > 0.0) {
127			indirect[i].sum += d * br;
128	greg	1.10	indirect[i].n += d;
129	greg	1.1	}
130			}
131			}
132
133
134			comp_thresh() /* compute glare threshold */
135			{
136			int h, v;
137			int nsamps;
138			double brsum, br;
139
140			if (verbose)
141			fprintf(stderr, "%s: computing glare threshold...\n",
142			progname);
143			brsum = 0.0;
144			nsamps = 0;
145	greg	1.12	for (v = vsize; v >= -vsize; v -= TSAMPSTEP) {
146	greg	1.1	for (h = -hsize; h <= hsize; h += TSAMPSTEP) {
147	greg	1.13	if ((br = getviewpix(h, v)) < 0.0)
148	greg	1.1	continue;
149			brsum += br;
150			nsamps++;
151			}
152	greg	1.12	}
153	greg	1.1	if (nsamps == 0) {
154			fprintf(stderr, "%s: no viewable scene!\n", progname);
155			exit(1);
156			}
157			threshold = GLAREBR * brsum / nsamps;
158			if (threshold <= FTINY) {
159			fprintf(stderr, "%s: threshold zero!\n", progname);
160			exit(1);
161			}
162	greg	1.2	if (verbose) {
163	greg	1.8	#ifdef DEBUG
164	greg	1.2	pict_stats();
165	greg	1.8	#endif
166	greg	1.1	fprintf(stderr,
167			"%s: threshold set to %f cd/m2 from %d samples\n",
168			progname, threshold, nsamps);
169	greg	1.2	}
170	greg	1.1	}
171
172
173			addsrcspan(nss) /* add new source span to our list */
174			struct srcspan *nss;
175			{
176			struct source last, cs, *this;
177			register struct srcspan *ss;
178
179			cs = NULL;
180			for (this = curlist; this != NULL; this = this->next) {
181			for (ss = this->first; ss != NULL; ss = ss->next) {
182			if (!vcont(nss->v, ss->v))
183			break;
184			if (hcont(ss, nss)) {
185			if (cs == NULL)
186			cs = this;
187			else {
188	greg	1.7	last->next = this->next;
189	greg	1.1	mergesource(cs, this);
190	greg	1.5	this = last;
191	greg	1.1	}
192			break;
193			}
194			}
195			last = this;
196			}
197			if (cs == NULL) {
198			cs = (struct source *)malloc(sizeof(struct source));
199			if (cs == NULL)
200			memerr("source records");
201	greg	1.7	cs->dom = 0.0;
202	greg	1.1	cs->first = NULL;
203			cs->next = curlist;
204			curlist = cs;
205			}
206			nss->next = cs->first;
207			cs->first = nss;
208			}
209
210
211			mergesource(sp, ap) /* merge source ap into source sp */
212			struct source sp, ap;
213			{
214			struct srcspan head;
215			register struct srcspan alp, prev, *tp;
216
217			head.next = sp->first;
218			prev = &head;
219			alp = ap->first;
220			while (alp != NULL && prev->next != NULL) {
221			if (prev->next->v > alp->v) {
222			tp = alp->next;
223			alp->next = prev->next;
224			prev->next = alp;
225			alp = tp;
226			}
227			prev = prev->next;
228			}
229			if (prev->next == NULL)
230			prev->next = alp;
231			sp->first = head.next;
232	greg	1.7	if (ap->dom > 0.0 && sp->dom > 0.0) { /* sources are done */
233			sp->dir[0] *= sp->dom;
234			sp->dir[1] *= sp->dom;
235			sp->dir[2] *= sp->dom;
236			fvsum(sp->dir, sp->dir, ap->dir, ap->dom);
237			normalize(sp->dir);
238			sp->brt = (sp->brtsp->dom + ap->brtap->dom)
239			/ (sp->dom + ap->dom);
240			}
241			free((char *)ap);
242	greg	1.1	}
243
244
245			close_sources(v) /* close sources above v */
246			int v;
247			{
248			struct source head;
249			register struct source last, this;
250
251			head.next = curlist;
252			last = &head;
253			for (this = curlist; this != NULL; this = this->next)
254			if (!vcont(v, this->first->v)) {
255			last->next = this->next;
256			donesource(this);
257			this = last;
258			} else
259			last = this;
260			curlist = head.next;
261			}
262
263
264			close_allsrcs() /* done with everything */
265			{
266	greg	1.5	register struct source this, next;
267	greg	1.1
268	greg	1.5	this = curlist;
269			while (this != NULL) {
270			next = this->next;
271	greg	1.1	donesource(this);
272	greg	1.5	this = next;
273	greg	1.4	}
274	greg	1.1	curlist = NULL;
275			}
276
277
278			donesource(sp) /* finished with this source */
279			register struct source *sp;
280			{
281			register struct srcspan *ss;
282			int h, n;
283	greg	1.13	double hsum, vsum, d;
284	greg	1.1
285			sp->dom = 0.0;
286	greg	1.13	hsum = vsum = 0.0;
287	greg	1.1	sp->brt = 0.0;
288	greg	1.3	n = 0;
289	greg	1.1	for (ss = sp->first; ss != NULL; ss = ss->next) {
290			sp->brt += ss->brsum;
291			n += ss->r - ss->l;
292	greg	1.14	for (h = ss->l; h < ss->r; h++) {
293	greg	1.13	d = pixsize(h, ss->v);
294			hsum += d*h;
295			vsum += d*ss->v;
296			sp->dom += d;
297			}
298	greg	1.1	free((char *)ss);
299			}
300			sp->first = NULL;
301			sp->brt /= (double)n;
302	greg	1.13	compdir(sp->dir, (int)(hsum/sp->dom), (int)(vsum/sp->dom));
303	greg	1.1	sp->next = donelist;
304			donelist = sp;
305			if (verbose)
306			fprintf(stderr,
307	greg	1.12	"%s: source at [%.3f,%.3f,%.3f], dw %.5f, br %.1f (%d samps)\n",
308	greg	1.2	progname, sp->dir[0], sp->dir[1], sp->dir[2],
309	greg	1.5	sp->dom, sp->brt, n);
310	greg	1.7	}
311
312
313			struct source *
314			findbuddy(s, l) /* find close enough source to s in l*/
315			register struct source s, l;
316			{
317			struct source *bestbuddy = NULL;
318			double d, r, mindist = MAXBUDDY;
319
320			r = sqrt(s->dom/PI);
321			for ( ; l != NULL; l = l->next) {
322			d = sqrt(dist2(l->dir, s->dir)) - sqrt(l->dom/PI) - r;
323			if (d < mindist) {
324			bestbuddy = l;
325			mindist = d;
326			}
327			}
328			return(bestbuddy);
329			}
330
331
332			absorb_specks() /* eliminate too-small sources */
333			{
334			struct source head, *buddy;
335			register struct source last, this;
336
337			if (verbose)
338			fprintf(stderr, "%s: absorbing small sources...\n", progname);
339			head.next = donelist;
340			last = &head;
341	greg	1.11	for (this = head.next; this != NULL; this = this->next)
342	greg	1.7	if (TOOSMALL(this)) {
343			last->next = this->next;
344	greg	1.11	buddy = findbuddy(this, head.next);
345	greg	1.7	if (buddy != NULL)
346			mergesource(buddy, this);
347			else
348			absorb(this);
349			this = last;
350			} else
351			last = this;
352			donelist = head.next;
353			}
354
355
356			absorb(s) /* absorb a source into indirect */
357			register struct source *s;
358			{
359			FVECT dir;
360	greg	1.10	double d;
361			register int i;
362	greg	1.7
363			for (i = 0; i < nglardirs; i++) {
364			spinvector(dir, ourview.vdir, ourview.vup, indirect[i].theta);
365	greg	1.10	d = DOT(dir,s->dir)s->dom(sampdens*sampdens);
366			if (d <= 0.0)
367	greg	1.7	continue;
368	greg	1.10	indirect[i].sum += d * s->brt;
369			indirect[i].n += d;
370	greg	1.7	}
371			for ( ; s->first != NULL; s->first = s->first->next)
372			free((char *)s->first);
373			free((char *)s);
374	greg	1.1	}