ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.10
Committed: Tue Nov 27 18:08:24 2018 UTC (5 years, 5 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +10 -4 lines
Log Message: 
Bug discovered by Randolph Fritz and fixed by Jan

File Contents

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