ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/evalglare.c
Revision: 3.07
Committed: Tue Sep 16 15:43:14 2025 UTC (12 days, 7 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 3.06: +3 -2 lines
Log Message: 
fix(evalglare): Compiler issue under Windows

File Contents

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