ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.8
Committed: Fri Aug 31 16:01:45 2018 UTC (5 years, 7 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R2
Changes since 2.7: +55 -54 lines
Log Message: 
New versions of evalglare and gendaylit from Jan Wienold

File Contents

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