ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.5
Committed: Tue Aug 2 16:35:01 2016 UTC (7 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +13 -6 lines
Log Message: 
Bug fix from Jan for output error spotted by David G-M

File Contents

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