ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.4
Committed: Thu Jul 28 16:46:10 2016 UTC (7 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.3: +5 -2 lines
Log Message:
Updated evalglare man page and small bug fix from J. Wienold

File Contents

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