ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/radiance/ray/src/util/evalglare.c

Comparing ray/src/util/evalglare.c (file contents):
Revision 2.2 by greg, Tue Aug 18 15:02:53 2015 UTC vs.
Revision 2.10 by greg, Tue Nov 27 18:08:24 2018 UTC

#	Line 1 | Line 1
1		#ifndef lint
2		static const char RCSid[] = "$Id$";
3		#endif
4	<	/* EVALGLARE V1.17
5	<	* Evalglare Software License, Version 1.0
4	>	/* EVALGLARE V2.08
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 23 | Line 23 | static const char RCSid[] = "$Id$";
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 31 | Line 31 | static const char RCSid[] = "$Id$";
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 48 | Line 48 | static const char RCSid[] = "$Id$";
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 278 | Line 278 | static const char RCSid[] = "$Id$";
278		remove of age factor due to non proven statistical evidence
279		*/
280
281	<	#define EVALGLARE
281	>	/* evalglare.c, v1.18 2015/11/10 add masking, UDP, PSGV, DGI_mod, UGR_exp and zoning. For zoning and band luminance, a statistical package from C. Reetz was added in order to derive median, std_dev, percentiles.
282	>	memory leakage was checked and is o.k.
283	>	*/
284
285	+	/* evalglare.c, v1.19 2015/12/04 sorting subroutines in statistical package changed by. C. Reetz, adding Position index weighted Luminance evaluation (mean and median of entire image, only in detailed output available)
286	+	*/
287	+	/* evalglare.c, v1.20 2015/12/21 removing -ansi gcc-option in makefile for the standalone version, adding DGR as output since VCP equation seems not to match with original publication
288	+	*/
289	+	/* evalglare.c, v1.21 2016/03/15 add a new pixel-overflow-correction option -N
290	+	*/
291	+	/* evalglare.c, v1.22 2016/03/22 fixed problem using -N option and angular distance calculation, when actual pixel equals disk-center (arccos returned nan)
292	+	*/
293	+	/* evalglare.c, v1.23 2016/04/18 fixed problem making black corners for -vth fish-eye images
294	+	*/
295	+	/* evalglare.c, v1.24 2016/04/26 significant speed improvement for 2nd glare source scan: for other glare sources will not be searched any more, when the 8 surrounding pixels have the same glare source number! Smaller speed improvement in the first glare source scan: remembering if pixel before was also a glare source, then no search for surrounding glare sources is applied.
296	+	changed masking threshold to 0.05 cd/m2
297	+	*/
298	+	/* evalglare.c, v1.25 2016/04/27 removed the change on the first glare source scan: in few cases glare sources are not merged together in the same way as before. Although the difference is marginal, this algorithm was removed in order to be consistent to the existing results.
299	+	*/
300	+	/* evalglare.c, v1.26 2016/04/28 removed the bug Lb getting nan in case all glare source pixels are above the peak extraction limit.
301	+	The calculation of the average source lumiance and solid angle was adapted to be more robust for extreme cases.
302	+	In case no glare source is found, the values of the glare metrics, relying only on glare sources, is set to zero (100 for VCP) to avoid nan and -inf in the output.
303	+	Changed glare section output when 0 glare source are found to have the same amount of columns than when glare sources are found
304	+	*/
305	+	/* evalglare.c, v1.27 2016/06/13 include a warning, if -vtv is in the header. Check, if the corners are black AND -vtv is used ->stopping (stopping can be avoided by using the forcing option -f ). Check, if an invalid exposure is in the header -> stopping in any case.
306	+	*/
307	+	/* evalglare.c, v1.28 2016/07/08 small code changes in the section of calculating the E_glare (disability) for each glare source.
308	+	*/
309	+	/* evalglare.c, v1.29 2016/07/12 change threshold for the black corner to 2cd/m2.
310	+	*/
311	+	/* evalglare.c, v1.30 2016/07/28 change zonal output: If just one zone is specified, only one zone shows up in the output (bug removal).
312	+	*/
313	+	/* evalglare.c, v1.31 2016/08/02  bug removal: default output did not calculate the amout of glare sources before and therefore no_glaresources was set to zero causing dgi,ugr being set to zero as well. Now renumbering of the glare sources and calculation of the amount of glare sources is done for all output versions.
314	+	*/
315	+	/* evalglare.c, v2.00 2016/11/15  add of a second fast calculation mode for annual calculations, activted by -2. Output: dgp,ugr
316	+	*/
317	+	/* evalglare.c, v2.01 2016/11/16  change of -2 option (now -2 dir_illum). External provision of the direct illuminance necessary, since the sun interpolation of daysim is causing problems in calculation of the background luminance.
318	+	*/
319	+	/* evalglare.c, v2.02 2017/02/28  change of warning message, when invalid exposure setting is found. Reason: tab removal is not in all cases the right measure - it depends which tool caused the invalid exposure entry   */
320	+
321	+	/* evalglare.c, v2.03 2017/08/04  ad of -O option - disk replacement by providing luminance, not documented
322	+	*/
323	+
324	+	/* evalglare.c, v2.04 2017/08/04  adding -q option: use of Ev/pi as background luminance. not documented. no combination with -n option!!!
325	+	*/
326	+
327	+	/* evalglare.c, v2.05 2018/08/28  change of the -q option for the choice of the background luminance calculation mode: 0: CIE method (Ev-Edir)/pi, 1: mathematical correct average background luminance, 2: Ev/pi.
328	+	change of default options:
329	+	- cosinus weighted calculation of the background luminance (according to CIE) is now default.
330	+	- absolute threshold for the glare source detection is now default (2000cd/m2), based on study of C. Pierson
331	+	*/
332	+
333	+	/* evalglare.c, v2.06 2018/08/29
334	+	change of default value of multiplier b to 5.0, if task options (-t or -T ) are activated AND -b NOT used. To be downward compatible when using the task method.
335	+	*/
336	+
337	+	/* evalglare.c, v2.07 2018/11/17
338	+	bugfix: correction of error in the equations of PGSV_con and PGSV_sat
339	+	all three PGSV equations are calculated now
340	+	illuminance from the masking area (E_v_mask) is also printed
341	+	bugfix: in VCPs error fuction equation, value of 6.347 replaced by 6.374
342	+	*/
343	+
344	+	/* evalglare.c, v2.08 2018/11/27
345	+	bugfix: checkroutine for same image size for the masking corrected
346	+	*/
347	+
348	+
349	+
350	+	#define EVALGLARE
351		#define PROGNAME "evalglare"
352	<	#define VERSION "1.17 release 15.07.2015 by EPFL, J.Wienold"
352	>	#define VERSION "2.08 release 27.11.2018 by EPFL, J.Wienold"
353		#define RELEASENAME PROGNAME " " VERSION
354
355	<	#include "rtio.h"
288	<	#include "platform.h"
355	>
356		#include "pictool.h"
357	+	#include "rtio.h"
358		#include <math.h>
359		#include <string.h>
360	+	#include "platform.h"
361	+	#include "muc_randvar.h"
362	+
363		char *progname;
364
365		/* subroutine to add a pixel to a glare source */
#	Line 564 | Line 635 | double get_task_lum(pict * p, int x, int y, float r, i
635		return av_lum;
636		}
637
567	–	/* subroutine for calculation of band luminance */
568	–	double get_band_lum(pict * p, float r, int task_color)
569	–	{
570	–	int x_min, x_max, y_min, y_max, ix, iy, y_mid;
571	–	double r_actual, av_lum, omega_sum, act_lum;
638
639
574	–	x_max = pict_get_xsize(p) - 1;
575	–	y_max = pict_get_ysize(p) - 1;
576	–	x_min = 0;
577	–	y_min = 0;
578	–	y_mid = rint(y_max/2);
640
641
581	–
582	–	av_lum = 0.0;
583	–	omega_sum = 0.0;
584	–
585	–	for (iy = y_min; iy <= y_max; iy++) {
586	–	for (ix = x_min; ix <= x_max; ix++) {
587	–
588	–	/*                      if (DOT(pict_get_cached_dir(p,ix,iy),p->view.vdir) < 0.0)
589	–	continue;*/
590	–	r_actual =
591	–	acos(DOT(pict_get_cached_dir(p, ix, y_mid),
592	–	pict_get_cached_dir(p, ix, iy))) ;
593	–	act_lum = luminance(pict_get_color(p, ix, iy));
594	–
595	–	if ((r_actual <= r) || (iy == y_mid) ) {
596	–	act_lum = luminance(pict_get_color(p, ix, iy));
597	–	av_lum += pict_get_omega(p, ix, iy) * act_lum;
598	–	omega_sum += pict_get_omega(p, ix, iy);
599	–	if (task_color == 1) {
600	–	pict_get_color(p, ix, iy)[RED] = 0.0;
601	–	pict_get_color(p, ix, iy)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
602	–	pict_get_color(p, ix, iy)[BLU] = 0.0;
603	–	}
604	–	}
605	–	}
606	–	}
607	–
608	–	av_lum = av_lum / omega_sum;
609	–	/*    printf("average luminance of band %f \n",av_lum);*/
610	–	/*      printf("solid angle of band %f \n",omega_sum);*/
611	–	return av_lum;
612	–	}
613	–
614	–
615	–
616	–
617	–
642		/* subroutine for coloring the glare sources
643		red is used only for explicit call of the subroutine with acol=0  , e.g. for disability glare
644		-1: set to grey*/
#	Line 626 | Line 650 | int setglcolor(pict * p, int x, int y, int acol, int u
650		l=u_r+u_g+u_b ;
651
652		pcol[RED][0] = 1.0 / CIE_rf;
653	<	pcol[GRN][0] = 0.;
654	<	pcol[BLU][0] = 0.;
653	>	pcol[GRN][0] = 0.0 / CIE_gf;
654	>	pcol[BLU][0] = 0.0 / CIE_bf;
655
656	<	pcol[RED][1] = 0.;
656	>	pcol[RED][1] = 0.0 / CIE_rf;
657		pcol[GRN][1] = 1.0 / CIE_gf;
658	<	pcol[BLU][1] = 0.;
658	>	pcol[BLU][1] = 0.0 / CIE_bf;
659
660	<	pcol[RED][2] = 0.;
661	<	pcol[GRN][2] = 0.;
662	<	pcol[BLU][2] = 1 / CIE_bf;
660	>	pcol[RED][2] = 0.15 / CIE_rf;
661	>	pcol[GRN][2] = 0.15 / CIE_gf;
662	>	pcol[BLU][2] = 0.7 / CIE_bf;
663
664	<	pcol[RED][3] = 1.0 / (1.0 - CIE_bf);
665	<	pcol[GRN][3] = 1.0 / (1.0 - CIE_bf);
666	<	pcol[BLU][3] = 0.;
664	>	pcol[RED][3] = .5 / CIE_rf;
665	>	pcol[GRN][3] = .5 / CIE_gf;
666	>	pcol[BLU][3] = 0.0 / CIE_bf;
667
668	<	pcol[RED][4] = 1.0 / (1.0 - CIE_gf);
669	<	pcol[GRN][4] = 0.;
670	<	pcol[BLU][4] = 1.0 / (1.0 - CIE_gf);
668	>	pcol[RED][4] = .5 / CIE_rf;
669	>	pcol[GRN][4] = .0 / CIE_gf;
670	>	pcol[BLU][4] = .5 / CIE_bf ;
671
672	<	pcol[RED][5] = 0.;
673	<	pcol[GRN][5] = 1.0 / (1.0 - CIE_rf);
674	<	pcol[BLU][5] = 1.0 / (1.0 - CIE_rf);
672	>	pcol[RED][5] = .0 / CIE_rf;
673	>	pcol[GRN][5] = .5 / CIE_gf;
674	>	pcol[BLU][5] = .5 / CIE_bf;
675
676	<	pcol[RED][6] = 0.;
677	<	pcol[GRN][6] = 1.0 / (1.0 - CIE_rf);
678	<	pcol[BLU][6] = 1.0 / (1.0 - CIE_rf);
676	>	pcol[RED][6] = 0.333 / CIE_rf;
677	>	pcol[GRN][6] = 0.333 / CIE_gf;
678	>	pcol[BLU][6] = 0.333 / CIE_bf;
679
680
681		pcol[RED][999] = 1.0 / WHTEFFICACY;
#	Line 662 | Line 686 | int setglcolor(pict * p, int x, int y, int acol, int u
686		pcol[RED][998] = u_r /(l* CIE_rf) ;
687		pcol[GRN][998] = u_g /(l* CIE_gf);
688		pcol[BLU][998] = u_b /(l* CIE_bf);
689	<	/*        printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l);*/
689	>	/*        printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l,WHTEFFICACY); */
690		icol = acol ;
691		if ( acol == -1) {icol=999;
692		}else{if (acol>0){icol = acol % 5 +1;
#	Line 675 | Line 699 | int setglcolor(pict * p, int x, int y, int acol, int u
699		pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY;
700		pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY;
701		pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY;
702	<
702	>	/*        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))); */
703		return icol;
704		}
705
#	Line 727 | Line 751 | void add_secondary_gs(pict * p, int x, int y, float r,
751
752		/* color pixel of gs */
753
754	<	icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b);
754	>	/*              icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b); */
755
756		}
757		}
#	Line 767 | Line 791 | void split_pixel_from_gs(pict * p, int x, int y, int n
791
792		/* color pixel of new gs */
793
794	<	icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b);
794	>	/*      icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b); */
795	>
796
772	–
797		}
798
799		/* subroutine for the calculation of the position index */
#	Line 800 | Line 824 | float get_posindex(pict * p, float x, float y, int pos
824		tau = tau * deg;
825		sigma = sigma * deg;
826
827	+	if (postype == 1) {
828	+	/* KIM  model */
829	+	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));
830	+	}else{
831	+	/* Guth model, equation from IES lighting handbook */
832		posindex =
833		exp((35.2 - 0.31889 * tau -
834		1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 +
#	Line 821 | Line 850 | float get_posindex(pict * p, float x, float y, int pos
850
851		posindex = 1 + fact * r;
852		}
853	<	if (posindex > 16)
853	>	if (posindex > 16)
854		posindex = 16;
855	+	}
856
857		return posindex;
858
#	Line 1035 | Line 1065 | return;
1065
1066		}
1067
1068	<	/* subroutine for the calculation of the daylight glare index */
1068	>	/* subroutine for the calculation of the daylight glare index DGI*/
1069
1070
1071		float get_dgi(pict * p, float lum_backg, int igs, int posindex_2)
#	Line 1062 | Line 1092 | float get_dgi(pict * p, float lum_backg, int igs, int
1092		return dgi;
1093
1094		}
1095	+	/* subroutine for the calculation of the modified daylight glare index DGI_mod (Fisekis2003)
1096	+	several glare sources are taken into account and summed up, original equation has no summary
1097	+	*/
1098
1099	+
1100	+	float get_dgi_mod(pict * p, float lum_a, int igs, int posindex_2)
1101	+	{
1102	+	float dgi_mod, sum_glare, omega_s;
1103	+	int i;
1104	+
1105	+
1106	+	sum_glare = 0;
1107	+	omega_s = 0;
1108	+
1109	+	for (i = 0; i <= igs; i++) {
1110	+
1111	+	if (pict_get_npix(p, i) > 0) {
1112	+	omega_s =
1113	+	pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1114	+	get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1115	+	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));
1116	+	/*    printf("i,sum_glare %i %f\n",i,sum_glare); */
1117	+	}
1118	+	}
1119	+	dgi_mod = 10 * log10(sum_glare);
1120	+
1121	+	return dgi_mod;
1122	+
1123	+	}
1124	+
1125		/* subroutine for the calculation of the daylight glare probability */
1126		double
1127		get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3,
#	Line 1103 | Line 1162 | get_dgp(pict * p, double E_v, int igs, double a1, doub
1162		}
1163
1164		/* subroutine for the calculation of the visual comfort probability */
1165	<	float get_vcp(pict * p, double lum_a, int igs, int posindex_2)
1165	>	float get_dgr(pict * p, double lum_a, int igs, int posindex_2)
1166		{
1167	<	float vcp;
1168	<	double sum_glare, dgr;
1167	>	float dgr;
1168	>	double sum_glare;
1169		int i, i_glare;
1170
1171
#	Line 1132 | Line 1191 | float get_vcp(pict * p, double lum_a, int igs, int pos
1191		}
1192		dgr = pow(sum_glare, pow(i_glare, -0.0914));
1193
1194	<	vcp = 50 * erf((6.347 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1194	>
1195	>
1196	>	return dgr;
1197	>
1198	>	}
1199	>
1200	>	float get_vcp(float dgr )
1201	>	{
1202	>	float vcp;
1203	>
1204	>	vcp = 50 * erf((6.374 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1205		if (dgr > 750) {
1206		vcp = 0;
1207		}
1208		if (dgr < 20) {
1209		vcp = 100;
1210		}
1142	–
1143	–
1211		return vcp;
1212
1213		}
1214
1215	+
1216		/* subroutine for the calculation of the unified glare rating */
1217		float get_ugr(pict * p, double lum_backg, int igs, int posindex_2)
1218		{
#	Line 1168 | Line 1236 | float get_ugr(pict * p, double lum_backg, int igs, int
1236		}
1237		}
1238		ugr = 8 * log10(0.25 / lum_backg * sum_glare);
1239	<
1239	>	if (sum_glare==0) {
1240	>	ugr=0.0;
1241	>	}
1242	>	if (lum_backg<=0) {
1243	>	ugr=-99.0;
1244	>	}
1245	>
1246		return ugr;
1247
1248		}
1249	+	/* subroutine for the calculation of the experimental unified glare rating (Fisekis 2003)*/
1250	+	float get_ugr_exp(pict * p, double lum_backg, double lum_a, int igs, int posindex_2)
1251	+	{
1252	+	float ugr_exp;
1253	+	double sum_glare;
1254	+	int i, i_glare;
1255
1256
1257	+	sum_glare = 0;
1258	+	i_glare = 0;
1259	+	for (i = 0; i <= igs; i++) {
1260	+
1261	+	if (pict_get_npix(p, i) > 0) {
1262	+	i_glare = i_glare + 1;
1263	+	sum_glare +=
1264	+	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);
1265	+
1266	+	}
1267	+	}
1268	+	ugr_exp =8 * log10(lum_a) + 8 * log10(sum_glare/lum_backg);
1269	+
1270	+	return ugr_exp;
1271	+
1272	+	}
1273	+	/* subroutine for the calculation of the unified glare probability (Hirning 2014)*/
1274	+	float get_ugp(pict * p, double lum_backg, int igs, int posindex_2)
1275	+	{
1276	+	float ugp;
1277	+	double sum_glare;
1278	+	int i, i_glare;
1279	+
1280	+
1281	+	sum_glare = 0;
1282	+	i_glare = 0;
1283	+	for (i = 0; i <= igs; i++) {
1284	+
1285	+	if (pict_get_npix(p, i) > 0) {
1286	+	i_glare = i_glare + 1;
1287	+	sum_glare +=
1288	+	pow(pict_get_av_lum(p, i) /
1289	+	get_posindex(p, pict_get_av_posx(p, i),
1290	+	pict_get_av_posy(p, i), posindex_2),
1291	+	2) * pict_get_av_omega(p, i);
1292	+
1293	+	}
1294	+	}
1295	+	ugp = 0.26 * log10(0.25 / lum_backg * sum_glare);
1296	+
1297	+	return ugp;
1298	+
1299	+	}
1300	+
1301	+
1302		/* subroutine for the calculation of the disability glare according to Poynter */
1303		float get_disability(pict * p, double lum_backg, int igs)
1304		{
#	Line 1221 | Line 1346 | float get_disability(pict * p, double lum_backg, int i
1346
1347		/* subroutine for the calculation of the cie glare index */
1348
1349	<	float
1225	<	get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1349	>	float get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1350		{
1351		float cgi;
1352		double sum_glare;
#	Line 1246 | Line 1370 | get_cgi(pict * p, double E_v, double E_v_dir, int igs,
1370		cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare);
1371
1372		return cgi;
1373	+	}
1374
1375	+
1376	+
1377	+	/* subroutine for the calculation of the PGSV_con; is only applied, when masking is done, since it refers only to the window. Important: masking area must be the window! */
1378	+	float get_pgsv_con(double E_v, double E_mask,double omega_mask,double lum_mask_av, double Lavg)
1379	+	{
1380	+	float pgsv_con;
1381	+	double Lb;
1382	+
1383	+	/*        Lb = (E_v-E_mask)/3.14159265359;  */
1384	+	/*        Lb = (2*E_v-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask); */
1385	+	Lb = (2*3.14159265359*Lavg-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask);
1386	+
1387	+
1388	+	pgsv_con =3.2*log10(lum_mask_av)-0.64*log10(omega_mask)+(0.79*log10(omega_mask)-0.61)*log10(Lb)-8.2 ;
1389	+
1390	+
1391	+	return pgsv_con;
1392		}
1393
1394	+	/* 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! */
1395	+	float get_pgsv_sat(double E_v)
1396	+	{
1397	+	float pgsv_sat;
1398
1399	+	pgsv_sat =3.3-(0.57+3.3)/(1+pow(E_v/3.14159265359/1250,1.7));
1400	+
1401	+
1402	+	return pgsv_sat;
1403	+	}
1404	+
1405	+	/* 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! */
1406	+
1407	+	float get_pgsv(double E_v, double E_mask,double omega_mask,double lum_mask_av,double Ltask, double Lavg)
1408	+	{
1409	+	float pgsv;
1410	+	double Lb;
1411	+
1412	+	/*        Lb = (E_v-E_mask)/3.14159265359;  */
1413	+	/*        Lb = (2*E_v-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask); */
1414	+	Lb = (2*3.14159265359*Lavg-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask);
1415	+
1416	+	if (Lb==0.0 ) {
1417	+	fprintf(stderr,  " warning: Background luminance is 0 or masking area = full image! pgsv cannot be calculated (set to -99)!!\n");
1418	+	pgsv=-99;
1419	+	}else{
1420	+	if ( (lum_mask_av/Lb) > (E_v/(3.14159265359*Ltask))) {
1421	+	pgsv=get_pgsv_con(E_v,E_mask,omega_mask,lum_mask_av, Lavg);
1422	+	}else{
1423	+	pgsv=get_pgsv_sat(E_v)  ;
1424	+	}}
1425	+	return pgsv;
1426	+
1427	+	}
1428	+
1429	+
1430	+
1431		#ifdef  EVALGLARE
1432
1433
#	Line 1261 | Line 1439 | int main(int argc, char **argv)
1439		{
1440		#define CLINEMAX 4095 /* memory allocated for command line string */
1441		pict *p = pict_create();
1442	<	int     skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, y1, fill, yfillmax, yfillmin,
1443	<	ext_vill, set_lum_max, set_lum_max2, task_color, i_splitstart,
1444	<	i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,
1445	<	igs, actual_igs, icol, xt, yt, change, before_igs, sgs, splithigh,uniform_gs,
1446	<	detail_out, posindex_picture, non_cos_lb, rx, ry, rmx, rmy,apply_disability,band_calc,band_color;
1447	<	double  lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma,
1448	<	E_vl_ext, lum_max, new_lum_max, r_center,
1449	<	search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,
1450	<	omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang,
1451	<	l_max, lum_backg, lum_backg_cos, omega_sources, lum_sources,
1452	<	lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum;
1453	<	float lum_task, lum_thres, dgi,  vcp, cgi, ugr, limit,
1442	>	pict *pm = pict_create();
1443	>	int     skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, x1,y1, fill, yfillmax, yfillmin,
1444	>	ext_vill, set_lum_max, set_lum_max2, img_corr,x_disk,y_disk,task_color, i_splitstart,zones,act_gsn,splitgs,
1445	>	i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,x_zone,y_zone, i_z1, i_z2, thres_activate,
1446	>	igs, actual_igs, lastpixelwas_gs, icol, xt, yt, change,checkpixels, before_igs, sgs, splithigh,uniform_gs,x_max, y_max,y_mid,
1447	>	detail_out, posindex_picture, non_cos_lb, rx, ry, rmx,rmy,apply_disability,band_calc,band_color,masking,i_mask,no_glaresources,force;
1448	>	double  LUM_replace,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,
1449	>	E_vl_ext, lum_max, new_lum_max, r_center, ugp, ugr_exp, dgi_mod,lum_a, E_v_mask,angle_disk,dist,n_corner_px,zero_corner_px,
1450	>	search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,r_actual,lum_actual,dir_ill,
1451	>	omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang, angle_z1, angle_z2,per_95_band,per_75_band,pos,
1452	>	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,
1453	>	lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum,
1454	>	lum_mask[1],lum_mask_av,lum_mask_std[1],lum_mask_median[1],omega_mask,bbox[2],
1455	>	lum_band[1],lum_band_av,lum_band_std[1],lum_band_median[1],omega_band,bbox_band[2],
1456	>	lum_z1[1],lum_z1_av,lum_z1_std[1],lum_z1_median[1],omega_z1,bbox_z1[2],
1457	>	lum_z2[1],lum_z2_av,lum_z2_std[1],lum_z2_median[1],omega_z2,bbox_z2[2],
1458	>	lum_pos[1],lum_nopos_median[1],lum_pos_median[1],lum_pos2_median[1],lum_pos_mean,lum_pos2_mean;
1459	>	float lum_task, lum_thres, dgi,  vcp, cgi, ugr, limit, dgr,pgsv ,pgsv_sat,pgsv_con,
1460		abs_max, Lveil;
1461	<	char file_out[500], file_out2[500], version[500];
1461	>	char maskfile[500],file_out[500], file_out2[500], version[500];
1462		char *cline;
1463		VIEW userview = STDVIEW;
1464		int gotuserview = 0;
1465	+	struct muc_rvar* s_mask;
1466	+	s_mask = muc_rvar_create();
1467	+	muc_rvar_set_dim(s_mask, 1);
1468	+	muc_rvar_clear(s_mask);
1469	+	struct muc_rvar* s_band;
1470	+	s_band = muc_rvar_create();
1471	+	muc_rvar_set_dim(s_band, 1);
1472	+	muc_rvar_clear(s_band);
1473	+	struct muc_rvar* s_z1;
1474	+	s_z1 = muc_rvar_create();
1475	+	muc_rvar_set_dim(s_z1, 1);
1476	+	muc_rvar_clear(s_z1);
1477
1478	+	struct muc_rvar* s_z2;
1479	+	s_z2 = muc_rvar_create();
1480	+	muc_rvar_set_dim(s_z2, 1);
1481	+	muc_rvar_clear(s_z2);
1482	+
1483	+	struct muc_rvar* s_noposweight;
1484	+	s_noposweight = muc_rvar_create();
1485	+	muc_rvar_set_dim(s_noposweight, 1);
1486	+	muc_rvar_clear(s_noposweight);
1487	+
1488	+	struct muc_rvar* s_posweight;
1489	+	s_posweight = muc_rvar_create();
1490	+	muc_rvar_set_dim(s_posweight, 1);
1491	+	muc_rvar_clear(s_posweight);
1492	+
1493	+	struct muc_rvar* s_posweight2;
1494	+	s_posweight2 = muc_rvar_create();
1495	+	muc_rvar_set_dim(s_posweight2, 1);
1496	+	muc_rvar_clear(s_posweight2);
1497	+
1498		/*set required user view parameters to invalid values*/
1499	<	uniform_gs = 0;
1499	>	dir_ill=0.0;
1500	>	delta_E=0.0;
1501	>	no_glaresources=0;
1502	>	n_corner_px=0;
1503	>	zero_corner_px=0;
1504	>	force=0;
1505	>	dist=0.0;
1506	>	u_r=0.33333;
1507	>	u_g=0.33333;
1508	>	u_b=0.33333;
1509	>	band_angle=0;
1510	>	angle_z1=0;
1511	>	angle_z2=0;
1512	>	x_zone=0;
1513	>	y_zone=0;
1514	>	per_75_z2=0;
1515	>	per_95_z2=0;
1516	>	lum_pos_mean=0;
1517	>	lum_pos2_mean=0;
1518	>	lum_band_av = 0.0;
1519	>	omega_band = 0.0;
1520	>	pgsv = 0.0 ;
1521	>	pgsv_con = 0.0 ;
1522	>	pgsv_sat = 0.0 ;
1523	>	E_v_mask = 0.0;
1524	>	lum_z1_av = 0.0;
1525	>	omega_z1 = 0.0;
1526	>	lum_z2_av = 0.0;
1527	>	omega_z2 = 0.0 ;
1528	>	i_z1 = 0;
1529	>	i_z2 = 0;
1530	>	zones = 0;
1531	>	lum_z2_av = 0.0;
1532	>	uniform_gs = 0;
1533		apply_disability = 0;
1534		disability_thresh = 0;
1535		Lveil_cie_sum=0.0;
#	Line 1300 | Line 1549 | int main(int argc, char **argv)
1549		omegat = 0.0;
1550		yt = 0;
1551		xt = 0;
1552	+	x_disk=0;
1553	+	y_disk=0;
1554	+	angle_disk=0;
1555		yfillmin = 0;
1556		yfillmax = 0;
1557		cut_view = 0;
#	Line 1308 | Line 1560 | int main(int argc, char **argv)
1560		omega_cos_contr = 0.0;
1561		lum_ideal = 0.0;
1562		max_angle = 0.2;
1563	<	lum_thres = 5.0;
1563	>	lum_thres = 2000.0;
1564	>	lum_task = 0.0;
1565		task_lum = 0;
1566		sgs = 0;
1567		splithigh = 1;
#	Line 1326 | Line 1579 | int main(int argc, char **argv)
1579		c1 = 5.87e-05;
1580		c2 = 0.092;
1581		c3 = 0.159;
1582	<	non_cos_lb = 1;
1582	>	non_cos_lb = 0;
1583		posindex_2 = 0;
1584		task_color = 0;
1585		limit = 50000.0;
1586		set_lum_max = 0;
1587		set_lum_max2 = 0;
1588	+	img_corr=0;
1589		abs_max = 0;
1590		progname = argv[0];
1591		E_v_contr = 0.0;
1592	<	strcpy(version, "1.15 release 05.02.2015 by J.Wienold");
1592	>	strcpy(version, "1.19 release 09.12.2015 by J.Wienold");
1593		low_light_corr=1.0;
1594		output = 0;
1595		calc_vill = 0;
1596		band_avlum = -99;
1597		band_calc = 0;
1598	+	masking = 0;
1599	+	lum_mask[0]=0.0;
1600	+	lum_mask_av=0.0;
1601	+	omega_mask=0.0;
1602	+	i_mask=0;
1603	+	actual_igs=0;
1604	+	LUM_replace=0;
1605	+	thres_activate=0;
1606		/* command line for output picture*/
1607
1608		cline = (char *) malloc(CLINEMAX+1);
#	Line 1382 | Line 1644 | int main(int argc, char **argv)
1644		printf("age factor not supported any more \n");
1645		exit(1);
1646		break;
1647	+	case 'A':
1648	+	masking = 1;
1649	+	detail_out = 1;
1650	+	strcpy(maskfile, argv[++i]);
1651	+	pict_read(pm,maskfile );
1652	+
1653	+	break;
1654		case 'b':
1655		lum_thres = atof(argv[++i]);
1656	+	thres_activate = 1;
1657		break;
1658		case 'c':
1659		checkfile = 1;
#	Line 1401 | Line 1671 | int main(int argc, char **argv)
1671		case 's':
1672		sgs = 1;
1673		break;
1674	+	case 'f':
1675	+	force = 1;
1676	+	break;
1677		case 'k':
1678		apply_disability = 1;
1679		disability_thresh = atof(argv[++i]);
#	Line 1409 | Line 1682 | int main(int argc, char **argv)
1682		case 'p':
1683		posindex_picture = 1;
1684		break;
1685	+	case 'P':
1686	+	posindex_2 = atoi(argv[++i]);
1687	+	break;
1688	+
1689
1690		case 'y':
1691		splithigh = 1;
#	Line 1439 | Line 1716 | int main(int argc, char **argv)
1716		detail_out2 = 1;
1717		break;
1718		case 'm':
1719	+	img_corr = 1;
1720		set_lum_max = 1;
1721		lum_max = atof(argv[++i]);
1722		new_lum_max = atof(argv[++i]);
#	Line 1447 | Line 1725 | int main(int argc, char **argv)
1725		break;
1726
1727		case 'M':
1728	+	img_corr = 1;
1729		set_lum_max2 = 1;
1730		lum_max = atof(argv[++i]);
1731		E_vl_ext = atof(argv[++i]);
1732		strcpy(file_out2, argv[++i]);
1733		/*                      printf("max lum set to %f\n",new_lum_max);*/
1734		break;
1735	<	case 'n':
1735	>	case 'N':
1736	>	img_corr = 1;
1737	>	set_lum_max2 = 2;
1738	>	x_disk = atoi(argv[++i]);
1739	>	y_disk = atoi(argv[++i]);
1740	>	angle_disk = atof(argv[++i]);
1741	>	E_vl_ext = atof(argv[++i]);
1742	>	strcpy(file_out2, argv[++i]);
1743	>	/*                      printf("max lum set to %f\n",new_lum_max);*/
1744	>	break;
1745	>	case 'O':
1746	>	img_corr = 1;
1747	>	set_lum_max2 = 3;
1748	>	x_disk = atoi(argv[++i]);
1749	>	y_disk = atoi(argv[++i]);
1750	>	angle_disk = atof(argv[++i]);
1751	>	LUM_replace = atof(argv[++i]);
1752	>	strcpy(file_out2, argv[++i]);
1753	>	/*                      printf("max lum set to %f\n",new_lum_max);*/
1754	>	break;
1755	>
1756	>
1757	>	/* deactivated          case 'n':
1758		non_cos_lb = 0;
1759		break;
1760	+	*/
1761	+	case 'q':
1762	+	non_cos_lb = atoi(argv[++i]);
1763	+	break;
1764
1765		case 't':
1766		task_lum = 1;
#	Line 1472 | Line 1777 | int main(int argc, char **argv)
1777		omegat = atof(argv[++i]);
1778		task_color = 1;
1779		break;
1780	+	case 'l':
1781	+	zones = 1;
1782	+	x_zone = atoi(argv[++i]);
1783	+	y_zone = atoi(argv[++i]);
1784	+	angle_z1 = atof(argv[++i]);
1785	+	angle_z2 = -1;
1786	+	break;
1787	+	case 'L':
1788	+	zones = 2;
1789	+	x_zone = atoi(argv[++i]);
1790	+	y_zone = atoi(argv[++i]);
1791	+	angle_z1 = atof(argv[++i]);
1792	+	angle_z2 = atof(argv[++i]);
1793	+	break;
1794	+
1795	+
1796		case 'B':
1797		band_calc = 1;
1798		band_angle = atof(argv[++i]);
#	Line 1508 | Line 1829 | int main(int argc, char **argv)
1829		case '1':
1830		output = 1;
1831		break;
1832	+	case '2':
1833	+	output = 2;
1834	+	dir_ill = atof(argv[++i]);
1835	+	break;
1836
1837		case 'v':
1838		if (argv[i][2] == '\0') {
#	Line 1537 | Line 1862 | int main(int argc, char **argv)
1862		}
1863		}
1864
1865	+	/* set multiplier for task method to 5, if not specified */
1866	+
1867	+	if ( task_lum == 1 && thres_activate == 0){
1868	+	lum_thres = 5.0;
1869	+	}
1870		/*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/
1871
1872	<	if (output == 1 && ext_vill == 1) {
1872	>	if (output == 1 && ext_vill == 1 ) {
1873		calcfast=1;
1874		}
1875	+
1876	+	if (output == 2 && ext_vill == 1 ) {
1877	+	calcfast=2;
1878	+	}
1879	+
1880	+	/*masking and zoning cannot be applied at the same time*/
1881
1882	+	if (masking ==1 && zones >0) {
1883	+	fprintf(stderr,  " masking and zoning cannot be activated at the same time!\n");
1884	+	exit(1);
1885	+	}
1886	+
1887		/* read picture file */
1888		if (i == argc) {
1889		SET_FILE_BINARY(stdin);
#	Line 1576 | Line 1917 | if (output == 1 && ext_vill == 1) {
1917		return EXIT_FAILURE;
1918		}
1919
1920	+
1921	+
1922		/*                fprintf(stderr,"Picture read!");*/
1923
1924		/* command line for output picture*/
1925
1926		pict_set_comment(p, cline);
1927
1928	+	/* several checks */
1929		if (p->view.type == VT_PAR) {
1930	<	fprintf(stderr, "wrong view type! must not be parallel ! \n");
1930	>	fprintf(stderr, "error: wrong view type! must not be parallel ! \n");
1931		exit(1);
1932		}
1933
1934	+	if (p->view.type == VT_PER) {
1935	+	fprintf(stderr, "warning: image has perspective view type specified in header ! \n");
1936	+	}
1937
1938	+	if (masking == 1) {
1939	+
1940	+	if (pict_get_xsize(p)!=pict_get_xsize(pm) || pict_get_ysize(p)!=pict_get_ysize(pm)) {
1941	+	fprintf(stderr, "error: masking image has other resolution than main image ! \n");
1942	+	fprintf(stderr, "size must be identical \n");
1943	+	printf("resolution main image : %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
1944	+	printf("resolution masking image : %dx%d\n",pict_get_xsize(pm),pict_get_ysize(pm));
1945	+	exit(1);
1946	+
1947	+	}
1948	+
1949	+	}
1950	+	/*      printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
1951	+
1952		/* Check size of search radius */
1953		rmx = (pict_get_xsize(p) / 2);
1954		rmy = (pict_get_ysize(p) / 2);
#	Line 1632 | Line 1993 | if (output == 1 && ext_vill == 1) {
1993		avlum = 0.0;
1994		pict_new_gli(p);
1995		igs = 0;
1996	+	pict_get_z1_gsn(p,igs) = 0;
1997	+	pict_get_z2_gsn(p,igs) = 0;
1998
1999	+
2000		/* cut out GUTH field of view and exit without glare evaluation */
2001		if (cut_view==2) {
2002		if (cut_view_type==1) {
#	Line 1699 | Line 2063 | if (cut_view==2) {
2063		if (calcfast == 0) {
2064		for (x = 0; x < pict_get_xsize(p); x++)
2065		for (y = 0; y < pict_get_ysize(p); y++) {
2066	+	lum = luminance(pict_get_color(p, x, y));
2067	+	dist=sqrt((x-rmx)*(x-rmx)+(y-rmy)*(y-rmy));
2068	+	if (dist>((rmx+rmy)/2)) {
2069	+	n_corner_px=n_corner_px+1;
2070	+	if (lum<7.0) {
2071	+	zero_corner_px=zero_corner_px+1;
2072	+	}
2073	+	}
2074	+
2075	+
2076		if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2077		if (pict_is_validpixel(p, x, y)) {
1704	–	lum = luminance(pict_get_color(p, x, y));
2078		if (fill == 1 && y >= yfillmax) {
2079		y1 = y - 1;
2080		lum = luminance(pict_get_color(p, x, y1));
#	Line 1714 | Line 2087 | if (cut_view==2) {
2087		abs_max = lum;
2088		}
2089		/* set luminance restriction, if -m is set */
2090	<	if (set_lum_max == 1 && lum >= lum_max) {
2091	<	pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
2092	<	pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
2093	<	pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
2094	<	/*                                    printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
2095	<	lum = luminance(pict_get_color(p, x, y));
2090	>	if (img_corr == 1 ) {
2091	>	if (set_lum_max == 1 && lum >= lum_max) {
2092	>	pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
2093	>	pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
2094	>	pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
2095	>	/*                                                      printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
2096	>	lum = luminance(pict_get_color(p, x, y));
2097	>	}
2098	>	if (set_lum_max2 == 1 && lum >= lum_max) {
2099
2100	<	}
2101	<	if (set_lum_max2 == 1 && lum >= lum_max) {
2100	>	E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2101	>	omega_cos_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2102	>	}
2103	>	if (set_lum_max2 == 2 ) {
2104	>	r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2105	>	if (x_disk == x && y_disk==y ) r_actual=0.0;
2106
2107	<	E_v_contr +=
2108	<	DOT(p->view.vdir,
2109	<	pict_get_cached_dir(p, x,
2110	<	y)) * pict_get_omega(p,
2111	<	x,
2112	<	y)
2113	<	* lum;
2114	<	omega_cos_contr +=
2115	<	DOT(p->view.vdir,
2116	<	pict_get_cached_dir(p, x,
1737	<	y)) * pict_get_omega(p,
1738	<	x,
1739	<	y)
1740	<	* 1;
2107	>	act_lum = luminance(pict_get_color(p, x, y));
2108	>
2109	>	if (r_actual <= angle_disk) {
2110	>	E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2111	>	omega_cos_contr += DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2112	>	}
2113	>
2114	>
2115	>
2116	>	}
2117		}
2118	+	om = pict_get_omega(p, x, y);
2119	+	sang += om;
2120	+	E_v += DOT(p->view.vdir, pict_get_cached_dir(p, x,y)) * om *lum;
2121	+	avlum += om * lum;
2122	+	pos=get_posindex(p, x, y,posindex_2);
2123
2124	+	lum_pos[0]=lum;
2125	+	muc_rvar_add_sample(s_noposweight,lum_pos);
2126	+	lum_pos[0]=lum/pos;
2127	+	lum_pos_mean+=lum_pos[0]*om;
2128	+	muc_rvar_add_sample(s_posweight, lum_pos);
2129	+	lum_pos[0]=lum_pos[0]/pos;
2130	+	lum_pos2_mean+=lum_pos[0]*om;
2131	+	muc_rvar_add_sample(s_posweight2,lum_pos);
2132
2133	<	sang += pict_get_omega(p, x, y);
2134	<	E_v +=
1746	<	DOT(p->view.vdir,
1747	<	pict_get_cached_dir(p, x,
1748	<	y)) * pict_get_omega(p, x,
1749	<	y) *
1750	<	lum;
1751	<	avlum +=
1752	<	pict_get_omega(p, x,
1753	<	y) * luminance(pict_get_color(p, x,
1754	<	y));
2133	>
2134	>
2135		} else {
2136		pict_get_color(p, x, y)[RED] = 0.0;
2137		pict_get_color(p, x, y)[GRN] = 0.0;
2138		pict_get_color(p, x, y)[BLU] = 0.0;
2139
2140		}
2141	<	}
2141	>	}else {
2142	>	pict_get_color(p, x, y)[RED] = 0.0;
2143	>	pict_get_color(p, x, y)[GRN] = 0.0;
2144	>	pict_get_color(p, x, y)[BLU] = 0.0;
2145	>
2146	>	}
2147		}
2148
2149	<	if (set_lum_max2 == 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) {
2149	>	/* check if image has black corners AND a perspective view */
2150
2151	+	if (zero_corner_px/n_corner_px > 0.70 && (p->view.type == VT_PER) ) {
2152	+	printf (" corner pixels are to  %f %% black! \n",zero_corner_px/n_corner_px*100);
2153	+	printf("error! The image has black corners AND header contains a perspective view type definition !!!\n") ;
2154	+
2155	+	if (force==0) {
2156	+	printf("stopping...!!!!\n") ;
2157	+
2158	+	exit(1);
2159	+	}
2160	+	}
2161	+
2162	+
2163	+
2164	+
2165	+	lum_pos_mean= lum_pos_mean/sang;
2166	+	lum_pos2_mean= lum_pos2_mean/sang;
2167	+
2168	+	if ((set_lum_max2 >= 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0 ) || set_lum_max2==3) {
2169	+
2170	+	if (set_lum_max2<3){
2171		lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
2172	+	if (set_lum_max2 == 2 && lum_ideal >= 2e9) {
2173	+	printf("warning! luminance of replacement pixels would be larger than 2e9 cd/m2. Value set to 2e9cd/m2!\n") ;
2174	+	}
2175	+
2176		if (lum_ideal > 0 && lum_ideal < setvalue) {
2177		setvalue = lum_ideal;
2178		}
2179	<	printf
1771	<	("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f  %f %f %f %f\n",
2179	>	printf("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f  %f %f %f %f\n",
2180		lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
2181	+	}else{setvalue=LUM_replace;
2182	+	}
2183
2184	<
2184	>
2185		for (x = 0; x < pict_get_xsize(p); x++)
2186		for (y = 0; y < pict_get_ysize(p); y++) {
2187		if (pict_get_hangle
2188		(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2189		if (pict_is_validpixel(p, x, y)) {
2190		lum = luminance(pict_get_color(p, x, y));
2191	+
2192	+
2193		if (set_lum_max2 == 1 && lum >= lum_max) {
2194
2195		pict_get_color(p, x, y)[RED] =
#	Line 1787 | Line 2199 | if (cut_view==2) {
2199		pict_get_color(p, x, y)[BLU] =
2200		setvalue / 179.0;
2201
2202	+	}else{ if(set_lum_max2 >1 ) {
2203	+	r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2204	+	if (x_disk == x && y_disk==y ) r_actual=0.0;
2205	+
2206	+	if (r_actual <= angle_disk) {
2207	+
2208	+	pict_get_color(p, x, y)[RED] =
2209	+	setvalue / 179.0;
2210	+	pict_get_color(p, x, y)[GRN] =
2211	+	setvalue / 179.0;
2212	+	pict_get_color(p, x, y)[BLU] =
2213	+	setvalue / 179.0;
2214	+
2215	+	}
2216		}
2217	+	}
2218		}
2219		}
2220		}
2221	+
2222
2223		pict_write(p, file_out2);
2224	<
2224	>	exit(1);
2225		}
2226		if (set_lum_max == 1) {
2227		pict_write(p, file_out2);
#	Line 1853 | Line 2281 | if (cut_view==1) {
2281		lum_source = lum_thres;
2282		}
2283
2284	<	/*     printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2284	>	/*     printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2285
2286		/* first glare source scan: find primary glare sources */
2287	<	for (x = 0; x < pict_get_xsize(p); x++)
2287	>	for (x = 0; x < pict_get_xsize(p); x++) {
2288	>	/*                lastpixelwas_gs=0; */
2289		for (y = 0; y < pict_get_ysize(p); y++) {
2290		if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2291		if (pict_is_validpixel(p, x, y)) {
#	Line 1865 | Line 2294 | if (cut_view==1) {
2294		if (act_lum >lum_total_max) {
2295		lum_total_max=act_lum;
2296		}
2297	<
2298	<	actual_igs =
2299	<	find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2297	>	/* speed improvement first scan: when last pixel was glare source, then use this glare source number instead of searching.
2298	>	has been deactivated with v1.25
2299	>
2300	>	if (lastpixelwas_gs==0 || search_pix <= 1.0 ) { */
2301	>	actual_igs = find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2302	>	/*  }*/
2303		if (actual_igs == 0) {
2304		igs = igs + 1;
2305		pict_new_gli(p);
2306		pict_get_lum_min(p, igs) = HUGE_VAL;
2307		pict_get_Eglare(p,igs) = 0.0;
2308		pict_get_Dglare(p,igs) = 0.0;
2309	+	pict_get_z1_gsn(p,igs) = 0;
2310	+	pict_get_z2_gsn(p,igs) = 0;
2311		actual_igs = igs;
2312	+
2313		}
2314	<	icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2314	>	/* no coloring any more here            icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2315		pict_get_gsn(p, x, y) = actual_igs;
2316		pict_get_pgs(p, x, y) = 1;
2317		add_pixel_to_gs(p, x, y, actual_igs);
2318	+	/*                                                lastpixelwas_gs=actual_igs; */
2319
2320		} else {
2321		pict_get_pgs(p, x, y) = 0;
2322	+	/*                                               lastpixelwas_gs=0;*/
2323		}
2324		}
2325		}
2326	<	}
2326	>	}
2327	>	}
2328		/*      pict_write(p,"firstscan.pic");   */
2329
2330	<	if (calcfast == 1 ) {
2330	>
2331	>
2332	>
2333	>	if (calcfast ==1 || search_pix <= 1.0 || calcfast == 2 ) {
2334		skip_second_scan=1;
2335		}
2336	+
2337
2338		/* second glare source scan: combine glare sources facing each other */
2339		change = 1;
2340	+	i = 0;
2341		while (change == 1 && skip_second_scan==0) {
2342		change = 0;
2343	<	for (x = 0; x < pict_get_xsize(p); x++)
2343	>	i = i+1;
2344	>	for (x = 0; x < pict_get_xsize(p); x++) {
2345		for (y = 0; y < pict_get_ysize(p); y++) {
1902	–	before_igs = pict_get_gsn(p, x, y);
2346		if (pict_get_hangle
2347		(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2348	<	if (pict_is_validpixel(p, x, y) && before_igs > 0) {
2348	>	checkpixels=1;
2349	>	before_igs = pict_get_gsn(p, x, y);
2350	>
2351	>	/* speed improvement: search for other nearby glare sources only, when at least one adjacend pixel has another glare source number  */
2352	>	if (x >1 && x<pict_get_xsize(p)-2 && y >1 && y<pict_get_ysize(p)-2 )   {
2353	>	x1=x-1;
2354	>	x2=x+1;
2355	>	y1=y-1;
2356	>	y2=y+1;
2357	>	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) ) {
2358	>	checkpixels = 0;
2359	>	actual_igs = before_igs;
2360	>	}  }
2361	>
2362	>
2363	>	if (pict_is_validpixel(p, x, y) && before_igs > 0 && checkpixels==1) {
2364		actual_igs =
2365		find_near_pgs(p, x, y, max_angle, before_igs,
2366		igs, 1);
#	Line 1911 | Line 2369 | if (calcfast == 1 ) {
2369		}
2370		if (before_igs > 0) {
2371		actual_igs = pict_get_gsn(p, x, y);
2372	<	icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2372	>	/*      icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2373		}
2374
2375		}
2376		}
2377		}
2378		/*      pict_write(p,"secondscan.pic");*/
2379	+	}}
2380
1922	–	}
1923	–
2381		/*smoothing: add secondary glare sources */
2382		i_max = igs;
2383		if (sgs == 1) {
#	Line 1977 | Line 2434 | if (calcfast == 1 ) {
2434		if (i_splitstart == (igs + 1)) {
2435		igs = igs + 1;
2436		pict_new_gli(p);
2437	+	pict_get_z1_gsn(p,igs) = 0;
2438	+	pict_get_z2_gsn(p,igs) = 0;
2439	+
2440		pict_get_Eglare(p,igs) = 0.0;
2441		pict_get_Dglare(p,igs) = 0.0;
2442		pict_get_lum_min(p, igs) =
#	Line 1990 | Line 2450 | if (calcfast == 1 ) {
2450		if (i_split == 0) {
2451		igs = igs + 1;
2452		pict_new_gli(p);
2453	+	pict_get_z1_gsn(p,igs) = 0;
2454	+	pict_get_z2_gsn(p,igs) = 0;
2455	+
2456		pict_get_Eglare(p,igs) = 0.0;
2457		pict_get_Dglare(p,igs) = 0.0;
2458		pict_get_lum_min(p, igs) =
#	Line 2026 | Line 2489 | if (calcfast == 1 ) {
2489		if (before_igs > 0) {
2490		actual_igs =
2491		pict_get_gsn(p, x, y);
2492	<	icol =
2492	>	/*                                                                                     icol =
2493		setglcolor(p, x, y,
2494	<	actual_igs, uniform_gs, u_r, u_g , u_b);
2494	>	actual_igs, uniform_gs, u_r, u_g , u_b);*/
2495		}
2496
2497		}
#	Line 2043 | Line 2506 | if (calcfast == 1 ) {
2506		}
2507		}
2508
2509	<	/* calculation of direct vertical illuminance for CGI and for disability glare*/
2509	>	/* calculation of direct vertical illuminance for CGI and for disability glare, coloring glare sources*/
2510
2511	<
2049	<	if (calcfast == 0) {
2511	>	if (calcfast == 0 || calcfast == 2) {
2512		for (x = 0; x < pict_get_xsize(p); x++)
2513		for (y = 0; y < pict_get_ysize(p); y++) {
2514		if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2515		if (pict_is_validpixel(p, x, y)) {
2516		if (pict_get_gsn(p, x, y) > 0) {
2517	<	pict_get_Eglare(p,pict_get_gsn(p, x, y)) +=
2518	<	DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2519	<	* pict_get_omega(p, x, y)
2520	<	* luminance(pict_get_color(p, x, y));
2521	<	E_v_dir +=
2060	<	DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2061	<	* pict_get_omega(p, x, y)
2062	<	* luminance(pict_get_color(p, x, y));
2517	>	actual_igs = pict_get_gsn(p, x, y);
2518	>	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));
2519	>	pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
2520	>	E_v_dir = E_v_dir + delta_E;
2521	>	setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2522		}
2523		}
2524		}
2525		}
2526		lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2527	<	lum_backg = lum_backg_cos;
2527	>
2528		}
2529	<	/* calc of band luminance if applied */
2529	>	/* calc of band luminance distribution if applied */
2530		if (band_calc == 1) {
2531	<	band_avlum=get_band_lum( p, band_angle, band_color);
2532	<	}
2531	>	x_max = pict_get_xsize(p) - 1;
2532	>	y_max = pict_get_ysize(p) - 1;
2533	>	y_mid = (int)(y_max/2);
2534	>	for (x = 0; x <= x_max; x++)
2535	>	for (y = 0; y <= y_max; y++) {
2536	>	if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2537	>	if (pict_is_validpixel(p, x, y)) {
2538
2539	+	r_actual = acos(DOT(pict_get_cached_dir(p, x, y_mid), pict_get_cached_dir(p, x, y))) ;
2540	+	act_lum = luminance(pict_get_color(p, x, y));
2541	+
2542	+	if ((r_actual <= band_angle) || (y == y_mid) ) {
2543	+	lum_band[0]=luminance(pict_get_color(p, x, y));
2544	+	muc_rvar_add_sample(s_band, lum_band);
2545	+	act_lum = luminance(pict_get_color(p, x, y));
2546	+	lum_band_av += pict_get_omega(p, x, y) * act_lum;
2547	+	omega_band += pict_get_omega(p, x, y);
2548	+	if (band_color == 1) {
2549	+	pict_get_color(p, x, y)[RED] = 0.0;
2550	+	pict_get_color(p, x, y)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
2551	+	pict_get_color(p, x, y)[BLU] = 0.0;
2552	+	}
2553	+	}
2554	+	}}}
2555	+	lum_band_av=lum_band_av/omega_band;
2556	+	muc_rvar_get_vx(s_band,lum_band_std);
2557	+	muc_rvar_get_percentile(s_band,lum_band_median,0.75);
2558	+	per_75_band=lum_band_median[0];
2559	+	muc_rvar_get_percentile(s_band,lum_band_median,0.95);
2560	+	per_95_band=lum_band_median[0];
2561	+	muc_rvar_get_median(s_band,lum_band_median);
2562	+	muc_rvar_get_bounding_box(s_band,bbox_band);
2563	+
2564	+	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] );
2565	+	/*      av_lum = av_lum / omega_sum; */
2566	+	/*    printf("average luminance of band %f \n",av_lum);*/
2567	+	/*      printf("solid angle of band %f \n",omega_sum);*/
2568	+	}
2569	+
2570		/*printf("total number of glare sources: %i\n",igs); */
2571		lum_sources = 0;
2572		omega_sources = 0;
2573	+	i=0;
2574		for (x = 0; x <= igs; x++) {
2575	+	if (pict_get_npix(p, x) > 0) {
2576		lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2577		omega_sources += pict_get_av_omega(p, x);
2578	+	i=i+1;
2579	+	}}
2580	+
2581	+	if (sang == omega_sources) {
2582	+	lum_backg =avlum;
2583	+	} else {
2584	+	lum_backg =(sang * avlum - lum_sources) / (sang - omega_sources);
2585		}
2586	<	if (non_cos_lb == 1) {
2587	<	lum_backg =
2588	<	(sang * avlum - lum_sources) / (sang - omega_sources);
2586	>
2587	>
2588	>	if (i == 0) {
2589	>	lum_sources =0.0;
2590	>	} else { lum_sources=lum_sources/omega_sources;
2591		}
2592	+
2593	+
2594	+
2595	+	if (non_cos_lb == 0) {
2596	+	lum_backg = lum_backg_cos;
2597	+	}
2598	+
2599	+	if (non_cos_lb == 2) {
2600	+	lum_backg = E_v / 3.1415927;
2601	+	}
2602	+
2603	+
2604	+	/* file writing NOT here
2605	+	if (checkfile == 1) {
2606	+	pict_write(p, file_out);
2607	+	}
2608	+	*/
2609	+
2610		/* print detailed output */
2611	+
2612	+
2613	+
2614	+	/* masking */
2615	+
2616	+	if ( masking == 1 ) {
2617	+
2618	+
2619	+	for (x = 0; x < pict_get_xsize(p); x++)
2620	+	for (y = 0; y < pict_get_ysize(p); y++) {
2621	+	if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2622	+	if (pict_is_validpixel(p, x, y)) {
2623	+	if (luminance(pict_get_color(pm, x, y))>0.1) {
2624	+	/*                                         printf ("hallo %i %i %f \n",x,y,luminance(pict_get_color(pm, x, y)));*/
2625	+	i_mask=i_mask+1;
2626	+	lum_mask[0]=luminance(pict_get_color(p, x, y));
2627	+	/* printf ("%f \n",lum_mask[0]);*/
2628	+	muc_rvar_add_sample(s_mask, lum_mask);
2629	+	omega_mask += pict_get_omega(p, x, y);
2630	+	lum_mask_av += pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2631	+	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));
2632	+
2633	+	pict_get_color(pm, x, y)[RED] = luminance(pict_get_color(p, x, y))/179.0;
2634	+	pict_get_color(pm, x, y)[GRN] = luminance(pict_get_color(p, x, y))/179.0;
2635	+	pict_get_color(pm, x, y)[BLU] = luminance(pict_get_color(p, x, y))/179.0;
2636	+
2637	+	} else {
2638	+	pict_get_color(p, x, y)[RED] = 0.0 ;
2639	+	pict_get_color(p, x, y)[GRN] = 0.0 ;
2640	+	pict_get_color(p, x, y)[BLU] = 0.0 ;
2641	+	}
2642	+	}
2643	+	}
2644	+	}
2645	+	/* Average luminance in masking area (weighted by solid angle) */
2646	+	lum_mask_av=lum_mask_av/omega_mask;
2647	+	muc_rvar_get_vx(s_mask,lum_mask_std);
2648	+	muc_rvar_get_percentile(s_mask,lum_mask_median,0.75);
2649	+	per_75_mask=lum_mask_median[0];
2650	+	muc_rvar_get_percentile(s_mask,lum_mask_median,0.95);
2651	+	per_95_mask=lum_mask_median[0];
2652	+	muc_rvar_get_median(s_mask,lum_mask_median);
2653	+	muc_rvar_get_bounding_box(s_mask,bbox);
2654	+	/* PSGV only why masking of window is applied! */
2655	+
2656	+
2657	+	if (task_lum == 0 || lum_task == 0.0 ) {
2658	+	fprintf(stderr, " warning: Task area not set or task luminance=0 ! pgsv cannot be calculated (set to -99)!!\n");
2659	+	pgsv = -99;
2660	+	} else {
2661	+	pgsv = get_pgsv(E_v, E_v_mask, omega_mask, lum_mask_av,lum_task,avlum);
2662	+	}
2663	+
2664	+	pgsv_con = get_pgsv_con(E_v, E_v_mask, omega_mask, lum_mask_av,avlum);
2665	+	pgsv_sat =get_pgsv_sat(E_v);
2666	+
2667		if (detail_out == 1) {
2668	+
2669	+	printf ("masking:no_pixels,omega,av_lum,median_lum,std_lum,perc_75,perc_95,lum_min,lum_max,pgsv_con,pgsv_sat,pgsv,Ev_mask: %i %f %f %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_con,pgsv_sat,pgsv,E_v_mask );
2670	+
2671	+	}
2672	+
2673	+	}
2674	+
2675	+	/* zones */
2676	+
2677	+	if ( zones > 0 ) {
2678	+
2679	+	for (x = 0; x < pict_get_xsize(p); x++)
2680	+	for (y = 0; y < pict_get_ysize(p); y++) {
2681	+	if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2682	+	if (pict_is_validpixel(p, x, y)) {
2683	+
2684	+
2685	+	r_actual =acos(DOT(pict_get_cached_dir(p, x, y), pict_get_cached_dir(p, x_zone, y_zone))) * 2;
2686	+	lum_actual = luminance(pict_get_color(p, x, y));
2687	+	act_gsn=pict_get_gsn(p, x, y);
2688	+	/*     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));*/
2689	+
2690	+	/*zone1 */
2691	+	if (r_actual <= angle_z1) {
2692	+	i_z1=i_z1+1;
2693	+	lum_z1[0]=lum_actual;
2694	+	muc_rvar_add_sample(s_z1, lum_z1);
2695	+	omega_z1 += pict_get_omega(p, x, y);
2696	+	lum_z1_av += pict_get_omega(p, x, y)* lum_actual;
2697	+	setglcolor(p,x,y,1,1 , 0.66, 0.01 ,0.33);
2698	+	/*check if separation gsn already exist */
2699	+
2700	+	if (act_gsn > 0){
2701	+	if (pict_get_z1_gsn(p,act_gsn) == 0) {
2702	+	pict_new_gli(p);
2703	+	igs = igs + 1;
2704	+	pict_get_z1_gsn(p,act_gsn) = igs*1.000;
2705	+	pict_get_z1_gsn(p,igs) = -1.0;
2706	+	pict_get_z2_gsn(p,igs) = -1.0;
2707	+	splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2708	+	/*                                                                                        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)); */
2709	+	}
2710	+	splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
2711	+	/*          printf("splitgs%i \n",splitgs);       */
2712	+	split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
2713	+	}
2714	+	}
2715	+	/*zone2 */
2716	+
2717	+	if (r_actual > angle_z1 && r_actual<= angle_z2 ) {
2718	+
2719	+	i_z2=i_z2+1;
2720	+	lum_z2[0]=lum_actual;
2721	+	muc_rvar_add_sample(s_z2, lum_z2);
2722	+	omega_z2 +=   pict_get_omega(p, x, y);
2723	+	lum_z2_av += pict_get_omega(p, x, y)* lum_actual;
2724	+	setglcolor(p,x,y,1,1 , 0.65, 0.33 ,0.02);
2725	+	/*                                                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));
2726	+	*/                                                if (act_gsn > 0){
2727	+	if (pict_get_z2_gsn(p,act_gsn) == 0) {
2728	+	pict_new_gli(p);
2729	+	igs = igs + 1;
2730	+	pict_get_z2_gsn(p,act_gsn) = igs*1.000;
2731	+	pict_get_z1_gsn(p,igs) = -2.0;
2732	+	pict_get_z2_gsn(p,igs) = -2.0;
2733	+	/*                                                                                         printf ("Glare source %i gets glare source %i in zone 2 \n",act_gsn,igs ); */
2734	+	}
2735	+	splitgs=(int)(pict_get_z2_gsn(p,act_gsn));
2736	+	/*                                              printf("splitgs %i \n",splitgs);*/
2737	+	split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
2738	+	}
2739	+	}
2740	+
2741	+	}
2742	+	} }
2743	+	/* Average luminance in zones (weighted by solid angle), calculate percentiles, median min and max in zones  */
2744	+	if (zones == 2 ) {
2745	+	lum_z2_av=lum_z2_av/omega_z2;
2746	+	muc_rvar_get_vx(s_z2,lum_z2_std);
2747	+	muc_rvar_get_percentile(s_z2,lum_z2_median,0.75);
2748	+	per_75_z2=lum_z2_median[0];
2749	+	muc_rvar_get_percentile(s_z2,lum_z2_median,0.95);
2750	+	per_95_z2=lum_z2_median[0];
2751	+	muc_rvar_get_median(s_z2,lum_z2_median);
2752	+	muc_rvar_get_bounding_box(s_z2,bbox_z2);
2753	+	}
2754	+	lum_z1_av=lum_z1_av/omega_z1;
2755	+	muc_rvar_get_vx(s_z1,lum_z1_std);
2756	+	muc_rvar_get_percentile(s_z1,lum_z1_median,0.75);
2757	+	per_75_z1=lum_z1_median[0];
2758	+	muc_rvar_get_percentile(s_z1,lum_z1_median,0.95);
2759	+	per_95_z1=lum_z1_median[0];
2760	+	muc_rvar_get_median(s_z1,lum_z1_median);
2761	+	muc_rvar_get_bounding_box(s_z1,bbox_z1);
2762	+	if (detail_out == 1) {
2763	+
2764	+	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] );
2765	+
2766	+	if (zones == 2 ) {
2767	+
2768	+	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] );
2769	+	} }
2770	+
2771	+	}
2772	+
2773		i = 0;
2774		for (x = 0; x <= igs; x++) {
2775		/* resorting glare source numbers */
#	Line 2092 | Line 2777 | if (calcfast == 1 ) {
2777		i = i + 1;
2778		}
2779		}
2780	+	no_glaresources=i;
2781
2782	+	/* glare sources */
2783	+	if (detail_out == 1) {
2784
2097	–
2785		printf
2786	<	("%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",
2786	>	("%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",
2787		i);
2788		if (i == 0) {
2789	<	printf("%i %f %f %f %f %.10f %f %f %f %f %f %f\n", i, 0.0, 0.0,
2790	<	0.0, 0.0, 0.0, 0.0, lum_backg, lum_task, E_v, E_v_dir,
2104	<	abs_max);
2789	>	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,
2790	>	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);
2791
2792		} else {
2793		i = 0;
#	Line 2121 | Line 2807 | if (calcfast == 1 ) {
2807		Lveil_cie =0 ;
2808		}
2809		Lveil_cie_sum =  Lveil_cie_sum + Lveil_cie;
2810	<	printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
2810	>	if (pict_get_z1_gsn(p,x)<0) {
2811	>	act_gsn=(int)(-pict_get_z1_gsn(p,x));
2812	>	}else{
2813	>	act_gsn=0;
2814	>	}
2815	>	printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %i \n",
2816		i, pict_get_npix(p, x), pict_get_av_posx(p, x),
2817		pict_get_ysize(p) - pict_get_av_posy(p, x),
2818		pict_get_av_lum(p, x), pict_get_av_omega(p, x),
#	Line 2129 | Line 2820 | if (calcfast == 1 ) {
2820		pict_get_av_posy(p, x),
2821		posindex_2), lum_backg, lum_task,
2822		E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
2823	<	,pict_get_cached_dir(p, x2,y2)[0],pict_get_cached_dir(p, x2,y2)[1],pict_get_cached_dir(p, x2,y2)[2],pict_get_Eglare(p,x),Lveil_cie,teta
2133	<
2134	<	);
2823	>	,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 );
2824		}
2825		}
2826		}
#	Line 2171 | Line 2860 | if (calcfast == 1 ) {
2860		}
2861
2862		if (calcfast == 0) {
2863	+	lum_a= E_v2/3.1415927;
2864		dgi = get_dgi(p, lum_backg, igs, posindex_2);
2865		ugr = get_ugr(p, lum_backg, igs, posindex_2);
2866	+	ugp = get_ugp (p,lum_backg, igs, posindex_2);
2867	+	ugr_exp = get_ugr_exp (p,lum_backg_cos,lum_a, igs, posindex_2);
2868	+	dgi_mod = get_dgi_mod(p, lum_a, igs, posindex_2);
2869		cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
2870	<	vcp = get_vcp(p, avlum, igs, posindex_2);
2870	>	dgr = get_dgr(p, avlum, igs, posindex_2);
2871	>	vcp = get_vcp(dgr);
2872		Lveil = get_disability(p, avlum, igs);
2873	+	if (no_glaresources == 0) {
2874	+	dgi = 0.0;
2875	+	ugr = 0.0;
2876	+	ugp = 0.0;
2877	+	ugr_exp = 0.0;
2878	+	dgi_mod = 0.0;
2879	+	cgi = 0.0;
2880	+	dgr = 0.0;
2881	+	vcp = 100.0;
2882	+	}
2883		}
2884		/* check dgp range */
2885		if (dgp <= 0.2) {
#	Line 2201 | Line 2905 | if (calcfast == 0) {
2905		dgp =low_light_corr*dgp;
2906		dgp =age_corr_factor*dgp;
2907		}
2908	+	muc_rvar_get_median(s_noposweight,lum_nopos_median);
2909	+	muc_rvar_get_median(s_posweight,lum_pos_median);
2910	+	muc_rvar_get_median(s_posweight2,lum_pos2_median);
2911	+
2912
2913	+
2914	+
2915		printf
2916	<	("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",
2916	>	("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",
2917		dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
2918	<	lum_sources / omega_sources, omega_sources, Lveil,Lveil_cie_sum,band_avlum);
2918	>	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]);
2919		} else {
2920		if (detail_out2 == 1) {
2921
#	Line 2236 | Line 2946 | if (calcfast == 0) {
2946		if (E_vl_ext < 1000) {
2947		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 ;}
2948		dgp =low_light_corr*dgp;
2949	<	dgp =age_corr_factor*dgp;
2950	<	printf("%f\n", dgp);
2949	>
2950	>	if (calcfast == 2) {
2951	>
2952	>	lum_backg_cos=(E_vl_ext-dir_ill)/3.1415927;
2953	>	ugr = get_ugr(p, lum_backg_cos, igs, posindex_2);
2954	>	printf("%f %f \n", dgp,ugr);
2955	>	}else{
2956	>	printf("%f\n", dgp);
2957	>	}
2958		}
2959
2960
#	Line 2268 | Line 2985 | has to be re-written from scratch....
2985
2986		/*output  */
2987		if (checkfile == 1) {
2988	+
2989	+
2990		pict_write(p, file_out);
2991		}
2992

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines