ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.7
Committed: Sat Aug 12 15:11:09 2017 UTC (6 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R1
Changes since 2.6: +71 -26 lines
Log Message:
Updates from Jan for next release

File Contents

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