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