ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.11
Committed: Tue Jun 30 15:57:30 2020 UTC (3 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R4, rad5R3, HEAD
Changes since 2.10: +385 -28 lines
Log Message:
feat(evalglare): new version of evalglare from Jan Wienold

File Contents

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