ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 2.14
Committed: Mon Feb 17 19:19:45 2020 UTC (4 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.13: +31 -25 lines
Log Message:
Hopeful fix to dgi calculation issues

File Contents

# User Rev Content
1 greg 1.1 #ifndef lint
2 greg 2.14 static const char RCSid[] = "$Id: glarendx.c,v 2.13 2018/08/02 18:33:50 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 greg 2.14 char *err;
168 greg 1.9 /* read header */
169     getheader(stdin, headline, NULL);
170     if (wrongformat) {
171     fprintf(stderr, "%s: bad input format\n", progname);
172     exit(1);
173     }
174     /* set view */
175 greg 2.14 if ((err = setview(&midview)) != NULL) {
176     fprintf(stderr, "%s: %s\n",
177     progname, err);
178     exit(1);
179     }
180     if (fabs(DOT(midview.vdir, midview.vup)) > 1e-4) {
181     fprintf(stderr,
182     "%s: Input view direction not perpendicular to up vector\n",
183 greg 1.9 progname);
184     exit(1);
185     }
186     }
187 schorsch 2.10
188    
189     static void
190     read_input(void) /* read glare sources from stdin */
191 greg 1.1 {
192     #define S_SEARCH 0
193     #define S_SOURCE 1
194     #define S_DIREC 2
195     int state = S_SEARCH;
196     char buf[128];
197 greg 2.11 struct glare_src *gs;
198     struct glare_dir *gd;
199 schorsch 2.10
200 greg 1.1 while (fgets(buf, sizeof(buf), stdin) != NULL)
201     switch (state) {
202     case S_SEARCH:
203     if (!strcmp(buf, "BEGIN glare source\n"))
204     state = S_SOURCE;
205     else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
206     state = S_DIREC;
207     break;
208     case S_SOURCE:
209     if (!strncmp(buf, "END", 3)) {
210     state = S_SEARCH;
211     break;
212     }
213     if ((gs = newp(struct glare_src)) == NULL)
214     goto memerr;
215 greg 2.11 if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
216     &gs->dir[2]) != 3 ||
217 greg 2.12 sscanf(sskip2(buf, 3), "%lf %lf",
218 greg 2.11 &gs->dom, &gs->lum) != 2)
219 greg 1.1 goto readerr;
220 greg 1.2 normalize(gs->dir);
221 greg 1.1 gs->next = all_srcs;
222     all_srcs = gs;
223     break;
224     case S_DIREC:
225     if (!strncmp(buf, "END", 3)) {
226     state = S_SEARCH;
227     break;
228     }
229     if ((gd = newp(struct glare_dir)) == NULL)
230     goto memerr;
231     if (sscanf(buf, "%lf %lf",
232     &gd->ang, &gd->indirect) != 2)
233     goto readerr;
234     gd->ang *= PI/180.0; /* convert to radians */
235     gd->next = all_dirs;
236     all_dirs = gd;
237     break;
238     }
239     return;
240     memerr:
241     perror(progname);
242     exit(1);
243     readerr:
244     fprintf(stderr, "%s: read error on input\n", progname);
245     exit(1);
246     #undef S_SEARCH
247     #undef S_SOURCE
248     #undef S_DIREC
249     }
250 schorsch 2.10
251    
252     static void
253     print_values( /* print out calculations */
254     gdfun *func
255     )
256 greg 1.1 {
257 greg 2.11 struct glare_dir *gd;
258 schorsch 2.10
259 greg 1.1 for (gd = all_dirs; gd != NULL; gd = gd->next)
260 schorsch 2.10 printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
261 greg 1.1 }
262 schorsch 2.10
263    
264     static double
265     direct( /* compute direct vertical illuminance */
266     struct glare_dir *gd
267     )
268 greg 1.1 {
269     FVECT mydir;
270     double d, dval;
271 greg 2.11 struct glare_src *gs;
272 schorsch 2.10
273 greg 1.1 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
274     dval = 0.0;
275     for (gs = all_srcs; gs != NULL; gs = gs->next) {
276     d = DOT(mydir,gs->dir);
277     if (d > FTINY)
278     dval += d * gs->dom * gs->lum;
279     }
280     return(dval);
281 greg 1.9 }
282 schorsch 2.10
283    
284     static double
285     indirect( /* return indirect vertical illuminance */
286     struct glare_dir *gd
287     )
288 greg 1.9 {
289     return(gd->indirect);
290     }
291 schorsch 2.10
292    
293     static double
294     total( /* return total vertical illuminance */
295     struct glare_dir *gd
296     )
297 greg 1.9 {
298     return(direct(gd)+gd->indirect);
299 greg 1.1 }
300 schorsch 2.10
301    
302 greg 1.1 /*
303     * posindex - compute glare position index from:
304     *
305     * Source Direction
306     * View Direction
307     * View Up Direction
308     *
309     * All vectors are assumed to be normalized.
310     * This function is an implementation of the method proposed by
311     * Robert Levin in his 1975 JIES article.
312 greg 1.5 * This calculation presumes the view direction and up vectors perpendicular.
313 greg 1.1 * We return a value less than zero for improper positions.
314     */
315 schorsch 2.10
316     static double
317     posindex( /* compute position index */
318     FVECT sd,
319     FVECT vd,
320     FVECT vu
321     )
322 greg 1.1 {
323     double sigma, tau;
324     double d;
325 schorsch 2.10
326 greg 1.1 d = DOT(sd,vd);
327     if (d <= 0.0)
328     return(-1.0);
329 greg 2.14 if (d >= 1.0-FTINY)
330 greg 1.5 return(1.0);
331 greg 1.1 sigma = acos(d) * (180./PI);
332 greg 2.4 d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
333 greg 2.14 tau = Acos(d) * (180./PI);
334 greg 1.1 return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
335     + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
336     ) );
337     }
338 schorsch 2.10
339    
340     static double
341     dgi( /* compute Daylight Glare Index */
342     struct glare_dir *gd
343     )
344 greg 2.4 {
345 greg 2.11 struct glare_src *gs;
346 greg 2.4 FVECT mydir,testdir[7],vhor;
347 greg 2.9 double r,posn,omega,p[7],sum;
348 greg 2.14 int i;
349 schorsch 2.10
350 greg 2.4 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
351 greg 2.14 sum = 0.0;
352 greg 2.4 for (gs = all_srcs; gs != NULL; gs = gs->next) {
353 greg 1.1
354 greg 2.4 /* compute 1/p^2 weighted solid angle of the source */
355 greg 2.14 r = 1. - gs->dom*(.5/PI);
356     r = sqrt(1. - r*r);
357 greg 2.4 fcross(vhor,gs->dir,midview.vup);
358 greg 2.14 if (normalize(vhor) == 0.)
359     continue;
360 greg 2.4 VCOPY(testdir[0],gs->dir);
361     fvsum(testdir[1],gs->dir,vhor,r);
362     fvsum(testdir[2],gs->dir,vhor,0.5*r);
363     fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
364     fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
365     fvsum(testdir[3],gs->dir,vhor,-r);
366     fvsum(testdir[4],gs->dir,vhor,-0.5*r);
367     fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
368     fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
369     for (i = 0; i < 7; i++) {
370 greg 2.14 p[i] = 0.0;
371     if (normalize(testdir[i]) == 0.)
372     continue;
373 greg 2.9 posn = posindex(testdir[i],mydir,midview.vup);
374 greg 2.14 if (posn > FTINY)
375 greg 2.9 p[i] = 1./(posn*posn);
376 greg 2.4 }
377 greg 2.14 r = 1. - gs->dom*(.5/PI);
378 greg 2.4 omega = gs->dom*p[0];
379 greg 2.14 omega += (r*PI*(1.+1./(r*r))-2.*PI)*(-p[0]+(p[1]+p[2])*0.5);
380     omega += (2.*PI-r*PI*(1.+1./(r*r)))*(-p[0]-0.1667*(p[1]+p[3])
381 greg 2.4 +0.3334*(p[2]+p[4]+p[5]+p[6]));
382 greg 2.14 if (omega <= 0.)
383     omega = 0.;
384 schorsch 2.10 sum += pow(gs->lum,1.6) * pow(omega,0.8) /
385 greg 2.14 (gd->indirect*(1./PI) + 0.07*sqrt(gs->dom)*gs->lum);
386 greg 2.4 }
387 greg 2.14 if (sum <= FTINY)
388 greg 2.4 return(0.0);
389 greg 2.14 return( 10.*log10(0.478*sum) );
390 greg 2.4 }
391    
392 schorsch 2.10
393     static double
394     brs_gi( /* compute BRS Glare Index */
395     struct glare_dir *gd
396     )
397 greg 2.4 {
398 greg 2.11 struct glare_src *gs;
399 greg 2.4 FVECT mydir;
400     double p;
401     double sum;
402 schorsch 2.10
403 greg 2.4 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
404     sum = 0.0;
405     for (gs = all_srcs; gs != NULL; gs = gs->next) {
406     p = posindex(gs->dir, mydir, midview.vup);
407     if (p <= FTINY)
408     continue;
409     sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
410     }
411     if (sum <= FTINY)
412     return(0.0);
413 greg 2.14 sum *= PI/gd->indirect;
414 greg 2.4 return(10*log10(0.478*sum));
415     }
416    
417    
418 schorsch 2.10 static double
419     guth_dgr( /* compute Guth discomfort glare rating */
420     struct glare_dir *gd
421     )
422 greg 1.1 {
423     #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
424 greg 2.11 struct glare_src *gs;
425 greg 1.7 FVECT mydir;
426 greg 1.1 double p;
427     double sum;
428 greg 1.3 double wtot, brsum;
429 greg 1.1 int n;
430 schorsch 2.10
431 greg 1.7 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
432 greg 1.3 sum = wtot = brsum = 0.0; n = 0;
433 greg 1.1 for (gs = all_srcs; gs != NULL; gs = gs->next) {
434 greg 1.7 p = posindex(gs->dir, mydir, midview.vup);
435 greg 1.1 if (p <= FTINY)
436     continue;
437     sum += gs->lum * q(gs->dom) / p;
438 greg 1.3 brsum += gs->lum * gs->dom;
439     wtot += gs->dom;
440 greg 1.1 n++;
441     }
442     if (n == 0)
443     return(0.0);
444 greg 1.8 return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
445 greg 1.1 pow((double)n, -.0914) ) );
446     #undef q
447     }
448 schorsch 2.10
449    
450 greg 1.1 #ifndef M_SQRT2
451     #define M_SQRT2 1.41421356237309504880
452     #endif
453 schorsch 2.10
454 greg 1.1 #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
455 schorsch 2.10
456    
457     static double
458     guth_vcp( /* compute Guth visual comfort probability */
459     struct glare_dir *gd
460     )
461 greg 1.1 {
462 greg 2.5 extern double erfc();
463 greg 1.7 double dgr;
464 schorsch 2.10
465 greg 1.7 dgr = guth_dgr(gd);
466     if (dgr <= FTINY)
467     return(100.0);
468     return(100.*norm_integral(6.374-1.3227*log(dgr)));
469 greg 1.8 }
470 schorsch 2.10
471    
472     static double
473     cie_cgi( /* compute CIE Glare Index */
474     struct glare_dir *gd
475     )
476 greg 1.8 {
477 greg 2.11 struct glare_src *gs;
478 greg 1.8 FVECT mydir;
479     double dillum;
480     double p;
481     double sum;
482 schorsch 2.10
483 greg 1.8 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
484     sum = 0.0;
485     for (gs = all_srcs; gs != NULL; gs = gs->next) {
486     p = posindex(gs->dir, mydir, midview.vup);
487     if (p <= FTINY)
488     continue;
489     sum += gs->lum*gs->lum * gs->dom / (p*p);
490     }
491     if (sum <= FTINY)
492     return(0.0);
493     dillum = direct(gd);
494     return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
495 greg 2.4 }
496    
497    
498 schorsch 2.10 static double
499     ugr( /* compute Unified Glare Rating */
500     struct glare_dir *gd
501     )
502 greg 2.4 {
503 greg 2.11 struct glare_src *gs;
504 greg 2.4 FVECT mydir;
505     double p;
506     double sum;
507 schorsch 2.10
508 greg 2.4 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
509     sum = 0.0;
510     for (gs = all_srcs; gs != NULL; gs = gs->next) {
511     p = posindex(gs->dir, mydir, midview.vup);
512     if (p <= FTINY)
513     continue;
514     sum += gs->lum*gs->lum * gs->dom / (p*p);
515     }
516     if (sum <= FTINY)
517     return(0.0);
518     return(8.*log10(0.25*sum*PI/gd->indirect));
519 greg 1.1 }