ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.10
Committed: Tue Nov 27 18:08:24 2018 UTC (5 years, 4 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.9: +10 -4 lines
Log Message: 
Bug discovered by Randolph Fritz and fixed by Jan

File Contents

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