ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.2
Committed: Tue Aug 18 15:02:53 2015 UTC (8 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R0
Changes since 2.1: +3 -0 lines
Log Message: 
Added RCCS id's to new source files (forgotten during initial check-in)

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id$";
3 #endif
4 /* EVALGLARE V1.17
5 * Evalglare Software License, Version 1.0
6 *
7 * Copyright (c) 1995 - 2015 Fraunhofer ISE, EPFL.
8 * All rights reserved.
9 *
10 *
11 * Redistribution and use in source and binary forms, WITHOUT
12 * modification, are permitted provided that the following all conditions
13 * are met:
14 *
15 * 1. Redistributions of source code must retain the above copyright
16 * notice, this list of conditions and the following disclaimer.
17 *
18 * 2. Redistributions in binary form must reproduce the above copyright
19 * notice, this list of conditions and the following disclaimer in
20 * the documentation and/or other materials provided with the
21 * distribution.
22 *
23 * 3. The end-user documentation included with the redistribution,
24 * if any, must include the following acknowledgments:
25 * "This product includes the evalglare software,
26 developed at Fraunhofer ISE by J. Wienold" and
27 * "This product includes Radiance software
28 * (http://radsite.lbl.gov/)
29 * developed by the Lawrence Berkeley National Laboratory
30 * (http://www.lbl.gov/)."
31 * Alternately, this acknowledgment may appear in the software itself,
32 * if and wherever such third-party acknowledgments normally appear.
33 *
34 * 4. The names "Evalglare," and "Fraunhofer ISE" must
35 * not be used to endorse or promote products derived from this
36 * software without prior written permission. For written
37 * permission, please contact [email protected]
38 *
39 * 5. Products derived from this software may not be called "evalglare",
40 * nor may "evalglare" appear in their name, without prior written
41 * permission of Fraunhofer ISE.
42 *
43 * Redistribution and use in source and binary forms, WITH
44 * modification, are permitted provided that the following conditions
45 * are met:
46 *
47 *
48 * conditions 1.-5. apply
49 *
50 * 6. In order to ensure scientific correctness, any modification in source code imply fulfilling all following comditions:
51 * a) A written permission from Fraunhofer ISE. For written permission, please contact [email protected].
52 * b) The permission can be applied via email and must contain the applied changes in source code and at least two example calculations,
53 * comparing the results of the modified version and the current version of evalglare.
54 * b) Any redistribution of a modified version of evalglare must contain following note:
55 * "This software uses a modified version of the source code of evalglare."
56 *
57 *
58 *
59 *
60 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
61 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
62 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
63 * DISCLAIMED. IN NO EVENT SHALL Fraunhofer ISE OR
64 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
65 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
66 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
67 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
68 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
69 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
70 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
71 * SUCH DAMAGE.
72 * ====================================================================
73 *
74 * This product includes Radiance software
75 * (http://radsite.lbl.gov/)
76 * developed by the Lawrence Berkeley National Laboratory
77 * (http://www.lbl.gov/).
78 *
79 * The Radiance Software License, Version 1.0
80 *
81 * Copyright (c) 1990 - 2009 The Regents of the University of California,
82 * through Lawrence Berkeley National Laboratory. All rights reserved.
83 *
84 * Redistribution and use in source and binary forms, with or without
85 * modification, are permitted provided that the following conditions
86 * are met:
87 *
88 * 1. Redistributions of source code must retain the above copyright
89 * notice, this list of conditions and the following disclaimer.
90 *
91 * 2. Redistributions in binary form must reproduce the above copyright
92 * notice, this list of conditions and the following disclaimer in
93 * the documentation and/or other materials provided with the
94 * distribution.
95 *
96 * 3. The end-user documentation included with the redistribution,
97 * if any, must include the following acknowledgment:
98 * "This product includes Radiance software
99 * (http://radsite.lbl.gov/)
100 * developed by the Lawrence Berkeley National Laboratory
101 * (http://www.lbl.gov/)."
102 * Alternately, this acknowledgment may appear in the software itself,
103 * if and wherever such third-party acknowledgments normally appear.
104 *
105 * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
106 * and "The Regents of the University of California" must
107 * not be used to endorse or promote products derived from this
108 * software without prior written permission. For written
109 * permission, please contact [email protected].
110 *
111 * 5. Products derived from this software may not be called "Radiance",
112 * nor may "Radiance" appear in their name, without prior written
113 * permission of Lawrence Berkeley National Laboratory.
114 *
115 * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
116 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
117 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
118 * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
119 * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
120 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
121 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
122 * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
123 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
124 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
125 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
126 * SUCH DAMAGE.
127 */
128
129 /* evalglare.c, v0.2 2005/08/21 18:00:00 wienold */
130 /* evalglare.c, v0.3 2006/06/01 16:20:00 wienold
131 changes to the v02 version:
132 -fix problem with non-square pictures
133 -set luminance values in the hemisphere behind the line of sight to 0
134 -fix error in vcp calculation */
135 /* evalglare.c, v0.4 2006/11/13 15:10:00 wienold
136 changes to the v03 version:
137 -fix problem with tabulator in picture header "cannot read view"
138 -fix missing pixels in task area
139 evalglare.c, v0.5 2006/11/29 9:30 wienold
140 changes to the v04 version:
141 - fix problem with glare sources at the edge of the picture
142 */
143 /* evalglare.c, v0.6 2006/11/29 9:30 wienold
144 changes to the v05 version:
145 - add luminance restriction parameters in order to take into account restrictions of measuring systems
146 */
147
148 /* evalglare.c, v0.7 2007/07/11 18:00 wienold
149 changes to the v06 version:
150 - add external provision of vertical illuminance
151 */
152
153 /* evalglare.c, v0.8 2007/12/04 1:35 wienold
154 changes to the v07 version:
155 - limit dgp to max 1.0
156 - fill up cutted pictures with last known values
157 - add second detailed output
158 */
159
160 /* evalglare.c, v0.9 2008/07/02 wienold
161 changes to the v08 version:
162 - add version number in picheader
163 - add option for showing version numer
164 */
165 /* evalglare.c, v0.9a 2008/09/20 wienold
166 changes to the v09 version:
167 - reduce fix value threshold from 500 to 100
168 */
169
170 /* evalglare.c, v0.9b 2008/11/12 wienold
171 changes to the v09a version:
172 - check view type vta, if wrong than exit
173 */
174 /* evalglare.c, v0.9c 2009/03/31 wienold
175 changes to the v09b version:
176 - peak extraction is default now (-y) , for deactivation use -x
177 */
178
179 /* evalglare.c, v0.9d 2009/06/24 wienold
180 changes to the v09c version:
181 - fixed memory problem while using strcat
182 */
183 /* evalglare.c, v0.9e 2009/10/10 wienold
184 changes to the v09d version:
185 - fixed problem while reading .pic of windows version
186 */
187 /* evalglare.c, v0.9f 2009/10/21 wienold
188 changes to the v09e version:
189 - piping of input pictures possible, allow also vtv and vtc viewtyps
190 */
191
192 /* evalglare.c, v0.9g 2009/10/21 wienold
193 changes to the v09f version:
194 - modified pictool.c
195 - added -V (calc only vertical illuminance)
196 */
197 /* evalglare.c, v0.9h 2011/10/10 wienold
198 changes to the v09g version:
199 - include disability glare for age factor of 1
200 -
201 */
202 /* evalglare.c, v0.9i 2011/10/17 wienold
203 changes to the v09h version:
204 - M option: Correct wrong (measured) luminance images by feeding Ev and a luminance value, which should be replaced
205
206 */
207 /* evalglare.c, v1.0 2012/02/08 wienold
208 changes to the v09h version:
209 - include all view types
210 - check view option in header
211 - add view options in command line
212 - option for cutting out GUTH field of view
213
214 */
215
216 /* evalglare.c, v1.02 2012/03/01 wienold,reetz,grobe
217 changes to the v1.0 version:
218 - fixed buffer overflow for string variable version[40]
219 - replaced various strings to handle release by #define
220 - removed most unused variables
221 - initialized variables
222 - removed nested functions
223 - compiles now with -ansi -pedantic
224
225 */
226
227 /* evalglare.c, v1.03 2012/04/17 wienold
228 - include low light correction
229 */
230
231
232 /* evalglare.c, v1.04 2012/04/23 wienold
233 - remove bug for gen_dgp_profile output
234 */
235
236 /* evalglare.c, v1.05 2012/05/29 wienold
237 - remove variable overflow of low-light-correction for large Ev
238 */
239 /* evalglare.c, v1.06 2012/05/29 wienold
240 - initiate low-light-correction-variable
241 */
242 /* evalglare.c, v1.07 2012/05/29 wienold
243 - remove edge pixels from evaluation, when center of pixel is behind view (stability)
244 */
245
246 /* evalglare.c, v1.08 2012/09/09 wienold
247 - add direction vector in detailed output for each glare source
248 - include age correction
249 */
250 /* evalglare.c, v1.09 2012/09/09 wienold
251 - faster calculation for the use of gen_dgp_profile: no vertical illuminance calculation, only dgp is calculated, second scan + split is deactivated, when no pixel >15000 is found
252 */
253 /* evalglare.c, v1.10 2012/09/09 wienold
254 - faster calculation for the use of gen_dgp_profile: no vertical illuminance calculation, only dgp is calculated, second scan + split is deactivated, when no pixel >15000 is found
255 */
256 /* evalglare.c, v1.11 2013/01/17 wienold
257 - fix output bug of dgp, when using -i or -I
258 */
259
260 /* evalglare.c, v1.12 2013/10/31 wienold
261 - include CIE equation for disability glare, Stiles-Holladay
262 */
263 /* evalglare.c, v1.13 2014/04/06 wienold
264 - remove bug: initialize Lveil_cie_sum
265 used for the CIE equation for disability glare
266 */
267
268 /* evalglare.c, v1.14 buggy changes... removed...
269 */
270 /* evalglare.c, v1.15 add option for uniform colored glare sources
271 */
272 /* evalglare.c, v1.16 2015/05/05 remove bugs: background luminance is now calculated not cos-weighted any more, to switch on cos-weighting use -n option
273 calculation of the background luminance is now based on the real solid angle and not hard coded 2*PI any more
274 fast calculation mode: total solid angle =2*PI (before PI)
275 */
276
277 /* evalglare.c, v1.17 2015/07/15 add option for calculating band luminance -B angle_of_band
278 remove of age factor due to non proven statistical evidence
279 */
280
281 #define EVALGLARE
282
283 #define PROGNAME "evalglare"
284 #define VERSION "1.17 release 15.07.2015 by EPFL, J.Wienold"
285 #define RELEASENAME PROGNAME " " VERSION
286
287 #include "rtio.h"
288 #include "platform.h"
289 #include "pictool.h"
290 #include <math.h>
291 #include <string.h>
292 char *progname;
293
294 /* subroutine to add a pixel to a glare source */
295 void add_pixel_to_gs(pict * p, int x, int y, int gsn)
296 {
297 double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
298 new_omega, act_lum;
299
300
301 pict_get_npix(p, gsn) = pict_get_npix(p, gsn) + 1;
302 old_av_posx = pict_get_av_posx(p, gsn);
303 old_av_posy = pict_get_av_posy(p, gsn);
304 old_av_lum = pict_get_av_lum(p, gsn);
305 old_omega = pict_get_av_omega(p, gsn);
306
307 act_omega = pict_get_omega(p, x, y);
308 act_lum = luminance(pict_get_color(p, x, y));
309 new_omega = old_omega + act_omega;
310 pict_get_av_posx(p, gsn) =
311 (old_av_posx * old_omega + x * act_omega) / new_omega;
312 pict_get_av_posy(p, gsn) =
313 (old_av_posy * old_omega + y * act_omega) / new_omega;
314 if (isnan((pict_get_av_posx(p, gsn))))
315 fprintf(stderr,"error in add_pixel_to_gs %d %d %f %f %f %f\n", x, y, old_av_posy, old_omega, act_omega, new_omega);
316
317 pict_get_av_lum(p, gsn) =
318 (old_av_lum * old_omega + act_lum * act_omega) / new_omega;
319 pict_get_av_omega(p, gsn) = new_omega;
320 pict_get_gsn(p, x, y) = gsn;
321 if (act_lum < pict_get_lum_min(p, gsn)) {
322 pict_get_lum_min(p, gsn) = act_lum;
323 }
324 if (act_lum > pict_get_lum_max(p, gsn)) {
325 pict_get_lum_max(p, gsn) = act_lum;
326 }
327
328 /* printf("gsn,x,y,av_posx,av_posy,av_lum %i %f %f %f %f %f\n",gsn,x,y,pict_get_av_posx(p, gsn),pict_get_av_posy(p, gsn),pict_get_av_lum(p, gsn)); */
329
330 }
331
332 /* subroutine for peak extraction */
333 int
334 find_split(pict * p, int x, int y, double r, int i_split_start,
335 int i_split_end)
336 {
337 int i_find_split, x_min, x_max, y_min, y_max, ix, iy, iix, iiy, dx, dy,
338 out_r;
339 double r_actual;
340
341 i_find_split = 0;
342
343 x_min = 0;
344 y_min = 0;
345 x_max = pict_get_ysize(p) - 1;
346
347 y_max = pict_get_ysize(p) - 1;
348
349 for (iiy = 1; iiy <= 2; iiy++) {
350 dy = iiy * 2 - 3;
351 if (dy == -1) {
352 iy = y;
353 } else {
354 iy = y + 1;
355 }
356 while (iy <= y_max && iy >= y_min) {
357 out_r = 0;
358 for (iix = 1; iix <= 2; iix++) {
359 dx = iix * 2 - 3;
360 if (dx == -1) {
361 ix = x;
362 } else {
363 ix = x + 1;
364 }
365 while (ix <= x_max && ix >= x_min && iy >= y_min) {
366
367 r_actual =
368 acos(DOT(pict_get_cached_dir(p, x, y),
369 pict_get_cached_dir(p, ix, iy))) * 2;
370 if (r_actual <= r) {
371 out_r = 1;
372 if (pict_get_gsn(p, ix, iy) >= i_split_start
373 && pict_get_gsn(p, ix, iy) <= i_split_end) {
374 i_find_split = pict_get_gsn(p, ix, iy);
375 }
376 } else {
377 ix = -99;
378 }
379 ix = ix + dx;
380 }
381 }
382 if (out_r == 0) {
383 iy = -99;
384 }
385 iy = iy + dy;
386
387 }
388 }
389
390
391
392 return i_find_split;
393 }
394
395 /*
396 static int icomp(int* a,int* b)
397 {
398 if (*a < *b)
399 return -1;
400 if (*a > *b)
401 return 1;
402 return 0;
403 }
404
405 */
406
407 /* subroutine to find nearby glare source pixels */
408
409 int
410 find_near_pgs(pict * p, int x, int y, float r, int act_gsn, int max_gsn,
411 int min_gsn)
412 {
413 int dx, dy, i_near_gs, xx, yy, x_min, x_max, y_min, y_max, ix, iy, iix,
414 iiy, old_gsn, new_gsn, find_gsn, change, stop_y_search,
415 stop_x_search;
416 double r_actual;
417 int ixm[3];
418
419 i_near_gs = 0;
420 stop_y_search = 0;
421 stop_x_search = 0;
422 x_min = 0;
423 y_min = 0;
424 if (act_gsn == 0) {
425 x_max = x;
426 } else {
427 x_max = pict_get_xsize(p) - 1;
428 }
429
430 y_max = pict_get_ysize(p) - 1;
431
432 old_gsn = pict_get_gsn(p, x, y);
433 new_gsn = old_gsn;
434 change = 0;
435 if (act_gsn > 0) {
436 i_near_gs = pict_get_gsn(p, x, y);
437 }
438 for (iiy = 1; iiy <= 2; iiy++) {
439 dy = iiy * 2 - 3;
440 if (dy == -1) {
441 iy = y;
442 } else {
443 iy = y + 1;
444 }
445 ixm[1] = x;
446 ixm[2] = x;
447 stop_y_search = 0;
448
449
450 while (iy <= y_max && iy >= y_min) {
451 for (iix = 1; iix <= 2; iix++) {
452 dx = iix * 2 - 3;
453 ix = (ixm[1] + ixm[2]) / 2;
454 stop_x_search = 0;
455 while (ix <= x_max && ix >= x_min && stop_x_search == 0
456 && stop_y_search == 0) {
457 /* printf(" dx,act_gsn, x,y,x_max, x_min, ix ,iy , ymax,ymin %i %i %i %i %i %i %i %i %i %i\n",dx,act_gsn,x,y,x_max,x_min,ix,iy,y_max,y_min);*/
458 r_actual =
459 acos(DOT(pict_get_cached_dir(p, x, y),
460 pict_get_cached_dir(p, ix, iy))) * 2;
461 if (r_actual <= r) {
462 if (pict_get_gsn(p, ix, iy) > 0) {
463 if (act_gsn == 0) {
464 i_near_gs = pict_get_gsn(p, ix, iy);
465 stop_x_search = 1;
466 stop_y_search = 1;
467 } else {
468 find_gsn = pict_get_gsn(p, ix, iy);
469 if (pict_get_av_omega(p, old_gsn) <
470 pict_get_av_omega(p, find_gsn)
471 && pict_get_av_omega(p,
472 find_gsn) >
473 pict_get_av_omega(p, new_gsn)
474 && find_gsn >= min_gsn
475 && find_gsn <= max_gsn) {
476 /* other primary glare source found with larger solid angle -> add together all */
477 new_gsn = find_gsn;
478 change = 1;
479 stop_x_search = 1;
480 stop_y_search = 1;
481 }
482 }
483 }
484 } else {
485 ixm[iix] = ix - dx;
486 stop_x_search = 1;
487 }
488 ix = ix + dx;
489 }
490 }
491 iy = iy + dy;
492 }
493 }
494 if (change > 0) {
495 pict_get_av_lum(p, old_gsn) = 0.;
496 pict_get_av_omega(p, old_gsn) = 0.;
497 pict_get_npix(p, old_gsn) = 0.;
498 pict_get_lum_max(p, old_gsn) = 0.;
499
500
501 i_near_gs = new_gsn;
502 /* printf(" changing gs no %i",old_gsn);
503 printf(" to %i\n",new_gsn);
504 */
505 for (xx = 0; xx < pict_get_xsize(p); xx++)
506 for (yy = 0; yy < pict_get_ysize(p); yy++) {
507 if (pict_is_validpixel(p, x, y)) {
508 if (pict_get_gsn(p, xx, yy) == old_gsn) {
509 add_pixel_to_gs(p, xx, yy, new_gsn);
510 }
511 }
512 }
513 }
514
515 return i_near_gs;
516
517
518 }
519
520 /* subroutine for calculation of task luminance */
521 double get_task_lum(pict * p, int x, int y, float r, int task_color)
522 {
523 int x_min, x_max, y_min, y_max, ix, iy;
524 double r_actual, av_lum, omega_sum, act_lum;
525
526
527 x_max = pict_get_xsize(p) - 1;
528 y_max = pict_get_ysize(p) - 1;
529 x_min = 0;
530 y_min = 0;
531
532
533
534 av_lum = 0.0;
535 omega_sum = 0.0;
536
537 for (iy = y_min; iy <= y_max; iy++) {
538 for (ix = x_min; ix <= x_max; ix++) {
539
540 /* if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
541 continue;*/
542 r_actual =
543 acos(DOT(pict_get_cached_dir(p, x, y),
544 pict_get_cached_dir(p, ix, iy))) * 2;
545 act_lum = luminance(pict_get_color(p, ix, iy));
546
547 if (r_actual <= r) {
548 act_lum = luminance(pict_get_color(p, ix, iy));
549 av_lum += pict_get_omega(p, ix, iy) * act_lum;
550 omega_sum += pict_get_omega(p, ix, iy);
551 if (task_color == 1) {
552 pict_get_color(p, ix, iy)[RED] = 0.0;
553 pict_get_color(p, ix, iy)[GRN] = 0.0;
554 pict_get_color(p, ix, iy)[BLU] =
555 act_lum / WHTEFFICACY / CIE_bf;
556 }
557 }
558 }
559 }
560
561 av_lum = av_lum / omega_sum;
562 /* printf("average luminance of task %f \n",av_lum);
563 printf("solid angle of task %f \n",omega_sum);*/
564 return av_lum;
565 }
566
567 /* subroutine for calculation of band luminance */
568 double get_band_lum(pict * p, float r, int task_color)
569 {
570 int x_min, x_max, y_min, y_max, ix, iy, y_mid;
571 double r_actual, av_lum, omega_sum, act_lum;
572
573
574 x_max = pict_get_xsize(p) - 1;
575 y_max = pict_get_ysize(p) - 1;
576 x_min = 0;
577 y_min = 0;
578 y_mid = rint(y_max/2);
579
580
581
582 av_lum = 0.0;
583 omega_sum = 0.0;
584
585 for (iy = y_min; iy <= y_max; iy++) {
586 for (ix = x_min; ix <= x_max; ix++) {
587
588 /* if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
589 continue;*/
590 r_actual =
591 acos(DOT(pict_get_cached_dir(p, ix, y_mid),
592 pict_get_cached_dir(p, ix, iy))) ;
593 act_lum = luminance(pict_get_color(p, ix, iy));
594
595 if ((r_actual <= r) || (iy == y_mid) ) {
596 act_lum = luminance(pict_get_color(p, ix, iy));
597 av_lum += pict_get_omega(p, ix, iy) * act_lum;
598 omega_sum += pict_get_omega(p, ix, iy);
599 if (task_color == 1) {
600 pict_get_color(p, ix, iy)[RED] = 0.0;
601 pict_get_color(p, ix, iy)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
602 pict_get_color(p, ix, iy)[BLU] = 0.0;
603 }
604 }
605 }
606 }
607
608 av_lum = av_lum / omega_sum;
609 /* printf("average luminance of band %f \n",av_lum);*/
610 /* printf("solid angle of band %f \n",omega_sum);*/
611 return av_lum;
612 }
613
614
615
616
617
618 /* subroutine for coloring the glare sources
619 red is used only for explicit call of the subroutine with acol=0 , e.g. for disability glare
620 -1: set to grey*/
621 int setglcolor(pict * p, int x, int y, int acol, int uniform_gs, double u_r, double u_g ,double u_b)
622 {
623 int icol;
624 double pcol[3][1000], act_lum, l;
625
626 l=u_r+u_g+u_b ;
627
628 pcol[RED][0] = 1.0 / CIE_rf;
629 pcol[GRN][0] = 0.;
630 pcol[BLU][0] = 0.;
631
632 pcol[RED][1] = 0.;
633 pcol[GRN][1] = 1.0 / CIE_gf;
634 pcol[BLU][1] = 0.;
635
636 pcol[RED][2] = 0.;
637 pcol[GRN][2] = 0.;
638 pcol[BLU][2] = 1 / CIE_bf;
639
640 pcol[RED][3] = 1.0 / (1.0 - CIE_bf);
641 pcol[GRN][3] = 1.0 / (1.0 - CIE_bf);
642 pcol[BLU][3] = 0.;
643
644 pcol[RED][4] = 1.0 / (1.0 - CIE_gf);
645 pcol[GRN][4] = 0.;
646 pcol[BLU][4] = 1.0 / (1.0 - CIE_gf);
647
648 pcol[RED][5] = 0.;
649 pcol[GRN][5] = 1.0 / (1.0 - CIE_rf);
650 pcol[BLU][5] = 1.0 / (1.0 - CIE_rf);
651
652 pcol[RED][6] = 0.;
653 pcol[GRN][6] = 1.0 / (1.0 - CIE_rf);
654 pcol[BLU][6] = 1.0 / (1.0 - CIE_rf);
655
656
657 pcol[RED][999] = 1.0 / WHTEFFICACY;
658 pcol[GRN][999] = 1.0 / WHTEFFICACY;
659 pcol[BLU][999] = 1.0 / WHTEFFICACY;
660
661
662 pcol[RED][998] = u_r /(l* CIE_rf) ;
663 pcol[GRN][998] = u_g /(l* CIE_gf);
664 pcol[BLU][998] = u_b /(l* CIE_bf);
665 /* printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l);*/
666 icol = acol ;
667 if ( acol == -1) {icol=999;
668 }else{if (acol>0){icol = acol % 5 +1;
669 }};
670 if ( uniform_gs == 1) {icol=998;
671 } ;
672
673 act_lum = luminance(pict_get_color(p, x, y));
674
675 pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY;
676 pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY;
677 pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY;
678
679 return icol;
680 }
681
682 /* this is the smooting subroutine */
683 void add_secondary_gs(pict * p, int x, int y, float r, int act_gs, int uniform_gs, double u_r, double u_g ,double u_b)
684 {
685 int x_min, x_max, y_min, y_max, ix, iy, icol;
686 double r_actual, omega_gs, omega_total, om;
687
688
689
690 omega_gs = 0.0;
691 omega_total = 0.0;
692 x_min = x - r;
693 if (x_min < 0) {
694 x_min = 0;
695 }
696 x_max = x + r;
697 if (x_max > pict_get_xsize(p) - 1) {
698 x_max = pict_get_xsize(p) - 2;
699 }
700 y_min = y - r;
701 if (y_min < 0) {
702 y_min = 0;
703 }
704 y_max = y + r;
705 if (y_max > pict_get_ysize(p) - 1) {
706 y_max = pict_get_ysize(p) - 2;
707 }
708
709 for (ix = x_min; ix <= x_max; ix++)
710 for (iy = y_min; iy <= y_max; iy++) {
711 r_actual = sqrt((x - ix) * (x - ix) + (y - iy) * (y - iy));
712 if (r_actual <= r) {
713 om = pict_get_omega(p, ix, iy);
714 omega_total += om;
715 if (pict_get_gsn(p, ix, iy) == act_gs
716 && pict_get_pgs(p, ix, iy) == 1) {
717 omega_gs = omega_gs + 1 * om;
718 }
719 }
720 }
721
722
723 if (omega_gs / omega_total > 0.2) {
724 /* add pixel to gs */
725
726 add_pixel_to_gs(p, x, y, act_gs);
727
728 /* color pixel of gs */
729
730 icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b);
731
732 }
733 }
734
735 /* subroutine for removing a pixel from a glare source */
736 void split_pixel_from_gs(pict * p, int x, int y, int new_gsn, int uniform_gs, double u_r, double u_g , double u_b)
737 {
738 int old_gsn, icol;
739 double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
740 new_omega, act_lum;
741
742
743 /* change existing gs */
744 old_gsn = pict_get_gsn(p, x, y);
745 pict_get_npix(p, old_gsn) = pict_get_npix(p, old_gsn) - 1;
746
747 act_omega = pict_get_omega(p, x, y);
748 old_av_posx = pict_get_av_posx(p, old_gsn);
749 old_av_posy = pict_get_av_posy(p, old_gsn);
750 old_omega = pict_get_av_omega(p, old_gsn);
751
752 new_omega = old_omega - act_omega;
753 pict_get_av_omega(p, old_gsn) = new_omega;
754 pict_get_av_posx(p, old_gsn) =
755 (old_av_posx * old_omega - x * act_omega) / new_omega;
756 pict_get_av_posy(p, old_gsn) =
757 (old_av_posy * old_omega - y * act_omega) / new_omega;
758
759
760 act_lum = luminance(pict_get_color(p, x, y));
761 old_av_lum = pict_get_av_lum(p, old_gsn);
762 pict_get_av_lum(p, old_gsn) =
763 (old_av_lum * old_omega - act_lum * act_omega) / new_omega;
764 /* add pixel to new gs */
765
766 add_pixel_to_gs(p, x, y, new_gsn);
767
768 /* color pixel of new gs */
769
770 icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b);
771
772
773 }
774
775 /* subroutine for the calculation of the position index */
776 float get_posindex(pict * p, float x, float y, int postype)
777 {
778 float posindex;
779 double teta, phi, sigma, tau, deg, d, s, r, fact;
780
781
782 pict_get_vangle(p, x, y, p->view.vdir, p->view.vup, &phi);
783 pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &teta);
784 pict_get_sigma(p, x, y, p->view.vdir, p->view.vup, &sigma);
785 pict_get_tau(p, x, y, p->view.vdir, p->view.vup, &tau);
786
787 /* return (phi+teta+2*3.1415927); */
788
789 deg = 180 / 3.1415927;
790 fact = 0.8;
791 if (phi == 0) {
792 phi = 0.00001;
793 }
794 if (sigma <= 0) {
795 sigma = -sigma;
796 }
797 if (teta == 0) {
798 teta = 0.0001;
799 }
800 tau = tau * deg;
801 sigma = sigma * deg;
802
803 posindex =
804 exp((35.2 - 0.31889 * tau -
805 1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 +
806 0.26667 * tau -
807 0.002963 * tau *
808 tau) / 100000 *
809 sigma * sigma);
810 /* below line of sight, using Iwata model */
811 if (phi < 0) {
812 d = 1 / tan(phi);
813 s = tan(teta) / tan(phi);
814 r = sqrt(1 / d * 1 / d + s * s / d / d);
815 if (r > 0.6)
816 fact = 1.2;
817 if (r > 3) {
818 fact = 1.2;
819 r = 3;
820 }
821
822 posindex = 1 + fact * r;
823 }
824 if (posindex > 16)
825 posindex = 16;
826
827 return posindex;
828
829 }
830
831
832
833 double get_upper_cut_2eyes(float teta)
834 {
835 double phi;
836
837 phi=pow(7.7458218+0.00057407915*teta-0.00021746318*teta*teta+8.5572726e-6*teta*teta*teta,2);
838
839 return phi;
840 }
841 double get_lower_cut_2eyes(float teta)
842 {
843 double phi;
844
845 phi=1/(-0.014699242-1.5541106e-5*teta+4.6898068e-6*teta*teta-5.1539687e-8*teta*teta*teta);
846
847 return phi;
848 }
849
850 double get_lower_cut_central(float teta)
851 {
852 double phi;
853
854 phi=(68.227109-2.9524084*teta+0.046674262*teta*teta)/(1-0.042317294*teta+0.00075698419*teta*teta-6.5364285e-7*teta*teta*teta);
855
856 if (teta>73) {
857 phi=60;
858 }
859
860 return phi;
861 }
862
863 /* subroutine for the cutting the total field of view */
864 void cut_view_1(pict*p)
865 {
866 int x,y;
867 double border,ang,teta,phi,phi2;
868 for(x=0;x<pict_get_xsize(p)-1;x++)
869 for(y=0;y<pict_get_ysize(p)-1;y++) {
870 if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
871 if (pict_is_validpixel(p, x, y)) {
872 pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
873 pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
874 pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
875 phi=phi*180/3.1415927;
876 phi2=phi2*180/3.1415927;
877 teta=teta*180/3.1415927;
878 if (teta<0) {
879 teta=-teta;
880 }
881 if(phi2>0){
882 border=get_upper_cut_2eyes(teta);
883
884 if (phi>border) {
885 pict_get_color(p,x,y)[RED]=0;
886 pict_get_color(p,x,y)[GRN]=0;
887 pict_get_color(p,x,y)[BLU]=0;
888 }
889 }else{
890 border=get_lower_cut_2eyes(180-teta);
891 if (-phi<border && teta>135) {
892 pict_get_color(p,x,y)[RED]=0;
893 pict_get_color(p,x,y)[GRN]=0;
894 pict_get_color(p,x,y)[BLU]=0;
895 }
896 }
897
898
899
900 }else{
901 pict_get_color(p,x,y)[RED]=0;
902 pict_get_color(p,x,y)[GRN]=0;
903 pict_get_color(p,x,y)[BLU]=0;
904
905
906 }
907
908 /* printf("teta,phi,border=%f %f %f\n",teta,phi,border);*/
909 }
910 }
911
912
913 return;
914
915 }
916 /* subroutine for the cutting the field of view seen by both eyes*/
917 void cut_view_2(pict*p)
918 {
919
920 int x,y;
921 double border,ang,teta,phi,phi2;
922 for(x=0;x<pict_get_xsize(p)-1;x++)
923 for(y=0;y<pict_get_ysize(p)-1;y++) {
924 if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
925 if (pict_is_validpixel(p, x, y)) {
926 pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
927 pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
928 pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
929 pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
930 phi=phi*180/3.1415927;
931 phi2=phi2*180/3.1415927;
932 teta=teta*180/3.1415927;
933 /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
934 if (teta<0) {
935 teta=-teta;
936 }
937 if(phi2>0){
938 border=60;
939
940 if (phi>border) {
941 pict_get_color(p,x,y)[RED]=0;
942 pict_get_color(p,x,y)[GRN]=0;
943 pict_get_color(p,x,y)[BLU]=0;
944 }
945 }else{
946 border=get_lower_cut_central(180-teta);
947 if (phi>border) {
948 pict_get_color(p,x,y)[RED]=0;
949 pict_get_color(p,x,y)[GRN]=0;
950 pict_get_color(p,x,y)[BLU]=0;
951 }
952 }
953
954
955
956 }else{
957 pict_get_color(p,x,y)[RED]=0;
958 pict_get_color(p,x,y)[GRN]=0;
959 pict_get_color(p,x,y)[BLU]=0;
960
961
962 }
963 }
964 }
965
966
967 return;
968
969 }
970
971 /* subroutine for the cutting the field of view seen by both eyes
972 luminance is modified by position index - just for experiments - undocumented */
973 void cut_view_3(pict*p)
974 {
975
976 int x,y;
977 double border,ang,teta,phi,phi2,lum,newlum;
978 for(x=0;x<pict_get_xsize(p)-1;x++)
979 for(y=0;y<pict_get_ysize(p)-1;y++) {
980 if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
981 if (DOT(pict_get_cached_dir(p,x,y),p->view.vdir) >= 0.0) {
982 pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
983 pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
984 pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
985 pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
986 phi=phi*180/3.1415927;
987 phi2=phi2*180/3.1415927;
988 teta=teta*180/3.1415927;
989 lum=luminance(pict_get_color(p,x,y));
990
991 newlum=lum/get_posindex(p,x,y,0);
992
993 pict_get_color(p,x,y)[RED]=newlum/WHTEFFICACY;
994 pict_get_color(p,x,y)[GRN]=newlum/WHTEFFICACY;
995 pict_get_color(p,x,y)[BLU]=newlum/WHTEFFICACY;
996
997
998
999
1000 /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
1001 if (teta<0) {
1002 teta=-teta;
1003 }
1004 if(phi2>0){
1005 border=60;
1006
1007 if (phi>border) {
1008 pict_get_color(p,x,y)[RED]=0;
1009 pict_get_color(p,x,y)[GRN]=0;
1010 pict_get_color(p,x,y)[BLU]=0;
1011 }
1012 }else{
1013 border=get_lower_cut_central(180-teta);
1014 if (phi>border) {
1015 pict_get_color(p,x,y)[RED]=0;
1016 pict_get_color(p,x,y)[GRN]=0;
1017 pict_get_color(p,x,y)[BLU]=0;
1018 }
1019 }
1020
1021
1022
1023 }else{
1024 pict_get_color(p,x,y)[RED]=0;
1025 pict_get_color(p,x,y)[GRN]=0;
1026 pict_get_color(p,x,y)[BLU]=0;
1027
1028
1029 }
1030 }
1031 }
1032
1033
1034 return;
1035
1036 }
1037
1038 /* subroutine for the calculation of the daylight glare index */
1039
1040
1041 float get_dgi(pict * p, float lum_backg, int igs, int posindex_2)
1042 {
1043 float dgi, sum_glare, omega_s;
1044 int i;
1045
1046
1047 sum_glare = 0;
1048 omega_s = 0;
1049
1050 for (i = 0; i <= igs; i++) {
1051
1052 if (pict_get_npix(p, i) > 0) {
1053 omega_s =
1054 pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1055 get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1056 sum_glare += 0.478 * pow(pict_get_av_lum(p, i), 1.6) * pow(omega_s,0.8) / (lum_backg + 0.07 * pow(pict_get_av_omega(p, i),0.5) * pict_get_av_lum(p, i));
1057 /* printf("i,sum_glare %i %f\n",i,sum_glare); */
1058 }
1059 }
1060 dgi = 10 * log10(sum_glare);
1061
1062 return dgi;
1063
1064 }
1065
1066 /* subroutine for the calculation of the daylight glare probability */
1067 double
1068 get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3,
1069 double a4, double a5, double c1, double c2, double c3,
1070 int posindex_2)
1071 {
1072 double dgp;
1073 double sum_glare;
1074 int i;
1075
1076 sum_glare = 0;
1077 if (igs > 0) {
1078 for (i = 0; i <= igs; i++) {
1079
1080 if (pict_get_npix(p, i) > 0) {
1081 sum_glare +=
1082 pow(pict_get_av_lum(p, i),
1083 a1) / pow(get_posindex(p, pict_get_av_posx(p, i),
1084 pict_get_av_posy(p, i),
1085 posindex_2),
1086 a4) * pow(pict_get_av_omega(p, i), a2);
1087 /* printf("i,sum_glare %i %f\n",i,sum_glare);*/
1088 }
1089 }
1090 dgp =
1091 c1 * pow(E_v,
1092 a5) + c3 + c2 * log10(1 + sum_glare / pow(E_v, a3));
1093 } else {
1094 dgp = c3 + c1 * pow(E_v, a5);
1095 }
1096
1097 if (dgp > 1) {
1098 dgp = 1;
1099 }
1100 /* printf("dgp %f\n",dgp);*/
1101 return dgp;
1102
1103 }
1104
1105 /* subroutine for the calculation of the visual comfort probability */
1106 float get_vcp(pict * p, double lum_a, int igs, int posindex_2)
1107 {
1108 float vcp;
1109 double sum_glare, dgr;
1110 int i, i_glare;
1111
1112
1113 sum_glare = 0;
1114 i_glare = 0;
1115 for (i = 0; i <= igs; i++) {
1116
1117 if (pict_get_npix(p, i) > 0) {
1118 i_glare = i_glare + 1;
1119 sum_glare +=
1120 (0.5 * pict_get_av_lum(p, i) *
1121 (20.4 * pict_get_av_omega(p, i) +
1122 1.52 * pow(pict_get_av_omega(p, i),
1123 0.2) - 0.075)) / (get_posindex(p,
1124 pict_get_av_posx
1125 (p, i),
1126 pict_get_av_posy
1127 (p, i),
1128 posindex_2) *
1129 pow(lum_a, 0.44));
1130
1131 }
1132 }
1133 dgr = pow(sum_glare, pow(i_glare, -0.0914));
1134
1135 vcp = 50 * erf((6.347 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1136 if (dgr > 750) {
1137 vcp = 0;
1138 }
1139 if (dgr < 20) {
1140 vcp = 100;
1141 }
1142
1143
1144 return vcp;
1145
1146 }
1147
1148 /* subroutine for the calculation of the unified glare rating */
1149 float get_ugr(pict * p, double lum_backg, int igs, int posindex_2)
1150 {
1151 float ugr;
1152 double sum_glare;
1153 int i, i_glare;
1154
1155
1156 sum_glare = 0;
1157 i_glare = 0;
1158 for (i = 0; i <= igs; i++) {
1159
1160 if (pict_get_npix(p, i) > 0) {
1161 i_glare = i_glare + 1;
1162 sum_glare +=
1163 pow(pict_get_av_lum(p, i) /
1164 get_posindex(p, pict_get_av_posx(p, i),
1165 pict_get_av_posy(p, i), posindex_2),
1166 2) * pict_get_av_omega(p, i);
1167
1168 }
1169 }
1170 ugr = 8 * log10(0.25 / lum_backg * sum_glare);
1171
1172 return ugr;
1173
1174 }
1175
1176
1177 /* subroutine for the calculation of the disability glare according to Poynter */
1178 float get_disability(pict * p, double lum_backg, int igs)
1179 {
1180 float disab;
1181 double sum_glare, sigma, deg;
1182 int i, i_glare;
1183
1184
1185 sum_glare = 0;
1186 i_glare = 0;
1187 deg = 180 / 3.1415927;
1188 for (i = 0; i <= igs; i++) {
1189
1190 if (pict_get_npix(p, i) > 0) {
1191 i_glare = i_glare + 1;
1192 pict_get_sigma(p, pict_get_av_posx(p, i),
1193 pict_get_av_posy(p, i), p->view.vdir,
1194 p->view.vup, &sigma);
1195
1196 sum_glare +=
1197 pict_get_av_lum(p,
1198 i) * cos(sigma +
1199 0.00000000001) *
1200 pict_get_av_omega(p, i)
1201 / (deg * sigma + 0.00000000001);
1202 if (isnan(sum_glare)) {
1203 printf("sigma for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1204 printf("omega for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1205 printf("avlum for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1206 printf("avlum for %f %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i), sigma);
1207 }
1208
1209 }
1210 }
1211 disab = 5 * sum_glare;
1212
1213 return disab;
1214
1215 }
1216
1217
1218
1219
1220
1221
1222 /* subroutine for the calculation of the cie glare index */
1223
1224 float
1225 get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1226 {
1227 float cgi;
1228 double sum_glare;
1229 int i, i_glare;
1230
1231
1232 sum_glare = 0;
1233 i_glare = 0;
1234 for (i = 0; i <= igs; i++) {
1235
1236 if (pict_get_npix(p, i) > 0) {
1237 i_glare = i_glare + 1;
1238 sum_glare +=
1239 pow(pict_get_av_lum(p, i) /
1240 get_posindex(p, pict_get_av_posx(p, i),
1241 pict_get_av_posy(p, i), posindex_2),
1242 2) * pict_get_av_omega(p, i);
1243
1244 }
1245 }
1246 cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare);
1247
1248 return cgi;
1249
1250 }
1251
1252
1253 #ifdef EVALGLARE
1254
1255
1256 /* main program
1257 ------------------------------------------------------------------------------------------------------------------*/
1258
1259
1260 int main(int argc, char **argv)
1261 {
1262 #define CLINEMAX 4095 /* memory allocated for command line string */
1263 pict *p = pict_create();
1264 int skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, y1, fill, yfillmax, yfillmin,
1265 ext_vill, set_lum_max, set_lum_max2, task_color, i_splitstart,
1266 i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,
1267 igs, actual_igs, icol, xt, yt, change, before_igs, sgs, splithigh,uniform_gs,
1268 detail_out, posindex_picture, non_cos_lb, rx, ry, rmx, rmy,apply_disability,band_calc,band_color;
1269 double lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma,
1270 E_vl_ext, lum_max, new_lum_max, r_center,
1271 search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,
1272 omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang,
1273 l_max, lum_backg, lum_backg_cos, omega_sources, lum_sources,
1274 lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum;
1275 float lum_task, lum_thres, dgi, vcp, cgi, ugr, limit,
1276 abs_max, Lveil;
1277 char file_out[500], file_out2[500], version[500];
1278 char *cline;
1279 VIEW userview = STDVIEW;
1280 int gotuserview = 0;
1281
1282 /*set required user view parameters to invalid values*/
1283 uniform_gs = 0;
1284 apply_disability = 0;
1285 disability_thresh = 0;
1286 Lveil_cie_sum=0.0;
1287 skip_second_scan=0;
1288 lum_total_max=0.0;
1289 calcfast=0;
1290 age_corr_factor = 1.0;
1291 age_corr = 0;
1292 age = 20.0;
1293 userview.horiz = 0;
1294 userview.vert = 0;
1295 userview.type = 0;
1296 dgp_ext = 0;
1297 E_vl_ext = 0.0;
1298 new_lum_max = 0.0;
1299 lum_max = 0.0;
1300 omegat = 0.0;
1301 yt = 0;
1302 xt = 0;
1303 yfillmin = 0;
1304 yfillmax = 0;
1305 cut_view = 0;
1306 cut_view_type = 0;
1307 setvalue = 2e09;
1308 omega_cos_contr = 0.0;
1309 lum_ideal = 0.0;
1310 max_angle = 0.2;
1311 lum_thres = 5.0;
1312 task_lum = 0;
1313 sgs = 0;
1314 splithigh = 1;
1315 detail_out = 0;
1316 detail_out2 = 0;
1317 posindex_picture = 0;
1318 checkfile = 0;
1319 ext_vill = 0;
1320 fill = 0;
1321 a1 = 2.0;
1322 a2 = 1.0;
1323 a3 = 1.87;
1324 a4 = 2.0;
1325 a5 = 1.0;
1326 c1 = 5.87e-05;
1327 c2 = 0.092;
1328 c3 = 0.159;
1329 non_cos_lb = 1;
1330 posindex_2 = 0;
1331 task_color = 0;
1332 limit = 50000.0;
1333 set_lum_max = 0;
1334 set_lum_max2 = 0;
1335 abs_max = 0;
1336 progname = argv[0];
1337 E_v_contr = 0.0;
1338 strcpy(version, "1.15 release 05.02.2015 by J.Wienold");
1339 low_light_corr=1.0;
1340 output = 0;
1341 calc_vill = 0;
1342 band_avlum = -99;
1343 band_calc = 0;
1344 /* command line for output picture*/
1345
1346 cline = (char *) malloc(CLINEMAX+1);
1347 cline[0] = '\0';
1348 for (i = 0; i < argc; i++) {
1349 /* fprintf(stderr, "%d %d \n",i,strlen(argv[i]));*/
1350 if (strlen(cline)+strlen(argv[i])+strlen(RELEASENAME)+2 >=CLINEMAX) {
1351 exit (-1);
1352 }
1353 else {
1354 strcat(cline, argv[i]);
1355 strcat(cline, " ");
1356 }
1357 }
1358 strcat(cline, "\n");
1359 strcat(cline, RELEASENAME);
1360 strcat(cline, "\n");
1361
1362
1363 /* program options */
1364
1365 for (i = 1; i < argc && argv[i][0] == '-'; i++) {
1366 /* expand arguments */
1367 while ((rval = expandarg(&argc, &argv, i)) > 0);
1368 if (rval < 0) {
1369 fprintf(stderr, "%s: cannot expand '%s'", argv[0], argv[i]);
1370 exit(1);
1371 }
1372 rval = getviewopt(&userview, argc - i, argv + i);
1373 if (rval >= 0) {
1374 i += rval;
1375 gotuserview++;
1376 continue;
1377 }
1378 switch (argv[i][1]) {
1379 case 'a':
1380 age = atof(argv[++i]);
1381 age_corr = 1;
1382 printf("age factor not supported any more \n");
1383 exit(1);
1384 break;
1385 case 'b':
1386 lum_thres = atof(argv[++i]);
1387 break;
1388 case 'c':
1389 checkfile = 1;
1390 strcpy(file_out, argv[++i]);
1391 break;
1392 case 'u':
1393 uniform_gs = 1;
1394 u_r = atof(argv[++i]);
1395 u_g = atof(argv[++i]);
1396 u_b = atof(argv[++i]);
1397 break;
1398 case 'r':
1399 max_angle = atof(argv[++i]);
1400 break;
1401 case 's':
1402 sgs = 1;
1403 break;
1404 case 'k':
1405 apply_disability = 1;
1406 disability_thresh = atof(argv[++i]);
1407 break;
1408
1409 case 'p':
1410 posindex_picture = 1;
1411 break;
1412
1413 case 'y':
1414 splithigh = 1;
1415 break;
1416 case 'x':
1417 splithigh = 0;
1418 break;
1419 case 'Y':
1420 splithigh = 1;
1421 limit = atof(argv[++i]);
1422 break;
1423
1424 case 'i':
1425 ext_vill = 1;
1426 E_vl_ext = atof(argv[++i]);
1427 break;
1428 case 'I':
1429 ext_vill = 1;
1430 fill = 1;
1431 E_vl_ext = atof(argv[++i]);
1432 yfillmax = atoi(argv[++i]);
1433 yfillmin = atoi(argv[++i]);
1434 break;
1435 case 'd':
1436 detail_out = 1;
1437 break;
1438 case 'D':
1439 detail_out2 = 1;
1440 break;
1441 case 'm':
1442 set_lum_max = 1;
1443 lum_max = atof(argv[++i]);
1444 new_lum_max = atof(argv[++i]);
1445 strcpy(file_out2, argv[++i]);
1446 /* printf("max lum set to %f\n",new_lum_max);*/
1447 break;
1448
1449 case 'M':
1450 set_lum_max2 = 1;
1451 lum_max = atof(argv[++i]);
1452 E_vl_ext = atof(argv[++i]);
1453 strcpy(file_out2, argv[++i]);
1454 /* printf("max lum set to %f\n",new_lum_max);*/
1455 break;
1456 case 'n':
1457 non_cos_lb = 0;
1458 break;
1459
1460 case 't':
1461 task_lum = 1;
1462 xt = atoi(argv[++i]);
1463 yt = atoi(argv[++i]);
1464 omegat = atof(argv[++i]);
1465 task_color = 0;
1466 break;
1467 case 'T':
1468 task_lum = 1;
1469 xt = atoi(argv[++i]);
1470 yt = atoi(argv[++i]);
1471 /* omegat= DEG2RAD(atof(argv[++i]));*/
1472 omegat = atof(argv[++i]);
1473 task_color = 1;
1474 break;
1475 case 'B':
1476 band_calc = 1;
1477 band_angle = atof(argv[++i]);
1478 band_color = 1;
1479 break;
1480
1481
1482 case 'w':
1483 a1 = atof(argv[++i]);
1484 a2 = atof(argv[++i]);
1485 a3 = atof(argv[++i]);
1486 a4 = atof(argv[++i]);
1487 a5 = atof(argv[++i]);
1488 c1 = atof(argv[++i]);
1489 c2 = atof(argv[++i]);
1490 c3 = atof(argv[++i]);
1491 break;
1492 case 'V':
1493 calc_vill = 1;
1494 break;
1495 case 'G':
1496 cut_view = 1;
1497 cut_view_type= atof(argv[++i]);
1498 break;
1499 case 'g':
1500 cut_view = 2;
1501 cut_view_type= atof(argv[++i]);
1502 break;
1503
1504 /*case 'v':
1505 printf("evalglare %s \n",version);
1506 exit(1); */
1507
1508 case '1':
1509 output = 1;
1510 break;
1511
1512 case 'v':
1513 if (argv[i][2] == '\0') {
1514 printf("%s \n",RELEASENAME);
1515 exit(1);
1516 }
1517 if (argv[i][2] != 'f')
1518 goto userr;
1519 rval = viewfile(argv[++i], &userview, NULL);
1520 if (rval < 0) {
1521 fprintf(stderr,
1522 "%s: cannot open view file \"%s\"\n",
1523 progname, argv[i]);
1524 exit(1);
1525 } else if (rval == 0) {
1526 fprintf(stderr,
1527 "%s: bad view file \"%s\"\n", progname, argv[i]);
1528 exit(1);
1529 } else {
1530 gotuserview++;
1531 }
1532 break;
1533
1534
1535 default:
1536 goto userr;
1537 }
1538 }
1539
1540 /*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/
1541
1542 if (output == 1 && ext_vill == 1) {
1543 calcfast=1;
1544 }
1545
1546 /* read picture file */
1547 if (i == argc) {
1548 SET_FILE_BINARY(stdin);
1549 FILE *fp = fdopen(fileno(stdin), "rb");
1550 if (!(fp)) {
1551 fprintf(stderr, "fdopen on stdin failed\n");
1552 return EXIT_FAILURE;
1553 }
1554 if (!(pict_read_fp(p, fp)))
1555 return EXIT_FAILURE;;
1556 } else {
1557 if (!pict_read(p, argv[i]))
1558 return EXIT_FAILURE;
1559 }
1560 if (gotuserview) {
1561 if (p->valid_view)
1562 fprintf(stderr,
1563 "warning: overriding image view by commandline argument\n");
1564 if ((userview.horiz == 0) || (userview.vert == 0) || (userview.type == 0)) {
1565 fprintf(stderr, "error: if specified, a view must at least contain -vt -vv and -vh\n");
1566 return EXIT_FAILURE;
1567 }
1568 p->view = userview;
1569 if (!(pict_update_view(p))) {
1570 fprintf(stderr, "error: invalid view specified");
1571 return EXIT_FAILURE;
1572 }
1573 pict_update_evalglare_caches(p);
1574 } else if (!(p->valid_view)) {
1575 fprintf(stderr, "error: no valid view specified\n");
1576 return EXIT_FAILURE;
1577 }
1578
1579 /* fprintf(stderr,"Picture read!");*/
1580
1581 /* command line for output picture*/
1582
1583 pict_set_comment(p, cline);
1584
1585 if (p->view.type == VT_PAR) {
1586 fprintf(stderr, "wrong view type! must not be parallel ! \n");
1587 exit(1);
1588 }
1589
1590
1591 /* Check size of search radius */
1592 rmx = (pict_get_xsize(p) / 2);
1593 rmy = (pict_get_ysize(p) / 2);
1594 rx = (pict_get_xsize(p) / 2 + 10);
1595 ry = (pict_get_ysize(p) / 2 + 10);
1596 r_center =
1597 acos(DOT(pict_get_cached_dir(p, rmx, rmy),
1598 pict_get_cached_dir(p, rx, ry))) * 2 / 10;
1599 search_pix = max_angle / r_center;
1600 if (search_pix < 1.0) {
1601 fprintf(stderr,
1602 "warning: search radius less than 1 pixel! deactivating smoothing and peak extraction...\n");
1603 splithigh = 0;
1604 sgs = 0;
1605
1606 } else {
1607 if (search_pix < 3.0) {
1608 fprintf(stderr,
1609 "warning: search radius less than 3 pixels! -> %f \n",
1610 search_pix);
1611
1612 }
1613 }
1614
1615 /* Check task position */
1616
1617 if (task_lum == 1) {
1618 if (xt >= pict_get_xsize(p) || yt >= pict_get_ysize(p) || xt < 0
1619 || yt < 0) {
1620 fprintf(stderr,
1621 "error: task position outside picture!! exit...");
1622 exit(1);
1623 }
1624 }
1625
1626
1627 /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1628
1629 sang = 0.0;
1630 E_v = 0.0;
1631 E_v_dir = 0.0;
1632 avlum = 0.0;
1633 pict_new_gli(p);
1634 igs = 0;
1635
1636 /* cut out GUTH field of view and exit without glare evaluation */
1637 if (cut_view==2) {
1638 if (cut_view_type==1) {
1639 cut_view_1(p);
1640 }else{
1641
1642 if (cut_view_type==2) {
1643 cut_view_2(p);
1644 }else{
1645 if (cut_view_type==3) {
1646 fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
1647 cut_view_3(p);
1648 }else{
1649 fprintf(stderr,"error: no valid option for view cutting!!");
1650 exit(1);
1651 }
1652 }}
1653 pict_write(p, file_out);
1654 exit(1); }
1655
1656
1657
1658
1659
1660
1661 /* write positionindex into checkfile and exit, activated by option -p */
1662 if (posindex_picture == 1) {
1663 for (x = 0; x < pict_get_xsize(p); x++)
1664 for (y = 0; y < pict_get_ysize(p); y++) {
1665 if (pict_get_hangle
1666 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1667 if (pict_is_validpixel(p, x, y)) {
1668 lum =
1669 get_posindex(p, x, y,
1670 posindex_2) / WHTEFFICACY;
1671
1672 pict_get_color(p, x, y)[RED] = lum;
1673 pict_get_color(p, x, y)[GRN] = lum;
1674 pict_get_color(p, x, y)[BLU] = lum;
1675 lum_task = luminance(pict_get_color(p, x, y));
1676 /*printf("x,y,posindex=%i %i %f %f\n",x,y,lum*WHTEFFICACY,lum_task);*/
1677 }
1678 }
1679 }
1680 pict_write(p, file_out);
1681 exit(1);
1682 }
1683
1684
1685 /* calculation of vertical illuminance, average luminance, in case of filling activated-> fill picture */
1686 /* fill, if necessary from 0 to yfillmin */
1687
1688 if (fill == 1) {
1689 for (x = 0; x < pict_get_xsize(p); x++)
1690 for (y = yfillmin; y > 0; y = y - 1) {
1691 y1 = y + 1;
1692 lum = luminance(pict_get_color(p, x, y1));
1693 pict_get_color(p, x, y)[RED] = lum / 179.0;
1694 pict_get_color(p, x, y)[GRN] = lum / 179.0;
1695 pict_get_color(p, x, y)[BLU] = lum / 179.0;
1696 }
1697 }
1698
1699 if (calcfast == 0) {
1700 for (x = 0; x < pict_get_xsize(p); x++)
1701 for (y = 0; y < pict_get_ysize(p); y++) {
1702 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
1703 if (pict_is_validpixel(p, x, y)) {
1704 lum = luminance(pict_get_color(p, x, y));
1705 if (fill == 1 && y >= yfillmax) {
1706 y1 = y - 1;
1707 lum = luminance(pict_get_color(p, x, y1));
1708 pict_get_color(p, x, y)[RED] = lum / 179.0;
1709 pict_get_color(p, x, y)[GRN] = lum / 179.0;
1710 pict_get_color(p, x, y)[BLU] = lum / 179.0;
1711 }
1712
1713 if (lum > abs_max) {
1714 abs_max = lum;
1715 }
1716 /* set luminance restriction, if -m is set */
1717 if (set_lum_max == 1 && lum >= lum_max) {
1718 pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
1719 pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
1720 pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
1721 /* printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
1722 lum = luminance(pict_get_color(p, x, y));
1723
1724 }
1725 if (set_lum_max2 == 1 && lum >= lum_max) {
1726
1727 E_v_contr +=
1728 DOT(p->view.vdir,
1729 pict_get_cached_dir(p, x,
1730 y)) * pict_get_omega(p,
1731 x,
1732 y)
1733 * lum;
1734 omega_cos_contr +=
1735 DOT(p->view.vdir,
1736 pict_get_cached_dir(p, x,
1737 y)) * pict_get_omega(p,
1738 x,
1739 y)
1740 * 1;
1741 }
1742
1743
1744 sang += pict_get_omega(p, x, y);
1745 E_v +=
1746 DOT(p->view.vdir,
1747 pict_get_cached_dir(p, x,
1748 y)) * pict_get_omega(p, x,
1749 y) *
1750 lum;
1751 avlum +=
1752 pict_get_omega(p, x,
1753 y) * luminance(pict_get_color(p, x,
1754 y));
1755 } else {
1756 pict_get_color(p, x, y)[RED] = 0.0;
1757 pict_get_color(p, x, y)[GRN] = 0.0;
1758 pict_get_color(p, x, y)[BLU] = 0.0;
1759
1760 }
1761 }
1762 }
1763
1764 if (set_lum_max2 == 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) {
1765
1766 lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
1767 if (lum_ideal > 0 && lum_ideal < setvalue) {
1768 setvalue = lum_ideal;
1769 }
1770 printf
1771 ("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f %f %f %f %f\n",
1772 lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
1773
1774
1775 for (x = 0; x < pict_get_xsize(p); x++)
1776 for (y = 0; y < pict_get_ysize(p); y++) {
1777 if (pict_get_hangle
1778 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1779 if (pict_is_validpixel(p, x, y)) {
1780 lum = luminance(pict_get_color(p, x, y));
1781 if (set_lum_max2 == 1 && lum >= lum_max) {
1782
1783 pict_get_color(p, x, y)[RED] =
1784 setvalue / 179.0;
1785 pict_get_color(p, x, y)[GRN] =
1786 setvalue / 179.0;
1787 pict_get_color(p, x, y)[BLU] =
1788 setvalue / 179.0;
1789
1790 }
1791 }
1792 }
1793 }
1794
1795 pict_write(p, file_out2);
1796
1797 }
1798 if (set_lum_max == 1) {
1799 pict_write(p, file_out2);
1800
1801 }
1802
1803 if (calc_vill == 1) {
1804 printf("%f\n", E_v);
1805 exit(1);
1806 }
1807 }else{
1808 /* in fast calculation mode: ev=ev_ext and sang=2*pi */
1809 sang=2*3.14159265359;
1810 lum_task =E_vl_ext/sang;
1811 E_v=E_vl_ext;
1812 /* printf("calc fast!! %f %f\n", sang,lum_task);*/
1813
1814
1815 }
1816
1817 /* cut out GUTH field of view for glare evaluation */
1818 if (cut_view==1) {
1819 if (cut_view_type==1) {
1820 cut_view_1(p);
1821 }else{
1822
1823 if (cut_view_type==2) {
1824 cut_view_2(p);
1825 }else{
1826 if (cut_view_type==3) {
1827 fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
1828 cut_view_3(p);
1829 }else{
1830 fprintf(stderr,"error: no valid option for view cutting!!");
1831 exit(1);
1832 }
1833 }}
1834 }
1835
1836 if (calcfast == 0) {
1837 avlum = avlum / sang;
1838 lum_task = avlum;
1839 }
1840 /* if (ext_vill==1) {
1841 E_v=E_vl_ext;
1842 avlum=E_v/3.1415927;
1843 } */
1844
1845 if (task_lum == 1) {
1846 lum_task = get_task_lum(p, xt, yt, omegat, task_color);
1847 }
1848 lum_source = lum_thres * lum_task;
1849 /* printf("Task Luminance=%f\n",lum_task);
1850 pict_write(p,"task.pic");*/
1851
1852 if (lum_thres > 100) {
1853 lum_source = lum_thres;
1854 }
1855
1856 /* printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
1857
1858 /* first glare source scan: find primary glare sources */
1859 for (x = 0; x < pict_get_xsize(p); x++)
1860 for (y = 0; y < pict_get_ysize(p); y++) {
1861 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
1862 if (pict_is_validpixel(p, x, y)) {
1863 act_lum = luminance(pict_get_color(p, x, y));
1864 if (act_lum > lum_source) {
1865 if (act_lum >lum_total_max) {
1866 lum_total_max=act_lum;
1867 }
1868
1869 actual_igs =
1870 find_near_pgs(p, x, y, max_angle, 0, igs, 1);
1871 if (actual_igs == 0) {
1872 igs = igs + 1;
1873 pict_new_gli(p);
1874 pict_get_lum_min(p, igs) = HUGE_VAL;
1875 pict_get_Eglare(p,igs) = 0.0;
1876 pict_get_Dglare(p,igs) = 0.0;
1877 actual_igs = igs;
1878 }
1879 icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
1880 pict_get_gsn(p, x, y) = actual_igs;
1881 pict_get_pgs(p, x, y) = 1;
1882 add_pixel_to_gs(p, x, y, actual_igs);
1883
1884 } else {
1885 pict_get_pgs(p, x, y) = 0;
1886 }
1887 }
1888 }
1889 }
1890 /* pict_write(p,"firstscan.pic"); */
1891
1892 if (calcfast == 1 ) {
1893 skip_second_scan=1;
1894 }
1895
1896 /* second glare source scan: combine glare sources facing each other */
1897 change = 1;
1898 while (change == 1 && skip_second_scan==0) {
1899 change = 0;
1900 for (x = 0; x < pict_get_xsize(p); x++)
1901 for (y = 0; y < pict_get_ysize(p); y++) {
1902 before_igs = pict_get_gsn(p, x, y);
1903 if (pict_get_hangle
1904 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1905 if (pict_is_validpixel(p, x, y) && before_igs > 0) {
1906 actual_igs =
1907 find_near_pgs(p, x, y, max_angle, before_igs,
1908 igs, 1);
1909 if (!(actual_igs == before_igs)) {
1910 change = 1;
1911 }
1912 if (before_igs > 0) {
1913 actual_igs = pict_get_gsn(p, x, y);
1914 icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
1915 }
1916
1917 }
1918 }
1919 }
1920 /* pict_write(p,"secondscan.pic");*/
1921
1922 }
1923
1924 /*smoothing: add secondary glare sources */
1925 i_max = igs;
1926 if (sgs == 1) {
1927
1928 /* simplified search radius: constant for entire picture, alway circle not an angle due to performance */
1929 x = (pict_get_xsize(p) / 2);
1930 y = (pict_get_ysize(p) / 2);
1931 rx = (pict_get_xsize(p) / 2 + 10);
1932 ry = (pict_get_ysize(p) / 2 + 10);
1933 r_center =
1934 acos(DOT(pict_get_cached_dir(p, x, y),
1935 pict_get_cached_dir(p, rx, ry))) * 2 / 10;
1936 search_pix = max_angle / r_center / 1.75;
1937
1938 for (x = 0; x < pict_get_xsize(p); x++) {
1939 for (y = 0; y < pict_get_ysize(p); y++) {
1940 if (pict_get_hangle
1941 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1942 if (pict_is_validpixel(p, x, y)
1943 && pict_get_gsn(p, x, y) == 0) {
1944 for (i = 1; i <= i_max; i++) {
1945
1946 if (pict_get_npix(p, i) > 0) {
1947 add_secondary_gs(p, x, y, search_pix, i, uniform_gs, u_r, u_g , u_b);
1948 }
1949 }
1950
1951 }
1952 }
1953 }
1954 }
1955
1956 }
1957
1958 /* extract extremes from glare sources to extra glare source */
1959 if (splithigh == 1 && lum_total_max>limit) {
1960
1961 r_split = max_angle / 2.0;
1962 for (i = 0; i <= i_max; i++) {
1963
1964 if (pict_get_npix(p, i) > 0) {
1965 l_max = pict_get_lum_max(p, i);
1966 i_splitstart = igs + 1;
1967 if (l_max >= limit) {
1968 for (x = 0; x < pict_get_xsize(p); x++)
1969 for (y = 0; y < pict_get_ysize(p); y++) {
1970 if (pict_get_hangle
1971 (p, x, y, p->view.vdir, p->view.vup, &ang))
1972 {
1973 if (pict_is_validpixel(p, x, y)
1974 && luminance(pict_get_color(p, x, y))
1975 >= limit
1976 && pict_get_gsn(p, x, y) == i) {
1977 if (i_splitstart == (igs + 1)) {
1978 igs = igs + 1;
1979 pict_new_gli(p);
1980 pict_get_Eglare(p,igs) = 0.0;
1981 pict_get_Dglare(p,igs) = 0.0;
1982 pict_get_lum_min(p, igs) =
1983 99999999999999.999;
1984 i_split = igs;
1985 } else {
1986 i_split =
1987 find_split(p, x, y, r_split,
1988 i_splitstart, igs);
1989 }
1990 if (i_split == 0) {
1991 igs = igs + 1;
1992 pict_new_gli(p);
1993 pict_get_Eglare(p,igs) = 0.0;
1994 pict_get_Dglare(p,igs) = 0.0;
1995 pict_get_lum_min(p, igs) =
1996 99999999999999.999;
1997 i_split = igs;
1998 }
1999 split_pixel_from_gs(p, x, y, i_split, uniform_gs, u_r, u_g , u_b);
2000
2001 }
2002 }
2003 }
2004
2005 }
2006 change = 1;
2007 while (change == 1) {
2008 change = 0;
2009 for (x = 0; x < pict_get_xsize(p); x++)
2010 for (y = 0; y < pict_get_ysize(p); y++) {
2011 before_igs = pict_get_gsn(p, x, y);
2012 if (before_igs >= i_splitstart) {
2013 if (pict_get_hangle
2014 (p, x, y, p->view.vdir, p->view.vup,
2015 &ang)) {
2016 if (pict_is_validpixel(p, x, y)
2017 && before_igs > 0) {
2018 actual_igs =
2019 find_near_pgs(p, x, y,
2020 max_angle,
2021 before_igs, igs,
2022 i_splitstart);
2023 if (!(actual_igs == before_igs)) {
2024 change = 1;
2025 }
2026 if (before_igs > 0) {
2027 actual_igs =
2028 pict_get_gsn(p, x, y);
2029 icol =
2030 setglcolor(p, x, y,
2031 actual_igs, uniform_gs, u_r, u_g , u_b);
2032 }
2033
2034 }
2035 }
2036 }
2037
2038 }
2039 }
2040
2041
2042 }
2043 }
2044 }
2045
2046 /* calculation of direct vertical illuminance for CGI and for disability glare*/
2047
2048
2049 if (calcfast == 0) {
2050 for (x = 0; x < pict_get_xsize(p); x++)
2051 for (y = 0; y < pict_get_ysize(p); y++) {
2052 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2053 if (pict_is_validpixel(p, x, y)) {
2054 if (pict_get_gsn(p, x, y) > 0) {
2055 pict_get_Eglare(p,pict_get_gsn(p, x, y)) +=
2056 DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2057 * pict_get_omega(p, x, y)
2058 * luminance(pict_get_color(p, x, y));
2059 E_v_dir +=
2060 DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2061 * pict_get_omega(p, x, y)
2062 * luminance(pict_get_color(p, x, y));
2063 }
2064 }
2065 }
2066 }
2067 lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2068 lum_backg = lum_backg_cos;
2069 }
2070 /* calc of band luminance if applied */
2071 if (band_calc == 1) {
2072 band_avlum=get_band_lum( p, band_angle, band_color);
2073 }
2074
2075 /*printf("total number of glare sources: %i\n",igs); */
2076 lum_sources = 0;
2077 omega_sources = 0;
2078 for (x = 0; x <= igs; x++) {
2079 lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2080 omega_sources += pict_get_av_omega(p, x);
2081 }
2082 if (non_cos_lb == 1) {
2083 lum_backg =
2084 (sang * avlum - lum_sources) / (sang - omega_sources);
2085 }
2086 /* print detailed output */
2087 if (detail_out == 1) {
2088 i = 0;
2089 for (x = 0; x <= igs; x++) {
2090 /* resorting glare source numbers */
2091 if (pict_get_npix(p, x) > 0) {
2092 i = i + 1;
2093 }
2094 }
2095
2096
2097
2098 printf
2099 ("%i No pixels x-pos y-pos L_s Omega_s Posindx L_b L_t E_vert Edir Max_Lum Sigma xdir ydir zdir Eglare_cie Lveil_cie \n",
2100 i);
2101 if (i == 0) {
2102 printf("%i %f %f %f %f %.10f %f %f %f %f %f %f\n", i, 0.0, 0.0,
2103 0.0, 0.0, 0.0, 0.0, lum_backg, lum_task, E_v, E_v_dir,
2104 abs_max);
2105
2106 } else {
2107 i = 0;
2108 for (x = 0; x <= igs; x++) {
2109 if (pict_get_npix(p, x) > 0) {
2110 i = i + 1;
2111 pict_get_sigma(p, pict_get_av_posx(p, x),
2112 pict_get_av_posy(p, x), p->view.vdir,
2113 p->view.vup, &sigma);
2114
2115 x2=pict_get_av_posx(p, x);
2116 y2=pict_get_av_posy(p, x);
2117 teta = 180.0 / 3.1415927 * acos(DOT(p->view.vdir, pict_get_cached_dir(p, x2, y2)));
2118 Lveil_cie = 10*pict_get_Eglare(p,x)/(teta*teta+0.0000000000000001);
2119
2120 if (apply_disability == 1 && Lveil_cie <=disability_thresh) {
2121 Lveil_cie =0 ;
2122 }
2123 Lveil_cie_sum = Lveil_cie_sum + Lveil_cie;
2124 printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
2125 i, pict_get_npix(p, x), pict_get_av_posx(p, x),
2126 pict_get_ysize(p) - pict_get_av_posy(p, x),
2127 pict_get_av_lum(p, x), pict_get_av_omega(p, x),
2128 get_posindex(p, pict_get_av_posx(p, x),
2129 pict_get_av_posy(p, x),
2130 posindex_2), lum_backg, lum_task,
2131 E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
2132 ,pict_get_cached_dir(p, x2,y2)[0],pict_get_cached_dir(p, x2,y2)[1],pict_get_cached_dir(p, x2,y2)[2],pict_get_Eglare(p,x),Lveil_cie,teta
2133
2134 );
2135 }
2136 }
2137 }
2138 }
2139
2140
2141
2142 /* calculation of indicees */
2143
2144 /* check vertical illuminance range */
2145 E_v2 = E_v;
2146
2147 if (E_v2 < 100) {
2148 fprintf(stderr,
2149 "Notice: Vertical illuminance is below 100 lux !!\n");
2150 }
2151 dgp =
2152 get_dgp(p, E_v2, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
2153 /* low light correction */
2154 if (E_v < 1000) {
2155 low_light_corr=1.0*exp(0.024*E_v-4)/(1+exp(0.024*E_v-4));} else {low_light_corr=1.0 ;}
2156 dgp =low_light_corr*dgp;
2157
2158 /* age correction */
2159
2160 if (age_corr == 1) {
2161 age_corr_factor=1.0/(1.1-0.5*age/100.0);
2162 }
2163 dgp =age_corr_factor*dgp;
2164
2165
2166 if (ext_vill == 1) {
2167 if (E_vl_ext < 100) {
2168 fprintf(stderr,
2169 "Notice: Vertical illuminance is below 100 lux !!\n");
2170 }
2171 }
2172
2173 if (calcfast == 0) {
2174 dgi = get_dgi(p, lum_backg, igs, posindex_2);
2175 ugr = get_ugr(p, lum_backg, igs, posindex_2);
2176 cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
2177 vcp = get_vcp(p, avlum, igs, posindex_2);
2178 Lveil = get_disability(p, avlum, igs);
2179 }
2180 /* check dgp range */
2181 if (dgp <= 0.2) {
2182 fprintf(stderr,
2183 "Notice: Low brightness scene. dgp below 0.2! dgp might underestimate glare sources\n");
2184 }
2185 if (E_v < 380) {
2186 fprintf(stderr,
2187 "Notice: Low brightness scene. Vertical illuminance less than 380 lux! dgp might underestimate glare sources\n");
2188 }
2189
2190
2191
2192 if (output == 0) {
2193 if (detail_out == 1) {
2194 if (ext_vill == 1) {
2195 dgp_ext =
2196 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2197 posindex_2);
2198 dgp = dgp_ext;
2199 if (E_vl_ext < 1000) {
2200 low_light_corr=1.0*exp(0.024*E_vl_ext-4)/(1+exp(0.024*E_vl_ext-4));} else {low_light_corr=1.0 ;}
2201 dgp =low_light_corr*dgp;
2202 dgp =age_corr_factor*dgp;
2203 }
2204
2205 printf
2206 ("dgp,av_lum,E_v,lum_backg,E_v_dir,dgi,ugr,vcp,cgi,lum_sources,omega_sources,Lveil,Lveil_cie,band_avlum: %f %f %f %f %f %f %f %f %f %f %f %f %f %f \n",
2207 dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
2208 lum_sources / omega_sources, omega_sources, Lveil,Lveil_cie_sum,band_avlum);
2209 } else {
2210 if (detail_out2 == 1) {
2211
2212 printf
2213 ("dgp,dgi,ugr,vcp,cgi,dgp_ext,Ev_calc,abs_max,Lveil: %f %f %f %f %f %f %f %f %f \n",
2214 dgp, dgi, ugr, vcp, cgi, dgp_ext, E_v, abs_max,
2215 Lveil);
2216
2217 } else {
2218 if (ext_vill == 1) {
2219 dgp_ext = get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,posindex_2);
2220
2221 if (E_vl_ext < 1000) {
2222 low_light_corr=1.0*exp(0.024*E_vl_ext-4)/(1+exp(0.024*E_vl_ext-4));} else {low_light_corr=1.0 ;}
2223 dgp =low_light_corr*dgp_ext;
2224 dgp =age_corr_factor*dgp;
2225 }
2226 printf("dgp,dgi,ugr,vcp,cgi,Lveil: %f %f %f %f %f %f \n",
2227 dgp, dgi, ugr, vcp, cgi, Lveil);
2228
2229 }
2230 }
2231 } else {
2232 dgp_ext =
2233 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2234 posindex_2);
2235 dgp = dgp_ext;
2236 if (E_vl_ext < 1000) {
2237 low_light_corr=1.0*exp(0.024*E_vl_ext-4)/(1+exp(0.024*E_vl_ext-4)); } else {low_light_corr=1.0 ;}
2238 dgp =low_light_corr*dgp;
2239 dgp =age_corr_factor*dgp;
2240 printf("%f\n", dgp);
2241 }
2242
2243
2244 /*printf("lowlight=%f\n", low_light_corr); */
2245
2246
2247 /* printf("hallo \n");
2248
2249
2250 apply_disability = 1;
2251 disability_thresh = atof(argv[++i]);
2252 coloring of disability glare sources red, remove all other colors
2253
2254
2255
2256 this function was removed because of bugs....
2257 has to be re-written from scratch....
2258
2259
2260 */
2261
2262
2263
2264
2265
2266
2267
2268
2269 /*output */
2270 if (checkfile == 1) {
2271 pict_write(p, file_out);
2272 }
2273
2274
2275
2276 return EXIT_SUCCESS;
2277 exit(0);
2278
2279 userr:
2280 fprintf(stderr,
2281 "Usage: %s [-s][-d][-c picture][-t xpos ypos angle] [-T xpos ypos angle] [-b fact] [-r angle] [-y] [-Y lum] [-i Ev] [-I Ev ymax ymin] [-v] picfile\n",
2282 progname);
2283 exit(1);
2284 }
2285
2286
2287
2288 #endif