ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.7
Committed: Sat Aug 12 15:11:09 2017 UTC (6 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R1
Changes since 2.6: +71 -26 lines
Log Message: 
Updates from Jan for next release

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