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