ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/glarendx.c
(Generate patch)

Comparing ray/src/util/glarendx.c (file contents):
Revision 2.5 by greg, Sun Jun 4 10:58:38 1995 UTC vs.
Revision 2.15 by greg, Mon Feb 17 19:35:24 2020 UTC

# Line 1 | Line 1
1 – /* Copyright (c) 1991 Regents of the University of California */
2 –
1   #ifndef lint
2 < static char SCCSid[] = "$SunId$ LBL";
2 > static const char       RCSid[] = "$Id$";
3   #endif
6 –
4   /*
5   * Compute Glare Index given by program name or -t option:
6   *
# Line 21 | Line 18 | static char SCCSid[] = "$SunId$ LBL";
18   *              12 April 1991   Greg Ward       EPFL
19   *              19 April 1993   R. Compagnon    EPFL (added dgi, brs_gi, ugr)
20   */
21 <
21 >
22 > #include <string.h>
23 >
24   #include "standard.h"
25   #include "view.h"
26 <
27 <
28 < double  posindex();
29 <
30 < double  direct(), guth_dgr(), guth_vcp(), cie_cgi(),
31 <        indirect(), total(), dgi(), brs_gi(), ugr();
32 <
26 >
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   struct named_func {
61          char    *name;
62 <        double  (*func)();
62 >        gdfun   *func;
63          char    *descrip;
64   } all_funcs[] = {
65          {"dgi", dgi, "Daylight Glare Index"},
# Line 47 | Line 73 | struct named_func {
73          {"vert_ind", indirect, "Indirect Vertical Illuminance"},
74          {NULL}
75   };
76 <
51 < struct glare_src {
52 <        FVECT   dir;            /* source direction */
53 <        double  dom;            /* solid angle */
54 <        double  lum;            /* average luminance */
55 <        struct glare_src        *next;
56 < } *all_srcs = NULL;
57 <
58 < struct glare_dir {
59 <        double  ang;            /* angle (in radians) */
60 <        double  indirect;       /* indirect illuminance */
61 <        struct glare_dir        *next;
62 < } *all_dirs = NULL;
63 <
76 >
77   #define newp(type)      (type *)malloc(sizeof(type))
78 <
78 >
79   char    *progname;
80   int     print_header = 1;
81 <
81 >
82   VIEW    midview = STDVIEW;
83 <
83 >
84   int     wrongformat = 0;
85 <
86 <
87 < main(argc, argv)
88 < int     argc;
89 < char    *argv[];
85 >
86 >
87 > int
88 > main(
89 >        int     argc,
90 >        char    *argv[]
91 > )
92   {
78 –        extern char     *rindex();
93          struct named_func       *funp;
94          char    *progtail;
95          int     i;
96                                          /* get program name */
97          progname = argv[0];
98 <        progtail = rindex(progname, '/');       /* final component */
98 >        progtail = strrchr(progname, '/');      /* final component */
99          if (progtail == NULL)
100                  progtail = progname;
101          else
# Line 125 | Line 139 | userr:
139                  fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
140          exit(1);
141   }
142 <
143 <
144 < headline(s)                     /* get line from header */
145 < char    *s;
142 >
143 >
144 > static int
145 > headline(                       /* get line from header */
146 >        char    *s,
147 >        void    *p
148 > )
149   {
150 <        char    fmt[32];
151 <
150 >        char    fmt[MAXFMTLEN];
151 >
152          if (print_header)               /* copy to output */
153                  fputs(s, stdout);
154          if (isview(s))
# Line 140 | Line 157 | char   *s;
157                  formatval(fmt, s);
158                  wrongformat = strcmp(fmt, "ascii");
159          }
160 +        return(0);
161   }
162 <
163 <
164 < init()                          /* initialize calculation */
162 >
163 >
164 > static void
165 > init(void)                              /* initialize calculation */
166   {
167 +        char    *err;
168                                          /* read header */
169          getheader(stdin, headline, NULL);
170          if (wrongformat) {
# Line 152 | Line 172 | init()                         /* initialize calculation */
172                  exit(1);
173          }
174                                          /* set view */
175 <        if (setview(&midview) != NULL) {
176 <                fprintf(stderr, "%s: bad view information in input\n",
175 >        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                                  progname);
184                  exit(1);
185          }
186   }
187 <
188 <
189 < read_input()                    /* read glare sources from stdin */
187 >
188 >
189 > static void
190 > read_input(void)                        /* read glare sources from stdin */
191   {
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 <        register struct glare_src       *gs;
198 <        register struct glare_dir       *gd;
199 <
197 >        struct glare_src        *gs;
198 >        struct glare_dir        *gd;
199 >
200          while (fgets(buf, sizeof(buf), stdin) != NULL)
201                  switch (state) {
202                  case S_SEARCH:
# Line 185 | Line 212 | read_input()                   /* read glare sources from stdin */
212                          }
213                          if ((gs = newp(struct glare_src)) == NULL)
214                                  goto memerr;
215 <                        if (sscanf(buf, "%lf %lf %lf %lf %lf",
216 <                                        &gs->dir[0], &gs->dir[1], &gs->dir[2],
217 <                                        &gs->dom, &gs->lum) != 5)
215 >                        if (sscanf(buf, FVFORMAT, &gs->dir[0], &gs->dir[1],
216 >                                                &gs->dir[2]) != 3 ||
217 >                                        sscanf(sskip2(buf, 3), "%lf %lf",
218 >                                                &gs->dom, &gs->lum) != 2)
219                                  goto readerr;
220                          normalize(gs->dir);
221                          gs->next = all_srcs;
# Line 219 | Line 247 | readerr:
247   #undef S_SOURCE
248   #undef S_DIREC
249   }
250 <
251 <
252 < print_values(funp)              /* print out calculations */
253 < double  (*funp)();
250 >
251 >
252 > static void
253 > print_values(           /* print out calculations */
254 >        gdfun *func
255 > )
256   {
257 <        register struct glare_dir       *gd;
258 <
257 >        struct glare_dir        *gd;
258 >
259          for (gd = all_dirs; gd != NULL; gd = gd->next)
260 <                printf("%f\t%f\n", gd->ang*(180.0/PI), (*funp)(gd));
260 >                printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
261   }
262 <
263 <
264 < double
265 < direct(gd)                      /* compute direct vertical illuminance */
266 < struct glare_dir        *gd;
262 >
263 >
264 > static double
265 > direct(                 /* compute direct vertical illuminance */
266 >        struct glare_dir        *gd
267 > )
268   {
269          FVECT   mydir;
270          double  d, dval;
271 <        register struct glare_src       *gs;
272 <
271 >        struct glare_src        *gs;
272 >
273          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
274          dval = 0.0;
275          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 248 | Line 279 | struct glare_dir       *gd;
279          }
280          return(dval);
281   }
282 <
283 <
284 < double
285 < indirect(gd)                    /* return indirect vertical illuminance */
286 < struct glare_dir        *gd;
282 >
283 >
284 > static double
285 > indirect(                       /* return indirect vertical illuminance */
286 >        struct glare_dir        *gd
287 > )
288   {
289          return(gd->indirect);
290   }
291 <
292 <
293 < double
294 < total(gd)                       /* return total vertical illuminance */
295 < struct glare_dir        *gd;
291 >
292 >
293 > static double
294 > total(                  /* return total vertical illuminance */
295 >        struct glare_dir        *gd
296 > )
297   {
298          return(direct(gd)+gd->indirect);
299   }
300 <
301 <
300 >
301 >
302   /*
303   * posindex -   compute glare position index from:
304   *
# Line 279 | Line 312 | struct glare_dir       *gd;
312   * This calculation presumes the view direction and up vectors perpendicular.
313   * We return a value less than zero for improper positions.
314   */
315 <
316 < double
317 < posindex(sd, vd, vu)                    /* compute position index */
318 < FVECT   sd, vd, vu;
315 >
316 > static double
317 > posindex(                       /* compute position index */
318 >        FVECT   sd,
319 >        FVECT   vd,
320 >        FVECT   vu
321 > )
322   {
323          double  sigma, tau;
324          double  d;
325 <
325 >
326          d = DOT(sd,vd);
327          if (d <= 0.0)
328                  return(-1.0);
329 <        if (d >= 1.0)
329 >        if (d >= 1.0-FTINY)
330                  return(1.0);
331          sigma = acos(d) * (180./PI);
332          d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
333 <        if (d >= 1.0)
298 <                tau = 0.0;
299 <        else
300 <                tau = acos(d) * (180./PI);
333 >        tau = Acos(d) * (180./PI);
334          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 <
339 <
340 < double
341 < dgi(gd)         /* compute Daylight Glare Index */
342 < struct glare_dir        *gd;
338 >
339 >
340 > static double
341 > dgi(            /* compute Daylight Glare Index */
342 >        struct glare_dir        *gd
343 > )
344   {
345 <        register struct glare_src       *gs;
345 >        struct glare_src        *gs;
346          FVECT   mydir,testdir[7],vhor;
347 <        double  r,omega,p[7],sum;
348 <        int     i,n;
349 <
347 >        double  r,posn,omega,p[7],sum;
348 >        int     i;
349 >
350          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
351 <        sum = 0.0; n = 0;
351 >        sum = 0.0;
352          for (gs = all_srcs; gs != NULL; gs = gs->next) {
353  
354                  /* compute 1/p^2 weighted solid angle of the source */
355 <                r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
355 >                r = 1. - gs->dom*(.5/PI);
356 >                r = sqrt(1. - r*r);
357                  fcross(vhor,gs->dir,midview.vup);
358 <                normalize(vhor);
358 >                if (normalize(vhor) == 0.)
359 >                        continue;
360                  VCOPY(testdir[0],gs->dir);
361                  fvsum(testdir[1],gs->dir,vhor,r);
362                  fvsum(testdir[2],gs->dir,vhor,0.5*r);
# Line 331 | Line 367 | struct glare_dir       *gd;
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 <                        normalize(testdir[i]);
371 <                        p[i] = pow(posindex(testdir[i],mydir,midview.vup),-2.0);
372 <                        if (p[i] <= FTINY) p[i] = 0.0;
370 >                        p[i] = 0.0;
371 >                        if (normalize(testdir[i]) == 0.)
372 >                                continue;
373 >                        posn = posindex(testdir[i],mydir,midview.vup);
374 >                        if (posn > FTINY)
375 >                                p[i] = 1./(posn*posn);
376                  }
377 <                r = 1-gs->dom/2./PI;
377 >                r = 1. - gs->dom*(.5/PI);
378                  omega = gs->dom*p[0];
379 <                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])
379 >                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                            +0.3334*(p[2]+p[4]+p[5]+p[6]));
382 <
383 <                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
384 <                       (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
385 <                n++;
382 >                if (omega <= 0.)
383 >                        continue;
384 >                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
385 >                       (gd->indirect*(1./PI) + 0.07*sqrt(gs->dom)*gs->lum);
386          }
387 <        if (n == 0)
387 >        if (sum <= FTINY)
388                  return(0.0);
389 <        return( 10*log10(0.478*sum) );
389 >        return( 10.*log10(0.478*sum) );
390   }
352 –
391  
392 < double
393 < brs_gi(gd)              /* compute BRS Glare Index */
394 < struct glare_dir        *gd;
392 >
393 > static double
394 > brs_gi(         /* compute BRS Glare Index */
395 >        struct glare_dir        *gd
396 > )
397   {
398 <        register struct glare_src       *gs;
398 >        struct glare_src        *gs;
399          FVECT   mydir;
400          double  p;
401          double  sum;
402 <
402 >
403          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
404          sum = 0.0;
405          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 370 | Line 410 | struct glare_dir       *gd;
410          }
411          if (sum <= FTINY)
412                  return(0.0);
413 <        sum /= gd->indirect/PI;
413 >        sum *= PI/gd->indirect;
414          return(10*log10(0.478*sum));
415   }
416  
417  
418 < double
419 < guth_dgr(gd)            /* compute Guth discomfort glare rating */
420 < struct glare_dir        *gd;
418 > static double
419 > guth_dgr(               /* compute Guth discomfort glare rating */
420 >        struct glare_dir        *gd
421 > )
422   {
423   #define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
424 <        register struct glare_src       *gs;
424 >        struct glare_src        *gs;
425          FVECT   mydir;
426          double  p;
427          double  sum;
428          double  wtot, brsum;
429          int     n;
430 <
430 >
431          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
432          sum = wtot = brsum = 0.0; n = 0;
433          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 404 | Line 445 | struct glare_dir       *gd;
445                          pow((double)n, -.0914) ) );
446   #undef q
447   }
448 <
449 <
448 >
449 >
450   #ifndef M_SQRT2
451   #define M_SQRT2 1.41421356237309504880
452   #endif
453 <
453 >
454   #define norm_integral(z)        (1.-.5*erfc((z)/M_SQRT2))
455 <
456 <
457 < double
458 < guth_vcp(gd)            /* compute Guth visual comfort probability */
459 < struct glare_dir        *gd;
455 >
456 >
457 > static double
458 > guth_vcp(               /* compute Guth visual comfort probability */
459 >        struct glare_dir        *gd
460 > )
461   {
462          extern double   erfc();
463          double  dgr;
464 <
464 >
465          dgr = guth_dgr(gd);
466          if (dgr <= FTINY)
467                  return(100.0);
468          return(100.*norm_integral(6.374-1.3227*log(dgr)));
469   }
470 <
471 <
472 < double
473 < cie_cgi(gd)             /* compute CIE Glare Index */
474 < struct glare_dir        *gd;
470 >
471 >
472 > static double
473 > cie_cgi(                /* compute CIE Glare Index */
474 >        struct glare_dir        *gd
475 > )
476   {
477 <        register struct glare_src       *gs;
477 >        struct glare_src        *gs;
478          FVECT   mydir;
479          double  dillum;
480          double  p;
481          double  sum;
482 <
482 >
483          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
484          sum = 0.0;
485          for (gs = all_srcs; gs != NULL; gs = gs->next) {
# Line 452 | Line 495 | struct glare_dir       *gd;
495   }
496  
497  
498 < double
499 < ugr(gd)         /* compute Unified Glare Rating */
500 < struct glare_dir        *gd;
498 > static double
499 > ugr(            /* compute Unified Glare Rating */
500 >        struct glare_dir        *gd
501 > )
502   {
503 <        register struct glare_src       *gs;
503 >        struct glare_src        *gs;
504          FVECT   mydir;
505          double  p;
506          double  sum;
507 <
507 >
508          spinvector(mydir, midview.vdir, midview.vup, gd->ang);
509          sum = 0.0;
510          for (gs = all_srcs; gs != NULL; gs = gs->next) {

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines