ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.2
Committed: Tue Aug 18 15:02:53 2015 UTC (8 years, 8 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad5R0
Changes since 2.1: +3 -0 lines
Log Message: 
Added RCCS id's to new source files (forgotten during initial check-in)

File Contents

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