ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/evalglare.c
Revision: 3.06
Committed: Mon Sep 15 20:21:24 2025 UTC (2 weeks, 2 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.05: +1 -3 lines
Log Message: 
fix(evalglare): Removed redundant progname definition

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: evalglare.c,v 3.05 2025/09/11 17:40:45 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 VIEW userview = STDVIEW;
1561 int gotuserview = 0;
1562
1563 struct muc_rvar* s_mask;
1564 s_mask = muc_rvar_create();
1565 muc_rvar_set_dim(s_mask, 1);
1566 muc_rvar_clear(s_mask);
1567 struct muc_rvar* s_band;
1568 s_band = muc_rvar_create();
1569 muc_rvar_set_dim(s_band, 1);
1570 muc_rvar_clear(s_band);
1571 struct muc_rvar* s_z1;
1572 s_z1 = muc_rvar_create();
1573 muc_rvar_set_dim(s_z1, 1);
1574 muc_rvar_clear(s_z1);
1575
1576 struct muc_rvar* s_z2;
1577 s_z2 = muc_rvar_create();
1578 muc_rvar_set_dim(s_z2, 1);
1579 muc_rvar_clear(s_z2);
1580
1581 struct muc_rvar* s_noposweight;
1582 s_noposweight = muc_rvar_create();
1583 muc_rvar_set_dim(s_noposweight, 1);
1584 muc_rvar_clear(s_noposweight);
1585
1586 struct muc_rvar* s_posweight;
1587 s_posweight = muc_rvar_create();
1588 muc_rvar_set_dim(s_posweight, 1);
1589 muc_rvar_clear(s_posweight);
1590
1591 struct muc_rvar* s_posweight2;
1592 s_posweight2 = muc_rvar_create();
1593 muc_rvar_set_dim(s_posweight2, 1);
1594 muc_rvar_clear(s_posweight2);
1595
1596 /*set required user view parameters to invalid values*/
1597 patchmode=0;
1598 patch_angle=25.0;
1599 patch_pixdistance_x=0;
1600 patch_pixdistance_y=0;
1601 lowlight=0;
1602 multi_image_mode=0;
1603 lastpixelwas_peak=0;
1604 num_images=0;
1605 dir_ill=0.0;
1606 delta_E=0.0;
1607 no_glaresources=0;
1608 n_corner_px=0;
1609 zero_corner_px=0;
1610 force=0;
1611 dist=0.0;
1612 u_r=0.33333;
1613 u_g=0.33333;
1614 u_b=0.33333;
1615 band_angle=0;
1616 angle_z1=0;
1617 angle_z2=0;
1618 x_zone=0;
1619 y_zone=0;
1620 per_75_z2=0;
1621 per_95_z2=0;
1622 lum_pos_mean=0;
1623 lum_pos2_mean=0;
1624 lum_band_av = 0.0;
1625 omega_band = 0.0;
1626 pgsv = 0.0 ;
1627 pgsv_con = 0.0 ;
1628 pgsv_sat = 0.0 ;
1629 E_v_mask = 0.0;
1630 Ez1 = 0.0;
1631 Ez2 = 0.0;
1632 lum_z1_av = 0.0;
1633 omega_z1 = 0.0;
1634 lum_z2_av = 0.0;
1635 omega_z2 = 0.0 ;
1636 i_z1 = 0;
1637 i_z2 = 0;
1638 zones = 0;
1639 lum_z2_av = 0.0;
1640 uniform_gs = 0;
1641 apply_disability = 0;
1642 disability_thresh = 0;
1643 Lveil_cie_sum=0.0;
1644 skip_second_scan=0;
1645 lum_total_max=0.0;
1646 calcfast=0;
1647 age_corr_factor = 1.0;
1648 age_corr = 0;
1649 age = 20.0;
1650 userview.horiz = 0;
1651 userview.vert = 0;
1652 userview.type = 0;
1653 dgp_ext = 0;
1654 E_vl_ext = 0.0;
1655 new_lum_max = 0.0;
1656 lum_max = 0.0;
1657 omegat = 0.0;
1658 yt = 0;
1659 xt = 0;
1660 x_disk=0;
1661 y_disk=0;
1662 angle_disk=0;
1663 yfillmin = 0;
1664 yfillmax = 0;
1665 cut_view = 0;
1666 cut_view_type = 0;
1667 setvalue = 2e09;
1668 omega_cos_contr = 0.0;
1669 lum_ideal = 0.0;
1670 max_angle = 0.2;
1671 lum_thres = 2000.0;
1672 lum_task = 0.0;
1673 task_lum = 0;
1674 sgs = 0;
1675 splithigh = 1;
1676 detail_out = 0;
1677 detail_out2 = 0;
1678 posindex_picture = 0;
1679 checkfile = 0;
1680 ext_vill = 0;
1681 fill = 0;
1682 a1 = 2.0;
1683 a2 = 1.0;
1684 a3 = 1.87;
1685 a4 = 2.0;
1686 a5 = 1.0;
1687 c1 = 5.87e-05;
1688 c2 = 0.092;
1689 c3 = 0.159;
1690 non_cos_lb = 0;
1691 posindex_2 = 0;
1692 task_color = 0;
1693 limit = 50000.0;
1694 set_lum_max = 0;
1695 set_lum_max2 = 0;
1696 img_corr=0;
1697 abs_max = 0;
1698 fixargv0(argv[0]);
1699 E_v_contr = 0.0;
1700 strcpy(version, "3.05 release 25.06.2024 by J.Wienold");
1701 low_light_corr=1.0;
1702 output = 0;
1703 calc_vill = 0;
1704 band_avlum = -99;
1705 band_calc = 0;
1706 masking = 0;
1707 lum_mask[0]=0.0;
1708 lum_mask_av=0.0;
1709 omega_mask=0.0;
1710 i_mask=0;
1711 actual_igs=0;
1712 LUM_replace=0;
1713 thres_activate=0;
1714 /* command line for output picture*/
1715
1716 cline = (char *) malloc(CLINEMAX+100);
1717 cline[0] = '\0';
1718 for (i = 0; i < argc; i++) {
1719 /* fprintf(stderr, "%d %d \n",i,strlen(argv[i]));*/
1720 if (strlen(cline)+strlen(argv[i])+strlen(RELEASENAME)+2 >=CLINEMAX) {
1721 exit (-1);
1722 }
1723 else {
1724 strcat(cline, argv[i]);
1725 strcat(cline, " ");
1726 }
1727 }
1728 strcat(cline, "\n");
1729 strcat(cline, RELEASENAME);
1730 strcat(cline, "\n");
1731
1732
1733 /* program options */
1734
1735 for (i = 1; i < argc && argv[i][0] == '-'; i++) {
1736 /* expand arguments */
1737 while ((rval = expandarg(&argc, &argv, i)) > 0);
1738 if (rval < 0) {
1739 fprintf(stderr, "%s: cannot expand '%s'", argv[0], argv[i]);
1740 exit(1);
1741 }
1742 rval = getviewopt(&userview, argc - i, argv + i);
1743 if (rval >= 0) {
1744 i += rval;
1745 gotuserview++;
1746 continue;
1747 }
1748 switch (argv[i][1]) {
1749 case 'a':
1750 age = atof(argv[++i]);
1751 age_corr = 1;
1752 printf("age factor not supported any more \n");
1753 exit(1);
1754 break;
1755 case 'A':
1756 masking = 1;
1757 detail_out = 1;
1758 strcpy(maskfile, argv[++i]);
1759 pict_read(pm,maskfile );
1760
1761 break;
1762 case 'b':
1763 lum_thres = atof(argv[++i]);
1764 lum_source =lum_thres;
1765 thres_activate = 1;
1766 break;
1767 case 'c':
1768 checkfile = 1;
1769 strcpy(file_out, argv[++i]);
1770 break;
1771 case 'u':
1772 uniform_gs = 1;
1773 u_r = atof(argv[++i]);
1774 u_g = atof(argv[++i]);
1775 u_b = atof(argv[++i]);
1776 break;
1777 case 'r':
1778 max_angle = atof(argv[++i]);
1779 break;
1780 case 'z':
1781 patch_angle = atof(argv[++i]);
1782 patchmode= atoi(argv[++i]);
1783 if ( patchmode == 3) { output=4;}
1784
1785 /* 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...) */
1786 break;
1787
1788 case 's':
1789 sgs = 1;
1790 break;
1791 case 'f':
1792 force = 1;
1793 break;
1794 case 'k':
1795 apply_disability = 1;
1796 disability_thresh = atof(argv[++i]);
1797 break;
1798
1799 case 'p':
1800 posindex_picture = 1;
1801 break;
1802 case 'P':
1803 posindex_2 = atoi(argv[++i]);
1804 break;
1805
1806
1807 case 'y':
1808 splithigh = 1;
1809 break;
1810 case 'x':
1811 splithigh = 0;
1812 break;
1813 case 'Y':
1814 splithigh = 1;
1815 limit = atof(argv[++i]);
1816 break;
1817
1818 case 'i':
1819 ext_vill = 1;
1820 E_vl_ext = atof(argv[++i]);
1821 break;
1822 case 'I':
1823 ext_vill = 1;
1824 fill = 1;
1825 E_vl_ext = atof(argv[++i]);
1826 yfillmax = atoi(argv[++i]);
1827 yfillmin = atoi(argv[++i]);
1828 break;
1829 case 'd':
1830 detail_out = 1;
1831 break;
1832 case 'D':
1833 detail_out2 = 1;
1834 break;
1835 case 'm':
1836 img_corr = 1;
1837 set_lum_max = 1;
1838 lum_max = atof(argv[++i]);
1839 new_lum_max = atof(argv[++i]);
1840 strcpy(file_out2, argv[++i]);
1841 /* printf("max lum set to %f\n",new_lum_max);*/
1842 break;
1843
1844 case 'M':
1845 img_corr = 1;
1846 set_lum_max2 = 1;
1847 lum_max = atof(argv[++i]);
1848 E_vl_ext = atof(argv[++i]);
1849 strcpy(file_out2, argv[++i]);
1850 /* printf("max lum set to %f\n",new_lum_max);*/
1851 break;
1852 case 'N':
1853 img_corr = 1;
1854 set_lum_max2 = 2;
1855 x_disk = atoi(argv[++i]);
1856 y_disk = atoi(argv[++i]);
1857 angle_disk = atof(argv[++i]);
1858 E_vl_ext = atof(argv[++i]);
1859 strcpy(file_out2, argv[++i]);
1860 /* printf("max lum set to %f\n",new_lum_max);*/
1861 break;
1862 case 'O':
1863 img_corr = 1;
1864 set_lum_max2 = 3;
1865 x_disk = atoi(argv[++i]);
1866 y_disk = atoi(argv[++i]);
1867 angle_disk = atof(argv[++i]);
1868 LUM_replace = atof(argv[++i]);
1869 strcpy(file_out2, argv[++i]);
1870 /* printf("max lum set to %f\n",new_lum_max);*/
1871 break;
1872
1873
1874 /* deactivated case 'n':
1875 non_cos_lb = 0;
1876 break;
1877 */
1878 case 'q':
1879 non_cos_lb = atoi(argv[++i]);
1880 break;
1881
1882 case 't':
1883 task_lum = 1;
1884 xt = atoi(argv[++i]);
1885 yt = atoi(argv[++i]);
1886 omegat = atof(argv[++i]);
1887 task_color = 0;
1888 break;
1889 case 'T':
1890 task_lum = 1;
1891 xt = atoi(argv[++i]);
1892 yt = atoi(argv[++i]);
1893 /* omegat= DEG2RAD(atof(argv[++i]));*/
1894 omegat = atof(argv[++i]);
1895 task_color = 1;
1896 break;
1897 case 'l':
1898 zones = 1;
1899 x_zone = atoi(argv[++i]);
1900 y_zone = atoi(argv[++i]);
1901 angle_z1 = atof(argv[++i]);
1902 angle_z2 = -1;
1903 break;
1904 case 'L':
1905 zones = 2;
1906 x_zone = atoi(argv[++i]);
1907 y_zone = atoi(argv[++i]);
1908 angle_z1 = atof(argv[++i]);
1909 angle_z2 = atof(argv[++i]);
1910 break;
1911
1912
1913 case 'B':
1914 band_calc = 1;
1915 band_angle = atof(argv[++i]);
1916 band_color = 1;
1917 break;
1918
1919
1920 case 'w':
1921 a1 = atof(argv[++i]);
1922 a2 = atof(argv[++i]);
1923 a3 = atof(argv[++i]);
1924 a4 = atof(argv[++i]);
1925 a5 = atof(argv[++i]);
1926 c1 = atof(argv[++i]);
1927 c2 = atof(argv[++i]);
1928 c3 = atof(argv[++i]);
1929 break;
1930 case 'V':
1931 calc_vill = 1;
1932 break;
1933 case 'G':
1934 cut_view = 1;
1935 cut_view_type= atof(argv[++i]);
1936 break;
1937 case 'g':
1938 cut_view = 2;
1939 cut_view_type= atof(argv[++i]);
1940 break;
1941
1942 /*case 'v':
1943 printf("evalglare %s \n",version);
1944 exit(1); */
1945 case 'C':
1946 strcpy(correction_type,argv[++i]);
1947
1948 if (!strcmp(correction_type,"l-") ){
1949 /* printf("low light off!\n"); */
1950 lowlight = 0; }
1951 if (!strcmp(correction_type,"l+") ){
1952 /* printf("low light on!\n"); */
1953 lowlight = 1; }
1954 if (!strcmp(correction_type,"0") ){
1955 /* printf("all corrections off!\n"); */
1956 lowlight = 0; }
1957
1958 break;
1959
1960 /*case 'v':
1961 printf("evalglare %s \n",version);
1962 exit(1); */
1963
1964 case '1':
1965 output = 1;
1966 break;
1967 case '2':
1968 output = 2;
1969 dir_ill = atof(argv[++i]);
1970 break;
1971 case '3':
1972 output = 3;
1973 break;
1974 case '4':
1975 lowlight = 0;
1976 break;
1977 case 'Q':
1978 multi_image_mode=1;
1979 output= 3;
1980 calcfast=1;
1981 num_images =atoi(argv[++i]);
1982 num_scans =atoi(argv[++i]);
1983 temp_image_name = malloc(sizeof(char*)*num_images);
1984
1985 x_temp_img=(int *) malloc(sizeof(int) * num_images);
1986 y_temp_img=(int *) malloc(sizeof(int) * num_images);
1987 scale_image_scans=(double *) malloc(sizeof(double) * (num_scans+1)*num_images);
1988 /* iterate through all images and allocate 256 characters to each: */
1989 for (j = 0; j < num_images; j++) {
1990 scale_image_scans[j*(num_scans+1)]=1.0;
1991 temp_image_name[j] = malloc(256*sizeof(char));
1992 strcpy(temp_image_name[j], argv[++i]);
1993 x_temp_img[j] = atoi(argv[++i]);
1994 y_temp_img[j] = atoi(argv[++i]);
1995 for (jj=1;jj<=num_scans;jj++){
1996 scale_image_scans[j*(num_scans+1)+jj]=atof(argv[++i]);}
1997 }
1998
1999
2000 break;
2001 case 'v':
2002 if (argv[i][2] == '\0') {
2003 printf("%s \n",RELEASENAME);
2004 exit(1);
2005 }
2006 if (argv[i][2] != 'f')
2007 goto userr;
2008 rval = viewfile(argv[++i], &userview, NULL);
2009 if (rval < 0) {
2010 fprintf(stderr,
2011 "%s: cannot open view file \"%s\"\n",
2012 progname, argv[i]);
2013 exit(1);
2014 } else if (rval == 0) {
2015 fprintf(stderr,
2016 "%s: bad view file \"%s\"\n", progname, argv[i]);
2017 exit(1);
2018 } else {
2019 gotuserview++;
2020 }
2021 break;
2022
2023
2024 default:
2025 goto userr;
2026 }
2027 }
2028
2029 /* set multiplier for task method to 5, if not specified */
2030
2031
2032
2033 if ( task_lum == 1 && thres_activate == 0){
2034 lum_thres = 5.0;
2035 }
2036 /*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/
2037
2038 if (output == 1 && ext_vill == 1 ) {
2039 calcfast=1;
2040 }
2041
2042 if (output == 2 && ext_vill == 1 ) {
2043 calcfast=2;
2044 }
2045
2046 /*masking and zoning cannot be applied at the same time*/
2047
2048 if (masking ==1 && zones >0) {
2049 fprintf(stderr, " masking and zoning cannot be activated at the same time!\n");
2050 exit(1);
2051 }
2052
2053 /* read picture file */
2054 if (i == argc) {
2055 SET_FILE_BINARY(stdin);
2056 FILE *fp = fdopen(fileno(stdin), "rb");
2057 if (!(fp)) {
2058 fprintf(stderr, "fdopen on stdin failed\n");
2059 return EXIT_FAILURE;
2060 }
2061 if (!(pict_read_fp(p, fp)))
2062 return EXIT_FAILURE;;
2063 } else {
2064 if (!pict_read(p, argv[i]))
2065 return EXIT_FAILURE;
2066 }
2067 if (gotuserview) {
2068 if (p->valid_view)
2069 fprintf(stderr,
2070 "warning: overriding image view by commandline argument\n");
2071 if ((userview.horiz == 0) || (userview.vert == 0) || (userview.type == 0)) {
2072 fprintf(stderr, "error: if specified, a view must at least contain -vt -vv and -vh\n");
2073 return EXIT_FAILURE;
2074 }
2075 p->view = userview;
2076 if (!(pict_update_view(p))) {
2077 fprintf(stderr, "error: invalid view specified");
2078 return EXIT_FAILURE;
2079 }
2080 pict_update_evalglare_caches(p);
2081 } else if (!(p->valid_view)) {
2082 fprintf(stderr, "error: no valid view specified\n");
2083 return EXIT_FAILURE;
2084 }
2085
2086
2087
2088 /* fprintf(stderr,"Picture read!");*/
2089
2090 /* command line for output picture*/
2091
2092 pict_set_comment(p, cline);
2093
2094 /* several checks */
2095 if (p->view.type == VT_PAR) {
2096 fprintf(stderr, "error: wrong view type! must not be parallel ! \n");
2097 exit(1);
2098 }
2099
2100 if ( patchmode > 0 && p->view.type != VT_ANG) {
2101
2102 fprintf(stderr, "error: Patchmode only possible with view type vta ! Stopping... \n");
2103 exit(1);
2104
2105 }
2106
2107
2108 if (p->view.type == VT_PER) {
2109 fprintf(stderr, "warning: image has perspective view type specified in header ! \n");
2110 }
2111
2112 if (masking == 1) {
2113
2114 if (pict_get_xsize(p)!=pict_get_xsize(pm) || pict_get_ysize(p)!=pict_get_ysize(pm)) {
2115 fprintf(stderr, "error: masking image has other resolution than main image ! \n");
2116 fprintf(stderr, "size must be identical \n");
2117 printf("resolution main image : %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
2118 printf("resolution masking image : %dx%d\n",pict_get_xsize(pm),pict_get_ysize(pm));
2119 exit(1);
2120
2121 }
2122
2123 }
2124 /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
2125
2126 /* Check size of search radius */
2127 rmx = (pict_get_xsize(p) / 2);
2128 rmy = (pict_get_ysize(p) / 2);
2129 rx = (pict_get_xsize(p) / 2 + 10);
2130 ry = (pict_get_ysize(p) / 2 + 10);
2131 r_center =
2132 acos(DOT(pict_get_cached_dir(p, rmx, rmy),
2133 pict_get_cached_dir(p, rx, ry))) * 2 / 10;
2134 search_pix = max_angle / r_center;
2135 if (search_pix < 1.0) {
2136 fprintf(stderr,
2137 "warning: search radius less than 1 pixel! deactivating smoothing and peak extraction...\n");
2138 splithigh = 0;
2139 sgs = 0;
2140
2141 } else {
2142 if (search_pix < 3.0) {
2143 fprintf(stderr,
2144 "warning: search radius less than 3 pixels! -> %f \n",
2145 search_pix);
2146
2147 }
2148 }
2149
2150
2151 /* Check task position */
2152
2153 if (task_lum == 1) {
2154 if (xt >= pict_get_xsize(p) || yt >= pict_get_ysize(p) || xt < 0
2155 || yt < 0) {
2156 fprintf(stderr,
2157 "error: task position outside picture!! exit...");
2158 exit(1);
2159 }
2160 }
2161
2162
2163 /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
2164
2165 sang = 0.0;
2166 E_v = 0.0;
2167 E_v_dir = 0.0;
2168 avlum = 0.0;
2169 pict_new_gli(p);
2170 igs = 0;
2171 pict_get_z1_gsn(p,igs) = 0;
2172 pict_get_z2_gsn(p,igs) = 0;
2173
2174 if (multi_image_mode<1) {
2175
2176
2177 /* cut out GUTH field of view and exit without glare evaluation */
2178 if (cut_view==2) {
2179 if (cut_view_type==1) {
2180 cut_view_1(p);
2181 }else{
2182
2183 if (cut_view_type==2) {
2184 cut_view_2(p);
2185 }else{
2186 if (cut_view_type==3) {
2187 fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
2188 cut_view_3(p);
2189 }else{
2190 fprintf(stderr,"error: no valid option for view cutting!!");
2191 exit(1);
2192 }
2193 }}
2194 pict_write(p, file_out);
2195 exit(1); }
2196
2197
2198
2199
2200
2201
2202 /* write positionindex into checkfile and exit, activated by option -p */
2203 if (posindex_picture == 1) {
2204 for (x = 0; x < pict_get_xsize(p); x++)
2205 for (y = 0; y < pict_get_ysize(p); y++) {
2206 if (pict_get_hangle
2207 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2208 if (pict_is_validpixel(p, x, y)) {
2209 lum =
2210 get_posindex(p, x, y,
2211 posindex_2) / WHTEFFICACY;
2212
2213 pict_get_color(p, x, y)[RED] = lum;
2214 pict_get_color(p, x, y)[GRN] = lum;
2215 pict_get_color(p, x, y)[BLU] = lum;
2216 lum_task = luminance(pict_get_color(p, x, y));
2217 /*printf("x,y,posindex=%i %i %f %f\n",x,y,lum*WHTEFFICACY,lum_task);*/
2218 }
2219 }
2220 }
2221 pict_write(p, file_out);
2222 exit(1);
2223 }
2224
2225
2226 /* calculation of vertical illuminance, average luminance, in case of filling activated-> fill picture */
2227 /* fill, if necessary from 0 to yfillmin */
2228
2229 if (fill == 1) {
2230 for (x = 0; x < pict_get_xsize(p); x++)
2231 for (y = yfillmin; y > 0; y = y - 1) {
2232 y1 = y + 1;
2233 lum = luminance(pict_get_color(p, x, y1));
2234 pict_get_color(p, x, y)[RED] = lum / 179.0;
2235 pict_get_color(p, x, y)[GRN] = lum / 179.0;
2236 pict_get_color(p, x, y)[BLU] = lum / 179.0;
2237 }
2238 }
2239
2240 if (calcfast == 0) {
2241 for (x = 0; x < pict_get_xsize(p); x++)
2242 for (y = 0; y < pict_get_ysize(p); y++) {
2243 lum = luminance(pict_get_color(p, x, y));
2244 dist=sqrt((x-rmx)*(x-rmx)+(y-rmy)*(y-rmy));
2245 if (dist>((rmx+rmy)/2)) {
2246 n_corner_px=n_corner_px+1;
2247 if (lum<7.0) {
2248 zero_corner_px=zero_corner_px+1;
2249 }
2250 }
2251
2252
2253 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2254 if (pict_is_validpixel(p, x, y)) {
2255 if (fill == 1 && y >= yfillmax) {
2256 y1 = y - 1;
2257 lum = luminance(pict_get_color(p, x, y1));
2258 pict_get_color(p, x, y)[RED] = lum / 179.0;
2259 pict_get_color(p, x, y)[GRN] = lum / 179.0;
2260 pict_get_color(p, x, y)[BLU] = lum / 179.0;
2261 }
2262
2263 if (lum > abs_max) {
2264 abs_max = lum;
2265 }
2266 /* set luminance restriction, if -m is set */
2267 if (img_corr == 1 ) {
2268 if (set_lum_max == 1 && lum >= lum_max) {
2269 pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
2270 pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
2271 pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
2272 /* printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
2273 lum = luminance(pict_get_color(p, x, y));
2274 }
2275 if (set_lum_max2 == 1 && lum >= lum_max) {
2276
2277 E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2278 omega_cos_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2279 }
2280 if (set_lum_max2 == 2 ) {
2281 r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2282 if (x_disk == x && y_disk==y ) r_actual=0.0;
2283
2284 act_lum = luminance(pict_get_color(p, x, y));
2285
2286 if (r_actual <= angle_disk) {
2287 E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2288 omega_cos_contr += DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2289 }
2290
2291
2292
2293 }
2294 }
2295 om = pict_get_omega(p, x, y);
2296 sang += om;
2297 E_v += DOT(p->view.vdir, pict_get_cached_dir(p, x,y)) * om *lum;
2298 avlum += om * lum;
2299 pos=get_posindex(p, x, y,posindex_2);
2300
2301 lum_pos[0]=lum;
2302 muc_rvar_add_sample(s_noposweight,lum_pos);
2303 lum_pos[0]=lum/pos;
2304 lum_pos_mean+=lum_pos[0]*om;
2305 muc_rvar_add_sample(s_posweight, lum_pos);
2306 lum_pos[0]=lum_pos[0]/pos;
2307 lum_pos2_mean+=lum_pos[0]*om;
2308 muc_rvar_add_sample(s_posweight2,lum_pos);
2309
2310
2311
2312 } else {
2313 pict_get_color(p, x, y)[RED] = 0.0;
2314 pict_get_color(p, x, y)[GRN] = 0.0;
2315 pict_get_color(p, x, y)[BLU] = 0.0;
2316
2317 }
2318 }else {
2319 pict_get_color(p, x, y)[RED] = 0.0;
2320 pict_get_color(p, x, y)[GRN] = 0.0;
2321 pict_get_color(p, x, y)[BLU] = 0.0;
2322
2323 }
2324 }
2325
2326 /* check if image has black corners AND a perspective view */
2327
2328 if (zero_corner_px/n_corner_px > 0.70 && (p->view.type == VT_PER) ) {
2329 printf (" corner pixels are to %f %% black! \n",zero_corner_px/n_corner_px*100);
2330 printf("error! The image has black corners AND header contains a perspective view type definition !!!\n") ;
2331
2332 if (force==0) {
2333 printf("stopping...!!!!\n") ;
2334
2335 exit(1);
2336 }
2337 }
2338
2339
2340
2341
2342 lum_pos_mean= lum_pos_mean/sang;
2343 lum_pos2_mean= lum_pos2_mean/sang;
2344
2345 if ((set_lum_max2 >= 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0 ) || set_lum_max2==3) {
2346
2347 if (set_lum_max2<3){
2348 lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
2349 if (set_lum_max2 == 2 && lum_ideal >= 2e9) {
2350 printf("warning! luminance of replacement pixels would be larger than 2e9 cd/m2. Value set to 2e9cd/m2!\n") ;
2351 }
2352
2353 if (lum_ideal > 0 && lum_ideal < setvalue) {
2354 setvalue = lum_ideal;
2355 }
2356 printf("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f %f %f %f %f\n",
2357 lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
2358 }else{setvalue=LUM_replace;
2359 }
2360
2361
2362 for (x = 0; x < pict_get_xsize(p); x++)
2363 for (y = 0; y < pict_get_ysize(p); y++) {
2364 if (pict_get_hangle
2365 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2366 if (pict_is_validpixel(p, x, y)) {
2367 lum = luminance(pict_get_color(p, x, y));
2368
2369
2370 if (set_lum_max2 == 1 && lum >= lum_max) {
2371
2372 pict_get_color(p, x, y)[RED] =
2373 setvalue / 179.0;
2374 pict_get_color(p, x, y)[GRN] =
2375 setvalue / 179.0;
2376 pict_get_color(p, x, y)[BLU] =
2377 setvalue / 179.0;
2378
2379 }else{ if(set_lum_max2 >1 ) {
2380 r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2381 if (x_disk == x && y_disk==y ) r_actual=0.0;
2382
2383 if (r_actual <= angle_disk) {
2384
2385 pict_get_color(p, x, y)[RED] =
2386 setvalue / 179.0;
2387 pict_get_color(p, x, y)[GRN] =
2388 setvalue / 179.0;
2389 pict_get_color(p, x, y)[BLU] =
2390 setvalue / 179.0;
2391
2392 }
2393 }
2394 }
2395 }
2396 }
2397 }
2398
2399
2400 pict_write(p, file_out2);
2401 exit(1);
2402 }
2403 if (set_lum_max == 1) {
2404 pict_write(p, file_out2);
2405
2406 }
2407
2408 if (calc_vill == 1) {
2409 printf("%f\n", E_v);
2410 exit(1);
2411 }
2412 }else{
2413 /* in fast calculation mode: ev=ev_ext and sang=2*pi */
2414 sang=2*3.14159265359;
2415 lum_task =E_vl_ext/sang;
2416 E_v=E_vl_ext;
2417 /* printf("calc fast!! %f %f\n", sang,lum_task);*/
2418
2419
2420 }
2421
2422 /* cut out GUTH field of view for glare evaluation */
2423 if (cut_view==1) {
2424 if (cut_view_type==1) {
2425 cut_view_1(p);
2426 }else{
2427
2428 if (cut_view_type==2) {
2429 cut_view_2(p);
2430 }else{
2431 if (cut_view_type==3) {
2432 fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
2433 cut_view_3(p);
2434 }else{
2435 fprintf(stderr,"error: no valid option for view cutting!!");
2436 exit(1);
2437 }
2438 }}
2439 }
2440
2441 if (calcfast == 0) {
2442 avlum = avlum / sang;
2443 lum_task = avlum;
2444 }
2445 /* if (ext_vill==1) {
2446 E_v=E_vl_ext;
2447 avlum=E_v/3.1415927;
2448 } */
2449
2450 if (task_lum == 1) {
2451 lum_task = get_task_lum(p, xt, yt, omegat, task_color);
2452 }
2453 lum_source = lum_thres * lum_task;
2454 /* printf("Task Luminance=%f\n",lum_task);
2455 pict_write(p,"task.pic");*/
2456
2457 if (lum_thres > 100) {
2458 lum_source = lum_thres;
2459 }
2460
2461 /* printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2462
2463 /* first glare source scan: find primary glare sources */
2464
2465 if (patchmode==0) {
2466 for (x = 0; x < pict_get_xsize(p); x++) {
2467 lastpixelwas_gs=0;
2468 /* lastpixelwas_peak=0; */
2469 for (y = 0; y < pict_get_ysize(p); y++) {
2470 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2471 if (pict_is_validpixel(p, x, y)) {
2472 act_lum = luminance(pict_get_color(p, x, y));
2473 if (act_lum > lum_source) {
2474 if (act_lum >lum_total_max) {
2475 lum_total_max=act_lum;
2476 }
2477 /* speed improvement first scan: when last pixel was glare source, then use this glare source number instead of searching.
2478 has been deactivated with v1.25, reactivated with v2.10 */
2479
2480 if (lastpixelwas_gs==0 || search_pix <= 1.0 ) {
2481 actual_igs = find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2482 }
2483 if (actual_igs == 0) {
2484 igs = igs + 1;
2485 pict_new_gli(p);
2486 pict_get_lum_min(p, igs) = HUGE_VAL;
2487 pict_get_Eglare(p,igs) = 0.0;
2488 pict_get_Dglare(p,igs) = 0.0;
2489 pict_get_z1_gsn(p,igs) = 0;
2490 pict_get_z2_gsn(p,igs) = 0;
2491 actual_igs = igs;
2492
2493 }
2494 /* no coloring any more here icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2495 pict_get_gsn(p, x, y) = actual_igs;
2496 pict_get_pgs(p, x, y) = 1;
2497 add_pixel_to_gs(p, x, y, actual_igs);
2498 lastpixelwas_gs=actual_igs;
2499
2500 } else {
2501 pict_get_pgs(p, x, y) = 0;
2502 lastpixelwas_gs=0;
2503 }
2504 }
2505 }
2506 }
2507 }
2508 } else {
2509 /* patchmode on!
2510 calculation only for angular projection!
2511
2512 */
2513
2514
2515 angle_v=p->view.vert;
2516 angle_h=p->view.horiz;
2517
2518
2519 /*printf ("angle_v,angle_h: %f %f \n",angle_v,angle_h) ;
2520
2521 setting up patches as glare sources */
2522
2523 patch_pixdistance_x=floor(pict_get_xsize(p)/angle_h*patch_angle);
2524 patch_pixdistance_y=floor(pict_get_ysize(p)/angle_v*patch_angle);
2525
2526 nx_patch=floor(angle_v/patch_angle)+1;
2527 ny_patch=floor(angle_h/patch_angle)+1;
2528
2529 y_offset=floor (patch_pixdistance_y/2);
2530 x_offset=floor (patch_pixdistance_x/2);
2531 /* 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) ; */
2532
2533 ix_offset=0;
2534 for (iy=1; iy<=ny_patch;iy++) {
2535
2536
2537
2538 for (ix=1; ix<=nx_patch;ix++) {
2539 igs = igs + 1;
2540 pict_new_gli(p);
2541
2542 pict_get_lum_min(p, igs) = HUGE_VAL;
2543 pict_get_Eglare(p,igs) = 0.0;
2544 pict_get_Dglare(p,igs) = 0.0;
2545 pict_get_z1_gsn(p,igs) = 0;
2546 pict_get_z2_gsn(p,igs) = 0;
2547 pict_get_dx_max(p,igs) = (x_offset+ix_offset*x_offset+(ix-1)*patch_pixdistance_x)*1.0;
2548 pict_get_dy_max(p,igs) = (y_offset+(iy-1)*patch_pixdistance_y)*1.0;
2549
2550 /* printf ("igs, x-patch, y-patch : %i %f %f \n",igs,pict_get_dx_max(p,igs),pict_get_dy_max(p,igs) ) ; */
2551
2552 }
2553 ix_offset=ix_offset+1;
2554 if (ix_offset==2) {
2555 ix_offset =0 ;
2556 }
2557
2558 }
2559 for (y = 0; y < pict_get_ysize(p); y++) {
2560 for (x = 0; x < pict_get_xsize(p); x++) {
2561
2562 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2563 if (pict_is_validpixel(p, x, y)) {
2564 act_lum = luminance(pict_get_color(p, x, y));
2565 if (act_lum > lum_source) {
2566 if (act_lum >lum_total_max) {
2567 lum_total_max=act_lum;
2568 }
2569
2570 y_lines = floor((y)/patch_pixdistance_y);
2571 xxx = (x+0.0)/(patch_pixdistance_x+0.0)-0.5*(y_lines % 2);
2572 if (xxx<0 ) { xxx=0.0 ;}
2573 i1 = y_lines*(nx_patch)+floor(xxx)+1;
2574 i2=0;
2575 add=0;
2576 if (y_lines % 2 == 1 ) {add=1;}
2577
2578 if (y >pict_get_dy_max(p,i1)) {
2579
2580 if (x > pict_get_dx_max(p,i1)) {
2581 i2=i1+nx_patch+add;
2582 }else {
2583 i2=i1+nx_patch-1+add;
2584 }
2585 }else {
2586
2587 if (x > pict_get_dx_max(p,i1)) {
2588 i2=i1-nx_patch+add;
2589 }else {
2590 i2=i1-nx_patch-1+add;
2591 }
2592 }
2593
2594
2595
2596
2597 if (i2 > igs || i2 < 1) {actual_igs=i1;}else{
2598 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))) ) {
2599 actual_igs=i1; }else {actual_igs=i2;}}
2600
2601
2602 pict_get_gsn(p, x, y) = actual_igs;
2603 pict_get_pgs(p, x, y) = 1;
2604 add_pixel_to_gs(p, x, y, actual_igs);
2605 /* setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2606 }
2607 }
2608 }
2609
2610
2611
2612
2613
2614 } }
2615
2616
2617
2618
2619 }
2620 /* pict_write(p,"firstscan.pic"); */
2621
2622
2623
2624
2625 if (calcfast ==1 || search_pix <= 1.0 || calcfast == 2 || patchmode > 0) {
2626 skip_second_scan=1;
2627 }
2628
2629
2630 /* second glare source scan: combine glare sources facing each other */
2631 change = 1;
2632 i = 0;
2633 while (change == 1 && skip_second_scan==0) {
2634 change = 0;
2635 i = i+1;
2636 for (x = 0; x < pict_get_xsize(p); x++) {
2637 for (y = 0; y < pict_get_ysize(p); y++) {
2638 if (pict_get_hangle
2639 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2640 checkpixels=1;
2641 before_igs = pict_get_gsn(p, x, y);
2642
2643 /* speed improvement: search for other nearby glare sources only, when at least one adjacend pixel has another glare source number */
2644 if (x >1 && x<pict_get_xsize(p)-2 && y >1 && y<pict_get_ysize(p)-2 ) {
2645 x1=x-1;
2646 x2=x+1;
2647 y1=y-1;
2648 y2=y+1;
2649 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) ) {
2650 checkpixels = 0;
2651 actual_igs = before_igs;
2652 } }
2653
2654
2655 if (pict_is_validpixel(p, x, y) && before_igs > 0 && checkpixels==1) {
2656 actual_igs =
2657 find_near_pgs(p, x, y, max_angle, before_igs,
2658 igs, 1);
2659 if (!(actual_igs == before_igs)) {
2660 change = 1;
2661 }
2662 if (before_igs > 0) {
2663 actual_igs = pict_get_gsn(p, x, y);
2664 /* icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2665 }
2666
2667 }
2668 }
2669 }
2670 /* pict_write(p,"secondscan.pic");*/
2671 }}
2672
2673 /*smoothing: add secondary glare sources */
2674 i_max = igs;
2675 if (sgs == 1) {
2676
2677 /* simplified search radius: constant for entire picture, alway circle not an angle due to performance */
2678 x = (pict_get_xsize(p) / 2);
2679 y = (pict_get_ysize(p) / 2);
2680 rx = (pict_get_xsize(p) / 2 + 10);
2681 ry = (pict_get_ysize(p) / 2 + 10);
2682 r_center =
2683 acos(DOT(pict_get_cached_dir(p, x, y),
2684 pict_get_cached_dir(p, rx, ry))) * 2 / 10;
2685 search_pix = max_angle / r_center / 1.75;
2686
2687 for (x = 0; x < pict_get_xsize(p); x++) {
2688 for (y = 0; y < pict_get_ysize(p); y++) {
2689 if (pict_get_hangle
2690 (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2691 if (pict_is_validpixel(p, x, y)
2692 && pict_get_gsn(p, x, y) == 0) {
2693 for (i = 1; i <= i_max; i++) {
2694
2695 if (pict_get_npix(p, i) > 0) {
2696 add_secondary_gs(p, x, y, search_pix, i, uniform_gs, u_r, u_g , u_b);
2697 }
2698 }
2699
2700 }
2701 }
2702 }
2703 }
2704
2705 }
2706
2707 /* extract extremes from glare sources to extra glare source */
2708 if (splithigh == 1 && lum_total_max>limit) {
2709 /* fprintf(stderr, " split of glare source!\n"); */
2710
2711 r_split = max_angle / 2.0;
2712 for (i = 0; i <= i_max; i++) {
2713
2714 if (pict_get_npix(p, i) > 0) {
2715 l_max = pict_get_lum_max(p, i);
2716 i_splitstart = igs + 1;
2717 if (l_max >= limit) {
2718 for (x = 0; x < pict_get_xsize(p); x++)
2719 for (y = 0; y < pict_get_ysize(p); y++) {
2720 if (pict_get_hangle
2721 (p, x, y, p->view.vdir, p->view.vup, &ang))
2722 {
2723 if (pict_is_validpixel(p, x, y)
2724 && luminance(pict_get_color(p, x, y))
2725 >= limit
2726 && pict_get_gsn(p, x, y) == i) {
2727 if (i_splitstart == (igs + 1)) {
2728 igs = igs + 1;
2729 pict_new_gli(p);
2730 pict_get_z1_gsn(p,igs) = 0;
2731 pict_get_z2_gsn(p,igs) = 0;
2732
2733 pict_get_Eglare(p,igs) = 0.0;
2734 pict_get_Dglare(p,igs) = 0.0;
2735 pict_get_lum_min(p, igs) =
2736 99999999999999.999;
2737 i_split = igs;
2738 } else {
2739 i_split =
2740 find_split(p, x, y, r_split,
2741 i_splitstart, igs);
2742 }
2743 if (i_split == 0) {
2744 igs = igs + 1;
2745 pict_new_gli(p);
2746 pict_get_z1_gsn(p,igs) = 0;
2747 pict_get_z2_gsn(p,igs) = 0;
2748
2749 pict_get_Eglare(p,igs) = 0.0;
2750 pict_get_Dglare(p,igs) = 0.0;
2751 pict_get_lum_min(p, igs) =
2752 99999999999999.999;
2753 i_split = igs;
2754 }
2755 split_pixel_from_gs(p, x, y, i_split, uniform_gs, u_r, u_g , u_b);
2756
2757 }
2758 }
2759 }
2760
2761 }
2762 change = 1;
2763 while (change == 1) {
2764 change = 0;
2765 for (x = 0; x < pict_get_xsize(p); x++)
2766 for (y = 0; y < pict_get_ysize(p); y++) {
2767 before_igs = pict_get_gsn(p, x, y);
2768 if (before_igs >= i_splitstart) {
2769 if (pict_get_hangle
2770 (p, x, y, p->view.vdir, p->view.vup,
2771 &ang)) {
2772 if (pict_is_validpixel(p, x, y)
2773 && before_igs > 0) {
2774 actual_igs =
2775 find_near_pgs(p, x, y,
2776 max_angle,
2777 before_igs, igs,
2778 i_splitstart);
2779 if (!(actual_igs == before_igs)) {
2780 change = 1;
2781 }
2782 if (before_igs > 0) {
2783 actual_igs =
2784 pict_get_gsn(p, x, y);
2785 /* icol =
2786 setglcolor(p, x, y,
2787 actual_igs, uniform_gs, u_r, u_g , u_b);*/
2788 }
2789
2790 }
2791 }
2792 }
2793
2794 }
2795 }
2796
2797
2798 }
2799 }
2800 }
2801
2802 /* calculation of direct vertical illuminance for CGI and for disability glare, coloring glare sources*/
2803
2804 if (calcfast == 0 || calcfast == 2) {
2805 for (x = 0; x < pict_get_xsize(p); x++)
2806 for (y = 0; y < pict_get_ysize(p); y++) {
2807 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2808 if (pict_is_validpixel(p, x, y)) {
2809 if (pict_get_gsn(p, x, y) > 0) {
2810 actual_igs = pict_get_gsn(p, x, y);
2811 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));
2812 pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
2813 E_v_dir = E_v_dir + delta_E;
2814 setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2815 }
2816 }
2817 }
2818 }
2819 lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2820
2821 }
2822 /* calc of band luminance distribution if applied */
2823 if (band_calc == 1) {
2824 x_max = pict_get_xsize(p) - 1;
2825 y_max = pict_get_ysize(p) - 1;
2826 y_mid = (int)(y_max/2);
2827 for (x = 0; x <= x_max; x++)
2828 for (y = 0; y <= y_max; y++) {
2829 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2830 if (pict_is_validpixel(p, x, y)) {
2831
2832 r_actual = acos(DOT(pict_get_cached_dir(p, x, y_mid), pict_get_cached_dir(p, x, y))) ;
2833 act_lum = luminance(pict_get_color(p, x, y));
2834
2835 if ((r_actual <= band_angle) || (y == y_mid) ) {
2836 lum_band[0]=luminance(pict_get_color(p, x, y));
2837 muc_rvar_add_sample(s_band, lum_band);
2838 act_lum = luminance(pict_get_color(p, x, y));
2839 lum_band_av += pict_get_omega(p, x, y) * act_lum;
2840 omega_band += pict_get_omega(p, x, y);
2841 if (band_color == 1) {
2842 pict_get_color(p, x, y)[RED] = 0.0;
2843 pict_get_color(p, x, y)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
2844 pict_get_color(p, x, y)[BLU] = 0.0;
2845 }
2846 }
2847 }}}
2848 lum_band_av=lum_band_av/omega_band;
2849 muc_rvar_get_vx(s_band,lum_band_std);
2850 muc_rvar_get_percentile(s_band,lum_band_median,0.75);
2851 per_75_band=lum_band_median[0];
2852 muc_rvar_get_percentile(s_band,lum_band_median,0.95);
2853 per_95_band=lum_band_median[0];
2854 muc_rvar_get_median(s_band,lum_band_median);
2855 muc_rvar_get_bounding_box(s_band,bbox_band);
2856
2857 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] );
2858 /* av_lum = av_lum / omega_sum; */
2859 /* printf("average luminance of band %f \n",av_lum);*/
2860 /* printf("solid angle of band %f \n",omega_sum);*/
2861 }
2862
2863 /*printf("total number of glare sources: %i\n",igs); */
2864 lum_sources = 0;
2865 omega_sources = 0;
2866 i=0;
2867 for (x = 0; x <= igs; x++) {
2868 if (pict_get_npix(p, x) > 0) {
2869 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);
2870 lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2871 omega_sources += pict_get_av_omega(p, x);
2872 i=i+1;
2873 }}
2874 sum_glare= c2*log10(1 + sum_glare / pow(E_v, a3));
2875 if (sang == omega_sources) {
2876 lum_backg =avlum;
2877 } else {
2878 lum_backg =(sang * avlum - lum_sources) / (sang - omega_sources);
2879 }
2880
2881
2882 if (i == 0) {
2883 lum_sources =0.0;
2884 } else { lum_sources=lum_sources/omega_sources;
2885 }
2886
2887
2888
2889 if (non_cos_lb == 0) {
2890 lum_backg = lum_backg_cos;
2891 }
2892
2893 if (non_cos_lb == 2) {
2894 lum_backg = E_v / 3.1415927;
2895 }
2896
2897
2898 /* file writing NOT here
2899 if (checkfile == 1) {
2900 pict_write(p, file_out);
2901 }
2902 */
2903
2904 /* print detailed output */
2905
2906
2907
2908 /* masking */
2909
2910 if ( masking == 1 ) {
2911
2912
2913 for (x = 0; x < pict_get_xsize(p); x++)
2914 for (y = 0; y < pict_get_ysize(p); y++) {
2915 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2916 if (pict_is_validpixel(p, x, y)) {
2917 if (luminance(pict_get_color(pm, x, y))>0.1) {
2918 /* printf ("hallo %i %i %f \n",x,y,luminance(pict_get_color(pm, x, y)));*/
2919 i_mask=i_mask+1;
2920 lum_mask[0]=luminance(pict_get_color(p, x, y));
2921 /* printf ("%f \n",lum_mask[0]);*/
2922 muc_rvar_add_sample(s_mask, lum_mask);
2923 omega_mask += pict_get_omega(p, x, y);
2924 lum_mask_av += pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2925 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));
2926
2927 pict_get_color(pm, x, y)[RED] = luminance(pict_get_color(p, x, y))/179.0;
2928 pict_get_color(pm, x, y)[GRN] = luminance(pict_get_color(p, x, y))/179.0;
2929 pict_get_color(pm, x, y)[BLU] = luminance(pict_get_color(p, x, y))/179.0;
2930
2931 } else {
2932 pict_get_color(p, x, y)[RED] = 0.0 ;
2933 pict_get_color(p, x, y)[GRN] = 0.0 ;
2934 pict_get_color(p, x, y)[BLU] = 0.0 ;
2935 }
2936 }
2937 }
2938 }
2939 /* Average luminance in masking area (weighted by solid angle) */
2940 lum_mask_av=lum_mask_av/omega_mask;
2941 muc_rvar_get_vx(s_mask,lum_mask_std);
2942 muc_rvar_get_percentile(s_mask,lum_mask_median,0.75);
2943 per_75_mask=lum_mask_median[0];
2944 muc_rvar_get_percentile(s_mask,lum_mask_median,0.95);
2945 per_95_mask=lum_mask_median[0];
2946 muc_rvar_get_median(s_mask,lum_mask_median);
2947 muc_rvar_get_bounding_box(s_mask,bbox);
2948 /* PSGV only why masking of window is applied! */
2949
2950
2951 if (task_lum == 0 || lum_task == 0.0 ) {
2952 fprintf(stderr, " warning: Task area not set or task luminance=0 ! pgsv cannot be calculated (set to -99)!!\n");
2953 pgsv = -99;
2954 } else {
2955 pgsv = get_pgsv(E_v, E_v_mask, omega_mask, lum_mask_av,lum_task,avlum);
2956 }
2957
2958 pgsv_con = get_pgsv_con(E_v, E_v_mask, omega_mask, lum_mask_av,avlum);
2959 pgsv_sat =get_pgsv_sat(E_v);
2960
2961 if (detail_out == 1) {
2962
2963 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 );
2964
2965 }
2966
2967 }
2968
2969 /* zones */
2970
2971 if ( zones > 0 ) {
2972
2973 for (x = 0; x < pict_get_xsize(p); x++)
2974 for (y = 0; y < pict_get_ysize(p); y++) {
2975 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2976 if (pict_is_validpixel(p, x, y)) {
2977
2978
2979 r_actual =acos(DOT(pict_get_cached_dir(p, x, y), pict_get_cached_dir(p, x_zone, y_zone))) * 2;
2980 lum_actual = luminance(pict_get_color(p, x, y));
2981 act_gsn=pict_get_gsn(p, x, y);
2982 /* 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));*/
2983
2984 /*zone1 */
2985 if (r_actual <= angle_z1) {
2986 i_z1=i_z1+1;
2987 lum_z1[0]=lum_actual;
2988 muc_rvar_add_sample(s_z1, lum_z1);
2989 omega_z1 += pict_get_omega(p, x, y);
2990 lum_z1_av += pict_get_omega(p, x, y)* lum_actual;
2991 Ez1 += DOT(p->view.vdir, pict_get_cached_dir(p, x, y))* pict_get_omega(p, x, y)* lum_actual;
2992 setglcolor(p,x,y,1,1 , 0.66, 0.01 ,0.33);
2993 /*check if separation gsn already exist */
2994
2995 if (act_gsn > 0){
2996 if (pict_get_z1_gsn(p,act_gsn) == 0) {
2997 pict_new_gli(p);
2998 igs = igs + 1;
2999 pict_get_z1_gsn(p,act_gsn) = igs*1.000;
3000 pict_get_z1_gsn(p,igs) = -1.0;
3001 pict_get_z2_gsn(p,igs) = -1.0;
3002 splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
3003 /* 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)); */
3004 }
3005 splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
3006 /* printf("splitgs%i \n",splitgs); */
3007 split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
3008 /* move direct illuminance contribution into zone -value */
3009 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));
3010 pict_get_Eglare(p,act_gsn ) = pict_get_Eglare(p,act_gsn ) - delta_E;
3011 pict_get_Eglare(p,igs ) = pict_get_Eglare(p,igs ) + delta_E;
3012
3013
3014 }
3015 }
3016 /*zone2 */
3017
3018 if (r_actual > angle_z1 && r_actual<= angle_z2 ) {
3019
3020 i_z2=i_z2+1;
3021 lum_z2[0]=lum_actual;
3022 muc_rvar_add_sample(s_z2, lum_z2);
3023 omega_z2 += pict_get_omega(p, x, y);
3024 lum_z2_av += pict_get_omega(p, x, y)* lum_actual;
3025 Ez2 += DOT(p->view.vdir, pict_get_cached_dir(p, x, y))* pict_get_omega(p, x, y)* lum_actual;
3026 setglcolor(p,x,y,1,1 , 0.65, 0.33 ,0.02);
3027 /* 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));
3028 */ if (act_gsn > 0){
3029 if (pict_get_z2_gsn(p,act_gsn) == 0) {
3030 pict_new_gli(p);
3031 igs = igs + 1;
3032 pict_get_z2_gsn(p,act_gsn) = igs*1.000;
3033 pict_get_z1_gsn(p,igs) = -2.0;
3034 pict_get_z2_gsn(p,igs) = -2.0;
3035 /* printf ("Glare source %i gets glare source %i in zone 2 \n",act_gsn,igs ); */
3036 }
3037 splitgs=(int)(pict_get_z2_gsn(p,act_gsn));
3038 /* printf("splitgs %i \n",splitgs);*/
3039 split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
3040 /* move direct illuminance contribution into zone -value */
3041 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));
3042 pict_get_Eglare(p,act_gsn ) = pict_get_Eglare(p,act_gsn ) - delta_E;
3043 pict_get_Eglare(p,igs ) = pict_get_Eglare(p,igs ) + delta_E;
3044
3045 }
3046 }
3047
3048 }
3049 } }
3050 /* Average luminance in zones (weighted by solid angle), calculate percentiles, median min and max in zones */
3051 if (zones == 2 ) {
3052 lum_z2_av=lum_z2_av/omega_z2;
3053 muc_rvar_get_vx(s_z2,lum_z2_std);
3054 muc_rvar_get_percentile(s_z2,lum_z2_median,0.75);
3055 per_75_z2=lum_z2_median[0];
3056 muc_rvar_get_percentile(s_z2,lum_z2_median,0.95);
3057 per_95_z2=lum_z2_median[0];
3058 muc_rvar_get_median(s_z2,lum_z2_median);
3059 muc_rvar_get_bounding_box(s_z2,bbox_z2);
3060 }
3061 lum_z1_av=lum_z1_av/omega_z1;
3062 muc_rvar_get_vx(s_z1,lum_z1_std);
3063 muc_rvar_get_percentile(s_z1,lum_z1_median,0.75);
3064 per_75_z1=lum_z1_median[0];
3065 muc_rvar_get_percentile(s_z1,lum_z1_median,0.95);
3066 per_95_z1=lum_z1_median[0];
3067 muc_rvar_get_median(s_z1,lum_z1_median);
3068 muc_rvar_get_bounding_box(s_z1,bbox_z1);
3069 if (detail_out == 1) {
3070
3071 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 );
3072
3073 if (zones == 2 ) {
3074
3075 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 );
3076 } }
3077
3078 }
3079
3080 int new_gs_number[igs+1];
3081
3082 i = 0;
3083 for (x = 0; x <= igs; x++) {
3084 /* resorting glare source numbers */
3085 if (pict_get_npix(p, x) > 0) {
3086 i = i + 1;
3087 pict_get_Dglare(p,x) = i;
3088 new_gs_number[x] = i;
3089 /* printf("%i %i %f %i \n", i,x,pict_get_Dglare(p,x),new_gs_number[x] ); */
3090 }
3091 }
3092 no_glaresources=i;
3093 /*printf("%i",no_glaresources );*/
3094 /* glare sources */
3095
3096 if (output == 4 ) {
3097
3098 i=0;
3099 for (x = 0; x <= igs; x++) {
3100 if (pict_get_npix(p, x) > 0) {
3101 i = i + 1;
3102
3103 x2=pict_get_av_posx(p, x);
3104 y2=pict_get_av_posy(p, x);
3105
3106 printf("%i %f %f %f %f %.10f %f %f %f %f \n",
3107 i, pict_get_npix(p, x), pict_get_av_posx(p, x), pict_get_av_posy(p, x),
3108 pict_get_av_lum(p, x), pict_get_av_omega(p, x),
3109 get_posindex(p, pict_get_av_posx(p, x),pict_get_av_posy(p, x),posindex_2),
3110 pict_get_cached_dir(p, x2,y2)[0],pict_get_cached_dir(p, x2,y2)[1],pict_get_cached_dir(p, x2,y2)[2]);
3111 }
3112 }
3113
3114
3115 for (y = 0; y < pict_get_ysize(p); y++) {
3116 for (x = 0; x < pict_get_xsize(p); x++) {
3117
3118 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
3119 if (pict_is_validpixel(p, x, y)) {
3120 if (pict_get_gsn(p,x,y) >0 ) {
3121 i=pict_get_gsn(p,x,y);
3122 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] );
3123
3124 }
3125
3126
3127
3128 }}}}
3129
3130
3131 }
3132
3133
3134
3135
3136 if (detail_out == 1 && output < 3) {
3137 dgp = get_dgp(p, E_v, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
3138
3139 printf
3140 ("%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",
3141 i);
3142 if (i == 0) {
3143 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,
3144 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);
3145
3146 } else {
3147 i = 0;
3148 for (x = 0; x <= igs; x++) {
3149 if (pict_get_npix(p, x) > 0) {
3150 i = i + 1;
3151 pict_get_sigma(p, pict_get_av_posx(p, x),
3152 pict_get_av_posy(p, x), p->view.vdir,
3153 p->view.vup, &sigma);
3154
3155 x2=pict_get_av_posx(p, x);
3156 y2=pict_get_av_posy(p, x);
3157 teta = 180.0 / 3.1415927 * acos(DOT(p->view.vdir, pict_get_cached_dir(p, x2, y2)));
3158 Lveil_cie = 10*pict_get_Eglare(p,x)/(teta*teta+0.0000000000000001);
3159
3160 if (apply_disability == 1 && Lveil_cie <=disability_thresh) {
3161 Lveil_cie =0 ;
3162 }
3163 Lveil_cie_sum = Lveil_cie_sum + Lveil_cie;
3164 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)) ;
3165 r_contrast=r_glare/sum_glare ;
3166 r_dgp=r_glare/dgp ;
3167 if (pict_get_z1_gsn(p,x)<0) {
3168 act_gsn=(int)(-pict_get_z1_gsn(p,x));
3169 }else{
3170 act_gsn=0;
3171 }
3172 printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %i %f %f \n",
3173 i, pict_get_npix(p, x), pict_get_av_posx(p, x),
3174 pict_get_ysize(p) - pict_get_av_posy(p, x),
3175 pict_get_av_lum(p, x), pict_get_av_omega(p, x),
3176 get_posindex(p, pict_get_av_posx(p, x),
3177 pict_get_av_posy(p, x),
3178 posindex_2), lum_backg, lum_task,
3179 E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
3180 ,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 );
3181 }
3182 }
3183 }
3184 }
3185
3186 if ( output < 3) {
3187
3188 /* calculation of indicees */
3189
3190 /* check vertical illuminance range */
3191 E_v2 = E_v;
3192
3193 if (E_v2 < 100) {
3194 fprintf(stderr,
3195 "Notice: Vertical illuminance is below 100 lux !!\n");
3196 }
3197 dgp =
3198 get_dgp(p, E_v2, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
3199 /* low light correction */
3200 if (lowlight ==1) {
3201 if (E_v < 1000) {
3202 low_light_corr=1.0*exp(0.024*E_v-4)/(1+exp(0.024*E_v-4));} else {low_light_corr=1.0 ;}
3203 dgp =low_light_corr*dgp;
3204 }
3205 /* age correction */
3206
3207 if (age_corr == 1) {
3208 age_corr_factor=1.0/(1.1-0.5*age/100.0);
3209 }
3210 dgp =age_corr_factor*dgp;
3211
3212
3213 if (ext_vill == 1) {
3214 if (E_vl_ext < 100) {
3215 fprintf(stderr,
3216 "Notice: Vertical illuminance is below 100 lux !!\n");
3217 }
3218 }
3219
3220 if (calcfast == 0) {
3221 lum_a= E_v2/3.1415927;
3222 dgi = get_dgi(p, lum_backg, igs, posindex_2);
3223 ugr = get_ugr(p, lum_backg, igs, posindex_2);
3224 ugp = get_ugp (p,lum_backg, igs, posindex_2);
3225 ugp2 = get_ugp2 (p,lum_backg, igs, posindex_2);
3226 ugr_exp = get_ugr_exp (p,lum_backg_cos,lum_a, igs, posindex_2);
3227 dgi_mod = get_dgi_mod(p, lum_a, igs, posindex_2);
3228 cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
3229 dgr = get_dgr(p, avlum, igs, posindex_2);
3230 vcp = get_vcp(dgr);
3231 Lveil = get_disability(p, avlum, igs);
3232 if (no_glaresources == 0) {
3233 dgi = 0.0;
3234 ugr = 0.0;
3235 ugp = 0.0;
3236 ugp2 =0.0;
3237 ugr_exp = 0.0;
3238 dgi_mod = 0.0;
3239 cgi = 0.0;
3240 dgr = 0.0;
3241 vcp = 100.0;
3242 }
3243 }
3244 /* check dgp range */
3245 if (dgp <= 0.2) {
3246 fprintf(stderr,
3247 "Notice: Low brightness scene. dgp below 0.2! dgp might underestimate glare sources\n");
3248 }
3249 if (E_v < 380) {
3250 fprintf(stderr,
3251 "Notice: Low brightness scene. Vertical illuminance less than 380 lux! dgp might underestimate glare sources\n");
3252 }
3253
3254
3255
3256 if (output == 0) {
3257 if (detail_out == 1) {
3258 if (ext_vill == 1) {
3259 dgp_ext =
3260 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
3261 posindex_2);
3262 dgp = dgp_ext;
3263 if (E_vl_ext < 1000 && lowlight ==1) {
3264 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 ;}
3265 dgp =low_light_corr*dgp;
3266 dgp =age_corr_factor*dgp;
3267 }
3268 muc_rvar_get_median(s_noposweight,lum_nopos_median);
3269 muc_rvar_get_median(s_posweight,lum_pos_median);
3270 muc_rvar_get_median(s_posweight2,lum_pos2_median);
3271
3272
3273
3274
3275 printf
3276 ("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",
3277 dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
3278 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);
3279 } else {
3280 if (detail_out2 == 1) {
3281
3282 printf
3283 ("dgp,dgi,ugr,vcp,cgi,dgp_ext,Ev_calc,abs_max,Lveil: %f %f %f %f %f %f %f %f %f \n",
3284 dgp, dgi, ugr, vcp, cgi, dgp_ext, E_v, abs_max,
3285 Lveil);
3286
3287 } else {
3288 if (ext_vill == 1) {
3289 dgp_ext = get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,posindex_2);
3290
3291 if (E_vl_ext < 1000 && lowlight ==1) {
3292 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 ;}
3293 dgp =low_light_corr*dgp_ext;
3294 dgp =age_corr_factor*dgp;
3295 }
3296 printf("dgp,dgi,ugr,vcp,cgi,Lveil: %f %f %f %f %f %f \n",
3297 dgp, dgi, ugr, vcp, cgi, Lveil);
3298
3299 }
3300 }
3301 } else {
3302 dgp_ext =
3303 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
3304 posindex_2);
3305 dgp = dgp_ext;
3306 if (E_vl_ext < 1000 && lowlight ==1) {
3307 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 ;}
3308 dgp =low_light_corr*dgp;
3309
3310 if (calcfast == 2) {
3311
3312 lum_backg_cos=(E_vl_ext-dir_ill)/3.1415927;
3313 ugr = get_ugr(p, lum_backg_cos, igs, posindex_2);
3314 printf("%f %f \n", dgp,ugr);
3315 }else{
3316 printf("%f\n", dgp);
3317 }
3318 }
3319 }
3320
3321 }else{
3322 /* only multiimagemode
3323
3324 */
3325
3326
3327 for (j = 0; j < num_images; j++) {
3328
3329
3330 /* loop over temporal images */
3331
3332 pict_read(pcoeff,temp_image_name[j] );
3333
3334 /*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));
3335 */
3336
3337
3338 /* loop over scans
3339 empty of */
3340 for (jj=0;jj<=num_scans;jj++){
3341
3342 /*printf("scale %i %i %i %f ",j,jj,num_scans,scale_image_scans[j*(num_scans+1)+jj]); */
3343
3344
3345 /* copy luminance value into big image and remove glare sources*/
3346 for (x = 0; x < pict_get_xsize(pcoeff); x++) {
3347 for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3348 xmap=x_temp_img[j]+x;
3349 ymap=y_temp_img[j]+y;
3350 if (xmap <0) { xmap=0;}
3351 if (ymap <0) { ymap=0;}
3352
3353 pict_get_color(p, xmap, ymap)[RED] = scale_image_scans[j*(num_scans+1)+jj]*pict_get_color(pcoeff, x, y)[RED];
3354 pict_get_color(p, xmap, ymap)[GRN] = scale_image_scans[j*(num_scans+1)+jj]*pict_get_color(pcoeff, x, y)[GRN];
3355 pict_get_color(p, xmap, ymap)[BLU] = scale_image_scans[j*(num_scans+1)+jj]*pict_get_color(pcoeff, x, y)[BLU];
3356 pict_get_gsn(p, xmap, ymap) = 0;
3357 pict_get_pgs(p, xmap, ymap) = 0;
3358 }}
3359
3360
3361
3362
3363
3364 actual_igs =0;
3365
3366 /* first glare source scan: find primary glare sources */
3367 for (x = 0; x < pict_get_xsize(pcoeff); x++) {
3368 lastpixelwas_gs=0;
3369 /* lastpixelwas_peak=0; */
3370 for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3371 xmap=x_temp_img[j]+x;
3372 ymap=y_temp_img[j]+y;
3373 if (pict_get_hangle(p, xmap, ymap, p->view.vdir, p->view.vup, &ang)) {
3374 if (pict_is_validpixel(p, xmap, ymap)) {
3375 act_lum = luminance(pict_get_color(p, xmap, ymap));
3376 if (act_lum> lum_source) {
3377 if (act_lum >lum_total_max) {
3378 lum_total_max=act_lum;
3379 }
3380
3381 if (lastpixelwas_gs==0 || search_pix <= 1.0 ) {
3382 actual_igs = find_near_pgs(p, xmap, ymap, max_angle, 0, igs, 1);
3383 }
3384 if (actual_igs == 0) {
3385 igs = igs + 1;
3386 pict_new_gli(p);
3387 pict_get_Eglare(p,igs) = 0.0;
3388 /* not necessary here pict_get_lum_min(p, igs) = HUGE_VAL;
3389 pict_get_Eglare(p,igs) = 0.0;
3390 pict_get_Dglare(p,igs) = 0.0;
3391 pict_get_z1_gsn(p,igs) = 0;
3392 pict_get_z2_gsn(p,igs) = 0; */
3393 actual_igs = igs;
3394
3395 }
3396 pict_get_gsn(p, xmap, ymap) = actual_igs;
3397 pict_get_pgs(p, xmap, ymap) = 1;
3398 add_pixel_to_gs(p, xmap, ymap, actual_igs);
3399 lastpixelwas_gs=actual_igs;
3400
3401
3402
3403 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));
3404 pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
3405
3406
3407
3408
3409 } else {
3410 pict_get_pgs(p, xmap, ymap) = 0;
3411 lastpixelwas_gs=0;
3412 }
3413 }
3414 }
3415 }
3416 }
3417
3418
3419 /* here should be peak extraction */
3420 i_max=igs;
3421 r_split = max_angle / 2.0;
3422 for (i = 0; i <= i_max; i++) {
3423
3424 if (pict_get_npix(p, i) > 0) {
3425 l_max = pict_get_lum_max(p, i);
3426 i_splitstart = igs + 1;
3427 if (l_max >= limit) {
3428 for (x = 0; x < pict_get_xsize(pcoeff); x++)
3429 for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3430 xmap=x_temp_img[j]+x;
3431 ymap=y_temp_img[j]+y;
3432
3433
3434 if (pict_get_hangle
3435 (p, xmap, ymap, p->view.vdir, p->view.vup, &ang))
3436 {
3437 if (pict_is_validpixel(p, xmap, ymap)
3438 && luminance(pict_get_color(p, xmap, ymap))
3439 >= limit
3440 && pict_get_gsn(p, xmap, ymap) == i) {
3441 if (i_splitstart == (igs + 1)) {
3442 igs = igs + 1;
3443 pict_new_gli(p);
3444 pict_get_z1_gsn(p,igs) = 0;
3445 pict_get_z2_gsn(p,igs) = 0;
3446
3447 pict_get_Eglare(p,igs) = 0.0;
3448 pict_get_Dglare(p,igs) = 0.0;
3449 pict_get_lum_min(p, igs) =
3450 99999999999999.999;
3451 i_split = igs;
3452 } else {
3453 i_split =
3454 find_split(p, xmap, ymap, r_split,
3455 i_splitstart, igs);
3456 }
3457 if (i_split == 0) {
3458 igs = igs + 1;
3459 pict_new_gli(p);
3460 pict_get_z1_gsn(p,igs) = 0;
3461 pict_get_z2_gsn(p,igs) = 0;
3462
3463 pict_get_Eglare(p,igs) = 0.0;
3464 pict_get_Dglare(p,igs) = 0.0;
3465 pict_get_lum_min(p, igs) =
3466 99999999999999.999;
3467 i_split = igs;
3468 }
3469 split_pixel_from_gs(p, xmap, ymap, i_split, uniform_gs, u_r, u_g , u_b);
3470
3471 }
3472 }
3473 }
3474
3475 }
3476 change = 1;
3477 while (change == 1) {
3478 change = 0;
3479 for (x = 0; x < pict_get_xsize(pcoeff); x++)
3480 for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3481 xmap=x_temp_img[j]+x;
3482 ymap=y_temp_img[j]+y;
3483 before_igs = pict_get_gsn(p, xmap, ymap);
3484 if (before_igs >= i_splitstart) {
3485 if (pict_get_hangle
3486 (p, xmap, ymap, p->view.vdir, p->view.vup,
3487 &ang)) {
3488 if (pict_is_validpixel(p, xmap, ymap)
3489 && before_igs > 0) {
3490 actual_igs =
3491 find_near_pgs(p, xmap, ymap,
3492 max_angle,
3493 before_igs, igs,
3494 i_splitstart);
3495 if (!(actual_igs == before_igs)) {
3496 change = 1;
3497 }
3498 if (before_igs > 0) {
3499 actual_igs =
3500 pict_get_gsn(p, xmap, ymap);
3501 /* icol =
3502 setglcolor(p, x, y,
3503 actual_igs, uniform_gs, u_r, u_g , u_b);*/
3504 }
3505
3506 }
3507 }
3508 }
3509
3510 }
3511 }
3512
3513
3514 }
3515 }
3516
3517 /* end peak extraction */
3518
3519
3520 /* calculation of direct vertical illuminance for th multi-image-mode */
3521
3522
3523 for (x = 0; x < pict_get_xsize(pcoeff); x++)
3524 for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3525 xmap=x_temp_img[j]+x;
3526 ymap=y_temp_img[j]+y;
3527 if (pict_get_hangle(p, xmap, ymap, p->view.vdir, p->view.vup, &ang)) {
3528 if (pict_is_validpixel(p, xmap, ymap)) {
3529 if (pict_get_gsn(p, xmap, ymap) > 0) {
3530 actual_igs = pict_get_gsn(p, xmap, ymap);
3531 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));
3532 pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
3533 }
3534 }
3535 }
3536 }
3537
3538
3539
3540
3541
3542
3543
3544
3545 i = 0;
3546 for (x = 0; x <= igs; x++) {
3547 if (pict_get_npix(p, x) > 0) {
3548 i = i + 1;
3549 }}
3550 no_glaresources=i;
3551
3552 /*
3553
3554 sprintf(file_out, "%s%i%s","ray2/img_",j,".hdr");
3555 pict_write(p, file_out);
3556 */
3557 printf("%i ",no_glaresources);
3558 i = 0;
3559 for (x = 0; x <= igs; x++) {
3560 if (pict_get_npix(p, x) > 0) {
3561 i = i + 1;
3562 pict_get_sigma(p, pict_get_av_posx(p, x), pict_get_av_posy(p, x), p->view.vdir, p->view.vup, &sigma);
3563
3564 x2=pict_get_av_posx(p, x);
3565 y2=pict_get_av_posy(p, x);
3566 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),
3567 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) );
3568 }
3569 pict_get_npix(p, x)=0;
3570 pict_get_av_lum(p, x)=0;
3571 pict_get_av_posy(p, x)=0;
3572 pict_get_av_posx(p, x)=0;
3573 pict_get_av_omega(p, x)=0;
3574 }
3575 printf("\n");
3576
3577
3578 /* empty big image and remove glare sources*/
3579 for (x = 0; x < pict_get_xsize(pcoeff); x++) {
3580 for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3581 xmap=x_temp_img[j]+x;
3582 ymap=y_temp_img[j]+y;
3583 pict_get_color(p, xmap, ymap)[RED] = 0;
3584 pict_get_color(p, xmap, ymap)[GRN] = 0;
3585 pict_get_color(p, xmap, ymap)[BLU] = 0;
3586 pict_get_gsn(p, xmap, ymap) = 0;
3587 pict_get_pgs(p, xmap, ymap) = 0;
3588 }}
3589 igs=0;
3590
3591
3592 }
3593
3594
3595 }
3596
3597 }
3598
3599 /* end multi-image-mode */
3600
3601 /*printf("lowlight=%f\n", low_light_corr); */
3602
3603
3604 /* printf("hallo \n");
3605
3606
3607 apply_disability = 1;
3608 disability_thresh = atof(argv[++i]);
3609 coloring of disability glare sources red, remove all other colors
3610
3611
3612
3613 this function was removed because of bugs....
3614 has to be re-written from scratch....
3615
3616
3617 */
3618
3619
3620
3621
3622
3623
3624
3625
3626 /*output */
3627 if (checkfile == 1) {
3628
3629
3630 pict_write(p, file_out);
3631 }
3632
3633
3634
3635 return EXIT_SUCCESS;
3636 exit(0);
3637
3638 userr:
3639 fprintf(stderr,
3640 "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",
3641 progname);
3642 exit(1);
3643 }
3644
3645
3646
3647 #endif