ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.11
Committed: Tue Jun 30 15:57:30 2020 UTC (3 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R3, HEAD
Changes since 2.10: +385 -28 lines
Log Message:
feat(evalglare): new version of evalglare from Jan Wienold

File Contents

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