ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/evalglare.c
Revision: 3.06
Committed: Mon Sep 15 20:21:24 2025 UTC (2 weeks, 2 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 3.05: +1 -3 lines
Log Message:
fix(evalglare): Removed redundant progname definition

File Contents

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