ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.1
Committed: Wed Aug 12 23:07:59 2015 UTC (8 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message: 
Added Jan Wienold's evalglare to distribution

File Contents

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