ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/evalglare.c
Revision: 3.07
Committed: Tue Sep 16 15:43:14 2025 UTC (11 days, 15 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 3.06: +3 -2 lines
Log Message:
fix(evalglare): Compiler issue under Windows

File Contents

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