ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/evalglare.c
Revision: 3.05
Committed: Thu Sep 11 17:40:45 2025 UTC (3 weeks, 2 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.12: +358 -71 lines
Log Message: 
Numerous changes and fixes from Jan Wienold

File Contents

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