ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.12
Committed: Tue Jun 3 21:31:51 2025 UTC (2 days, 23 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.11: +2 -4 lines
Log Message: 
refactor: More consistent use of global char * progname and fixargv0()

File Contents

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