ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.9
Committed: Fri Nov 23 19:32:17 2018 UTC (5 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.8: +59 -14 lines
Log Message:
Bug fixes and man page updates from Jan Wienold

File Contents

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