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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines