ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/evalglare.c
(Generate patch)

Comparing ray/src/util/evalglare.c (file contents):
Revision 2.1 by greg, Wed Aug 12 23:07:59 2015 UTC vs.
Revision 2.5 by greg, Tue Aug 2 16:35:01 2016 UTC

# Line 1 | Line 1
1 < /* EVALGLARE V1.17
2 < * Evalglare Software License, Version 1.0
1 > #ifndef lint
2 > static const char RCSid[] = "$Id$";
3 > #endif
4 > /* EVALGLARE V1.31
5 > * Evalglare Software License, Version 2.0
6   *
7 < * Copyright (c) 1995 - 2015 Fraunhofer ISE, EPFL.
7 > * Copyright (c) 1995 - 2016 Fraunhofer ISE, EPFL.
8   * All rights reserved.
9   *
10   *
# Line 20 | Line 23
23   * 3. The end-user documentation included with the redistribution,
24   *           if any, must include the following acknowledgments:
25   *             "This product includes the evalglare software,
26 <                developed at Fraunhofer ISE by J. Wienold" and
26 >                developed at Fraunhofer ISE and EPFL by J. Wienold" and
27   *             "This product includes Radiance software
28   *                 (http://radsite.lbl.gov/)
29   *                 developed by the Lawrence Berkeley National Laboratory
# Line 28 | Line 31
31   *       Alternately, this acknowledgment may appear in the software itself,
32   *       if and wherever such third-party acknowledgments normally appear.
33   *
34 < * 4. The names "Evalglare," and "Fraunhofer ISE" must
34 > * 4. The names "Evalglare," "EPFL" and "Fraunhofer ISE" must
35   *       not be used to endorse or promote products derived from this
36   *       software without prior written permission. For written
37   *       permission, please contact [email protected]
38   *
39   * 5. Products derived from this software may not be called "evalglare",
40   *       nor may "evalglare" appear in their name, without prior written
41 < *       permission of Fraunhofer ISE.
41 > *       permission of EPFL.
42   *
43   * Redistribution and use in source and binary forms, WITH
44   * modification, are permitted provided that the following conditions
# Line 45 | Line 48
48   * conditions 1.-5. apply
49   *
50   * 6. In order to ensure scientific correctness, any modification in source code imply fulfilling all following comditions:
51 < *     a) A written permission from Fraunhofer ISE. For written permission, please contact [email protected].
51 > *     a) A written permission from EPFL. For written permission, please contact [email protected].
52   *     b) The permission can be applied via email and must contain the applied changes in source code and at least two example calculations,
53   *        comparing the results of the modified version and the current version of evalglare.
54   *     b) Any redistribution of a modified version of evalglare must contain following note:
# Line 275 | Line 278
278     remove of age factor due to non proven statistical evidence
279     */
280  
281 < #define EVALGLARE
281 > /* 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 + /* 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 + /* 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 + #define EVALGLARE
316   #define PROGNAME "evalglare"
317 < #define VERSION "1.17 release 15.07.2015 by EPFL, J.Wienold"
317 > #define VERSION "1.31 release 02.08.2016 by EPFL, J.Wienold"
318   #define RELEASENAME PROGNAME " " VERSION
319  
320 < #include "rtio.h"
285 < #include "platform.h"
320 >
321   #include "pictool.h"
322 + #include "rtio.h"
323   #include <math.h>
324   #include <string.h>
325 + #include "platform.h"
326 + #include "muc_randvar.h"
327 +
328   char *progname;
329  
330   /* subroutine to add a pixel to a glare source */
# Line 561 | Line 600 | double get_task_lum(pict * p, int x, int y, float r, i
600          return av_lum;
601   }
602  
564 /* subroutine for calculation of band luminance */
565 double get_band_lum(pict * p, float r, int task_color)
566 {
567        int x_min, x_max, y_min, y_max, ix, iy, y_mid;
568        double r_actual, av_lum, omega_sum, act_lum;
603  
604  
571        x_max = pict_get_xsize(p) - 1;
572        y_max = pict_get_ysize(p) - 1;
573        x_min = 0;
574        y_min = 0;
575        y_mid = rint(y_max/2);
605  
606  
578
579        av_lum = 0.0;
580        omega_sum = 0.0;
581
582        for (iy = y_min; iy <= y_max; iy++) {
583                for (ix = x_min; ix <= x_max; ix++) {
584
585 /*                      if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
586                                continue;*/
587                        r_actual =
588                                acos(DOT(pict_get_cached_dir(p, ix, y_mid),
589                                                 pict_get_cached_dir(p, ix, iy))) ;
590                        act_lum = luminance(pict_get_color(p, ix, iy));
591
592                        if ((r_actual <= r) || (iy == y_mid) ) {
593                                act_lum = luminance(pict_get_color(p, ix, iy));
594                                av_lum += pict_get_omega(p, ix, iy) * act_lum;
595                                omega_sum += pict_get_omega(p, ix, iy);
596                                if (task_color == 1) {
597                                        pict_get_color(p, ix, iy)[RED] = 0.0;
598                                        pict_get_color(p, ix, iy)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
599                                        pict_get_color(p, ix, iy)[BLU] = 0.0;
600                                }
601                        }
602                }
603        }
604
605        av_lum = av_lum / omega_sum;
606 /*    printf("average luminance of band %f \n",av_lum);*/
607 /*      printf("solid angle of band %f \n",omega_sum);*/
608        return av_lum;
609 }
610
611
612
613
614
607   /* subroutine for coloring the glare sources
608       red is used only for explicit call of the subroutine with acol=0  , e.g. for disability glare
609       -1: set to grey*/
# Line 623 | Line 615 | int setglcolor(pict * p, int x, int y, int acol, int u
615          l=u_r+u_g+u_b ;
616          
617          pcol[RED][0] = 1.0 / CIE_rf;
618 <        pcol[GRN][0] = 0.;
619 <        pcol[BLU][0] = 0.;
618 >        pcol[GRN][0] = 0.0 / CIE_gf;
619 >        pcol[BLU][0] = 0.0 / CIE_bf;
620  
621 <        pcol[RED][1] = 0.;
621 >        pcol[RED][1] = 0.0 / CIE_rf;
622          pcol[GRN][1] = 1.0 / CIE_gf;
623 <        pcol[BLU][1] = 0.;
623 >        pcol[BLU][1] = 0.0 / CIE_bf;
624  
625 <        pcol[RED][2] = 0.;
626 <        pcol[GRN][2] = 0.;
627 <        pcol[BLU][2] = 1 / CIE_bf;
625 >        pcol[RED][2] = 0.15 / CIE_rf;
626 >        pcol[GRN][2] = 0.15 / CIE_gf;
627 >        pcol[BLU][2] = 0.7 / CIE_bf;
628  
629 <        pcol[RED][3] = 1.0 / (1.0 - CIE_bf);
630 <        pcol[GRN][3] = 1.0 / (1.0 - CIE_bf);
631 <        pcol[BLU][3] = 0.;
629 >        pcol[RED][3] = .5 / CIE_rf;
630 >        pcol[GRN][3] = .5 / CIE_gf;
631 >        pcol[BLU][3] = 0.0 / CIE_bf;
632  
633 <        pcol[RED][4] = 1.0 / (1.0 - CIE_gf);
634 <        pcol[GRN][4] = 0.;
635 <        pcol[BLU][4] = 1.0 / (1.0 - CIE_gf);
633 >        pcol[RED][4] = .5 / CIE_rf;
634 >        pcol[GRN][4] = .0 / CIE_gf;
635 >        pcol[BLU][4] = .5 / CIE_bf ;
636  
637 <        pcol[RED][5] = 0.;
638 <        pcol[GRN][5] = 1.0 / (1.0 - CIE_rf);
639 <        pcol[BLU][5] = 1.0 / (1.0 - CIE_rf);
637 >        pcol[RED][5] = .0 / CIE_rf;
638 >        pcol[GRN][5] = .5 / CIE_gf;
639 >        pcol[BLU][5] = .5 / CIE_bf;
640  
641 <        pcol[RED][6] = 0.;
642 <        pcol[GRN][6] = 1.0 / (1.0 - CIE_rf);
643 <        pcol[BLU][6] = 1.0 / (1.0 - CIE_rf);
641 >        pcol[RED][6] = 0.333 / CIE_rf;
642 >        pcol[GRN][6] = 0.333 / CIE_gf;
643 >        pcol[BLU][6] = 0.333 / CIE_bf;
644  
645  
646          pcol[RED][999] = 1.0 / WHTEFFICACY;
# Line 659 | Line 651 | int setglcolor(pict * p, int x, int y, int acol, int u
651          pcol[RED][998] = u_r /(l* CIE_rf) ;
652          pcol[GRN][998] = u_g /(l* CIE_gf);
653          pcol[BLU][998] = u_b /(l* CIE_bf);
654 < /*        printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l);*/
654 > /*        printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l,WHTEFFICACY); */
655          icol = acol ;
656          if ( acol == -1) {icol=999;
657                                    }else{if (acol>0){icol = acol % 5 +1;
# Line 672 | Line 664 | int setglcolor(pict * p, int x, int y, int acol, int u
664          pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY;
665          pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY;
666          pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY;
667 <
667 > /*        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))); */
668          return icol;
669   }
670  
# Line 724 | Line 716 | void add_secondary_gs(pict * p, int x, int y, float r,
716  
717   /* color pixel of gs */
718  
719 <                icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b);
719 > /*              icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b); */
720  
721          }
722   }
# Line 764 | Line 756 | void split_pixel_from_gs(pict * p, int x, int y, int n
756  
757   /* color pixel of new gs */
758  
759 <        icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b);
759 > /*      icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b); */
760 >        
761  
769
762   }
763  
764   /* subroutine for the calculation of the position index */
# Line 797 | Line 789 | float get_posindex(pict * p, float x, float y, int pos
789          tau = tau * deg;
790          sigma = sigma * deg;
791  
792 +        if (postype == 1) {
793 + /* KIM  model */        
794 +        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));
795 +        }else{
796 + /* Guth model, equation from IES lighting handbook */
797          posindex =
798                  exp((35.2 - 0.31889 * tau -
799                           1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 +
# Line 818 | Line 815 | float get_posindex(pict * p, float x, float y, int pos
815  
816                  posindex = 1 + fact * r;
817          }
818 <        if (posindex > 16)
818 >                if (posindex > 16)
819                  posindex = 16;
820 + }
821  
822          return posindex;
823  
# Line 1032 | Line 1030 | return;
1030  
1031   }
1032  
1033 < /* subroutine for the calculation of the daylight glare index */
1033 > /* subroutine for the calculation of the daylight glare index DGI*/
1034  
1035  
1036   float get_dgi(pict * p, float lum_backg, int igs, int posindex_2)
# Line 1059 | Line 1057 | float get_dgi(pict * p, float lum_backg, int igs, int
1057          return dgi;
1058  
1059   }
1060 + /* subroutine for the calculation of the modified daylight glare index DGI_mod (Fisekis2003)
1061 +   several glare sources are taken into account and summed up, original equation has no summary
1062 + */
1063  
1064 +
1065 + float get_dgi_mod(pict * p, float lum_a, int igs, int posindex_2)
1066 + {
1067 +        float dgi_mod, sum_glare, omega_s;
1068 +        int i;
1069 +
1070 +
1071 +        sum_glare = 0;
1072 +        omega_s = 0;
1073 +
1074 +        for (i = 0; i <= igs; i++) {
1075 +
1076 +                if (pict_get_npix(p, i) > 0) {
1077 +                        omega_s =
1078 +                                pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1079 +                                get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1080 +                        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));
1081 +                    /*    printf("i,sum_glare %i %f\n",i,sum_glare); */
1082 +                }
1083 +        }
1084 +        dgi_mod = 10 * log10(sum_glare);
1085 +
1086 +        return dgi_mod;
1087 +
1088 + }
1089 +
1090   /* subroutine for the calculation of the daylight glare probability */
1091   double
1092   get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3,
# Line 1100 | Line 1127 | get_dgp(pict * p, double E_v, int igs, double a1, doub
1127   }
1128  
1129   /* subroutine for the calculation of the visual comfort probability */
1130 < float get_vcp(pict * p, double lum_a, int igs, int posindex_2)
1130 > float get_dgr(pict * p, double lum_a, int igs, int posindex_2)
1131   {
1132 <        float vcp;
1133 <        double sum_glare, dgr;
1132 >        float dgr;
1133 >        double sum_glare;
1134          int i, i_glare;
1135  
1136  
# Line 1129 | Line 1156 | float get_vcp(pict * p, double lum_a, int igs, int pos
1156          }
1157          dgr = pow(sum_glare, pow(i_glare, -0.0914));
1158  
1159 +
1160 +
1161 +        return dgr;
1162 +
1163 + }
1164 +
1165 + float get_vcp(float dgr )
1166 + {
1167 +        float vcp;
1168 +
1169          vcp = 50 * erf((6.347 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1170          if (dgr > 750) {
1171                  vcp = 0;
# Line 1136 | Line 1173 | float get_vcp(pict * p, double lum_a, int igs, int pos
1173          if (dgr < 20) {
1174                  vcp = 100;
1175          }
1139
1140
1176          return vcp;
1177  
1178   }
1179  
1180 +
1181   /* subroutine for the calculation of the unified glare rating */
1182   float get_ugr(pict * p, double lum_backg, int igs, int posindex_2)
1183   {
# Line 1169 | Line 1205 | float get_ugr(pict * p, double lum_backg, int igs, int
1205          return ugr;
1206  
1207   }
1208 + /* subroutine for the calculation of the experimental unified glare rating (Fisekis 2003)*/
1209 + float get_ugr_exp(pict * p, double lum_backg, double lum_a, int igs, int posindex_2)
1210 + {
1211 +        float ugr_exp;
1212 +        double sum_glare;
1213 +        int i, i_glare;
1214  
1215  
1216 +        sum_glare = 0;
1217 +        i_glare = 0;
1218 +        for (i = 0; i <= igs; i++) {
1219 +
1220 +                if (pict_get_npix(p, i) > 0) {
1221 +                        i_glare = i_glare + 1;
1222 +                        sum_glare +=
1223 +                                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);
1224 +
1225 +                }
1226 +        }
1227 +        ugr_exp =8 * log10(lum_a) + 8 * log10(sum_glare/lum_backg);
1228 +
1229 +        return ugr_exp;
1230 +
1231 + }
1232 + /* subroutine for the calculation of the unified glare probability (Hirning 2014)*/
1233 + float get_ugp(pict * p, double lum_backg, int igs, int posindex_2)
1234 + {
1235 +        float ugp;
1236 +        double sum_glare;
1237 +        int i, i_glare;
1238 +
1239 +
1240 +        sum_glare = 0;
1241 +        i_glare = 0;
1242 +        for (i = 0; i <= igs; i++) {
1243 +
1244 +                if (pict_get_npix(p, i) > 0) {
1245 +                        i_glare = i_glare + 1;
1246 +                        sum_glare +=
1247 +                                pow(pict_get_av_lum(p, i) /
1248 +                                        get_posindex(p, pict_get_av_posx(p, i),
1249 +                                                                 pict_get_av_posy(p, i), posindex_2),
1250 +                                        2) * pict_get_av_omega(p, i);
1251 +
1252 +                }
1253 +        }
1254 +        ugp = 0.26 * log10(0.25 / lum_backg * sum_glare);
1255 +
1256 +        return ugp;
1257 +
1258 + }
1259 +
1260 +
1261   /* subroutine for the calculation of the disability glare according to Poynter */
1262   float get_disability(pict * p, double lum_backg, int igs)
1263   {
# Line 1218 | Line 1305 | float get_disability(pict * p, double lum_backg, int i
1305  
1306   /* subroutine for the calculation of the cie glare index */
1307  
1308 < float
1222 < get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1308 > float get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1309   {
1310          float cgi;
1311          double sum_glare;
# Line 1243 | Line 1329 | get_cgi(pict * p, double E_v, double E_v_dir, int igs,
1329          cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare);
1330  
1331          return cgi;
1332 + }      
1333  
1334 + /* 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! */
1335 + float get_pgsv(double E_v, double E_mask,double omega_mask,double lum_mask_av)
1336 + {
1337 +        float pgsv;
1338 +        double Lb;
1339 +
1340 +        Lb = (E_v-E_mask)/1.414213562373;
1341 +
1342 +        pgsv =3.2*log10(lum_mask_av)-0.64*log10(omega_mask)+(0.79*log10(omega_mask)-0.61)*log10(Lb)-8.2 ;
1343 +
1344 +
1345 +        return pgsv;
1346   }
1347  
1348 + /* 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! */
1349 + float get_pgsv_sat(double E_v)
1350 + {
1351 +        float pgsv_sat;
1352  
1353 +        pgsv_sat =3.3-(0.57+3.3)/pow((1+E_v/1.414213562373/1250),1.7);
1354 +
1355 +
1356 +        return pgsv_sat;
1357 + }
1358 +
1359 +
1360 +
1361 +
1362   #ifdef  EVALGLARE
1363  
1364  
# Line 1258 | Line 1370 | int main(int argc, char **argv)
1370   {
1371          #define CLINEMAX 4095 /* memory allocated for command line string */
1372          pict *p = pict_create();
1373 <        int     skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, y1, fill, yfillmax, yfillmin,
1374 <                ext_vill, set_lum_max, set_lum_max2, task_color, i_splitstart,
1375 <                i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,
1376 <                igs, actual_igs, icol, xt, yt, change, before_igs, sgs, splithigh,uniform_gs,
1377 <                detail_out, posindex_picture, non_cos_lb, rx, ry, rmx, rmy,apply_disability,band_calc,band_color;
1378 <        double  lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma,
1379 <                E_vl_ext, lum_max, new_lum_max, r_center,
1380 <                search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,
1381 <                omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang,
1382 <                l_max, lum_backg, lum_backg_cos, omega_sources, lum_sources,
1383 <                lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum;
1384 <        float lum_task, lum_thres, dgi,  vcp, cgi, ugr, limit,
1373 >        pict *pm = pict_create();
1374 >        int     skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, x1,y1, fill, yfillmax, yfillmin,
1375 >                ext_vill, set_lum_max, set_lum_max2, img_corr,x_disk,y_disk,task_color, i_splitstart,zones,act_gsn,splitgs,
1376 >                i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,x_zone,y_zone, i_z1, i_z2,
1377 >                igs, actual_igs, lastpixelwas_gs, icol, xt, yt, change,checkpixels, before_igs, sgs, splithigh,uniform_gs,x_max, y_max,y_mid,
1378 >                detail_out, posindex_picture, non_cos_lb, rx, ry, rmx,rmy,apply_disability,band_calc,band_color,masking,i_mask,no_glaresources,force;
1379 >        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,
1380 >                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,
1381 >                search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,r_actual,lum_actual,
1382 >                omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang, angle_z1, angle_z2,per_95_band,per_75_band,pos,
1383 >                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,
1384 >                lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum,
1385 >                lum_mask[1],lum_mask_av,lum_mask_std[1],lum_mask_median[1],omega_mask,bbox[2],
1386 >                lum_band[1],lum_band_av,lum_band_std[1],lum_band_median[1],omega_band,bbox_band[2],
1387 >                lum_z1[1],lum_z1_av,lum_z1_std[1],lum_z1_median[1],omega_z1,bbox_z1[2],
1388 >                lum_z2[1],lum_z2_av,lum_z2_std[1],lum_z2_median[1],omega_z2,bbox_z2[2],
1389 >                lum_pos[1],lum_nopos_median[1],lum_pos_median[1],lum_pos2_median[1],lum_pos_mean,lum_pos2_mean;
1390 >        float lum_task, lum_thres, dgi,  vcp, cgi, ugr, limit, dgr,
1391                  abs_max, Lveil;
1392 <        char file_out[500], file_out2[500], version[500];
1392 >        char maskfile[500],file_out[500], file_out2[500], version[500];
1393          char *cline;
1394          VIEW userview = STDVIEW;
1395          int gotuserview = 0;
1396 +        struct muc_rvar* s_mask;
1397 +        s_mask = muc_rvar_create();
1398 +        muc_rvar_set_dim(s_mask, 1);
1399 +        muc_rvar_clear(s_mask);
1400 +        struct muc_rvar* s_band;
1401 +        s_band = muc_rvar_create();
1402 +        muc_rvar_set_dim(s_band, 1);
1403 +        muc_rvar_clear(s_band);
1404 +        struct muc_rvar* s_z1;
1405 +        s_z1 = muc_rvar_create();
1406 +        muc_rvar_set_dim(s_z1, 1);
1407 +        muc_rvar_clear(s_z1);
1408  
1409 +        struct muc_rvar* s_z2;
1410 +        s_z2 = muc_rvar_create();
1411 +        muc_rvar_set_dim(s_z2, 1);
1412 +        muc_rvar_clear(s_z2);
1413 +
1414 +        struct muc_rvar* s_noposweight;
1415 +        s_noposweight = muc_rvar_create();
1416 +        muc_rvar_set_dim(s_noposweight, 1);
1417 +        muc_rvar_clear(s_noposweight);
1418 +
1419 +        struct muc_rvar* s_posweight;
1420 +        s_posweight = muc_rvar_create();
1421 +        muc_rvar_set_dim(s_posweight, 1);
1422 +        muc_rvar_clear(s_posweight);
1423 +
1424 +        struct muc_rvar* s_posweight2;
1425 +        s_posweight2 = muc_rvar_create();
1426 +        muc_rvar_set_dim(s_posweight2, 1);
1427 +        muc_rvar_clear(s_posweight2);
1428 +
1429          /*set required user view parameters to invalid values*/
1430 <        uniform_gs = 0;
1430 >        delta_E=0.0;
1431 >        no_glaresources=0;
1432 >        n_corner_px=0;
1433 >        zero_corner_px=0;
1434 >        force=0;
1435 >        dist=0.0;
1436 >        u_r=0.33333;
1437 >        u_g=0.33333;
1438 >        u_b=0.33333;
1439 >        band_angle=0;
1440 >        angle_z1=0;
1441 >        angle_z2=0;
1442 >        x_zone=0;
1443 >        y_zone=0;
1444 >        per_75_z2=0;
1445 >        per_95_z2=0;
1446 >        lum_pos_mean=0;
1447 >        lum_pos2_mean=0;
1448 >        lum_band_av = 0.0;
1449 >        omega_band = 0.0;
1450 >        pgsv = 0.0 ;
1451 >        pgsv_sat = 0.0 ;
1452 >        E_v_mask = 0.0;
1453 >        lum_z1_av = 0.0;
1454 >        omega_z1 = 0.0;
1455 >        lum_z2_av = 0.0;
1456 >        omega_z2 = 0.0 ;
1457 >        i_z1 = 0;
1458 >        i_z2 = 0;        
1459 >        zones = 0;
1460 >        lum_z2_av = 0.0;
1461 >        uniform_gs = 0;
1462          apply_disability = 0;
1463          disability_thresh = 0;
1464          Lveil_cie_sum=0.0;
# Line 1297 | Line 1478 | int main(int argc, char **argv)
1478          omegat = 0.0;
1479          yt = 0;
1480          xt = 0;
1481 +        x_disk=0;
1482 +        y_disk=0;
1483 +        angle_disk=0;
1484          yfillmin = 0;
1485          yfillmax = 0;
1486          cut_view = 0;
# Line 1329 | Line 1513 | int main(int argc, char **argv)
1513          limit = 50000.0;
1514          set_lum_max = 0;
1515          set_lum_max2 = 0;
1516 +        img_corr=0;
1517          abs_max = 0;
1518          progname = argv[0];
1519          E_v_contr = 0.0;
1520 <        strcpy(version, "1.15 release 05.02.2015 by J.Wienold");
1520 >        strcpy(version, "1.19 release 09.12.2015 by J.Wienold");
1521          low_light_corr=1.0;
1522          output = 0;
1523          calc_vill = 0;
1524          band_avlum = -99;
1525          band_calc = 0;
1526 +        masking = 0;
1527 +        lum_mask[0]=0.0;
1528 +        lum_mask_av=0.0;
1529 +        omega_mask=0.0;
1530 +        i_mask=0;
1531 +        actual_igs=0;
1532   /* command line for output picture*/
1533  
1534          cline = (char *) malloc(CLINEMAX+1);
# Line 1379 | Line 1570 | int main(int argc, char **argv)
1570                          printf("age factor not supported any more \n");
1571                          exit(1);
1572                          break;
1573 +                case 'A':
1574 +                        masking = 1;
1575 +                        detail_out = 1;
1576 +                        strcpy(maskfile, argv[++i]);
1577 +                        pict_read(pm,maskfile );
1578 +
1579 +                        break;
1580                  case 'b':
1581                          lum_thres = atof(argv[++i]);
1582                          break;
# Line 1398 | Line 1596 | int main(int argc, char **argv)
1596                  case 's':
1597                          sgs = 1;
1598                          break;
1599 +                case 'f':
1600 +                        force = 1;
1601 +                        break;
1602                  case 'k':
1603                          apply_disability = 1;
1604                          disability_thresh = atof(argv[++i]);
# Line 1406 | Line 1607 | int main(int argc, char **argv)
1607                  case 'p':
1608                          posindex_picture = 1;
1609                          break;
1610 +                case 'P':
1611 +                        posindex_2 = atoi(argv[++i]);
1612 +                        break;
1613 +                        
1614  
1615                  case 'y':
1616                          splithigh = 1;
# Line 1436 | Line 1641 | int main(int argc, char **argv)
1641                          detail_out2 = 1;
1642                          break;
1643                  case 'm':
1644 +                        img_corr = 1;
1645                          set_lum_max = 1;
1646                          lum_max = atof(argv[++i]);
1647                          new_lum_max = atof(argv[++i]);
# Line 1444 | Line 1650 | int main(int argc, char **argv)
1650                          break;
1651  
1652                  case 'M':
1653 +                        img_corr = 1;
1654                          set_lum_max2 = 1;
1655                          lum_max = atof(argv[++i]);
1656                          E_vl_ext = atof(argv[++i]);
1657                          strcpy(file_out2, argv[++i]);
1658   /*                      printf("max lum set to %f\n",new_lum_max);*/
1659                          break;
1660 +                case 'N':
1661 +                        img_corr = 1;
1662 +                        set_lum_max2 = 2;
1663 +                        x_disk = atoi(argv[++i]);
1664 +                        y_disk = atoi(argv[++i]);
1665 +                        angle_disk = atof(argv[++i]);
1666 +                        E_vl_ext = atof(argv[++i]);
1667 +                        strcpy(file_out2, argv[++i]);
1668 + /*                      printf("max lum set to %f\n",new_lum_max);*/
1669 +                        break;
1670 +
1671                  case 'n':
1672                          non_cos_lb = 0;
1673                          break;
# Line 1469 | Line 1687 | int main(int argc, char **argv)
1687                          omegat = atof(argv[++i]);
1688                          task_color = 1;
1689                          break;
1690 +                case 'l':
1691 +                        zones = 1;
1692 +                        x_zone = atoi(argv[++i]);
1693 +                        y_zone = atoi(argv[++i]);
1694 +                        angle_z1 = atof(argv[++i]);
1695 +                        angle_z2 = -1;
1696 +                        break;
1697 +                case 'L':
1698 +                        zones = 2;
1699 +                        x_zone = atoi(argv[++i]);
1700 +                        y_zone = atoi(argv[++i]);
1701 +                        angle_z1 = atof(argv[++i]);
1702 +                        angle_z2 = atof(argv[++i]);
1703 +                        break;
1704 +                        
1705 +                        
1706                  case 'B':
1707                          band_calc = 1;
1708                          band_angle = atof(argv[++i]);
# Line 1539 | Line 1773 | int main(int argc, char **argv)
1773   if (output == 1 && ext_vill == 1) {
1774                         calcfast=1;
1775                         }
1776 + /*masking and zoning cannot be applied at the same time*/
1777  
1778 + if (masking ==1 && zones >0) {
1779 +               fprintf(stderr,  " masking and zoning cannot be activated at the same time!\n");
1780 +               exit(1);
1781 + }
1782 +
1783   /* read picture file */
1784          if (i == argc) {
1785                  SET_FILE_BINARY(stdin);
# Line 1573 | Line 1813 | if (output == 1 && ext_vill == 1) {
1813                  return EXIT_FAILURE;
1814          }
1815  
1816 +
1817 +
1818   /*                fprintf(stderr,"Picture read!");*/
1819  
1820   /* command line for output picture*/
1821  
1822          pict_set_comment(p, cline);
1823  
1824 + /* several checks */
1825          if (p->view.type == VT_PAR) {
1826 <                fprintf(stderr, "wrong view type! must not be parallel ! \n");
1826 >                fprintf(stderr, "error: wrong view type! must not be parallel ! \n");
1827                  exit(1);
1828          }
1829  
1830 +        if (p->view.type == VT_PER) {
1831 +                fprintf(stderr, "warning: image has perspective view type specified in header ! \n");
1832 +        }
1833  
1834 +        if (masking == 1) {
1835 +
1836 +                if (!pict_get_xsize(p)==pict_get_xsize(pm) || !pict_get_ysize(p)==pict_get_ysize(pm)) {
1837 +                fprintf(stderr, "error: masking image has other resolution than main image ! \n");
1838 +                fprintf(stderr, "size must be identical \n");
1839 +                printf("resolution main image : %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
1840 +                printf("resolution masking image : %dx%d\n",pict_get_xsize(pm),pict_get_ysize(pm));
1841 +                exit(1);
1842 +                
1843 +                }
1844 +
1845 +        }
1846 + /*      printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1847 +
1848   /* Check size of search radius */
1849          rmx = (pict_get_xsize(p) / 2);
1850          rmy = (pict_get_ysize(p) / 2);
# Line 1629 | Line 1889 | if (output == 1 && ext_vill == 1) {
1889          avlum = 0.0;
1890          pict_new_gli(p);
1891          igs = 0;
1892 +        pict_get_z1_gsn(p,igs) = 0;
1893 +        pict_get_z2_gsn(p,igs) = 0;
1894  
1895 +
1896   /* cut out GUTH field of view and exit without glare evaluation */
1897   if (cut_view==2) {
1898          if (cut_view_type==1) {
# Line 1696 | Line 1959 | if (cut_view==2) {
1959          if (calcfast == 0) {
1960          for (x = 0; x < pict_get_xsize(p); x++)
1961                  for (y = 0; y < pict_get_ysize(p); y++) {
1962 +                   lum = luminance(pict_get_color(p, x, y));              
1963 +                   dist=sqrt((x-rmx)*(x-rmx)+(y-rmy)*(y-rmy));
1964 +                  if (dist>((rmx+rmy)/2)) {
1965 +                     n_corner_px=n_corner_px+1;
1966 +                     if (lum<7.0) {
1967 +                         zero_corner_px=zero_corner_px+1;
1968 +                         }
1969 +                     }
1970 +                
1971 +                
1972                          if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
1973                                  if (pict_is_validpixel(p, x, y)) {
1701                                        lum = luminance(pict_get_color(p, x, y));
1974                                          if (fill == 1 && y >= yfillmax) {
1975                                                  y1 = y - 1;
1976                                                  lum = luminance(pict_get_color(p, x, y1));
# Line 1711 | Line 1983 | if (cut_view==2) {
1983                                                  abs_max = lum;
1984                                          }
1985   /* set luminance restriction, if -m is set */
1986 <                                        if (set_lum_max == 1 && lum >= lum_max) {
1987 <                                                pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
1988 <                                                pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
1989 <                                                pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
1990 < /*                                    printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
1991 <                                                lum = luminance(pict_get_color(p, x, y));
1986 >                                        if (img_corr == 1 ) {
1987 >                                                if (set_lum_max == 1 && lum >= lum_max) {
1988 >                                                        pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
1989 >                                                        pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
1990 >                                                        pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
1991 > /*                                                      printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
1992 >                                                        lum = luminance(pict_get_color(p, x, y));
1993 >                                                        }
1994 >                                                if (set_lum_max2 == 1 && lum >= lum_max) {
1995  
1996 <                                        }
1997 <                                        if (set_lum_max2 == 1 && lum >= lum_max) {
1996 >                                                        E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
1997 >                                                        omega_cos_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
1998 >                                                        }
1999 >                                                if (set_lum_max2 == 2 ) {
2000 >                                                        r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2001 >                                                        if (x_disk == x && y_disk==y ) r_actual=0.0;
2002  
2003 <                                                E_v_contr +=
2004 <                                                        DOT(p->view.vdir,
2005 <                                                                pict_get_cached_dir(p, x,
2006 <                                                                                                        y)) * pict_get_omega(p,
2007 <                                                                                                                                                 x,
2008 <                                                                                                                                                 y)
2009 <                                                        * lum;
2010 <                                                omega_cos_contr +=
2011 <                                                        DOT(p->view.vdir,
2012 <                                                                pict_get_cached_dir(p, x,
1734 <                                                                                                        y)) * pict_get_omega(p,
1735 <                                                                                                                                                 x,
1736 <                                                                                                                                                 y)
1737 <                                                        * 1;
2003 >                                                        act_lum = luminance(pict_get_color(p, x, y));
2004 >
2005 >                                                       if (r_actual <= angle_disk) {
2006 >                                                              E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2007 >                                                              omega_cos_contr += DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2008 >                                                                }
2009 >
2010 >                                                
2011 >                                                
2012 >                                                        }
2013                                          }
2014 +                                        om = pict_get_omega(p, x, y);
2015 +                                        sang += om;
2016 +                                        E_v += DOT(p->view.vdir, pict_get_cached_dir(p, x,y)) * om *lum;
2017 +                                        avlum += om * lum;
2018 +                                        pos=get_posindex(p, x, y,posindex_2);
2019  
2020 +                                        lum_pos[0]=lum;
2021 +                                        muc_rvar_add_sample(s_noposweight,lum_pos);
2022 +                                        lum_pos[0]=lum/pos;
2023 +                                        lum_pos_mean+=lum_pos[0]*om;
2024 +                                        muc_rvar_add_sample(s_posweight, lum_pos);
2025 +                                        lum_pos[0]=lum_pos[0]/pos;
2026 +                                        lum_pos2_mean+=lum_pos[0]*om;
2027 +                                        muc_rvar_add_sample(s_posweight2,lum_pos);
2028  
2029 <                                        sang += pict_get_omega(p, x, y);
2030 <                                        E_v +=
1743 <                                                DOT(p->view.vdir,
1744 <                                                        pict_get_cached_dir(p, x,
1745 <                                                                                                y)) * pict_get_omega(p, x,
1746 <                                                                                                                                         y) *
1747 <                                                lum;
1748 <                                        avlum +=
1749 <                                                pict_get_omega(p, x,
1750 <                                                                           y) * luminance(pict_get_color(p, x,
1751 <                                                                                                                                         y));
2029 >
2030 >
2031                                  } else {
2032                                          pict_get_color(p, x, y)[RED] = 0.0;
2033                                          pict_get_color(p, x, y)[GRN] = 0.0;
2034                                          pict_get_color(p, x, y)[BLU] = 0.0;
2035  
2036                                  }
2037 <                        }
2037 >                        }else {
2038 >                                        pict_get_color(p, x, y)[RED] = 0.0;
2039 >                                        pict_get_color(p, x, y)[GRN] = 0.0;
2040 >                                        pict_get_color(p, x, y)[BLU] = 0.0;
2041 >
2042 >                                }
2043                  }
2044  
2045 <        if (set_lum_max2 == 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) {
2045 > /* check if image has black corners AND a perspective view */
2046  
2047 +       if (zero_corner_px/n_corner_px > 0.70 && (p->view.type == VT_PER) ) {
2048 +       printf (" corner pixels are to  %f %% black! \n",zero_corner_px/n_corner_px*100);
2049 +                printf("error! The image has black corners AND header contains a perspective view type definition !!!\n") ;
2050 +                
2051 +                if (force==0) {
2052 +                                printf("stopping...!!!!\n") ;
2053 +
2054 +                      exit(1);
2055 +                 }
2056 +       }
2057 +
2058 +
2059 +
2060 +
2061 +        lum_pos_mean= lum_pos_mean/sang;
2062 +        lum_pos2_mean= lum_pos2_mean/sang;
2063 +
2064 +        if (set_lum_max2 >= 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) {
2065 +
2066                  lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
2067 +                if (set_lum_max2 == 2 && lum_ideal >= 2e9) {
2068 +                printf("warning! luminance of replacement pixels would be larger than 2e9 cd/m2. Value set to 2e9cd/m2!\n") ;
2069 +                }
2070 +
2071                  if (lum_ideal > 0 && lum_ideal < setvalue) {
2072                          setvalue = lum_ideal;
2073                  }
2074 <                printf
1768 <                        ("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f  %f %f %f %f\n",
2074 >                printf("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f  %f %f %f %f\n",
2075                           lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
2076  
2077 <
2077 >            
2078                  for (x = 0; x < pict_get_xsize(p); x++)
2079                          for (y = 0; y < pict_get_ysize(p); y++) {
2080                                  if (pict_get_hangle
2081                                          (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2082                                          if (pict_is_validpixel(p, x, y)) {
2083                                                  lum = luminance(pict_get_color(p, x, y));
2084 +                                                
2085 +                                                
2086                                                  if (set_lum_max2 == 1 && lum >= lum_max) {
2087  
2088                                                          pict_get_color(p, x, y)[RED] =
# Line 1784 | Line 2092 | if (cut_view==2) {
2092                                                          pict_get_color(p, x, y)[BLU] =
2093                                                                  setvalue / 179.0;
2094  
2095 +                                                }else{ if(set_lum_max2 ==2 ) {
2096 +                                                        r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2097 +                                                        if (x_disk == x && y_disk==y ) r_actual=0.0;
2098 +
2099 +                                                       if (r_actual <= angle_disk) {
2100 +                                                      
2101 +                                                        pict_get_color(p, x, y)[RED] =
2102 +                                                                setvalue / 179.0;
2103 +                                                        pict_get_color(p, x, y)[GRN] =
2104 +                                                                setvalue / 179.0;
2105 +                                                        pict_get_color(p, x, y)[BLU] =
2106 +                                                                setvalue / 179.0;
2107 +                                                      
2108 +                                                       }
2109 +                                                                
2110 +                                                
2111                                                  }
2112 +                                                }
2113                                          }
2114                                  }
2115                          }
2116 +                        
2117  
2118                  pict_write(p, file_out2);
2119 <
2119 >        exit(1);
2120          }
2121          if (set_lum_max == 1) {
2122                  pict_write(p, file_out2);
# Line 1850 | Line 2176 | if (cut_view==1) {
2176                  lum_source = lum_thres;
2177          }
2178  
2179 <        /*     printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2179 > /*     printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2180  
2181   /* first glare source scan: find primary glare sources */
2182 <        for (x = 0; x < pict_get_xsize(p); x++)
2182 >        for (x = 0; x < pict_get_xsize(p); x++) {
2183 > /*                lastpixelwas_gs=0; */
2184                  for (y = 0; y < pict_get_ysize(p); y++) {
2185                          if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2186                                  if (pict_is_validpixel(p, x, y)) {
# Line 1862 | Line 2189 | if (cut_view==1) {
2189                                                  if (act_lum >lum_total_max) {
2190                                                                               lum_total_max=act_lum;
2191                                                                                 }
2192 <                                                
2193 <                                                actual_igs =
2194 <                                                        find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2192 > /* speed improvement first scan: when last pixel was glare source, then use this glare source number instead of searching.
2193 >   has been deactivated with v1.25
2194 >                      
2195 >                                                if (lastpixelwas_gs==0 || search_pix <= 1.0 ) { */
2196 >                                                actual_igs = find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2197 > /*  }*/
2198                                                  if (actual_igs == 0) {
2199                                                          igs = igs + 1;
2200                                                          pict_new_gli(p);
2201                                                          pict_get_lum_min(p, igs) = HUGE_VAL;
2202                                                          pict_get_Eglare(p,igs) = 0.0;
2203                                                          pict_get_Dglare(p,igs) = 0.0;
2204 +                                                        pict_get_z1_gsn(p,igs) = 0;
2205 +                                                        pict_get_z2_gsn(p,igs) = 0;
2206                                                          actual_igs = igs;
2207 +                                                        
2208                                                  }
2209 <                                                icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2209 > /* no coloring any more here            icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2210                                                  pict_get_gsn(p, x, y) = actual_igs;
2211                                                  pict_get_pgs(p, x, y) = 1;
2212                                                  add_pixel_to_gs(p, x, y, actual_igs);
2213 + /*                                                lastpixelwas_gs=actual_igs; */
2214  
2215                                          } else {
2216                                                  pict_get_pgs(p, x, y) = 0;
2217 + /*                                               lastpixelwas_gs=0;*/
2218                                          }
2219                                  }
2220                          }
2221 <                }
2221 >                }
2222 >             }
2223   /*      pict_write(p,"firstscan.pic");   */
2224  
2225 < if (calcfast == 1 ) {
2225 >
2226 > if (calcfast == 1 || search_pix <= 1.0) {
2227     skip_second_scan=1;
2228     }
2229  
2230   /* second glare source scan: combine glare sources facing each other */
2231          change = 1;
2232 +        i = 0;
2233          while (change == 1 && skip_second_scan==0) {
2234                  change = 0;
2235 <                for (x = 0; x < pict_get_xsize(p); x++)
2235 >                i = i+1;
2236 >                for (x = 0; x < pict_get_xsize(p); x++) {
2237                          for (y = 0; y < pict_get_ysize(p); y++) {
1899                                before_igs = pict_get_gsn(p, x, y);
2238                                  if (pict_get_hangle
2239                                          (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2240 <                                        if (pict_is_validpixel(p, x, y) && before_igs > 0) {
2240 >                                        checkpixels=1;
2241 >                                        before_igs = pict_get_gsn(p, x, y);
2242 >
2243 > /* speed improvement: search for other nearby glare sources only, when at least one adjacend pixel has another glare source number  */
2244 >                                        if (x >1 && x<pict_get_xsize(p)-2 && y >1 && y<pict_get_ysize(p)-2 )   {
2245 >                                                x1=x-1;
2246 >                                                x2=x+1;
2247 >                                                y1=y-1;
2248 >                                                y2=y+1;
2249 >                                            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) ) {
2250 >                                            checkpixels = 0;
2251 >                                             actual_igs = before_igs;
2252 >                                             }  }
2253 >
2254 >
2255 >                                        if (pict_is_validpixel(p, x, y) && before_igs > 0 && checkpixels==1) {
2256                                                  actual_igs =
2257                                                          find_near_pgs(p, x, y, max_angle, before_igs,
2258                                                                                    igs, 1);
# Line 1908 | Line 2261 | if (calcfast == 1 ) {
2261                                                  }
2262                                                  if (before_igs > 0) {
2263                                                          actual_igs = pict_get_gsn(p, x, y);
2264 <                                                        icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2264 >                                                /*      icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2265                                                  }
2266  
2267                                          }
2268                                  }
2269                          }
2270   /*      pict_write(p,"secondscan.pic");*/
2271 +        }}
2272  
1919        }
1920
2273   /*smoothing: add secondary glare sources */
2274          i_max = igs;
2275          if (sgs == 1) {
# Line 1974 | Line 2326 | if (calcfast == 1 ) {
2326                                                                          if (i_splitstart == (igs + 1)) {
2327                                                                                  igs = igs + 1;
2328                                                                                  pict_new_gli(p);
2329 +                                                                                pict_get_z1_gsn(p,igs) = 0;
2330 +                                                                                pict_get_z2_gsn(p,igs) = 0;
2331 +
2332                                                                                  pict_get_Eglare(p,igs) = 0.0;
2333                                                                                  pict_get_Dglare(p,igs) = 0.0;
2334                                                                                  pict_get_lum_min(p, igs) =
# Line 1987 | Line 2342 | if (calcfast == 1 ) {
2342                                                                          if (i_split == 0) {
2343                                                                                  igs = igs + 1;
2344                                                                                  pict_new_gli(p);
2345 +                                                                                pict_get_z1_gsn(p,igs) = 0;
2346 +                                                                                pict_get_z2_gsn(p,igs) = 0;
2347 +
2348                                                                                  pict_get_Eglare(p,igs) = 0.0;
2349                                                                                  pict_get_Dglare(p,igs) = 0.0;
2350                                                                                  pict_get_lum_min(p, igs) =
# Line 2023 | Line 2381 | if (calcfast == 1 ) {
2381                                                                                  if (before_igs > 0) {
2382                                                                                          actual_igs =
2383                                                                                                  pict_get_gsn(p, x, y);
2384 <                                                                                        icol =
2384 > /*                                                                                     icol =
2385                                                                                                  setglcolor(p, x, y,
2386 <                                                                                                                   actual_igs, uniform_gs, u_r, u_g , u_b);
2386 >                                                                                                                   actual_igs, uniform_gs, u_r, u_g , u_b);*/
2387                                                                                  }
2388  
2389                                                                          }
# Line 2040 | Line 2398 | if (calcfast == 1 ) {
2398                  }
2399          }
2400  
2401 < /* calculation of direct vertical illuminance for CGI and for disability glare*/
2401 > /* calculation of direct vertical illuminance for CGI and for disability glare, coloring glare sources*/
2402  
2045
2403          if (calcfast == 0) {
2404          for (x = 0; x < pict_get_xsize(p); x++)
2405                  for (y = 0; y < pict_get_ysize(p); y++) {
2406                          if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2407                                  if (pict_is_validpixel(p, x, y)) {
2408                                          if (pict_get_gsn(p, x, y) > 0) {
2409 <                                                pict_get_Eglare(p,pict_get_gsn(p, x, y)) +=
2410 <                                                        DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2411 <                                                        * pict_get_omega(p, x, y)
2412 <                                                        * luminance(pict_get_color(p, x, y));
2413 <                                                E_v_dir +=
2057 <                                                        DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2058 <                                                        * pict_get_omega(p, x, y)
2059 <                                                        * luminance(pict_get_color(p, x, y));
2409 >                                                actual_igs = pict_get_gsn(p, x, y);
2410 >                                                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));
2411 >                                                pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
2412 >                                                E_v_dir = E_v_dir + delta_E;
2413 >                                                setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2414                                          }
2415                                  }
2416                          }
2417                  }
2418          lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2419 <        lum_backg = lum_backg_cos;
2419 >
2420           }
2421 < /* calc of band luminance if applied */
2421 > /* calc of band luminance distribution if applied */
2422          if (band_calc == 1) {
2423 <        band_avlum=get_band_lum( p, band_angle, band_color);
2424 <        }
2423 >                x_max = pict_get_xsize(p) - 1;
2424 >                y_max = pict_get_ysize(p) - 1;
2425 >                y_mid = (int)(y_max/2);
2426 >                for (x = 0; x <= x_max; x++)
2427 >                for (y = 0; y <= y_max; y++) {
2428 >                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2429 >                                if (pict_is_validpixel(p, x, y)) {
2430  
2431 +                        r_actual = acos(DOT(pict_get_cached_dir(p, x, y_mid), pict_get_cached_dir(p, x, y))) ;
2432 +                        act_lum = luminance(pict_get_color(p, x, y));
2433 +
2434 +                        if ((r_actual <= band_angle) || (y == y_mid) ) {
2435 +                                lum_band[0]=luminance(pict_get_color(p, x, y));
2436 +                                muc_rvar_add_sample(s_band, lum_band);
2437 +                                act_lum = luminance(pict_get_color(p, x, y));
2438 +                                lum_band_av += pict_get_omega(p, x, y) * act_lum;
2439 +                                omega_band += pict_get_omega(p, x, y);
2440 +                                if (band_color == 1) {
2441 +                                        pict_get_color(p, x, y)[RED] = 0.0;
2442 +                                        pict_get_color(p, x, y)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
2443 +                                        pict_get_color(p, x, y)[BLU] = 0.0;
2444 +                                }
2445 +                        }
2446 +                }}}
2447 +                lum_band_av=lum_band_av/omega_band;
2448 +                muc_rvar_get_vx(s_band,lum_band_std);
2449 +                muc_rvar_get_percentile(s_band,lum_band_median,0.75);
2450 +                per_75_band=lum_band_median[0];
2451 +                muc_rvar_get_percentile(s_band,lum_band_median,0.95);
2452 +                per_95_band=lum_band_median[0];
2453 +                muc_rvar_get_median(s_band,lum_band_median);
2454 +                muc_rvar_get_bounding_box(s_band,bbox_band);
2455 +        
2456 +                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] );
2457 + /*      av_lum = av_lum / omega_sum; */
2458 + /*    printf("average luminance of band %f \n",av_lum);*/
2459 + /*      printf("solid angle of band %f \n",omega_sum);*/
2460 +        }
2461 +
2462   /*printf("total number of glare sources: %i\n",igs); */
2463          lum_sources = 0;
2464          omega_sources = 0;
2465 +        i=0;
2466          for (x = 0; x <= igs; x++) {
2467 +        if (pict_get_npix(p, x) > 0) {
2468                  lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2469                  omega_sources += pict_get_av_omega(p, x);
2470 +                i=i+1;
2471 +        }}
2472 +      
2473 +        if (sang == omega_sources) {
2474 +               lum_backg =avlum;
2475 +        } else {
2476 +               lum_backg =(sang * avlum - lum_sources) / (sang - omega_sources);
2477          }
2478 <        if (non_cos_lb == 1) {
2479 <                lum_backg =
2480 <                        (sang * avlum - lum_sources) / (sang - omega_sources);
2478 >
2479 >
2480 >        if (i == 0) {
2481 >               lum_sources =0.0;
2482 >        } else { lum_sources=lum_sources/omega_sources;
2483          }
2484 +
2485 +
2486 +
2487 +        if (non_cos_lb == 0) {
2488 +                        lum_backg = lum_backg_cos;
2489 +        }
2490 +
2491 + /* file writing NOT here
2492 +        if (checkfile == 1) {
2493 +                pict_write(p, file_out);
2494 +        }
2495 + */
2496 +
2497   /* print detailed output */
2498 +        
2499 +        
2500 +
2501 + /* masking */
2502 +
2503 +        if ( masking == 1 ) {
2504 +        
2505 +        
2506 +                for (x = 0; x < pict_get_xsize(p); x++)
2507 +                for (y = 0; y < pict_get_ysize(p); y++) {
2508 +                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2509 +                                if (pict_is_validpixel(p, x, y)) {
2510 +                                        if (luminance(pict_get_color(pm, x, y))>0.1) {
2511 + /*                                         printf ("hallo %i %i %f \n",x,y,luminance(pict_get_color(pm, x, y)));*/
2512 +                                                  i_mask=i_mask+1;
2513 +                                                  lum_mask[0]=luminance(pict_get_color(p, x, y));
2514 +                                                 /* printf ("%f \n",lum_mask[0]);*/
2515 +                                                  muc_rvar_add_sample(s_mask, lum_mask);
2516 +                                                  omega_mask += pict_get_omega(p, x, y);
2517 +                                                  lum_mask_av += pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2518 +                                                  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));
2519 +
2520 +                                                  pict_get_color(pm, x, y)[RED] = luminance(pict_get_color(p, x, y))/179.0;
2521 +                                                  pict_get_color(pm, x, y)[GRN] = luminance(pict_get_color(p, x, y))/179.0;
2522 +                                                  pict_get_color(pm, x, y)[BLU] = luminance(pict_get_color(p, x, y))/179.0;
2523 +  
2524 +                                           } else {
2525 +                                           pict_get_color(p, x, y)[RED] = 0.0 ;
2526 +                                           pict_get_color(p, x, y)[GRN] = 0.0 ;
2527 +                                           pict_get_color(p, x, y)[BLU] = 0.0 ;
2528 +                                           }
2529 +                                }
2530 +                        }
2531 +                }
2532 + /* Average luminance in masking area (weighted by solid angle) */
2533 +                lum_mask_av=lum_mask_av/omega_mask;
2534 +                muc_rvar_get_vx(s_mask,lum_mask_std);
2535 +                muc_rvar_get_percentile(s_mask,lum_mask_median,0.75);
2536 +                per_75_mask=lum_mask_median[0];
2537 +                muc_rvar_get_percentile(s_mask,lum_mask_median,0.95);
2538 +                per_95_mask=lum_mask_median[0];
2539 +                muc_rvar_get_median(s_mask,lum_mask_median);
2540 +                muc_rvar_get_bounding_box(s_mask,bbox);
2541 + /* PSGV only why masking of window is applied! */
2542 +                 pgsv = get_pgsv(E_v, E_v_mask, omega_mask, lum_mask_av);
2543 +                 pgsv_sat =get_pgsv_sat(E_v);
2544 +
2545          if (detail_out == 1) {
2546 +
2547 +                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 );
2548 +
2549 +        }      
2550 +                
2551 +        }
2552 +
2553 + /* zones */
2554 +
2555 +        if ( zones > 0 ) {
2556 +        
2557 +                for (x = 0; x < pict_get_xsize(p); x++)
2558 +                for (y = 0; y < pict_get_ysize(p); y++) {
2559 +                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2560 +                                if (pict_is_validpixel(p, x, y)) {
2561 +
2562 +
2563 +                                 r_actual =acos(DOT(pict_get_cached_dir(p, x, y), pict_get_cached_dir(p, x_zone, y_zone))) * 2;
2564 +                                 lum_actual = luminance(pict_get_color(p, x, y));
2565 +                                 act_gsn=pict_get_gsn(p, x, y);
2566 +                            /*     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));*/
2567 +                                
2568 + /*zone1 */
2569 +                        if (r_actual <= angle_z1) {
2570 +                                                  i_z1=i_z1+1;
2571 +                                                  lum_z1[0]=lum_actual;
2572 +                                                  muc_rvar_add_sample(s_z1, lum_z1);
2573 +                                                  omega_z1 += pict_get_omega(p, x, y);
2574 +                                                  lum_z1_av += pict_get_omega(p, x, y)* lum_actual;
2575 +                                                  setglcolor(p,x,y,1,1 , 0.66, 0.01 ,0.33);
2576 + /*check if separation gsn already exist */
2577 +
2578 +                                                   if (act_gsn > 0){
2579 +                                                   if (pict_get_z1_gsn(p,act_gsn) == 0) {
2580 +                                                                                          pict_new_gli(p);
2581 +                                                                                          igs = igs + 1;
2582 +                                                                                          pict_get_z1_gsn(p,act_gsn) = igs*1.000;
2583 +                                                                                          pict_get_z1_gsn(p,igs) = -1.0;
2584 +                                                                                          pict_get_z2_gsn(p,igs) = -1.0;
2585 +                                                                                          splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2586 + /*                                                                                        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)); */
2587 +                                                                                         }
2588 +                                                    splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2589 +                                        /*          printf("splitgs%i \n",splitgs);       */              
2590 +                                                    split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
2591 +                                                }                                
2592 +                                                  }
2593 + /*zone2 */
2594 +
2595 +                        if (r_actual > angle_z1 && r_actual<= angle_z2 ) {
2596 +
2597 +                                                  i_z2=i_z2+1;
2598 +                                                  lum_z2[0]=lum_actual;
2599 +                                                  muc_rvar_add_sample(s_z2, lum_z2);
2600 +                                                  omega_z2 +=   pict_get_omega(p, x, y);
2601 +                                                  lum_z2_av += pict_get_omega(p, x, y)* lum_actual;
2602 +                                                  setglcolor(p,x,y,1,1 , 0.65, 0.33 ,0.02);
2603 + /*                                                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));
2604 + */                                                if (act_gsn > 0){
2605 +                                                   if (pict_get_z2_gsn(p,act_gsn) == 0) {
2606 +                                                                                          pict_new_gli(p);
2607 +                                                                                          igs = igs + 1;
2608 +                                                                                          pict_get_z2_gsn(p,act_gsn) = igs*1.000;
2609 +                                                                                          pict_get_z1_gsn(p,igs) = -2.0;
2610 +                                                                                          pict_get_z2_gsn(p,igs) = -2.0;
2611 + /*                                                                                         printf ("Glare source %i gets glare source %i in zone 2 \n",act_gsn,igs ); */
2612 +                                                                                         }
2613 +                                                splitgs=(int)(pict_get_z2_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 +
2619 +                        }
2620 +                } }
2621 + /* Average luminance in zones (weighted by solid angle), calculate percentiles, median min and max in zones  */
2622 +               if (zones == 2 ) {
2623 +                lum_z2_av=lum_z2_av/omega_z2;
2624 +                muc_rvar_get_vx(s_z2,lum_z2_std);
2625 +                muc_rvar_get_percentile(s_z2,lum_z2_median,0.75);
2626 +                per_75_z2=lum_z2_median[0];
2627 +                muc_rvar_get_percentile(s_z2,lum_z2_median,0.95);
2628 +                per_95_z2=lum_z2_median[0];
2629 +                muc_rvar_get_median(s_z2,lum_z2_median);
2630 +                muc_rvar_get_bounding_box(s_z2,bbox_z2);
2631 +                }
2632 +                lum_z1_av=lum_z1_av/omega_z1;
2633 +                muc_rvar_get_vx(s_z1,lum_z1_std);
2634 +                muc_rvar_get_percentile(s_z1,lum_z1_median,0.75);
2635 +                per_75_z1=lum_z1_median[0];
2636 +                muc_rvar_get_percentile(s_z1,lum_z1_median,0.95);
2637 +                per_95_z1=lum_z1_median[0];
2638 +                muc_rvar_get_median(s_z1,lum_z1_median);
2639 +                muc_rvar_get_bounding_box(s_z1,bbox_z1);
2640 +        if (detail_out == 1) {
2641 +
2642 +                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] );
2643 +
2644 +               if (zones == 2 ) {
2645 +
2646 +                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] );
2647 + } }            
2648 +                
2649 +        }
2650 +
2651                  i = 0;
2652                  for (x = 0; x <= igs; x++) {
2653   /* resorting glare source numbers */
# Line 2089 | Line 2655 | if (calcfast == 1 ) {
2655                                  i = i + 1;
2656                          }
2657                  }
2658 +                no_glaresources=i;
2659  
2660 + /* glare sources */
2661 +        if (detail_out == 1) {
2662  
2094
2663                  printf
2664 <                        ("%i No pixels x-pos y-pos L_s Omega_s Posindx L_b L_t E_vert Edir Max_Lum Sigma xdir ydir zdir Eglare_cie Lveil_cie \n",
2664 >                        ("%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",
2665                           i);
2666                  if (i == 0) {
2667 <                        printf("%i %f %f %f %f %.10f %f %f %f %f %f %f\n", i, 0.0, 0.0,
2668 <                                   0.0, 0.0, 0.0, 0.0, lum_backg, lum_task, E_v, E_v_dir,
2101 <                                   abs_max);
2667 >                        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,
2668 >                                   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);
2669  
2670                  } else {
2671                          i = 0;
# Line 2118 | Line 2685 | if (calcfast == 1 ) {
2685                                                                       Lveil_cie =0 ;
2686                                                                     }
2687                                          Lveil_cie_sum =  Lveil_cie_sum + Lveil_cie;
2688 <                                        printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
2688 >                                        if (pict_get_z1_gsn(p,x)<0) {
2689 >                                        act_gsn=(int)(-pict_get_z1_gsn(p,x));
2690 >                                        }else{
2691 >                                        act_gsn=0;
2692 >                                        }
2693 >                                        printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %i \n",
2694                                                     i, pict_get_npix(p, x), pict_get_av_posx(p, x),
2695                                                     pict_get_ysize(p) - pict_get_av_posy(p, x),
2696                                                     pict_get_av_lum(p, x), pict_get_av_omega(p, x),
# Line 2126 | Line 2698 | if (calcfast == 1 ) {
2698                                                                                  pict_get_av_posy(p, x),
2699                                                                                  posindex_2), lum_backg, lum_task,
2700                                                     E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
2701 <                                                   ,pict_get_cached_dir(p, x2,y2)[0],pict_get_cached_dir(p, x2,y2)[1],pict_get_cached_dir(p, x2,y2)[2],pict_get_Eglare(p,x),Lveil_cie,teta
2130 <                                                  
2131 <                                                   );
2701 >                                                   ,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 );
2702                                  }
2703                          }
2704                  }
# Line 2168 | Line 2738 | if (calcfast == 1 ) {
2738          }
2739  
2740   if (calcfast == 0) {
2741 +        lum_a= E_v2/3.1415927;
2742          dgi = get_dgi(p, lum_backg, igs, posindex_2);
2743          ugr = get_ugr(p, lum_backg, igs, posindex_2);
2744 +        ugp = get_ugp (p,lum_backg, igs, posindex_2);
2745 +        ugr_exp = get_ugr_exp (p,lum_backg_cos,lum_a, igs, posindex_2);
2746 +        dgi_mod = get_dgi_mod(p, lum_a, igs, posindex_2);
2747          cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
2748 <        vcp = get_vcp(p, avlum, igs, posindex_2);
2748 >        dgr = get_dgr(p, avlum, igs, posindex_2);
2749 >        vcp = get_vcp(dgr);
2750          Lveil = get_disability(p, avlum, igs);
2751 +        if (no_glaresources == 0) {
2752 +                dgi = 0.0;
2753 +                ugr = 0.0;
2754 +                ugp = 0.0;
2755 +                ugr_exp = 0.0;
2756 +                dgi_mod = 0.0;
2757 +                cgi = 0.0;
2758 +                dgr = 0.0;
2759 +                vcp = 100.0;
2760 +        }
2761   }
2762   /* check dgp range */
2763          if (dgp <= 0.2) {
# Line 2198 | Line 2783 | if (calcfast == 0) {
2783                                  dgp =low_light_corr*dgp;
2784                                  dgp =age_corr_factor*dgp;
2785                          }
2786 +                muc_rvar_get_median(s_noposweight,lum_nopos_median);
2787 +                muc_rvar_get_median(s_posweight,lum_pos_median);
2788 +                muc_rvar_get_median(s_posweight2,lum_pos2_median);
2789 +
2790  
2791 +
2792 +        
2793                          printf
2794 <                                ("dgp,av_lum,E_v,lum_backg,E_v_dir,dgi,ugr,vcp,cgi,lum_sources,omega_sources,Lveil,Lveil_cie,band_avlum: %f %f %f %f %f %f %f %f %f %f %f %f %f %f \n",
2794 >                                ("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",
2795                                   dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
2796 <                                 lum_sources / omega_sources, omega_sources, Lveil,Lveil_cie_sum,band_avlum);
2796 >                                 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]);
2797                  } else {
2798                          if (detail_out2 == 1) {
2799  
# Line 2265 | Line 2856 | has to be re-written from scratch....
2856  
2857   /*output  */
2858          if (checkfile == 1) {
2859 +                        
2860 +        
2861                  pict_write(p, file_out);
2862          }
2863  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines