ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.13
Committed: Thu Aug 2 18:33:50 2018 UTC (5 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2
Changes since 2.12: +2 -2 lines
Log Message:
Created MAXFMTLEN to guard against buffer overrun attacks in header input

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.13 static const char RCSid[] = "$Id: glarendx.c,v 2.12 2018/04/02 16:45:41 greg Exp $";
3 greg 1.1 #endif
4     /*
5 greg 1.9 * Compute Glare Index given by program name or -t option:
6 greg 1.1 *
7 greg 2.4 * dgi - Daylight Glare Index
8     * brs_gi - Building Research Station Glare Index (Petherbridge
9     * & Hopkinson)
10     * ugr - Unified Glare Rating System (Fischer)
11 greg 1.2 * guth_dgr - Guth discomfort glare rating
12     * guth_vcp - Guth visual comfort probability
13 greg 1.9 * cie_cgi - CIE Glare Index (1983, due to Einhorn)
14     * vert_dir - Direct vertical illuminance
15     * vert_ind - Indirect vertical illuminance (from input)
16     * vert_ill - Total vertical illuminance
17 greg 1.1 *
18     * 12 April 1991 Greg Ward EPFL
19 greg 2.4 * 19 April 1993 R. Compagnon EPFL (added dgi, brs_gi, ugr)
20 greg 1.1 */
21 schorsch 2.10
22 schorsch 2.8 #include <string.h>
23    
24 greg 1.1 #include "standard.h"
25     #include "view.h"
26 schorsch 2.10
27    
28     struct glare_src {
29     FVECT dir; /* source direction */
30     double dom; /* solid angle */
31     double lum; /* average luminance */
32     struct glare_src *next;
33     } *all_srcs = NULL;
34    
35     struct glare_dir {
36     double ang; /* angle (in radians) */
37     double indirect; /* indirect illuminance */
38     struct glare_dir *next;
39     } *all_dirs = NULL;
40    
41     typedef double gdfun(struct glare_dir *gd);
42    
43     static gdfun dgi;
44     static gdfun brs_gi;
45     static gdfun ugr;
46     static gdfun guth_dgr;
47     static gdfun guth_vcp;
48     static gdfun cie_cgi;
49     static gdfun direct;
50     static gdfun indirect;
51     static gdfun total;
52    
53     static gethfunc headline;
54     static void init(void);
55     static void read_input(void);
56     static void print_values(gdfun *func);
57     static double posindex(FVECT sd, FVECT vd, FVECT vu);
58    
59    
60 greg 1.1 struct named_func {
61     char *name;
62 schorsch 2.10 gdfun *func;
63 greg 1.9 char *descrip;
64 greg 1.1 } all_funcs[] = {
65 greg 2.4 {"dgi", dgi, "Daylight Glare Index"},
66     {"brs_gi", brs_gi, "BRS Glare Index"},
67     {"ugr", ugr, "Unified Glare Rating"},
68 greg 1.9 {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
69     {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
70     {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
71     {"vert_dir", direct, "Direct Vertical Illuminance"},
72     {"vert_ill", total, "Total Vertical Illuminance"},
73     {"vert_ind", indirect, "Indirect Vertical Illuminance"},
74 greg 1.1 {NULL}
75     };
76 schorsch 2.10
77 greg 1.1 #define newp(type) (type *)malloc(sizeof(type))
78 schorsch 2.10
79 greg 1.1 char *progname;
80     int print_header = 1;
81 schorsch 2.10
82 greg 1.1 VIEW midview = STDVIEW;
83 schorsch 2.10
84 greg 1.4 int wrongformat = 0;
85 schorsch 2.10
86    
87     int
88     main(
89     int argc,
90     char *argv[]
91     )
92 greg 1.1 {
93     struct named_func *funp;
94     char *progtail;
95     int i;
96     /* get program name */
97     progname = argv[0];
98 schorsch 2.8 progtail = strrchr(progname, '/'); /* final component */
99 greg 1.1 if (progtail == NULL)
100     progtail = progname;
101     else
102     progtail++;
103     /* get options */
104     for (i = 1; i < argc && argv[i][0] == '-'; i++)
105     switch (argv[i][1]) {
106     case 't':
107     progtail = argv[++i];
108     break;
109     case 'h':
110     print_header = 0;
111     break;
112     default:
113     goto userr;
114     }
115     if (i < argc-1)
116     goto userr;
117     if (i == argc-1) /* open file */
118     if (freopen(argv[i], "r", stdin) == NULL) {
119     perror(argv[i]);
120     exit(1);
121     }
122 greg 1.9 /* find and run calculation */
123 greg 1.1 for (funp = all_funcs; funp->name != NULL; funp++)
124     if (!strcmp(funp->name, progtail)) {
125 greg 1.9 init();
126     read_input();
127     if (print_header) {
128     printargs(i, argv, stdout);
129     putchar('\n');
130     }
131 greg 1.1 print_values(funp->func);
132     exit(0); /* we're done */
133     }
134     /* invalid function */
135     userr:
136 greg 1.9 fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
137     fprintf(stderr, "\twhere 'type' is one of the following:\n");
138     for (funp = all_funcs; funp->name != NULL; funp++)
139     fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
140 greg 1.1 exit(1);
141     }
142 schorsch 2.10
143    
144     static int
145     headline( /* get line from header */
146     char *s,
147     void *p
148     )
149 greg 1.1 {
150 greg 2.13 char fmt[MAXFMTLEN];
151 schorsch 2.10
152 greg 1.1 if (print_header) /* copy to output */
153     fputs(s, stdout);
154 greg 2.2 if (isview(s))
155     sscanview(&midview, s);
156 greg 1.4 else if (isformat(s)) {
157     formatval(fmt, s);
158 greg 1.6 wrongformat = strcmp(fmt, "ascii");
159 greg 1.4 }
160 gwlarson 2.6 return(0);
161 greg 1.1 }
162 schorsch 2.10
163    
164     static void
165     init(void) /* initialize calculation */
166 greg 1.9 {
167     /* read header */
168     getheader(stdin, headline, NULL);
169     if (wrongformat) {
170     fprintf(stderr, "%s: bad input format\n", progname);
171     exit(1);
172     }
173     /* set view */
174     if (setview(&midview) != NULL) {
175     fprintf(stderr, "%s: bad view information in input\n",
176     progname);
177     exit(1);
178     }
179     }
180 schorsch 2.10
181    
182     static void
183     read_input(void) /* read glare sources from stdin */
184 greg 1.1 {
185     #define S_SEARCH 0
186     #define S_SOURCE 1
187     #define S_DIREC 2
188     int state = S_SEARCH;
189     char buf[128];
190 greg 2.11 struct glare_src *gs;
191     struct glare_dir *gd;
192 schorsch 2.10
193 greg 1.1 while (fgets(buf, sizeof(buf), stdin) != NULL)
194     switch (state) {
195     case S_SEARCH:
196     if (!strcmp(buf, "BEGIN glare source\n"))
197     state = S_SOURCE;
198     else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
199     state = S_DIREC;
200     break;
201     case S_SOURCE:
202     if (!strncmp(buf, "END", 3)) {
203     state = S_SEARCH;
204     break;
205     }
206     if ((gs = newp(struct glare_src)) == NULL)
207     goto memerr;
208 greg 2.11 if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
209     &gs->dir[2]) != 3 ||
210 greg 2.12 sscanf(sskip2(buf, 3), "%lf %lf",
211 greg 2.11 &gs->dom, &gs->lum) != 2)
212 greg 1.1 goto readerr;
213 greg 1.2 normalize(gs->dir);
214 greg 1.1 gs->next = all_srcs;
215     all_srcs = gs;
216     break;
217     case S_DIREC:
218     if (!strncmp(buf, "END", 3)) {
219     state = S_SEARCH;
220     break;
221     }
222     if ((gd = newp(struct glare_dir)) == NULL)
223     goto memerr;
224     if (sscanf(buf, "%lf %lf",
225     &gd->ang, &gd->indirect) != 2)
226     goto readerr;
227     gd->ang *= PI/180.0; /* convert to radians */
228     gd->next = all_dirs;
229     all_dirs = gd;
230     break;
231     }
232     return;
233     memerr:
234     perror(progname);
235     exit(1);
236     readerr:
237     fprintf(stderr, "%s: read error on input\n", progname);
238     exit(1);
239     #undef S_SEARCH
240     #undef S_SOURCE
241     #undef S_DIREC
242     }
243 schorsch 2.10
244    
245     static void
246     print_values( /* print out calculations */
247     gdfun *func
248     )
249 greg 1.1 {
250 greg 2.11 struct glare_dir *gd;
251 schorsch 2.10
252 greg 1.1 for (gd = all_dirs; gd != NULL; gd = gd->next)
253 schorsch 2.10 printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
254 greg 1.1 }
255 schorsch 2.10
256    
257     static double
258     direct( /* compute direct vertical illuminance */
259     struct glare_dir *gd
260     )
261 greg 1.1 {
262     FVECT mydir;
263     double d, dval;
264 greg 2.11 struct glare_src *gs;
265 schorsch 2.10
266 greg 1.1 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
267     dval = 0.0;
268     for (gs = all_srcs; gs != NULL; gs = gs->next) {
269     d = DOT(mydir,gs->dir);
270     if (d > FTINY)
271     dval += d * gs->dom * gs->lum;
272     }
273     return(dval);
274 greg 1.9 }
275 schorsch 2.10
276    
277     static double
278     indirect( /* return indirect vertical illuminance */
279     struct glare_dir *gd
280     )
281 greg 1.9 {
282     return(gd->indirect);
283     }
284 schorsch 2.10
285    
286     static double
287     total( /* return total vertical illuminance */
288     struct glare_dir *gd
289     )
290 greg 1.9 {
291     return(direct(gd)+gd->indirect);
292 greg 1.1 }
293 schorsch 2.10
294    
295 greg 1.1 /*
296     * posindex - compute glare position index from:
297     *
298     * Source Direction
299     * View Direction
300     * View Up Direction
301     *
302     * All vectors are assumed to be normalized.
303     * This function is an implementation of the method proposed by
304     * Robert Levin in his 1975 JIES article.
305 greg 1.5 * This calculation presumes the view direction and up vectors perpendicular.
306 greg 1.1 * We return a value less than zero for improper positions.
307     */
308 schorsch 2.10
309     static double
310     posindex( /* compute position index */
311     FVECT sd,
312     FVECT vd,
313     FVECT vu
314     )
315 greg 1.1 {
316     double sigma, tau;
317     double d;
318 schorsch 2.10
319 greg 1.1 d = DOT(sd,vd);
320     if (d <= 0.0)
321     return(-1.0);
322 greg 1.5 if (d >= 1.0)
323     return(1.0);
324 greg 1.1 sigma = acos(d) * (180./PI);
325 greg 2.4 d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
326 greg 1.5 if (d >= 1.0)
327     tau = 0.0;
328     else
329     tau = acos(d) * (180./PI);
330 greg 1.1 return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
331     + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
332     ) );
333     }
334 schorsch 2.10
335    
336     static double
337     dgi( /* compute Daylight Glare Index */
338     struct glare_dir *gd
339     )
340 greg 2.4 {
341 greg 2.11 struct glare_src *gs;
342 greg 2.4 FVECT mydir,testdir[7],vhor;
343 greg 2.9 double r,posn,omega,p[7],sum;
344 greg 2.4 int i,n;
345 schorsch 2.10
346 greg 2.4 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
347     sum = 0.0; n = 0;
348     for (gs = all_srcs; gs != NULL; gs = gs->next) {
349 greg 1.1
350 greg 2.4 /* compute 1/p^2 weighted solid angle of the source */
351     r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
352     fcross(vhor,gs->dir,midview.vup);
353     normalize(vhor);
354     VCOPY(testdir[0],gs->dir);
355     fvsum(testdir[1],gs->dir,vhor,r);
356     fvsum(testdir[2],gs->dir,vhor,0.5*r);
357     fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
358     fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
359     fvsum(testdir[3],gs->dir,vhor,-r);
360     fvsum(testdir[4],gs->dir,vhor,-0.5*r);
361     fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
362     fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
363     for (i = 0; i < 7; i++) {
364     normalize(testdir[i]);
365 greg 2.9 posn = posindex(testdir[i],mydir,midview.vup);
366     if (posn <= FTINY)
367     p[i] = 0.0;
368     else
369     p[i] = 1./(posn*posn);
370 greg 2.4 }
371     r = 1-gs->dom/2./PI;
372     omega = gs->dom*p[0];
373     omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
374     omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
375     +0.3334*(p[2]+p[4]+p[5]+p[6]));
376 greg 1.1
377 schorsch 2.10 sum += pow(gs->lum,1.6) * pow(omega,0.8) /
378 greg 2.4 (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
379     n++;
380     }
381     if (n == 0)
382     return(0.0);
383     return( 10*log10(0.478*sum) );
384     }
385    
386 schorsch 2.10
387     static double
388     brs_gi( /* compute BRS Glare Index */
389     struct glare_dir *gd
390     )
391 greg 2.4 {
392 greg 2.11 struct glare_src *gs;
393 greg 2.4 FVECT mydir;
394     double p;
395     double sum;
396 schorsch 2.10
397 greg 2.4 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
398     sum = 0.0;
399     for (gs = all_srcs; gs != NULL; gs = gs->next) {
400     p = posindex(gs->dir, mydir, midview.vup);
401     if (p <= FTINY)
402     continue;
403     sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
404     }
405     if (sum <= FTINY)
406     return(0.0);
407     sum /= gd->indirect/PI;
408     return(10*log10(0.478*sum));
409     }
410    
411    
412 schorsch 2.10 static double
413     guth_dgr( /* compute Guth discomfort glare rating */
414     struct glare_dir *gd
415     )
416 greg 1.1 {
417     #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
418 greg 2.11 struct glare_src *gs;
419 greg 1.7 FVECT mydir;
420 greg 1.1 double p;
421     double sum;
422 greg 1.3 double wtot, brsum;
423 greg 1.1 int n;
424 schorsch 2.10
425 greg 1.7 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
426 greg 1.3 sum = wtot = brsum = 0.0; n = 0;
427 greg 1.1 for (gs = all_srcs; gs != NULL; gs = gs->next) {
428 greg 1.7 p = posindex(gs->dir, mydir, midview.vup);
429 greg 1.1 if (p <= FTINY)
430     continue;
431     sum += gs->lum * q(gs->dom) / p;
432 greg 1.3 brsum += gs->lum * gs->dom;
433     wtot += gs->dom;
434 greg 1.1 n++;
435     }
436     if (n == 0)
437     return(0.0);
438 greg 1.8 return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
439 greg 1.1 pow((double)n, -.0914) ) );
440     #undef q
441     }
442 schorsch 2.10
443    
444 greg 1.1 #ifndef M_SQRT2
445     #define M_SQRT2 1.41421356237309504880
446     #endif
447 schorsch 2.10
448 greg 1.1 #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
449 schorsch 2.10
450    
451     static double
452     guth_vcp( /* compute Guth visual comfort probability */
453     struct glare_dir *gd
454     )
455 greg 1.1 {
456 greg 2.5 extern double erfc();
457 greg 1.7 double dgr;
458 schorsch 2.10
459 greg 1.7 dgr = guth_dgr(gd);
460     if (dgr <= FTINY)
461     return(100.0);
462     return(100.*norm_integral(6.374-1.3227*log(dgr)));
463 greg 1.8 }
464 schorsch 2.10
465    
466     static double
467     cie_cgi( /* compute CIE Glare Index */
468     struct glare_dir *gd
469     )
470 greg 1.8 {
471 greg 2.11 struct glare_src *gs;
472 greg 1.8 FVECT mydir;
473     double dillum;
474     double p;
475     double sum;
476 schorsch 2.10
477 greg 1.8 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
478     sum = 0.0;
479     for (gs = all_srcs; gs != NULL; gs = gs->next) {
480     p = posindex(gs->dir, mydir, midview.vup);
481     if (p <= FTINY)
482     continue;
483     sum += gs->lum*gs->lum * gs->dom / (p*p);
484     }
485     if (sum <= FTINY)
486     return(0.0);
487     dillum = direct(gd);
488     return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
489 greg 2.4 }
490    
491    
492 schorsch 2.10 static double
493     ugr( /* compute Unified Glare Rating */
494     struct glare_dir *gd
495     )
496 greg 2.4 {
497 greg 2.11 struct glare_src *gs;
498 greg 2.4 FVECT mydir;
499     double p;
500     double sum;
501 schorsch 2.10
502 greg 2.4 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
503     sum = 0.0;
504     for (gs = all_srcs; gs != NULL; gs = gs->next) {
505     p = posindex(gs->dir, mydir, midview.vup);
506     if (p <= FTINY)
507     continue;
508     sum += gs->lum*gs->lum * gs->dom / (p*p);
509     }
510     if (sum <= FTINY)
511     return(0.0);
512     return(8.*log10(0.25*sum*PI/gd->indirect));
513 greg 1.1 }