ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/src/util/evalglare.c
Revision: 3.05
Committed: Thu Sep 11 17:40:45 2025 UTC (3 weeks, 2 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.12: +358 -71 lines
Log Message: 
Numerous changes and fixes from Jan Wienold

File Contents

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