--- ray/src/util/evalglare.c 2015/08/18 15:02:53 2.2 +++ ray/src/util/evalglare.c 2017/05/18 02:25:27 2.6 @@ -1,10 +1,10 @@ #ifndef lint -static const char RCSid[] = "$Id: evalglare.c,v 2.2 2015/08/18 15:02:53 greg Exp $"; +static const char RCSid[] = "$Id: evalglare.c,v 2.6 2017/05/18 02:25:27 greg Exp $"; #endif -/* EVALGLARE V1.17 - * Evalglare Software License, Version 1.0 +/* EVALGLARE V2.00 + * Evalglare Software License, Version 2.0 * - * Copyright (c) 1995 - 2015 Fraunhofer ISE, EPFL. + * Copyright (c) 1995 - 2016 Fraunhofer ISE, EPFL. * All rights reserved. * * @@ -23,7 +23,7 @@ static const char RCSid[] = "$Id: evalglare.c,v 2.2 20 * 3. The end-user documentation included with the redistribution, * if any, must include the following acknowledgments: * "This product includes the evalglare software, - developed at Fraunhofer ISE by J. Wienold" and + developed at Fraunhofer ISE and EPFL by J. Wienold" and * "This product includes Radiance software * (http://radsite.lbl.gov/) * developed by the Lawrence Berkeley National Laboratory @@ -31,14 +31,14 @@ static const char RCSid[] = "$Id: evalglare.c,v 2.2 20 * Alternately, this acknowledgment may appear in the software itself, * if and wherever such third-party acknowledgments normally appear. * - * 4. The names "Evalglare," and "Fraunhofer ISE" must + * 4. The names "Evalglare," "EPFL" and "Fraunhofer ISE" must * not be used to endorse or promote products derived from this * software without prior written permission. For written * permission, please contact jan.wienold@epfl.ch * * 5. Products derived from this software may not be called "evalglare", * nor may "evalglare" appear in their name, without prior written - * permission of Fraunhofer ISE. + * permission of EPFL. * * Redistribution and use in source and binary forms, WITH * modification, are permitted provided that the following conditions @@ -48,7 +48,7 @@ static const char RCSid[] = "$Id: evalglare.c,v 2.2 20 * conditions 1.-5. apply * * 6. In order to ensure scientific correctness, any modification in source code imply fulfilling all following comditions: - * a) A written permission from Fraunhofer ISE. For written permission, please contact jan.wienold@ise.fraunhofer.de. + * a) A written permission from EPFL. For written permission, please contact jan.wienold@epfl.ch. * b) The permission can be applied via email and must contain the applied changes in source code and at least two example calculations, * comparing the results of the modified version and the current version of evalglare. * b) Any redistribution of a modified version of evalglare must contain following note: @@ -278,17 +278,59 @@ static const char RCSid[] = "$Id: evalglare.c,v 2.2 20 remove of age factor due to non proven statistical evidence */ -#define EVALGLARE +/* 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. +memory leakage was checked and is o.k. + */ +/* 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) + */ +/* 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 + */ +/* evalglare.c, v1.21 2016/03/15 add a new pixel-overflow-correction option -N + */ +/* 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) + */ +/* evalglare.c, v1.23 2016/04/18 fixed problem making black corners for -vth fish-eye images + */ +/* 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. +changed masking threshold to 0.05 cd/m2 + */ +/* 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. + */ +/* 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. + The calculation of the average source lumiance and solid angle was adapted to be more robust for extreme cases. + 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. + Changed glare section output when 0 glare source are found to have the same amount of columns than when glare sources are found + */ +/* 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. + */ +/* evalglare.c, v1.28 2016/07/08 small code changes in the section of calculating the E_glare (disability) for each glare source. + */ +/* evalglare.c, v1.29 2016/07/12 change threshold for the black corner to 2cd/m2. + */ +/* 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). + */ +/* 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. + */ +/* evalglare.c, v2.00 2016/11/15 add of a second fast calculation mode for annual calculations, activted by -2. Output: dgp,ugr + */ +/* evalglare.c, v2.01 2016/11/16 change of -2 option (now -2 dir_illum). External provision of the direct illuminance necessary, since the sun interpolation of daysim is causing problems in calculation of the background luminance. + */ +/* evalglare.c, v2.02 2017/02/28 change of warning message, when invalid exposure setting is found. Reason: tab removal is not in all cases the right measure - it depends which tool caused the invalid exposure entry */ + +#define EVALGLARE #define PROGNAME "evalglare" -#define VERSION "1.17 release 15.07.2015 by EPFL, J.Wienold" +#define VERSION "2.02 release 28.02.2017 by EPFL, J.Wienold" #define RELEASENAME PROGNAME " " VERSION -#include "rtio.h" -#include "platform.h" + #include "pictool.h" +#include "rtio.h" #include #include +#include "platform.h" +#include "muc_randvar.h" + char *progname; /* subroutine to add a pixel to a glare source */ @@ -564,57 +606,10 @@ double get_task_lum(pict * p, int x, int y, float r, i return av_lum; } -/* subroutine for calculation of band luminance */ -double get_band_lum(pict * p, float r, int task_color) -{ - int x_min, x_max, y_min, y_max, ix, iy, y_mid; - double r_actual, av_lum, omega_sum, act_lum; - x_max = pict_get_xsize(p) - 1; - y_max = pict_get_ysize(p) - 1; - x_min = 0; - y_min = 0; - y_mid = rint(y_max/2); - - av_lum = 0.0; - omega_sum = 0.0; - - for (iy = y_min; iy <= y_max; iy++) { - for (ix = x_min; ix <= x_max; ix++) { - -/* if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0) - continue;*/ - r_actual = - acos(DOT(pict_get_cached_dir(p, ix, y_mid), - pict_get_cached_dir(p, ix, iy))) ; - act_lum = luminance(pict_get_color(p, ix, iy)); - - if ((r_actual <= r) || (iy == y_mid) ) { - act_lum = luminance(pict_get_color(p, ix, iy)); - av_lum += pict_get_omega(p, ix, iy) * act_lum; - omega_sum += pict_get_omega(p, ix, iy); - if (task_color == 1) { - pict_get_color(p, ix, iy)[RED] = 0.0; - pict_get_color(p, ix, iy)[GRN] = act_lum / WHTEFFICACY / CIE_gf; - pict_get_color(p, ix, iy)[BLU] = 0.0; - } - } - } - } - - av_lum = av_lum / omega_sum; -/* printf("average luminance of band %f \n",av_lum);*/ -/* printf("solid angle of band %f \n",omega_sum);*/ - return av_lum; -} - - - - - /* subroutine for coloring the glare sources red is used only for explicit call of the subroutine with acol=0 , e.g. for disability glare -1: set to grey*/ @@ -626,32 +621,32 @@ int setglcolor(pict * p, int x, int y, int acol, int u l=u_r+u_g+u_b ; pcol[RED][0] = 1.0 / CIE_rf; - pcol[GRN][0] = 0.; - pcol[BLU][0] = 0.; + pcol[GRN][0] = 0.0 / CIE_gf; + pcol[BLU][0] = 0.0 / CIE_bf; - pcol[RED][1] = 0.; + pcol[RED][1] = 0.0 / CIE_rf; pcol[GRN][1] = 1.0 / CIE_gf; - pcol[BLU][1] = 0.; + pcol[BLU][1] = 0.0 / CIE_bf; - pcol[RED][2] = 0.; - pcol[GRN][2] = 0.; - pcol[BLU][2] = 1 / CIE_bf; + pcol[RED][2] = 0.15 / CIE_rf; + pcol[GRN][2] = 0.15 / CIE_gf; + pcol[BLU][2] = 0.7 / CIE_bf; - pcol[RED][3] = 1.0 / (1.0 - CIE_bf); - pcol[GRN][3] = 1.0 / (1.0 - CIE_bf); - pcol[BLU][3] = 0.; + pcol[RED][3] = .5 / CIE_rf; + pcol[GRN][3] = .5 / CIE_gf; + pcol[BLU][3] = 0.0 / CIE_bf; - pcol[RED][4] = 1.0 / (1.0 - CIE_gf); - pcol[GRN][4] = 0.; - pcol[BLU][4] = 1.0 / (1.0 - CIE_gf); + pcol[RED][4] = .5 / CIE_rf; + pcol[GRN][4] = .0 / CIE_gf; + pcol[BLU][4] = .5 / CIE_bf ; - pcol[RED][5] = 0.; - pcol[GRN][5] = 1.0 / (1.0 - CIE_rf); - pcol[BLU][5] = 1.0 / (1.0 - CIE_rf); + pcol[RED][5] = .0 / CIE_rf; + pcol[GRN][5] = .5 / CIE_gf; + pcol[BLU][5] = .5 / CIE_bf; - pcol[RED][6] = 0.; - pcol[GRN][6] = 1.0 / (1.0 - CIE_rf); - pcol[BLU][6] = 1.0 / (1.0 - CIE_rf); + pcol[RED][6] = 0.333 / CIE_rf; + pcol[GRN][6] = 0.333 / CIE_gf; + pcol[BLU][6] = 0.333 / CIE_bf; pcol[RED][999] = 1.0 / WHTEFFICACY; @@ -662,7 +657,7 @@ int setglcolor(pict * p, int x, int y, int acol, int u pcol[RED][998] = u_r /(l* CIE_rf) ; pcol[GRN][998] = u_g /(l* CIE_gf); pcol[BLU][998] = u_b /(l* CIE_bf); -/* printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l);*/ +/* printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l,WHTEFFICACY); */ icol = acol ; if ( acol == -1) {icol=999; }else{if (acol>0){icol = acol % 5 +1; @@ -675,7 +670,7 @@ int setglcolor(pict * p, int x, int y, int acol, int u pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY; pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY; pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY; - +/* 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))); */ return icol; } @@ -727,7 +722,7 @@ void add_secondary_gs(pict * p, int x, int y, float r, /* color pixel of gs */ - icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b); +/* icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b); */ } } @@ -767,9 +762,9 @@ void split_pixel_from_gs(pict * p, int x, int y, int n /* color pixel of new gs */ - icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b); +/* icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b); */ + - } /* subroutine for the calculation of the position index */ @@ -800,6 +795,11 @@ float get_posindex(pict * p, float x, float y, int pos tau = tau * deg; sigma = sigma * deg; + if (postype == 1) { +/* KIM model */ + 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)); + }else{ +/* Guth model, equation from IES lighting handbook */ posindex = exp((35.2 - 0.31889 * tau - 1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 + @@ -821,8 +821,9 @@ float get_posindex(pict * p, float x, float y, int pos posindex = 1 + fact * r; } - if (posindex > 16) + if (posindex > 16) posindex = 16; +} return posindex; @@ -1035,7 +1036,7 @@ return; } -/* subroutine for the calculation of the daylight glare index */ +/* subroutine for the calculation of the daylight glare index DGI*/ float get_dgi(pict * p, float lum_backg, int igs, int posindex_2) @@ -1062,7 +1063,36 @@ float get_dgi(pict * p, float lum_backg, int igs, int return dgi; } +/* subroutine for the calculation of the modified daylight glare index DGI_mod (Fisekis2003) + several glare sources are taken into account and summed up, original equation has no summary +*/ + +float get_dgi_mod(pict * p, float lum_a, int igs, int posindex_2) +{ + float dgi_mod, sum_glare, omega_s; + int i; + + + sum_glare = 0; + omega_s = 0; + + for (i = 0; i <= igs; i++) { + + if (pict_get_npix(p, i) > 0) { + omega_s = + pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) / + get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2); + 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)); + /* printf("i,sum_glare %i %f\n",i,sum_glare); */ + } + } + dgi_mod = 10 * log10(sum_glare); + + return dgi_mod; + +} + /* subroutine for the calculation of the daylight glare probability */ double get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3, @@ -1103,10 +1133,10 @@ get_dgp(pict * p, double E_v, int igs, double a1, doub } /* subroutine for the calculation of the visual comfort probability */ -float get_vcp(pict * p, double lum_a, int igs, int posindex_2) +float get_dgr(pict * p, double lum_a, int igs, int posindex_2) { - float vcp; - double sum_glare, dgr; + float dgr; + double sum_glare; int i, i_glare; @@ -1132,6 +1162,16 @@ float get_vcp(pict * p, double lum_a, int igs, int pos } dgr = pow(sum_glare, pow(i_glare, -0.0914)); + + + return dgr; + +} + +float get_vcp(float dgr ) +{ + float vcp; + vcp = 50 * erf((6.347 - 1.3227 * log(dgr)) / 1.414213562373) + 50; if (dgr > 750) { vcp = 0; @@ -1139,12 +1179,11 @@ float get_vcp(pict * p, double lum_a, int igs, int pos if (dgr < 20) { vcp = 100; } - - return vcp; } + /* subroutine for the calculation of the unified glare rating */ float get_ugr(pict * p, double lum_backg, int igs, int posindex_2) { @@ -1168,12 +1207,69 @@ float get_ugr(pict * p, double lum_backg, int igs, int } } ugr = 8 * log10(0.25 / lum_backg * sum_glare); - + if (sum_glare==0) { + ugr=0.0; + } + if (lum_backg<=0) { + ugr=-99.0; + } + return ugr; } +/* subroutine for the calculation of the experimental unified glare rating (Fisekis 2003)*/ +float get_ugr_exp(pict * p, double lum_backg, double lum_a, int igs, int posindex_2) +{ + float ugr_exp; + double sum_glare; + int i, i_glare; + sum_glare = 0; + i_glare = 0; + for (i = 0; i <= igs; i++) { + + if (pict_get_npix(p, i) > 0) { + i_glare = i_glare + 1; + sum_glare += + 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); + + } + } + ugr_exp =8 * log10(lum_a) + 8 * log10(sum_glare/lum_backg); + + return ugr_exp; + +} +/* subroutine for the calculation of the unified glare probability (Hirning 2014)*/ +float get_ugp(pict * p, double lum_backg, int igs, int posindex_2) +{ + float ugp; + double sum_glare; + int i, i_glare; + + + sum_glare = 0; + i_glare = 0; + for (i = 0; i <= igs; i++) { + + if (pict_get_npix(p, i) > 0) { + i_glare = i_glare + 1; + sum_glare += + pow(pict_get_av_lum(p, i) / + get_posindex(p, pict_get_av_posx(p, i), + pict_get_av_posy(p, i), posindex_2), + 2) * pict_get_av_omega(p, i); + + } + } + ugp = 0.26 * log10(0.25 / lum_backg * sum_glare); + + return ugp; + +} + + /* subroutine for the calculation of the disability glare according to Poynter */ float get_disability(pict * p, double lum_backg, int igs) { @@ -1221,8 +1317,7 @@ float get_disability(pict * p, double lum_backg, int i /* subroutine for the calculation of the cie glare index */ -float -get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2) +float get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2) { float cgi; double sum_glare; @@ -1246,10 +1341,36 @@ get_cgi(pict * p, double E_v, double E_v_dir, int igs, cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare); return cgi; +} +/* 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! */ +float get_pgsv(double E_v, double E_mask,double omega_mask,double lum_mask_av) +{ + float pgsv; + double Lb; + + Lb = (E_v-E_mask)/1.414213562373; + + pgsv =3.2*log10(lum_mask_av)-0.64*log10(omega_mask)+(0.79*log10(omega_mask)-0.61)*log10(Lb)-8.2 ; + + + return pgsv; } +/* 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! */ +float get_pgsv_sat(double E_v) +{ + float pgsv_sat; + pgsv_sat =3.3-(0.57+3.3)/pow((1+E_v/1.414213562373/1250),1.7); + + + return pgsv_sat; +} + + + + #ifdef EVALGLARE @@ -1261,26 +1382,96 @@ int main(int argc, char **argv) { #define CLINEMAX 4095 /* memory allocated for command line string */ pict *p = pict_create(); - int skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, y1, fill, yfillmax, yfillmin, - ext_vill, set_lum_max, set_lum_max2, task_color, i_splitstart, - i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2, - igs, actual_igs, icol, xt, yt, change, before_igs, sgs, splithigh,uniform_gs, - detail_out, posindex_picture, non_cos_lb, rx, ry, rmx, rmy,apply_disability,band_calc,band_color; - double lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma, - E_vl_ext, lum_max, new_lum_max, r_center, - search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle, - omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang, - l_max, lum_backg, lum_backg_cos, omega_sources, lum_sources, - lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum; - float lum_task, lum_thres, dgi, vcp, cgi, ugr, limit, + pict *pm = pict_create(); + int skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, x1,y1, fill, yfillmax, yfillmin, + ext_vill, set_lum_max, set_lum_max2, img_corr,x_disk,y_disk,task_color, i_splitstart,zones,act_gsn,splitgs, + i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,x_zone,y_zone, i_z1, i_z2, + igs, actual_igs, lastpixelwas_gs, icol, xt, yt, change,checkpixels, before_igs, sgs, splithigh,uniform_gs,x_max, y_max,y_mid, + detail_out, posindex_picture, non_cos_lb, rx, ry, rmx,rmy,apply_disability,band_calc,band_color,masking,i_mask,no_glaresources,force; + 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, + 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, + search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,r_actual,lum_actual,dir_ill, + omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang, angle_z1, angle_z2,per_95_band,per_75_band,pos, + 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, + lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum, + lum_mask[1],lum_mask_av,lum_mask_std[1],lum_mask_median[1],omega_mask,bbox[2], + lum_band[1],lum_band_av,lum_band_std[1],lum_band_median[1],omega_band,bbox_band[2], + lum_z1[1],lum_z1_av,lum_z1_std[1],lum_z1_median[1],omega_z1,bbox_z1[2], + lum_z2[1],lum_z2_av,lum_z2_std[1],lum_z2_median[1],omega_z2,bbox_z2[2], + lum_pos[1],lum_nopos_median[1],lum_pos_median[1],lum_pos2_median[1],lum_pos_mean,lum_pos2_mean; + float lum_task, lum_thres, dgi, vcp, cgi, ugr, limit, dgr, abs_max, Lveil; - char file_out[500], file_out2[500], version[500]; + char maskfile[500],file_out[500], file_out2[500], version[500]; char *cline; VIEW userview = STDVIEW; int gotuserview = 0; + struct muc_rvar* s_mask; + s_mask = muc_rvar_create(); + muc_rvar_set_dim(s_mask, 1); + muc_rvar_clear(s_mask); + struct muc_rvar* s_band; + s_band = muc_rvar_create(); + muc_rvar_set_dim(s_band, 1); + muc_rvar_clear(s_band); + struct muc_rvar* s_z1; + s_z1 = muc_rvar_create(); + muc_rvar_set_dim(s_z1, 1); + muc_rvar_clear(s_z1); + struct muc_rvar* s_z2; + s_z2 = muc_rvar_create(); + muc_rvar_set_dim(s_z2, 1); + muc_rvar_clear(s_z2); + + struct muc_rvar* s_noposweight; + s_noposweight = muc_rvar_create(); + muc_rvar_set_dim(s_noposweight, 1); + muc_rvar_clear(s_noposweight); + + struct muc_rvar* s_posweight; + s_posweight = muc_rvar_create(); + muc_rvar_set_dim(s_posweight, 1); + muc_rvar_clear(s_posweight); + + struct muc_rvar* s_posweight2; + s_posweight2 = muc_rvar_create(); + muc_rvar_set_dim(s_posweight2, 1); + muc_rvar_clear(s_posweight2); + /*set required user view parameters to invalid values*/ - uniform_gs = 0; + dir_ill=0.0; + delta_E=0.0; + no_glaresources=0; + n_corner_px=0; + zero_corner_px=0; + force=0; + dist=0.0; + u_r=0.33333; + u_g=0.33333; + u_b=0.33333; + band_angle=0; + angle_z1=0; + angle_z2=0; + x_zone=0; + y_zone=0; + per_75_z2=0; + per_95_z2=0; + lum_pos_mean=0; + lum_pos2_mean=0; + lum_band_av = 0.0; + omega_band = 0.0; + pgsv = 0.0 ; + pgsv_sat = 0.0 ; + E_v_mask = 0.0; + lum_z1_av = 0.0; + omega_z1 = 0.0; + lum_z2_av = 0.0; + omega_z2 = 0.0 ; + i_z1 = 0; + i_z2 = 0; + zones = 0; + lum_z2_av = 0.0; + uniform_gs = 0; apply_disability = 0; disability_thresh = 0; Lveil_cie_sum=0.0; @@ -1300,6 +1491,9 @@ int main(int argc, char **argv) omegat = 0.0; yt = 0; xt = 0; + x_disk=0; + y_disk=0; + angle_disk=0; yfillmin = 0; yfillmax = 0; cut_view = 0; @@ -1332,15 +1526,22 @@ int main(int argc, char **argv) limit = 50000.0; set_lum_max = 0; set_lum_max2 = 0; + img_corr=0; abs_max = 0; progname = argv[0]; E_v_contr = 0.0; - strcpy(version, "1.15 release 05.02.2015 by J.Wienold"); + strcpy(version, "1.19 release 09.12.2015 by J.Wienold"); low_light_corr=1.0; output = 0; calc_vill = 0; band_avlum = -99; band_calc = 0; + masking = 0; + lum_mask[0]=0.0; + lum_mask_av=0.0; + omega_mask=0.0; + i_mask=0; + actual_igs=0; /* command line for output picture*/ cline = (char *) malloc(CLINEMAX+1); @@ -1382,6 +1583,13 @@ int main(int argc, char **argv) printf("age factor not supported any more \n"); exit(1); break; + case 'A': + masking = 1; + detail_out = 1; + strcpy(maskfile, argv[++i]); + pict_read(pm,maskfile ); + + break; case 'b': lum_thres = atof(argv[++i]); break; @@ -1401,6 +1609,9 @@ int main(int argc, char **argv) case 's': sgs = 1; break; + case 'f': + force = 1; + break; case 'k': apply_disability = 1; disability_thresh = atof(argv[++i]); @@ -1409,6 +1620,10 @@ int main(int argc, char **argv) case 'p': posindex_picture = 1; break; + case 'P': + posindex_2 = atoi(argv[++i]); + break; + case 'y': splithigh = 1; @@ -1439,6 +1654,7 @@ int main(int argc, char **argv) detail_out2 = 1; break; case 'm': + img_corr = 1; set_lum_max = 1; lum_max = atof(argv[++i]); new_lum_max = atof(argv[++i]); @@ -1447,12 +1663,24 @@ int main(int argc, char **argv) break; case 'M': + img_corr = 1; set_lum_max2 = 1; lum_max = atof(argv[++i]); E_vl_ext = atof(argv[++i]); strcpy(file_out2, argv[++i]); /* printf("max lum set to %f\n",new_lum_max);*/ break; + case 'N': + img_corr = 1; + set_lum_max2 = 2; + x_disk = atoi(argv[++i]); + y_disk = atoi(argv[++i]); + angle_disk = atof(argv[++i]); + E_vl_ext = atof(argv[++i]); + strcpy(file_out2, argv[++i]); +/* printf("max lum set to %f\n",new_lum_max);*/ + break; + case 'n': non_cos_lb = 0; break; @@ -1472,6 +1700,22 @@ int main(int argc, char **argv) omegat = atof(argv[++i]); task_color = 1; break; + case 'l': + zones = 1; + x_zone = atoi(argv[++i]); + y_zone = atoi(argv[++i]); + angle_z1 = atof(argv[++i]); + angle_z2 = -1; + break; + case 'L': + zones = 2; + x_zone = atoi(argv[++i]); + y_zone = atoi(argv[++i]); + angle_z1 = atof(argv[++i]); + angle_z2 = atof(argv[++i]); + break; + + case 'B': band_calc = 1; band_angle = atof(argv[++i]); @@ -1508,6 +1752,10 @@ int main(int argc, char **argv) case '1': output = 1; break; + case '2': + output = 2; + dir_ill = atof(argv[++i]); + break; case 'v': if (argv[i][2] == '\0') { @@ -1539,10 +1787,21 @@ int main(int argc, char **argv) /*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/ -if (output == 1 && ext_vill == 1) { +if (output == 1 && ext_vill == 1 ) { calcfast=1; } + +if (output == 2 && ext_vill == 1 ) { + calcfast=2; + } + +/*masking and zoning cannot be applied at the same time*/ +if (masking ==1 && zones >0) { + fprintf(stderr, " masking and zoning cannot be activated at the same time!\n"); + exit(1); +} + /* read picture file */ if (i == argc) { SET_FILE_BINARY(stdin); @@ -1576,18 +1835,38 @@ if (output == 1 && ext_vill == 1) { return EXIT_FAILURE; } + + /* fprintf(stderr,"Picture read!");*/ /* command line for output picture*/ pict_set_comment(p, cline); +/* several checks */ if (p->view.type == VT_PAR) { - fprintf(stderr, "wrong view type! must not be parallel ! \n"); + fprintf(stderr, "error: wrong view type! must not be parallel ! \n"); exit(1); } + if (p->view.type == VT_PER) { + fprintf(stderr, "warning: image has perspective view type specified in header ! \n"); + } + if (masking == 1) { + + if (!pict_get_xsize(p)==pict_get_xsize(pm) || !pict_get_ysize(p)==pict_get_ysize(pm)) { + fprintf(stderr, "error: masking image has other resolution than main image ! \n"); + fprintf(stderr, "size must be identical \n"); + printf("resolution main image : %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); + printf("resolution masking image : %dx%d\n",pict_get_xsize(pm),pict_get_ysize(pm)); + exit(1); + + } + + } +/* printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */ + /* Check size of search radius */ rmx = (pict_get_xsize(p) / 2); rmy = (pict_get_ysize(p) / 2); @@ -1632,7 +1911,10 @@ if (output == 1 && ext_vill == 1) { avlum = 0.0; pict_new_gli(p); igs = 0; + pict_get_z1_gsn(p,igs) = 0; + pict_get_z2_gsn(p,igs) = 0; + /* cut out GUTH field of view and exit without glare evaluation */ if (cut_view==2) { if (cut_view_type==1) { @@ -1699,9 +1981,18 @@ if (cut_view==2) { if (calcfast == 0) { for (x = 0; x < pict_get_xsize(p); x++) for (y = 0; y < pict_get_ysize(p); y++) { + lum = luminance(pict_get_color(p, x, y)); + dist=sqrt((x-rmx)*(x-rmx)+(y-rmy)*(y-rmy)); + if (dist>((rmx+rmy)/2)) { + n_corner_px=n_corner_px+1; + if (lum<7.0) { + zero_corner_px=zero_corner_px+1; + } + } + + if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) { if (pict_is_validpixel(p, x, y)) { - lum = luminance(pict_get_color(p, x, y)); if (fill == 1 && y >= yfillmax) { y1 = y - 1; lum = luminance(pict_get_color(p, x, y1)); @@ -1714,70 +2005,106 @@ if (cut_view==2) { abs_max = lum; } /* set luminance restriction, if -m is set */ - if (set_lum_max == 1 && lum >= lum_max) { - pict_get_color(p, x, y)[RED] = new_lum_max / 179.0; - pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0; - pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0; -/* printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */ - lum = luminance(pict_get_color(p, x, y)); + if (img_corr == 1 ) { + if (set_lum_max == 1 && lum >= lum_max) { + pict_get_color(p, x, y)[RED] = new_lum_max / 179.0; + pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0; + pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0; +/* printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */ + lum = luminance(pict_get_color(p, x, y)); + } + if (set_lum_max2 == 1 && lum >= lum_max) { - } - if (set_lum_max2 == 1 && lum >= lum_max) { + E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum; + omega_cos_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1; + } + if (set_lum_max2 == 2 ) { + r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2; + if (x_disk == x && y_disk==y ) r_actual=0.0; - E_v_contr += - DOT(p->view.vdir, - pict_get_cached_dir(p, x, - y)) * pict_get_omega(p, - x, - y) - * lum; - omega_cos_contr += - DOT(p->view.vdir, - pict_get_cached_dir(p, x, - y)) * pict_get_omega(p, - x, - y) - * 1; + act_lum = luminance(pict_get_color(p, x, y)); + + if (r_actual <= angle_disk) { + E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum; + omega_cos_contr += DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1; + } + + + + } } + om = pict_get_omega(p, x, y); + sang += om; + E_v += DOT(p->view.vdir, pict_get_cached_dir(p, x,y)) * om *lum; + avlum += om * lum; + pos=get_posindex(p, x, y,posindex_2); + lum_pos[0]=lum; + muc_rvar_add_sample(s_noposweight,lum_pos); + lum_pos[0]=lum/pos; + lum_pos_mean+=lum_pos[0]*om; + muc_rvar_add_sample(s_posweight, lum_pos); + lum_pos[0]=lum_pos[0]/pos; + lum_pos2_mean+=lum_pos[0]*om; + muc_rvar_add_sample(s_posweight2,lum_pos); - sang += pict_get_omega(p, x, y); - E_v += - DOT(p->view.vdir, - pict_get_cached_dir(p, x, - y)) * pict_get_omega(p, x, - y) * - lum; - avlum += - pict_get_omega(p, x, - y) * luminance(pict_get_color(p, x, - y)); + + } else { pict_get_color(p, x, y)[RED] = 0.0; pict_get_color(p, x, y)[GRN] = 0.0; pict_get_color(p, x, y)[BLU] = 0.0; } - } + }else { + pict_get_color(p, x, y)[RED] = 0.0; + pict_get_color(p, x, y)[GRN] = 0.0; + pict_get_color(p, x, y)[BLU] = 0.0; + + } } - if (set_lum_max2 == 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) { +/* check if image has black corners AND a perspective view */ + if (zero_corner_px/n_corner_px > 0.70 && (p->view.type == VT_PER) ) { + printf (" corner pixels are to %f %% black! \n",zero_corner_px/n_corner_px*100); + printf("error! The image has black corners AND header contains a perspective view type definition !!!\n") ; + + if (force==0) { + printf("stopping...!!!!\n") ; + + exit(1); + } + } + + + + + lum_pos_mean= lum_pos_mean/sang; + lum_pos2_mean= lum_pos2_mean/sang; + + if (set_lum_max2 >= 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) { + lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr; + if (set_lum_max2 == 2 && lum_ideal >= 2e9) { + printf("warning! luminance of replacement pixels would be larger than 2e9 cd/m2. Value set to 2e9cd/m2!\n") ; + } + if (lum_ideal > 0 && lum_ideal < setvalue) { setvalue = lum_ideal; } - printf - ("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f %f %f %f %f\n", + printf("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f %f %f %f %f\n", lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr); - + for (x = 0; x < pict_get_xsize(p); x++) for (y = 0; y < pict_get_ysize(p); y++) { if (pict_get_hangle (p, x, y, p->view.vdir, p->view.vup, &ang)) { if (pict_is_validpixel(p, x, y)) { lum = luminance(pict_get_color(p, x, y)); + + if (set_lum_max2 == 1 && lum >= lum_max) { pict_get_color(p, x, y)[RED] = @@ -1787,13 +2114,31 @@ if (cut_view==2) { pict_get_color(p, x, y)[BLU] = setvalue / 179.0; + }else{ if(set_lum_max2 ==2 ) { + r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2; + if (x_disk == x && y_disk==y ) r_actual=0.0; + + if (r_actual <= angle_disk) { + + pict_get_color(p, x, y)[RED] = + setvalue / 179.0; + pict_get_color(p, x, y)[GRN] = + setvalue / 179.0; + pict_get_color(p, x, y)[BLU] = + setvalue / 179.0; + + } + + } + } } } } + pict_write(p, file_out2); - + exit(1); } if (set_lum_max == 1) { pict_write(p, file_out2); @@ -1853,10 +2198,11 @@ if (cut_view==1) { lum_source = lum_thres; } - /* printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */ +/* printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */ /* first glare source scan: find primary glare sources */ - for (x = 0; x < pict_get_xsize(p); x++) + for (x = 0; x < pict_get_xsize(p); x++) { +/* lastpixelwas_gs=0; */ for (y = 0; y < pict_get_ysize(p); y++) { if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) { if (pict_is_validpixel(p, x, y)) { @@ -1865,44 +2211,73 @@ if (cut_view==1) { if (act_lum >lum_total_max) { lum_total_max=act_lum; } - - actual_igs = - find_near_pgs(p, x, y, max_angle, 0, igs, 1); +/* speed improvement first scan: when last pixel was glare source, then use this glare source number instead of searching. + has been deactivated with v1.25 + + if (lastpixelwas_gs==0 || search_pix <= 1.0 ) { */ + actual_igs = find_near_pgs(p, x, y, max_angle, 0, igs, 1); +/* }*/ if (actual_igs == 0) { igs = igs + 1; pict_new_gli(p); pict_get_lum_min(p, igs) = HUGE_VAL; pict_get_Eglare(p,igs) = 0.0; pict_get_Dglare(p,igs) = 0.0; + pict_get_z1_gsn(p,igs) = 0; + pict_get_z2_gsn(p,igs) = 0; actual_igs = igs; + } - icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); +/* no coloring any more here icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */ pict_get_gsn(p, x, y) = actual_igs; pict_get_pgs(p, x, y) = 1; add_pixel_to_gs(p, x, y, actual_igs); +/* lastpixelwas_gs=actual_igs; */ } else { pict_get_pgs(p, x, y) = 0; + /* lastpixelwas_gs=0;*/ } } } - } + } + } /* pict_write(p,"firstscan.pic"); */ -if (calcfast == 1 ) { + + + +if (calcfast ==1 || search_pix <= 1.0 || calcfast == 2 ) { skip_second_scan=1; } + /* second glare source scan: combine glare sources facing each other */ change = 1; + i = 0; while (change == 1 && skip_second_scan==0) { change = 0; - for (x = 0; x < pict_get_xsize(p); x++) + i = i+1; + for (x = 0; x < pict_get_xsize(p); x++) { for (y = 0; y < pict_get_ysize(p); y++) { - before_igs = pict_get_gsn(p, x, y); if (pict_get_hangle (p, x, y, p->view.vdir, p->view.vup, &ang)) { - if (pict_is_validpixel(p, x, y) && before_igs > 0) { + checkpixels=1; + before_igs = pict_get_gsn(p, x, y); + +/* speed improvement: search for other nearby glare sources only, when at least one adjacend pixel has another glare source number */ + if (x >1 && x1 && y 0 && checkpixels==1) { actual_igs = find_near_pgs(p, x, y, max_angle, before_igs, igs, 1); @@ -1911,16 +2286,15 @@ if (calcfast == 1 ) { } if (before_igs > 0) { actual_igs = pict_get_gsn(p, x, y); - icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); + /* icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */ } } } } /* pict_write(p,"secondscan.pic");*/ + }} - } - /*smoothing: add secondary glare sources */ i_max = igs; if (sgs == 1) { @@ -1977,6 +2351,9 @@ if (calcfast == 1 ) { if (i_splitstart == (igs + 1)) { igs = igs + 1; pict_new_gli(p); + pict_get_z1_gsn(p,igs) = 0; + pict_get_z2_gsn(p,igs) = 0; + pict_get_Eglare(p,igs) = 0.0; pict_get_Dglare(p,igs) = 0.0; pict_get_lum_min(p, igs) = @@ -1990,6 +2367,9 @@ if (calcfast == 1 ) { if (i_split == 0) { igs = igs + 1; pict_new_gli(p); + pict_get_z1_gsn(p,igs) = 0; + pict_get_z2_gsn(p,igs) = 0; + pict_get_Eglare(p,igs) = 0.0; pict_get_Dglare(p,igs) = 0.0; pict_get_lum_min(p, igs) = @@ -2026,9 +2406,9 @@ if (calcfast == 1 ) { if (before_igs > 0) { actual_igs = pict_get_gsn(p, x, y); - icol = + /* icol = setglcolor(p, x, y, - actual_igs, uniform_gs, u_r, u_g , u_b); + actual_igs, uniform_gs, u_r, u_g , u_b);*/ } } @@ -2043,48 +2423,256 @@ if (calcfast == 1 ) { } } -/* calculation of direct vertical illuminance for CGI and for disability glare*/ +/* calculation of direct vertical illuminance for CGI and for disability glare, coloring glare sources*/ - - if (calcfast == 0) { + if (calcfast == 0 || calcfast == 2) { for (x = 0; x < pict_get_xsize(p); x++) for (y = 0; y < pict_get_ysize(p); y++) { if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) { if (pict_is_validpixel(p, x, y)) { if (pict_get_gsn(p, x, y) > 0) { - pict_get_Eglare(p,pict_get_gsn(p, x, y)) += - DOT(p->view.vdir, pict_get_cached_dir(p, x, y)) - * pict_get_omega(p, x, y) - * luminance(pict_get_color(p, x, y)); - E_v_dir += - DOT(p->view.vdir, pict_get_cached_dir(p, x, y)) - * pict_get_omega(p, x, y) - * luminance(pict_get_color(p, x, y)); + actual_igs = pict_get_gsn(p, x, y); + 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)); + pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E; + E_v_dir = E_v_dir + delta_E; + setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); } } } } lum_backg_cos = (E_v - E_v_dir) / 3.1415927; - lum_backg = lum_backg_cos; + } -/* calc of band luminance if applied */ +/* calc of band luminance distribution if applied */ if (band_calc == 1) { - band_avlum=get_band_lum( p, band_angle, band_color); - } + x_max = pict_get_xsize(p) - 1; + y_max = pict_get_ysize(p) - 1; + y_mid = (int)(y_max/2); + for (x = 0; x <= x_max; x++) + for (y = 0; y <= y_max; y++) { + if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) { + if (pict_is_validpixel(p, x, y)) { + r_actual = acos(DOT(pict_get_cached_dir(p, x, y_mid), pict_get_cached_dir(p, x, y))) ; + act_lum = luminance(pict_get_color(p, x, y)); + + if ((r_actual <= band_angle) || (y == y_mid) ) { + lum_band[0]=luminance(pict_get_color(p, x, y)); + muc_rvar_add_sample(s_band, lum_band); + act_lum = luminance(pict_get_color(p, x, y)); + lum_band_av += pict_get_omega(p, x, y) * act_lum; + omega_band += pict_get_omega(p, x, y); + if (band_color == 1) { + pict_get_color(p, x, y)[RED] = 0.0; + pict_get_color(p, x, y)[GRN] = act_lum / WHTEFFICACY / CIE_gf; + pict_get_color(p, x, y)[BLU] = 0.0; + } + } + }}} + lum_band_av=lum_band_av/omega_band; + muc_rvar_get_vx(s_band,lum_band_std); + muc_rvar_get_percentile(s_band,lum_band_median,0.75); + per_75_band=lum_band_median[0]; + muc_rvar_get_percentile(s_band,lum_band_median,0.95); + per_95_band=lum_band_median[0]; + muc_rvar_get_median(s_band,lum_band_median); + muc_rvar_get_bounding_box(s_band,bbox_band); + + 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] ); +/* av_lum = av_lum / omega_sum; */ +/* printf("average luminance of band %f \n",av_lum);*/ +/* printf("solid angle of band %f \n",omega_sum);*/ + } + /*printf("total number of glare sources: %i\n",igs); */ lum_sources = 0; omega_sources = 0; + i=0; for (x = 0; x <= igs; x++) { + if (pict_get_npix(p, x) > 0) { lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x); omega_sources += pict_get_av_omega(p, x); + i=i+1; + }} + + if (sang == omega_sources) { + lum_backg =avlum; + } else { + lum_backg =(sang * avlum - lum_sources) / (sang - omega_sources); } - if (non_cos_lb == 1) { - lum_backg = - (sang * avlum - lum_sources) / (sang - omega_sources); + + + if (i == 0) { + lum_sources =0.0; + } else { lum_sources=lum_sources/omega_sources; } + + + + if (non_cos_lb == 0) { + lum_backg = lum_backg_cos; + } + +/* file writing NOT here + if (checkfile == 1) { + pict_write(p, file_out); + } +*/ + /* print detailed output */ + + + +/* masking */ + + if ( masking == 1 ) { + + + for (x = 0; x < pict_get_xsize(p); x++) + for (y = 0; y < pict_get_ysize(p); y++) { + if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) { + if (pict_is_validpixel(p, x, y)) { + if (luminance(pict_get_color(pm, x, y))>0.1) { +/* printf ("hallo %i %i %f \n",x,y,luminance(pict_get_color(pm, x, y)));*/ + i_mask=i_mask+1; + lum_mask[0]=luminance(pict_get_color(p, x, y)); + /* printf ("%f \n",lum_mask[0]);*/ + muc_rvar_add_sample(s_mask, lum_mask); + omega_mask += pict_get_omega(p, x, y); + lum_mask_av += pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y)); + 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)); + + pict_get_color(pm, x, y)[RED] = luminance(pict_get_color(p, x, y))/179.0; + pict_get_color(pm, x, y)[GRN] = luminance(pict_get_color(p, x, y))/179.0; + pict_get_color(pm, x, y)[BLU] = luminance(pict_get_color(p, x, y))/179.0; + + } else { + pict_get_color(p, x, y)[RED] = 0.0 ; + pict_get_color(p, x, y)[GRN] = 0.0 ; + pict_get_color(p, x, y)[BLU] = 0.0 ; + } + } + } + } +/* Average luminance in masking area (weighted by solid angle) */ + lum_mask_av=lum_mask_av/omega_mask; + muc_rvar_get_vx(s_mask,lum_mask_std); + muc_rvar_get_percentile(s_mask,lum_mask_median,0.75); + per_75_mask=lum_mask_median[0]; + muc_rvar_get_percentile(s_mask,lum_mask_median,0.95); + per_95_mask=lum_mask_median[0]; + muc_rvar_get_median(s_mask,lum_mask_median); + muc_rvar_get_bounding_box(s_mask,bbox); +/* PSGV only why masking of window is applied! */ + pgsv = get_pgsv(E_v, E_v_mask, omega_mask, lum_mask_av); + pgsv_sat =get_pgsv_sat(E_v); + if (detail_out == 1) { + + 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 ); + + } + + } + +/* zones */ + + if ( zones > 0 ) { + + for (x = 0; x < pict_get_xsize(p); x++) + for (y = 0; y < pict_get_ysize(p); y++) { + if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) { + if (pict_is_validpixel(p, x, y)) { + + + r_actual =acos(DOT(pict_get_cached_dir(p, x, y), pict_get_cached_dir(p, x_zone, y_zone))) * 2; + lum_actual = luminance(pict_get_color(p, x, y)); + act_gsn=pict_get_gsn(p, x, y); + /* 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));*/ + +/*zone1 */ + if (r_actual <= angle_z1) { + i_z1=i_z1+1; + lum_z1[0]=lum_actual; + muc_rvar_add_sample(s_z1, lum_z1); + omega_z1 += pict_get_omega(p, x, y); + lum_z1_av += pict_get_omega(p, x, y)* lum_actual; + setglcolor(p,x,y,1,1 , 0.66, 0.01 ,0.33); +/*check if separation gsn already exist */ + + if (act_gsn > 0){ + if (pict_get_z1_gsn(p,act_gsn) == 0) { + pict_new_gli(p); + igs = igs + 1; + pict_get_z1_gsn(p,act_gsn) = igs*1.000; + pict_get_z1_gsn(p,igs) = -1.0; + pict_get_z2_gsn(p,igs) = -1.0; + splitgs=(int)(pict_get_z1_gsn(p,act_gsn)); +/* 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)); */ + } + splitgs=(int)(pict_get_z1_gsn(p,act_gsn)); + /* printf("splitgs%i \n",splitgs); */ + split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b); + } + } +/*zone2 */ + + if (r_actual > angle_z1 && r_actual<= angle_z2 ) { + + i_z2=i_z2+1; + lum_z2[0]=lum_actual; + muc_rvar_add_sample(s_z2, lum_z2); + omega_z2 += pict_get_omega(p, x, y); + lum_z2_av += pict_get_omega(p, x, y)* lum_actual; + setglcolor(p,x,y,1,1 , 0.65, 0.33 ,0.02); +/* 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)); +*/ if (act_gsn > 0){ + if (pict_get_z2_gsn(p,act_gsn) == 0) { + pict_new_gli(p); + igs = igs + 1; + pict_get_z2_gsn(p,act_gsn) = igs*1.000; + pict_get_z1_gsn(p,igs) = -2.0; + pict_get_z2_gsn(p,igs) = -2.0; +/* printf ("Glare source %i gets glare source %i in zone 2 \n",act_gsn,igs ); */ + } + splitgs=(int)(pict_get_z2_gsn(p,act_gsn)); +/* printf("splitgs %i \n",splitgs);*/ + split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b); + } + } + + } + } } +/* Average luminance in zones (weighted by solid angle), calculate percentiles, median min and max in zones */ + if (zones == 2 ) { + lum_z2_av=lum_z2_av/omega_z2; + muc_rvar_get_vx(s_z2,lum_z2_std); + muc_rvar_get_percentile(s_z2,lum_z2_median,0.75); + per_75_z2=lum_z2_median[0]; + muc_rvar_get_percentile(s_z2,lum_z2_median,0.95); + per_95_z2=lum_z2_median[0]; + muc_rvar_get_median(s_z2,lum_z2_median); + muc_rvar_get_bounding_box(s_z2,bbox_z2); + } + lum_z1_av=lum_z1_av/omega_z1; + muc_rvar_get_vx(s_z1,lum_z1_std); + muc_rvar_get_percentile(s_z1,lum_z1_median,0.75); + per_75_z1=lum_z1_median[0]; + muc_rvar_get_percentile(s_z1,lum_z1_median,0.95); + per_95_z1=lum_z1_median[0]; + muc_rvar_get_median(s_z1,lum_z1_median); + muc_rvar_get_bounding_box(s_z1,bbox_z1); + if (detail_out == 1) { + + 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] ); + + if (zones == 2 ) { + + 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] ); + } } + + } + i = 0; for (x = 0; x <= igs; x++) { /* resorting glare source numbers */ @@ -2092,16 +2680,17 @@ if (calcfast == 1 ) { i = i + 1; } } + no_glaresources=i; +/* glare sources */ + if (detail_out == 1) { - printf - ("%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", + ("%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", i); if (i == 0) { - printf("%i %f %f %f %f %.10f %f %f %f %f %f %f\n", i, 0.0, 0.0, - 0.0, 0.0, 0.0, 0.0, lum_backg, lum_task, E_v, E_v_dir, - abs_max); + 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, + 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); } else { i = 0; @@ -2121,7 +2710,12 @@ if (calcfast == 1 ) { Lveil_cie =0 ; } Lveil_cie_sum = Lveil_cie_sum + Lveil_cie; - printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", + if (pict_get_z1_gsn(p,x)<0) { + act_gsn=(int)(-pict_get_z1_gsn(p,x)); + }else{ + act_gsn=0; + } + printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %i \n", i, pict_get_npix(p, x), pict_get_av_posx(p, x), pict_get_ysize(p) - pict_get_av_posy(p, x), pict_get_av_lum(p, x), pict_get_av_omega(p, x), @@ -2129,9 +2723,7 @@ if (calcfast == 1 ) { pict_get_av_posy(p, x), posindex_2), lum_backg, lum_task, E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927 - ,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 - - ); + ,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 ); } } } @@ -2171,11 +2763,26 @@ if (calcfast == 1 ) { } if (calcfast == 0) { + lum_a= E_v2/3.1415927; dgi = get_dgi(p, lum_backg, igs, posindex_2); ugr = get_ugr(p, lum_backg, igs, posindex_2); + ugp = get_ugp (p,lum_backg, igs, posindex_2); + ugr_exp = get_ugr_exp (p,lum_backg_cos,lum_a, igs, posindex_2); + dgi_mod = get_dgi_mod(p, lum_a, igs, posindex_2); cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2); - vcp = get_vcp(p, avlum, igs, posindex_2); + dgr = get_dgr(p, avlum, igs, posindex_2); + vcp = get_vcp(dgr); Lveil = get_disability(p, avlum, igs); + if (no_glaresources == 0) { + dgi = 0.0; + ugr = 0.0; + ugp = 0.0; + ugr_exp = 0.0; + dgi_mod = 0.0; + cgi = 0.0; + dgr = 0.0; + vcp = 100.0; + } } /* check dgp range */ if (dgp <= 0.2) { @@ -2201,11 +2808,17 @@ if (calcfast == 0) { dgp =low_light_corr*dgp; dgp =age_corr_factor*dgp; } + muc_rvar_get_median(s_noposweight,lum_nopos_median); + muc_rvar_get_median(s_posweight,lum_pos_median); + muc_rvar_get_median(s_posweight2,lum_pos2_median); + + + printf - ("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", + ("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", dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi, - lum_sources / omega_sources, omega_sources, Lveil,Lveil_cie_sum,band_avlum); + 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]); } else { if (detail_out2 == 1) { @@ -2236,8 +2849,15 @@ if (calcfast == 0) { if (E_vl_ext < 1000) { low_light_corr=1.0*exp(0.024*E_vl_ext-4)/(1+exp(0.024*E_vl_ext-4)); } else {low_light_corr=1.0 ;} dgp =low_light_corr*dgp; - dgp =age_corr_factor*dgp; - printf("%f\n", dgp); + + if (calcfast == 2) { + + lum_backg_cos=(E_vl_ext-dir_ill)/3.1415927; + ugr = get_ugr(p, lum_backg_cos, igs, posindex_2); + printf("%f %f \n", dgp,ugr); + }else{ + printf("%f\n", dgp); + } } @@ -2268,6 +2888,8 @@ has to be re-written from scratch.... /*output */ if (checkfile == 1) { + + pict_write(p, file_out); }