ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/util/evalglare.c
Revision: 3.08
Committed: Wed Oct 1 17:48:15 2025 UTC (3 days, 21 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 3.07: +16 -10 lines
Log Message: 
fix(evalglare): Jan fixed issue with border masking

File Contents

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