ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
Revision: 1.8
Committed: Thu May 2 13:17:41 1991 UTC (32 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 1.7: +30 -4 lines
Log Message:
added CGI (Einhorn) formula

File Contents

# 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     * Compute Glare Index given by program name:
9     *
10 greg 1.2 * guth_dgr - Guth discomfort glare rating
11     * guth_vcp - Guth visual comfort probability
12 greg 1.8 * cie_cgi - CIE Glare Index (1983)
13 greg 1.1 *
14     * 12 April 1991 Greg Ward EPFL
15     */
16    
17     #include "standard.h"
18     #include "view.h"
19    
20     extern double erfc();
21    
22     double posindex();
23     int headline();
24    
25 greg 1.8 double direct(), guth_dgr(), guth_vcp(), cie_cgi();
26 greg 1.1
27     struct named_func {
28     char *name;
29     double (*func)();
30     } all_funcs[] = {
31     {"direct", direct},
32 greg 1.2 {"guth_dgr", guth_dgr},
33     {"guth_vcp", guth_vcp},
34 greg 1.8 {"cie_cgi", cie_cgi},
35 greg 1.1 {NULL}
36     };
37    
38     struct glare_src {
39     FVECT dir; /* source direction */
40     double dom; /* solid angle */
41     double lum; /* average luminance */
42     struct glare_src *next;
43     } *all_srcs = NULL;
44    
45     struct glare_dir {
46     double ang; /* angle (in radians) */
47     double indirect; /* indirect illuminance */
48     struct glare_dir *next;
49     } *all_dirs = NULL;
50    
51     #define newp(type) (type *)malloc(sizeof(type))
52    
53     char *progname;
54     int print_header = 1;
55    
56     VIEW midview = STDVIEW;
57    
58 greg 1.4 int wrongformat = 0;
59 greg 1.1
60 greg 1.4
61 greg 1.1 main(argc, argv)
62     int argc;
63     char *argv[];
64     {
65     extern char *rindex();
66     struct named_func *funp;
67     char *progtail;
68     int i;
69     /* get program name */
70     progname = argv[0];
71     progtail = rindex(progname, '/'); /* final component */
72     if (progtail == NULL)
73     progtail = progname;
74     else
75     progtail++;
76     /* get options */
77     for (i = 1; i < argc && argv[i][0] == '-'; i++)
78     switch (argv[i][1]) {
79     case 't':
80     progtail = argv[++i];
81     break;
82     case 'h':
83     print_header = 0;
84     break;
85     default:
86     goto userr;
87     }
88     if (i < argc-1)
89     goto userr;
90     if (i == argc-1) /* open file */
91     if (freopen(argv[i], "r", stdin) == NULL) {
92     perror(argv[i]);
93     exit(1);
94     }
95     /* read header */
96 greg 1.4 getheader(stdin, headline, NULL);
97     if (wrongformat) {
98     fprintf(stderr, "%s: bad input format\n", progname);
99     exit(1);
100     }
101 greg 1.1 if (print_header) { /* add to header */
102     printargs(i, argv, stdout);
103     putchar('\n');
104     }
105     /* set view */
106     if (setview(&midview) != NULL) {
107 greg 1.4 fprintf(stderr, "%s: bad view information in input\n",
108     progname);
109 greg 1.1 exit(1);
110     }
111     /* get findglare data */
112     read_input();
113     /* find calculation */
114     for (funp = all_funcs; funp->name != NULL; funp++)
115     if (!strcmp(funp->name, progtail)) {
116     print_values(funp->func);
117     exit(0); /* we're done */
118     }
119     /* invalid function */
120     fprintf(stderr, "%s: unknown function!\n", progtail);
121     exit(1);
122     userr:
123     fprintf(stderr, "Usage: %s [-t type][-h] [input]\n", progname);
124     exit(1);
125     }
126    
127    
128     headline(s) /* get line from header */
129     char *s;
130     {
131 greg 1.4 char fmt[32];
132    
133 greg 1.1 if (print_header) /* copy to output */
134     fputs(s, stdout);
135     if (!strncmp(s, VIEWSTR, VIEWSTRL))
136     sscanview(&midview, s+VIEWSTRL);
137 greg 1.4 else if (isformat(s)) {
138     formatval(fmt, s);
139 greg 1.6 wrongformat = strcmp(fmt, "ascii");
140 greg 1.4 }
141 greg 1.1 }
142    
143    
144     read_input() /* read glare sources from stdin */
145     {
146     #define S_SEARCH 0
147     #define S_SOURCE 1
148     #define S_DIREC 2
149     int state = S_SEARCH;
150     char buf[128];
151     register struct glare_src *gs;
152     register struct glare_dir *gd;
153    
154     while (fgets(buf, sizeof(buf), stdin) != NULL)
155     switch (state) {
156     case S_SEARCH:
157     if (!strcmp(buf, "BEGIN glare source\n"))
158     state = S_SOURCE;
159     else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
160     state = S_DIREC;
161     break;
162     case S_SOURCE:
163     if (!strncmp(buf, "END", 3)) {
164     state = S_SEARCH;
165     break;
166     }
167     if ((gs = newp(struct glare_src)) == NULL)
168     goto memerr;
169     if (sscanf(buf, "%lf %lf %lf %lf %lf",
170     &gs->dir[0], &gs->dir[1], &gs->dir[2],
171     &gs->dom, &gs->lum) != 5)
172     goto readerr;
173 greg 1.2 normalize(gs->dir);
174 greg 1.1 gs->next = all_srcs;
175     all_srcs = gs;
176     break;
177     case S_DIREC:
178     if (!strncmp(buf, "END", 3)) {
179     state = S_SEARCH;
180     break;
181     }
182     if ((gd = newp(struct glare_dir)) == NULL)
183     goto memerr;
184     if (sscanf(buf, "%lf %lf",
185     &gd->ang, &gd->indirect) != 2)
186     goto readerr;
187     gd->ang *= PI/180.0; /* convert to radians */
188     gd->next = all_dirs;
189     all_dirs = gd;
190     break;
191     }
192     return;
193     memerr:
194     perror(progname);
195     exit(1);
196     readerr:
197     fprintf(stderr, "%s: read error on input\n", progname);
198     exit(1);
199     #undef S_SEARCH
200     #undef S_SOURCE
201     #undef S_DIREC
202     }
203    
204    
205     print_values(funp) /* print out calculations */
206     double (*funp)();
207     {
208     register struct glare_dir *gd;
209    
210     for (gd = all_dirs; gd != NULL; gd = gd->next)
211     printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
212     }
213    
214    
215     double
216     direct(gd) /* compute direct illuminance */
217     struct glare_dir *gd;
218     {
219     FVECT mydir;
220     double d, dval;
221     register struct glare_src *gs;
222    
223     spinvector(mydir, midview.vdir, midview.vup, gd->ang);
224     dval = 0.0;
225     for (gs = all_srcs; gs != NULL; gs = gs->next) {
226     d = DOT(mydir,gs->dir);
227     if (d > FTINY)
228     dval += d * gs->dom * gs->lum;
229     }
230     return(dval);
231     }
232    
233    
234     /*
235     * posindex - compute glare position index from:
236     *
237     * Source Direction
238     * View Direction
239     * View Up Direction
240     *
241     * All vectors are assumed to be normalized.
242     * This function is an implementation of the method proposed by
243     * Robert Levin in his 1975 JIES article.
244 greg 1.5 * This calculation presumes the view direction and up vectors perpendicular.
245 greg 1.1 * We return a value less than zero for improper positions.
246     */
247    
248     double
249     posindex(sd, vd, vu) /* compute position index */
250     FVECT sd, vd, vu;
251     {
252     double sigma, tau;
253     double d;
254    
255     d = DOT(sd,vd);
256     if (d <= 0.0)
257     return(-1.0);
258 greg 1.5 if (d >= 1.0)
259     return(1.0);
260 greg 1.1 sigma = acos(d) * (180./PI);
261 greg 1.5 d = DOT(sd,vu)/sqrt(1.0-d*d);
262     if (d >= 1.0)
263     tau = 0.0;
264     else
265     tau = acos(d) * (180./PI);
266 greg 1.1 return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
267     + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
268     ) );
269     }
270    
271    
272     double
273 greg 1.2 guth_dgr(gd) /* compute Guth discomfort glare rating */
274 greg 1.1 struct glare_dir *gd;
275     {
276     #define q(w) (20.4*w+1.52*pow(w,.2)-.075)
277     register struct glare_src *gs;
278 greg 1.7 FVECT mydir;
279 greg 1.1 double p;
280     double sum;
281 greg 1.3 double wtot, brsum;
282 greg 1.1 int n;
283    
284 greg 1.7 spinvector(mydir, midview.vdir, midview.vup, gd->ang);
285 greg 1.3 sum = wtot = brsum = 0.0; n = 0;
286 greg 1.1 for (gs = all_srcs; gs != NULL; gs = gs->next) {
287 greg 1.7 p = posindex(gs->dir, mydir, midview.vup);
288 greg 1.1 if (p <= FTINY)
289     continue;
290     sum += gs->lum * q(gs->dom) / p;
291 greg 1.3 brsum += gs->lum * gs->dom;
292     wtot += gs->dom;
293 greg 1.1 n++;
294     }
295     if (n == 0)
296     return(0.0);
297 greg 1.8
298     return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
299 greg 1.1 pow((double)n, -.0914) ) );
300     #undef q
301     }
302    
303    
304     extern double erf(), erfc();
305    
306     #ifndef M_SQRT2
307     #define M_SQRT2 1.41421356237309504880
308     #endif
309    
310     #define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
311    
312    
313     double
314 greg 1.2 guth_vcp(gd) /* compute Guth visual comfort probability */
315 greg 1.1 struct glare_dir *gd;
316     {
317 greg 1.7 double dgr;
318    
319     dgr = guth_dgr(gd);
320     if (dgr <= FTINY)
321     return(100.0);
322     return(100.*norm_integral(6.374-1.3227*log(dgr)));
323 greg 1.8 }
324    
325    
326     double
327     cie_cgi(gd) /* compute CIE Glare Index */
328     struct glare_dir *gd;
329     {
330     register struct glare_src *gs;
331     FVECT mydir;
332     double dillum;
333     double p;
334     double sum;
335    
336     spinvector(mydir, midview.vdir, midview.vup, gd->ang);
337     sum = 0.0;
338     for (gs = all_srcs; gs != NULL; gs = gs->next) {
339     p = posindex(gs->dir, mydir, midview.vup);
340     if (p <= FTINY)
341     continue;
342     sum += gs->lum*gs->lum * gs->dom / (p*p);
343     }
344     if (sum <= FTINY)
345     return(0.0);
346     dillum = direct(gd);
347     return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
348 greg 1.1 }