ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.9
Committed: Fri Nov 23 19:32:17 2018 UTC (5 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.8: +59 -14 lines
Log Message: 
Bug fixes and man page updates from Jan Wienold

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id$";
3 #endif
4 /* EVALGLARE V2.07
5 * Evalglare Software License, Version 2.0
6 *
7 * Copyright (c) 1995 - 2016 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 and EPFL 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," "EPFL" 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 EPFL.
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 EPFL. 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 /* evalglare.c, v1.18 2015/11/10 add masking, UDP, PSGV, DGI_mod, UGR_exp and zoning. For zoning and band luminance, a statistical package from C. Reetz was added in order to derive median, std_dev, percentiles.
282 memory leakage was checked and is o.k.
283 */
284
285 /* evalglare.c, v1.19 2015/12/04 sorting subroutines in statistical package changed by. C. Reetz, adding Position index weighted Luminance evaluation (mean and median of entire image, only in detailed output available)
286 */
287 /* evalglare.c, v1.20 2015/12/21 removing -ansi gcc-option in makefile for the standalone version, adding DGR as output since VCP equation seems not to match with original publication
288 */
289 /* evalglare.c, v1.21 2016/03/15 add a new pixel-overflow-correction option -N
290 */
291 /* evalglare.c, v1.22 2016/03/22 fixed problem using -N option and angular distance calculation, when actual pixel equals disk-center (arccos returned nan)
292 */
293 /* evalglare.c, v1.23 2016/04/18 fixed problem making black corners for -vth fish-eye images
294 */
295 /* evalglare.c, v1.24 2016/04/26 significant speed improvement for 2nd glare source scan: for other glare sources will not be searched any more, when the 8 surrounding pixels have the same glare source number! Smaller speed improvement in the first glare source scan: remembering if pixel before was also a glare source, then no search for surrounding glare sources is applied.
296 changed masking threshold to 0.05 cd/m2
297 */
298 /* evalglare.c, v1.25 2016/04/27 removed the change on the first glare source scan: in few cases glare sources are not merged together in the same way as before. Although the difference is marginal, this algorithm was removed in order to be consistent to the existing results.
299 */
300 /* evalglare.c, v1.26 2016/04/28 removed the bug Lb getting nan in case all glare source pixels are above the peak extraction limit.
301 The calculation of the average source lumiance and solid angle was adapted to be more robust for extreme cases.
302 In case no glare source is found, the values of the glare metrics, relying only on glare sources, is set to zero (100 for VCP) to avoid nan and -inf in the output.
303 Changed glare section output when 0 glare source are found to have the same amount of columns than when glare sources are found
304 */
305 /* evalglare.c, v1.27 2016/06/13 include a warning, if -vtv is in the header. Check, if the corners are black AND -vtv is used ->stopping (stopping can be avoided by using the forcing option -f ). Check, if an invalid exposure is in the header -> stopping in any case.
306 */
307 /* evalglare.c, v1.28 2016/07/08 small code changes in the section of calculating the E_glare (disability) for each glare source.
308 */
309 /* evalglare.c, v1.29 2016/07/12 change threshold for the black corner to 2cd/m2.
310 */
311 /* evalglare.c, v1.30 2016/07/28 change zonal output: If just one zone is specified, only one zone shows up in the output (bug removal).
312 */
313 /* evalglare.c, v1.31 2016/08/02 bug removal: default output did not calculate the amout of glare sources before and therefore no_glaresources was set to zero causing dgi,ugr being set to zero as well. Now renumbering of the glare sources and calculation of the amount of glare sources is done for all output versions.
314 */
315 /* evalglare.c, v2.00 2016/11/15 add of a second fast calculation mode for annual calculations, activted by -2. Output: dgp,ugr
316 */
317 /* evalglare.c, v2.01 2016/11/16 change of -2 option (now -2 dir_illum). External provision of the direct illuminance necessary, since the sun interpolation of daysim is causing problems in calculation of the background luminance.
318 */
319 /* evalglare.c, v2.02 2017/02/28 change of warning message, when invalid exposure setting is found. Reason: tab removal is not in all cases the right measure - it depends which tool caused the invalid exposure entry */
320
321 /* evalglare.c, v2.03 2017/08/04 ad of -O option - disk replacement by providing luminance, not documented
322 */
323
324 /* evalglare.c, v2.04 2017/08/04 adding -q option: use of Ev/pi as background luminance. not documented. no combination with -n option!!!
325 */
326
327 /* evalglare.c, v2.05 2018/08/28 change of the -q option for the choice of the background luminance calculation mode: 0: CIE method (Ev-Edir)/pi, 1: mathematical correct average background luminance, 2: Ev/pi.
328 change of default options:
329 - cosinus weighted calculation of the background luminance (according to CIE) is now default.
330 - absolute threshold for the glare source detection is now default (2000cd/m2), based on study of C. Pierson
331 */
332
333 /* evalglare.c, v2.06 2018/08/29
334 change of default value of multiplier b to 5.0, if task options (-t or -T ) are activated AND -b NOT used. To be downward compatible when using the task method.
335 */
336
337 /* evalglare.c, v2.07 2018/11/06
338 bugfix: correction of error in the equations of PGSV_con and PGSV_sat
339 all three PGSV equations are calculated now
340 illuminance from the masking area (E_v_mask) is also printed
341 bugfix: in VCPs error fuction equation, value of 6.347 replaced by 6.374
342 */
343
344 #define EVALGLARE
345 #define PROGNAME "evalglare"
346 #define VERSION "2.07 release 17.11.2018 by EPFL, J.Wienold"
347 #define RELEASENAME PROGNAME " " VERSION
348
349
350 #include "pictool.h"
351 #include "rtio.h"
352 #include <math.h>
353 #include <string.h>
354 #include "platform.h"
355 #include "muc_randvar.h"
356
357 char *progname;
358
359 /* subroutine to add a pixel to a glare source */
360 void add_pixel_to_gs(pict * p, int x, int y, int gsn)
361 {
362 double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
363 new_omega, act_lum;
364
365
366 pict_get_npix(p, gsn) = pict_get_npix(p, gsn) + 1;
367 old_av_posx = pict_get_av_posx(p, gsn);
368 old_av_posy = pict_get_av_posy(p, gsn);
369 old_av_lum = pict_get_av_lum(p, gsn);
370 old_omega = pict_get_av_omega(p, gsn);
371
372 act_omega = pict_get_omega(p, x, y);
373 act_lum = luminance(pict_get_color(p, x, y));
374 new_omega = old_omega + act_omega;
375 pict_get_av_posx(p, gsn) =
376 (old_av_posx * old_omega + x * act_omega) / new_omega;
377 pict_get_av_posy(p, gsn) =
378 (old_av_posy * old_omega + y * act_omega) / new_omega;
379 if (isnan((pict_get_av_posx(p, gsn))))
380 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);
381
382 pict_get_av_lum(p, gsn) =
383 (old_av_lum * old_omega + act_lum * act_omega) / new_omega;
384 pict_get_av_omega(p, gsn) = new_omega;
385 pict_get_gsn(p, x, y) = gsn;
386 if (act_lum < pict_get_lum_min(p, gsn)) {
387 pict_get_lum_min(p, gsn) = act_lum;
388 }
389 if (act_lum > pict_get_lum_max(p, gsn)) {
390 pict_get_lum_max(p, gsn) = act_lum;
391 }
392
393 /* 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)); */
394
395 }
396
397 /* subroutine for peak extraction */
398 int
399 find_split(pict * p, int x, int y, double r, int i_split_start,
400 int i_split_end)
401 {
402 int i_find_split, x_min, x_max, y_min, y_max, ix, iy, iix, iiy, dx, dy,
403 out_r;
404 double r_actual;
405
406 i_find_split = 0;
407
408 x_min = 0;
409 y_min = 0;
410 x_max = pict_get_ysize(p) - 1;
411
412 y_max = pict_get_ysize(p) - 1;
413
414 for (iiy = 1; iiy <= 2; iiy++) {
415 dy = iiy * 2 - 3;
416 if (dy == -1) {
417 iy = y;
418 } else {
419 iy = y + 1;
420 }
421 while (iy <= y_max && iy >= y_min) {
422 out_r = 0;
423 for (iix = 1; iix <= 2; iix++) {
424 dx = iix * 2 - 3;
425 if (dx == -1) {
426 ix = x;
427 } else {
428 ix = x + 1;
429 }
430 while (ix <= x_max && ix >= x_min && iy >= y_min) {
431
432 r_actual =
433 acos(DOT(pict_get_cached_dir(p, x, y),
434 pict_get_cached_dir(p, ix, iy))) * 2;
435 if (r_actual <= r) {
436 out_r = 1;
437 if (pict_get_gsn(p, ix, iy) >= i_split_start
438 && pict_get_gsn(p, ix, iy) <= i_split_end) {
439 i_find_split = pict_get_gsn(p, ix, iy);
440 }
441 } else {
442 ix = -99;
443 }
444 ix = ix + dx;
445 }
446 }
447 if (out_r == 0) {
448 iy = -99;
449 }
450 iy = iy + dy;
451
452 }
453 }
454
455
456
457 return i_find_split;
458 }
459
460 /*
461 static int icomp(int* a,int* b)
462 {
463 if (*a < *b)
464 return -1;
465 if (*a > *b)
466 return 1;
467 return 0;
468 }
469
470 */
471
472 /* subroutine to find nearby glare source pixels */
473
474 int
475 find_near_pgs(pict * p, int x, int y, float r, int act_gsn, int max_gsn,
476 int min_gsn)
477 {
478 int dx, dy, i_near_gs, xx, yy, x_min, x_max, y_min, y_max, ix, iy, iix,
479 iiy, old_gsn, new_gsn, find_gsn, change, stop_y_search,
480 stop_x_search;
481 double r_actual;
482 int ixm[3];
483
484 i_near_gs = 0;
485 stop_y_search = 0;
486 stop_x_search = 0;
487 x_min = 0;
488 y_min = 0;
489 if (act_gsn == 0) {
490 x_max = x;
491 } else {
492 x_max = pict_get_xsize(p) - 1;
493 }
494
495 y_max = pict_get_ysize(p) - 1;
496
497 old_gsn = pict_get_gsn(p, x, y);
498 new_gsn = old_gsn;
499 change = 0;
500 if (act_gsn > 0) {
501 i_near_gs = pict_get_gsn(p, x, y);
502 }
503 for (iiy = 1; iiy <= 2; iiy++) {
504 dy = iiy * 2 - 3;
505 if (dy == -1) {
506 iy = y;
507 } else {
508 iy = y + 1;
509 }
510 ixm[1] = x;
511 ixm[2] = x;
512 stop_y_search = 0;
513
514
515 while (iy <= y_max && iy >= y_min) {
516 for (iix = 1; iix <= 2; iix++) {
517 dx = iix * 2 - 3;
518 ix = (ixm[1] + ixm[2]) / 2;
519 stop_x_search = 0;
520 while (ix <= x_max && ix >= x_min && stop_x_search == 0
521 && stop_y_search == 0) {
522 /* 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);*/
523 r_actual =
524 acos(DOT(pict_get_cached_dir(p, x, y),
525 pict_get_cached_dir(p, ix, iy))) * 2;
526 if (r_actual <= r) {
527 if (pict_get_gsn(p, ix, iy) > 0) {
528 if (act_gsn == 0) {
529 i_near_gs = pict_get_gsn(p, ix, iy);
530 stop_x_search = 1;
531 stop_y_search = 1;
532 } else {
533 find_gsn = pict_get_gsn(p, ix, iy);
534 if (pict_get_av_omega(p, old_gsn) <
535 pict_get_av_omega(p, find_gsn)
536 && pict_get_av_omega(p,
537 find_gsn) >
538 pict_get_av_omega(p, new_gsn)
539 && find_gsn >= min_gsn
540 && find_gsn <= max_gsn) {
541 /* other primary glare source found with larger solid angle -> add together all */
542 new_gsn = find_gsn;
543 change = 1;
544 stop_x_search = 1;
545 stop_y_search = 1;
546 }
547 }
548 }
549 } else {
550 ixm[iix] = ix - dx;
551 stop_x_search = 1;
552 }
553 ix = ix + dx;
554 }
555 }
556 iy = iy + dy;
557 }
558 }
559 if (change > 0) {
560 pict_get_av_lum(p, old_gsn) = 0.;
561 pict_get_av_omega(p, old_gsn) = 0.;
562 pict_get_npix(p, old_gsn) = 0.;
563 pict_get_lum_max(p, old_gsn) = 0.;
564
565
566 i_near_gs = new_gsn;
567 /* printf(" changing gs no %i",old_gsn);
568 printf(" to %i\n",new_gsn);
569 */
570 for (xx = 0; xx < pict_get_xsize(p); xx++)
571 for (yy = 0; yy < pict_get_ysize(p); yy++) {
572 if (pict_is_validpixel(p, x, y)) {
573 if (pict_get_gsn(p, xx, yy) == old_gsn) {
574 add_pixel_to_gs(p, xx, yy, new_gsn);
575 }
576 }
577 }
578 }
579
580 return i_near_gs;
581
582
583 }
584
585 /* subroutine for calculation of task luminance */
586 double get_task_lum(pict * p, int x, int y, float r, int task_color)
587 {
588 int x_min, x_max, y_min, y_max, ix, iy;
589 double r_actual, av_lum, omega_sum, act_lum;
590
591
592 x_max = pict_get_xsize(p) - 1;
593 y_max = pict_get_ysize(p) - 1;
594 x_min = 0;
595 y_min = 0;
596
597
598
599 av_lum = 0.0;
600 omega_sum = 0.0;
601
602 for (iy = y_min; iy <= y_max; iy++) {
603 for (ix = x_min; ix <= x_max; ix++) {
604
605 /* if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
606 continue;*/
607 r_actual =
608 acos(DOT(pict_get_cached_dir(p, x, y),
609 pict_get_cached_dir(p, ix, iy))) * 2;
610 act_lum = luminance(pict_get_color(p, ix, iy));
611
612 if (r_actual <= r) {
613 act_lum = luminance(pict_get_color(p, ix, iy));
614 av_lum += pict_get_omega(p, ix, iy) * act_lum;
615 omega_sum += pict_get_omega(p, ix, iy);
616 if (task_color == 1) {
617 pict_get_color(p, ix, iy)[RED] = 0.0;
618 pict_get_color(p, ix, iy)[GRN] = 0.0;
619 pict_get_color(p, ix, iy)[BLU] =
620 act_lum / WHTEFFICACY / CIE_bf;
621 }
622 }
623 }
624 }
625
626 av_lum = av_lum / omega_sum;
627 /* printf("average luminance of task %f \n",av_lum);
628 printf("solid angle of task %f \n",omega_sum);*/
629 return av_lum;
630 }
631
632
633
634
635
636 /* subroutine for coloring the glare sources
637 red is used only for explicit call of the subroutine with acol=0 , e.g. for disability glare
638 -1: set to grey*/
639 int setglcolor(pict * p, int x, int y, int acol, int uniform_gs, double u_r, double u_g ,double u_b)
640 {
641 int icol;
642 double pcol[3][1000], act_lum, l;
643
644 l=u_r+u_g+u_b ;
645
646 pcol[RED][0] = 1.0 / CIE_rf;
647 pcol[GRN][0] = 0.0 / CIE_gf;
648 pcol[BLU][0] = 0.0 / CIE_bf;
649
650 pcol[RED][1] = 0.0 / CIE_rf;
651 pcol[GRN][1] = 1.0 / CIE_gf;
652 pcol[BLU][1] = 0.0 / CIE_bf;
653
654 pcol[RED][2] = 0.15 / CIE_rf;
655 pcol[GRN][2] = 0.15 / CIE_gf;
656 pcol[BLU][2] = 0.7 / CIE_bf;
657
658 pcol[RED][3] = .5 / CIE_rf;
659 pcol[GRN][3] = .5 / CIE_gf;
660 pcol[BLU][3] = 0.0 / CIE_bf;
661
662 pcol[RED][4] = .5 / CIE_rf;
663 pcol[GRN][4] = .0 / CIE_gf;
664 pcol[BLU][4] = .5 / CIE_bf ;
665
666 pcol[RED][5] = .0 / CIE_rf;
667 pcol[GRN][5] = .5 / CIE_gf;
668 pcol[BLU][5] = .5 / CIE_bf;
669
670 pcol[RED][6] = 0.333 / CIE_rf;
671 pcol[GRN][6] = 0.333 / CIE_gf;
672 pcol[BLU][6] = 0.333 / CIE_bf;
673
674
675 pcol[RED][999] = 1.0 / WHTEFFICACY;
676 pcol[GRN][999] = 1.0 / WHTEFFICACY;
677 pcol[BLU][999] = 1.0 / WHTEFFICACY;
678
679
680 pcol[RED][998] = u_r /(l* CIE_rf) ;
681 pcol[GRN][998] = u_g /(l* CIE_gf);
682 pcol[BLU][998] = u_b /(l* CIE_bf);
683 /* printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l,WHTEFFICACY); */
684 icol = acol ;
685 if ( acol == -1) {icol=999;
686 }else{if (acol>0){icol = acol % 5 +1;
687 }};
688 if ( uniform_gs == 1) {icol=998;
689 } ;
690
691 act_lum = luminance(pict_get_color(p, x, y));
692
693 pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY;
694 pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY;
695 pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY;
696 /* printf("x,y,lum_before,lum_after,diff %i %i %f %f %f\n",x,y, act_lum,luminance(pict_get_color(p, x, y)), act_lum-luminance(pict_get_color(p, x, y))); */
697 return icol;
698 }
699
700 /* this is the smooting subroutine */
701 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)
702 {
703 int x_min, x_max, y_min, y_max, ix, iy, icol;
704 double r_actual, omega_gs, omega_total, om;
705
706
707
708 omega_gs = 0.0;
709 omega_total = 0.0;
710 x_min = x - r;
711 if (x_min < 0) {
712 x_min = 0;
713 }
714 x_max = x + r;
715 if (x_max > pict_get_xsize(p) - 1) {
716 x_max = pict_get_xsize(p) - 2;
717 }
718 y_min = y - r;
719 if (y_min < 0) {
720 y_min = 0;
721 }
722 y_max = y + r;
723 if (y_max > pict_get_ysize(p) - 1) {
724 y_max = pict_get_ysize(p) - 2;
725 }
726
727 for (ix = x_min; ix <= x_max; ix++)
728 for (iy = y_min; iy <= y_max; iy++) {
729 r_actual = sqrt((x - ix) * (x - ix) + (y - iy) * (y - iy));
730 if (r_actual <= r) {
731 om = pict_get_omega(p, ix, iy);
732 omega_total += om;
733 if (pict_get_gsn(p, ix, iy) == act_gs
734 && pict_get_pgs(p, ix, iy) == 1) {
735 omega_gs = omega_gs + 1 * om;
736 }
737 }
738 }
739
740
741 if (omega_gs / omega_total > 0.2) {
742 /* add pixel to gs */
743
744 add_pixel_to_gs(p, x, y, act_gs);
745
746 /* color pixel of gs */
747
748 /* icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b); */
749
750 }
751 }
752
753 /* subroutine for removing a pixel from a glare source */
754 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)
755 {
756 int old_gsn, icol;
757 double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
758 new_omega, act_lum;
759
760
761 /* change existing gs */
762 old_gsn = pict_get_gsn(p, x, y);
763 pict_get_npix(p, old_gsn) = pict_get_npix(p, old_gsn) - 1;
764
765 act_omega = pict_get_omega(p, x, y);
766 old_av_posx = pict_get_av_posx(p, old_gsn);
767 old_av_posy = pict_get_av_posy(p, old_gsn);
768 old_omega = pict_get_av_omega(p, old_gsn);
769
770 new_omega = old_omega - act_omega;
771 pict_get_av_omega(p, old_gsn) = new_omega;
772 pict_get_av_posx(p, old_gsn) =
773 (old_av_posx * old_omega - x * act_omega) / new_omega;
774 pict_get_av_posy(p, old_gsn) =
775 (old_av_posy * old_omega - y * act_omega) / new_omega;
776
777
778 act_lum = luminance(pict_get_color(p, x, y));
779 old_av_lum = pict_get_av_lum(p, old_gsn);
780 pict_get_av_lum(p, old_gsn) =
781 (old_av_lum * old_omega - act_lum * act_omega) / new_omega;
782 /* add pixel to new gs */
783
784 add_pixel_to_gs(p, x, y, new_gsn);
785
786 /* color pixel of new gs */
787
788 /* icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b); */
789
790
791 }
792
793 /* subroutine for the calculation of the position index */
794 float get_posindex(pict * p, float x, float y, int postype)
795 {
796 float posindex;
797 double teta, phi, sigma, tau, deg, d, s, r, fact;
798
799
800 pict_get_vangle(p, x, y, p->view.vdir, p->view.vup, &phi);
801 pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &teta);
802 pict_get_sigma(p, x, y, p->view.vdir, p->view.vup, &sigma);
803 pict_get_tau(p, x, y, p->view.vdir, p->view.vup, &tau);
804
805 /* return (phi+teta+2*3.1415927); */
806
807 deg = 180 / 3.1415927;
808 fact = 0.8;
809 if (phi == 0) {
810 phi = 0.00001;
811 }
812 if (sigma <= 0) {
813 sigma = -sigma;
814 }
815 if (teta == 0) {
816 teta = 0.0001;
817 }
818 tau = tau * deg;
819 sigma = sigma * deg;
820
821 if (postype == 1) {
822 /* KIM model */
823 posindex = exp ((sigma-(-0.000009*tau*tau*tau+0.0014*tau*tau+0.0866*tau+21.633))/(-0.000009*tau*tau*tau+0.0013*tau*tau+0.0853*tau+8.772));
824 }else{
825 /* Guth model, equation from IES lighting handbook */
826 posindex =
827 exp((35.2 - 0.31889 * tau -
828 1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 +
829 0.26667 * tau -
830 0.002963 * tau *
831 tau) / 100000 *
832 sigma * sigma);
833 /* below line of sight, using Iwata model */
834 if (phi < 0) {
835 d = 1 / tan(phi);
836 s = tan(teta) / tan(phi);
837 r = sqrt(1 / d * 1 / d + s * s / d / d);
838 if (r > 0.6)
839 fact = 1.2;
840 if (r > 3) {
841 fact = 1.2;
842 r = 3;
843 }
844
845 posindex = 1 + fact * r;
846 }
847 if (posindex > 16)
848 posindex = 16;
849 }
850
851 return posindex;
852
853 }
854
855
856
857 double get_upper_cut_2eyes(float teta)
858 {
859 double phi;
860
861 phi=pow(7.7458218+0.00057407915*teta-0.00021746318*teta*teta+8.5572726e-6*teta*teta*teta,2);
862
863 return phi;
864 }
865 double get_lower_cut_2eyes(float teta)
866 {
867 double phi;
868
869 phi=1/(-0.014699242-1.5541106e-5*teta+4.6898068e-6*teta*teta-5.1539687e-8*teta*teta*teta);
870
871 return phi;
872 }
873
874 double get_lower_cut_central(float teta)
875 {
876 double phi;
877
878 phi=(68.227109-2.9524084*teta+0.046674262*teta*teta)/(1-0.042317294*teta+0.00075698419*teta*teta-6.5364285e-7*teta*teta*teta);
879
880 if (teta>73) {
881 phi=60;
882 }
883
884 return phi;
885 }
886
887 /* subroutine for the cutting the total field of view */
888 void cut_view_1(pict*p)
889 {
890 int x,y;
891 double border,ang,teta,phi,phi2;
892 for(x=0;x<pict_get_xsize(p)-1;x++)
893 for(y=0;y<pict_get_ysize(p)-1;y++) {
894 if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
895 if (pict_is_validpixel(p, x, y)) {
896 pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
897 pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
898 pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
899 phi=phi*180/3.1415927;
900 phi2=phi2*180/3.1415927;
901 teta=teta*180/3.1415927;
902 if (teta<0) {
903 teta=-teta;
904 }
905 if(phi2>0){
906 border=get_upper_cut_2eyes(teta);
907
908 if (phi>border) {
909 pict_get_color(p,x,y)[RED]=0;
910 pict_get_color(p,x,y)[GRN]=0;
911 pict_get_color(p,x,y)[BLU]=0;
912 }
913 }else{
914 border=get_lower_cut_2eyes(180-teta);
915 if (-phi<border && teta>135) {
916 pict_get_color(p,x,y)[RED]=0;
917 pict_get_color(p,x,y)[GRN]=0;
918 pict_get_color(p,x,y)[BLU]=0;
919 }
920 }
921
922
923
924 }else{
925 pict_get_color(p,x,y)[RED]=0;
926 pict_get_color(p,x,y)[GRN]=0;
927 pict_get_color(p,x,y)[BLU]=0;
928
929
930 }
931
932 /* printf("teta,phi,border=%f %f %f\n",teta,phi,border);*/
933 }
934 }
935
936
937 return;
938
939 }
940 /* subroutine for the cutting the field of view seen by both eyes*/
941 void cut_view_2(pict*p)
942 {
943
944 int x,y;
945 double border,ang,teta,phi,phi2;
946 for(x=0;x<pict_get_xsize(p)-1;x++)
947 for(y=0;y<pict_get_ysize(p)-1;y++) {
948 if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
949 if (pict_is_validpixel(p, x, y)) {
950 pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
951 pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
952 pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
953 pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
954 phi=phi*180/3.1415927;
955 phi2=phi2*180/3.1415927;
956 teta=teta*180/3.1415927;
957 /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
958 if (teta<0) {
959 teta=-teta;
960 }
961 if(phi2>0){
962 border=60;
963
964 if (phi>border) {
965 pict_get_color(p,x,y)[RED]=0;
966 pict_get_color(p,x,y)[GRN]=0;
967 pict_get_color(p,x,y)[BLU]=0;
968 }
969 }else{
970 border=get_lower_cut_central(180-teta);
971 if (phi>border) {
972 pict_get_color(p,x,y)[RED]=0;
973 pict_get_color(p,x,y)[GRN]=0;
974 pict_get_color(p,x,y)[BLU]=0;
975 }
976 }
977
978
979
980 }else{
981 pict_get_color(p,x,y)[RED]=0;
982 pict_get_color(p,x,y)[GRN]=0;
983 pict_get_color(p,x,y)[BLU]=0;
984
985
986 }
987 }
988 }
989
990
991 return;
992
993 }
994
995 /* subroutine for the cutting the field of view seen by both eyes
996 luminance is modified by position index - just for experiments - undocumented */
997 void cut_view_3(pict*p)
998 {
999
1000 int x,y;
1001 double border,ang,teta,phi,phi2,lum,newlum;
1002 for(x=0;x<pict_get_xsize(p)-1;x++)
1003 for(y=0;y<pict_get_ysize(p)-1;y++) {
1004 if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
1005 if (DOT(pict_get_cached_dir(p,x,y),p->view.vdir) >= 0.0) {
1006 pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
1007 pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
1008 pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
1009 pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
1010 phi=phi*180/3.1415927;
1011 phi2=phi2*180/3.1415927;
1012 teta=teta*180/3.1415927;
1013 lum=luminance(pict_get_color(p,x,y));
1014
1015 newlum=lum/get_posindex(p,x,y,0);
1016
1017 pict_get_color(p,x,y)[RED]=newlum/WHTEFFICACY;
1018 pict_get_color(p,x,y)[GRN]=newlum/WHTEFFICACY;
1019 pict_get_color(p,x,y)[BLU]=newlum/WHTEFFICACY;
1020
1021
1022
1023
1024 /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
1025 if (teta<0) {
1026 teta=-teta;
1027 }
1028 if(phi2>0){
1029 border=60;
1030
1031 if (phi>border) {
1032 pict_get_color(p,x,y)[RED]=0;
1033 pict_get_color(p,x,y)[GRN]=0;
1034 pict_get_color(p,x,y)[BLU]=0;
1035 }
1036 }else{
1037 border=get_lower_cut_central(180-teta);
1038 if (phi>border) {
1039 pict_get_color(p,x,y)[RED]=0;
1040 pict_get_color(p,x,y)[GRN]=0;
1041 pict_get_color(p,x,y)[BLU]=0;
1042 }
1043 }
1044
1045
1046
1047 }else{
1048 pict_get_color(p,x,y)[RED]=0;
1049 pict_get_color(p,x,y)[GRN]=0;
1050 pict_get_color(p,x,y)[BLU]=0;
1051
1052
1053 }
1054 }
1055 }
1056
1057
1058 return;
1059
1060 }
1061
1062 /* subroutine for the calculation of the daylight glare index DGI*/
1063
1064
1065 float get_dgi(pict * p, float lum_backg, int igs, int posindex_2)
1066 {
1067 float dgi, sum_glare, omega_s;
1068 int i;
1069
1070
1071 sum_glare = 0;
1072 omega_s = 0;
1073
1074 for (i = 0; i <= igs; i++) {
1075
1076 if (pict_get_npix(p, i) > 0) {
1077 omega_s =
1078 pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1079 get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1080 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));
1081 /* printf("i,sum_glare %i %f\n",i,sum_glare); */
1082 }
1083 }
1084 dgi = 10 * log10(sum_glare);
1085
1086 return dgi;
1087
1088 }
1089 /* subroutine for the calculation of the modified daylight glare index DGI_mod (Fisekis2003)
1090 several glare sources are taken into account and summed up, original equation has no summary
1091 */
1092
1093
1094 float get_dgi_mod(pict * p, float lum_a, int igs, int posindex_2)
1095 {
1096 float dgi_mod, sum_glare, omega_s;
1097 int i;
1098
1099
1100 sum_glare = 0;
1101 omega_s = 0;
1102
1103 for (i = 0; i <= igs; i++) {
1104
1105 if (pict_get_npix(p, i) > 0) {
1106 omega_s =
1107 pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1108 get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1109 sum_glare += 0.478 * pow(pict_get_av_lum(p, i), 1.6) * pow(omega_s,0.8) / (pow(lum_a,0.85) + 0.07 * pow(pict_get_av_omega(p, i),0.5) * pict_get_av_lum(p, i));
1110 /* printf("i,sum_glare %i %f\n",i,sum_glare); */
1111 }
1112 }
1113 dgi_mod = 10 * log10(sum_glare);
1114
1115 return dgi_mod;
1116
1117 }
1118
1119 /* subroutine for the calculation of the daylight glare probability */
1120 double
1121 get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3,
1122 double a4, double a5, double c1, double c2, double c3,
1123 int posindex_2)
1124 {
1125 double dgp;
1126 double sum_glare;
1127 int i;
1128
1129 sum_glare = 0;
1130 if (igs > 0) {
1131 for (i = 0; i <= igs; i++) {
1132
1133 if (pict_get_npix(p, i) > 0) {
1134 sum_glare +=
1135 pow(pict_get_av_lum(p, i),
1136 a1) / pow(get_posindex(p, pict_get_av_posx(p, i),
1137 pict_get_av_posy(p, i),
1138 posindex_2),
1139 a4) * pow(pict_get_av_omega(p, i), a2);
1140 /* printf("i,sum_glare %i %f\n",i,sum_glare);*/
1141 }
1142 }
1143 dgp =
1144 c1 * pow(E_v,
1145 a5) + c3 + c2 * log10(1 + sum_glare / pow(E_v, a3));
1146 } else {
1147 dgp = c3 + c1 * pow(E_v, a5);
1148 }
1149
1150 if (dgp > 1) {
1151 dgp = 1;
1152 }
1153 /* printf("dgp %f\n",dgp);*/
1154 return dgp;
1155
1156 }
1157
1158 /* subroutine for the calculation of the visual comfort probability */
1159 float get_dgr(pict * p, double lum_a, int igs, int posindex_2)
1160 {
1161 float dgr;
1162 double sum_glare;
1163 int i, i_glare;
1164
1165
1166 sum_glare = 0;
1167 i_glare = 0;
1168 for (i = 0; i <= igs; i++) {
1169
1170 if (pict_get_npix(p, i) > 0) {
1171 i_glare = i_glare + 1;
1172 sum_glare +=
1173 (0.5 * pict_get_av_lum(p, i) *
1174 (20.4 * pict_get_av_omega(p, i) +
1175 1.52 * pow(pict_get_av_omega(p, i),
1176 0.2) - 0.075)) / (get_posindex(p,
1177 pict_get_av_posx
1178 (p, i),
1179 pict_get_av_posy
1180 (p, i),
1181 posindex_2) *
1182 pow(lum_a, 0.44));
1183
1184 }
1185 }
1186 dgr = pow(sum_glare, pow(i_glare, -0.0914));
1187
1188
1189
1190 return dgr;
1191
1192 }
1193
1194 float get_vcp(float dgr )
1195 {
1196 float vcp;
1197
1198 vcp = 50 * erf((6.374 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1199 if (dgr > 750) {
1200 vcp = 0;
1201 }
1202 if (dgr < 20) {
1203 vcp = 100;
1204 }
1205 return vcp;
1206
1207 }
1208
1209
1210 /* subroutine for the calculation of the unified glare rating */
1211 float get_ugr(pict * p, double lum_backg, int igs, int posindex_2)
1212 {
1213 float ugr;
1214 double sum_glare;
1215 int i, i_glare;
1216
1217
1218 sum_glare = 0;
1219 i_glare = 0;
1220 for (i = 0; i <= igs; i++) {
1221
1222 if (pict_get_npix(p, i) > 0) {
1223 i_glare = i_glare + 1;
1224 sum_glare +=
1225 pow(pict_get_av_lum(p, i) /
1226 get_posindex(p, pict_get_av_posx(p, i),
1227 pict_get_av_posy(p, i), posindex_2),
1228 2) * pict_get_av_omega(p, i);
1229
1230 }
1231 }
1232 ugr = 8 * log10(0.25 / lum_backg * sum_glare);
1233 if (sum_glare==0) {
1234 ugr=0.0;
1235 }
1236 if (lum_backg<=0) {
1237 ugr=-99.0;
1238 }
1239
1240 return ugr;
1241
1242 }
1243 /* subroutine for the calculation of the experimental unified glare rating (Fisekis 2003)*/
1244 float get_ugr_exp(pict * p, double lum_backg, double lum_a, int igs, int posindex_2)
1245 {
1246 float ugr_exp;
1247 double sum_glare;
1248 int i, i_glare;
1249
1250
1251 sum_glare = 0;
1252 i_glare = 0;
1253 for (i = 0; i <= igs; i++) {
1254
1255 if (pict_get_npix(p, i) > 0) {
1256 i_glare = i_glare + 1;
1257 sum_glare +=
1258 pow(1 /get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2),2) *pict_get_av_lum(p, i)* pict_get_av_omega(p, i);
1259
1260 }
1261 }
1262 ugr_exp =8 * log10(lum_a) + 8 * log10(sum_glare/lum_backg);
1263
1264 return ugr_exp;
1265
1266 }
1267 /* subroutine for the calculation of the unified glare probability (Hirning 2014)*/
1268 float get_ugp(pict * p, double lum_backg, int igs, int posindex_2)
1269 {
1270 float ugp;
1271 double sum_glare;
1272 int i, i_glare;
1273
1274
1275 sum_glare = 0;
1276 i_glare = 0;
1277 for (i = 0; i <= igs; i++) {
1278
1279 if (pict_get_npix(p, i) > 0) {
1280 i_glare = i_glare + 1;
1281 sum_glare +=
1282 pow(pict_get_av_lum(p, i) /
1283 get_posindex(p, pict_get_av_posx(p, i),
1284 pict_get_av_posy(p, i), posindex_2),
1285 2) * pict_get_av_omega(p, i);
1286
1287 }
1288 }
1289 ugp = 0.26 * log10(0.25 / lum_backg * sum_glare);
1290
1291 return ugp;
1292
1293 }
1294
1295
1296 /* subroutine for the calculation of the disability glare according to Poynter */
1297 float get_disability(pict * p, double lum_backg, int igs)
1298 {
1299 float disab;
1300 double sum_glare, sigma, deg;
1301 int i, i_glare;
1302
1303
1304 sum_glare = 0;
1305 i_glare = 0;
1306 deg = 180 / 3.1415927;
1307 for (i = 0; i <= igs; i++) {
1308
1309 if (pict_get_npix(p, i) > 0) {
1310 i_glare = i_glare + 1;
1311 pict_get_sigma(p, pict_get_av_posx(p, i),
1312 pict_get_av_posy(p, i), p->view.vdir,
1313 p->view.vup, &sigma);
1314
1315 sum_glare +=
1316 pict_get_av_lum(p,
1317 i) * cos(sigma +
1318 0.00000000001) *
1319 pict_get_av_omega(p, i)
1320 / (deg * sigma + 0.00000000001);
1321 if (isnan(sum_glare)) {
1322 printf("sigma for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1323 printf("omega for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1324 printf("avlum for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1325 printf("avlum for %f %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i), sigma);
1326 }
1327
1328 }
1329 }
1330 disab = 5 * sum_glare;
1331
1332 return disab;
1333
1334 }
1335
1336
1337
1338
1339
1340
1341 /* subroutine for the calculation of the cie glare index */
1342
1343 float get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1344 {
1345 float cgi;
1346 double sum_glare;
1347 int i, i_glare;
1348
1349
1350 sum_glare = 0;
1351 i_glare = 0;
1352 for (i = 0; i <= igs; i++) {
1353
1354 if (pict_get_npix(p, i) > 0) {
1355 i_glare = i_glare + 1;
1356 sum_glare +=
1357 pow(pict_get_av_lum(p, i) /
1358 get_posindex(p, pict_get_av_posx(p, i),
1359 pict_get_av_posy(p, i), posindex_2),
1360 2) * pict_get_av_omega(p, i);
1361
1362 }
1363 }
1364 cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare);
1365
1366 return cgi;
1367 }
1368
1369
1370
1371 /* subroutine for the calculation of the PGSV_con; is only applied, when masking is done, since it refers only to the window. Important: masking area must be the window! */
1372 float get_pgsv_con(double E_v, double E_mask,double omega_mask,double lum_mask_av, double Lavg)
1373 {
1374 float pgsv_con;
1375 double Lb;
1376
1377 /* Lb = (E_v-E_mask)/3.14159265359; */
1378 /* Lb = (2*E_v-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask); */
1379 Lb = (2*3.14159265359*Lavg-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask);
1380
1381
1382 pgsv_con =3.2*log10(lum_mask_av)-0.64*log10(omega_mask)+(0.79*log10(omega_mask)-0.61)*log10(Lb)-8.2 ;
1383
1384
1385 return pgsv_con;
1386 }
1387
1388 /* subroutine for the calculation of the PGSV_saturation; is only applied, when masking is done, since it refers only to the window. Important: masking area must be the window! */
1389 float get_pgsv_sat(double E_v)
1390 {
1391 float pgsv_sat;
1392
1393 pgsv_sat =3.3-(0.57+3.3)/(1+pow(E_v/3.14159265359/1250,1.7));
1394
1395
1396 return pgsv_sat;
1397 }
1398
1399 /* subroutine for the calculation of the PGSV; is only applied, when masking is done, since it refers only to the window. Important: masking area must be the window! */
1400
1401 float get_pgsv(double E_v, double E_mask,double omega_mask,double lum_mask_av,double Ltask, double Lavg)
1402 {
1403 float pgsv;
1404 double Lb;
1405
1406 /* Lb = (E_v-E_mask)/3.14159265359; */
1407 /* Lb = (2*E_v-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask); */
1408 Lb = (2*3.14159265359*Lavg-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask);
1409
1410 if (Lb==0.0 ) {
1411 fprintf(stderr, " warning: Background luminance is 0 or masking area = full image! pgsv cannot be calculated (set to -99)!!\n");
1412 pgsv=-99;
1413 }else{
1414 if ( (lum_mask_av/Lb) > (E_v/(3.14159265359*Ltask))) {
1415 pgsv=get_pgsv_con(E_v,E_mask,omega_mask,lum_mask_av, Lavg);
1416 }else{
1417 pgsv=get_pgsv_sat(E_v) ;
1418 }}
1419 return pgsv;
1420
1421 }
1422
1423
1424
1425 #ifdef EVALGLARE
1426
1427
1428 /* main program
1429 ------------------------------------------------------------------------------------------------------------------*/
1430
1431
1432 int main(int argc, char **argv)
1433 {
1434 #define CLINEMAX 4095 /* memory allocated for command line string */
1435 pict *p = pict_create();
1436 pict *pm = pict_create();
1437 int skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, x1,y1, fill, yfillmax, yfillmin,
1438 ext_vill, set_lum_max, set_lum_max2, img_corr,x_disk,y_disk,task_color, i_splitstart,zones,act_gsn,splitgs,
1439 i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,x_zone,y_zone, i_z1, i_z2, thres_activate,
1440 igs, actual_igs, lastpixelwas_gs, icol, xt, yt, change,checkpixels, before_igs, sgs, splithigh,uniform_gs,x_max, y_max,y_mid,
1441 detail_out, posindex_picture, non_cos_lb, rx, ry, rmx,rmy,apply_disability,band_calc,band_color,masking,i_mask,no_glaresources,force;
1442 double LUM_replace,lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma,om,delta_E,
1443 E_vl_ext, lum_max, new_lum_max, r_center, ugp, ugr_exp, dgi_mod,lum_a, E_v_mask,angle_disk,dist,n_corner_px,zero_corner_px,
1444 search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,r_actual,lum_actual,dir_ill,
1445 omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang, angle_z1, angle_z2,per_95_band,per_75_band,pos,
1446 l_max, lum_backg, lum_backg_cos, omega_sources, lum_sources,per_75_mask,per_95_mask,per_75_z1,per_95_z1,per_75_z2,per_95_z2,
1447 lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum,
1448 lum_mask[1],lum_mask_av,lum_mask_std[1],lum_mask_median[1],omega_mask,bbox[2],
1449 lum_band[1],lum_band_av,lum_band_std[1],lum_band_median[1],omega_band,bbox_band[2],
1450 lum_z1[1],lum_z1_av,lum_z1_std[1],lum_z1_median[1],omega_z1,bbox_z1[2],
1451 lum_z2[1],lum_z2_av,lum_z2_std[1],lum_z2_median[1],omega_z2,bbox_z2[2],
1452 lum_pos[1],lum_nopos_median[1],lum_pos_median[1],lum_pos2_median[1],lum_pos_mean,lum_pos2_mean;
1453 float lum_task, lum_thres, dgi, vcp, cgi, ugr, limit, dgr,pgsv ,pgsv_sat,pgsv_con,
1454 abs_max, Lveil;
1455 char maskfile[500],file_out[500], file_out2[500], version[500];
1456 char *cline;
1457 VIEW userview = STDVIEW;
1458 int gotuserview = 0;
1459 struct muc_rvar* s_mask;
1460 s_mask = muc_rvar_create();
1461 muc_rvar_set_dim(s_mask, 1);
1462 muc_rvar_clear(s_mask);
1463 struct muc_rvar* s_band;
1464 s_band = muc_rvar_create();
1465 muc_rvar_set_dim(s_band, 1);
1466 muc_rvar_clear(s_band);
1467 struct muc_rvar* s_z1;
1468 s_z1 = muc_rvar_create();
1469 muc_rvar_set_dim(s_z1, 1);
1470 muc_rvar_clear(s_z1);
1471
1472 struct muc_rvar* s_z2;
1473 s_z2 = muc_rvar_create();
1474 muc_rvar_set_dim(s_z2, 1);
1475 muc_rvar_clear(s_z2);
1476
1477 struct muc_rvar* s_noposweight;
1478 s_noposweight = muc_rvar_create();
1479 muc_rvar_set_dim(s_noposweight, 1);
1480 muc_rvar_clear(s_noposweight);
1481
1482 struct muc_rvar* s_posweight;
1483 s_posweight = muc_rvar_create();
1484 muc_rvar_set_dim(s_posweight, 1);
1485 muc_rvar_clear(s_posweight);
1486
1487 struct muc_rvar* s_posweight2;
1488 s_posweight2 = muc_rvar_create();
1489 muc_rvar_set_dim(s_posweight2, 1);
1490 muc_rvar_clear(s_posweight2);
1491
1492 /*set required user view parameters to invalid values*/
1493 dir_ill=0.0;
1494 delta_E=0.0;
1495 no_glaresources=0;
1496 n_corner_px=0;
1497 zero_corner_px=0;
1498 force=0;
1499 dist=0.0;
1500 u_r=0.33333;
1501 u_g=0.33333;
1502 u_b=0.33333;
1503 band_angle=0;
1504 angle_z1=0;
1505 angle_z2=0;
1506 x_zone=0;
1507 y_zone=0;
1508 per_75_z2=0;
1509 per_95_z2=0;
1510 lum_pos_mean=0;
1511 lum_pos2_mean=0;
1512 lum_band_av = 0.0;
1513 omega_band = 0.0;
1514 pgsv = 0.0 ;
1515 pgsv_con = 0.0 ;
1516 pgsv_sat = 0.0 ;
1517 E_v_mask = 0.0;
1518 lum_z1_av = 0.0;
1519 omega_z1 = 0.0;
1520 lum_z2_av = 0.0;
1521 omega_z2 = 0.0 ;
1522 i_z1 = 0;
1523 i_z2 = 0;
1524 zones = 0;
1525 lum_z2_av = 0.0;
1526 uniform_gs = 0;
1527 apply_disability = 0;
1528 disability_thresh = 0;
1529 Lveil_cie_sum=0.0;
1530 skip_second_scan=0;
1531 lum_total_max=0.0;
1532 calcfast=0;
1533 age_corr_factor = 1.0;
1534 age_corr = 0;
1535 age = 20.0;
1536 userview.horiz = 0;
1537 userview.vert = 0;
1538 userview.type = 0;
1539 dgp_ext = 0;
1540 E_vl_ext = 0.0;
1541 new_lum_max = 0.0;
1542 lum_max = 0.0;
1543 omegat = 0.0;
1544 yt = 0;
1545 xt = 0;
1546 x_disk=0;
1547 y_disk=0;
1548 angle_disk=0;
1549 yfillmin = 0;
1550 yfillmax = 0;
1551 cut_view = 0;
1552 cut_view_type = 0;
1553 setvalue = 2e09;
1554 omega_cos_contr = 0.0;
1555 lum_ideal = 0.0;
1556 max_angle = 0.2;
1557 lum_thres = 2000.0;
1558 lum_task = 0.0;
1559 task_lum = 0;
1560 sgs = 0;
1561 splithigh = 1;
1562 detail_out = 0;
1563 detail_out2 = 0;
1564 posindex_picture = 0;
1565 checkfile = 0;
1566 ext_vill = 0;
1567 fill = 0;
1568 a1 = 2.0;
1569 a2 = 1.0;
1570 a3 = 1.87;
1571 a4 = 2.0;
1572 a5 = 1.0;
1573 c1 = 5.87e-05;
1574 c2 = 0.092;
1575 c3 = 0.159;
1576 non_cos_lb = 0;
1577 posindex_2 = 0;
1578 task_color = 0;
1579 limit = 50000.0;
1580 set_lum_max = 0;
1581 set_lum_max2 = 0;
1582 img_corr=0;
1583 abs_max = 0;
1584 progname = argv[0];
1585 E_v_contr = 0.0;
1586 strcpy(version, "1.19 release 09.12.2015 by J.Wienold");
1587 low_light_corr=1.0;
1588 output = 0;
1589 calc_vill = 0;
1590 band_avlum = -99;
1591 band_calc = 0;
1592 masking = 0;
1593 lum_mask[0]=0.0;
1594 lum_mask_av=0.0;
1595 omega_mask=0.0;
1596 i_mask=0;
1597 actual_igs=0;
1598 LUM_replace=0;
1599 thres_activate=0;
1600 /* command line for output picture*/
1601
1602 cline = (char *) malloc(CLINEMAX+1);
1603 cline[0] = '\0';
1604 for (i = 0; i < argc; i++) {
1605 /* fprintf(stderr, "%d %d \n",i,strlen(argv[i]));*/
1606 if (strlen(cline)+strlen(argv[i])+strlen(RELEASENAME)+2 >=CLINEMAX) {
1607 exit (-1);
1608 }
1609 else {
1610 strcat(cline, argv[i]);
1611 strcat(cline, " ");
1612 }
1613 }
1614 strcat(cline, "\n");
1615 strcat(cline, RELEASENAME);
1616 strcat(cline, "\n");
1617
1618
1619 /* program options */
1620
1621 for (i = 1; i < argc && argv[i][0] == '-'; i++) {
1622 /* expand arguments */
1623 while ((rval = expandarg(&argc, &argv, i)) > 0);
1624 if (rval < 0) {
1625 fprintf(stderr, "%s: cannot expand '%s'", argv[0], argv[i]);
1626 exit(1);
1627 }
1628 rval = getviewopt(&userview, argc - i, argv + i);
1629 if (rval >= 0) {
1630 i += rval;
1631 gotuserview++;
1632 continue;
1633 }
1634 switch (argv[i][1]) {
1635 case 'a':
1636 age = atof(argv[++i]);
1637 age_corr = 1;
1638 printf("age factor not supported any more \n");
1639 exit(1);
1640 break;
1641 case 'A':
1642 masking = 1;
1643 detail_out = 1;
1644 strcpy(maskfile, argv[++i]);
1645 pict_read(pm,maskfile );
1646
1647 break;
1648 case 'b':
1649 lum_thres = atof(argv[++i]);
1650 thres_activate = 1;
1651 break;
1652 case 'c':
1653 checkfile = 1;
1654 strcpy(file_out, argv[++i]);
1655 break;
1656 case 'u':
1657 uniform_gs = 1;
1658 u_r = atof(argv[++i]);
1659 u_g = atof(argv[++i]);
1660 u_b = atof(argv[++i]);
1661 break;
1662 case 'r':
1663 max_angle = atof(argv[++i]);
1664 break;
1665 case 's':
1666 sgs = 1;
1667 break;
1668 case 'f':
1669 force = 1;
1670 break;
1671 case 'k':
1672 apply_disability = 1;
1673 disability_thresh = atof(argv[++i]);
1674 break;
1675
1676 case 'p':
1677 posindex_picture = 1;
1678 break;
1679 case 'P':
1680 posindex_2 = atoi(argv[++i]);
1681 break;
1682
1683
1684 case 'y':
1685 splithigh = 1;
1686 break;
1687 case 'x':
1688 splithigh = 0;
1689 break;
1690 case 'Y':
1691 splithigh = 1;
1692 limit = atof(argv[++i]);
1693 break;
1694
1695 case 'i':
1696 ext_vill = 1;
1697 E_vl_ext = atof(argv[++i]);
1698 break;
1699 case 'I':
1700 ext_vill = 1;
1701 fill = 1;
1702 E_vl_ext = atof(argv[++i]);
1703 yfillmax = atoi(argv[++i]);
1704 yfillmin = atoi(argv[++i]);
1705 break;
1706 case 'd':
1707 detail_out = 1;
1708 break;
1709 case 'D':
1710 detail_out2 = 1;
1711 break;
1712 case 'm':
1713 img_corr = 1;
1714 set_lum_max = 1;
1715 lum_max = atof(argv[++i]);
1716 new_lum_max = atof(argv[++i]);
1717 strcpy(file_out2, argv[++i]);
1718 /* printf("max lum set to %f\n",new_lum_max);*/
1719 break;
1720
1721 case 'M':
1722 img_corr = 1;
1723 set_lum_max2 = 1;
1724 lum_max = atof(argv[++i]);
1725 E_vl_ext = atof(argv[++i]);
1726 strcpy(file_out2, argv[++i]);
1727 /* printf("max lum set to %f\n",new_lum_max);*/
1728 break;
1729 case 'N':
1730 img_corr = 1;
1731 set_lum_max2 = 2;
1732 x_disk = atoi(argv[++i]);
1733 y_disk = atoi(argv[++i]);
1734 angle_disk = atof(argv[++i]);
1735 E_vl_ext = atof(argv[++i]);
1736 strcpy(file_out2, argv[++i]);
1737 /* printf("max lum set to %f\n",new_lum_max);*/
1738 break;
1739 case 'O':
1740 img_corr = 1;
1741 set_lum_max2 = 3;
1742 x_disk = atoi(argv[++i]);
1743 y_disk = atoi(argv[++i]);
1744 angle_disk = atof(argv[++i]);
1745 LUM_replace = atof(argv[++i]);
1746 strcpy(file_out2, argv[++i]);
1747 /* printf("max lum set to %f\n",new_lum_max);*/
1748 break;
1749
1750
1751 /* deactivated case 'n':
1752 non_cos_lb = 0;
1753 break;
1754 */
1755 case 'q':
1756 non_cos_lb = atoi(argv[++i]);
1757 break;
1758
1759 case 't':
1760 task_lum = 1;
1761 xt = atoi(argv[++i]);
1762 yt = atoi(argv[++i]);
1763 omegat = atof(argv[++i]);
1764 task_color = 0;
1765 break;
1766 case 'T':
1767 task_lum = 1;
1768 xt = atoi(argv[++i]);
1769 yt = atoi(argv[++i]);
1770 /* omegat= DEG2RAD(atof(argv[++i]));*/
1771 omegat = atof(argv[++i]);
1772 task_color = 1;
1773 break;
1774 case 'l':
1775 zones = 1;
1776 x_zone = atoi(argv[++i]);
1777 y_zone = atoi(argv[++i]);
1778 angle_z1 = atof(argv[++i]);
1779 angle_z2 = -1;
1780 break;
1781 case 'L':
1782 zones = 2;
1783 x_zone = atoi(argv[++i]);
1784 y_zone = atoi(argv[++i]);
1785 angle_z1 = atof(argv[++i]);
1786 angle_z2 = atof(argv[++i]);
1787 break;
1788
1789
1790 case 'B':
1791 band_calc = 1;
1792 band_angle = atof(argv[++i]);
1793 band_color = 1;
1794 break;
1795
1796
1797 case 'w':
1798 a1 = atof(argv[++i]);
1799 a2 = atof(argv[++i]);
1800 a3 = atof(argv[++i]);
1801 a4 = atof(argv[++i]);
1802 a5 = atof(argv[++i]);
1803 c1 = atof(argv[++i]);
1804 c2 = atof(argv[++i]);
1805 c3 = atof(argv[++i]);
1806 break;
1807 case 'V':
1808 calc_vill = 1;
1809 break;
1810 case 'G':
1811 cut_view = 1;
1812 cut_view_type= atof(argv[++i]);
1813 break;
1814 case 'g':
1815 cut_view = 2;
1816 cut_view_type= atof(argv[++i]);
1817 break;
1818
1819 /*case 'v':
1820 printf("evalglare %s \n",version);
1821 exit(1); */
1822
1823 case '1':
1824 output = 1;
1825 break;
1826 case '2':
1827 output = 2;
1828 dir_ill = atof(argv[++i]);
1829 break;
1830
1831 case 'v':
1832 if (argv[i][2] == '\0') {
1833 printf("%s \n",RELEASENAME);
1834 exit(1);
1835 }
1836 if (argv[i][2] != 'f')
1837 goto userr;
1838 rval = viewfile(argv[++i], &userview, NULL);
1839 if (rval < 0) {
1840 fprintf(stderr,
1841 "%s: cannot open view file \"%s\"\n",
1842 progname, argv[i]);
1843 exit(1);
1844 } else if (rval == 0) {
1845 fprintf(stderr,
1846 "%s: bad view file \"%s\"\n", progname, argv[i]);
1847 exit(1);
1848 } else {
1849 gotuserview++;
1850 }
1851 break;
1852
1853
1854 default:
1855 goto userr;
1856 }
1857 }
1858
1859 /* set multiplier for task method to 5, if not specified */
1860
1861 if ( task_lum == 1 && thres_activate == 0){
1862 lum_thres = 5.0;
1863 }
1864 /*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/
1865
1866 if (output == 1 && ext_vill == 1 ) {
1867 calcfast=1;
1868 }
1869
1870 if (output == 2 && ext_vill == 1 ) {
1871 calcfast=2;
1872 }
1873
1874 /*masking and zoning cannot be applied at the same time*/
1875
1876 if (masking ==1 && zones >0) {
1877 fprintf(stderr, " masking and zoning cannot be activated at the same time!\n");
1878 exit(1);
1879 }
1880
1881 /* read picture file */
1882 if (i == argc) {
1883 SET_FILE_BINARY(stdin);
1884 FILE *fp = fdopen(fileno(stdin), "rb");
1885 if (!(fp)) {
1886 fprintf(stderr, "fdopen on stdin failed\n");
1887 return EXIT_FAILURE;
1888 }
1889 if (!(pict_read_fp(p, fp)))
1890 return EXIT_FAILURE;;
1891 } else {
1892 if (!pict_read(p, argv[i]))
1893 return EXIT_FAILURE;
1894 }
1895 if (gotuserview) {
1896 if (p->valid_view)
1897 fprintf(stderr,
1898 "warning: overriding image view by commandline argument\n");
1899 if ((userview.horiz == 0) || (userview.vert == 0) || (userview.type == 0)) {
1900 fprintf(stderr, "error: if specified, a view must at least contain -vt -vv and -vh\n");
1901 return EXIT_FAILURE;
1902 }
1903 p->view = userview;
1904 if (!(pict_update_view(p))) {
1905 fprintf(stderr, "error: invalid view specified");
1906 return EXIT_FAILURE;
1907 }
1908 pict_update_evalglare_caches(p);
1909 } else if (!(p->valid_view)) {
1910 fprintf(stderr, "error: no valid view specified\n");
1911 return EXIT_FAILURE;
1912 }
1913
1914
1915
1916 /* fprintf(stderr,"Picture read!");*/
1917
1918 /* command line for output picture*/
1919
1920 pict_set_comment(p, cline);
1921
1922 /* several checks */
1923 if (p->view.type == VT_PAR) {
1924 fprintf(stderr, "error: wrong view type! must not be parallel ! \n");
1925 exit(1);
1926 }
1927
1928 if (p->view.type == VT_PER) {
1929 fprintf(stderr, "warning: image has perspective view type specified in header ! \n");
1930 }
1931
1932 if (masking == 1) {
1933
1934 if (!pict_get_xsize(p)==pict_get_xsize(pm) || !pict_get_ysize(p)==pict_get_ysize(pm)) {
1935 fprintf(stderr, "error: masking image has other resolution than main image ! \n");
1936 fprintf(stderr, "size must be identical \n");
1937 printf("resolution main image : %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
1938 printf("resolution masking image : %dx%d\n",pict_get_xsize(pm),pict_get_ysize(pm));
1939 exit(1);
1940
1941 }
1942
1943 }
1944 /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1945
1946 /* Check size of search radius */
1947 rmx = (pict_get_xsize(p) / 2);
1948 rmy = (pict_get_ysize(p) / 2);
1949 rx = (pict_get_xsize(p) / 2 + 10);
1950 ry = (pict_get_ysize(p) / 2 + 10);
1951 r_center =
1952 acos(DOT(pict_get_cached_dir(p, rmx, rmy),
1953 pict_get_cached_dir(p, rx, ry))) * 2 / 10;
1954 search_pix = max_angle / r_center;
1955 if (search_pix < 1.0) {
1956 fprintf(stderr,
1957 "warning: search radius less than 1 pixel! deactivating smoothing and peak extraction...\n");
1958 splithigh = 0;
1959 sgs = 0;
1960
1961 } else {
1962 if (search_pix < 3.0) {
1963 fprintf(stderr,
1964 "warning: search radius less than 3 pixels! -> %f \n",
1965 search_pix);
1966
1967 }
1968 }
1969
1970 /* Check task position */
1971
1972 if (task_lum == 1) {
1973 if (xt >= pict_get_xsize(p) || yt >= pict_get_ysize(p) || xt < 0
1974 || yt < 0) {
1975 fprintf(stderr,
1976 "error: task position outside picture!! exit...");
1977 exit(1);
1978 }
1979 }
1980
1981
1982 /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1983
1984 sang = 0.0;
1985 E_v = 0.0;
1986 E_v_dir = 0.0;
1987 avlum = 0.0;
1988 pict_new_gli(p);
1989 igs = 0;
1990 pict_get_z1_gsn(p,igs) = 0;
1991 pict_get_z2_gsn(p,igs) = 0;
1992
1993
1994 /* cut out GUTH field of view and exit without glare evaluation */
1995 if (cut_view==2) {
1996 if (cut_view_type==1) {
1997 cut_view_1(p);
1998 }else{
1999
2000 if (cut_view_type==2) {
2001 cut_view_2(p);
2002 }else{
2003 if (cut_view_type==3) {
2004 fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
2005 cut_view_3(p);
2006 }else{
2007 fprintf(stderr,"error: no valid option for view cutting!!");
2008 exit(1);
2009 }
2010 }}
2011 pict_write(p, file_out);
2012 exit(1); }
2013
2014
2015
2016
2017
2018
2019 /* write positionindex into checkfile and exit, activated by option -p */
2020 if (posindex_picture == 1) {
2021 for (x = 0; x < pict_get_xsize(p); x++)
2022 for (y = 0; y < pict_get_ysize(p); y++) {
2023 if (pict_get_hangle
2024 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2025 if (pict_is_validpixel(p, x, y)) {
2026 lum =
2027 get_posindex(p, x, y,
2028 posindex_2) / WHTEFFICACY;
2029
2030 pict_get_color(p, x, y)[RED] = lum;
2031 pict_get_color(p, x, y)[GRN] = lum;
2032 pict_get_color(p, x, y)[BLU] = lum;
2033 lum_task = luminance(pict_get_color(p, x, y));
2034 /*printf("x,y,posindex=%i %i %f %f\n",x,y,lum*WHTEFFICACY,lum_task);*/
2035 }
2036 }
2037 }
2038 pict_write(p, file_out);
2039 exit(1);
2040 }
2041
2042
2043 /* calculation of vertical illuminance, average luminance, in case of filling activated-> fill picture */
2044 /* fill, if necessary from 0 to yfillmin */
2045
2046 if (fill == 1) {
2047 for (x = 0; x < pict_get_xsize(p); x++)
2048 for (y = yfillmin; y > 0; y = y - 1) {
2049 y1 = y + 1;
2050 lum = luminance(pict_get_color(p, x, y1));
2051 pict_get_color(p, x, y)[RED] = lum / 179.0;
2052 pict_get_color(p, x, y)[GRN] = lum / 179.0;
2053 pict_get_color(p, x, y)[BLU] = lum / 179.0;
2054 }
2055 }
2056
2057 if (calcfast == 0) {
2058 for (x = 0; x < pict_get_xsize(p); x++)
2059 for (y = 0; y < pict_get_ysize(p); y++) {
2060 lum = luminance(pict_get_color(p, x, y));
2061 dist=sqrt((x-rmx)*(x-rmx)+(y-rmy)*(y-rmy));
2062 if (dist>((rmx+rmy)/2)) {
2063 n_corner_px=n_corner_px+1;
2064 if (lum<7.0) {
2065 zero_corner_px=zero_corner_px+1;
2066 }
2067 }
2068
2069
2070 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2071 if (pict_is_validpixel(p, x, y)) {
2072 if (fill == 1 && y >= yfillmax) {
2073 y1 = y - 1;
2074 lum = luminance(pict_get_color(p, x, y1));
2075 pict_get_color(p, x, y)[RED] = lum / 179.0;
2076 pict_get_color(p, x, y)[GRN] = lum / 179.0;
2077 pict_get_color(p, x, y)[BLU] = lum / 179.0;
2078 }
2079
2080 if (lum > abs_max) {
2081 abs_max = lum;
2082 }
2083 /* set luminance restriction, if -m is set */
2084 if (img_corr == 1 ) {
2085 if (set_lum_max == 1 && lum >= lum_max) {
2086 pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
2087 pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
2088 pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
2089 /* printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
2090 lum = luminance(pict_get_color(p, x, y));
2091 }
2092 if (set_lum_max2 == 1 && lum >= lum_max) {
2093
2094 E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2095 omega_cos_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2096 }
2097 if (set_lum_max2 == 2 ) {
2098 r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2099 if (x_disk == x && y_disk==y ) r_actual=0.0;
2100
2101 act_lum = luminance(pict_get_color(p, x, y));
2102
2103 if (r_actual <= angle_disk) {
2104 E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2105 omega_cos_contr += DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2106 }
2107
2108
2109
2110 }
2111 }
2112 om = pict_get_omega(p, x, y);
2113 sang += om;
2114 E_v += DOT(p->view.vdir, pict_get_cached_dir(p, x,y)) * om *lum;
2115 avlum += om * lum;
2116 pos=get_posindex(p, x, y,posindex_2);
2117
2118 lum_pos[0]=lum;
2119 muc_rvar_add_sample(s_noposweight,lum_pos);
2120 lum_pos[0]=lum/pos;
2121 lum_pos_mean+=lum_pos[0]*om;
2122 muc_rvar_add_sample(s_posweight, lum_pos);
2123 lum_pos[0]=lum_pos[0]/pos;
2124 lum_pos2_mean+=lum_pos[0]*om;
2125 muc_rvar_add_sample(s_posweight2,lum_pos);
2126
2127
2128
2129 } else {
2130 pict_get_color(p, x, y)[RED] = 0.0;
2131 pict_get_color(p, x, y)[GRN] = 0.0;
2132 pict_get_color(p, x, y)[BLU] = 0.0;
2133
2134 }
2135 }else {
2136 pict_get_color(p, x, y)[RED] = 0.0;
2137 pict_get_color(p, x, y)[GRN] = 0.0;
2138 pict_get_color(p, x, y)[BLU] = 0.0;
2139
2140 }
2141 }
2142
2143 /* check if image has black corners AND a perspective view */
2144
2145 if (zero_corner_px/n_corner_px > 0.70 && (p->view.type == VT_PER) ) {
2146 printf (" corner pixels are to %f %% black! \n",zero_corner_px/n_corner_px*100);
2147 printf("error! The image has black corners AND header contains a perspective view type definition !!!\n") ;
2148
2149 if (force==0) {
2150 printf("stopping...!!!!\n") ;
2151
2152 exit(1);
2153 }
2154 }
2155
2156
2157
2158
2159 lum_pos_mean= lum_pos_mean/sang;
2160 lum_pos2_mean= lum_pos2_mean/sang;
2161
2162 if ((set_lum_max2 >= 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0 ) || set_lum_max2==3) {
2163
2164 if (set_lum_max2<3){
2165 lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
2166 if (set_lum_max2 == 2 && lum_ideal >= 2e9) {
2167 printf("warning! luminance of replacement pixels would be larger than 2e9 cd/m2. Value set to 2e9cd/m2!\n") ;
2168 }
2169
2170 if (lum_ideal > 0 && lum_ideal < setvalue) {
2171 setvalue = lum_ideal;
2172 }
2173 printf("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f %f %f %f %f\n",
2174 lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
2175 }else{setvalue=LUM_replace;
2176 }
2177
2178
2179 for (x = 0; x < pict_get_xsize(p); x++)
2180 for (y = 0; y < pict_get_ysize(p); y++) {
2181 if (pict_get_hangle
2182 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2183 if (pict_is_validpixel(p, x, y)) {
2184 lum = luminance(pict_get_color(p, x, y));
2185
2186
2187 if (set_lum_max2 == 1 && lum >= lum_max) {
2188
2189 pict_get_color(p, x, y)[RED] =
2190 setvalue / 179.0;
2191 pict_get_color(p, x, y)[GRN] =
2192 setvalue / 179.0;
2193 pict_get_color(p, x, y)[BLU] =
2194 setvalue / 179.0;
2195
2196 }else{ if(set_lum_max2 >1 ) {
2197 r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2198 if (x_disk == x && y_disk==y ) r_actual=0.0;
2199
2200 if (r_actual <= angle_disk) {
2201
2202 pict_get_color(p, x, y)[RED] =
2203 setvalue / 179.0;
2204 pict_get_color(p, x, y)[GRN] =
2205 setvalue / 179.0;
2206 pict_get_color(p, x, y)[BLU] =
2207 setvalue / 179.0;
2208
2209 }
2210 }
2211 }
2212 }
2213 }
2214 }
2215
2216
2217 pict_write(p, file_out2);
2218 exit(1);
2219 }
2220 if (set_lum_max == 1) {
2221 pict_write(p, file_out2);
2222
2223 }
2224
2225 if (calc_vill == 1) {
2226 printf("%f\n", E_v);
2227 exit(1);
2228 }
2229 }else{
2230 /* in fast calculation mode: ev=ev_ext and sang=2*pi */
2231 sang=2*3.14159265359;
2232 lum_task =E_vl_ext/sang;
2233 E_v=E_vl_ext;
2234 /* printf("calc fast!! %f %f\n", sang,lum_task);*/
2235
2236
2237 }
2238
2239 /* cut out GUTH field of view for glare evaluation */
2240 if (cut_view==1) {
2241 if (cut_view_type==1) {
2242 cut_view_1(p);
2243 }else{
2244
2245 if (cut_view_type==2) {
2246 cut_view_2(p);
2247 }else{
2248 if (cut_view_type==3) {
2249 fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
2250 cut_view_3(p);
2251 }else{
2252 fprintf(stderr,"error: no valid option for view cutting!!");
2253 exit(1);
2254 }
2255 }}
2256 }
2257
2258 if (calcfast == 0) {
2259 avlum = avlum / sang;
2260 lum_task = avlum;
2261 }
2262 /* if (ext_vill==1) {
2263 E_v=E_vl_ext;
2264 avlum=E_v/3.1415927;
2265 } */
2266
2267 if (task_lum == 1) {
2268 lum_task = get_task_lum(p, xt, yt, omegat, task_color);
2269 }
2270 lum_source = lum_thres * lum_task;
2271 /* printf("Task Luminance=%f\n",lum_task);
2272 pict_write(p,"task.pic");*/
2273
2274 if (lum_thres > 100) {
2275 lum_source = lum_thres;
2276 }
2277
2278 /* printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2279
2280 /* first glare source scan: find primary glare sources */
2281 for (x = 0; x < pict_get_xsize(p); x++) {
2282 /* lastpixelwas_gs=0; */
2283 for (y = 0; y < pict_get_ysize(p); y++) {
2284 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2285 if (pict_is_validpixel(p, x, y)) {
2286 act_lum = luminance(pict_get_color(p, x, y));
2287 if (act_lum > lum_source) {
2288 if (act_lum >lum_total_max) {
2289 lum_total_max=act_lum;
2290 }
2291 /* speed improvement first scan: when last pixel was glare source, then use this glare source number instead of searching.
2292 has been deactivated with v1.25
2293
2294 if (lastpixelwas_gs==0 || search_pix <= 1.0 ) { */
2295 actual_igs = find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2296 /* }*/
2297 if (actual_igs == 0) {
2298 igs = igs + 1;
2299 pict_new_gli(p);
2300 pict_get_lum_min(p, igs) = HUGE_VAL;
2301 pict_get_Eglare(p,igs) = 0.0;
2302 pict_get_Dglare(p,igs) = 0.0;
2303 pict_get_z1_gsn(p,igs) = 0;
2304 pict_get_z2_gsn(p,igs) = 0;
2305 actual_igs = igs;
2306
2307 }
2308 /* no coloring any more here icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2309 pict_get_gsn(p, x, y) = actual_igs;
2310 pict_get_pgs(p, x, y) = 1;
2311 add_pixel_to_gs(p, x, y, actual_igs);
2312 /* lastpixelwas_gs=actual_igs; */
2313
2314 } else {
2315 pict_get_pgs(p, x, y) = 0;
2316 /* lastpixelwas_gs=0;*/
2317 }
2318 }
2319 }
2320 }
2321 }
2322 /* pict_write(p,"firstscan.pic"); */
2323
2324
2325
2326
2327 if (calcfast ==1 || search_pix <= 1.0 || calcfast == 2 ) {
2328 skip_second_scan=1;
2329 }
2330
2331
2332 /* second glare source scan: combine glare sources facing each other */
2333 change = 1;
2334 i = 0;
2335 while (change == 1 && skip_second_scan==0) {
2336 change = 0;
2337 i = i+1;
2338 for (x = 0; x < pict_get_xsize(p); x++) {
2339 for (y = 0; y < pict_get_ysize(p); y++) {
2340 if (pict_get_hangle
2341 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2342 checkpixels=1;
2343 before_igs = pict_get_gsn(p, x, y);
2344
2345 /* speed improvement: search for other nearby glare sources only, when at least one adjacend pixel has another glare source number */
2346 if (x >1 && x<pict_get_xsize(p)-2 && y >1 && y<pict_get_ysize(p)-2 ) {
2347 x1=x-1;
2348 x2=x+1;
2349 y1=y-1;
2350 y2=y+1;
2351 if (before_igs == pict_get_gsn(p, x1, y) && before_igs == pict_get_gsn(p, x2, y) && before_igs == pict_get_gsn(p, x, y1) && before_igs == pict_get_gsn(p, x, y2) && before_igs == pict_get_gsn(p, x1, y1) && before_igs == pict_get_gsn(p, x2, y1) && before_igs == pict_get_gsn(p, x1, y2) && before_igs == pict_get_gsn(p, x2, y2) ) {
2352 checkpixels = 0;
2353 actual_igs = before_igs;
2354 } }
2355
2356
2357 if (pict_is_validpixel(p, x, y) && before_igs > 0 && checkpixels==1) {
2358 actual_igs =
2359 find_near_pgs(p, x, y, max_angle, before_igs,
2360 igs, 1);
2361 if (!(actual_igs == before_igs)) {
2362 change = 1;
2363 }
2364 if (before_igs > 0) {
2365 actual_igs = pict_get_gsn(p, x, y);
2366 /* icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2367 }
2368
2369 }
2370 }
2371 }
2372 /* pict_write(p,"secondscan.pic");*/
2373 }}
2374
2375 /*smoothing: add secondary glare sources */
2376 i_max = igs;
2377 if (sgs == 1) {
2378
2379 /* simplified search radius: constant for entire picture, alway circle not an angle due to performance */
2380 x = (pict_get_xsize(p) / 2);
2381 y = (pict_get_ysize(p) / 2);
2382 rx = (pict_get_xsize(p) / 2 + 10);
2383 ry = (pict_get_ysize(p) / 2 + 10);
2384 r_center =
2385 acos(DOT(pict_get_cached_dir(p, x, y),
2386 pict_get_cached_dir(p, rx, ry))) * 2 / 10;
2387 search_pix = max_angle / r_center / 1.75;
2388
2389 for (x = 0; x < pict_get_xsize(p); x++) {
2390 for (y = 0; y < pict_get_ysize(p); y++) {
2391 if (pict_get_hangle
2392 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2393 if (pict_is_validpixel(p, x, y)
2394 && pict_get_gsn(p, x, y) == 0) {
2395 for (i = 1; i <= i_max; i++) {
2396
2397 if (pict_get_npix(p, i) > 0) {
2398 add_secondary_gs(p, x, y, search_pix, i, uniform_gs, u_r, u_g , u_b);
2399 }
2400 }
2401
2402 }
2403 }
2404 }
2405 }
2406
2407 }
2408
2409 /* extract extremes from glare sources to extra glare source */
2410 if (splithigh == 1 && lum_total_max>limit) {
2411
2412 r_split = max_angle / 2.0;
2413 for (i = 0; i <= i_max; i++) {
2414
2415 if (pict_get_npix(p, i) > 0) {
2416 l_max = pict_get_lum_max(p, i);
2417 i_splitstart = igs + 1;
2418 if (l_max >= limit) {
2419 for (x = 0; x < pict_get_xsize(p); x++)
2420 for (y = 0; y < pict_get_ysize(p); y++) {
2421 if (pict_get_hangle
2422 (p, x, y, p->view.vdir, p->view.vup, &ang))
2423 {
2424 if (pict_is_validpixel(p, x, y)
2425 && luminance(pict_get_color(p, x, y))
2426 >= limit
2427 && pict_get_gsn(p, x, y) == i) {
2428 if (i_splitstart == (igs + 1)) {
2429 igs = igs + 1;
2430 pict_new_gli(p);
2431 pict_get_z1_gsn(p,igs) = 0;
2432 pict_get_z2_gsn(p,igs) = 0;
2433
2434 pict_get_Eglare(p,igs) = 0.0;
2435 pict_get_Dglare(p,igs) = 0.0;
2436 pict_get_lum_min(p, igs) =
2437 99999999999999.999;
2438 i_split = igs;
2439 } else {
2440 i_split =
2441 find_split(p, x, y, r_split,
2442 i_splitstart, igs);
2443 }
2444 if (i_split == 0) {
2445 igs = igs + 1;
2446 pict_new_gli(p);
2447 pict_get_z1_gsn(p,igs) = 0;
2448 pict_get_z2_gsn(p,igs) = 0;
2449
2450 pict_get_Eglare(p,igs) = 0.0;
2451 pict_get_Dglare(p,igs) = 0.0;
2452 pict_get_lum_min(p, igs) =
2453 99999999999999.999;
2454 i_split = igs;
2455 }
2456 split_pixel_from_gs(p, x, y, i_split, uniform_gs, u_r, u_g , u_b);
2457
2458 }
2459 }
2460 }
2461
2462 }
2463 change = 1;
2464 while (change == 1) {
2465 change = 0;
2466 for (x = 0; x < pict_get_xsize(p); x++)
2467 for (y = 0; y < pict_get_ysize(p); y++) {
2468 before_igs = pict_get_gsn(p, x, y);
2469 if (before_igs >= i_splitstart) {
2470 if (pict_get_hangle
2471 (p, x, y, p->view.vdir, p->view.vup,
2472 &ang)) {
2473 if (pict_is_validpixel(p, x, y)
2474 && before_igs > 0) {
2475 actual_igs =
2476 find_near_pgs(p, x, y,
2477 max_angle,
2478 before_igs, igs,
2479 i_splitstart);
2480 if (!(actual_igs == before_igs)) {
2481 change = 1;
2482 }
2483 if (before_igs > 0) {
2484 actual_igs =
2485 pict_get_gsn(p, x, y);
2486 /* icol =
2487 setglcolor(p, x, y,
2488 actual_igs, uniform_gs, u_r, u_g , u_b);*/
2489 }
2490
2491 }
2492 }
2493 }
2494
2495 }
2496 }
2497
2498
2499 }
2500 }
2501 }
2502
2503 /* calculation of direct vertical illuminance for CGI and for disability glare, coloring glare sources*/
2504
2505 if (calcfast == 0 || calcfast == 2) {
2506 for (x = 0; x < pict_get_xsize(p); x++)
2507 for (y = 0; y < pict_get_ysize(p); y++) {
2508 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2509 if (pict_is_validpixel(p, x, y)) {
2510 if (pict_get_gsn(p, x, y) > 0) {
2511 actual_igs = pict_get_gsn(p, x, y);
2512 delta_E=DOT(p->view.vdir, pict_get_cached_dir(p, x, y))* pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2513 pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
2514 E_v_dir = E_v_dir + delta_E;
2515 setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2516 }
2517 }
2518 }
2519 }
2520 lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2521
2522 }
2523 /* calc of band luminance distribution if applied */
2524 if (band_calc == 1) {
2525 x_max = pict_get_xsize(p) - 1;
2526 y_max = pict_get_ysize(p) - 1;
2527 y_mid = (int)(y_max/2);
2528 for (x = 0; x <= x_max; x++)
2529 for (y = 0; y <= y_max; y++) {
2530 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2531 if (pict_is_validpixel(p, x, y)) {
2532
2533 r_actual = acos(DOT(pict_get_cached_dir(p, x, y_mid), pict_get_cached_dir(p, x, y))) ;
2534 act_lum = luminance(pict_get_color(p, x, y));
2535
2536 if ((r_actual <= band_angle) || (y == y_mid) ) {
2537 lum_band[0]=luminance(pict_get_color(p, x, y));
2538 muc_rvar_add_sample(s_band, lum_band);
2539 act_lum = luminance(pict_get_color(p, x, y));
2540 lum_band_av += pict_get_omega(p, x, y) * act_lum;
2541 omega_band += pict_get_omega(p, x, y);
2542 if (band_color == 1) {
2543 pict_get_color(p, x, y)[RED] = 0.0;
2544 pict_get_color(p, x, y)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
2545 pict_get_color(p, x, y)[BLU] = 0.0;
2546 }
2547 }
2548 }}}
2549 lum_band_av=lum_band_av/omega_band;
2550 muc_rvar_get_vx(s_band,lum_band_std);
2551 muc_rvar_get_percentile(s_band,lum_band_median,0.75);
2552 per_75_band=lum_band_median[0];
2553 muc_rvar_get_percentile(s_band,lum_band_median,0.95);
2554 per_95_band=lum_band_median[0];
2555 muc_rvar_get_median(s_band,lum_band_median);
2556 muc_rvar_get_bounding_box(s_band,bbox_band);
2557
2558 printf ("band:band_omega,band_av_lum,band_median_lum,band_std_lum,band_perc_75,band_perc_95,band_lum_min,band_lum_max: %f %f %f %f %f %f %f %f\n",omega_band,lum_band_av,lum_band_median[0],sqrt(lum_band_std[0]),per_75_band,per_95_band,bbox_band[0],bbox_band[1] );
2559 /* av_lum = av_lum / omega_sum; */
2560 /* printf("average luminance of band %f \n",av_lum);*/
2561 /* printf("solid angle of band %f \n",omega_sum);*/
2562 }
2563
2564 /*printf("total number of glare sources: %i\n",igs); */
2565 lum_sources = 0;
2566 omega_sources = 0;
2567 i=0;
2568 for (x = 0; x <= igs; x++) {
2569 if (pict_get_npix(p, x) > 0) {
2570 lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2571 omega_sources += pict_get_av_omega(p, x);
2572 i=i+1;
2573 }}
2574
2575 if (sang == omega_sources) {
2576 lum_backg =avlum;
2577 } else {
2578 lum_backg =(sang * avlum - lum_sources) / (sang - omega_sources);
2579 }
2580
2581
2582 if (i == 0) {
2583 lum_sources =0.0;
2584 } else { lum_sources=lum_sources/omega_sources;
2585 }
2586
2587
2588
2589 if (non_cos_lb == 0) {
2590 lum_backg = lum_backg_cos;
2591 }
2592
2593 if (non_cos_lb == 2) {
2594 lum_backg = E_v / 3.1415927;
2595 }
2596
2597
2598 /* file writing NOT here
2599 if (checkfile == 1) {
2600 pict_write(p, file_out);
2601 }
2602 */
2603
2604 /* print detailed output */
2605
2606
2607
2608 /* masking */
2609
2610 if ( masking == 1 ) {
2611
2612
2613 for (x = 0; x < pict_get_xsize(p); x++)
2614 for (y = 0; y < pict_get_ysize(p); y++) {
2615 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2616 if (pict_is_validpixel(p, x, y)) {
2617 if (luminance(pict_get_color(pm, x, y))>0.1) {
2618 /* printf ("hallo %i %i %f \n",x,y,luminance(pict_get_color(pm, x, y)));*/
2619 i_mask=i_mask+1;
2620 lum_mask[0]=luminance(pict_get_color(p, x, y));
2621 /* printf ("%f \n",lum_mask[0]);*/
2622 muc_rvar_add_sample(s_mask, lum_mask);
2623 omega_mask += pict_get_omega(p, x, y);
2624 lum_mask_av += pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2625 E_v_mask +=DOT(p->view.vdir, pict_get_cached_dir(p, x, y))*pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2626
2627 pict_get_color(pm, x, y)[RED] = luminance(pict_get_color(p, x, y))/179.0;
2628 pict_get_color(pm, x, y)[GRN] = luminance(pict_get_color(p, x, y))/179.0;
2629 pict_get_color(pm, x, y)[BLU] = luminance(pict_get_color(p, x, y))/179.0;
2630
2631 } else {
2632 pict_get_color(p, x, y)[RED] = 0.0 ;
2633 pict_get_color(p, x, y)[GRN] = 0.0 ;
2634 pict_get_color(p, x, y)[BLU] = 0.0 ;
2635 }
2636 }
2637 }
2638 }
2639 /* Average luminance in masking area (weighted by solid angle) */
2640 lum_mask_av=lum_mask_av/omega_mask;
2641 muc_rvar_get_vx(s_mask,lum_mask_std);
2642 muc_rvar_get_percentile(s_mask,lum_mask_median,0.75);
2643 per_75_mask=lum_mask_median[0];
2644 muc_rvar_get_percentile(s_mask,lum_mask_median,0.95);
2645 per_95_mask=lum_mask_median[0];
2646 muc_rvar_get_median(s_mask,lum_mask_median);
2647 muc_rvar_get_bounding_box(s_mask,bbox);
2648 /* PSGV only why masking of window is applied! */
2649
2650
2651 if (task_lum == 0 || lum_task == 0.0 ) {
2652 fprintf(stderr, " warning: Task area not set or task luminance=0 ! pgsv cannot be calculated (set to -99)!!\n");
2653 pgsv = -99;
2654 } else {
2655 pgsv = get_pgsv(E_v, E_v_mask, omega_mask, lum_mask_av,lum_task,avlum);
2656 }
2657
2658 pgsv_con = get_pgsv_con(E_v, E_v_mask, omega_mask, lum_mask_av,avlum);
2659 pgsv_sat =get_pgsv_sat(E_v);
2660
2661 if (detail_out == 1) {
2662
2663 printf ("masking:no_pixels,omega,av_lum,median_lum,std_lum,perc_75,perc_95,lum_min,lum_max,pgsv_con,pgsv_sat,pgsv,Ev_mask: %i %f %f %f %f %f %f %f %f %f %f %f %f\n",i_mask,omega_mask,lum_mask_av,lum_mask_median[0],sqrt(lum_mask_std[0]),per_75_mask,per_95_mask,bbox[0],bbox[1],pgsv_con,pgsv_sat,pgsv,E_v_mask );
2664
2665 }
2666
2667 }
2668
2669 /* zones */
2670
2671 if ( zones > 0 ) {
2672
2673 for (x = 0; x < pict_get_xsize(p); x++)
2674 for (y = 0; y < pict_get_ysize(p); y++) {
2675 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2676 if (pict_is_validpixel(p, x, y)) {
2677
2678
2679 r_actual =acos(DOT(pict_get_cached_dir(p, x, y), pict_get_cached_dir(p, x_zone, y_zone))) * 2;
2680 lum_actual = luminance(pict_get_color(p, x, y));
2681 act_gsn=pict_get_gsn(p, x, y);
2682 /* printf("1act_gsn,pict_get_z1_gsn,pict_get_z2_gsn %i %f %f \n",act_gsn,pict_get_z1_gsn(p,act_gsn),pict_get_z2_gsn(p,act_gsn));*/
2683
2684 /*zone1 */
2685 if (r_actual <= angle_z1) {
2686 i_z1=i_z1+1;
2687 lum_z1[0]=lum_actual;
2688 muc_rvar_add_sample(s_z1, lum_z1);
2689 omega_z1 += pict_get_omega(p, x, y);
2690 lum_z1_av += pict_get_omega(p, x, y)* lum_actual;
2691 setglcolor(p,x,y,1,1 , 0.66, 0.01 ,0.33);
2692 /*check if separation gsn already exist */
2693
2694 if (act_gsn > 0){
2695 if (pict_get_z1_gsn(p,act_gsn) == 0) {
2696 pict_new_gli(p);
2697 igs = igs + 1;
2698 pict_get_z1_gsn(p,act_gsn) = igs*1.000;
2699 pict_get_z1_gsn(p,igs) = -1.0;
2700 pict_get_z2_gsn(p,igs) = -1.0;
2701 splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2702 /* printf ("Glare source %i gets glare source %i in zone 1 : %i %f %f \n",act_gsn,igs,splitgs,pict_get_z1_gsn(p,act_gsn),pict_get_z1_gsn(p,igs)); */
2703 }
2704 splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2705 /* printf("splitgs%i \n",splitgs); */
2706 split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
2707 }
2708 }
2709 /*zone2 */
2710
2711 if (r_actual > angle_z1 && r_actual<= angle_z2 ) {
2712
2713 i_z2=i_z2+1;
2714 lum_z2[0]=lum_actual;
2715 muc_rvar_add_sample(s_z2, lum_z2);
2716 omega_z2 += pict_get_omega(p, x, y);
2717 lum_z2_av += pict_get_omega(p, x, y)* lum_actual;
2718 setglcolor(p,x,y,1,1 , 0.65, 0.33 ,0.02);
2719 /* printf("zone 2 x y act_gsn pict_get_z1_gsn(p,act_gsn) pict_get_z2_gsn(p,act_gsn): %i %i %f 1:%f 2:%f \n",x,y,act_gsn,pict_get_z1_gsn(p,act_gsn),pict_get_z2_gsn(p,act_gsn));
2720 */ if (act_gsn > 0){
2721 if (pict_get_z2_gsn(p,act_gsn) == 0) {
2722 pict_new_gli(p);
2723 igs = igs + 1;
2724 pict_get_z2_gsn(p,act_gsn) = igs*1.000;
2725 pict_get_z1_gsn(p,igs) = -2.0;
2726 pict_get_z2_gsn(p,igs) = -2.0;
2727 /* printf ("Glare source %i gets glare source %i in zone 2 \n",act_gsn,igs ); */
2728 }
2729 splitgs=(int)(pict_get_z2_gsn(p,act_gsn));
2730 /* printf("splitgs %i \n",splitgs);*/
2731 split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
2732 }
2733 }
2734
2735 }
2736 } }
2737 /* Average luminance in zones (weighted by solid angle), calculate percentiles, median min and max in zones */
2738 if (zones == 2 ) {
2739 lum_z2_av=lum_z2_av/omega_z2;
2740 muc_rvar_get_vx(s_z2,lum_z2_std);
2741 muc_rvar_get_percentile(s_z2,lum_z2_median,0.75);
2742 per_75_z2=lum_z2_median[0];
2743 muc_rvar_get_percentile(s_z2,lum_z2_median,0.95);
2744 per_95_z2=lum_z2_median[0];
2745 muc_rvar_get_median(s_z2,lum_z2_median);
2746 muc_rvar_get_bounding_box(s_z2,bbox_z2);
2747 }
2748 lum_z1_av=lum_z1_av/omega_z1;
2749 muc_rvar_get_vx(s_z1,lum_z1_std);
2750 muc_rvar_get_percentile(s_z1,lum_z1_median,0.75);
2751 per_75_z1=lum_z1_median[0];
2752 muc_rvar_get_percentile(s_z1,lum_z1_median,0.95);
2753 per_95_z1=lum_z1_median[0];
2754 muc_rvar_get_median(s_z1,lum_z1_median);
2755 muc_rvar_get_bounding_box(s_z1,bbox_z1);
2756 if (detail_out == 1) {
2757
2758 printf ("zoning:z1_omega,z1_av_lum,z1_median_lum,z1_std_lum,z1_perc_75,z1_perc_95,z1_lum_min,z1_lum_max: %f %f %f %f %f %f %f %f\n",omega_z1,lum_z1_av,lum_z1_median[0],sqrt(lum_z1_std[0]),per_75_z1,per_95_z1,bbox_z1[0],bbox_z1[1] );
2759
2760 if (zones == 2 ) {
2761
2762 printf ("zoning:z2_omega,z2_av_lum,z2_median_lum,z2_std_lum,z2_perc_75,z2_perc_95,z2_lum_min,z2_lum_max: %f %f %f %f %f %f %f %f\n",omega_z2,lum_z2_av,lum_z2_median[0],sqrt(lum_z2_std[0]),per_75_z2,per_95_z2,bbox_z2[0],bbox_z2[1] );
2763 } }
2764
2765 }
2766
2767 i = 0;
2768 for (x = 0; x <= igs; x++) {
2769 /* resorting glare source numbers */
2770 if (pict_get_npix(p, x) > 0) {
2771 i = i + 1;
2772 }
2773 }
2774 no_glaresources=i;
2775
2776 /* glare sources */
2777 if (detail_out == 1) {
2778
2779 printf
2780 ("%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 teta glare_zone\n",
2781 i);
2782 if (i == 0) {
2783 printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", i, 0.0, 0.0,
2784 0.0, 0.0, 0.0, 0.0, lum_backg, lum_task, E_v, E_v_dir,abs_max,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
2785
2786 } else {
2787 i = 0;
2788 for (x = 0; x <= igs; x++) {
2789 if (pict_get_npix(p, x) > 0) {
2790 i = i + 1;
2791 pict_get_sigma(p, pict_get_av_posx(p, x),
2792 pict_get_av_posy(p, x), p->view.vdir,
2793 p->view.vup, &sigma);
2794
2795 x2=pict_get_av_posx(p, x);
2796 y2=pict_get_av_posy(p, x);
2797 teta = 180.0 / 3.1415927 * acos(DOT(p->view.vdir, pict_get_cached_dir(p, x2, y2)));
2798 Lveil_cie = 10*pict_get_Eglare(p,x)/(teta*teta+0.0000000000000001);
2799
2800 if (apply_disability == 1 && Lveil_cie <=disability_thresh) {
2801 Lveil_cie =0 ;
2802 }
2803 Lveil_cie_sum = Lveil_cie_sum + Lveil_cie;
2804 if (pict_get_z1_gsn(p,x)<0) {
2805 act_gsn=(int)(-pict_get_z1_gsn(p,x));
2806 }else{
2807 act_gsn=0;
2808 }
2809 printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %i \n",
2810 i, pict_get_npix(p, x), pict_get_av_posx(p, x),
2811 pict_get_ysize(p) - pict_get_av_posy(p, x),
2812 pict_get_av_lum(p, x), pict_get_av_omega(p, x),
2813 get_posindex(p, pict_get_av_posx(p, x),
2814 pict_get_av_posy(p, x),
2815 posindex_2), lum_backg, lum_task,
2816 E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
2817 ,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,act_gsn );
2818 }
2819 }
2820 }
2821 }
2822
2823
2824
2825 /* calculation of indicees */
2826
2827 /* check vertical illuminance range */
2828 E_v2 = E_v;
2829
2830 if (E_v2 < 100) {
2831 fprintf(stderr,
2832 "Notice: Vertical illuminance is below 100 lux !!\n");
2833 }
2834 dgp =
2835 get_dgp(p, E_v2, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
2836 /* low light correction */
2837 if (E_v < 1000) {
2838 low_light_corr=1.0*exp(0.024*E_v-4)/(1+exp(0.024*E_v-4));} else {low_light_corr=1.0 ;}
2839 dgp =low_light_corr*dgp;
2840
2841 /* age correction */
2842
2843 if (age_corr == 1) {
2844 age_corr_factor=1.0/(1.1-0.5*age/100.0);
2845 }
2846 dgp =age_corr_factor*dgp;
2847
2848
2849 if (ext_vill == 1) {
2850 if (E_vl_ext < 100) {
2851 fprintf(stderr,
2852 "Notice: Vertical illuminance is below 100 lux !!\n");
2853 }
2854 }
2855
2856 if (calcfast == 0) {
2857 lum_a= E_v2/3.1415927;
2858 dgi = get_dgi(p, lum_backg, igs, posindex_2);
2859 ugr = get_ugr(p, lum_backg, igs, posindex_2);
2860 ugp = get_ugp (p,lum_backg, igs, posindex_2);
2861 ugr_exp = get_ugr_exp (p,lum_backg_cos,lum_a, igs, posindex_2);
2862 dgi_mod = get_dgi_mod(p, lum_a, igs, posindex_2);
2863 cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
2864 dgr = get_dgr(p, avlum, igs, posindex_2);
2865 vcp = get_vcp(dgr);
2866 Lveil = get_disability(p, avlum, igs);
2867 if (no_glaresources == 0) {
2868 dgi = 0.0;
2869 ugr = 0.0;
2870 ugp = 0.0;
2871 ugr_exp = 0.0;
2872 dgi_mod = 0.0;
2873 cgi = 0.0;
2874 dgr = 0.0;
2875 vcp = 100.0;
2876 }
2877 }
2878 /* check dgp range */
2879 if (dgp <= 0.2) {
2880 fprintf(stderr,
2881 "Notice: Low brightness scene. dgp below 0.2! dgp might underestimate glare sources\n");
2882 }
2883 if (E_v < 380) {
2884 fprintf(stderr,
2885 "Notice: Low brightness scene. Vertical illuminance less than 380 lux! dgp might underestimate glare sources\n");
2886 }
2887
2888
2889
2890 if (output == 0) {
2891 if (detail_out == 1) {
2892 if (ext_vill == 1) {
2893 dgp_ext =
2894 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2895 posindex_2);
2896 dgp = dgp_ext;
2897 if (E_vl_ext < 1000) {
2898 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 ;}
2899 dgp =low_light_corr*dgp;
2900 dgp =age_corr_factor*dgp;
2901 }
2902 muc_rvar_get_median(s_noposweight,lum_nopos_median);
2903 muc_rvar_get_median(s_posweight,lum_pos_median);
2904 muc_rvar_get_median(s_posweight2,lum_pos2_median);
2905
2906
2907
2908
2909 printf
2910 ("dgp,av_lum,E_v,lum_backg,E_v_dir,dgi,ugr,vcp,cgi,lum_sources,omega_sources,Lveil,Lveil_cie,dgr,ugp,ugr_exp,dgi_mod,av_lum_pos,av_lum_pos2,med_lum,med_lum_pos,med_lum_pos2: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
2911 dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
2912 lum_sources, omega_sources, Lveil,Lveil_cie_sum,dgr,ugp,ugr_exp,dgi_mod,lum_pos_mean,lum_pos2_mean/sang,lum_nopos_median[0],lum_pos_median[0],lum_pos2_median[0]);
2913 } else {
2914 if (detail_out2 == 1) {
2915
2916 printf
2917 ("dgp,dgi,ugr,vcp,cgi,dgp_ext,Ev_calc,abs_max,Lveil: %f %f %f %f %f %f %f %f %f \n",
2918 dgp, dgi, ugr, vcp, cgi, dgp_ext, E_v, abs_max,
2919 Lveil);
2920
2921 } else {
2922 if (ext_vill == 1) {
2923 dgp_ext = get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,posindex_2);
2924
2925 if (E_vl_ext < 1000) {
2926 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 ;}
2927 dgp =low_light_corr*dgp_ext;
2928 dgp =age_corr_factor*dgp;
2929 }
2930 printf("dgp,dgi,ugr,vcp,cgi,Lveil: %f %f %f %f %f %f \n",
2931 dgp, dgi, ugr, vcp, cgi, Lveil);
2932
2933 }
2934 }
2935 } else {
2936 dgp_ext =
2937 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2938 posindex_2);
2939 dgp = dgp_ext;
2940 if (E_vl_ext < 1000) {
2941 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 ;}
2942 dgp =low_light_corr*dgp;
2943
2944 if (calcfast == 2) {
2945
2946 lum_backg_cos=(E_vl_ext-dir_ill)/3.1415927;
2947 ugr = get_ugr(p, lum_backg_cos, igs, posindex_2);
2948 printf("%f %f \n", dgp,ugr);
2949 }else{
2950 printf("%f\n", dgp);
2951 }
2952 }
2953
2954
2955 /*printf("lowlight=%f\n", low_light_corr); */
2956
2957
2958 /* printf("hallo \n");
2959
2960
2961 apply_disability = 1;
2962 disability_thresh = atof(argv[++i]);
2963 coloring of disability glare sources red, remove all other colors
2964
2965
2966
2967 this function was removed because of bugs....
2968 has to be re-written from scratch....
2969
2970
2971 */
2972
2973
2974
2975
2976
2977
2978
2979
2980 /*output */
2981 if (checkfile == 1) {
2982
2983
2984 pict_write(p, file_out);
2985 }
2986
2987
2988
2989 return EXIT_SUCCESS;
2990 exit(0);
2991
2992 userr:
2993 fprintf(stderr,
2994 "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",
2995 progname);
2996 exit(1);
2997 }
2998
2999
3000
3001 #endif