ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/util/evalglare.c
Revision: 3.08
Committed: Wed Oct 1 17:48:15 2025 UTC (3 days, 21 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 3.07: +16 -10 lines
Log Message: 
fix(evalglare): Jan fixed issue with border masking

File Contents

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