ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.1
Committed: Wed Aug 12 23:07:59 2015 UTC (8 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
Log Message: 
Added Jan Wienold's evalglare to distribution

File Contents

# User Rev Content
1 greg 2.1 /* EVALGLARE V1.17
2     * Evalglare Software License, Version 1.0
3     *
4     * Copyright (c) 1995 - 2015 Fraunhofer ISE, EPFL.
5     * All rights reserved.
6     *
7     *
8     * Redistribution and use in source and binary forms, WITHOUT
9     * modification, are permitted provided that the following all conditions
10     * are met:
11     *
12     * 1. Redistributions of source code must retain the above copyright
13     * notice, this list of conditions and the following disclaimer.
14     *
15     * 2. Redistributions in binary form must reproduce the above copyright
16     * notice, this list of conditions and the following disclaimer in
17     * the documentation and/or other materials provided with the
18     * distribution.
19     *
20     * 3. The end-user documentation included with the redistribution,
21     * if any, must include the following acknowledgments:
22     * "This product includes the evalglare software,
23     developed at Fraunhofer ISE by J. Wienold" and
24     * "This product includes Radiance software
25     * (http://radsite.lbl.gov/)
26     * developed by the Lawrence Berkeley National Laboratory
27     * (http://www.lbl.gov/)."
28     * Alternately, this acknowledgment may appear in the software itself,
29     * if and wherever such third-party acknowledgments normally appear.
30     *
31     * 4. The names "Evalglare," and "Fraunhofer ISE" must
32     * not be used to endorse or promote products derived from this
33     * software without prior written permission. For written
34     * permission, please contact [email protected]
35     *
36     * 5. Products derived from this software may not be called "evalglare",
37     * nor may "evalglare" appear in their name, without prior written
38     * permission of Fraunhofer ISE.
39     *
40     * Redistribution and use in source and binary forms, WITH
41     * modification, are permitted provided that the following conditions
42     * are met:
43     *
44     *
45     * conditions 1.-5. apply
46     *
47     * 6. In order to ensure scientific correctness, any modification in source code imply fulfilling all following comditions:
48     * a) A written permission from Fraunhofer ISE. For written permission, please contact [email protected].
49     * b) The permission can be applied via email and must contain the applied changes in source code and at least two example calculations,
50     * comparing the results of the modified version and the current version of evalglare.
51     * b) Any redistribution of a modified version of evalglare must contain following note:
52     * "This software uses a modified version of the source code of evalglare."
53     *
54     *
55     *
56     *
57     * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
58     * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
59     * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
60     * DISCLAIMED. IN NO EVENT SHALL Fraunhofer ISE OR
61     * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
62     * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
63     * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
64     * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
65     * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
66     * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
67     * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
68     * SUCH DAMAGE.
69     * ====================================================================
70     *
71     * This product includes Radiance software
72     * (http://radsite.lbl.gov/)
73     * developed by the Lawrence Berkeley National Laboratory
74     * (http://www.lbl.gov/).
75     *
76     * The Radiance Software License, Version 1.0
77     *
78     * Copyright (c) 1990 - 2009 The Regents of the University of California,
79     * through Lawrence Berkeley National Laboratory. All rights reserved.
80     *
81     * Redistribution and use in source and binary forms, with or without
82     * modification, are permitted provided that the following conditions
83     * are met:
84     *
85     * 1. Redistributions of source code must retain the above copyright
86     * notice, this list of conditions and the following disclaimer.
87     *
88     * 2. Redistributions in binary form must reproduce the above copyright
89     * notice, this list of conditions and the following disclaimer in
90     * the documentation and/or other materials provided with the
91     * distribution.
92     *
93     * 3. The end-user documentation included with the redistribution,
94     * if any, must include the following acknowledgment:
95     * "This product includes Radiance software
96     * (http://radsite.lbl.gov/)
97     * developed by the Lawrence Berkeley National Laboratory
98     * (http://www.lbl.gov/)."
99     * Alternately, this acknowledgment may appear in the software itself,
100     * if and wherever such third-party acknowledgments normally appear.
101     *
102     * 4. The names "Radiance," "Lawrence Berkeley National Laboratory"
103     * and "The Regents of the University of California" must
104     * not be used to endorse or promote products derived from this
105     * software without prior written permission. For written
106     * permission, please contact [email protected].
107     *
108     * 5. Products derived from this software may not be called "Radiance",
109     * nor may "Radiance" appear in their name, without prior written
110     * permission of Lawrence Berkeley National Laboratory.
111     *
112     * THIS SOFTWARE IS PROVIDED ``AS IS'' AND ANY EXPRESSED OR IMPLIED
113     * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
114     * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
115     * DISCLAIMED. IN NO EVENT SHALL Lawrence Berkeley National Laboratory OR
116     * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
117     * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
118     * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
119     * USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
120     * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
121     * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
122     * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
123     * SUCH DAMAGE.
124     */
125    
126     /* evalglare.c, v0.2 2005/08/21 18:00:00 wienold */
127     /* evalglare.c, v0.3 2006/06/01 16:20:00 wienold
128     changes to the v02 version:
129     -fix problem with non-square pictures
130     -set luminance values in the hemisphere behind the line of sight to 0
131     -fix error in vcp calculation */
132     /* evalglare.c, v0.4 2006/11/13 15:10:00 wienold
133     changes to the v03 version:
134     -fix problem with tabulator in picture header "cannot read view"
135     -fix missing pixels in task area
136     evalglare.c, v0.5 2006/11/29 9:30 wienold
137     changes to the v04 version:
138     - fix problem with glare sources at the edge of the picture
139     */
140     /* evalglare.c, v0.6 2006/11/29 9:30 wienold
141     changes to the v05 version:
142     - add luminance restriction parameters in order to take into account restrictions of measuring systems
143     */
144    
145     /* evalglare.c, v0.7 2007/07/11 18:00 wienold
146     changes to the v06 version:
147     - add external provision of vertical illuminance
148     */
149    
150     /* evalglare.c, v0.8 2007/12/04 1:35 wienold
151     changes to the v07 version:
152     - limit dgp to max 1.0
153     - fill up cutted pictures with last known values
154     - add second detailed output
155     */
156    
157     /* evalglare.c, v0.9 2008/07/02 wienold
158     changes to the v08 version:
159     - add version number in picheader
160     - add option for showing version numer
161     */
162     /* evalglare.c, v0.9a 2008/09/20 wienold
163     changes to the v09 version:
164     - reduce fix value threshold from 500 to 100
165     */
166    
167     /* evalglare.c, v0.9b 2008/11/12 wienold
168     changes to the v09a version:
169     - check view type vta, if wrong than exit
170     */
171     /* evalglare.c, v0.9c 2009/03/31 wienold
172     changes to the v09b version:
173     - peak extraction is default now (-y) , for deactivation use -x
174     */
175    
176     /* evalglare.c, v0.9d 2009/06/24 wienold
177     changes to the v09c version:
178     - fixed memory problem while using strcat
179     */
180     /* evalglare.c, v0.9e 2009/10/10 wienold
181     changes to the v09d version:
182     - fixed problem while reading .pic of windows version
183     */
184     /* evalglare.c, v0.9f 2009/10/21 wienold
185     changes to the v09e version:
186     - piping of input pictures possible, allow also vtv and vtc viewtyps
187     */
188    
189     /* evalglare.c, v0.9g 2009/10/21 wienold
190     changes to the v09f version:
191     - modified pictool.c
192     - added -V (calc only vertical illuminance)
193     */
194     /* evalglare.c, v0.9h 2011/10/10 wienold
195     changes to the v09g version:
196     - include disability glare for age factor of 1
197     -
198     */
199     /* evalglare.c, v0.9i 2011/10/17 wienold
200     changes to the v09h version:
201     - M option: Correct wrong (measured) luminance images by feeding Ev and a luminance value, which should be replaced
202    
203     */
204     /* evalglare.c, v1.0 2012/02/08 wienold
205     changes to the v09h version:
206     - include all view types
207     - check view option in header
208     - add view options in command line
209     - option for cutting out GUTH field of view
210    
211     */
212    
213     /* evalglare.c, v1.02 2012/03/01 wienold,reetz,grobe
214     changes to the v1.0 version:
215     - fixed buffer overflow for string variable version[40]
216     - replaced various strings to handle release by #define
217     - removed most unused variables
218     - initialized variables
219     - removed nested functions
220     - compiles now with -ansi -pedantic
221    
222     */
223    
224     /* evalglare.c, v1.03 2012/04/17 wienold
225     - include low light correction
226     */
227    
228    
229     /* evalglare.c, v1.04 2012/04/23 wienold
230     - remove bug for gen_dgp_profile output
231     */
232    
233     /* evalglare.c, v1.05 2012/05/29 wienold
234     - remove variable overflow of low-light-correction for large Ev
235     */
236     /* evalglare.c, v1.06 2012/05/29 wienold
237     - initiate low-light-correction-variable
238     */
239     /* evalglare.c, v1.07 2012/05/29 wienold
240     - remove edge pixels from evaluation, when center of pixel is behind view (stability)
241     */
242    
243     /* evalglare.c, v1.08 2012/09/09 wienold
244     - add direction vector in detailed output for each glare source
245     - include age correction
246     */
247     /* evalglare.c, v1.09 2012/09/09 wienold
248     - 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
249     */
250     /* evalglare.c, v1.10 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.11 2013/01/17 wienold
254     - fix output bug of dgp, when using -i or -I
255     */
256    
257     /* evalglare.c, v1.12 2013/10/31 wienold
258     - include CIE equation for disability glare, Stiles-Holladay
259     */
260     /* evalglare.c, v1.13 2014/04/06 wienold
261     - remove bug: initialize Lveil_cie_sum
262     used for the CIE equation for disability glare
263     */
264    
265     /* evalglare.c, v1.14 buggy changes... removed...
266     */
267     /* evalglare.c, v1.15 add option for uniform colored glare sources
268     */
269     /* 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
270     calculation of the background luminance is now based on the real solid angle and not hard coded 2*PI any more
271     fast calculation mode: total solid angle =2*PI (before PI)
272     */
273    
274     /* evalglare.c, v1.17 2015/07/15 add option for calculating band luminance -B angle_of_band
275     remove of age factor due to non proven statistical evidence
276     */
277    
278     #define EVALGLARE
279    
280     #define PROGNAME "evalglare"
281     #define VERSION "1.17 release 15.07.2015 by EPFL, J.Wienold"
282     #define RELEASENAME PROGNAME " " VERSION
283    
284     #include "rtio.h"
285     #include "platform.h"
286     #include "pictool.h"
287     #include <math.h>
288     #include <string.h>
289     char *progname;
290    
291     /* subroutine to add a pixel to a glare source */
292     void add_pixel_to_gs(pict * p, int x, int y, int gsn)
293     {
294     double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
295     new_omega, act_lum;
296    
297    
298     pict_get_npix(p, gsn) = pict_get_npix(p, gsn) + 1;
299     old_av_posx = pict_get_av_posx(p, gsn);
300     old_av_posy = pict_get_av_posy(p, gsn);
301     old_av_lum = pict_get_av_lum(p, gsn);
302     old_omega = pict_get_av_omega(p, gsn);
303    
304     act_omega = pict_get_omega(p, x, y);
305     act_lum = luminance(pict_get_color(p, x, y));
306     new_omega = old_omega + act_omega;
307     pict_get_av_posx(p, gsn) =
308     (old_av_posx * old_omega + x * act_omega) / new_omega;
309     pict_get_av_posy(p, gsn) =
310     (old_av_posy * old_omega + y * act_omega) / new_omega;
311     if (isnan((pict_get_av_posx(p, gsn))))
312     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);
313    
314     pict_get_av_lum(p, gsn) =
315     (old_av_lum * old_omega + act_lum * act_omega) / new_omega;
316     pict_get_av_omega(p, gsn) = new_omega;
317     pict_get_gsn(p, x, y) = gsn;
318     if (act_lum < pict_get_lum_min(p, gsn)) {
319     pict_get_lum_min(p, gsn) = act_lum;
320     }
321     if (act_lum > pict_get_lum_max(p, gsn)) {
322     pict_get_lum_max(p, gsn) = act_lum;
323     }
324    
325     /* 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)); */
326    
327     }
328    
329     /* subroutine for peak extraction */
330     int
331     find_split(pict * p, int x, int y, double r, int i_split_start,
332     int i_split_end)
333     {
334     int i_find_split, x_min, x_max, y_min, y_max, ix, iy, iix, iiy, dx, dy,
335     out_r;
336     double r_actual;
337    
338     i_find_split = 0;
339    
340     x_min = 0;
341     y_min = 0;
342     x_max = pict_get_ysize(p) - 1;
343    
344     y_max = pict_get_ysize(p) - 1;
345    
346     for (iiy = 1; iiy <= 2; iiy++) {
347     dy = iiy * 2 - 3;
348     if (dy == -1) {
349     iy = y;
350     } else {
351     iy = y + 1;
352     }
353     while (iy <= y_max && iy >= y_min) {
354     out_r = 0;
355     for (iix = 1; iix <= 2; iix++) {
356     dx = iix * 2 - 3;
357     if (dx == -1) {
358     ix = x;
359     } else {
360     ix = x + 1;
361     }
362     while (ix <= x_max && ix >= x_min && iy >= y_min) {
363    
364     r_actual =
365     acos(DOT(pict_get_cached_dir(p, x, y),
366     pict_get_cached_dir(p, ix, iy))) * 2;
367     if (r_actual <= r) {
368     out_r = 1;
369     if (pict_get_gsn(p, ix, iy) >= i_split_start
370     && pict_get_gsn(p, ix, iy) <= i_split_end) {
371     i_find_split = pict_get_gsn(p, ix, iy);
372     }
373     } else {
374     ix = -99;
375     }
376     ix = ix + dx;
377     }
378     }
379     if (out_r == 0) {
380     iy = -99;
381     }
382     iy = iy + dy;
383    
384     }
385     }
386    
387    
388    
389     return i_find_split;
390     }
391    
392     /*
393     static int icomp(int* a,int* b)
394     {
395     if (*a < *b)
396     return -1;
397     if (*a > *b)
398     return 1;
399     return 0;
400     }
401    
402     */
403    
404     /* subroutine to find nearby glare source pixels */
405    
406     int
407     find_near_pgs(pict * p, int x, int y, float r, int act_gsn, int max_gsn,
408     int min_gsn)
409     {
410     int dx, dy, i_near_gs, xx, yy, x_min, x_max, y_min, y_max, ix, iy, iix,
411     iiy, old_gsn, new_gsn, find_gsn, change, stop_y_search,
412     stop_x_search;
413     double r_actual;
414     int ixm[3];
415    
416     i_near_gs = 0;
417     stop_y_search = 0;
418     stop_x_search = 0;
419     x_min = 0;
420     y_min = 0;
421     if (act_gsn == 0) {
422     x_max = x;
423     } else {
424     x_max = pict_get_xsize(p) - 1;
425     }
426    
427     y_max = pict_get_ysize(p) - 1;
428    
429     old_gsn = pict_get_gsn(p, x, y);
430     new_gsn = old_gsn;
431     change = 0;
432     if (act_gsn > 0) {
433     i_near_gs = pict_get_gsn(p, x, y);
434     }
435     for (iiy = 1; iiy <= 2; iiy++) {
436     dy = iiy * 2 - 3;
437     if (dy == -1) {
438     iy = y;
439     } else {
440     iy = y + 1;
441     }
442     ixm[1] = x;
443     ixm[2] = x;
444     stop_y_search = 0;
445    
446    
447     while (iy <= y_max && iy >= y_min) {
448     for (iix = 1; iix <= 2; iix++) {
449     dx = iix * 2 - 3;
450     ix = (ixm[1] + ixm[2]) / 2;
451     stop_x_search = 0;
452     while (ix <= x_max && ix >= x_min && stop_x_search == 0
453     && stop_y_search == 0) {
454     /* 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);*/
455     r_actual =
456     acos(DOT(pict_get_cached_dir(p, x, y),
457     pict_get_cached_dir(p, ix, iy))) * 2;
458     if (r_actual <= r) {
459     if (pict_get_gsn(p, ix, iy) > 0) {
460     if (act_gsn == 0) {
461     i_near_gs = pict_get_gsn(p, ix, iy);
462     stop_x_search = 1;
463     stop_y_search = 1;
464     } else {
465     find_gsn = pict_get_gsn(p, ix, iy);
466     if (pict_get_av_omega(p, old_gsn) <
467     pict_get_av_omega(p, find_gsn)
468     && pict_get_av_omega(p,
469     find_gsn) >
470     pict_get_av_omega(p, new_gsn)
471     && find_gsn >= min_gsn
472     && find_gsn <= max_gsn) {
473     /* other primary glare source found with larger solid angle -> add together all */
474     new_gsn = find_gsn;
475     change = 1;
476     stop_x_search = 1;
477     stop_y_search = 1;
478     }
479     }
480     }
481     } else {
482     ixm[iix] = ix - dx;
483     stop_x_search = 1;
484     }
485     ix = ix + dx;
486     }
487     }
488     iy = iy + dy;
489     }
490     }
491     if (change > 0) {
492     pict_get_av_lum(p, old_gsn) = 0.;
493     pict_get_av_omega(p, old_gsn) = 0.;
494     pict_get_npix(p, old_gsn) = 0.;
495     pict_get_lum_max(p, old_gsn) = 0.;
496    
497    
498     i_near_gs = new_gsn;
499     /* printf(" changing gs no %i",old_gsn);
500     printf(" to %i\n",new_gsn);
501     */
502     for (xx = 0; xx < pict_get_xsize(p); xx++)
503     for (yy = 0; yy < pict_get_ysize(p); yy++) {
504     if (pict_is_validpixel(p, x, y)) {
505     if (pict_get_gsn(p, xx, yy) == old_gsn) {
506     add_pixel_to_gs(p, xx, yy, new_gsn);
507     }
508     }
509     }
510     }
511    
512     return i_near_gs;
513    
514    
515     }
516    
517     /* subroutine for calculation of task luminance */
518     double get_task_lum(pict * p, int x, int y, float r, int task_color)
519     {
520     int x_min, x_max, y_min, y_max, ix, iy;
521     double r_actual, av_lum, omega_sum, act_lum;
522    
523    
524     x_max = pict_get_xsize(p) - 1;
525     y_max = pict_get_ysize(p) - 1;
526     x_min = 0;
527     y_min = 0;
528    
529    
530    
531     av_lum = 0.0;
532     omega_sum = 0.0;
533    
534     for (iy = y_min; iy <= y_max; iy++) {
535     for (ix = x_min; ix <= x_max; ix++) {
536    
537     /* if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
538     continue;*/
539     r_actual =
540     acos(DOT(pict_get_cached_dir(p, x, y),
541     pict_get_cached_dir(p, ix, iy))) * 2;
542     act_lum = luminance(pict_get_color(p, ix, iy));
543    
544     if (r_actual <= r) {
545     act_lum = luminance(pict_get_color(p, ix, iy));
546     av_lum += pict_get_omega(p, ix, iy) * act_lum;
547     omega_sum += pict_get_omega(p, ix, iy);
548     if (task_color == 1) {
549     pict_get_color(p, ix, iy)[RED] = 0.0;
550     pict_get_color(p, ix, iy)[GRN] = 0.0;
551     pict_get_color(p, ix, iy)[BLU] =
552     act_lum / WHTEFFICACY / CIE_bf;
553     }
554     }
555     }
556     }
557    
558     av_lum = av_lum / omega_sum;
559     /* printf("average luminance of task %f \n",av_lum);
560     printf("solid angle of task %f \n",omega_sum);*/
561     return av_lum;
562     }
563    
564     /* subroutine for calculation of band luminance */
565     double get_band_lum(pict * p, float r, int task_color)
566     {
567     int x_min, x_max, y_min, y_max, ix, iy, y_mid;
568     double r_actual, av_lum, omega_sum, act_lum;
569    
570    
571     x_max = pict_get_xsize(p) - 1;
572     y_max = pict_get_ysize(p) - 1;
573     x_min = 0;
574     y_min = 0;
575     y_mid = rint(y_max/2);
576    
577    
578    
579     av_lum = 0.0;
580     omega_sum = 0.0;
581    
582     for (iy = y_min; iy <= y_max; iy++) {
583     for (ix = x_min; ix <= x_max; ix++) {
584    
585     /* if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
586     continue;*/
587     r_actual =
588     acos(DOT(pict_get_cached_dir(p, ix, y_mid),
589     pict_get_cached_dir(p, ix, iy))) ;
590     act_lum = luminance(pict_get_color(p, ix, iy));
591    
592     if ((r_actual <= r) || (iy == y_mid) ) {
593     act_lum = luminance(pict_get_color(p, ix, iy));
594     av_lum += pict_get_omega(p, ix, iy) * act_lum;
595     omega_sum += pict_get_omega(p, ix, iy);
596     if (task_color == 1) {
597     pict_get_color(p, ix, iy)[RED] = 0.0;
598     pict_get_color(p, ix, iy)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
599     pict_get_color(p, ix, iy)[BLU] = 0.0;
600     }
601     }
602     }
603     }
604    
605     av_lum = av_lum / omega_sum;
606     /* printf("average luminance of band %f \n",av_lum);*/
607     /* printf("solid angle of band %f \n",omega_sum);*/
608     return av_lum;
609     }
610    
611    
612    
613    
614    
615     /* subroutine for coloring the glare sources
616     red is used only for explicit call of the subroutine with acol=0 , e.g. for disability glare
617     -1: set to grey*/
618     int setglcolor(pict * p, int x, int y, int acol, int uniform_gs, double u_r, double u_g ,double u_b)
619     {
620     int icol;
621     double pcol[3][1000], act_lum, l;
622    
623     l=u_r+u_g+u_b ;
624    
625     pcol[RED][0] = 1.0 / CIE_rf;
626     pcol[GRN][0] = 0.;
627     pcol[BLU][0] = 0.;
628    
629     pcol[RED][1] = 0.;
630     pcol[GRN][1] = 1.0 / CIE_gf;
631     pcol[BLU][1] = 0.;
632    
633     pcol[RED][2] = 0.;
634     pcol[GRN][2] = 0.;
635     pcol[BLU][2] = 1 / CIE_bf;
636    
637     pcol[RED][3] = 1.0 / (1.0 - CIE_bf);
638     pcol[GRN][3] = 1.0 / (1.0 - CIE_bf);
639     pcol[BLU][3] = 0.;
640    
641     pcol[RED][4] = 1.0 / (1.0 - CIE_gf);
642     pcol[GRN][4] = 0.;
643     pcol[BLU][4] = 1.0 / (1.0 - CIE_gf);
644    
645     pcol[RED][5] = 0.;
646     pcol[GRN][5] = 1.0 / (1.0 - CIE_rf);
647     pcol[BLU][5] = 1.0 / (1.0 - CIE_rf);
648    
649     pcol[RED][6] = 0.;
650     pcol[GRN][6] = 1.0 / (1.0 - CIE_rf);
651     pcol[BLU][6] = 1.0 / (1.0 - CIE_rf);
652    
653    
654     pcol[RED][999] = 1.0 / WHTEFFICACY;
655     pcol[GRN][999] = 1.0 / WHTEFFICACY;
656     pcol[BLU][999] = 1.0 / WHTEFFICACY;
657    
658    
659     pcol[RED][998] = u_r /(l* CIE_rf) ;
660     pcol[GRN][998] = u_g /(l* CIE_gf);
661     pcol[BLU][998] = u_b /(l* CIE_bf);
662     /* printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l);*/
663     icol = acol ;
664     if ( acol == -1) {icol=999;
665     }else{if (acol>0){icol = acol % 5 +1;
666     }};
667     if ( uniform_gs == 1) {icol=998;
668     } ;
669    
670     act_lum = luminance(pict_get_color(p, x, y));
671    
672     pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY;
673     pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY;
674     pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY;
675    
676     return icol;
677     }
678    
679     /* this is the smooting subroutine */
680     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)
681     {
682     int x_min, x_max, y_min, y_max, ix, iy, icol;
683     double r_actual, omega_gs, omega_total, om;
684    
685    
686    
687     omega_gs = 0.0;
688     omega_total = 0.0;
689     x_min = x - r;
690     if (x_min < 0) {
691     x_min = 0;
692     }
693     x_max = x + r;
694     if (x_max > pict_get_xsize(p) - 1) {
695     x_max = pict_get_xsize(p) - 2;
696     }
697     y_min = y - r;
698     if (y_min < 0) {
699     y_min = 0;
700     }
701     y_max = y + r;
702     if (y_max > pict_get_ysize(p) - 1) {
703     y_max = pict_get_ysize(p) - 2;
704     }
705    
706     for (ix = x_min; ix <= x_max; ix++)
707     for (iy = y_min; iy <= y_max; iy++) {
708     r_actual = sqrt((x - ix) * (x - ix) + (y - iy) * (y - iy));
709     if (r_actual <= r) {
710     om = pict_get_omega(p, ix, iy);
711     omega_total += om;
712     if (pict_get_gsn(p, ix, iy) == act_gs
713     && pict_get_pgs(p, ix, iy) == 1) {
714     omega_gs = omega_gs + 1 * om;
715     }
716     }
717     }
718    
719    
720     if (omega_gs / omega_total > 0.2) {
721     /* add pixel to gs */
722    
723     add_pixel_to_gs(p, x, y, act_gs);
724    
725     /* color pixel of gs */
726    
727     icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b);
728    
729     }
730     }
731    
732     /* subroutine for removing a pixel from a glare source */
733     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)
734     {
735     int old_gsn, icol;
736     double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
737     new_omega, act_lum;
738    
739    
740     /* change existing gs */
741     old_gsn = pict_get_gsn(p, x, y);
742     pict_get_npix(p, old_gsn) = pict_get_npix(p, old_gsn) - 1;
743    
744     act_omega = pict_get_omega(p, x, y);
745     old_av_posx = pict_get_av_posx(p, old_gsn);
746     old_av_posy = pict_get_av_posy(p, old_gsn);
747     old_omega = pict_get_av_omega(p, old_gsn);
748    
749     new_omega = old_omega - act_omega;
750     pict_get_av_omega(p, old_gsn) = new_omega;
751     pict_get_av_posx(p, old_gsn) =
752     (old_av_posx * old_omega - x * act_omega) / new_omega;
753     pict_get_av_posy(p, old_gsn) =
754     (old_av_posy * old_omega - y * act_omega) / new_omega;
755    
756    
757     act_lum = luminance(pict_get_color(p, x, y));
758     old_av_lum = pict_get_av_lum(p, old_gsn);
759     pict_get_av_lum(p, old_gsn) =
760     (old_av_lum * old_omega - act_lum * act_omega) / new_omega;
761     /* add pixel to new gs */
762    
763     add_pixel_to_gs(p, x, y, new_gsn);
764    
765     /* color pixel of new gs */
766    
767     icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b);
768    
769    
770     }
771    
772     /* subroutine for the calculation of the position index */
773     float get_posindex(pict * p, float x, float y, int postype)
774     {
775     float posindex;
776     double teta, phi, sigma, tau, deg, d, s, r, fact;
777    
778    
779     pict_get_vangle(p, x, y, p->view.vdir, p->view.vup, &phi);
780     pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &teta);
781     pict_get_sigma(p, x, y, p->view.vdir, p->view.vup, &sigma);
782     pict_get_tau(p, x, y, p->view.vdir, p->view.vup, &tau);
783    
784     /* return (phi+teta+2*3.1415927); */
785    
786     deg = 180 / 3.1415927;
787     fact = 0.8;
788     if (phi == 0) {
789     phi = 0.00001;
790     }
791     if (sigma <= 0) {
792     sigma = -sigma;
793     }
794     if (teta == 0) {
795     teta = 0.0001;
796     }
797     tau = tau * deg;
798     sigma = sigma * deg;
799    
800     posindex =
801     exp((35.2 - 0.31889 * tau -
802     1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 +
803     0.26667 * tau -
804     0.002963 * tau *
805     tau) / 100000 *
806     sigma * sigma);
807     /* below line of sight, using Iwata model */
808     if (phi < 0) {
809     d = 1 / tan(phi);
810     s = tan(teta) / tan(phi);
811     r = sqrt(1 / d * 1 / d + s * s / d / d);
812     if (r > 0.6)
813     fact = 1.2;
814     if (r > 3) {
815     fact = 1.2;
816     r = 3;
817     }
818    
819     posindex = 1 + fact * r;
820     }
821     if (posindex > 16)
822     posindex = 16;
823    
824     return posindex;
825    
826     }
827    
828    
829    
830     double get_upper_cut_2eyes(float teta)
831     {
832     double phi;
833    
834     phi=pow(7.7458218+0.00057407915*teta-0.00021746318*teta*teta+8.5572726e-6*teta*teta*teta,2);
835    
836     return phi;
837     }
838     double get_lower_cut_2eyes(float teta)
839     {
840     double phi;
841    
842     phi=1/(-0.014699242-1.5541106e-5*teta+4.6898068e-6*teta*teta-5.1539687e-8*teta*teta*teta);
843    
844     return phi;
845     }
846    
847     double get_lower_cut_central(float teta)
848     {
849     double phi;
850    
851     phi=(68.227109-2.9524084*teta+0.046674262*teta*teta)/(1-0.042317294*teta+0.00075698419*teta*teta-6.5364285e-7*teta*teta*teta);
852    
853     if (teta>73) {
854     phi=60;
855     }
856    
857     return phi;
858     }
859    
860     /* subroutine for the cutting the total field of view */
861     void cut_view_1(pict*p)
862     {
863     int x,y;
864     double border,ang,teta,phi,phi2;
865     for(x=0;x<pict_get_xsize(p)-1;x++)
866     for(y=0;y<pict_get_ysize(p)-1;y++) {
867     if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
868     if (pict_is_validpixel(p, x, y)) {
869     pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
870     pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
871     pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
872     phi=phi*180/3.1415927;
873     phi2=phi2*180/3.1415927;
874     teta=teta*180/3.1415927;
875     if (teta<0) {
876     teta=-teta;
877     }
878     if(phi2>0){
879     border=get_upper_cut_2eyes(teta);
880    
881     if (phi>border) {
882     pict_get_color(p,x,y)[RED]=0;
883     pict_get_color(p,x,y)[GRN]=0;
884     pict_get_color(p,x,y)[BLU]=0;
885     }
886     }else{
887     border=get_lower_cut_2eyes(180-teta);
888     if (-phi<border && teta>135) {
889     pict_get_color(p,x,y)[RED]=0;
890     pict_get_color(p,x,y)[GRN]=0;
891     pict_get_color(p,x,y)[BLU]=0;
892     }
893     }
894    
895    
896    
897     }else{
898     pict_get_color(p,x,y)[RED]=0;
899     pict_get_color(p,x,y)[GRN]=0;
900     pict_get_color(p,x,y)[BLU]=0;
901    
902    
903     }
904    
905     /* printf("teta,phi,border=%f %f %f\n",teta,phi,border);*/
906     }
907     }
908    
909    
910     return;
911    
912     }
913     /* subroutine for the cutting the field of view seen by both eyes*/
914     void cut_view_2(pict*p)
915     {
916    
917     int x,y;
918     double border,ang,teta,phi,phi2;
919     for(x=0;x<pict_get_xsize(p)-1;x++)
920     for(y=0;y<pict_get_ysize(p)-1;y++) {
921     if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
922     if (pict_is_validpixel(p, x, y)) {
923     pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
924     pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
925     pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
926     pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
927     phi=phi*180/3.1415927;
928     phi2=phi2*180/3.1415927;
929     teta=teta*180/3.1415927;
930     /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
931     if (teta<0) {
932     teta=-teta;
933     }
934     if(phi2>0){
935     border=60;
936    
937     if (phi>border) {
938     pict_get_color(p,x,y)[RED]=0;
939     pict_get_color(p,x,y)[GRN]=0;
940     pict_get_color(p,x,y)[BLU]=0;
941     }
942     }else{
943     border=get_lower_cut_central(180-teta);
944     if (phi>border) {
945     pict_get_color(p,x,y)[RED]=0;
946     pict_get_color(p,x,y)[GRN]=0;
947     pict_get_color(p,x,y)[BLU]=0;
948     }
949     }
950    
951    
952    
953     }else{
954     pict_get_color(p,x,y)[RED]=0;
955     pict_get_color(p,x,y)[GRN]=0;
956     pict_get_color(p,x,y)[BLU]=0;
957    
958    
959     }
960     }
961     }
962    
963    
964     return;
965    
966     }
967    
968     /* subroutine for the cutting the field of view seen by both eyes
969     luminance is modified by position index - just for experiments - undocumented */
970     void cut_view_3(pict*p)
971     {
972    
973     int x,y;
974     double border,ang,teta,phi,phi2,lum,newlum;
975     for(x=0;x<pict_get_xsize(p)-1;x++)
976     for(y=0;y<pict_get_ysize(p)-1;y++) {
977     if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
978     if (DOT(pict_get_cached_dir(p,x,y),p->view.vdir) >= 0.0) {
979     pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
980     pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&teta);
981     pict_get_sigma(p,x,y,p->view.vdir,p->view.vup,&phi);
982     pict_get_tau(p,x,y,p->view.vdir,p->view.vup,&teta);
983     phi=phi*180/3.1415927;
984     phi2=phi2*180/3.1415927;
985     teta=teta*180/3.1415927;
986     lum=luminance(pict_get_color(p,x,y));
987    
988     newlum=lum/get_posindex(p,x,y,0);
989    
990     pict_get_color(p,x,y)[RED]=newlum/WHTEFFICACY;
991     pict_get_color(p,x,y)[GRN]=newlum/WHTEFFICACY;
992     pict_get_color(p,x,y)[BLU]=newlum/WHTEFFICACY;
993    
994    
995    
996    
997     /* printf("x,y,teta,phi,border=%i %i %f %f %f\n",x,y,teta,phi,border);*/
998     if (teta<0) {
999     teta=-teta;
1000     }
1001     if(phi2>0){
1002     border=60;
1003    
1004     if (phi>border) {
1005     pict_get_color(p,x,y)[RED]=0;
1006     pict_get_color(p,x,y)[GRN]=0;
1007     pict_get_color(p,x,y)[BLU]=0;
1008     }
1009     }else{
1010     border=get_lower_cut_central(180-teta);
1011     if (phi>border) {
1012     pict_get_color(p,x,y)[RED]=0;
1013     pict_get_color(p,x,y)[GRN]=0;
1014     pict_get_color(p,x,y)[BLU]=0;
1015     }
1016     }
1017    
1018    
1019    
1020     }else{
1021     pict_get_color(p,x,y)[RED]=0;
1022     pict_get_color(p,x,y)[GRN]=0;
1023     pict_get_color(p,x,y)[BLU]=0;
1024    
1025    
1026     }
1027     }
1028     }
1029    
1030    
1031     return;
1032    
1033     }
1034    
1035     /* subroutine for the calculation of the daylight glare index */
1036    
1037    
1038     float get_dgi(pict * p, float lum_backg, int igs, int posindex_2)
1039     {
1040     float dgi, sum_glare, omega_s;
1041     int i;
1042    
1043    
1044     sum_glare = 0;
1045     omega_s = 0;
1046    
1047     for (i = 0; i <= igs; i++) {
1048    
1049     if (pict_get_npix(p, i) > 0) {
1050     omega_s =
1051     pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1052     get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1053     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));
1054     /* printf("i,sum_glare %i %f\n",i,sum_glare); */
1055     }
1056     }
1057     dgi = 10 * log10(sum_glare);
1058    
1059     return dgi;
1060    
1061     }
1062    
1063     /* subroutine for the calculation of the daylight glare probability */
1064     double
1065     get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3,
1066     double a4, double a5, double c1, double c2, double c3,
1067     int posindex_2)
1068     {
1069     double dgp;
1070     double sum_glare;
1071     int i;
1072    
1073     sum_glare = 0;
1074     if (igs > 0) {
1075     for (i = 0; i <= igs; i++) {
1076    
1077     if (pict_get_npix(p, i) > 0) {
1078     sum_glare +=
1079     pow(pict_get_av_lum(p, i),
1080     a1) / pow(get_posindex(p, pict_get_av_posx(p, i),
1081     pict_get_av_posy(p, i),
1082     posindex_2),
1083     a4) * pow(pict_get_av_omega(p, i), a2);
1084     /* printf("i,sum_glare %i %f\n",i,sum_glare);*/
1085     }
1086     }
1087     dgp =
1088     c1 * pow(E_v,
1089     a5) + c3 + c2 * log10(1 + sum_glare / pow(E_v, a3));
1090     } else {
1091     dgp = c3 + c1 * pow(E_v, a5);
1092     }
1093    
1094     if (dgp > 1) {
1095     dgp = 1;
1096     }
1097     /* printf("dgp %f\n",dgp);*/
1098     return dgp;
1099    
1100     }
1101    
1102     /* subroutine for the calculation of the visual comfort probability */
1103     float get_vcp(pict * p, double lum_a, int igs, int posindex_2)
1104     {
1105     float vcp;
1106     double sum_glare, dgr;
1107     int i, i_glare;
1108    
1109    
1110     sum_glare = 0;
1111     i_glare = 0;
1112     for (i = 0; i <= igs; i++) {
1113    
1114     if (pict_get_npix(p, i) > 0) {
1115     i_glare = i_glare + 1;
1116     sum_glare +=
1117     (0.5 * pict_get_av_lum(p, i) *
1118     (20.4 * pict_get_av_omega(p, i) +
1119     1.52 * pow(pict_get_av_omega(p, i),
1120     0.2) - 0.075)) / (get_posindex(p,
1121     pict_get_av_posx
1122     (p, i),
1123     pict_get_av_posy
1124     (p, i),
1125     posindex_2) *
1126     pow(lum_a, 0.44));
1127    
1128     }
1129     }
1130     dgr = pow(sum_glare, pow(i_glare, -0.0914));
1131    
1132     vcp = 50 * erf((6.347 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1133     if (dgr > 750) {
1134     vcp = 0;
1135     }
1136     if (dgr < 20) {
1137     vcp = 100;
1138     }
1139    
1140    
1141     return vcp;
1142    
1143     }
1144    
1145     /* subroutine for the calculation of the unified glare rating */
1146     float get_ugr(pict * p, double lum_backg, int igs, int posindex_2)
1147     {
1148     float ugr;
1149     double sum_glare;
1150     int i, i_glare;
1151    
1152    
1153     sum_glare = 0;
1154     i_glare = 0;
1155     for (i = 0; i <= igs; i++) {
1156    
1157     if (pict_get_npix(p, i) > 0) {
1158     i_glare = i_glare + 1;
1159     sum_glare +=
1160     pow(pict_get_av_lum(p, i) /
1161     get_posindex(p, pict_get_av_posx(p, i),
1162     pict_get_av_posy(p, i), posindex_2),
1163     2) * pict_get_av_omega(p, i);
1164    
1165     }
1166     }
1167     ugr = 8 * log10(0.25 / lum_backg * sum_glare);
1168    
1169     return ugr;
1170    
1171     }
1172    
1173    
1174     /* subroutine for the calculation of the disability glare according to Poynter */
1175     float get_disability(pict * p, double lum_backg, int igs)
1176     {
1177     float disab;
1178     double sum_glare, sigma, deg;
1179     int i, i_glare;
1180    
1181    
1182     sum_glare = 0;
1183     i_glare = 0;
1184     deg = 180 / 3.1415927;
1185     for (i = 0; i <= igs; i++) {
1186    
1187     if (pict_get_npix(p, i) > 0) {
1188     i_glare = i_glare + 1;
1189     pict_get_sigma(p, pict_get_av_posx(p, i),
1190     pict_get_av_posy(p, i), p->view.vdir,
1191     p->view.vup, &sigma);
1192    
1193     sum_glare +=
1194     pict_get_av_lum(p,
1195     i) * cos(sigma +
1196     0.00000000001) *
1197     pict_get_av_omega(p, i)
1198     / (deg * sigma + 0.00000000001);
1199     if (isnan(sum_glare)) {
1200     printf("sigma for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1201     printf("omega for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1202     printf("avlum for %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i));
1203     printf("avlum for %f %f %f\n", pict_get_av_posx(p, i), pict_get_av_posy(p, i), sigma);
1204     }
1205    
1206     }
1207     }
1208     disab = 5 * sum_glare;
1209    
1210     return disab;
1211    
1212     }
1213    
1214    
1215    
1216    
1217    
1218    
1219     /* subroutine for the calculation of the cie glare index */
1220    
1221     float
1222     get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1223     {
1224     float cgi;
1225     double sum_glare;
1226     int i, i_glare;
1227    
1228    
1229     sum_glare = 0;
1230     i_glare = 0;
1231     for (i = 0; i <= igs; i++) {
1232    
1233     if (pict_get_npix(p, i) > 0) {
1234     i_glare = i_glare + 1;
1235     sum_glare +=
1236     pow(pict_get_av_lum(p, i) /
1237     get_posindex(p, pict_get_av_posx(p, i),
1238     pict_get_av_posy(p, i), posindex_2),
1239     2) * pict_get_av_omega(p, i);
1240    
1241     }
1242     }
1243     cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare);
1244    
1245     return cgi;
1246    
1247     }
1248    
1249    
1250     #ifdef EVALGLARE
1251    
1252    
1253     /* main program
1254     ------------------------------------------------------------------------------------------------------------------*/
1255    
1256    
1257     int main(int argc, char **argv)
1258     {
1259     #define CLINEMAX 4095 /* memory allocated for command line string */
1260     pict *p = pict_create();
1261     int skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, y1, fill, yfillmax, yfillmin,
1262     ext_vill, set_lum_max, set_lum_max2, task_color, i_splitstart,
1263     i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,
1264     igs, actual_igs, icol, xt, yt, change, before_igs, sgs, splithigh,uniform_gs,
1265     detail_out, posindex_picture, non_cos_lb, rx, ry, rmx, rmy,apply_disability,band_calc,band_color;
1266     double lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma,
1267     E_vl_ext, lum_max, new_lum_max, r_center,
1268     search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,
1269     omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang,
1270     l_max, lum_backg, lum_backg_cos, omega_sources, lum_sources,
1271     lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum;
1272     float lum_task, lum_thres, dgi, vcp, cgi, ugr, limit,
1273     abs_max, Lveil;
1274     char file_out[500], file_out2[500], version[500];
1275     char *cline;
1276     VIEW userview = STDVIEW;
1277     int gotuserview = 0;
1278    
1279     /*set required user view parameters to invalid values*/
1280     uniform_gs = 0;
1281     apply_disability = 0;
1282     disability_thresh = 0;
1283     Lveil_cie_sum=0.0;
1284     skip_second_scan=0;
1285     lum_total_max=0.0;
1286     calcfast=0;
1287     age_corr_factor = 1.0;
1288     age_corr = 0;
1289     age = 20.0;
1290     userview.horiz = 0;
1291     userview.vert = 0;
1292     userview.type = 0;
1293     dgp_ext = 0;
1294     E_vl_ext = 0.0;
1295     new_lum_max = 0.0;
1296     lum_max = 0.0;
1297     omegat = 0.0;
1298     yt = 0;
1299     xt = 0;
1300     yfillmin = 0;
1301     yfillmax = 0;
1302     cut_view = 0;
1303     cut_view_type = 0;
1304     setvalue = 2e09;
1305     omega_cos_contr = 0.0;
1306     lum_ideal = 0.0;
1307     max_angle = 0.2;
1308     lum_thres = 5.0;
1309     task_lum = 0;
1310     sgs = 0;
1311     splithigh = 1;
1312     detail_out = 0;
1313     detail_out2 = 0;
1314     posindex_picture = 0;
1315     checkfile = 0;
1316     ext_vill = 0;
1317     fill = 0;
1318     a1 = 2.0;
1319     a2 = 1.0;
1320     a3 = 1.87;
1321     a4 = 2.0;
1322     a5 = 1.0;
1323     c1 = 5.87e-05;
1324     c2 = 0.092;
1325     c3 = 0.159;
1326     non_cos_lb = 1;
1327     posindex_2 = 0;
1328     task_color = 0;
1329     limit = 50000.0;
1330     set_lum_max = 0;
1331     set_lum_max2 = 0;
1332     abs_max = 0;
1333     progname = argv[0];
1334     E_v_contr = 0.0;
1335     strcpy(version, "1.15 release 05.02.2015 by J.Wienold");
1336     low_light_corr=1.0;
1337     output = 0;
1338     calc_vill = 0;
1339     band_avlum = -99;
1340     band_calc = 0;
1341     /* command line for output picture*/
1342    
1343     cline = (char *) malloc(CLINEMAX+1);
1344     cline[0] = '\0';
1345     for (i = 0; i < argc; i++) {
1346     /* fprintf(stderr, "%d %d \n",i,strlen(argv[i]));*/
1347     if (strlen(cline)+strlen(argv[i])+strlen(RELEASENAME)+2 >=CLINEMAX) {
1348     exit (-1);
1349     }
1350     else {
1351     strcat(cline, argv[i]);
1352     strcat(cline, " ");
1353     }
1354     }
1355     strcat(cline, "\n");
1356     strcat(cline, RELEASENAME);
1357     strcat(cline, "\n");
1358    
1359    
1360     /* program options */
1361    
1362     for (i = 1; i < argc && argv[i][0] == '-'; i++) {
1363     /* expand arguments */
1364     while ((rval = expandarg(&argc, &argv, i)) > 0);
1365     if (rval < 0) {
1366     fprintf(stderr, "%s: cannot expand '%s'", argv[0], argv[i]);
1367     exit(1);
1368     }
1369     rval = getviewopt(&userview, argc - i, argv + i);
1370     if (rval >= 0) {
1371     i += rval;
1372     gotuserview++;
1373     continue;
1374     }
1375     switch (argv[i][1]) {
1376     case 'a':
1377     age = atof(argv[++i]);
1378     age_corr = 1;
1379     printf("age factor not supported any more \n");
1380     exit(1);
1381     break;
1382     case 'b':
1383     lum_thres = atof(argv[++i]);
1384     break;
1385     case 'c':
1386     checkfile = 1;
1387     strcpy(file_out, argv[++i]);
1388     break;
1389     case 'u':
1390     uniform_gs = 1;
1391     u_r = atof(argv[++i]);
1392     u_g = atof(argv[++i]);
1393     u_b = atof(argv[++i]);
1394     break;
1395     case 'r':
1396     max_angle = atof(argv[++i]);
1397     break;
1398     case 's':
1399     sgs = 1;
1400     break;
1401     case 'k':
1402     apply_disability = 1;
1403     disability_thresh = atof(argv[++i]);
1404     break;
1405    
1406     case 'p':
1407     posindex_picture = 1;
1408     break;
1409    
1410     case 'y':
1411     splithigh = 1;
1412     break;
1413     case 'x':
1414     splithigh = 0;
1415     break;
1416     case 'Y':
1417     splithigh = 1;
1418     limit = atof(argv[++i]);
1419     break;
1420    
1421     case 'i':
1422     ext_vill = 1;
1423     E_vl_ext = atof(argv[++i]);
1424     break;
1425     case 'I':
1426     ext_vill = 1;
1427     fill = 1;
1428     E_vl_ext = atof(argv[++i]);
1429     yfillmax = atoi(argv[++i]);
1430     yfillmin = atoi(argv[++i]);
1431     break;
1432     case 'd':
1433     detail_out = 1;
1434     break;
1435     case 'D':
1436     detail_out2 = 1;
1437     break;
1438     case 'm':
1439     set_lum_max = 1;
1440     lum_max = atof(argv[++i]);
1441     new_lum_max = atof(argv[++i]);
1442     strcpy(file_out2, argv[++i]);
1443     /* printf("max lum set to %f\n",new_lum_max);*/
1444     break;
1445    
1446     case 'M':
1447     set_lum_max2 = 1;
1448     lum_max = atof(argv[++i]);
1449     E_vl_ext = atof(argv[++i]);
1450     strcpy(file_out2, argv[++i]);
1451     /* printf("max lum set to %f\n",new_lum_max);*/
1452     break;
1453     case 'n':
1454     non_cos_lb = 0;
1455     break;
1456    
1457     case 't':
1458     task_lum = 1;
1459     xt = atoi(argv[++i]);
1460     yt = atoi(argv[++i]);
1461     omegat = atof(argv[++i]);
1462     task_color = 0;
1463     break;
1464     case 'T':
1465     task_lum = 1;
1466     xt = atoi(argv[++i]);
1467     yt = atoi(argv[++i]);
1468     /* omegat= DEG2RAD(atof(argv[++i]));*/
1469     omegat = atof(argv[++i]);
1470     task_color = 1;
1471     break;
1472     case 'B':
1473     band_calc = 1;
1474     band_angle = atof(argv[++i]);
1475     band_color = 1;
1476     break;
1477    
1478    
1479     case 'w':
1480     a1 = atof(argv[++i]);
1481     a2 = atof(argv[++i]);
1482     a3 = atof(argv[++i]);
1483     a4 = atof(argv[++i]);
1484     a5 = atof(argv[++i]);
1485     c1 = atof(argv[++i]);
1486     c2 = atof(argv[++i]);
1487     c3 = atof(argv[++i]);
1488     break;
1489     case 'V':
1490     calc_vill = 1;
1491     break;
1492     case 'G':
1493     cut_view = 1;
1494     cut_view_type= atof(argv[++i]);
1495     break;
1496     case 'g':
1497     cut_view = 2;
1498     cut_view_type= atof(argv[++i]);
1499     break;
1500    
1501     /*case 'v':
1502     printf("evalglare %s \n",version);
1503     exit(1); */
1504    
1505     case '1':
1506     output = 1;
1507     break;
1508    
1509     case 'v':
1510     if (argv[i][2] == '\0') {
1511     printf("%s \n",RELEASENAME);
1512     exit(1);
1513     }
1514     if (argv[i][2] != 'f')
1515     goto userr;
1516     rval = viewfile(argv[++i], &userview, NULL);
1517     if (rval < 0) {
1518     fprintf(stderr,
1519     "%s: cannot open view file \"%s\"\n",
1520     progname, argv[i]);
1521     exit(1);
1522     } else if (rval == 0) {
1523     fprintf(stderr,
1524     "%s: bad view file \"%s\"\n", progname, argv[i]);
1525     exit(1);
1526     } else {
1527     gotuserview++;
1528     }
1529     break;
1530    
1531    
1532     default:
1533     goto userr;
1534     }
1535     }
1536    
1537     /*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/
1538    
1539     if (output == 1 && ext_vill == 1) {
1540     calcfast=1;
1541     }
1542    
1543     /* read picture file */
1544     if (i == argc) {
1545     SET_FILE_BINARY(stdin);
1546     FILE *fp = fdopen(fileno(stdin), "rb");
1547     if (!(fp)) {
1548     fprintf(stderr, "fdopen on stdin failed\n");
1549     return EXIT_FAILURE;
1550     }
1551     if (!(pict_read_fp(p, fp)))
1552     return EXIT_FAILURE;;
1553     } else {
1554     if (!pict_read(p, argv[i]))
1555     return EXIT_FAILURE;
1556     }
1557     if (gotuserview) {
1558     if (p->valid_view)
1559     fprintf(stderr,
1560     "warning: overriding image view by commandline argument\n");
1561     if ((userview.horiz == 0) || (userview.vert == 0) || (userview.type == 0)) {
1562     fprintf(stderr, "error: if specified, a view must at least contain -vt -vv and -vh\n");
1563     return EXIT_FAILURE;
1564     }
1565     p->view = userview;
1566     if (!(pict_update_view(p))) {
1567     fprintf(stderr, "error: invalid view specified");
1568     return EXIT_FAILURE;
1569     }
1570     pict_update_evalglare_caches(p);
1571     } else if (!(p->valid_view)) {
1572     fprintf(stderr, "error: no valid view specified\n");
1573     return EXIT_FAILURE;
1574     }
1575    
1576     /* fprintf(stderr,"Picture read!");*/
1577    
1578     /* command line for output picture*/
1579    
1580     pict_set_comment(p, cline);
1581    
1582     if (p->view.type == VT_PAR) {
1583     fprintf(stderr, "wrong view type! must not be parallel ! \n");
1584     exit(1);
1585     }
1586    
1587    
1588     /* Check size of search radius */
1589     rmx = (pict_get_xsize(p) / 2);
1590     rmy = (pict_get_ysize(p) / 2);
1591     rx = (pict_get_xsize(p) / 2 + 10);
1592     ry = (pict_get_ysize(p) / 2 + 10);
1593     r_center =
1594     acos(DOT(pict_get_cached_dir(p, rmx, rmy),
1595     pict_get_cached_dir(p, rx, ry))) * 2 / 10;
1596     search_pix = max_angle / r_center;
1597     if (search_pix < 1.0) {
1598     fprintf(stderr,
1599     "warning: search radius less than 1 pixel! deactivating smoothing and peak extraction...\n");
1600     splithigh = 0;
1601     sgs = 0;
1602    
1603     } else {
1604     if (search_pix < 3.0) {
1605     fprintf(stderr,
1606     "warning: search radius less than 3 pixels! -> %f \n",
1607     search_pix);
1608    
1609     }
1610     }
1611    
1612     /* Check task position */
1613    
1614     if (task_lum == 1) {
1615     if (xt >= pict_get_xsize(p) || yt >= pict_get_ysize(p) || xt < 0
1616     || yt < 0) {
1617     fprintf(stderr,
1618     "error: task position outside picture!! exit...");
1619     exit(1);
1620     }
1621     }
1622    
1623    
1624     /* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1625    
1626     sang = 0.0;
1627     E_v = 0.0;
1628     E_v_dir = 0.0;
1629     avlum = 0.0;
1630     pict_new_gli(p);
1631     igs = 0;
1632    
1633     /* cut out GUTH field of view and exit without glare evaluation */
1634     if (cut_view==2) {
1635     if (cut_view_type==1) {
1636     cut_view_1(p);
1637     }else{
1638    
1639     if (cut_view_type==2) {
1640     cut_view_2(p);
1641     }else{
1642     if (cut_view_type==3) {
1643     fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
1644     cut_view_3(p);
1645     }else{
1646     fprintf(stderr,"error: no valid option for view cutting!!");
1647     exit(1);
1648     }
1649     }}
1650     pict_write(p, file_out);
1651     exit(1); }
1652    
1653    
1654    
1655    
1656    
1657    
1658     /* write positionindex into checkfile and exit, activated by option -p */
1659     if (posindex_picture == 1) {
1660     for (x = 0; x < pict_get_xsize(p); x++)
1661     for (y = 0; y < pict_get_ysize(p); y++) {
1662     if (pict_get_hangle
1663     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1664     if (pict_is_validpixel(p, x, y)) {
1665     lum =
1666     get_posindex(p, x, y,
1667     posindex_2) / WHTEFFICACY;
1668    
1669     pict_get_color(p, x, y)[RED] = lum;
1670     pict_get_color(p, x, y)[GRN] = lum;
1671     pict_get_color(p, x, y)[BLU] = lum;
1672     lum_task = luminance(pict_get_color(p, x, y));
1673     /*printf("x,y,posindex=%i %i %f %f\n",x,y,lum*WHTEFFICACY,lum_task);*/
1674     }
1675     }
1676     }
1677     pict_write(p, file_out);
1678     exit(1);
1679     }
1680    
1681    
1682     /* calculation of vertical illuminance, average luminance, in case of filling activated-> fill picture */
1683     /* fill, if necessary from 0 to yfillmin */
1684    
1685     if (fill == 1) {
1686     for (x = 0; x < pict_get_xsize(p); x++)
1687     for (y = yfillmin; y > 0; y = y - 1) {
1688     y1 = y + 1;
1689     lum = luminance(pict_get_color(p, x, y1));
1690     pict_get_color(p, x, y)[RED] = lum / 179.0;
1691     pict_get_color(p, x, y)[GRN] = lum / 179.0;
1692     pict_get_color(p, x, y)[BLU] = lum / 179.0;
1693     }
1694     }
1695    
1696     if (calcfast == 0) {
1697     for (x = 0; x < pict_get_xsize(p); x++)
1698     for (y = 0; y < pict_get_ysize(p); y++) {
1699     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
1700     if (pict_is_validpixel(p, x, y)) {
1701     lum = luminance(pict_get_color(p, x, y));
1702     if (fill == 1 && y >= yfillmax) {
1703     y1 = y - 1;
1704     lum = luminance(pict_get_color(p, x, y1));
1705     pict_get_color(p, x, y)[RED] = lum / 179.0;
1706     pict_get_color(p, x, y)[GRN] = lum / 179.0;
1707     pict_get_color(p, x, y)[BLU] = lum / 179.0;
1708     }
1709    
1710     if (lum > abs_max) {
1711     abs_max = lum;
1712     }
1713     /* set luminance restriction, if -m is set */
1714     if (set_lum_max == 1 && lum >= lum_max) {
1715     pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
1716     pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
1717     pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
1718     /* printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
1719     lum = luminance(pict_get_color(p, x, y));
1720    
1721     }
1722     if (set_lum_max2 == 1 && lum >= lum_max) {
1723    
1724     E_v_contr +=
1725     DOT(p->view.vdir,
1726     pict_get_cached_dir(p, x,
1727     y)) * pict_get_omega(p,
1728     x,
1729     y)
1730     * lum;
1731     omega_cos_contr +=
1732     DOT(p->view.vdir,
1733     pict_get_cached_dir(p, x,
1734     y)) * pict_get_omega(p,
1735     x,
1736     y)
1737     * 1;
1738     }
1739    
1740    
1741     sang += pict_get_omega(p, x, y);
1742     E_v +=
1743     DOT(p->view.vdir,
1744     pict_get_cached_dir(p, x,
1745     y)) * pict_get_omega(p, x,
1746     y) *
1747     lum;
1748     avlum +=
1749     pict_get_omega(p, x,
1750     y) * luminance(pict_get_color(p, x,
1751     y));
1752     } else {
1753     pict_get_color(p, x, y)[RED] = 0.0;
1754     pict_get_color(p, x, y)[GRN] = 0.0;
1755     pict_get_color(p, x, y)[BLU] = 0.0;
1756    
1757     }
1758     }
1759     }
1760    
1761     if (set_lum_max2 == 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) {
1762    
1763     lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
1764     if (lum_ideal > 0 && lum_ideal < setvalue) {
1765     setvalue = lum_ideal;
1766     }
1767     printf
1768     ("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f %f %f %f %f\n",
1769     lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
1770    
1771    
1772     for (x = 0; x < pict_get_xsize(p); x++)
1773     for (y = 0; y < pict_get_ysize(p); y++) {
1774     if (pict_get_hangle
1775     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1776     if (pict_is_validpixel(p, x, y)) {
1777     lum = luminance(pict_get_color(p, x, y));
1778     if (set_lum_max2 == 1 && lum >= lum_max) {
1779    
1780     pict_get_color(p, x, y)[RED] =
1781     setvalue / 179.0;
1782     pict_get_color(p, x, y)[GRN] =
1783     setvalue / 179.0;
1784     pict_get_color(p, x, y)[BLU] =
1785     setvalue / 179.0;
1786    
1787     }
1788     }
1789     }
1790     }
1791    
1792     pict_write(p, file_out2);
1793    
1794     }
1795     if (set_lum_max == 1) {
1796     pict_write(p, file_out2);
1797    
1798     }
1799    
1800     if (calc_vill == 1) {
1801     printf("%f\n", E_v);
1802     exit(1);
1803     }
1804     }else{
1805     /* in fast calculation mode: ev=ev_ext and sang=2*pi */
1806     sang=2*3.14159265359;
1807     lum_task =E_vl_ext/sang;
1808     E_v=E_vl_ext;
1809     /* printf("calc fast!! %f %f\n", sang,lum_task);*/
1810    
1811    
1812     }
1813    
1814     /* cut out GUTH field of view for glare evaluation */
1815     if (cut_view==1) {
1816     if (cut_view_type==1) {
1817     cut_view_1(p);
1818     }else{
1819    
1820     if (cut_view_type==2) {
1821     cut_view_2(p);
1822     }else{
1823     if (cut_view_type==3) {
1824     fprintf(stderr,"warning: pixel luminance is weighted by position index - do not use image for glare evaluations!!");
1825     cut_view_3(p);
1826     }else{
1827     fprintf(stderr,"error: no valid option for view cutting!!");
1828     exit(1);
1829     }
1830     }}
1831     }
1832    
1833     if (calcfast == 0) {
1834     avlum = avlum / sang;
1835     lum_task = avlum;
1836     }
1837     /* if (ext_vill==1) {
1838     E_v=E_vl_ext;
1839     avlum=E_v/3.1415927;
1840     } */
1841    
1842     if (task_lum == 1) {
1843     lum_task = get_task_lum(p, xt, yt, omegat, task_color);
1844     }
1845     lum_source = lum_thres * lum_task;
1846     /* printf("Task Luminance=%f\n",lum_task);
1847     pict_write(p,"task.pic");*/
1848    
1849     if (lum_thres > 100) {
1850     lum_source = lum_thres;
1851     }
1852    
1853     /* printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
1854    
1855     /* first glare source scan: find primary glare sources */
1856     for (x = 0; x < pict_get_xsize(p); x++)
1857     for (y = 0; y < pict_get_ysize(p); y++) {
1858     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
1859     if (pict_is_validpixel(p, x, y)) {
1860     act_lum = luminance(pict_get_color(p, x, y));
1861     if (act_lum > lum_source) {
1862     if (act_lum >lum_total_max) {
1863     lum_total_max=act_lum;
1864     }
1865    
1866     actual_igs =
1867     find_near_pgs(p, x, y, max_angle, 0, igs, 1);
1868     if (actual_igs == 0) {
1869     igs = igs + 1;
1870     pict_new_gli(p);
1871     pict_get_lum_min(p, igs) = HUGE_VAL;
1872     pict_get_Eglare(p,igs) = 0.0;
1873     pict_get_Dglare(p,igs) = 0.0;
1874     actual_igs = igs;
1875     }
1876     icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
1877     pict_get_gsn(p, x, y) = actual_igs;
1878     pict_get_pgs(p, x, y) = 1;
1879     add_pixel_to_gs(p, x, y, actual_igs);
1880    
1881     } else {
1882     pict_get_pgs(p, x, y) = 0;
1883     }
1884     }
1885     }
1886     }
1887     /* pict_write(p,"firstscan.pic"); */
1888    
1889     if (calcfast == 1 ) {
1890     skip_second_scan=1;
1891     }
1892    
1893     /* second glare source scan: combine glare sources facing each other */
1894     change = 1;
1895     while (change == 1 && skip_second_scan==0) {
1896     change = 0;
1897     for (x = 0; x < pict_get_xsize(p); x++)
1898     for (y = 0; y < pict_get_ysize(p); y++) {
1899     before_igs = pict_get_gsn(p, x, y);
1900     if (pict_get_hangle
1901     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1902     if (pict_is_validpixel(p, x, y) && before_igs > 0) {
1903     actual_igs =
1904     find_near_pgs(p, x, y, max_angle, before_igs,
1905     igs, 1);
1906     if (!(actual_igs == before_igs)) {
1907     change = 1;
1908     }
1909     if (before_igs > 0) {
1910     actual_igs = pict_get_gsn(p, x, y);
1911     icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
1912     }
1913    
1914     }
1915     }
1916     }
1917     /* pict_write(p,"secondscan.pic");*/
1918    
1919     }
1920    
1921     /*smoothing: add secondary glare sources */
1922     i_max = igs;
1923     if (sgs == 1) {
1924    
1925     /* simplified search radius: constant for entire picture, alway circle not an angle due to performance */
1926     x = (pict_get_xsize(p) / 2);
1927     y = (pict_get_ysize(p) / 2);
1928     rx = (pict_get_xsize(p) / 2 + 10);
1929     ry = (pict_get_ysize(p) / 2 + 10);
1930     r_center =
1931     acos(DOT(pict_get_cached_dir(p, x, y),
1932     pict_get_cached_dir(p, rx, ry))) * 2 / 10;
1933     search_pix = max_angle / r_center / 1.75;
1934    
1935     for (x = 0; x < pict_get_xsize(p); x++) {
1936     for (y = 0; y < pict_get_ysize(p); y++) {
1937     if (pict_get_hangle
1938     (p, x, y, p->view.vdir, p->view.vup, &ang)) {
1939     if (pict_is_validpixel(p, x, y)
1940     && pict_get_gsn(p, x, y) == 0) {
1941     for (i = 1; i <= i_max; i++) {
1942    
1943     if (pict_get_npix(p, i) > 0) {
1944     add_secondary_gs(p, x, y, search_pix, i, uniform_gs, u_r, u_g , u_b);
1945     }
1946     }
1947    
1948     }
1949     }
1950     }
1951     }
1952    
1953     }
1954    
1955     /* extract extremes from glare sources to extra glare source */
1956     if (splithigh == 1 && lum_total_max>limit) {
1957    
1958     r_split = max_angle / 2.0;
1959     for (i = 0; i <= i_max; i++) {
1960    
1961     if (pict_get_npix(p, i) > 0) {
1962     l_max = pict_get_lum_max(p, i);
1963     i_splitstart = igs + 1;
1964     if (l_max >= limit) {
1965     for (x = 0; x < pict_get_xsize(p); x++)
1966     for (y = 0; y < pict_get_ysize(p); y++) {
1967     if (pict_get_hangle
1968     (p, x, y, p->view.vdir, p->view.vup, &ang))
1969     {
1970     if (pict_is_validpixel(p, x, y)
1971     && luminance(pict_get_color(p, x, y))
1972     >= limit
1973     && pict_get_gsn(p, x, y) == i) {
1974     if (i_splitstart == (igs + 1)) {
1975     igs = igs + 1;
1976     pict_new_gli(p);
1977     pict_get_Eglare(p,igs) = 0.0;
1978     pict_get_Dglare(p,igs) = 0.0;
1979     pict_get_lum_min(p, igs) =
1980     99999999999999.999;
1981     i_split = igs;
1982     } else {
1983     i_split =
1984     find_split(p, x, y, r_split,
1985     i_splitstart, igs);
1986     }
1987     if (i_split == 0) {
1988     igs = igs + 1;
1989     pict_new_gli(p);
1990     pict_get_Eglare(p,igs) = 0.0;
1991     pict_get_Dglare(p,igs) = 0.0;
1992     pict_get_lum_min(p, igs) =
1993     99999999999999.999;
1994     i_split = igs;
1995     }
1996     split_pixel_from_gs(p, x, y, i_split, uniform_gs, u_r, u_g , u_b);
1997    
1998     }
1999     }
2000     }
2001    
2002     }
2003     change = 1;
2004     while (change == 1) {
2005     change = 0;
2006     for (x = 0; x < pict_get_xsize(p); x++)
2007     for (y = 0; y < pict_get_ysize(p); y++) {
2008     before_igs = pict_get_gsn(p, x, y);
2009     if (before_igs >= i_splitstart) {
2010     if (pict_get_hangle
2011     (p, x, y, p->view.vdir, p->view.vup,
2012     &ang)) {
2013     if (pict_is_validpixel(p, x, y)
2014     && before_igs > 0) {
2015     actual_igs =
2016     find_near_pgs(p, x, y,
2017     max_angle,
2018     before_igs, igs,
2019     i_splitstart);
2020     if (!(actual_igs == before_igs)) {
2021     change = 1;
2022     }
2023     if (before_igs > 0) {
2024     actual_igs =
2025     pict_get_gsn(p, x, y);
2026     icol =
2027     setglcolor(p, x, y,
2028     actual_igs, uniform_gs, u_r, u_g , u_b);
2029     }
2030    
2031     }
2032     }
2033     }
2034    
2035     }
2036     }
2037    
2038    
2039     }
2040     }
2041     }
2042    
2043     /* calculation of direct vertical illuminance for CGI and for disability glare*/
2044    
2045    
2046     if (calcfast == 0) {
2047     for (x = 0; x < pict_get_xsize(p); x++)
2048     for (y = 0; y < pict_get_ysize(p); y++) {
2049     if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2050     if (pict_is_validpixel(p, x, y)) {
2051     if (pict_get_gsn(p, x, y) > 0) {
2052     pict_get_Eglare(p,pict_get_gsn(p, x, y)) +=
2053     DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2054     * pict_get_omega(p, x, y)
2055     * luminance(pict_get_color(p, x, y));
2056     E_v_dir +=
2057     DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2058     * pict_get_omega(p, x, y)
2059     * luminance(pict_get_color(p, x, y));
2060     }
2061     }
2062     }
2063     }
2064     lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2065     lum_backg = lum_backg_cos;
2066     }
2067     /* calc of band luminance if applied */
2068     if (band_calc == 1) {
2069     band_avlum=get_band_lum( p, band_angle, band_color);
2070     }
2071    
2072     /*printf("total number of glare sources: %i\n",igs); */
2073     lum_sources = 0;
2074     omega_sources = 0;
2075     for (x = 0; x <= igs; x++) {
2076     lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2077     omega_sources += pict_get_av_omega(p, x);
2078     }
2079     if (non_cos_lb == 1) {
2080     lum_backg =
2081     (sang * avlum - lum_sources) / (sang - omega_sources);
2082     }
2083     /* print detailed output */
2084     if (detail_out == 1) {
2085     i = 0;
2086     for (x = 0; x <= igs; x++) {
2087     /* resorting glare source numbers */
2088     if (pict_get_npix(p, x) > 0) {
2089     i = i + 1;
2090     }
2091     }
2092    
2093    
2094    
2095     printf
2096     ("%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 \n",
2097     i);
2098     if (i == 0) {
2099     printf("%i %f %f %f %f %.10f %f %f %f %f %f %f\n", i, 0.0, 0.0,
2100     0.0, 0.0, 0.0, 0.0, lum_backg, lum_task, E_v, E_v_dir,
2101     abs_max);
2102    
2103     } else {
2104     i = 0;
2105     for (x = 0; x <= igs; x++) {
2106     if (pict_get_npix(p, x) > 0) {
2107     i = i + 1;
2108     pict_get_sigma(p, pict_get_av_posx(p, x),
2109     pict_get_av_posy(p, x), p->view.vdir,
2110     p->view.vup, &sigma);
2111    
2112     x2=pict_get_av_posx(p, x);
2113     y2=pict_get_av_posy(p, x);
2114     teta = 180.0 / 3.1415927 * acos(DOT(p->view.vdir, pict_get_cached_dir(p, x2, y2)));
2115     Lveil_cie = 10*pict_get_Eglare(p,x)/(teta*teta+0.0000000000000001);
2116    
2117     if (apply_disability == 1 && Lveil_cie <=disability_thresh) {
2118     Lveil_cie =0 ;
2119     }
2120     Lveil_cie_sum = Lveil_cie_sum + Lveil_cie;
2121     printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
2122     i, pict_get_npix(p, x), pict_get_av_posx(p, x),
2123     pict_get_ysize(p) - pict_get_av_posy(p, x),
2124     pict_get_av_lum(p, x), pict_get_av_omega(p, x),
2125     get_posindex(p, pict_get_av_posx(p, x),
2126     pict_get_av_posy(p, x),
2127     posindex_2), lum_backg, lum_task,
2128     E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
2129     ,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
2130    
2131     );
2132     }
2133     }
2134     }
2135     }
2136    
2137    
2138    
2139     /* calculation of indicees */
2140    
2141     /* check vertical illuminance range */
2142     E_v2 = E_v;
2143    
2144     if (E_v2 < 100) {
2145     fprintf(stderr,
2146     "Notice: Vertical illuminance is below 100 lux !!\n");
2147     }
2148     dgp =
2149     get_dgp(p, E_v2, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
2150     /* low light correction */
2151     if (E_v < 1000) {
2152     low_light_corr=1.0*exp(0.024*E_v-4)/(1+exp(0.024*E_v-4));} else {low_light_corr=1.0 ;}
2153     dgp =low_light_corr*dgp;
2154    
2155     /* age correction */
2156    
2157     if (age_corr == 1) {
2158     age_corr_factor=1.0/(1.1-0.5*age/100.0);
2159     }
2160     dgp =age_corr_factor*dgp;
2161    
2162    
2163     if (ext_vill == 1) {
2164     if (E_vl_ext < 100) {
2165     fprintf(stderr,
2166     "Notice: Vertical illuminance is below 100 lux !!\n");
2167     }
2168     }
2169    
2170     if (calcfast == 0) {
2171     dgi = get_dgi(p, lum_backg, igs, posindex_2);
2172     ugr = get_ugr(p, lum_backg, igs, posindex_2);
2173     cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
2174     vcp = get_vcp(p, avlum, igs, posindex_2);
2175     Lveil = get_disability(p, avlum, igs);
2176     }
2177     /* check dgp range */
2178     if (dgp <= 0.2) {
2179     fprintf(stderr,
2180     "Notice: Low brightness scene. dgp below 0.2! dgp might underestimate glare sources\n");
2181     }
2182     if (E_v < 380) {
2183     fprintf(stderr,
2184     "Notice: Low brightness scene. Vertical illuminance less than 380 lux! dgp might underestimate glare sources\n");
2185     }
2186    
2187    
2188    
2189     if (output == 0) {
2190     if (detail_out == 1) {
2191     if (ext_vill == 1) {
2192     dgp_ext =
2193     get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2194     posindex_2);
2195     dgp = dgp_ext;
2196     if (E_vl_ext < 1000) {
2197     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 ;}
2198     dgp =low_light_corr*dgp;
2199     dgp =age_corr_factor*dgp;
2200     }
2201    
2202     printf
2203     ("dgp,av_lum,E_v,lum_backg,E_v_dir,dgi,ugr,vcp,cgi,lum_sources,omega_sources,Lveil,Lveil_cie,band_avlum: %f %f %f %f %f %f %f %f %f %f %f %f %f %f \n",
2204     dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
2205     lum_sources / omega_sources, omega_sources, Lveil,Lveil_cie_sum,band_avlum);
2206     } else {
2207     if (detail_out2 == 1) {
2208    
2209     printf
2210     ("dgp,dgi,ugr,vcp,cgi,dgp_ext,Ev_calc,abs_max,Lveil: %f %f %f %f %f %f %f %f %f \n",
2211     dgp, dgi, ugr, vcp, cgi, dgp_ext, E_v, abs_max,
2212     Lveil);
2213    
2214     } else {
2215     if (ext_vill == 1) {
2216     dgp_ext = get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,posindex_2);
2217    
2218     if (E_vl_ext < 1000) {
2219     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 ;}
2220     dgp =low_light_corr*dgp_ext;
2221     dgp =age_corr_factor*dgp;
2222     }
2223     printf("dgp,dgi,ugr,vcp,cgi,Lveil: %f %f %f %f %f %f \n",
2224     dgp, dgi, ugr, vcp, cgi, Lveil);
2225    
2226     }
2227     }
2228     } else {
2229     dgp_ext =
2230     get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
2231     posindex_2);
2232     dgp = dgp_ext;
2233     if (E_vl_ext < 1000) {
2234     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 ;}
2235     dgp =low_light_corr*dgp;
2236     dgp =age_corr_factor*dgp;
2237     printf("%f\n", dgp);
2238     }
2239    
2240    
2241     /*printf("lowlight=%f\n", low_light_corr); */
2242    
2243    
2244     /* printf("hallo \n");
2245    
2246    
2247     apply_disability = 1;
2248     disability_thresh = atof(argv[++i]);
2249     coloring of disability glare sources red, remove all other colors
2250    
2251    
2252    
2253     this function was removed because of bugs....
2254     has to be re-written from scratch....
2255    
2256    
2257     */
2258    
2259    
2260    
2261    
2262    
2263    
2264    
2265    
2266     /*output */
2267     if (checkfile == 1) {
2268     pict_write(p, file_out);
2269     }
2270    
2271    
2272    
2273     return EXIT_SUCCESS;
2274     exit(0);
2275    
2276     userr:
2277     fprintf(stderr,
2278     "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",
2279     progname);
2280     exit(1);
2281     }
2282    
2283    
2284    
2285     #endif