ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.6
Committed: Thu May 18 02:25:27 2017 UTC (6 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +41 -9 lines
Log Message: 
Bug fixes and other changes from Jan Wienold

File Contents

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