ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
Revision: 2.6
Committed: Thu May 18 02:25:27 2017 UTC (6 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.5: +41 -9 lines
Log Message: 
Bug fixes and other changes from Jan Wienold

File Contents

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