ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.4
Committed: Thu Jul 28 16:46:10 2016 UTC (7 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.3: +5 -2 lines
Log Message:
Updated evalglare man page and small bug fix from J. Wienold

File Contents

# User Rev Content
1 greg 2.2 #ifndef lint
2     static const char RCSid[] = "$Id$";
3     #endif
4 greg 2.3 /* EVALGLARE V1.29
5     * Evalglare Software License, Version 2.0
6 greg 2.1 *
7 greg 2.3 * Copyright (c) 1995 - 2016 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.1 #define EVALGLARE
314     #define PROGNAME "evalglare"
315 greg 2.4 #define VERSION "1.30 release 29.07.2016 by EPFL, J.Wienold"
316 greg 2.1 #define RELEASENAME PROGNAME " " VERSION
317    
318 greg 2.3
319     #include "pictool.h"
320 greg 2.1 #include "rtio.h"
321     #include <math.h>
322     #include <string.h>
323 greg 2.3 #include "platform.h"
324     #include "muc_randvar.h"
325    
326 greg 2.1 char *progname;
327    
328     /* subroutine to add a pixel to a glare source */
329     void add_pixel_to_gs(pict * p, int x, int y, int gsn)
330     {
331     double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
332     new_omega, act_lum;
333    
334    
335     pict_get_npix(p, gsn) = pict_get_npix(p, gsn) + 1;
336     old_av_posx = pict_get_av_posx(p, gsn);
337     old_av_posy = pict_get_av_posy(p, gsn);
338     old_av_lum = pict_get_av_lum(p, gsn);
339     old_omega = pict_get_av_omega(p, gsn);
340    
341     act_omega = pict_get_omega(p, x, y);
342     act_lum = luminance(pict_get_color(p, x, y));
343     new_omega = old_omega + act_omega;
344     pict_get_av_posx(p, gsn) =
345     (old_av_posx * old_omega + x * act_omega) / new_omega;
346     pict_get_av_posy(p, gsn) =
347     (old_av_posy * old_omega + y * act_omega) / new_omega;
348     if (isnan((pict_get_av_posx(p, gsn))))
349     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);
350    
351     pict_get_av_lum(p, gsn) =
352     (old_av_lum * old_omega + act_lum * act_omega) / new_omega;
353     pict_get_av_omega(p, gsn) = new_omega;
354     pict_get_gsn(p, x, y) = gsn;
355     if (act_lum < pict_get_lum_min(p, gsn)) {
356     pict_get_lum_min(p, gsn) = act_lum;
357     }
358     if (act_lum > pict_get_lum_max(p, gsn)) {
359     pict_get_lum_max(p, gsn) = act_lum;
360     }
361    
362     /* 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)); */
363    
364     }
365    
366     /* subroutine for peak extraction */
367     int
368     find_split(pict * p, int x, int y, double r, int i_split_start,
369     int i_split_end)
370     {
371     int i_find_split, x_min, x_max, y_min, y_max, ix, iy, iix, iiy, dx, dy,
372     out_r;
373     double r_actual;
374    
375     i_find_split = 0;
376    
377     x_min = 0;
378     y_min = 0;
379     x_max = pict_get_ysize(p) - 1;
380    
381     y_max = pict_get_ysize(p) - 1;
382    
383     for (iiy = 1; iiy <= 2; iiy++) {
384     dy = iiy * 2 - 3;
385     if (dy == -1) {
386     iy = y;
387     } else {
388     iy = y + 1;
389     }
390     while (iy <= y_max && iy >= y_min) {
391     out_r = 0;
392     for (iix = 1; iix <= 2; iix++) {
393     dx = iix * 2 - 3;
394     if (dx == -1) {
395     ix = x;
396     } else {
397     ix = x + 1;
398     }
399     while (ix <= x_max && ix >= x_min && iy >= y_min) {
400    
401     r_actual =
402     acos(DOT(pict_get_cached_dir(p, x, y),
403     pict_get_cached_dir(p, ix, iy))) * 2;
404     if (r_actual <= r) {
405     out_r = 1;
406     if (pict_get_gsn(p, ix, iy) >= i_split_start
407     && pict_get_gsn(p, ix, iy) <= i_split_end) {
408     i_find_split = pict_get_gsn(p, ix, iy);
409     }
410     } else {
411     ix = -99;
412     }
413     ix = ix + dx;
414     }
415     }
416     if (out_r == 0) {
417     iy = -99;
418     }
419     iy = iy + dy;
420    
421     }
422     }
423    
424    
425    
426     return i_find_split;
427     }
428    
429     /*
430     static int icomp(int* a,int* b)
431     {
432     if (*a < *b)
433     return -1;
434     if (*a > *b)
435     return 1;
436     return 0;
437     }
438    
439     */
440    
441     /* subroutine to find nearby glare source pixels */
442    
443     int
444     find_near_pgs(pict * p, int x, int y, float r, int act_gsn, int max_gsn,
445     int min_gsn)
446     {
447     int dx, dy, i_near_gs, xx, yy, x_min, x_max, y_min, y_max, ix, iy, iix,
448     iiy, old_gsn, new_gsn, find_gsn, change, stop_y_search,
449     stop_x_search;
450     double r_actual;
451     int ixm[3];
452    
453     i_near_gs = 0;
454     stop_y_search = 0;
455     stop_x_search = 0;
456     x_min = 0;
457     y_min = 0;
458     if (act_gsn == 0) {
459     x_max = x;
460     } else {
461     x_max = pict_get_xsize(p) - 1;
462     }
463    
464     y_max = pict_get_ysize(p) - 1;
465    
466     old_gsn = pict_get_gsn(p, x, y);
467     new_gsn = old_gsn;
468     change = 0;
469     if (act_gsn > 0) {
470     i_near_gs = pict_get_gsn(p, x, y);
471     }
472     for (iiy = 1; iiy <= 2; iiy++) {
473     dy = iiy * 2 - 3;
474     if (dy == -1) {
475     iy = y;
476     } else {
477     iy = y + 1;
478     }
479     ixm[1] = x;
480     ixm[2] = x;
481     stop_y_search = 0;
482    
483    
484     while (iy <= y_max && iy >= y_min) {
485     for (iix = 1; iix <= 2; iix++) {
486     dx = iix * 2 - 3;
487     ix = (ixm[1] + ixm[2]) / 2;
488     stop_x_search = 0;
489     while (ix <= x_max && ix >= x_min && stop_x_search == 0
490     && stop_y_search == 0) {
491     /* 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);*/
492     r_actual =
493     acos(DOT(pict_get_cached_dir(p, x, y),
494     pict_get_cached_dir(p, ix, iy))) * 2;
495     if (r_actual <= r) {
496     if (pict_get_gsn(p, ix, iy) > 0) {
497     if (act_gsn == 0) {
498     i_near_gs = pict_get_gsn(p, ix, iy);
499     stop_x_search = 1;
500     stop_y_search = 1;
501     } else {
502     find_gsn = pict_get_gsn(p, ix, iy);
503     if (pict_get_av_omega(p, old_gsn) <
504     pict_get_av_omega(p, find_gsn)
505     && pict_get_av_omega(p,
506     find_gsn) >
507     pict_get_av_omega(p, new_gsn)
508     && find_gsn >= min_gsn
509     && find_gsn <= max_gsn) {
510     /* other primary glare source found with larger solid angle -> add together all */
511     new_gsn = find_gsn;
512     change = 1;
513     stop_x_search = 1;
514     stop_y_search = 1;
515     }
516     }
517     }
518     } else {
519     ixm[iix] = ix - dx;
520     stop_x_search = 1;
521     }
522     ix = ix + dx;
523     }
524     }
525     iy = iy + dy;
526     }
527     }
528     if (change > 0) {
529     pict_get_av_lum(p, old_gsn) = 0.;
530     pict_get_av_omega(p, old_gsn) = 0.;
531     pict_get_npix(p, old_gsn) = 0.;
532     pict_get_lum_max(p, old_gsn) = 0.;
533    
534    
535     i_near_gs = new_gsn;
536     /* printf(" changing gs no %i",old_gsn);
537     printf(" to %i\n",new_gsn);
538     */
539     for (xx = 0; xx < pict_get_xsize(p); xx++)
540     for (yy = 0; yy < pict_get_ysize(p); yy++) {
541     if (pict_is_validpixel(p, x, y)) {
542     if (pict_get_gsn(p, xx, yy) == old_gsn) {
543     add_pixel_to_gs(p, xx, yy, new_gsn);
544     }
545     }
546     }
547     }
548    
549     return i_near_gs;
550    
551    
552     }
553    
554     /* subroutine for calculation of task luminance */
555     double get_task_lum(pict * p, int x, int y, float r, int task_color)
556     {
557     int x_min, x_max, y_min, y_max, ix, iy;
558     double r_actual, av_lum, omega_sum, act_lum;
559    
560    
561     x_max = pict_get_xsize(p) - 1;
562     y_max = pict_get_ysize(p) - 1;
563     x_min = 0;
564     y_min = 0;
565    
566    
567    
568     av_lum = 0.0;
569     omega_sum = 0.0;
570    
571     for (iy = y_min; iy <= y_max; iy++) {
572     for (ix = x_min; ix <= x_max; ix++) {
573    
574     /* if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
575     continue;*/
576     r_actual =
577     acos(DOT(pict_get_cached_dir(p, x, y),
578     pict_get_cached_dir(p, ix, iy))) * 2;
579     act_lum = luminance(pict_get_color(p, ix, iy));
580    
581     if (r_actual <= r) {
582     act_lum = luminance(pict_get_color(p, ix, iy));
583     av_lum += pict_get_omega(p, ix, iy) * act_lum;
584     omega_sum += pict_get_omega(p, ix, iy);
585     if (task_color == 1) {
586     pict_get_color(p, ix, iy)[RED] = 0.0;
587     pict_get_color(p, ix, iy)[GRN] = 0.0;
588     pict_get_color(p, ix, iy)[BLU] =
589     act_lum / WHTEFFICACY / CIE_bf;
590     }
591     }
592     }
593     }
594    
595     av_lum = av_lum / omega_sum;
596     /* printf("average luminance of task %f \n",av_lum);
597     printf("solid angle of task %f \n",omega_sum);*/
598     return av_lum;
599     }
600    
601    
602    
603    
604    
605     /* subroutine for coloring the glare sources
606     red is used only for explicit call of the subroutine with acol=0 , e.g. for disability glare
607     -1: set to grey*/
608     int setglcolor(pict * p, int x, int y, int acol, int uniform_gs, double u_r, double u_g ,double u_b)
609     {
610     int icol;
611     double pcol[3][1000], act_lum, l;
612    
613     l=u_r+u_g+u_b ;
614    
615     pcol[RED][0] = 1.0 / CIE_rf;
616 greg 2.3 pcol[GRN][0] = 0.0 / CIE_gf;
617     pcol[BLU][0] = 0.0 / CIE_bf;
618 greg 2.1
619 greg 2.3 pcol[RED][1] = 0.0 / CIE_rf;
620 greg 2.1 pcol[GRN][1] = 1.0 / CIE_gf;
621 greg 2.3 pcol[BLU][1] = 0.0 / CIE_bf;
622 greg 2.1
623 greg 2.3 pcol[RED][2] = 0.15 / CIE_rf;
624     pcol[GRN][2] = 0.15 / CIE_gf;
625     pcol[BLU][2] = 0.7 / CIE_bf;
626    
627     pcol[RED][3] = .5 / CIE_rf;
628     pcol[GRN][3] = .5 / CIE_gf;
629     pcol[BLU][3] = 0.0 / CIE_bf;
630    
631     pcol[RED][4] = .5 / CIE_rf;
632     pcol[GRN][4] = .0 / CIE_gf;
633     pcol[BLU][4] = .5 / CIE_bf ;
634    
635     pcol[RED][5] = .0 / CIE_rf;
636     pcol[GRN][5] = .5 / CIE_gf;
637     pcol[BLU][5] = .5 / CIE_bf;
638    
639     pcol[RED][6] = 0.333 / CIE_rf;
640     pcol[GRN][6] = 0.333 / CIE_gf;
641     pcol[BLU][6] = 0.333 / CIE_bf;
642 greg 2.1
643    
644     pcol[RED][999] = 1.0 / WHTEFFICACY;
645     pcol[GRN][999] = 1.0 / WHTEFFICACY;
646     pcol[BLU][999] = 1.0 / WHTEFFICACY;
647    
648    
649     pcol[RED][998] = u_r /(l* CIE_rf) ;
650     pcol[GRN][998] = u_g /(l* CIE_gf);
651     pcol[BLU][998] = u_b /(l* CIE_bf);
652 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); */
653 greg 2.1 icol = acol ;
654     if ( acol == -1) {icol=999;
655     }else{if (acol>0){icol = acol % 5 +1;
656     }};
657     if ( uniform_gs == 1) {icol=998;
658     } ;
659    
660     act_lum = luminance(pict_get_color(p, x, y));
661    
662     pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY;
663     pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY;
664     pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY;
665 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))); */
666 greg 2.1 return icol;
667     }
668    
669     /* this is the smooting subroutine */
670     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)
671     {
672     int x_min, x_max, y_min, y_max, ix, iy, icol;
673     double r_actual, omega_gs, omega_total, om;
674    
675    
676    
677     omega_gs = 0.0;
678     omega_total = 0.0;
679     x_min = x - r;
680     if (x_min < 0) {
681     x_min = 0;
682     }
683     x_max = x + r;
684     if (x_max > pict_get_xsize(p) - 1) {
685     x_max = pict_get_xsize(p) - 2;
686     }
687     y_min = y - r;
688     if (y_min < 0) {
689     y_min = 0;
690     }
691     y_max = y + r;
692     if (y_max > pict_get_ysize(p) - 1) {
693     y_max = pict_get_ysize(p) - 2;
694     }
695    
696     for (ix = x_min; ix <= x_max; ix++)
697     for (iy = y_min; iy <= y_max; iy++) {
698     r_actual = sqrt((x - ix) * (x - ix) + (y - iy) * (y - iy));
699     if (r_actual <= r) {
700     om = pict_get_omega(p, ix, iy);
701     omega_total += om;
702     if (pict_get_gsn(p, ix, iy) == act_gs
703     && pict_get_pgs(p, ix, iy) == 1) {
704     omega_gs = omega_gs + 1 * om;
705     }
706     }
707     }
708    
709    
710     if (omega_gs / omega_total > 0.2) {
711     /* add pixel to gs */
712    
713     add_pixel_to_gs(p, x, y, act_gs);
714    
715     /* color pixel of gs */
716    
717 greg 2.3 /* icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b); */
718 greg 2.1
719     }
720     }
721    
722     /* subroutine for removing a pixel from a glare source */
723     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)
724     {
725     int old_gsn, icol;
726     double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
727     new_omega, act_lum;
728    
729    
730     /* change existing gs */
731     old_gsn = pict_get_gsn(p, x, y);
732     pict_get_npix(p, old_gsn) = pict_get_npix(p, old_gsn) - 1;
733    
734     act_omega = pict_get_omega(p, x, y);
735     old_av_posx = pict_get_av_posx(p, old_gsn);
736     old_av_posy = pict_get_av_posy(p, old_gsn);
737     old_omega = pict_get_av_omega(p, old_gsn);
738    
739     new_omega = old_omega - act_omega;
740     pict_get_av_omega(p, old_gsn) = new_omega;
741     pict_get_av_posx(p, old_gsn) =
742     (old_av_posx * old_omega - x * act_omega) / new_omega;
743     pict_get_av_posy(p, old_gsn) =
744     (old_av_posy * old_omega - y * act_omega) / new_omega;
745    
746    
747     act_lum = luminance(pict_get_color(p, x, y));
748     old_av_lum = pict_get_av_lum(p, old_gsn);
749     pict_get_av_lum(p, old_gsn) =
750     (old_av_lum * old_omega - act_lum * act_omega) / new_omega;
751     /* add pixel to new gs */
752    
753     add_pixel_to_gs(p, x, y, new_gsn);
754    
755     /* color pixel of new gs */
756    
757 greg 2.3 /* icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b); */
758    
759 greg 2.1
760     }
761    
762     /* subroutine for the calculation of the position index */
763     float get_posindex(pict * p, float x, float y, int postype)
764     {
765     float posindex;
766     double teta, phi, sigma, tau, deg, d, s, r, fact;
767    
768    
769     pict_get_vangle(p, x, y, p->view.vdir, p->view.vup, &phi);
770     pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &teta);
771     pict_get_sigma(p, x, y, p->view.vdir, p->view.vup, &sigma);
772     pict_get_tau(p, x, y, p->view.vdir, p->view.vup, &tau);
773    
774     /* return (phi+teta+2*3.1415927); */
775    
776     deg = 180 / 3.1415927;
777     fact = 0.8;
778     if (phi == 0) {
779     phi = 0.00001;
780     }
781     if (sigma <= 0) {
782     sigma = -sigma;
783     }
784     if (teta == 0) {
785     teta = 0.0001;
786     }
787     tau = tau * deg;
788     sigma = sigma * deg;
789    
790 greg 2.3 if (postype == 1) {
791     /* KIM model */
792     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));
793     }else{
794     /* Guth model, equation from IES lighting handbook */
795 greg 2.1 posindex =
796     exp((35.2 - 0.31889 * tau -
797     1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 +
798     0.26667 * tau -
799     0.002963 * tau *
800     tau) / 100000 *
801     sigma * sigma);
802     /* below line of sight, using Iwata model */
803     if (phi < 0) {
804     d = 1 / tan(phi);
805     s = tan(teta) / tan(phi);
806     r = sqrt(1 / d * 1 / d + s * s / d / d);
807     if (r > 0.6)
808     fact = 1.2;
809     if (r > 3) {
810     fact = 1.2;
811     r = 3;
812     }
813    
814     posindex = 1 + fact * r;
815     }
816 greg 2.3 if (posindex > 16)
817 greg 2.1 posindex = 16;
818 greg 2.3 }
819 greg 2.1
820     return posindex;
821    
822     }
823    
824    
825    
826     double get_upper_cut_2eyes(float teta)
827     {
828     double phi;
829    
830     phi=pow(7.7458218+0.00057407915*teta-0.00021746318*teta*teta+8.5572726e-6*teta*teta*teta,2);
831    
832     return phi;
833     }
834     double get_lower_cut_2eyes(float teta)
835     {
836     double phi;
837    
838     phi=1/(-0.014699242-1.5541106e-5*teta+4.6898068e-6*teta*teta-5.1539687e-8*teta*teta*teta);
839    
840     return phi;
841     }
842    
843     double get_lower_cut_central(float teta)
844     {
845     double phi;
846    
847     phi=(68.227109-2.9524084*teta+0.046674262*teta*teta)/(1-0.042317294*teta+0.00075698419*teta*teta-6.5364285e-7*teta*teta*teta);
848    
849     if (teta>73) {
850     phi=60;
851     }
852    
853     return phi;
854     }
855    
856     /* subroutine for the cutting the total field of view */
857     void cut_view_1(pict*p)
858     {
859     int x,y;
860     double border,ang,teta,phi,phi2;
861     for(x=0;x<pict_get_xsize(p)-1;x++)
862     for(y=0;y<pict_get_ysize(p)-1;y++) {
863     if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
864     if (pict_is_validpixel(p, x, y)) {
865     pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
866     pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
867     pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
868     phi=phi*180/3.1415927;
869     phi2=phi2*180/3.1415927;
870     teta=teta*180/3.1415927;
871     if (teta<0) {
872     teta=-teta;
873     }
874     if(phi2>0){
875     border=get_upper_cut_2eyes(teta);
876    
877     if (phi>border) {
878     pict_get_color(p,x,y)[RED]=0;
879     pict_get_color(p,x,y)[GRN]=0;
880     pict_get_color(p,x,y)[BLU]=0;
881     }
882     }else{
883     border=get_lower_cut_2eyes(180-teta);
884     if (-phi<border && teta>135) {
885     pict_get_color(p,x,y)[RED]=0;
886     pict_get_color(p,x,y)[GRN]=0;
887     pict_get_color(p,x,y)[BLU]=0;
888     }
889     }
890    
891    
892    
893     }else{
894     pict_get_color(p,x,y)[RED]=0;
895     pict_get_color(p,x,y)[GRN]=0;
896     pict_get_color(p,x,y)[BLU]=0;
897    
898    
899     }
900    
901     /* printf("teta,phi,border=%f %f %f\n",teta,phi,border);*/
902     }
903     }
904    
905    
906     return;
907    
908     }
909     /* subroutine for the cutting the field of view seen by both eyes*/
910     void cut_view_2(pict*p)
911     {
912    
913     int x,y;
914     double border,ang,teta,phi,phi2;
915     for(x=0;x<pict_get_xsize(p)-1;x++)
916     for(y=0;y<pict_get_ysize(p)-1;y++) {
917     if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
918     if (pict_is_validpixel(p, x, y)) {
919     pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
920     pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
921     pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
922     pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
923     phi=phi*180/3.1415927;
924     phi2=phi2*180/3.1415927;
925     teta=teta*180/3.1415927;
926     /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
927     if (teta<0) {
928     teta=-teta;
929     }
930     if(phi2>0){
931     border=60;
932    
933     if (phi>border) {
934     pict_get_color(p,x,y)[RED]=0;
935     pict_get_color(p,x,y)[GRN]=0;
936     pict_get_color(p,x,y)[BLU]=0;
937     }
938     }else{
939     border=get_lower_cut_central(180-teta);
940     if (phi>border) {
941     pict_get_color(p,x,y)[RED]=0;
942     pict_get_color(p,x,y)[GRN]=0;
943     pict_get_color(p,x,y)[BLU]=0;
944     }
945     }
946    
947    
948    
949     }else{
950     pict_get_color(p,x,y)[RED]=0;
951     pict_get_color(p,x,y)[GRN]=0;
952     pict_get_color(p,x,y)[BLU]=0;
953    
954    
955     }
956     }
957     }
958    
959    
960     return;
961    
962     }
963    
964     /* subroutine for the cutting the field of view seen by both eyes
965     luminance is modified by position index - just for experiments - undocumented */
966     void cut_view_3(pict*p)
967     {
968    
969     int x,y;
970     double border,ang,teta,phi,phi2,lum,newlum;
971     for(x=0;x<pict_get_xsize(p)-1;x++)
972     for(y=0;y<pict_get_ysize(p)-1;y++) {
973     if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
974     if (DOT(pict_get_cached_dir(p,x,y),p->view.vdir) >= 0.0) {
975     pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
976     pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
977     pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
978     pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
979     phi=phi*180/3.1415927;
980     phi2=phi2*180/3.1415927;
981     teta=teta*180/3.1415927;
982     lum=luminance(pict_get_color(p,x,y));
983    
984     newlum=lum/get_posindex(p,x,y,0);
985    
986     pict_get_color(p,x,y)[RED]=newlum/WHTEFFICACY;
987     pict_get_color(p,x,y)[GRN]=newlum/WHTEFFICACY;
988     pict_get_color(p,x,y)[BLU]=newlum/WHTEFFICACY;
989    
990    
991    
992    
993     /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
994     if (teta<0) {
995     teta=-teta;
996     }
997     if(phi2>0){
998     border=60;
999    
1000     if (phi>border) {
1001     pict_get_color(p,x,y)[RED]=0;
1002     pict_get_color(p,x,y)[GRN]=0;
1003     pict_get_color(p,x,y)[BLU]=0;
1004     }
1005     }else{
1006     border=get_lower_cut_central(180-teta);
1007     if (phi>border) {
1008     pict_get_color(p,x,y)[RED]=0;
1009     pict_get_color(p,x,y)[GRN]=0;
1010     pict_get_color(p,x,y)[BLU]=0;
1011     }
1012     }
1013    
1014    
1015    
1016     }else{
1017     pict_get_color(p,x,y)[RED]=0;
1018     pict_get_color(p,x,y)[GRN]=0;
1019     pict_get_color(p,x,y)[BLU]=0;
1020    
1021    
1022     }
1023     }
1024     }
1025    
1026    
1027     return;
1028    
1029     }
1030    
1031 greg 2.3 /* subroutine for the calculation of the daylight glare index DGI*/
1032 greg 2.1
1033    
1034     float get_dgi(pict * p, float lum_backg, int igs, int posindex_2)
1035     {
1036     float dgi, sum_glare, omega_s;
1037     int i;
1038    
1039    
1040     sum_glare = 0;
1041     omega_s = 0;
1042    
1043     for (i = 0; i <= igs; i++) {
1044    
1045     if (pict_get_npix(p, i) > 0) {
1046     omega_s =
1047     pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1048     get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1049     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));
1050     /* printf("i,sum_glare %i %f\n",i,sum_glare); */
1051     }
1052     }
1053     dgi = 10 * log10(sum_glare);
1054    
1055     return dgi;
1056    
1057     }
1058 greg 2.3 /* subroutine for the calculation of the modified daylight glare index DGI_mod (Fisekis2003)
1059     several glare sources are taken into account and summed up, original equation has no summary
1060     */
1061    
1062    
1063     float get_dgi_mod(pict * p, float lum_a, int igs, int posindex_2)
1064     {
1065     float dgi_mod, sum_glare, omega_s;
1066     int i;
1067    
1068    
1069     sum_glare = 0;
1070     omega_s = 0;
1071    
1072     for (i = 0; i <= igs; i++) {
1073    
1074     if (pict_get_npix(p, i) > 0) {
1075     omega_s =
1076     pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1077     get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1078     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));
1079     /* printf("i,sum_glare %i %f\n",i,sum_glare); */
1080     }
1081     }
1082     dgi_mod = 10 * log10(sum_glare);
1083    
1084     return dgi_mod;
1085    
1086     }
1087 greg 2.1
1088     /* subroutine for the calculation of the daylight glare probability */
1089     double
1090     get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3,
1091     double a4, double a5, double c1, double c2, double c3,
1092     int posindex_2)
1093     {
1094     double dgp;
1095     double sum_glare;
1096     int i;
1097    
1098     sum_glare = 0;
1099     if (igs > 0) {
1100     for (i = 0; i <= igs; i++) {
1101    
1102     if (pict_get_npix(p, i) > 0) {
1103     sum_glare +=
1104     pow(pict_get_av_lum(p, i),
1105     a1) / pow(get_posindex(p, pict_get_av_posx(p, i),
1106     pict_get_av_posy(p, i),
1107     posindex_2),
1108     a4) * pow(pict_get_av_omega(p, i), a2);
1109     /* printf("i,sum_glare %i %f\n",i,sum_glare);*/
1110     }
1111     }
1112     dgp =
1113     c1 * pow(E_v,
1114     a5) + c3 + c2 * log10(1 + sum_glare / pow(E_v, a3));
1115     } else {
1116     dgp = c3 + c1 * pow(E_v, a5);
1117     }
1118    
1119     if (dgp > 1) {
1120     dgp = 1;
1121     }
1122     /* printf("dgp %f\n",dgp);*/
1123     return dgp;
1124    
1125     }
1126    
1127     /* subroutine for the calculation of the visual comfort probability */
1128 greg 2.3 float get_dgr(pict * p, double lum_a, int igs, int posindex_2)
1129 greg 2.1 {
1130 greg 2.3 float dgr;
1131     double sum_glare;
1132 greg 2.1 int i, i_glare;
1133    
1134    
1135     sum_glare = 0;
1136     i_glare = 0;
1137     for (i = 0; i <= igs; i++) {
1138    
1139     if (pict_get_npix(p, i) > 0) {
1140     i_glare = i_glare + 1;
1141     sum_glare +=
1142     (0.5 * pict_get_av_lum(p, i) *
1143     (20.4 * pict_get_av_omega(p, i) +
1144     1.52 * pow(pict_get_av_omega(p, i),
1145     0.2) - 0.075)) / (get_posindex(p,
1146     pict_get_av_posx
1147     (p, i),
1148     pict_get_av_posy
1149     (p, i),
1150     posindex_2) *
1151     pow(lum_a, 0.44));
1152    
1153     }
1154     }
1155     dgr = pow(sum_glare, pow(i_glare, -0.0914));
1156    
1157 greg 2.3
1158    
1159     return dgr;
1160    
1161     }
1162    
1163     float get_vcp(float dgr )
1164     {
1165     float vcp;
1166    
1167 greg 2.1 vcp = 50 * erf((6.347 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1168     if (dgr > 750) {
1169     vcp = 0;
1170     }
1171     if (dgr < 20) {
1172     vcp = 100;
1173     }
1174     return vcp;
1175    
1176     }
1177    
1178 greg 2.3
1179 greg 2.1 /* subroutine for the calculation of the unified glare rating */
1180     float get_ugr(pict * p, double lum_backg, int igs, int posindex_2)
1181     {
1182     float ugr;
1183     double sum_glare;
1184     int i, i_glare;
1185    
1186    
1187     sum_glare = 0;
1188     i_glare = 0;
1189     for (i = 0; i <= igs; i++) {
1190    
1191     if (pict_get_npix(p, i) > 0) {
1192     i_glare = i_glare + 1;
1193     sum_glare +=
1194     pow(pict_get_av_lum(p, i) /
1195     get_posindex(p, pict_get_av_posx(p, i),
1196     pict_get_av_posy(p, i), posindex_2),
1197     2) * pict_get_av_omega(p, i);
1198    
1199     }
1200     }
1201     ugr = 8 * log10(0.25 / lum_backg * sum_glare);
1202    
1203     return ugr;
1204    
1205     }
1206 greg 2.3 /* subroutine for the calculation of the experimental unified glare rating (Fisekis 2003)*/
1207     float get_ugr_exp(pict * p, double lum_backg, double lum_a, int igs, int posindex_2)
1208     {
1209     float ugr_exp;
1210     double sum_glare;
1211     int i, i_glare;
1212    
1213    
1214     sum_glare = 0;
1215     i_glare = 0;
1216     for (i = 0; i <= igs; i++) {
1217    
1218     if (pict_get_npix(p, i) > 0) {
1219     i_glare = i_glare + 1;
1220     sum_glare +=
1221     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);
1222    
1223     }
1224     }
1225     ugr_exp =8 * log10(lum_a) + 8 * log10(sum_glare/lum_backg);
1226    
1227     return ugr_exp;
1228    
1229     }
1230     /* subroutine for the calculation of the unified glare probability (Hirning 2014)*/
1231     float get_ugp(pict * p, double lum_backg, int igs, int posindex_2)
1232     {
1233     float ugp;
1234     double sum_glare;
1235     int i, i_glare;
1236    
1237    
1238     sum_glare = 0;
1239     i_glare = 0;
1240     for (i = 0; i <= igs; i++) {
1241    
1242     if (pict_get_npix(p, i) > 0) {
1243     i_glare = i_glare + 1;
1244     sum_glare +=
1245     pow(pict_get_av_lum(p, i) /
1246     get_posindex(p, pict_get_av_posx(p, i),
1247     pict_get_av_posy(p, i), posindex_2),
1248     2) * pict_get_av_omega(p, i);
1249    
1250     }
1251     }
1252     ugp = 0.26 * log10(0.25 / lum_backg * sum_glare);
1253    
1254     return ugp;
1255    
1256     }
1257 greg 2.1
1258    
1259     /* subroutine for the calculation of the disability glare according to Poynter */
1260     float get_disability(pict * p, double lum_backg, int igs)
1261     {
1262     float disab;
1263     double sum_glare, sigma, deg;
1264     int i, i_glare;
1265    
1266    
1267     sum_glare = 0;
1268     i_glare = 0;
1269     deg = 180 / 3.1415927;
1270     for (i = 0; i <= igs; i++) {
1271    
1272     if (pict_get_npix(p, i) > 0) {
1273     i_glare = i_glare + 1;
1274     pict_get_sigma(p, pict_get_av_posx(p, i),
1275     pict_get_av_posy(p, i), p->view.vdir,
1276     p->view.vup, &sigma);
1277    
1278     sum_glare +=
1279     pict_get_av_lum(p,
1280     i) * cos(sigma +
1281     0.00000000001) *
1282     pict_get_av_omega(p, i)
1283     / (deg * sigma + 0.00000000001);
1284     if (isnan(sum_glare)) {
1285     printf("sigma for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1286     printf("omega for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1287     printf("avlum for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1288     printf("avlum for %f %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i), sigma);
1289     }
1290    
1291     }
1292     }
1293     disab = 5 * sum_glare;
1294    
1295     return disab;
1296    
1297     }
1298    
1299    
1300    
1301    
1302    
1303    
1304     /* subroutine for the calculation of the cie glare index */
1305    
1306 greg 2.3 float get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1307 greg 2.1 {
1308     float cgi;
1309     double sum_glare;
1310     int i, i_glare;
1311    
1312    
1313     sum_glare = 0;
1314     i_glare = 0;
1315     for (i = 0; i <= igs; i++) {
1316    
1317     if (pict_get_npix(p, i) > 0) {
1318     i_glare = i_glare + 1;
1319     sum_glare +=
1320     pow(pict_get_av_lum(p, i) /
1321     get_posindex(p, pict_get_av_posx(p, i),
1322     pict_get_av_posy(p, i), posindex_2),
1323     2) * pict_get_av_omega(p, i);
1324    
1325     }
1326     }
1327     cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare);
1328    
1329     return cgi;
1330 greg 2.3 }
1331    
1332     /* 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! */
1333     float get_pgsv(double E_v, double E_mask,double omega_mask,double lum_mask_av)
1334     {
1335     float pgsv;
1336     double Lb;
1337    
1338     Lb = (E_v-E_mask)/1.414213562373;
1339    
1340     pgsv =3.2*log10(lum_mask_av)-0.64*log10(omega_mask)+(0.79*log10(omega_mask)-0.61)*log10(Lb)-8.2 ;
1341    
1342    
1343     return pgsv;
1344     }
1345    
1346     /* 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! */
1347     float get_pgsv_sat(double E_v)
1348     {
1349     float pgsv_sat;
1350    
1351     pgsv_sat =3.3-(0.57+3.3)/pow((1+E_v/1.414213562373/1250),1.7);
1352    
1353 greg 2.1
1354 greg 2.3 return pgsv_sat;
1355 greg 2.1 }
1356    
1357    
1358 greg 2.3
1359    
1360 greg 2.1 #ifdef EVALGLARE
1361    
1362    
1363     /* main program
1364     ------------------------------------------------------------------------------------------------------------------*/
1365    
1366    
1367     int main(int argc, char **argv)
1368     {
1369     #define CLINEMAX 4095 /* memory allocated for command line string */
1370     pict *p = pict_create();
1371 greg 2.3 pict *pm = pict_create();
1372     int skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, x1,y1, fill, yfillmax, yfillmin,
1373     ext_vill, set_lum_max, set_lum_max2, img_corr,x_disk,y_disk,task_color, i_splitstart,zones,act_gsn,splitgs,
1374     i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,x_zone,y_zone, i_z1, i_z2,
1375     igs, actual_igs, lastpixelwas_gs, icol, xt, yt, change,checkpixels, before_igs, sgs, splithigh,uniform_gs,x_max, y_max,y_mid,
1376     detail_out, posindex_picture, non_cos_lb, rx, ry, rmx,rmy,apply_disability,band_calc,band_color,masking,i_mask,no_glaresources,force;
1377     double lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma,om,delta_E,
1378     E_vl_ext, lum_max, new_lum_max, r_center, ugp, ugr_exp, dgi_mod,lum_a, pgsv,E_v_mask,pgsv_sat,angle_disk,dist,n_corner_px,zero_corner_px,
1379     search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,r_actual,lum_actual,
1380     omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang, angle_z1, angle_z2,per_95_band,per_75_band,pos,
1381     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,
1382     lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum,
1383     lum_mask[1],lum_mask_av,lum_mask_std[1],lum_mask_median[1],omega_mask,bbox[2],
1384     lum_band[1],lum_band_av,lum_band_std[1],lum_band_median[1],omega_band,bbox_band[2],
1385     lum_z1[1],lum_z1_av,lum_z1_std[1],lum_z1_median[1],omega_z1,bbox_z1[2],
1386     lum_z2[1],lum_z2_av,lum_z2_std[1],lum_z2_median[1],omega_z2,bbox_z2[2],
1387     lum_pos[1],lum_nopos_median[1],lum_pos_median[1],lum_pos2_median[1],lum_pos_mean,lum_pos2_mean;
1388     float lum_task, lum_thres, dgi, vcp, cgi, ugr, limit, dgr,
1389 greg 2.1 abs_max, Lveil;
1390 greg 2.3 char maskfile[500],file_out[500], file_out2[500], version[500];
1391 greg 2.1 char *cline;
1392     VIEW userview = STDVIEW;
1393     int gotuserview = 0;
1394 greg 2.3 struct muc_rvar* s_mask;
1395     s_mask = muc_rvar_create();
1396     muc_rvar_set_dim(s_mask, 1);
1397     muc_rvar_clear(s_mask);
1398     struct muc_rvar* s_band;
1399     s_band = muc_rvar_create();
1400     muc_rvar_set_dim(s_band, 1);
1401     muc_rvar_clear(s_band);
1402     struct muc_rvar* s_z1;
1403     s_z1 = muc_rvar_create();
1404     muc_rvar_set_dim(s_z1, 1);
1405     muc_rvar_clear(s_z1);
1406    
1407     struct muc_rvar* s_z2;
1408     s_z2 = muc_rvar_create();
1409     muc_rvar_set_dim(s_z2, 1);
1410     muc_rvar_clear(s_z2);
1411    
1412     struct muc_rvar* s_noposweight;
1413     s_noposweight = muc_rvar_create();
1414     muc_rvar_set_dim(s_noposweight, 1);
1415     muc_rvar_clear(s_noposweight);
1416    
1417     struct muc_rvar* s_posweight;
1418     s_posweight = muc_rvar_create();
1419     muc_rvar_set_dim(s_posweight, 1);
1420     muc_rvar_clear(s_posweight);
1421    
1422     struct muc_rvar* s_posweight2;
1423     s_posweight2 = muc_rvar_create();
1424     muc_rvar_set_dim(s_posweight2, 1);
1425     muc_rvar_clear(s_posweight2);
1426 greg 2.1
1427     /*set required user view parameters to invalid values*/
1428 greg 2.3 delta_E=0.0;
1429     no_glaresources=0;
1430     n_corner_px=0;
1431     zero_corner_px=0;
1432     force=0;
1433     dist=0.0;
1434     u_r=0.33333;
1435     u_g=0.33333;
1436     u_b=0.33333;
1437     band_angle=0;
1438     angle_z1=0;
1439     angle_z2=0;
1440     x_zone=0;
1441     y_zone=0;
1442     per_75_z2=0;
1443     per_95_z2=0;
1444     lum_pos_mean=0;
1445     lum_pos2_mean=0;
1446     lum_band_av = 0.0;
1447     omega_band = 0.0;
1448     pgsv = 0.0 ;
1449     pgsv_sat = 0.0 ;
1450     E_v_mask = 0.0;
1451     lum_z1_av = 0.0;
1452     omega_z1 = 0.0;
1453     lum_z2_av = 0.0;
1454     omega_z2 = 0.0 ;
1455     i_z1 = 0;
1456     i_z2 = 0;
1457     zones = 0;
1458     lum_z2_av = 0.0;
1459     uniform_gs = 0;
1460 greg 2.1 apply_disability = 0;
1461     disability_thresh = 0;
1462     Lveil_cie_sum=0.0;
1463     skip_second_scan=0;
1464     lum_total_max=0.0;
1465     calcfast=0;
1466     age_corr_factor = 1.0;
1467     age_corr = 0;
1468     age = 20.0;
1469     userview.horiz = 0;
1470     userview.vert = 0;
1471     userview.type = 0;
1472     dgp_ext = 0;
1473     E_vl_ext = 0.0;
1474     new_lum_max = 0.0;
1475     lum_max = 0.0;
1476     omegat = 0.0;
1477     yt = 0;
1478     xt = 0;
1479 greg 2.3 x_disk=0;
1480     y_disk=0;
1481     angle_disk=0;
1482 greg 2.1 yfillmin = 0;
1483     yfillmax = 0;
1484     cut_view = 0;
1485     cut_view_type = 0;
1486     setvalue = 2e09;
1487     omega_cos_contr = 0.0;
1488     lum_ideal = 0.0;
1489     max_angle = 0.2;
1490     lum_thres = 5.0;
1491     task_lum = 0;
1492     sgs = 0;
1493     splithigh = 1;
1494     detail_out = 0;
1495     detail_out2 = 0;
1496     posindex_picture = 0;
1497     checkfile = 0;
1498     ext_vill = 0;
1499     fill = 0;
1500     a1 = 2.0;
1501     a2 = 1.0;
1502     a3 = 1.87;
1503     a4 = 2.0;
1504     a5 = 1.0;
1505     c1 = 5.87e-05;
1506     c2 = 0.092;
1507     c3 = 0.159;
1508     non_cos_lb = 1;
1509     posindex_2 = 0;
1510     task_color = 0;
1511     limit = 50000.0;
1512     set_lum_max = 0;
1513     set_lum_max2 = 0;
1514 greg 2.3 img_corr=0;
1515 greg 2.1 abs_max = 0;
1516     progname = argv[0];
1517     E_v_contr = 0.0;
1518 greg 2.3 strcpy(version, "1.19 release 09.12.2015 by J.Wienold");
1519 greg 2.1 low_light_corr=1.0;
1520     output = 0;
1521     calc_vill = 0;
1522     band_avlum = -99;
1523     band_calc = 0;
1524 greg 2.3 masking = 0;
1525     lum_mask[0]=0.0;
1526     lum_mask_av=0.0;
1527     omega_mask=0.0;
1528     i_mask=0;
1529     actual_igs=0;
1530 greg 2.1 /* command line for output picture*/
1531    
1532     cline = (char *) malloc(CLINEMAX+1);
1533     cline[0] = '\0';
1534     for (i = 0; i < argc; i++) {
1535     /* fprintf(stderr, "%d %d \n",i,strlen(argv[i]));*/
1536     if (strlen(cline)+strlen(argv[i])+strlen(RELEASENAME)+2 >=CLINEMAX) {
1537     exit (-1);
1538     }
1539     else {
1540     strcat(cline, argv[i]);
1541     strcat(cline, " ");
1542     }
1543     }
1544     strcat(cline, "\n");
1545     strcat(cline, RELEASENAME);
1546     strcat(cline, "\n");
1547    
1548    
1549     /* program options */
1550    
1551     for (i = 1; i < argc && argv[i][0] == '-'; i++) {
1552     /* expand arguments */
1553     while ((rval = expandarg(&argc, &argv, i)) > 0);
1554     if (rval < 0) {
1555     fprintf(stderr, "%s: cannot expand '%s'", argv[0], argv[i]);
1556     exit(1);
1557     }
1558     rval = getviewopt(&userview, argc - i, argv + i);
1559     if (rval >= 0) {
1560     i += rval;
1561     gotuserview++;
1562     continue;
1563     }
1564     switch (argv[i][1]) {
1565     case 'a':
1566     age = atof(argv[++i]);
1567     age_corr = 1;
1568     printf("age factor not supported any more \n");
1569     exit(1);
1570     break;
1571 greg 2.3 case 'A':
1572     masking = 1;
1573     detail_out = 1;
1574     strcpy(maskfile, argv[++i]);
1575     pict_read(pm,maskfile );
1576    
1577     break;
1578 greg 2.1 case 'b':
1579     lum_thres = atof(argv[++i]);
1580     break;
1581     case 'c':
1582     checkfile = 1;
1583     strcpy(file_out, argv[++i]);
1584     break;
1585     case 'u':
1586     uniform_gs = 1;
1587     u_r = atof(argv[++i]);
1588     u_g = atof(argv[++i]);
1589     u_b = atof(argv[++i]);
1590     break;
1591     case 'r':
1592     max_angle = atof(argv[++i]);
1593     break;
1594     case 's':
1595     sgs = 1;
1596     break;
1597 greg 2.3 case 'f':
1598     force = 1;
1599     break;
1600 greg 2.1 case 'k':
1601     apply_disability = 1;
1602     disability_thresh = atof(argv[++i]);
1603     break;
1604    
1605     case 'p':
1606     posindex_picture = 1;
1607     break;
1608 greg 2.3 case 'P':
1609     posindex_2 = atoi(argv[++i]);
1610     break;
1611    
1612 greg 2.1
1613     case 'y':
1614     splithigh = 1;
1615     break;
1616     case 'x':
1617     splithigh = 0;
1618     break;
1619     case 'Y':
1620     splithigh = 1;
1621     limit = atof(argv[++i]);
1622     break;
1623    
1624     case 'i':
1625     ext_vill = 1;
1626     E_vl_ext = atof(argv[++i]);
1627     break;
1628     case 'I':
1629     ext_vill = 1;
1630     fill = 1;
1631     E_vl_ext = atof(argv[++i]);
1632     yfillmax = atoi(argv[++i]);
1633     yfillmin = atoi(argv[++i]);
1634     break;
1635     case 'd':
1636     detail_out = 1;
1637     break;
1638     case 'D':
1639     detail_out2 = 1;
1640     break;
1641     case 'm':
1642 greg 2.3 img_corr = 1;
1643 greg 2.1 set_lum_max = 1;
1644     lum_max = atof(argv[++i]);
1645     new_lum_max = atof(argv[++i]);
1646     strcpy(file_out2, argv[++i]);
1647     /* printf("max lum set to %f\n",new_lum_max);*/
1648     break;
1649    
1650     case 'M':
1651 greg 2.3 img_corr = 1;
1652 greg 2.1 set_lum_max2 = 1;
1653     lum_max = atof(argv[++i]);
1654     E_vl_ext = atof(argv[++i]);
1655     strcpy(file_out2, argv[++i]);
1656     /* printf("max lum set to %f\n",new_lum_max);*/
1657     break;
1658 greg 2.3 case 'N':
1659     img_corr = 1;
1660     set_lum_max2 = 2;
1661     x_disk = atoi(argv[++i]);
1662     y_disk = atoi(argv[++i]);
1663     angle_disk = atof(argv[++i]);
1664     E_vl_ext = atof(argv[++i]);
1665     strcpy(file_out2, argv[++i]);
1666     /* printf("max lum set to %f\n",new_lum_max);*/
1667     break;
1668    
1669 greg 2.1 case 'n':
1670     non_cos_lb = 0;
1671     break;
1672    
1673     case 't':
1674     task_lum = 1;
1675     xt = atoi(argv[++i]);
1676     yt = atoi(argv[++i]);
1677     omegat = atof(argv[++i]);
1678     task_color = 0;
1679     break;
1680     case 'T':
1681     task_lum = 1;
1682     xt = atoi(argv[++i]);
1683     yt = atoi(argv[++i]);
1684     /* omegat= DEG2RAD(atof(argv[++i]));*/
1685     omegat = atof(argv[++i]);
1686     task_color = 1;
1687     break;
1688 greg 2.3 case 'l':
1689     zones = 1;
1690     x_zone = atoi(argv[++i]);
1691     y_zone = atoi(argv[++i]);
1692     angle_z1 = atof(argv[++i]);
1693     angle_z2 = -1;
1694     break;
1695     case 'L':
1696     zones = 2;
1697     x_zone = atoi(argv[++i]);
1698     y_zone = atoi(argv[++i]);
1699     angle_z1 = atof(argv[++i]);
1700     angle_z2 = atof(argv[++i]);
1701     break;
1702    
1703    
1704 greg 2.1 case 'B':
1705     band_calc = 1;
1706     band_angle = atof(argv[++i]);
1707     band_color = 1;
1708     break;
1709    
1710    
1711     case 'w':
1712     a1 = atof(argv[++i]);
1713     a2 = atof(argv[++i]);
1714     a3 = atof(argv[++i]);
1715     a4 = atof(argv[++i]);
1716     a5 = atof(argv[++i]);
1717     c1 = atof(argv[++i]);
1718     c2 = atof(argv[++i]);
1719     c3 = atof(argv[++i]);
1720     break;
1721     case 'V':
1722     calc_vill = 1;
1723     break;
1724     case 'G':
1725     cut_view = 1;
1726     cut_view_type= atof(argv[++i]);
1727     break;
1728     case 'g':
1729     cut_view = 2;
1730     cut_view_type= atof(argv[++i]);
1731     break;
1732    
1733     /*case 'v':
1734     printf("evalglare %s \n",version);
1735     exit(1); */
1736    
1737     case '1':
1738     output = 1;
1739     break;
1740    
1741     case 'v':
1742     if (argv[i][2] == '\0') {
1743     printf("%s \n",RELEASENAME);
1744     exit(1);
1745     }
1746     if (argv[i][2] != 'f')
1747     goto userr;
1748     rval = viewfile(argv[++i], &userview, NULL);
1749     if (rval < 0) {
1750     fprintf(stderr,
1751     "%s: cannot open view file \"%s\"\n",
1752     progname, argv[i]);
1753     exit(1);
1754     } else if (rval == 0) {
1755     fprintf(stderr,
1756     "%s: bad view file \"%s\"\n", progname, argv[i]);
1757     exit(1);
1758     } else {
1759     gotuserview++;
1760     }
1761     break;
1762    
1763    
1764     default:
1765     goto userr;
1766     }
1767     }
1768    
1769     /*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/
1770    
1771     if (output == 1 && ext_vill == 1) {
1772     calcfast=1;
1773     }
1774 greg 2.3 /*masking and zoning cannot be applied at the same time*/
1775    
1776     if (masking ==1 && zones >0) {
1777     fprintf(stderr, " masking and zoning cannot be activated at the same time!\n");
1778     exit(1);
1779     }
1780 greg 2.1
1781     /* read picture file */
1782     if (i == argc) {
1783     SET_FILE_BINARY(stdin);
1784     FILE *fp = fdopen(fileno(stdin), "rb");
1785     if (!(fp)) {
1786     fprintf(stderr, "fdopen on stdin failed\n");
1787     return EXIT_FAILURE;
1788     }
1789     if (!(pict_read_fp(p, fp)))
1790     return EXIT_FAILURE;;
1791     } else {
1792     if (!pict_read(p, argv[i]))
1793     return EXIT_FAILURE;
1794     }
1795     if (gotuserview) {
1796     if (p->valid_view)
1797     fprintf(stderr,
1798     "warning: overriding image view by commandline argument\n");
1799     if ((userview.horiz == 0) || (userview.vert == 0) || (userview.type == 0)) {
1800     fprintf(stderr, "error: if specified, a view must at least contain -vt -vv and -vh\n");
1801     return EXIT_FAILURE;
1802     }
1803     p->view = userview;
1804     if (!(pict_update_view(p))) {
1805     fprintf(stderr, "error: invalid view specified");
1806     return EXIT_FAILURE;
1807     }
1808     pict_update_evalglare_caches(p);
1809     } else if (!(p->valid_view)) {
1810     fprintf(stderr, "error: no valid view specified\n");
1811     return EXIT_FAILURE;
1812     }
1813    
1814 greg 2.3
1815    
1816 greg 2.1 /* fprintf(stderr,"Picture read!");*/
1817    
1818     /* command line for output picture*/
1819    
1820     pict_set_comment(p, cline);
1821    
1822 greg 2.3 /* several checks */
1823 greg 2.1 if (p->view.type == VT_PAR) {
1824 greg 2.3 fprintf(stderr, "error: wrong view type! must not be parallel ! \n");
1825 greg 2.1 exit(1);
1826     }
1827    
1828 greg 2.3 if (p->view.type == VT_PER) {
1829     fprintf(stderr, "warning: image has perspective view type specified in header ! \n");
1830     }
1831    
1832     if (masking == 1) {
1833    
1834     if (!pict_get_xsize(p)==pict_get_xsize(pm) || !pict_get_ysize(p)==pict_get_ysize(pm)) {
1835     fprintf(stderr, "error: masking image has other resolution than main image ! \n");
1836     fprintf(stderr, "size must be identical \n");
1837     printf("resolution main image : %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
1838     printf("resolution masking image : %dx%d\n",pict_get_xsize(pm),pict_get_ysize(pm));
1839     exit(1);
1840    
1841     }
1842    
1843     }
1844     /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1845 greg 2.1
1846     /* Check size of search radius */
1847     rmx = (pict_get_xsize(p) / 2);
1848     rmy = (pict_get_ysize(p) / 2);
1849     rx = (pict_get_xsize(p) / 2 + 10);
1850     ry = (pict_get_ysize(p) / 2 + 10);
1851     r_center =
1852     acos(DOT(pict_get_cached_dir(p, rmx, rmy),
1853     pict_get_cached_dir(p, rx, ry))) * 2 / 10;
1854     search_pix = max_angle / r_center;
1855     if (search_pix < 1.0) {
1856     fprintf(stderr,
1857     "warning: search radius less than 1 pixel! deactivating smoothing and peak extraction...\n");
1858     splithigh = 0;
1859     sgs = 0;
1860    
1861     } else {
1862     if (search_pix < 3.0) {
1863     fprintf(stderr,
1864     "warning: search radius less than 3 pixels! -> %f \n",
1865     search_pix);
1866    
1867     }
1868     }
1869    
1870     /* Check task position */
1871    
1872     if (task_lum == 1) {
1873     if (xt >= pict_get_xsize(p) || yt >= pict_get_ysize(p) || xt < 0
1874     || yt < 0) {
1875     fprintf(stderr,
1876     "error: task position outside picture!! exit...");
1877     exit(1);
1878     }
1879     }
1880    
1881    
1882     /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1883    
1884     sang = 0.0;
1885     E_v = 0.0;
1886     E_v_dir = 0.0;
1887     avlum = 0.0;
1888     pict_new_gli(p);
1889     igs = 0;
1890 greg 2.3 pict_get_z1_gsn(p,igs) = 0;
1891     pict_get_z2_gsn(p,igs) = 0;
1892    
1893 greg 2.1
1894     /* cut out GUTH field of view and exit without glare evaluation */
1895     if (cut_view==2) {
1896     if (cut_view_type==1) {
1897     cut_view_1(p);
1898     }else{
1899    
1900     if (cut_view_type==2) {
1901     cut_view_2(p);
1902     }else{
1903     if (cut_view_type==3) {
1904     fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
1905     cut_view_3(p);
1906     }else{
1907     fprintf(stderr,"error: no valid option for view cutting!!");
1908     exit(1);
1909     }
1910     }}
1911     pict_write(p, file_out);
1912     exit(1); }
1913    
1914    
1915    
1916    
1917    
1918    
1919     /* write positionindex into checkfile and exit, activated by option -p */
1920     if (posindex_picture == 1) {
1921     for (x = 0; x < pict_get_xsize(p); x++)
1922     for (y = 0; y < pict_get_ysize(p); y++) {
1923     if (pict_get_hangle
1924     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1925     if (pict_is_validpixel(p, x, y)) {
1926     lum =
1927     get_posindex(p, x, y,
1928     posindex_2) / WHTEFFICACY;
1929    
1930     pict_get_color(p, x, y)[RED] = lum;
1931     pict_get_color(p, x, y)[GRN] = lum;
1932     pict_get_color(p, x, y)[BLU] = lum;
1933     lum_task = luminance(pict_get_color(p, x, y));
1934     /*printf("x,y,posindex=%i %i %f %f\n",x,y,lum*WHTEFFICACY,lum_task);*/
1935     }
1936     }
1937     }
1938     pict_write(p, file_out);
1939     exit(1);
1940     }
1941    
1942    
1943     /* calculation of vertical illuminance, average luminance, in case of filling activated-> fill picture */
1944     /* fill, if necessary from 0 to yfillmin */
1945    
1946     if (fill == 1) {
1947     for (x = 0; x < pict_get_xsize(p); x++)
1948     for (y = yfillmin; y > 0; y = y - 1) {
1949     y1 = y + 1;
1950     lum = luminance(pict_get_color(p, x, y1));
1951     pict_get_color(p, x, y)[RED] = lum / 179.0;
1952     pict_get_color(p, x, y)[GRN] = lum / 179.0;
1953     pict_get_color(p, x, y)[BLU] = lum / 179.0;
1954     }
1955     }
1956    
1957     if (calcfast == 0) {
1958     for (x = 0; x < pict_get_xsize(p); x++)
1959     for (y = 0; y < pict_get_ysize(p); y++) {
1960 greg 2.3 lum = luminance(pict_get_color(p, x, y));
1961     dist=sqrt((x-rmx)*(x-rmx)+(y-rmy)*(y-rmy));
1962     if (dist>((rmx+rmy)/2)) {
1963     n_corner_px=n_corner_px+1;
1964     if (lum<7.0) {
1965     zero_corner_px=zero_corner_px+1;
1966     }
1967     }
1968    
1969    
1970 greg 2.1 if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
1971     if (pict_is_validpixel(p, x, y)) {
1972     if (fill == 1 && y >= yfillmax) {
1973     y1 = y - 1;
1974     lum = luminance(pict_get_color(p, x, y1));
1975     pict_get_color(p, x, y)[RED] = lum / 179.0;
1976     pict_get_color(p, x, y)[GRN] = lum / 179.0;
1977     pict_get_color(p, x, y)[BLU] = lum / 179.0;
1978     }
1979    
1980     if (lum > abs_max) {
1981     abs_max = lum;
1982     }
1983     /* set luminance restriction, if -m is set */
1984 greg 2.3 if (img_corr == 1 ) {
1985     if (set_lum_max == 1 && lum >= lum_max) {
1986     pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
1987     pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
1988     pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
1989     /* printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
1990     lum = luminance(pict_get_color(p, x, y));
1991     }
1992     if (set_lum_max2 == 1 && lum >= lum_max) {
1993    
1994     E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
1995     omega_cos_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
1996     }
1997     if (set_lum_max2 == 2 ) {
1998     r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
1999     if (x_disk == x && y_disk==y ) r_actual=0.0;
2000    
2001     act_lum = luminance(pict_get_color(p, x, y));
2002    
2003     if (r_actual <= angle_disk) {
2004     E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2005     omega_cos_contr += DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2006     }
2007 greg 2.1
2008 greg 2.3
2009    
2010     }
2011 greg 2.1 }
2012 greg 2.3 om = pict_get_omega(p, x, y);
2013     sang += om;
2014     E_v += DOT(p->view.vdir, pict_get_cached_dir(p, x,y)) * om *lum;
2015     avlum += om * lum;
2016     pos=get_posindex(p, x, y,posindex_2);
2017    
2018     lum_pos[0]=lum;
2019     muc_rvar_add_sample(s_noposweight,lum_pos);
2020     lum_pos[0]=lum/pos;
2021     lum_pos_mean+=lum_pos[0]*om;
2022     muc_rvar_add_sample(s_posweight, lum_pos);
2023     lum_pos[0]=lum_pos[0]/pos;
2024     lum_pos2_mean+=lum_pos[0]*om;
2025     muc_rvar_add_sample(s_posweight2,lum_pos);
2026 greg 2.1
2027    
2028    
2029     } else {
2030     pict_get_color(p, x, y)[RED] = 0.0;
2031     pict_get_color(p, x, y)[GRN] = 0.0;
2032     pict_get_color(p, x, y)[BLU] = 0.0;
2033    
2034     }
2035 greg 2.3 }else {
2036     pict_get_color(p, x, y)[RED] = 0.0;
2037     pict_get_color(p, x, y)[GRN] = 0.0;
2038     pict_get_color(p, x, y)[BLU] = 0.0;
2039    
2040     }
2041 greg 2.1 }
2042    
2043 greg 2.3 /* check if image has black corners AND a perspective view */
2044    
2045     if (zero_corner_px/n_corner_px > 0.70 && (p->view.type == VT_PER) ) {
2046     printf (" corner pixels are to %f %% black! \n",zero_corner_px/n_corner_px*100);
2047     printf("error! The image has black corners AND header contains a perspective view type definition !!!\n") ;
2048    
2049     if (force==0) {
2050     printf("stopping...!!!!\n") ;
2051    
2052     exit(1);
2053     }
2054     }
2055    
2056    
2057    
2058    
2059     lum_pos_mean= lum_pos_mean/sang;
2060     lum_pos2_mean= lum_pos2_mean/sang;
2061    
2062     if (set_lum_max2 >= 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) {
2063 greg 2.1
2064     lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
2065 greg 2.3 if (set_lum_max2 == 2 && lum_ideal >= 2e9) {
2066     printf("warning! luminance of replacement pixels would be larger than 2e9 cd/m2. Value set to 2e9cd/m2!\n") ;
2067     }
2068    
2069 greg 2.1 if (lum_ideal > 0 && lum_ideal < setvalue) {
2070     setvalue = lum_ideal;
2071     }
2072 greg 2.3 printf("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f %f %f %f %f\n",
2073 greg 2.1 lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
2074    
2075 greg 2.3
2076 greg 2.1 for (x = 0; x < pict_get_xsize(p); x++)
2077     for (y = 0; y < pict_get_ysize(p); y++) {
2078     if (pict_get_hangle
2079     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2080     if (pict_is_validpixel(p, x, y)) {
2081     lum = luminance(pict_get_color(p, x, y));
2082 greg 2.3
2083    
2084 greg 2.1 if (set_lum_max2 == 1 && lum >= lum_max) {
2085    
2086     pict_get_color(p, x, y)[RED] =
2087     setvalue / 179.0;
2088     pict_get_color(p, x, y)[GRN] =
2089     setvalue / 179.0;
2090     pict_get_color(p, x, y)[BLU] =
2091     setvalue / 179.0;
2092    
2093 greg 2.3 }else{ if(set_lum_max2 ==2 ) {
2094     r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2095     if (x_disk == x && y_disk==y ) r_actual=0.0;
2096    
2097     if (r_actual <= angle_disk) {
2098    
2099     pict_get_color(p, x, y)[RED] =
2100     setvalue / 179.0;
2101     pict_get_color(p, x, y)[GRN] =
2102     setvalue / 179.0;
2103     pict_get_color(p, x, y)[BLU] =
2104     setvalue / 179.0;
2105    
2106     }
2107    
2108    
2109     }
2110 greg 2.1 }
2111     }
2112     }
2113     }
2114 greg 2.3
2115 greg 2.1
2116     pict_write(p, file_out2);
2117 greg 2.3 exit(1);
2118 greg 2.1 }
2119     if (set_lum_max == 1) {
2120     pict_write(p, file_out2);
2121    
2122     }
2123    
2124     if (calc_vill == 1) {
2125     printf("%f\n", E_v);
2126     exit(1);
2127     }
2128     }else{
2129     /* in fast calculation mode: ev=ev_ext and sang=2*pi */
2130     sang=2*3.14159265359;
2131     lum_task =E_vl_ext/sang;
2132     E_v=E_vl_ext;
2133     /* printf("calc fast!! %f %f\n", sang,lum_task);*/
2134    
2135    
2136     }
2137    
2138     /* cut out GUTH field of view for glare evaluation */
2139     if (cut_view==1) {
2140     if (cut_view_type==1) {
2141     cut_view_1(p);
2142     }else{
2143    
2144     if (cut_view_type==2) {
2145     cut_view_2(p);
2146     }else{
2147     if (cut_view_type==3) {
2148     fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
2149     cut_view_3(p);
2150     }else{
2151     fprintf(stderr,"error: no valid option for view cutting!!");
2152     exit(1);
2153     }
2154     }}
2155     }
2156    
2157     if (calcfast == 0) {
2158     avlum = avlum / sang;
2159     lum_task = avlum;
2160     }
2161     /* if (ext_vill==1) {
2162     E_v=E_vl_ext;
2163     avlum=E_v/3.1415927;
2164     } */
2165    
2166     if (task_lum == 1) {
2167     lum_task = get_task_lum(p, xt, yt, omegat, task_color);
2168     }
2169     lum_source = lum_thres * lum_task;
2170     /* printf("Task Luminance=%f\n",lum_task);
2171     pict_write(p,"task.pic");*/
2172    
2173     if (lum_thres > 100) {
2174     lum_source = lum_thres;
2175     }
2176    
2177 greg 2.3 /* printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2178 greg 2.1
2179     /* first glare source scan: find primary glare sources */
2180 greg 2.3 for (x = 0; x < pict_get_xsize(p); x++) {
2181     /* lastpixelwas_gs=0; */
2182 greg 2.1 for (y = 0; y < pict_get_ysize(p); y++) {
2183     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2184     if (pict_is_validpixel(p, x, y)) {
2185     act_lum = luminance(pict_get_color(p, x, y));
2186     if (act_lum > lum_source) {
2187     if (act_lum >lum_total_max) {
2188     lum_total_max=act_lum;
2189     }
2190 greg 2.3 /* speed improvement first scan: when last pixel was glare source, then use this glare source number instead of searching.
2191     has been deactivated with v1.25
2192    
2193     if (lastpixelwas_gs==0 || search_pix <= 1.0 ) { */
2194     actual_igs = find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2195     /* }*/
2196 greg 2.1 if (actual_igs == 0) {
2197     igs = igs + 1;
2198     pict_new_gli(p);
2199     pict_get_lum_min(p, igs) = HUGE_VAL;
2200     pict_get_Eglare(p,igs) = 0.0;
2201     pict_get_Dglare(p,igs) = 0.0;
2202 greg 2.3 pict_get_z1_gsn(p,igs) = 0;
2203     pict_get_z2_gsn(p,igs) = 0;
2204 greg 2.1 actual_igs = igs;
2205 greg 2.3
2206 greg 2.1 }
2207 greg 2.3 /* no coloring any more here icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2208 greg 2.1 pict_get_gsn(p, x, y) = actual_igs;
2209     pict_get_pgs(p, x, y) = 1;
2210     add_pixel_to_gs(p, x, y, actual_igs);
2211 greg 2.3 /* lastpixelwas_gs=actual_igs; */
2212 greg 2.1
2213     } else {
2214     pict_get_pgs(p, x, y) = 0;
2215 greg 2.3 /* lastpixelwas_gs=0;*/
2216 greg 2.1 }
2217     }
2218     }
2219 greg 2.3 }
2220     }
2221 greg 2.1 /* pict_write(p,"firstscan.pic"); */
2222    
2223 greg 2.3
2224     if (calcfast == 1 || search_pix <= 1.0) {
2225 greg 2.1 skip_second_scan=1;
2226     }
2227    
2228     /* second glare source scan: combine glare sources facing each other */
2229     change = 1;
2230 greg 2.3 i = 0;
2231 greg 2.1 while (change == 1 && skip_second_scan==0) {
2232     change = 0;
2233 greg 2.3 i = i+1;
2234     for (x = 0; x < pict_get_xsize(p); x++) {
2235 greg 2.1 for (y = 0; y < pict_get_ysize(p); y++) {
2236     if (pict_get_hangle
2237     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2238 greg 2.3 checkpixels=1;
2239     before_igs = pict_get_gsn(p, x, y);
2240    
2241     /* speed improvement: search for other nearby glare sources only, when at least one adjacend pixel has another glare source number */
2242     if (x >1 && x<pict_get_xsize(p)-2 && y >1 && y<pict_get_ysize(p)-2 ) {
2243     x1=x-1;
2244     x2=x+1;
2245     y1=y-1;
2246     y2=y+1;
2247     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) ) {
2248     checkpixels = 0;
2249     actual_igs = before_igs;
2250     } }
2251    
2252    
2253     if (pict_is_validpixel(p, x, y) && before_igs > 0 && checkpixels==1) {
2254 greg 2.1 actual_igs =
2255     find_near_pgs(p, x, y, max_angle, before_igs,
2256     igs, 1);
2257     if (!(actual_igs == before_igs)) {
2258     change = 1;
2259     }
2260     if (before_igs > 0) {
2261     actual_igs = pict_get_gsn(p, x, y);
2262 greg 2.3 /* icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2263 greg 2.1 }
2264    
2265     }
2266     }
2267     }
2268     /* pict_write(p,"secondscan.pic");*/
2269 greg 2.3 }}
2270 greg 2.1
2271     /*smoothing: add secondary glare sources */
2272     i_max = igs;
2273     if (sgs == 1) {
2274    
2275     /* simplified search radius: constant for entire picture, alway circle not an angle due to performance */
2276     x = (pict_get_xsize(p) / 2);
2277     y = (pict_get_ysize(p) / 2);
2278     rx = (pict_get_xsize(p) / 2 + 10);
2279     ry = (pict_get_ysize(p) / 2 + 10);
2280     r_center =
2281     acos(DOT(pict_get_cached_dir(p, x, y),
2282     pict_get_cached_dir(p, rx, ry))) * 2 / 10;
2283     search_pix = max_angle / r_center / 1.75;
2284    
2285     for (x = 0; x < pict_get_xsize(p); x++) {
2286     for (y = 0; y < pict_get_ysize(p); y++) {
2287     if (pict_get_hangle
2288     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2289     if (pict_is_validpixel(p, x, y)
2290     && pict_get_gsn(p, x, y) == 0) {
2291     for (i = 1; i <= i_max; i++) {
2292    
2293     if (pict_get_npix(p, i) > 0) {
2294     add_secondary_gs(p, x, y, search_pix, i, uniform_gs, u_r, u_g , u_b);
2295     }
2296     }
2297    
2298     }
2299     }
2300     }
2301     }
2302    
2303     }
2304    
2305     /* extract extremes from glare sources to extra glare source */
2306     if (splithigh == 1 && lum_total_max>limit) {
2307    
2308     r_split = max_angle / 2.0;
2309     for (i = 0; i <= i_max; i++) {
2310    
2311     if (pict_get_npix(p, i) > 0) {
2312     l_max = pict_get_lum_max(p, i);
2313     i_splitstart = igs + 1;
2314     if (l_max >= limit) {
2315     for (x = 0; x < pict_get_xsize(p); x++)
2316     for (y = 0; y < pict_get_ysize(p); y++) {
2317     if (pict_get_hangle
2318     (p, x, y, p->view.vdir, p->view.vup, &ang))
2319     {
2320     if (pict_is_validpixel(p, x, y)
2321     && luminance(pict_get_color(p, x, y))
2322     >= limit
2323     && pict_get_gsn(p, x, y) == i) {
2324     if (i_splitstart == (igs + 1)) {
2325     igs = igs + 1;
2326     pict_new_gli(p);
2327 greg 2.3 pict_get_z1_gsn(p,igs) = 0;
2328     pict_get_z2_gsn(p,igs) = 0;
2329    
2330 greg 2.1 pict_get_Eglare(p,igs) = 0.0;
2331     pict_get_Dglare(p,igs) = 0.0;
2332     pict_get_lum_min(p, igs) =
2333     99999999999999.999;
2334     i_split = igs;
2335     } else {
2336     i_split =
2337     find_split(p, x, y, r_split,
2338     i_splitstart, igs);
2339     }
2340     if (i_split == 0) {
2341     igs = igs + 1;
2342     pict_new_gli(p);
2343 greg 2.3 pict_get_z1_gsn(p,igs) = 0;
2344     pict_get_z2_gsn(p,igs) = 0;
2345    
2346 greg 2.1 pict_get_Eglare(p,igs) = 0.0;
2347     pict_get_Dglare(p,igs) = 0.0;
2348     pict_get_lum_min(p, igs) =
2349     99999999999999.999;
2350     i_split = igs;
2351     }
2352     split_pixel_from_gs(p, x, y, i_split, uniform_gs, u_r, u_g , u_b);
2353    
2354     }
2355     }
2356     }
2357    
2358     }
2359     change = 1;
2360     while (change == 1) {
2361     change = 0;
2362     for (x = 0; x < pict_get_xsize(p); x++)
2363     for (y = 0; y < pict_get_ysize(p); y++) {
2364     before_igs = pict_get_gsn(p, x, y);
2365     if (before_igs >= i_splitstart) {
2366     if (pict_get_hangle
2367     (p, x, y, p->view.vdir, p->view.vup,
2368     &ang)) {
2369     if (pict_is_validpixel(p, x, y)
2370     && before_igs > 0) {
2371     actual_igs =
2372     find_near_pgs(p, x, y,
2373     max_angle,
2374     before_igs, igs,
2375     i_splitstart);
2376     if (!(actual_igs == before_igs)) {
2377     change = 1;
2378     }
2379     if (before_igs > 0) {
2380     actual_igs =
2381     pict_get_gsn(p, x, y);
2382 greg 2.3 /* icol =
2383 greg 2.1 setglcolor(p, x, y,
2384 greg 2.3 actual_igs, uniform_gs, u_r, u_g , u_b);*/
2385 greg 2.1 }
2386    
2387     }
2388     }
2389     }
2390    
2391     }
2392     }
2393    
2394    
2395     }
2396     }
2397     }
2398    
2399 greg 2.3 /* calculation of direct vertical illuminance for CGI and for disability glare, coloring glare sources*/
2400 greg 2.1
2401     if (calcfast == 0) {
2402     for (x = 0; x < pict_get_xsize(p); x++)
2403     for (y = 0; y < pict_get_ysize(p); y++) {
2404     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2405     if (pict_is_validpixel(p, x, y)) {
2406     if (pict_get_gsn(p, x, y) > 0) {
2407 greg 2.3 actual_igs = pict_get_gsn(p, x, y);
2408     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));
2409     pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
2410     E_v_dir = E_v_dir + delta_E;
2411     setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2412 greg 2.1 }
2413     }
2414     }
2415     }
2416     lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2417 greg 2.3
2418 greg 2.1 }
2419 greg 2.3 /* calc of band luminance distribution if applied */
2420 greg 2.1 if (band_calc == 1) {
2421 greg 2.3 x_max = pict_get_xsize(p) - 1;
2422     y_max = pict_get_ysize(p) - 1;
2423     y_mid = (int)(y_max/2);
2424     for (x = 0; x <= x_max; x++)
2425     for (y = 0; y <= y_max; y++) {
2426     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2427     if (pict_is_validpixel(p, x, y)) {
2428    
2429     r_actual = acos(DOT(pict_get_cached_dir(p, x, y_mid), pict_get_cached_dir(p, x, y))) ;
2430     act_lum = luminance(pict_get_color(p, x, y));
2431    
2432     if ((r_actual <= band_angle) || (y == y_mid) ) {
2433     lum_band[0]=luminance(pict_get_color(p, x, y));
2434     muc_rvar_add_sample(s_band, lum_band);
2435     act_lum = luminance(pict_get_color(p, x, y));
2436     lum_band_av += pict_get_omega(p, x, y) * act_lum;
2437     omega_band += pict_get_omega(p, x, y);
2438     if (band_color == 1) {
2439     pict_get_color(p, x, y)[RED] = 0.0;
2440     pict_get_color(p, x, y)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
2441     pict_get_color(p, x, y)[BLU] = 0.0;
2442     }
2443     }
2444     }}}
2445     lum_band_av=lum_band_av/omega_band;
2446     muc_rvar_get_vx(s_band,lum_band_std);
2447     muc_rvar_get_percentile(s_band,lum_band_median,0.75);
2448     per_75_band=lum_band_median[0];
2449     muc_rvar_get_percentile(s_band,lum_band_median,0.95);
2450     per_95_band=lum_band_median[0];
2451     muc_rvar_get_median(s_band,lum_band_median);
2452     muc_rvar_get_bounding_box(s_band,bbox_band);
2453    
2454     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] );
2455     /* av_lum = av_lum / omega_sum; */
2456     /* printf("average luminance of band %f \n",av_lum);*/
2457     /* printf("solid angle of band %f \n",omega_sum);*/
2458     }
2459 greg 2.1
2460     /*printf("total number of glare sources: %i\n",igs); */
2461     lum_sources = 0;
2462     omega_sources = 0;
2463 greg 2.3 i=0;
2464 greg 2.1 for (x = 0; x <= igs; x++) {
2465 greg 2.3 if (pict_get_npix(p, x) > 0) {
2466 greg 2.1 lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2467     omega_sources += pict_get_av_omega(p, x);
2468 greg 2.3 i=i+1;
2469     }}
2470    
2471     if (sang == omega_sources) {
2472     lum_backg =avlum;
2473     } else {
2474     lum_backg =(sang * avlum - lum_sources) / (sang - omega_sources);
2475 greg 2.1 }
2476 greg 2.3
2477    
2478     if (i == 0) {
2479     lum_sources =0.0;
2480     } else { lum_sources=lum_sources/omega_sources;
2481     }
2482    
2483    
2484    
2485     if (non_cos_lb == 0) {
2486     lum_backg = lum_backg_cos;
2487 greg 2.1 }
2488 greg 2.3
2489     /* file writing NOT here
2490     if (checkfile == 1) {
2491     pict_write(p, file_out);
2492     }
2493     */
2494    
2495 greg 2.1 /* print detailed output */
2496 greg 2.3
2497    
2498 greg 2.1 if (detail_out == 1) {
2499 greg 2.3
2500     /* masking */
2501    
2502     if ( masking == 1 ) {
2503    
2504    
2505     for (x = 0; x < pict_get_xsize(p); x++)
2506     for (y = 0; y < pict_get_ysize(p); y++) {
2507     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2508     if (pict_is_validpixel(p, x, y)) {
2509     if (luminance(pict_get_color(pm, x, y))>0.1) {
2510     /* printf ("hallo %i %i %f \n",x,y,luminance(pict_get_color(pm, x, y)));*/
2511     i_mask=i_mask+1;
2512     lum_mask[0]=luminance(pict_get_color(p, x, y));
2513     /* printf ("%f \n",lum_mask[0]);*/
2514     muc_rvar_add_sample(s_mask, lum_mask);
2515     omega_mask += pict_get_omega(p, x, y);
2516     lum_mask_av += pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2517     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));
2518    
2519     pict_get_color(pm, x, y)[RED] = luminance(pict_get_color(p, x, y))/179.0;
2520     pict_get_color(pm, x, y)[GRN] = luminance(pict_get_color(p, x, y))/179.0;
2521     pict_get_color(pm, x, y)[BLU] = luminance(pict_get_color(p, x, y))/179.0;
2522    
2523     } else {
2524     pict_get_color(p, x, y)[RED] = 0.0 ;
2525     pict_get_color(p, x, y)[GRN] = 0.0 ;
2526     pict_get_color(p, x, y)[BLU] = 0.0 ;
2527     }
2528     }
2529     }
2530     }
2531     /* Average luminance in masking area (weighted by solid angle) */
2532     lum_mask_av=lum_mask_av/omega_mask;
2533     muc_rvar_get_vx(s_mask,lum_mask_std);
2534     muc_rvar_get_percentile(s_mask,lum_mask_median,0.75);
2535     per_75_mask=lum_mask_median[0];
2536     muc_rvar_get_percentile(s_mask,lum_mask_median,0.95);
2537     per_95_mask=lum_mask_median[0];
2538     muc_rvar_get_median(s_mask,lum_mask_median);
2539     muc_rvar_get_bounding_box(s_mask,bbox);
2540     /* PSGV only why masking of window is applied! */
2541     pgsv = get_pgsv(E_v, E_v_mask, omega_mask, lum_mask_av);
2542     pgsv_sat =get_pgsv_sat(E_v);
2543     printf ("masking:no_pixels,omega,av_lum,median_lum,std_lum,perc_75,perc_95,lum_min,lum_max,pgsv,pgsv_sat: %i %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,pgsv_sat );
2544    
2545    
2546    
2547     }
2548    
2549     /* zones */
2550    
2551     if ( zones > 0 ) {
2552    
2553     for (x = 0; x < pict_get_xsize(p); x++)
2554     for (y = 0; y < pict_get_ysize(p); y++) {
2555     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2556     if (pict_is_validpixel(p, x, y)) {
2557    
2558    
2559     r_actual =acos(DOT(pict_get_cached_dir(p, x, y), pict_get_cached_dir(p, x_zone, y_zone))) * 2;
2560     lum_actual = luminance(pict_get_color(p, x, y));
2561     act_gsn=pict_get_gsn(p, x, y);
2562     /* 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));*/
2563    
2564     /*zone1 */
2565     if (r_actual <= angle_z1) {
2566     i_z1=i_z1+1;
2567     lum_z1[0]=lum_actual;
2568     muc_rvar_add_sample(s_z1, lum_z1);
2569     omega_z1 += pict_get_omega(p, x, y);
2570     lum_z1_av += pict_get_omega(p, x, y)* lum_actual;
2571     setglcolor(p,x,y,1,1 , 0.66, 0.01 ,0.33);
2572     /*check if separation gsn already exist */
2573    
2574     if (act_gsn > 0){
2575     if (pict_get_z1_gsn(p,act_gsn) == 0) {
2576     pict_new_gli(p);
2577     igs = igs + 1;
2578     pict_get_z1_gsn(p,act_gsn) = igs*1.000;
2579     pict_get_z1_gsn(p,igs) = -1.0;
2580     pict_get_z2_gsn(p,igs) = -1.0;
2581     splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2582     /* 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)); */
2583     }
2584     splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2585     /* printf("splitgs%i \n",splitgs); */
2586     split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
2587     }
2588     }
2589     /*zone2 */
2590    
2591     if (r_actual > angle_z1 && r_actual<= angle_z2 ) {
2592    
2593     i_z2=i_z2+1;
2594     lum_z2[0]=lum_actual;
2595     muc_rvar_add_sample(s_z2, lum_z2);
2596     omega_z2 += pict_get_omega(p, x, y);
2597     lum_z2_av += pict_get_omega(p, x, y)* lum_actual;
2598     setglcolor(p,x,y,1,1 , 0.65, 0.33 ,0.02);
2599     /* 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));
2600     */ if (act_gsn > 0){
2601     if (pict_get_z2_gsn(p,act_gsn) == 0) {
2602     pict_new_gli(p);
2603     igs = igs + 1;
2604     pict_get_z2_gsn(p,act_gsn) = igs*1.000;
2605     pict_get_z1_gsn(p,igs) = -2.0;
2606     pict_get_z2_gsn(p,igs) = -2.0;
2607     /* printf ("Glare source %i gets glare source %i in zone 2 \n",act_gsn,igs ); */
2608     }
2609     splitgs=(int)(pict_get_z2_gsn(p,act_gsn));
2610     /* printf("splitgs %i \n",splitgs);*/
2611     split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
2612     }
2613     }
2614    
2615     }
2616     } }
2617     /* Average luminance in zones (weighted by solid angle), calculate percentiles, median min and max in zones */
2618     if (zones == 2 ) {
2619     lum_z2_av=lum_z2_av/omega_z2;
2620     muc_rvar_get_vx(s_z2,lum_z2_std);
2621     muc_rvar_get_percentile(s_z2,lum_z2_median,0.75);
2622     per_75_z2=lum_z2_median[0];
2623     muc_rvar_get_percentile(s_z2,lum_z2_median,0.95);
2624     per_95_z2=lum_z2_median[0];
2625     muc_rvar_get_median(s_z2,lum_z2_median);
2626     muc_rvar_get_bounding_box(s_z2,bbox_z2);
2627     }
2628     lum_z1_av=lum_z1_av/omega_z1;
2629     muc_rvar_get_vx(s_z1,lum_z1_std);
2630     muc_rvar_get_percentile(s_z1,lum_z1_median,0.75);
2631     per_75_z1=lum_z1_median[0];
2632     muc_rvar_get_percentile(s_z1,lum_z1_median,0.95);
2633     per_95_z1=lum_z1_median[0];
2634     muc_rvar_get_median(s_z1,lum_z1_median);
2635     muc_rvar_get_bounding_box(s_z1,bbox_z1);
2636    
2637     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: %f %f %f %f %f %f %f %f\n",omega_z1,lum_z1_av,lum_z1_median[0],sqrt(lum_z1_std[0]),per_75_z1,per_95_z1,bbox_z1[0],bbox_z1[1] );
2638    
2639 greg 2.4 if (zones == 2 ) {
2640 greg 2.3
2641     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: %f %f %f %f %f %f %f %f\n",omega_z2,lum_z2_av,lum_z2_median[0],sqrt(lum_z2_std[0]),per_75_z2,per_95_z2,bbox_z2[0],bbox_z2[1] );
2642 greg 2.4 }
2643 greg 2.3
2644     }
2645    
2646 greg 2.1 i = 0;
2647     for (x = 0; x <= igs; x++) {
2648     /* resorting glare source numbers */
2649     if (pict_get_npix(p, x) > 0) {
2650     i = i + 1;
2651     }
2652     }
2653 greg 2.3 no_glaresources=i;
2654 greg 2.1
2655 greg 2.3 /* glare sources */
2656 greg 2.1 printf
2657 greg 2.3 ("%i No pixels x-pos y-pos L_s Omega_s Posindx L_b L_t E_vert Edir Max_Lum Sigma xdir ydir zdir Eglare_cie Lveil_cie teta glare_zone\n",
2658 greg 2.1 i);
2659     if (i == 0) {
2660 greg 2.3 printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", i, 0.0, 0.0,
2661     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);
2662 greg 2.1
2663     } else {
2664     i = 0;
2665     for (x = 0; x <= igs; x++) {
2666     if (pict_get_npix(p, x) > 0) {
2667     i = i + 1;
2668     pict_get_sigma(p, pict_get_av_posx(p, x),
2669     pict_get_av_posy(p, x), p->view.vdir,
2670     p->view.vup, &sigma);
2671    
2672     x2=pict_get_av_posx(p, x);
2673     y2=pict_get_av_posy(p, x);
2674     teta = 180.0 / 3.1415927 * acos(DOT(p->view.vdir, pict_get_cached_dir(p, x2, y2)));
2675     Lveil_cie = 10*pict_get_Eglare(p,x)/(teta*teta+0.0000000000000001);
2676    
2677     if (apply_disability == 1 && Lveil_cie <=disability_thresh) {
2678     Lveil_cie =0 ;
2679     }
2680     Lveil_cie_sum = Lveil_cie_sum + Lveil_cie;
2681 greg 2.3 if (pict_get_z1_gsn(p,x)<0) {
2682     act_gsn=(int)(-pict_get_z1_gsn(p,x));
2683     }else{
2684     act_gsn=0;
2685     }
2686     printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %i \n",
2687 greg 2.1 i, pict_get_npix(p, x), pict_get_av_posx(p, x),
2688     pict_get_ysize(p) - pict_get_av_posy(p, x),
2689     pict_get_av_lum(p, x), pict_get_av_omega(p, x),
2690     get_posindex(p, pict_get_av_posx(p, x),
2691     pict_get_av_posy(p, x),
2692     posindex_2), lum_backg, lum_task,
2693     E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
2694 greg 2.3 ,pict_get_cached_dir(p, x2,y2)[0],pict_get_cached_dir(p, x2,y2)[1],pict_get_cached_dir(p, x2,y2)[2],pict_get_Eglare(p,x),Lveil_cie,teta,act_gsn );
2695 greg 2.1 }
2696     }
2697     }
2698     }
2699    
2700    
2701    
2702     /* calculation of indicees */
2703    
2704     /* check vertical illuminance range */
2705     E_v2 = E_v;
2706    
2707     if (E_v2 < 100) {
2708     fprintf(stderr,
2709     "Notice: Vertical illuminance is below 100 lux !!\n");
2710     }
2711     dgp =
2712     get_dgp(p, E_v2, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
2713     /* low light correction */
2714     if (E_v < 1000) {
2715     low_light_corr=1.0*exp(0.024*E_v-4)/(1+exp(0.024*E_v-4));} else {low_light_corr=1.0 ;}
2716     dgp =low_light_corr*dgp;
2717    
2718     /* age correction */
2719    
2720     if (age_corr == 1) {
2721     age_corr_factor=1.0/(1.1-0.5*age/100.0);
2722     }
2723     dgp =age_corr_factor*dgp;
2724    
2725    
2726     if (ext_vill == 1) {
2727     if (E_vl_ext < 100) {
2728     fprintf(stderr,
2729     "Notice: Vertical illuminance is below 100 lux !!\n");
2730     }
2731     }
2732    
2733     if (calcfast == 0) {
2734 greg 2.3 lum_a= E_v2/3.1415927;
2735 greg 2.1 dgi = get_dgi(p, lum_backg, igs, posindex_2);
2736     ugr = get_ugr(p, lum_backg, igs, posindex_2);
2737 greg 2.3 ugp = get_ugp (p,lum_backg, igs, posindex_2);
2738     ugr_exp = get_ugr_exp (p,lum_backg_cos,lum_a, igs, posindex_2);
2739     dgi_mod = get_dgi_mod(p, lum_a, igs, posindex_2);
2740 greg 2.1 cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
2741 greg 2.3 dgr = get_dgr(p, avlum, igs, posindex_2);
2742     vcp = get_vcp(dgr);
2743 greg 2.1 Lveil = get_disability(p, avlum, igs);
2744 greg 2.3 if (no_glaresources == 0) {
2745     dgi = 0.0;
2746     ugr = 0.0;
2747     ugp = 0.0;
2748     ugr_exp = 0.0;
2749     dgi_mod = 0.0;
2750     cgi = 0.0;
2751     dgr = 0.0;
2752     vcp = 100.0;
2753     }
2754 greg 2.1 }
2755     /* check dgp range */
2756     if (dgp <= 0.2) {
2757     fprintf(stderr,
2758     "Notice: Low brightness scene. dgp below 0.2! dgp might underestimate glare sources\n");
2759     }
2760     if (E_v < 380) {
2761     fprintf(stderr,
2762     "Notice: Low brightness scene. Vertical illuminance less than 380 lux! dgp might underestimate glare sources\n");
2763     }
2764    
2765    
2766    
2767     if (output == 0) {
2768     if (detail_out == 1) {
2769     if (ext_vill == 1) {
2770     dgp_ext =
2771     get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2772     posindex_2);
2773     dgp = dgp_ext;
2774     if (E_vl_ext < 1000) {
2775     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 ;}
2776     dgp =low_light_corr*dgp;
2777     dgp =age_corr_factor*dgp;
2778     }
2779 greg 2.3 muc_rvar_get_median(s_noposweight,lum_nopos_median);
2780     muc_rvar_get_median(s_posweight,lum_pos_median);
2781     muc_rvar_get_median(s_posweight2,lum_pos2_median);
2782    
2783    
2784 greg 2.1
2785 greg 2.3
2786 greg 2.1 printf
2787 greg 2.3 ("dgp,av_lum,E_v,lum_backg,E_v_dir,dgi,ugr,vcp,cgi,lum_sources,omega_sources,Lveil,Lveil_cie,dgr,ugp,ugr_exp,dgi_mod,av_lum_pos,av_lum_pos2,med_lum,med_lum_pos,med_lum_pos2: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
2788 greg 2.1 dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
2789 greg 2.3 lum_sources, omega_sources, Lveil,Lveil_cie_sum,dgr,ugp,ugr_exp,dgi_mod,lum_pos_mean,lum_pos2_mean/sang,lum_nopos_median[0],lum_pos_median[0],lum_pos2_median[0]);
2790 greg 2.1 } else {
2791     if (detail_out2 == 1) {
2792    
2793     printf
2794     ("dgp,dgi,ugr,vcp,cgi,dgp_ext,Ev_calc,abs_max,Lveil: %f %f %f %f %f %f %f %f %f \n",
2795     dgp, dgi, ugr, vcp, cgi, dgp_ext, E_v, abs_max,
2796     Lveil);
2797    
2798     } else {
2799     if (ext_vill == 1) {
2800     dgp_ext = get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,posindex_2);
2801    
2802     if (E_vl_ext < 1000) {
2803     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 ;}
2804     dgp =low_light_corr*dgp_ext;
2805     dgp =age_corr_factor*dgp;
2806     }
2807     printf("dgp,dgi,ugr,vcp,cgi,Lveil: %f %f %f %f %f %f \n",
2808     dgp, dgi, ugr, vcp, cgi, Lveil);
2809    
2810     }
2811     }
2812     } else {
2813     dgp_ext =
2814     get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2815     posindex_2);
2816     dgp = dgp_ext;
2817     if (E_vl_ext < 1000) {
2818     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 ;}
2819     dgp =low_light_corr*dgp;
2820     dgp =age_corr_factor*dgp;
2821     printf("%f\n", dgp);
2822     }
2823    
2824    
2825     /*printf("lowlight=%f\n", low_light_corr); */
2826    
2827    
2828     /* printf("hallo \n");
2829    
2830    
2831     apply_disability = 1;
2832     disability_thresh = atof(argv[++i]);
2833     coloring of disability glare sources red, remove all other colors
2834    
2835    
2836    
2837     this function was removed because of bugs....
2838     has to be re-written from scratch....
2839    
2840    
2841     */
2842    
2843    
2844    
2845    
2846    
2847    
2848    
2849    
2850     /*output */
2851     if (checkfile == 1) {
2852 greg 2.3
2853    
2854 greg 2.1 pict_write(p, file_out);
2855     }
2856    
2857    
2858    
2859     return EXIT_SUCCESS;
2860     exit(0);
2861    
2862     userr:
2863     fprintf(stderr,
2864     "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",
2865     progname);
2866     exit(1);
2867     }
2868    
2869    
2870    
2871     #endif