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

Comparing ray/src/util/evalglare.c (file contents):
Revision 2.2 by greg, Tue Aug 18 15:02:53 2015 UTC vs.
Revision 3.08 by greg, Wed Oct 1 17:48:15 2025 UTC

# Line 1 | Line 1
1   #ifndef lint
2 < static const char RCSid[] = "$Id$";
2 > static const char RCSid[] = "$Id$";
3   #endif
4 < /* EVALGLARE V1.17
5 < * Evalglare Software License, Version 1.0
4 > /* EVALGLARE V3.06
5 > * Evalglare Software License, Version 3.0x
6   *
7 < * Copyright (c) 1995 - 2015 Fraunhofer ISE, EPFL.
7 > * Copyright (c) 1995 - 2022 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 + /* evalglare.c, v2.09 2019/01/18
349 + calculate "illuminance-contribution of zones"
350 + -switch to turn off low-light correction: 4
351 +  */
352 +
353 + /* evalglare.c, v2.10 2019/07/30, 2020/06/25
354 + -add multiimage-mode for annual glare evaluation -Q
355 + -extra output for multiimage mode
356 + 2020/06/22 - added Ev contribution from each glare source as output to the multiimage mode
357 + -switch for correction modes (-C), value of 0 switches off all corrections
358 +  */
359 +
360 +
361 + /* evalglare.c, v2.11  2020/07/06
362 + -add extra output for DGP-contribution in -d mode
363 + -fix (remove) 1 pixel offset of the mulltiimage-mode
364 +  */
365 +
366 + /* evalglare.c, v2.12  2021/07/19
367 + -bugfix: center of glare source calculation now based on l*omega weighting, before it was just on omega.
368 +  */
369 +
370 + /* evalglare.c, v2.13  2021/07/21
371 + - center of glare source calculation now based on l*l*omega weighting.
372 + - bugfix: Position-index below line of sight is now corrected and based on the modified equation from Iwata and Osterhaus 2010
373 +  */
374 +
375 + /* evalglare.c, v2.14  2021/10/20
376 + - center of glare source calculation now based on l*omega weighting.
377 + - adding UGP2 equation to the output
378 + - new default setting: low-light correction off !
379 +  */
380 +
381 + /* evalglare.c, v3.0  2021/12/24
382 + - adding patch-mode (-z patchangle combineflag ) as glare source grouping method
383 +  */
384 +
385 + /* evalglare.c, v3.01  2022/01/05
386 + - change of multi-image-mode. Adding extra scans with scaled image
387 + syntax now: -Q n_images n_extrascans_per_image n_images*( imagename x y  n_extrascans*(scale) )
388 +  */
389 + /* evalglare.c, v3.02  2022/02/11
390 + - correction of position index below line of sight after Takuro Kikuchi's comments in radiance-discourse from Dec 2021.  
391 +  */
392 +
393 + /* evalglare.c, v3.03  2022/03/18, 2024/06/25
394 + - add no pixels output for zonal analysis.
395 +
396 + /* evalglare.c, v3.04  2023/09/01
397 + - test version to include peripherical dependency, changes not (yet) used for following versions
398 +
399 + /* evalglare.c, v3.05  2024/06/25
400 + - fix bug of still having lowlight correction even if turned off
401 + - changing abs_max,Lveil from float to double
402 +  */
403 +
404 + /* evalglare.c, v3.06  2024/06/25
405 + - fix bug missing border pixel when cutting with -G
406 +  */
407 +
408 +
409 +
410 +
411 + #define EVALGLARE
412   #define PROGNAME "evalglare"
413 < #define VERSION "1.17 release 15.07.2015 by EPFL, J.Wienold"
413 > #define VERSION "3.06 release 01.10.2025 by J.Wienold, EPFL"
414   #define RELEASENAME PROGNAME " " VERSION
415  
416 < #include "rtio.h"
288 < #include "platform.h"
416 >
417   #include "pictool.h"
418 + #include "rtio.h"
419   #include <math.h>
420   #include <string.h>
421 < char *progname;
421 > #include "platform.h"
422 > #include "muc_randvar.h"
423  
424   /* subroutine to add a pixel to a glare source */
425   void add_pixel_to_gs(pict * p, int x, int y, int gsn)
426   {
427          double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
428 <                new_omega, act_lum;
428 >                new_omega, act_lum,temp_av_posx,temp_av_posy,av_lum;
429  
430  
431          pict_get_npix(p, gsn) = pict_get_npix(p, gsn) + 1;
# Line 307 | Line 437 | void add_pixel_to_gs(pict * p, int x, int y, int gsn)
437          act_omega = pict_get_omega(p, x, y);
438          act_lum = luminance(pict_get_color(p, x, y));
439          new_omega = old_omega + act_omega;
440 <        pict_get_av_posx(p, gsn) =
441 <                (old_av_posx * old_omega + x * act_omega) / new_omega;
442 <        pict_get_av_posy(p, gsn) =
443 <                (old_av_posy * old_omega + y * act_omega) / new_omega;
440 >        pict_get_av_lum(p, gsn) =
441 >                (old_av_lum * old_omega + act_lum * act_omega) / new_omega;
442 >
443 >        av_lum=pict_get_av_lum(p, gsn);
444 >        temp_av_posx =
445 >                (old_av_posx *  old_omega* old_av_lum + x * act_omega*act_lum) / (old_av_lum*old_omega + act_lum* act_omega );
446 >        pict_get_av_posx(p, gsn) = temp_av_posx;
447 >        temp_av_posy =
448 >                (old_av_posy *  old_omega* old_av_lum + y * act_omega*act_lum) / (old_av_lum*old_omega + act_lum* act_omega );
449 >        
450 >        pict_get_av_posy(p, gsn) = temp_av_posy;
451          if (isnan((pict_get_av_posx(p, gsn))))
452                  fprintf(stderr,"error in add_pixel_to_gs %d %d %f %f %f %f\n", x, y, old_av_posy, old_omega, act_omega, new_omega);
453  
317        pict_get_av_lum(p, gsn) =
318                (old_av_lum * old_omega + act_lum * act_omega) / new_omega;
454          pict_get_av_omega(p, gsn) = new_omega;
455          pict_get_gsn(p, x, y) = gsn;
456          if (act_lum < pict_get_lum_min(p, gsn)) {
# Line 564 | Line 699 | double get_task_lum(pict * p, int x, int y, float r, i
699          return av_lum;
700   }
701  
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;
702  
703  
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);
704  
705  
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
706   /* subroutine for coloring the glare sources
707       red is used only for explicit call of the subroutine with acol=0  , e.g. for disability glare
708       -1: set to grey*/
# Line 626 | Line 714 | int setglcolor(pict * p, int x, int y, int acol, int u
714          l=u_r+u_g+u_b ;
715          
716          pcol[RED][0] = 1.0 / CIE_rf;
717 <        pcol[GRN][0] = 0.;
718 <        pcol[BLU][0] = 0.;
717 >        pcol[GRN][0] = 0.0 / CIE_gf;
718 >        pcol[BLU][0] = 0.0 / CIE_bf;
719  
720 <        pcol[RED][1] = 0.;
720 >        pcol[RED][1] = 0.0 / CIE_rf;
721          pcol[GRN][1] = 1.0 / CIE_gf;
722 <        pcol[BLU][1] = 0.;
722 >        pcol[BLU][1] = 0.0 / CIE_bf;
723  
724 <        pcol[RED][2] = 0.;
725 <        pcol[GRN][2] = 0.;
726 <        pcol[BLU][2] = 1 / CIE_bf;
724 >        pcol[RED][2] = 0.15 / CIE_rf;
725 >        pcol[GRN][2] = 0.15 / CIE_gf;
726 >        pcol[BLU][2] = 0.7 / CIE_bf;
727  
728 <        pcol[RED][3] = 1.0 / (1.0 - CIE_bf);
729 <        pcol[GRN][3] = 1.0 / (1.0 - CIE_bf);
730 <        pcol[BLU][3] = 0.;
728 >        pcol[RED][3] = .5 / CIE_rf;
729 >        pcol[GRN][3] = .5 / CIE_gf;
730 >        pcol[BLU][3] = 0.0 / CIE_bf;
731  
732 <        pcol[RED][4] = 1.0 / (1.0 - CIE_gf);
733 <        pcol[GRN][4] = 0.;
734 <        pcol[BLU][4] = 1.0 / (1.0 - CIE_gf);
732 >        pcol[RED][4] = .5 / CIE_rf;
733 >        pcol[GRN][4] = .0 / CIE_gf;
734 >        pcol[BLU][4] = .5 / CIE_bf ;
735  
736 <        pcol[RED][5] = 0.;
737 <        pcol[GRN][5] = 1.0 / (1.0 - CIE_rf);
738 <        pcol[BLU][5] = 1.0 / (1.0 - CIE_rf);
736 >        pcol[RED][5] = .0 / CIE_rf;
737 >        pcol[GRN][5] = .5 / CIE_gf;
738 >        pcol[BLU][5] = .5 / CIE_bf;
739  
740 <        pcol[RED][6] = 0.;
741 <        pcol[GRN][6] = 1.0 / (1.0 - CIE_rf);
742 <        pcol[BLU][6] = 1.0 / (1.0 - CIE_rf);
740 >        pcol[RED][6] = 0.333 / CIE_rf;
741 >        pcol[GRN][6] = 0.333 / CIE_gf;
742 >        pcol[BLU][6] = 0.333 / CIE_bf;
743  
744  
745          pcol[RED][999] = 1.0 / WHTEFFICACY;
# Line 662 | Line 750 | int setglcolor(pict * p, int x, int y, int acol, int u
750          pcol[RED][998] = u_r /(l* CIE_rf) ;
751          pcol[GRN][998] = u_g /(l* CIE_gf);
752          pcol[BLU][998] = u_b /(l* CIE_bf);
753 < /*        printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l);*/
753 > /*        printf("CIE_rf,CIE_gf,CIE_bf,l=%f %f %f %f %f\n",CIE_rf,CIE_gf,CIE_bf,l,WHTEFFICACY); */
754          icol = acol ;
755          if ( acol == -1) {icol=999;
756                                    }else{if (acol>0){icol = acol % 5 +1;
# Line 675 | Line 763 | int setglcolor(pict * p, int x, int y, int acol, int u
763          pict_get_color(p, x, y)[RED] = pcol[RED][icol] * act_lum / WHTEFFICACY;
764          pict_get_color(p, x, y)[GRN] = pcol[GRN][icol] * act_lum / WHTEFFICACY;
765          pict_get_color(p, x, y)[BLU] = pcol[BLU][icol] * act_lum / WHTEFFICACY;
766 <
766 > /*        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))); */
767          return icol;
768   }
769  
# Line 727 | Line 815 | void add_secondary_gs(pict * p, int x, int y, float r,
815  
816   /* color pixel of gs */
817  
818 <                icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b);
818 > /*              icol = setglcolor(p, x, y, act_gs, uniform_gs, u_r, u_g , u_b); */
819  
820          }
821   }
# Line 736 | Line 824 | void add_secondary_gs(pict * p, int x, int y, float r,
824   void split_pixel_from_gs(pict * p, int x, int y, int new_gsn, int uniform_gs, double u_r, double u_g , double u_b)
825   {
826          int old_gsn, icol;
827 <        double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega,
828 <                new_omega, act_lum;
827 >        double old_av_posx, old_av_posy, old_av_lum, old_omega, act_omega, delta_E,
828 >                new_omega, act_lum,temp_av_posx,temp_av_posy,av_lum;
829  
830  
831   /* change existing gs */
# Line 748 | Line 836 | void split_pixel_from_gs(pict * p, int x, int y, int n
836          old_av_posx = pict_get_av_posx(p, old_gsn);
837          old_av_posy = pict_get_av_posy(p, old_gsn);
838          old_omega = pict_get_av_omega(p, old_gsn);
839 +        
840 +        
841 +        
842  
843          new_omega = old_omega - act_omega;
844          pict_get_av_omega(p, old_gsn) = new_omega;
754        pict_get_av_posx(p, old_gsn) =
755                (old_av_posx * old_omega - x * act_omega) / new_omega;
756        pict_get_av_posy(p, old_gsn) =
757                (old_av_posy * old_omega - y * act_omega) / new_omega;
845  
759
846          act_lum = luminance(pict_get_color(p, x, y));
847          old_av_lum = pict_get_av_lum(p, old_gsn);
848          pict_get_av_lum(p, old_gsn) =
849                  (old_av_lum * old_omega - act_lum * act_omega) / new_omega;
850 +
851 +        av_lum = pict_get_av_lum(p, old_gsn);
852 +        temp_av_posx =
853 +                (old_av_posx *old_av_lum* old_omega  - x *act_lum* act_omega ) / (old_av_lum*old_omega  - act_lum* act_omega);
854 +
855 +        pict_get_av_posx(p, old_gsn) = temp_av_posx;
856 +        temp_av_posy =
857 +                (old_av_posy *old_av_lum* old_omega - y *act_lum* act_omega ) / (old_av_lum*old_omega  - act_lum* act_omega);
858 +        pict_get_av_posy(p, old_gsn) = temp_av_posy;
859 +
860          /* add pixel to new  gs */
861  
862          add_pixel_to_gs(p, x, y, new_gsn);
863  
768 /* color pixel of new gs */
864  
770        icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b);
865  
866  
867 + /* color pixel of new gs */
868 +
869 + /*      icol = setglcolor(p, x, y, new_gsn, uniform_gs, u_r, u_g , u_b); */
870 +        
871 +
872   }
873  
874   /* subroutine for the calculation of the position index */
875   float get_posindex(pict * p, float x, float y, int postype)
876   {
877          float posindex;
878 <        double teta, phi, sigma, tau, deg, d, s, r, fact;
878 >        double teta, beta, phi, sigma, tau, deg, d, s, r, fact;
879  
880  
881          pict_get_vangle(p, x, y, p->view.vdir, p->view.vup, &phi);
# Line 800 | Line 899 | float get_posindex(pict * p, float x, float y, int pos
899          tau = tau * deg;
900          sigma = sigma * deg;
901  
902 +        if (postype == 1) {
903 + /* KIM  model */        
904 +        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));
905 +        }else{
906 +
907 + /* Guth model, equation from IES lighting handbook */
908          posindex =
909                  exp((35.2 - 0.31889 * tau -
910                           1.22 * exp(-2 * tau / 9)) / 1000 * sigma + (21 +
# Line 807 | Line 912 | float get_posindex(pict * p, float x, float y, int pos
912                                                                                                                   0.002963 * tau *
913                                                                                                                   tau) / 100000 *
914                          sigma * sigma);
810 /* below line of sight, using Iwata model */
811        if (phi < 0) {
812                d = 1 / tan(phi);
813                s = tan(teta) / tan(phi);
814                r = sqrt(1 / d * 1 / d + s * s / d / d);
815                if (r > 0.6)
816                        fact = 1.2;
817                if (r > 3) {
818                        fact = 1.2;
819                        r = 3;
820                }
915  
916 <                posindex = 1 + fact * r;
916 > /* below line of sight, using Iwata model, CIE2010, converted coordinate system according to Takuro Kikuchi */
917 >
918 >        if (phi < 0) {
919 >          beta = atan(tan(sigma/deg)* sqrt(1 + 0.3225 * pow(cos(tau/deg),2))) * deg;
920 >          posindex = exp(6.49 / 1000 * beta + 21.0 / 100000 * beta * beta);
921 >
922          }
923 <        if (posindex > 16)
923 >
924 >                if (posindex > 16)
925                  posindex = 16;
926 + }
927  
928          return posindex;
929  
# Line 865 | Line 966 | void cut_view_1(pict*p)
966   {
967   int x,y;
968   double border,ang,teta,phi,phi2;
969 <                  for(x=0;x<pict_get_xsize(p)-1;x++)
970 <                   for(y=0;y<pict_get_ysize(p)-1;y++) {
969 >                  for(x=0;x<pict_get_xsize(p);x++)
970 >                   for(y=0;y<pict_get_ysize(p);y++) {
971                          if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
972                                  if (pict_is_validpixel(p, x, y)) {
973                                          pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
# Line 919 | Line 1020 | void cut_view_2(pict*p)
1020  
1021   int x,y;
1022   double border,ang,teta,phi,phi2;
1023 <                  for(x=0;x<pict_get_xsize(p)-1;x++)
1024 <                   for(y=0;y<pict_get_ysize(p)-1;y++) {
1023 >                  for(x=0;x<pict_get_xsize(p);x++)
1024 >                   for(y=0;y<pict_get_ysize(p);y++) {
1025                          if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
1026                                  if (pict_is_validpixel(p, x, y)) {
1027                                          pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
# Line 975 | Line 1076 | void cut_view_3(pict*p)
1076  
1077   int x,y;
1078   double border,ang,teta,phi,phi2,lum,newlum;
1079 <                  for(x=0;x<pict_get_xsize(p)-1;x++)
1080 <                   for(y=0;y<pict_get_ysize(p)-1;y++) {
1079 >                  for(x=0;x<pict_get_xsize(p);x++)
1080 >                   for(y=0;y<pict_get_ysize(p);y++) {
1081                          if (pict_get_hangle(p,x,y,p->view.vdir,p->view.vup,&ang)) {
1082                               if (DOT(pict_get_cached_dir(p,x,y),p->view.vdir) >= 0.0) {
1083                                          pict_get_vangle(p,x,y,p->view.vdir,p->view.vup,&phi2);
# Line 1035 | Line 1136 | return;
1136  
1137   }
1138  
1139 < /* subroutine for the calculation of the daylight glare index */
1139 > /* subroutine for the calculation of the daylight glare index DGI*/
1140  
1141  
1142   float get_dgi(pict * p, float lum_backg, int igs, int posindex_2)
# Line 1062 | Line 1163 | float get_dgi(pict * p, float lum_backg, int igs, int
1163          return dgi;
1164  
1165   }
1166 + /* subroutine for the calculation of the modified daylight glare index DGI_mod (Fisekis2003)
1167 +   several glare sources are taken into account and summed up, original equation has no summary
1168 + */
1169  
1170 +
1171 + float get_dgi_mod(pict * p, float lum_a, int igs, int posindex_2)
1172 + {
1173 +        float dgi_mod, sum_glare, omega_s;
1174 +        int i;
1175 +
1176 +
1177 +        sum_glare = 0;
1178 +        omega_s = 0;
1179 +
1180 +        for (i = 0; i <= igs; i++) {
1181 +
1182 +                if (pict_get_npix(p, i) > 0) {
1183 +                        omega_s =
1184 +                                pict_get_av_omega(p, i) / get_posindex(p,pict_get_av_posx(p,i), pict_get_av_posy(p,i),posindex_2) /
1185 +                                get_posindex(p, pict_get_av_posx(p, i), pict_get_av_posy(p, i), posindex_2);
1186 +                        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));
1187 +                    /*    printf("i,sum_glare %i %f\n",i,sum_glare); */
1188 +                }
1189 +        }
1190 +        dgi_mod = 10 * log10(sum_glare);
1191 +
1192 +        return dgi_mod;
1193 +
1194 + }
1195 +
1196   /* subroutine for the calculation of the daylight glare probability */
1197   double
1198   get_dgp(pict * p, double E_v, int igs, double a1, double a2, double a3,
# Line 1084 | Line 1214 | get_dgp(pict * p, double E_v, int igs, double a1, doub
1214                                                                                             pict_get_av_posy(p, i),
1215                                                                                             posindex_2),
1216                                                                    a4) * pow(pict_get_av_omega(p, i), a2);
1217 < /*                      printf("i,sum_glare %i %f\n",i,sum_glare);*/                                      
1217 >                        /*printf("i,sum_glare %i %f\n",i,sum_glare);    */                                
1218                          }
1219                  }
1220                  dgp =
# Line 1103 | Line 1233 | get_dgp(pict * p, double E_v, int igs, double a1, doub
1233   }
1234  
1235   /* subroutine for the calculation of the visual comfort probability */
1236 < float get_vcp(pict * p, double lum_a, int igs, int posindex_2)
1236 > float get_dgr(pict * p, double lum_a, int igs, int posindex_2)
1237   {
1238 <        float vcp;
1239 <        double sum_glare, dgr;
1238 >        float dgr;
1239 >        double sum_glare;
1240          int i, i_glare;
1241  
1242  
# Line 1132 | Line 1262 | float get_vcp(pict * p, double lum_a, int igs, int pos
1262          }
1263          dgr = pow(sum_glare, pow(i_glare, -0.0914));
1264  
1265 <        vcp = 50 * erf((6.347 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1265 >
1266 >
1267 >        return dgr;
1268 >
1269 > }
1270 >
1271 > float get_vcp(float dgr )
1272 > {
1273 >        float vcp;
1274 >
1275 >        vcp = 50 * erf((6.374 - 1.3227 * log(dgr)) / 1.414213562373) + 50;
1276          if (dgr > 750) {
1277                  vcp = 0;
1278          }
1279          if (dgr < 20) {
1280                  vcp = 100;
1281          }
1142
1143
1282          return vcp;
1283  
1284   }
1285  
1286 +
1287   /* subroutine for the calculation of the unified glare rating */
1288   float get_ugr(pict * p, double lum_backg, int igs, int posindex_2)
1289   {
# Line 1168 | Line 1307 | float get_ugr(pict * p, double lum_backg, int igs, int
1307                  }
1308          }
1309          ugr = 8 * log10(0.25 / lum_backg * sum_glare);
1310 <
1310 >        if (sum_glare==0) {
1311 >        ugr=0.0;
1312 >        }
1313 >        if (lum_backg<=0) {
1314 >        ugr=-99.0;
1315 >        }
1316 >        
1317          return ugr;
1318  
1319   }
1320 + /* subroutine for the calculation of the experimental unified glare rating (Fisekis 2003)*/
1321 + float get_ugr_exp(pict * p, double lum_backg, double lum_a, int igs, int posindex_2)
1322 + {
1323 +        float ugr_exp;
1324 +        double sum_glare;
1325 +        int i, i_glare;
1326  
1327  
1328 +        sum_glare = 0;
1329 +        i_glare = 0;
1330 +        for (i = 0; i <= igs; i++) {
1331 +
1332 +                if (pict_get_npix(p, i) > 0) {
1333 +                        i_glare = i_glare + 1;
1334 +                        sum_glare +=
1335 +                                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);
1336 +
1337 +                }
1338 +        }
1339 +        ugr_exp =8 * log10(lum_a) + 8 * log10(sum_glare/lum_backg);
1340 +
1341 +        return ugr_exp;
1342 +
1343 + }
1344 + /* subroutine for the calculation of the unified glare probability (Hirning 2014)*/
1345 + float get_ugp(pict * p, double lum_backg, int igs, int posindex_2)
1346 + {
1347 +        float ugp;
1348 +        double sum_glare;
1349 +        int i, i_glare;
1350 +
1351 +
1352 +        sum_glare = 0;
1353 +        i_glare = 0;
1354 +        for (i = 0; i <= igs; i++) {
1355 +
1356 +                if (pict_get_npix(p, i) > 0) {
1357 +                        i_glare = i_glare + 1;
1358 +                        sum_glare +=
1359 +                                pow(pict_get_av_lum(p, i) /
1360 +                                        get_posindex(p, pict_get_av_posx(p, i),
1361 +                                                                 pict_get_av_posy(p, i), posindex_2),
1362 +                                        2) * pict_get_av_omega(p, i);
1363 +
1364 +                }
1365 +        }
1366 +        ugp = 0.26 * log10(0.25 / lum_backg * sum_glare);
1367 +
1368 +        return ugp;
1369 +
1370 + }
1371 + /* subroutine for the calculation of the updated unified glare probability (Hirning 2016)*/
1372 + float get_ugp2(pict * p, double lum_backg, int igs, int posindex_2)
1373 + {
1374 +        float ugp;
1375 +        double sum_glare,ugp2;
1376 +        int i, i_glare;
1377 +
1378 +
1379 +        sum_glare = 0;
1380 +        i_glare = 0;
1381 +        for (i = 0; i <= igs; i++) {
1382 +
1383 +                if (pict_get_npix(p, i) > 0) {
1384 +                        i_glare = i_glare + 1;
1385 +                        sum_glare +=
1386 +                                pow(pict_get_av_lum(p, i) /
1387 +                                        get_posindex(p, pict_get_av_posx(p, i),
1388 +                                                                 pict_get_av_posy(p, i), posindex_2),
1389 +                                        2) * pict_get_av_omega(p, i);
1390 +
1391 +                }
1392 +        }
1393 +        ugp2 = 1/pow(1.0+2.0/7.0*pow(sum_glare/lum_backg,-0.2),10.0);
1394 +
1395 +        return ugp2;
1396 +
1397 + }
1398 +
1399 +
1400   /* subroutine for the calculation of the disability glare according to Poynter */
1401   float get_disability(pict * p, double lum_backg, int igs)
1402   {
# Line 1221 | Line 1444 | float get_disability(pict * p, double lum_backg, int i
1444  
1445   /* subroutine for the calculation of the cie glare index */
1446  
1447 < float
1225 < get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1447 > float get_cgi(pict * p, double E_v, double E_v_dir, int igs, int posindex_2)
1448   {
1449          float cgi;
1450          double sum_glare;
# Line 1246 | Line 1468 | get_cgi(pict * p, double E_v, double E_v_dir, int igs,
1468          cgi = 8 * log10((2 * (1 + E_v_dir / 500) / E_v) * sum_glare);
1469  
1470          return cgi;
1471 + }      
1472  
1473 +
1474 +
1475 + /* 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! */
1476 + float get_pgsv_con(double E_v, double E_mask,double omega_mask,double lum_mask_av, double Lavg)
1477 + {
1478 +        float pgsv_con;
1479 +        double Lb;
1480 +
1481 + /*        Lb = (E_v-E_mask)/3.14159265359;  */
1482 + /*        Lb = (2*E_v-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask); */
1483 +          Lb = (2*3.14159265359*Lavg-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask);
1484 +
1485 +
1486 +        pgsv_con =3.2*log10(lum_mask_av)-0.64*log10(omega_mask)+(0.79*log10(omega_mask)-0.61)*log10(Lb)-8.2 ;
1487 +
1488 +
1489 +        return pgsv_con;
1490   }
1491  
1492 + /* 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! */
1493 + float get_pgsv_sat(double E_v)
1494 + {
1495 +        float pgsv_sat;
1496  
1497 +        pgsv_sat =3.3-(0.57+3.3)/(1+pow(E_v/3.14159265359/1250,1.7));
1498 +
1499 +
1500 +        return pgsv_sat;
1501 + }
1502 +
1503 + /* 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! */
1504 +
1505 + float get_pgsv(double E_v, double E_mask,double omega_mask,double lum_mask_av,double Ltask, double Lavg)
1506 + {
1507 +        float pgsv;
1508 +        double Lb;
1509 +
1510 + /*        Lb = (E_v-E_mask)/3.14159265359;  */
1511 + /*        Lb = (2*E_v-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask); */
1512 +          Lb = (2*3.14159265359*Lavg-lum_mask_av*omega_mask)/(2*3.14159265359-omega_mask);
1513 +        
1514 +        if (Lb==0.0 ) {
1515 +               fprintf(stderr,  " warning: Background luminance is 0 or masking area = full image! pgsv cannot be calculated (set to -99)!!\n");
1516 +                pgsv=-99;
1517 +                        }else{
1518 +                if ( (lum_mask_av/Lb) > (E_v/(3.14159265359*Ltask))) {
1519 +                        pgsv=get_pgsv_con(E_v,E_mask,omega_mask,lum_mask_av, Lavg);
1520 +                }else{
1521 +                        pgsv=get_pgsv_sat(E_v)  ;
1522 +                        }}
1523 +        return pgsv;
1524 +
1525 + }
1526 +
1527 +
1528 +
1529   #ifdef  EVALGLARE
1530  
1531  
# Line 1259 | Line 1535 | get_cgi(pict * p, double E_v, double E_v_dir, int igs,
1535  
1536   int main(int argc, char **argv)
1537   {
1538 <        #define CLINEMAX 4095 /* memory allocated for command line string */
1538 >        #define CLINEMAX 999999999 /* memory allocated for command line string */
1539          pict *p = pict_create();
1540 <        int     skip_second_scan,calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, y1, fill, yfillmax, yfillmin,
1541 <                ext_vill, set_lum_max, set_lum_max2, task_color, i_splitstart,
1542 <                i_split, posindex_2, task_lum, checkfile, rval, i, i_max, x, y,x2,y2,
1543 <                igs, actual_igs, icol, xt, yt, change, before_igs, sgs, splithigh,uniform_gs,
1544 <                detail_out, posindex_picture, non_cos_lb, rx, ry, rmx, rmy,apply_disability,band_calc,band_color;
1545 <        double  lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue, lum_ideal, E_v_contr, sigma,
1546 <                E_vl_ext, lum_max, new_lum_max, r_center,
1547 <                search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,
1548 <                omegat, sang, E_v, E_v2, E_v_dir, avlum, act_lum, ang,
1549 <                l_max, lum_backg, lum_backg_cos, omega_sources, lum_sources,
1550 <                lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum;
1551 <        float lum_task, lum_thres, dgi,  vcp, cgi, ugr, limit,
1552 <                abs_max, Lveil;
1553 <        char file_out[500], file_out2[500], version[500];
1540 >        pict *pm = pict_create();
1541 >        pict *pcoeff = pict_create();
1542 >        int     add,i1,i2,y_lines,ix_offset,x_offset,y_offset, patchmode,ix, iy, patch_pixdistance_x, patch_pixdistance_y,nx_patch,ny_patch, lowlight,skip_second_scan,
1543 >                calcfast,age_corr,cut_view,cut_view_type,calc_vill, output, detail_out2, x1,y1, fill, yfillmax, yfillmin,
1544 >                ext_vill, set_lum_max, set_lum_max2, img_corr,x_disk,y_disk,task_color, i_splitstart,zones,act_gsn,splitgs,j,jj,multi_image_mode,
1545 >                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,num_images,xmap,ymap,
1546 >                igs, actual_igs, lastpixelwas_gs, icol, xt, yt, change,checkpixels, before_igs, sgs, splithigh,uniform_gs,x_max, y_max,y_mid,
1547 >                detail_out, posindex_picture, non_cos_lb, rx, ry,lastpixelwas_peak, rmx,rmy,apply_disability,band_calc,band_color,masking,i_mask,no_glaresources,force,num_scans;
1548 >        double  abs_max,Lveil,xxx,angle_v,angle_h,patch_angle, r_contrast,r_dgp,r_glare,sum_glare, LUM_replace,lum_total_max,age_corr_factor,age,dgp_ext,dgp,low_light_corr,omega_cos_contr, setvalue,
1549 >                lum_ideal, E_v_contr, sigma,om,delta_E,
1550 >                E_vl_ext, lum_max, new_lum_max, r_center, ugp,ugp2, ugr_exp, dgi_mod,lum_a, E_v_mask,angle_disk,dist,n_corner_px,zero_corner_px,
1551 >                search_pix, a1, a2, a3, a4, a5, c3, c1, c2, r_split, max_angle,r_actual,lum_actual,dir_ill,
1552 >                omegat, sang,Ez1,Ez2, E_v, E_v2, E_v_dir, avlum, act_lum, ang, angle_z1, angle_z2,per_95_band,per_75_band,pos,
1553 >                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,
1554 >                lum, lum_source,teta,Lveil_cie,Lveil_cie_sum,disability_thresh,u_r,u_g,u_b,band_angle,band_avlum,
1555 >                lum_mask[1],lum_mask_av,lum_mask_std[1],lum_mask_median[1],omega_mask,bbox[2],
1556 >                lum_band[1],lum_band_av,lum_band_std[1],lum_band_median[1],omega_band,bbox_band[2],
1557 >                lum_z1[1],lum_z1_av,lum_z1_std[1],lum_z1_median[1],omega_z1,bbox_z1[2],
1558 >                lum_z2[1],lum_z2_av,lum_z2_std[1],lum_z2_median[1],omega_z2,bbox_z2[2],
1559 >                lum_pos[1],lum_nopos_median[1],lum_pos_median[1],lum_pos2_median[1],lum_pos_mean,lum_pos2_mean;
1560 >        float lum_task, lum_thres, dgi,  vcp, cgi, ugr, limit, dgr,pgsv ,pgsv_sat,pgsv_con;
1561 >        char maskfile[500],file_out[500], file_out2[500], version[500],correction_type[500];
1562          char *cline;
1563 +        char** temp_image_name;
1564 +        int *x_temp_img, *y_temp_img;
1565 +        double *scale_image_scans;
1566 +        int *new_gs_number = NULL;
1567          VIEW userview = STDVIEW;
1568          int gotuserview = 0;
1569 +        
1570 +        struct muc_rvar* s_mask;
1571 +        s_mask = muc_rvar_create();
1572 +        muc_rvar_set_dim(s_mask, 1);
1573 +        muc_rvar_clear(s_mask);
1574 +        struct muc_rvar* s_band;
1575 +        s_band = muc_rvar_create();
1576 +        muc_rvar_set_dim(s_band, 1);
1577 +        muc_rvar_clear(s_band);
1578 +        struct muc_rvar* s_z1;
1579 +        s_z1 = muc_rvar_create();
1580 +        muc_rvar_set_dim(s_z1, 1);
1581 +        muc_rvar_clear(s_z1);
1582  
1583 +        struct muc_rvar* s_z2;
1584 +        s_z2 = muc_rvar_create();
1585 +        muc_rvar_set_dim(s_z2, 1);
1586 +        muc_rvar_clear(s_z2);
1587 +
1588 +        struct muc_rvar* s_noposweight;
1589 +        s_noposweight = muc_rvar_create();
1590 +        muc_rvar_set_dim(s_noposweight, 1);
1591 +        muc_rvar_clear(s_noposweight);
1592 +
1593 +        struct muc_rvar* s_posweight;
1594 +        s_posweight = muc_rvar_create();
1595 +        muc_rvar_set_dim(s_posweight, 1);
1596 +        muc_rvar_clear(s_posweight);
1597 +
1598 +        struct muc_rvar* s_posweight2;
1599 +        s_posweight2 = muc_rvar_create();
1600 +        muc_rvar_set_dim(s_posweight2, 1);
1601 +        muc_rvar_clear(s_posweight2);
1602 +
1603          /*set required user view parameters to invalid values*/
1604 <        uniform_gs = 0;
1604 >        patchmode=0;
1605 >        patch_angle=25.0;
1606 >        patch_pixdistance_x=0;
1607 >        patch_pixdistance_y=0;
1608 >        lowlight=0;
1609 >        multi_image_mode=0;
1610 >        lastpixelwas_peak=0;    
1611 >        num_images=0;        
1612 >        dir_ill=0.0;
1613 >        delta_E=0.0;
1614 >        no_glaresources=0;
1615 >        n_corner_px=0;
1616 >        zero_corner_px=0;
1617 >        force=0;
1618 >        dist=0.0;
1619 >        u_r=0.33333;
1620 >        u_g=0.33333;
1621 >        u_b=0.33333;
1622 >        band_angle=0;
1623 >        angle_z1=0;
1624 >        angle_z2=0;
1625 >        x_zone=0;
1626 >        y_zone=0;
1627 >        per_75_z2=0;
1628 >        per_95_z2=0;
1629 >        lum_pos_mean=0;
1630 >        lum_pos2_mean=0;
1631 >        lum_band_av = 0.0;
1632 >        omega_band = 0.0;
1633 >        pgsv = 0.0 ;
1634 >        pgsv_con = 0.0 ;
1635 >        pgsv_sat = 0.0 ;
1636 >        E_v_mask = 0.0;
1637 >        Ez1 = 0.0;
1638 >        Ez2 = 0.0;
1639 >        lum_z1_av = 0.0;
1640 >        omega_z1 = 0.0;
1641 >        lum_z2_av = 0.0;
1642 >        omega_z2 = 0.0 ;
1643 >        i_z1 = 0;
1644 >        i_z2 = 0;        
1645 >        zones = 0;
1646 >        lum_z2_av = 0.0;
1647 >        uniform_gs = 0;
1648          apply_disability = 0;
1649          disability_thresh = 0;
1650          Lveil_cie_sum=0.0;
# Line 1300 | Line 1664 | int main(int argc, char **argv)
1664          omegat = 0.0;
1665          yt = 0;
1666          xt = 0;
1667 +        x_disk=0;
1668 +        y_disk=0;
1669 +        angle_disk=0;
1670          yfillmin = 0;
1671          yfillmax = 0;
1672          cut_view = 0;
# Line 1308 | Line 1675 | int main(int argc, char **argv)
1675          omega_cos_contr = 0.0;
1676          lum_ideal = 0.0;
1677          max_angle = 0.2;
1678 <        lum_thres = 5.0;
1678 >        lum_thres = 2000.0;
1679 >        lum_task = 0.0;
1680          task_lum = 0;
1681          sgs = 0;
1682          splithigh = 1;
# Line 1326 | Line 1694 | int main(int argc, char **argv)
1694          c1 = 5.87e-05;
1695          c2 = 0.092;
1696          c3 = 0.159;
1697 <        non_cos_lb = 1;
1697 >        non_cos_lb = 0;
1698          posindex_2 = 0;
1699          task_color = 0;
1700          limit = 50000.0;
1701          set_lum_max = 0;
1702          set_lum_max2 = 0;
1703 +        img_corr=0;
1704          abs_max = 0;
1705 <        progname = argv[0];
1705 >        fixargv0(argv[0]);
1706          E_v_contr = 0.0;
1707 <        strcpy(version, "1.15 release 05.02.2015 by J.Wienold");
1707 >        strcpy(version, "3.06 release 01.10.2025 by J.Wienold");
1708          low_light_corr=1.0;
1709          output = 0;
1710          calc_vill = 0;
1711          band_avlum = -99;
1712          band_calc = 0;
1713 +        masking = 0;
1714 +        lum_mask[0]=0.0;
1715 +        lum_mask_av=0.0;
1716 +        omega_mask=0.0;
1717 +        i_mask=0;
1718 +        actual_igs=0;
1719 +        LUM_replace=0;
1720 +        thres_activate=0;
1721   /* command line for output picture*/
1722  
1723 <        cline = (char *) malloc(CLINEMAX+1);
1723 >        cline = (char *) malloc(CLINEMAX+100);
1724          cline[0] = '\0';
1725          for (i = 0; i < argc; i++) {
1726   /*       fprintf(stderr, "%d %d \n",i,strlen(argv[i]));*/
# Line 1382 | Line 1759 | int main(int argc, char **argv)
1759                          printf("age factor not supported any more \n");
1760                          exit(1);
1761                          break;
1762 +                case 'A':
1763 +                        masking = 1;
1764 +                        detail_out = 1;
1765 +                        strcpy(maskfile, argv[++i]);
1766 +                        pict_read(pm,maskfile );
1767 +
1768 +                        break;
1769                  case 'b':
1770                          lum_thres = atof(argv[++i]);
1771 +                        lum_source =lum_thres;
1772 +                        thres_activate = 1;
1773                          break;
1774                  case 'c':
1775                          checkfile = 1;
# Line 1398 | Line 1784 | int main(int argc, char **argv)
1784                  case 'r':
1785                          max_angle = atof(argv[++i]);
1786                          break;
1787 +                case 'z':
1788 +                        patch_angle = atof(argv[++i]);
1789 +                        patchmode= atoi(argv[++i]);
1790 +                        if ( patchmode == 3) { output=4;}
1791 +                        
1792 +                        /* patchmode = 0 : deactivated; patchmode =1: patch without combining, normal output; patchmode =3: patch without combining, coefficient output; patchmode =2 patch with combining (not yet implemented...) */
1793 +                        break;
1794 +
1795                  case 's':
1796                          sgs = 1;
1797                          break;
1798 +                case 'f':
1799 +                        force = 1;
1800 +                        break;
1801                  case 'k':
1802                          apply_disability = 1;
1803                          disability_thresh = atof(argv[++i]);
# Line 1409 | Line 1806 | int main(int argc, char **argv)
1806                  case 'p':
1807                          posindex_picture = 1;
1808                          break;
1809 +                case 'P':
1810 +                        posindex_2 = atoi(argv[++i]);
1811 +                        break;
1812 +                        
1813  
1814                  case 'y':
1815                          splithigh = 1;
# Line 1439 | Line 1840 | int main(int argc, char **argv)
1840                          detail_out2 = 1;
1841                          break;
1842                  case 'm':
1843 +                        img_corr = 1;
1844                          set_lum_max = 1;
1845                          lum_max = atof(argv[++i]);
1846                          new_lum_max = atof(argv[++i]);
# Line 1447 | Line 1849 | int main(int argc, char **argv)
1849                          break;
1850  
1851                  case 'M':
1852 +                        img_corr = 1;
1853                          set_lum_max2 = 1;
1854                          lum_max = atof(argv[++i]);
1855                          E_vl_ext = atof(argv[++i]);
1856                          strcpy(file_out2, argv[++i]);
1857   /*                      printf("max lum set to %f\n",new_lum_max);*/
1858                          break;
1859 <                case 'n':
1859 >                case 'N':
1860 >                        img_corr = 1;
1861 >                        set_lum_max2 = 2;
1862 >                        x_disk = atoi(argv[++i]);
1863 >                        y_disk = atoi(argv[++i]);
1864 >                        angle_disk = atof(argv[++i]);
1865 >                        E_vl_ext = atof(argv[++i]);
1866 >                        strcpy(file_out2, argv[++i]);
1867 > /*                      printf("max lum set to %f\n",new_lum_max);*/
1868 >                        break;
1869 >                case 'O':
1870 >                        img_corr = 1;
1871 >                        set_lum_max2 = 3;
1872 >                        x_disk = atoi(argv[++i]);
1873 >                        y_disk = atoi(argv[++i]);
1874 >                        angle_disk = atof(argv[++i]);
1875 >                        LUM_replace = atof(argv[++i]);
1876 >                        strcpy(file_out2, argv[++i]);
1877 > /*                      printf("max lum set to %f\n",new_lum_max);*/
1878 >                        break;
1879 >
1880 >
1881 > /* deactivated          case 'n':
1882                          non_cos_lb = 0;
1883                          break;
1884 + */
1885 +                case 'q':
1886 +                        non_cos_lb = atoi(argv[++i]);
1887 +                        break;
1888  
1889                  case 't':
1890                          task_lum = 1;
# Line 1472 | Line 1901 | int main(int argc, char **argv)
1901                          omegat = atof(argv[++i]);
1902                          task_color = 1;
1903                          break;
1904 +                case 'l':
1905 +                        zones = 1;
1906 +                        x_zone = atoi(argv[++i]);
1907 +                        y_zone = atoi(argv[++i]);
1908 +                        angle_z1 = atof(argv[++i]);
1909 +                        angle_z2 = -1;
1910 +                        break;
1911 +                case 'L':
1912 +                        zones = 2;
1913 +                        x_zone = atoi(argv[++i]);
1914 +                        y_zone = atoi(argv[++i]);
1915 +                        angle_z1 = atof(argv[++i]);
1916 +                        angle_z2 = atof(argv[++i]);
1917 +                        break;
1918 +                        
1919 +                        
1920                  case 'B':
1921                          band_calc = 1;
1922                          band_angle = atof(argv[++i]);
# Line 1504 | Line 1949 | int main(int argc, char **argv)
1949                          /*case 'v':
1950                          printf("evalglare  %s \n",version);
1951                          exit(1); */
1952 +                case 'C':
1953 +                        strcpy(correction_type,argv[++i]);
1954 +                        
1955 +                        if (!strcmp(correction_type,"l-")  ){
1956 +                /*      printf("low light off!\n"); */
1957 +                                                                                        lowlight = 0; }
1958 +                        if (!strcmp(correction_type,"l+")  ){
1959 +                /*      printf("low light on!\n"); */
1960 +                                                                                        lowlight = 1; }
1961 +                        if (!strcmp(correction_type,"0")  ){
1962 +                /*      printf("all corrections off!\n"); */
1963 +                                                                                        lowlight = 0; }
1964 +                        
1965 +                        break;
1966  
1967 +                        /*case 'v':
1968 +                        printf("evalglare  %s \n",version);
1969 +                        exit(1); */
1970 +
1971                  case '1':
1972                          output = 1;
1973                          break;
1974 <
1974 >                case '2':
1975 >                        output = 2;
1976 >                        dir_ill = atof(argv[++i]);
1977 >                        break;
1978 >                case '3':
1979 >                        output = 3;
1980 >                        break;
1981 >                case '4':
1982 >                        lowlight = 0;
1983 >                        break;
1984 >                case 'Q':
1985 >                        multi_image_mode=1;
1986 >                        output= 3;                      
1987 >                        calcfast=1;
1988 >                        num_images =atoi(argv[++i]);
1989 >                        num_scans =atoi(argv[++i]);
1990 >                        temp_image_name = malloc(sizeof(char*)*num_images);
1991 >                        
1992 >                        x_temp_img=(int *) malloc(sizeof(int) * num_images);
1993 >                        y_temp_img=(int *) malloc(sizeof(int) * num_images);
1994 >                        scale_image_scans=(double *) malloc(sizeof(double) * (num_scans+1)*num_images);
1995 >        /* iterate through all images and allocate 256 characters to each: */
1996 >                        for (j = 0; j < num_images; j++) {
1997 >                                scale_image_scans[j*(num_scans+1)]=1.0;
1998 >                                temp_image_name[j] = malloc(256*sizeof(char));
1999 >                                strcpy(temp_image_name[j], argv[++i]);
2000 >                                x_temp_img[j] = atoi(argv[++i]);
2001 >                                y_temp_img[j] = atoi(argv[++i]);
2002 >                                for (jj=1;jj<=num_scans;jj++){
2003 >                                      scale_image_scans[j*(num_scans+1)+jj]=atof(argv[++i]);}
2004 >                        }
2005 >                
2006 >                
2007 >                        break;
2008                  case 'v':
2009                          if (argv[i][2] == '\0') {
2010                               printf("%s \n",RELEASENAME);                              
# Line 1537 | Line 2033 | int main(int argc, char **argv)
2033                  }
2034          }
2035  
2036 + /* set multiplier for task method to 5, if not specified */
2037 +                
2038 +
2039 +
2040 + if ( task_lum == 1 && thres_activate == 0){
2041 +                lum_thres = 5.0;
2042 + }
2043   /*fast calculation, if gendgp_profile is used: No Vertical illuminance calculation, only dgp is calculated*/
2044  
2045 < if (output == 1 && ext_vill == 1) {
2045 > if (output == 1 && ext_vill == 1 ) {
2046                         calcfast=1;
2047                         }
2048 +                      
2049 + if (output == 2 && ext_vill == 1 ) {
2050 +                       calcfast=2;
2051 +                       }
2052 +                      
2053 + /*masking and zoning cannot be applied at the same time*/
2054  
2055 + if (masking ==1 && zones >0) {
2056 +               fprintf(stderr,  " masking and zoning cannot be activated at the same time!\n");
2057 +               exit(1);
2058 + }
2059 +
2060   /* read picture file */
2061          if (i == argc) {
2062                  SET_FILE_BINARY(stdin);
# Line 1576 | Line 2090 | if (output == 1 && ext_vill == 1) {
2090                  return EXIT_FAILURE;
2091          }
2092  
2093 +
2094 +
2095   /*                fprintf(stderr,"Picture read!");*/
2096  
2097   /* command line for output picture*/
2098  
2099          pict_set_comment(p, cline);
2100  
2101 + /* several checks */
2102          if (p->view.type == VT_PAR) {
2103 <                fprintf(stderr, "wrong view type! must not be parallel ! \n");
2103 >                fprintf(stderr, "error: wrong view type! must not be parallel ! \n");
2104                  exit(1);
2105          }
2106  
2107 +        if ( patchmode > 0 && p->view.type != VT_ANG) {
2108 +        
2109 +                fprintf(stderr, "error: Patchmode only possible with view type vta  !  Stopping... \n");
2110 +                exit(1);
2111 +        
2112 +        }
2113 +        
2114 +        
2115 +        if (p->view.type == VT_PER) {
2116 +                fprintf(stderr, "warning: image has perspective view type specified in header ! \n");
2117 +        }
2118  
2119 +        if (masking == 1) {
2120 +
2121 +                if (pict_get_xsize(p)!=pict_get_xsize(pm) || pict_get_ysize(p)!=pict_get_ysize(pm)) {
2122 +                fprintf(stderr, "error: masking image has other resolution than main image ! \n");
2123 +                fprintf(stderr, "size must be identical \n");
2124 +                printf("resolution main image : %dx%d\n",pict_get_xsize(p),pict_get_ysize(p));
2125 +                printf("resolution masking image : %dx%d\n",pict_get_xsize(pm),pict_get_ysize(pm));
2126 +                exit(1);
2127 +                
2128 +                }
2129 +
2130 +        }
2131 + /*      printf("resolution: %dx%d\n",pict_get_xsize(p),pict_get_ysize(p)); */
2132 +
2133   /* Check size of search radius */
2134          rmx = (pict_get_xsize(p) / 2);
2135          rmy = (pict_get_ysize(p) / 2);
# Line 1611 | Line 2153 | if (output == 1 && ext_vill == 1) {
2153  
2154                  }
2155          }
2156 <
2156 >        
2157 >
2158   /* Check task position  */
2159  
2160          if (task_lum == 1) {
# Line 1632 | Line 2175 | if (output == 1 && ext_vill == 1) {
2175          avlum = 0.0;
2176          pict_new_gli(p);
2177          igs = 0;
2178 +        pict_get_z1_gsn(p,igs) = 0;
2179 +        pict_get_z2_gsn(p,igs) = 0;
2180  
2181 + if (multi_image_mode<1) {
2182 +
2183 +
2184   /* cut out GUTH field of view and exit without glare evaluation */
2185   if (cut_view==2) {
2186          if (cut_view_type==1) {
# Line 1699 | Line 2247 | if (cut_view==2) {
2247          if (calcfast == 0) {
2248          for (x = 0; x < pict_get_xsize(p); x++)
2249                  for (y = 0; y < pict_get_ysize(p); y++) {
2250 +                   lum = luminance(pict_get_color(p, x, y));              
2251 +                   dist=sqrt((x-rmx)*(x-rmx)+(y-rmy)*(y-rmy));
2252 +                  if (dist>((rmx+rmy)/2)) {
2253 +                     n_corner_px=n_corner_px+1;
2254 +                     if (lum<7.0) {
2255 +                         zero_corner_px=zero_corner_px+1;
2256 +                         }
2257 +                     }
2258 +                
2259 +                
2260                          if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2261                                  if (pict_is_validpixel(p, x, y)) {
1704                                        lum = luminance(pict_get_color(p, x, y));
2262                                          if (fill == 1 && y >= yfillmax) {
2263                                                  y1 = y - 1;
2264                                                  lum = luminance(pict_get_color(p, x, y1));
# Line 1714 | Line 2271 | if (cut_view==2) {
2271                                                  abs_max = lum;
2272                                          }
2273   /* set luminance restriction, if -m is set */
2274 <                                        if (set_lum_max == 1 && lum >= lum_max) {
2275 <                                                pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
2276 <                                                pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
2277 <                                                pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
2278 < /*                                    printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
2279 <                                                lum = luminance(pict_get_color(p, x, y));
2274 >                                        if (img_corr == 1 ) {
2275 >                                                if (set_lum_max == 1 && lum >= lum_max) {
2276 >                                                        pict_get_color(p, x, y)[RED] = new_lum_max / 179.0;
2277 >                                                        pict_get_color(p, x, y)[GRN] = new_lum_max / 179.0;
2278 >                                                        pict_get_color(p, x, y)[BLU] = new_lum_max / 179.0;
2279 > /*                                                      printf("old lum, new lum %f %f \n",lum,luminance(pict_get_color(p,x,y))); */
2280 >                                                        lum = luminance(pict_get_color(p, x, y));
2281 >                                                        }
2282 >                                                if (set_lum_max2 == 1 && lum >= lum_max) {
2283  
2284 <                                        }
2285 <                                        if (set_lum_max2 == 1 && lum >= lum_max) {
2284 >                                                        E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2285 >                                                        omega_cos_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2286 >                                                        }
2287 >                                                if (set_lum_max2 == 2 ) {
2288 >                                                        r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2289 >                                                        if (x_disk == x && y_disk==y ) r_actual=0.0;
2290  
2291 <                                                E_v_contr +=
2292 <                                                        DOT(p->view.vdir,
2293 <                                                                pict_get_cached_dir(p, x,
2294 <                                                                                                        y)) * pict_get_omega(p,
2295 <                                                                                                                                                 x,
2296 <                                                                                                                                                 y)
2297 <                                                        * lum;
2298 <                                                omega_cos_contr +=
2299 <                                                        DOT(p->view.vdir,
2300 <                                                                pict_get_cached_dir(p, x,
1737 <                                                                                                        y)) * pict_get_omega(p,
1738 <                                                                                                                                                 x,
1739 <                                                                                                                                                 y)
1740 <                                                        * 1;
2291 >                                                        act_lum = luminance(pict_get_color(p, x, y));
2292 >
2293 >                                                       if (r_actual <= angle_disk) {
2294 >                                                              E_v_contr +=DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* lum;
2295 >                                                              omega_cos_contr += DOT(p->view.vdir,pict_get_cached_dir(p, x,y)) * pict_get_omega(p,x,y)* 1;
2296 >                                                                }
2297 >
2298 >                                                
2299 >                                                
2300 >                                                        }
2301                                          }
2302 +                                        om = pict_get_omega(p, x, y);
2303 +                                        sang += om;
2304 +                                        E_v += DOT(p->view.vdir, pict_get_cached_dir(p, x,y)) * om *lum;
2305 +                                        avlum += om * lum;
2306 +                                        pos=get_posindex(p, x, y,posindex_2);
2307  
2308 +                                        lum_pos[0]=lum;
2309 +                                        muc_rvar_add_sample(s_noposweight,lum_pos);
2310 +                                        lum_pos[0]=lum/pos;
2311 +                                        lum_pos_mean+=lum_pos[0]*om;
2312 +                                        muc_rvar_add_sample(s_posweight, lum_pos);
2313 +                                        lum_pos[0]=lum_pos[0]/pos;
2314 +                                        lum_pos2_mean+=lum_pos[0]*om;
2315 +                                        muc_rvar_add_sample(s_posweight2,lum_pos);
2316  
2317 <                                        sang += pict_get_omega(p, x, y);
2318 <                                        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));
2317 >
2318 >
2319                                  } else {
2320                                          pict_get_color(p, x, y)[RED] = 0.0;
2321                                          pict_get_color(p, x, y)[GRN] = 0.0;
2322                                          pict_get_color(p, x, y)[BLU] = 0.0;
2323  
2324                                  }
2325 <                        }
2325 >                        }else {
2326 >                                        pict_get_color(p, x, y)[RED] = 0.0;
2327 >                                        pict_get_color(p, x, y)[GRN] = 0.0;
2328 >                                        pict_get_color(p, x, y)[BLU] = 0.0;
2329 >
2330 >                                }
2331                  }
2332  
2333 <        if (set_lum_max2 == 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0) {
2333 > /* check if image has black corners AND a perspective view */
2334  
2335 +       if (zero_corner_px/n_corner_px > 0.70 && (p->view.type == VT_PER) ) {
2336 +       printf (" corner pixels are to  %f %% black! \n",zero_corner_px/n_corner_px*100);
2337 +                printf("error! The image has black corners AND header contains a perspective view type definition !!!\n") ;
2338 +                
2339 +                if (force==0) {
2340 +                                printf("stopping...!!!!\n") ;
2341 +
2342 +                      exit(1);
2343 +                 }
2344 +       }
2345 +
2346 +
2347 +
2348 +
2349 +        lum_pos_mean= lum_pos_mean/sang;
2350 +        lum_pos2_mean= lum_pos2_mean/sang;
2351 +
2352 +        if ((set_lum_max2 >= 1 && E_v_contr > 0 && (E_vl_ext - E_v) > 0 ) || set_lum_max2==3) {
2353 +
2354 +                if (set_lum_max2<3){
2355                  lum_ideal = (E_vl_ext - E_v + E_v_contr) / omega_cos_contr;
2356 +                if (set_lum_max2 == 2 && lum_ideal >= 2e9) {
2357 +                printf("warning! luminance of replacement pixels would be larger than 2e9 cd/m2. Value set to 2e9cd/m2!\n") ;
2358 +                }
2359 +
2360                  if (lum_ideal > 0 && lum_ideal < setvalue) {
2361                          setvalue = lum_ideal;
2362                  }
2363 <                printf
1771 <                        ("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f  %f %f %f %f\n",
2363 >                printf("change luminance values!! lum_ideal,setvalue,E_vl_ext,E_v,E_v_contr %f  %f %f %f %f\n",
2364                           lum_ideal, setvalue, E_vl_ext, E_v, E_v_contr);
2365 +                          }else{setvalue=LUM_replace;
2366 +                         }
2367  
2368 <
2368 >            
2369                  for (x = 0; x < pict_get_xsize(p); x++)
2370                          for (y = 0; y < pict_get_ysize(p); y++) {
2371                                  if (pict_get_hangle
2372                                          (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2373                                          if (pict_is_validpixel(p, x, y)) {
2374                                                  lum = luminance(pict_get_color(p, x, y));
2375 +                                                
2376 +                                                
2377                                                  if (set_lum_max2 == 1 && lum >= lum_max) {
2378  
2379                                                          pict_get_color(p, x, y)[RED] =
# Line 1787 | Line 2383 | if (cut_view==2) {
2383                                                          pict_get_color(p, x, y)[BLU] =
2384                                                                  setvalue / 179.0;
2385  
2386 +                                                }else{ if(set_lum_max2 >1 ) {
2387 +                                                        r_actual =acos(DOT(pict_get_cached_dir(p, x_disk, y_disk), pict_get_cached_dir(p, x, y))) * 2;
2388 +                                                        if (x_disk == x && y_disk==y ) r_actual=0.0;
2389 +
2390 +                                                       if (r_actual <= angle_disk) {
2391 +                                                      
2392 +                                                        pict_get_color(p, x, y)[RED] =
2393 +                                                                setvalue / 179.0;
2394 +                                                        pict_get_color(p, x, y)[GRN] =
2395 +                                                                setvalue / 179.0;
2396 +                                                        pict_get_color(p, x, y)[BLU] =
2397 +                                                                setvalue / 179.0;
2398 +                                                      
2399 +                                                       }                                                
2400                                                  }
2401 +                                                }
2402                                          }
2403                                  }
2404                          }
2405 +                        
2406  
2407                  pict_write(p, file_out2);
2408 <
2408 >        exit(1);
2409          }
2410          if (set_lum_max == 1) {
2411                  pict_write(p, file_out2);
# Line 1853 | Line 2465 | if (cut_view==1) {
2465                  lum_source = lum_thres;
2466          }
2467  
2468 <        /*     printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2468 > /*     printf("Ev, avlum, Omega_tot%f %f %f \n",E_v, avlum, sang); */
2469  
2470   /* first glare source scan: find primary glare sources */
2471 <        for (x = 0; x < pict_get_xsize(p); x++)
2471 >
2472 > if (patchmode==0) {
2473 >        for (x = 0; x < pict_get_xsize(p); x++) {
2474 >                lastpixelwas_gs=0;
2475 > /*              lastpixelwas_peak=0; */
2476                  for (y = 0; y < pict_get_ysize(p); y++) {
2477                          if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2478                                  if (pict_is_validpixel(p, x, y)) {
# Line 1865 | Line 2481 | if (cut_view==1) {
2481                                                  if (act_lum >lum_total_max) {
2482                                                                               lum_total_max=act_lum;
2483                                                                                 }
2484 <                                                
2485 <                                                actual_igs =
2486 <                                                        find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2484 > /* speed improvement first scan: when last pixel was glare source, then use this glare source number instead of searching.
2485 >   has been deactivated with v1.25, reactivated with v2.10 */
2486 >                      
2487 >                                                if (lastpixelwas_gs==0 || search_pix <= 1.0 ) {
2488 >                                                actual_igs = find_near_pgs(p, x, y, max_angle, 0, igs, 1);
2489 >  }
2490                                                  if (actual_igs == 0) {
2491                                                          igs = igs + 1;
2492                                                          pict_new_gli(p);
2493                                                          pict_get_lum_min(p, igs) = HUGE_VAL;
2494                                                          pict_get_Eglare(p,igs) = 0.0;
2495                                                          pict_get_Dglare(p,igs) = 0.0;
2496 +                                                        pict_get_z1_gsn(p,igs) = 0;
2497 +                                                        pict_get_z2_gsn(p,igs) = 0;
2498                                                          actual_igs = igs;
2499 +                                                        
2500                                                  }
2501 <                                                icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2501 > /* no coloring any more here            icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2502                                                  pict_get_gsn(p, x, y) = actual_igs;
2503                                                  pict_get_pgs(p, x, y) = 1;
2504                                                  add_pixel_to_gs(p, x, y, actual_igs);
2505 +                                                lastpixelwas_gs=actual_igs;
2506  
2507                                          } else {
2508                                                  pict_get_pgs(p, x, y) = 0;
2509 +                                                lastpixelwas_gs=0;
2510                                          }
2511                                  }
2512                          }
2513 <                }
2514 < /*      pict_write(p,"firstscan.pic");   */
2513 >                }
2514 >             }
2515 >         } else {    
2516 > /* patchmode on!
2517 > calculation only for angular projection!
2518  
2519 < if (calcfast == 1 ) {
2519 > */
2520 >
2521 >
2522 > angle_v=p->view.vert;
2523 > angle_h=p->view.horiz;
2524 >
2525 >
2526 > /*printf ("angle_v,angle_h: %f %f \n",angle_v,angle_h) ;
2527 >
2528 >  setting up patches as glare sources */
2529 >
2530 > patch_pixdistance_x=floor(pict_get_xsize(p)/angle_h*patch_angle);
2531 > patch_pixdistance_y=floor(pict_get_ysize(p)/angle_v*patch_angle);
2532 >
2533 > nx_patch=floor(angle_v/patch_angle)+1;
2534 > ny_patch=floor(angle_h/patch_angle)+1;
2535 >
2536 > y_offset=floor (patch_pixdistance_y/2);
2537 > x_offset=floor (patch_pixdistance_x/2);
2538 > /* printf ("nx_patch,ny_patch,x_offset,y_offset,patch_pixdistance_x,patch_pixdistance_y   %i %i %i %i %i %i\n",nx_patch,ny_patch,x_offset,y_offset,patch_pixdistance_x,patch_pixdistance_y) ; */
2539 >
2540 > ix_offset=0;
2541 > for (iy=1; iy<=ny_patch;iy++) {
2542 >
2543 >
2544 >
2545 > for (ix=1; ix<=nx_patch;ix++) {
2546 >                                                        igs = igs + 1;
2547 >                                                        pict_new_gli(p);
2548 >
2549 >                                                        pict_get_lum_min(p, igs) = HUGE_VAL;
2550 >                                                        pict_get_Eglare(p,igs) = 0.0;
2551 >                                                        pict_get_Dglare(p,igs) = 0.0;
2552 >                                                        pict_get_z1_gsn(p,igs) = 0;
2553 >                                                        pict_get_z2_gsn(p,igs) = 0;
2554 >                                                        pict_get_dx_max(p,igs) = (x_offset+ix_offset*x_offset+(ix-1)*patch_pixdistance_x)*1.0;
2555 >                                                        pict_get_dy_max(p,igs) = (y_offset+(iy-1)*patch_pixdistance_y)*1.0;
2556 >
2557 > /* printf ("igs, x-patch, y-patch : %i %f %f \n",igs,pict_get_dx_max(p,igs),pict_get_dy_max(p,igs) ) ; */
2558 >
2559 > }
2560 > ix_offset=ix_offset+1;
2561 > if (ix_offset==2) {
2562 >                        ix_offset =0 ;
2563 >                        }
2564 >
2565 > }
2566 >                for (y = 0; y < pict_get_ysize(p); y++) {
2567 >                        for (x = 0; x < pict_get_xsize(p); x++) {
2568 >
2569 >                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2570 >                                if (pict_is_validpixel(p, x, y)) {
2571 >                                        act_lum = luminance(pict_get_color(p, x, y));
2572 >                                        if (act_lum > lum_source) {
2573 >                                                if (act_lum >lum_total_max) {
2574 >                                                                             lum_total_max=act_lum;
2575 >                                                                               }
2576 >                                        
2577 >                                                y_lines = floor((y)/patch_pixdistance_y);
2578 >                                                xxx = (x+0.0)/(patch_pixdistance_x+0.0)-0.5*(y_lines % 2);
2579 >                                                if (xxx<0 ) { xxx=0.0 ;}
2580 >                                                i1 = y_lines*(nx_patch)+floor(xxx)+1;
2581 >                                                i2=0;
2582 >                                                add=0;
2583 >                                                if (y_lines % 2 == 1 ) {add=1;}
2584 >                                                
2585 >                                                if (y >pict_get_dy_max(p,i1)) {
2586 >                                                        
2587 >                                                        if (x > pict_get_dx_max(p,i1)) {
2588 >                                                                i2=i1+nx_patch+add;
2589 >                                                                }else {
2590 >                                                                i2=i1+nx_patch-1+add;
2591 >                                                                }
2592 >                                                        }else {
2593 >                                                
2594 >                                                        if (x > pict_get_dx_max(p,i1)) {
2595 >                                                                i2=i1-nx_patch+add;
2596 >                                                                }else {
2597 >                                                                i2=i1-nx_patch-1+add;
2598 >                                                                }
2599 >                                                        }
2600 >                                                
2601 >                                                
2602 >        
2603 >                                                
2604 >                                                if (i2 > igs || i2 < 1) {actual_igs=i1;}else{
2605 >                                                if ( ((x-pict_get_dx_max(p,i1))*(x-pict_get_dx_max(p,i1))+(y-pict_get_dy_max(p,i1))*(y-pict_get_dy_max(p,i1))) < ((x-pict_get_dx_max(p,i2))*(x-pict_get_dx_max(p,i2))+(y-pict_get_dy_max(p,i2))*(y-pict_get_dy_max(p,i2))) ) {
2606 >                                                 actual_igs=i1; }else {actual_igs=i2;}}
2607 >                                                
2608 >        
2609 >                                                pict_get_gsn(p, x, y) = actual_igs;
2610 >                                                pict_get_pgs(p, x, y) = 1;
2611 >                                                add_pixel_to_gs(p, x, y, actual_igs);
2612 >                        /*                      setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2613 >                                }
2614 >                                }
2615 >                                }
2616 >
2617 >
2618 >
2619 >
2620 >
2621 > }               }
2622 >
2623 >
2624 >
2625 >            
2626 >             }
2627 >  /*                    pict_write(p,"firstscan.pic");   */
2628 >
2629 >
2630 >
2631 >
2632 > if (calcfast ==1 || search_pix <= 1.0 || calcfast == 2 || patchmode > 0) {
2633     skip_second_scan=1;
2634     }
2635 +  
2636  
2637   /* second glare source scan: combine glare sources facing each other */
2638          change = 1;
2639 +        i = 0;
2640          while (change == 1 && skip_second_scan==0) {
2641                  change = 0;
2642 <                for (x = 0; x < pict_get_xsize(p); x++)
2642 >                i = i+1;
2643 >                for (x = 0; x < pict_get_xsize(p); x++) {
2644                          for (y = 0; y < pict_get_ysize(p); y++) {
1902                                before_igs = pict_get_gsn(p, x, y);
2645                                  if (pict_get_hangle
2646                                          (p, x, y, p->view.vdir, p->view.vup, &ang)) {
2647 <                                        if (pict_is_validpixel(p, x, y) && before_igs > 0) {
2647 >                                        checkpixels=1;
2648 >                                        before_igs = pict_get_gsn(p, x, y);
2649 >
2650 > /* speed improvement: search for other nearby glare sources only, when at least one adjacend pixel has another glare source number  */
2651 >                                        if (x >1 && x<pict_get_xsize(p)-2 && y >1 && y<pict_get_ysize(p)-2 )   {
2652 >                                                x1=x-1;
2653 >                                                x2=x+1;
2654 >                                                y1=y-1;
2655 >                                                y2=y+1;
2656 >                                            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) ) {
2657 >                                            checkpixels = 0;
2658 >                                             actual_igs = before_igs;
2659 >                                             }  }
2660 >
2661 >
2662 >                                        if (pict_is_validpixel(p, x, y) && before_igs > 0 && checkpixels==1) {
2663                                                  actual_igs =
2664                                                          find_near_pgs(p, x, y, max_angle, before_igs,
2665                                                                                    igs, 1);
# Line 1911 | Line 2668 | if (calcfast == 1 ) {
2668                                                  }
2669                                                  if (before_igs > 0) {
2670                                                          actual_igs = pict_get_gsn(p, x, y);
2671 <                                                        icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2671 >                                                /*      icol = setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b); */
2672                                                  }
2673  
2674                                          }
2675                                  }
2676                          }
2677   /*      pict_write(p,"secondscan.pic");*/
2678 +        }}
2679  
1922        }
1923
2680   /*smoothing: add secondary glare sources */
2681          i_max = igs;
2682          if (sgs == 1) {
# Line 1957 | Line 2713 | if (calcfast == 1 ) {
2713  
2714   /* extract extremes from glare sources to extra glare source */
2715          if (splithigh == 1 && lum_total_max>limit) {
2716 + /*             fprintf(stderr,  " split of glare source!\n"); */
2717  
2718                  r_split = max_angle / 2.0;
2719                  for (i = 0; i <= i_max; i++) {
# Line 1977 | Line 2734 | if (calcfast == 1 ) {
2734                                                                          if (i_splitstart == (igs + 1)) {
2735                                                                                  igs = igs + 1;
2736                                                                                  pict_new_gli(p);
2737 +                                                                                pict_get_z1_gsn(p,igs) = 0;
2738 +                                                                                pict_get_z2_gsn(p,igs) = 0;
2739 +
2740                                                                                  pict_get_Eglare(p,igs) = 0.0;
2741                                                                                  pict_get_Dglare(p,igs) = 0.0;
2742                                                                                  pict_get_lum_min(p, igs) =
# Line 1990 | Line 2750 | if (calcfast == 1 ) {
2750                                                                          if (i_split == 0) {
2751                                                                                  igs = igs + 1;
2752                                                                                  pict_new_gli(p);
2753 +                                                                                pict_get_z1_gsn(p,igs) = 0;
2754 +                                                                                pict_get_z2_gsn(p,igs) = 0;
2755 +
2756                                                                                  pict_get_Eglare(p,igs) = 0.0;
2757                                                                                  pict_get_Dglare(p,igs) = 0.0;
2758                                                                                  pict_get_lum_min(p, igs) =
# Line 2026 | Line 2789 | if (calcfast == 1 ) {
2789                                                                                  if (before_igs > 0) {
2790                                                                                          actual_igs =
2791                                                                                                  pict_get_gsn(p, x, y);
2792 <                                                                                        icol =
2792 > /*                                                                                     icol =
2793                                                                                                  setglcolor(p, x, y,
2794 <                                                                                                                   actual_igs, uniform_gs, u_r, u_g , u_b);
2794 >                                                                                                                   actual_igs, uniform_gs, u_r, u_g , u_b);*/
2795                                                                                  }
2796  
2797                                                                          }
# Line 2043 | Line 2806 | if (calcfast == 1 ) {
2806                  }
2807          }
2808  
2809 < /* calculation of direct vertical illuminance for CGI and for disability glare*/
2809 > /* calculation of direct vertical illuminance for CGI and for disability glare, coloring glare sources*/
2810  
2811 <
2049 <        if (calcfast == 0) {
2811 >        if (calcfast == 0 || calcfast == 2) {
2812          for (x = 0; x < pict_get_xsize(p); x++)
2813                  for (y = 0; y < pict_get_ysize(p); y++) {
2814                          if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2815                                  if (pict_is_validpixel(p, x, y)) {
2816                                          if (pict_get_gsn(p, x, y) > 0) {
2817 <                                                pict_get_Eglare(p,pict_get_gsn(p, x, y)) +=
2818 <                                                        DOT(p->view.vdir, pict_get_cached_dir(p, x, y))
2819 <                                                        * pict_get_omega(p, x, y)
2820 <                                                        * luminance(pict_get_color(p, x, y));
2821 <                                                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));
2817 >                                                actual_igs = pict_get_gsn(p, x, y);
2818 >                                                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));
2819 >                                                pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
2820 >                                                E_v_dir = E_v_dir + delta_E;
2821 >                                                setglcolor(p, x, y, actual_igs, uniform_gs, u_r, u_g , u_b);
2822                                          }
2823                                  }
2824                          }
2825                  }
2826          lum_backg_cos = (E_v - E_v_dir) / 3.1415927;
2827 <        lum_backg = lum_backg_cos;
2827 >
2828           }
2829 < /* calc of band luminance if applied */
2829 > /* calc of band luminance distribution if applied */
2830          if (band_calc == 1) {
2831 <        band_avlum=get_band_lum( p, band_angle, band_color);
2832 <        }
2831 >                x_max = pict_get_xsize(p) - 1;
2832 >                y_max = pict_get_ysize(p) - 1;
2833 >                y_mid = (int)(y_max/2);
2834 >                for (x = 0; x <= x_max; x++)
2835 >                for (y = 0; y <= y_max; y++) {
2836 >                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2837 >                                if (pict_is_validpixel(p, x, y)) {
2838  
2839 +                        r_actual = acos(DOT(pict_get_cached_dir(p, x, y_mid), pict_get_cached_dir(p, x, y))) ;
2840 +                        act_lum = luminance(pict_get_color(p, x, y));
2841 +
2842 +                        if ((r_actual <= band_angle) || (y == y_mid) ) {
2843 +                                lum_band[0]=luminance(pict_get_color(p, x, y));
2844 +                                muc_rvar_add_sample(s_band, lum_band);
2845 +                                act_lum = luminance(pict_get_color(p, x, y));
2846 +                                lum_band_av += pict_get_omega(p, x, y) * act_lum;
2847 +                                omega_band += pict_get_omega(p, x, y);
2848 +                                if (band_color == 1) {
2849 +                                        pict_get_color(p, x, y)[RED] = 0.0;
2850 +                                        pict_get_color(p, x, y)[GRN] = act_lum / WHTEFFICACY / CIE_gf;
2851 +                                        pict_get_color(p, x, y)[BLU] = 0.0;
2852 +                                }
2853 +                        }
2854 +                }}}
2855 +                lum_band_av=lum_band_av/omega_band;
2856 +                muc_rvar_get_vx(s_band,lum_band_std);
2857 +                muc_rvar_get_percentile(s_band,lum_band_median,0.75);
2858 +                per_75_band=lum_band_median[0];
2859 +                muc_rvar_get_percentile(s_band,lum_band_median,0.95);
2860 +                per_95_band=lum_band_median[0];
2861 +                muc_rvar_get_median(s_band,lum_band_median);
2862 +                muc_rvar_get_bounding_box(s_band,bbox_band);
2863 +        
2864 +                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] );
2865 + /*      av_lum = av_lum / omega_sum; */
2866 + /*    printf("average luminance of band %f \n",av_lum);*/
2867 + /*      printf("solid angle of band %f \n",omega_sum);*/
2868 +        }
2869 +
2870   /*printf("total number of glare sources: %i\n",igs); */
2871          lum_sources = 0;
2872          omega_sources = 0;
2873 +        i=0;
2874          for (x = 0; x <= igs; x++) {
2875 +        if (pict_get_npix(p, x) > 0) {
2876 +                sum_glare += pow(pict_get_av_lum(p, x),2.0)* pict_get_av_omega(p, x)/pow(get_posindex(p, pict_get_av_posx(p, x), pict_get_av_posy(p, x),posindex_2), 2.0);
2877                  lum_sources += pict_get_av_lum(p, x) * pict_get_av_omega(p, x);
2878                  omega_sources += pict_get_av_omega(p, x);
2879 +                i=i+1;
2880 +        }}
2881 +        sum_glare= c2*log10(1 + sum_glare / pow(E_v, a3));
2882 +        if (sang == omega_sources) {
2883 +               lum_backg =avlum;
2884 +        } else {
2885 +               lum_backg =(sang * avlum - lum_sources) / (sang - omega_sources);
2886          }
2887 <        if (non_cos_lb == 1) {
2888 <                lum_backg =
2889 <                        (sang * avlum - lum_sources) / (sang - omega_sources);
2887 >
2888 >
2889 >        if (i == 0) {
2890 >               lum_sources =0.0;
2891 >        } else { lum_sources=lum_sources/omega_sources;
2892          }
2893 +
2894 +
2895 +
2896 +        if (non_cos_lb == 0) {
2897 +                        lum_backg = lum_backg_cos;
2898 +        }
2899 +
2900 +        if (non_cos_lb == 2) {
2901 +                        lum_backg = E_v / 3.1415927;
2902 +        }
2903 +
2904 +
2905 + /* file writing NOT here
2906 +        if (checkfile == 1) {
2907 +                pict_write(p, file_out);
2908 +        }
2909 + */
2910 +
2911   /* print detailed output */
2912 +        
2913 +        
2914 +
2915 + /* masking */
2916 +
2917 +        if ( masking == 1 ) {
2918 +        
2919 +        
2920 +                for (x = 0; x < pict_get_xsize(p); x++)
2921 +                for (y = 0; y < pict_get_ysize(p); y++) {
2922 +                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2923 +                                if (pict_is_validpixel(p, x, y)) {
2924 +                                        if (luminance(pict_get_color(pm, x, y))>0.1) {
2925 + /*                                         printf ("hallo %i %i %f \n",x,y,luminance(pict_get_color(pm, x, y)));*/
2926 +                                                  i_mask=i_mask+1;
2927 +                                                  lum_mask[0]=luminance(pict_get_color(p, x, y));
2928 +                                                 /* printf ("%f \n",lum_mask[0]);*/
2929 +                                                  muc_rvar_add_sample(s_mask, lum_mask);
2930 +                                                  omega_mask += pict_get_omega(p, x, y);
2931 +                                                  lum_mask_av += pict_get_omega(p, x, y)* luminance(pict_get_color(p, x, y));
2932 +                                                  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));
2933 +
2934 +                                                  pict_get_color(pm, x, y)[RED] = luminance(pict_get_color(p, x, y))/179.0;
2935 +                                                  pict_get_color(pm, x, y)[GRN] = luminance(pict_get_color(p, x, y))/179.0;
2936 +                                                  pict_get_color(pm, x, y)[BLU] = luminance(pict_get_color(p, x, y))/179.0;
2937 +  
2938 +                                           } else {
2939 +                                           pict_get_color(p, x, y)[RED] = 0.0 ;
2940 +                                           pict_get_color(p, x, y)[GRN] = 0.0 ;
2941 +                                           pict_get_color(p, x, y)[BLU] = 0.0 ;
2942 +                                           }
2943 +                                }
2944 +                        }
2945 +                }
2946 + /* Average luminance in masking area (weighted by solid angle) */
2947 +                lum_mask_av=lum_mask_av/omega_mask;
2948 +                muc_rvar_get_vx(s_mask,lum_mask_std);
2949 +                muc_rvar_get_percentile(s_mask,lum_mask_median,0.75);
2950 +                per_75_mask=lum_mask_median[0];
2951 +                muc_rvar_get_percentile(s_mask,lum_mask_median,0.95);
2952 +                per_95_mask=lum_mask_median[0];
2953 +                muc_rvar_get_median(s_mask,lum_mask_median);
2954 +                muc_rvar_get_bounding_box(s_mask,bbox);
2955 + /* PSGV only why masking of window is applied! */
2956 +
2957 +        
2958 +        if (task_lum == 0 || lum_task == 0.0 ) {
2959 +                        fprintf(stderr, " warning: Task area not set or task luminance=0 ! pgsv cannot be calculated (set to -99)!!\n");
2960 +                        pgsv = -99;
2961 +                } else {
2962 +                        pgsv = get_pgsv(E_v, E_v_mask, omega_mask, lum_mask_av,lum_task,avlum);
2963 +                        }
2964 +
2965 +                 pgsv_con = get_pgsv_con(E_v, E_v_mask, omega_mask, lum_mask_av,avlum);
2966 +                 pgsv_sat =get_pgsv_sat(E_v);
2967 +
2968          if (detail_out == 1) {
2969 +
2970 +                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 );
2971 +
2972 +        }      
2973 +                
2974 +        }
2975 +
2976 + /* zones */
2977 +
2978 +        if ( zones > 0 ) {
2979 +        
2980 +                for (x = 0; x < pict_get_xsize(p); x++)
2981 +                for (y = 0; y < pict_get_ysize(p); y++) {
2982 +                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
2983 +                                if (pict_is_validpixel(p, x, y)) {
2984 +
2985 +
2986 +                                 r_actual =acos(DOT(pict_get_cached_dir(p, x, y), pict_get_cached_dir(p, x_zone, y_zone))) * 2;
2987 +                                 lum_actual = luminance(pict_get_color(p, x, y));
2988 +                                 act_gsn=pict_get_gsn(p, x, y);
2989 +                            /*     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));*/
2990 +                                
2991 + /*zone1 */
2992 +                        if (r_actual <= angle_z1) {
2993 +                                                  i_z1=i_z1+1;
2994 +                                                  lum_z1[0]=lum_actual;
2995 +                                                  muc_rvar_add_sample(s_z1, lum_z1);
2996 +                                                  omega_z1 += pict_get_omega(p, x, y);
2997 +                                                  lum_z1_av += pict_get_omega(p, x, y)* lum_actual;
2998 +                                                  Ez1 += DOT(p->view.vdir, pict_get_cached_dir(p, x, y))* pict_get_omega(p, x, y)* lum_actual;
2999 +                                                  setglcolor(p,x,y,1,1 , 0.66, 0.01 ,0.33);
3000 + /*check if separation gsn already exist */
3001 +
3002 +                                                   if (act_gsn > 0){
3003 +                                                   if (pict_get_z1_gsn(p,act_gsn) == 0) {
3004 +                                                                                          pict_new_gli(p);
3005 +                                                                                          igs = igs + 1;
3006 +                                                                                          pict_get_z1_gsn(p,act_gsn) = igs*1.000;
3007 +                                                                                          pict_get_z1_gsn(p,igs) = -1.0;
3008 +                                                                                          pict_get_z2_gsn(p,igs) = -1.0;
3009 +                                                                                          splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
3010 + /*                                                                                        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)); */
3011 +                                                                                         }
3012 +                                                    splitgs=(int)(pict_get_z1_gsn(p,act_gsn));
3013 +                                        /*          printf("splitgs%i \n",splitgs);       */              
3014 +                                                    split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
3015 +                                        /* move direct illuminance contribution into  zone -value           */
3016 +                                                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));
3017 +                                                pict_get_Eglare(p,act_gsn ) = pict_get_Eglare(p,act_gsn ) - delta_E;
3018 +                                                pict_get_Eglare(p,igs ) = pict_get_Eglare(p,igs ) + delta_E;
3019 +            
3020 +                                                    
3021 +                                                }                                
3022 +                                                  }
3023 + /*zone2 */
3024 +
3025 +                        if (r_actual > angle_z1 && r_actual<= angle_z2 ) {
3026 +
3027 +                                                  i_z2=i_z2+1;
3028 +                                                  lum_z2[0]=lum_actual;
3029 +                                                  muc_rvar_add_sample(s_z2, lum_z2);
3030 +                                                  omega_z2 +=   pict_get_omega(p, x, y);
3031 +                                                  lum_z2_av += pict_get_omega(p, x, y)* lum_actual;
3032 +                                                  Ez2 += DOT(p->view.vdir, pict_get_cached_dir(p, x, y))* pict_get_omega(p, x, y)* lum_actual;
3033 +                                                  setglcolor(p,x,y,1,1 , 0.65, 0.33 ,0.02);
3034 + /*                                                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));
3035 + */                                                if (act_gsn > 0){
3036 +                                                   if (pict_get_z2_gsn(p,act_gsn) == 0) {
3037 +                                                                                          pict_new_gli(p);
3038 +                                                                                          igs = igs + 1;
3039 +                                                                                          pict_get_z2_gsn(p,act_gsn) = igs*1.000;
3040 +                                                                                          pict_get_z1_gsn(p,igs) = -2.0;
3041 +                                                                                          pict_get_z2_gsn(p,igs) = -2.0;
3042 + /*                                                                                         printf ("Glare source %i gets glare source %i in zone 2 \n",act_gsn,igs ); */
3043 +                                                                                         }
3044 +                                                splitgs=(int)(pict_get_z2_gsn(p,act_gsn));
3045 + /*                                              printf("splitgs %i \n",splitgs);*/                              
3046 +                                                    split_pixel_from_gs(p, x, y, splitgs, uniform_gs, u_r, u_g , u_b);
3047 +                                        /* move direct illuminance contribution into  zone -value           */
3048 +                                                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));
3049 +                                                pict_get_Eglare(p,act_gsn ) = pict_get_Eglare(p,act_gsn ) - delta_E;
3050 +                                                pict_get_Eglare(p,igs ) = pict_get_Eglare(p,igs ) + delta_E;
3051 +
3052 +                                           }
3053 +                                }
3054 +
3055 +                        }
3056 +                } }
3057 + /* Average luminance in zones (weighted by solid angle), calculate percentiles, median min and max in zones  */
3058 +               if (zones == 2 ) {
3059 +                lum_z2_av=lum_z2_av/omega_z2;
3060 +                muc_rvar_get_vx(s_z2,lum_z2_std);
3061 +                muc_rvar_get_percentile(s_z2,lum_z2_median,0.75);
3062 +                per_75_z2=lum_z2_median[0];
3063 +                muc_rvar_get_percentile(s_z2,lum_z2_median,0.95);
3064 +                per_95_z2=lum_z2_median[0];
3065 +                muc_rvar_get_median(s_z2,lum_z2_median);
3066 +                muc_rvar_get_bounding_box(s_z2,bbox_z2);
3067 +                }
3068 +                lum_z1_av=lum_z1_av/omega_z1;
3069 +                muc_rvar_get_vx(s_z1,lum_z1_std);
3070 +                muc_rvar_get_percentile(s_z1,lum_z1_median,0.75);
3071 +                per_75_z1=lum_z1_median[0];
3072 +                muc_rvar_get_percentile(s_z1,lum_z1_median,0.95);
3073 +                per_95_z1=lum_z1_median[0];
3074 +                muc_rvar_get_median(s_z1,lum_z1_median);
3075 +                muc_rvar_get_bounding_box(s_z1,bbox_z1);
3076 +        if (detail_out == 1) {
3077 +
3078 +                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,Ez1,i_z1: %f %f %f %f %f %f %f %f %f %i\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],Ez1,i_z1 );
3079 +
3080 +               if (zones == 2 ) {
3081 +
3082 +                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,Ez2,i_z2:  %f %f %f %f %f %f %f %f %f %i\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],Ez2,i_z1 );
3083 + } }            
3084 +                
3085 +        }
3086 +
3087 + new_gs_number = (int *)malloc(sizeof(int)*(igs+1));
3088 +
3089                  i = 0;
3090                  for (x = 0; x <= igs; x++) {
3091   /* resorting glare source numbers */
3092                          if (pict_get_npix(p, x) > 0) {
3093                                  i = i + 1;
3094 +                                pict_get_Dglare(p,x) = i;
3095 +                                new_gs_number[x] = i;
3096 + /*                      printf("%i %i %f %i \n", i,x,pict_get_Dglare(p,x),new_gs_number[x] ); */
3097                          }
3098                  }
3099 +                no_glaresources=i;
3100 + /*printf("%i",no_glaresources );*/
3101 + /* glare sources */
3102  
3103 +        if (output == 4 ) {
3104 +
3105 +        i=0;
3106 +        for (x = 0; x <= igs; x++) {
3107 +                                if (pict_get_npix(p, x) > 0) {
3108 +                                        i = i + 1;
3109  
3110 +                                        x2=pict_get_av_posx(p, x);
3111 +                                        y2=pict_get_av_posy(p, x);
3112  
3113 +                                        printf("%i %f %f %f %f %.10f %f %f %f %f \n",
3114 +                                                   i, pict_get_npix(p, x), pict_get_av_posx(p, x), pict_get_av_posy(p, x),
3115 +                                                   pict_get_av_lum(p, x), pict_get_av_omega(p, x),
3116 +                                                   get_posindex(p, pict_get_av_posx(p, x),pict_get_av_posy(p, x),posindex_2),
3117 +                                                   pict_get_cached_dir(p, x2,y2)[0],pict_get_cached_dir(p, x2,y2)[1],pict_get_cached_dir(p, x2,y2)[2]);
3118 +        }
3119 +        }
3120 +
3121 +
3122 +                for (y = 0; y < pict_get_ysize(p); y++) {
3123 +                        for (x = 0; x < pict_get_xsize(p); x++) {
3124 +
3125 +                        if (pict_get_hangle(p, x, y, p->view.vdir, p->view.vup, &ang)) {
3126 +                                if (pict_is_validpixel(p, x, y)) {
3127 +                        if (pict_get_gsn(p,x,y) >0 ) {
3128 +                        i=pict_get_gsn(p,x,y);
3129 +                        printf("%i %i %f %f %f \n", i, new_gs_number[i], pict_get_cached_dir(p, x,y)[0],pict_get_cached_dir(p, x,y)[1],pict_get_cached_dir(p, x,y)[2] );
3130 +                        
3131 +                        }
3132 +
3133 +
3134 +
3135 + }}}}
3136 +
3137 +
3138 + }
3139 +
3140 +
3141 +
3142 +
3143 +        if (detail_out == 1 && output < 3) {
3144 +        dgp = get_dgp(p, E_v, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
3145 +
3146                  printf
3147 <                        ("%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",
3147 >                        ("%i No pixels x-pos y-pos L_s Omega_s Posindx L_b L_t E_v Edir Max_Lum Sigma xdir ydir zdir Eglare Lveil_cie teta glare_zone r_contrast r_dgp\n",
3148                           i);
3149                  if (i == 0) {
3150 <                        printf("%i %f %f %f %f %.10f %f %f %f %f %f %f\n", i, 0.0, 0.0,
3151 <                                   0.0, 0.0, 0.0, 0.0, lum_backg, lum_task, E_v, E_v_dir,
2104 <                                   abs_max);
3150 >                        printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n", i, 0.0, 0.0,
3151 >                                   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,0.0,0.0);
3152  
3153                  } else {
3154                          i = 0;
# Line 2121 | Line 3168 | if (calcfast == 1 ) {
3168                                                                       Lveil_cie =0 ;
3169                                                                     }
3170                                          Lveil_cie_sum =  Lveil_cie_sum + Lveil_cie;
3171 <                                        printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
3171 >                                        r_glare= c2*log10(1 + pow(pict_get_av_lum(p, x),2)* pict_get_av_omega(p, x)/pow(get_posindex(p, pict_get_av_posx(p, x), pict_get_av_posy(p, x),posindex_2), 2.0) / pow(E_v, a3))  ;
3172 >                                        r_contrast=r_glare/sum_glare  ;
3173 >                                        r_dgp=r_glare/dgp ;
3174 >                                        if (pict_get_z1_gsn(p,x)<0) {
3175 >                                        act_gsn=(int)(-pict_get_z1_gsn(p,x));
3176 >                                        }else{
3177 >                                        act_gsn=0;
3178 >                                        }
3179 >                                        printf("%i %f %f %f %f %.10f %f %f %f %f %f %f %f %f %f %f %f %f %f %i %f %f \n",
3180                                                     i, pict_get_npix(p, x), pict_get_av_posx(p, x),
3181                                                     pict_get_ysize(p) - pict_get_av_posy(p, x),
3182                                                     pict_get_av_lum(p, x), pict_get_av_omega(p, x),
# Line 2129 | Line 3184 | if (calcfast == 1 ) {
3184                                                                                  pict_get_av_posy(p, x),
3185                                                                                  posindex_2), lum_backg, lum_task,
3186                                                     E_v, E_v_dir, abs_max, sigma * 180 / 3.1415927
3187 <                                                   ,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 <                                                   );
3187 >                                                   ,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,r_contrast,r_dgp );
3188                                  }
3189                          }
3190                  }
3191          }
3192  
3193 + if ( output < 3) {
3194  
2141
3195   /* calculation of indicees */
3196  
3197   /* check vertical illuminance range */
# Line 2151 | Line 3204 | if (calcfast == 1 ) {
3204          dgp =
3205                  get_dgp(p, E_v2, igs, a1, a2, a3, a4, a5, c1, c2, c3, posindex_2);
3206   /* low light correction */
3207 +     if (lowlight ==1) {
3208         if (E_v < 1000) {
3209         low_light_corr=1.0*exp(0.024*E_v-4)/(1+exp(0.024*E_v-4));} else {low_light_corr=1.0 ;}
3210         dgp =low_light_corr*dgp;
3211 <      
3211 >       }
3212   /* age correction */
3213        
3214          if (age_corr == 1) {
# Line 2171 | Line 3225 | if (calcfast == 1 ) {
3225          }
3226  
3227   if (calcfast == 0) {
3228 +        lum_a= E_v2/3.1415927;
3229          dgi = get_dgi(p, lum_backg, igs, posindex_2);
3230          ugr = get_ugr(p, lum_backg, igs, posindex_2);
3231 +        ugp = get_ugp (p,lum_backg, igs, posindex_2);
3232 +        ugp2 = get_ugp2 (p,lum_backg, igs, posindex_2);
3233 +        ugr_exp = get_ugr_exp (p,lum_backg_cos,lum_a, igs, posindex_2);
3234 +        dgi_mod = get_dgi_mod(p, lum_a, igs, posindex_2);
3235          cgi = get_cgi(p, E_v, E_v_dir, igs, posindex_2);
3236 <        vcp = get_vcp(p, avlum, igs, posindex_2);
3236 >        dgr = get_dgr(p, avlum, igs, posindex_2);
3237 >        vcp = get_vcp(dgr);
3238          Lveil = get_disability(p, avlum, igs);
3239 +        if (no_glaresources == 0) {
3240 +                dgi = 0.0;
3241 +                ugr = 0.0;
3242 +                ugp = 0.0;
3243 +                ugp2 =0.0;
3244 +                ugr_exp = 0.0;
3245 +                dgi_mod = 0.0;
3246 +                cgi = 0.0;
3247 +                dgr = 0.0;
3248 +                vcp = 100.0;
3249 +        }
3250   }
3251   /* check dgp range */
3252          if (dgp <= 0.2) {
# Line 2196 | Line 3267 | if (calcfast == 0) {
3267                                 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
3268                                  posindex_2);
3269                                  dgp = dgp_ext;
3270 <                                if (E_vl_ext < 1000) {
3270 >                                if (E_vl_ext < 1000 && lowlight ==1) {
3271                                  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 ;}
3272                                  dgp =low_light_corr*dgp;
3273                                  dgp =age_corr_factor*dgp;
3274                          }
3275 +                muc_rvar_get_median(s_noposweight,lum_nopos_median);
3276 +                muc_rvar_get_median(s_posweight,lum_pos_median);
3277 +                muc_rvar_get_median(s_posweight2,lum_pos2_median);
3278 +
3279  
3280 +
3281 +        
3282                          printf
3283 <                                ("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",
3283 >                                ("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,ugp2: %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
3284                                   dgp, avlum, E_v, lum_backg, E_v_dir, dgi, ugr, vcp, cgi,
3285 <                                 lum_sources / omega_sources, omega_sources, Lveil,Lveil_cie_sum,band_avlum);
3285 >                                 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],ugp2);
3286                  } else {
3287                          if (detail_out2 == 1) {
3288  
# Line 2218 | Line 3295 | if (calcfast == 0) {
3295                                  if (ext_vill == 1) {
3296                  dgp_ext = get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,posindex_2);
3297                                  
3298 <                                if (E_vl_ext < 1000) {
3298 >                                if (E_vl_ext < 1000 && lowlight ==1) {
3299                                  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 ;}
3300                                          dgp =low_light_corr*dgp_ext;
3301                                          dgp =age_corr_factor*dgp;
# Line 2233 | Line 3310 | if (calcfast == 0) {
3310                                 get_dgp(p, E_vl_ext, igs, a1, a2, a3, a4, a5, c1, c2, c3,
3311                                  posindex_2);
3312                                  dgp = dgp_ext;
3313 <                                if (E_vl_ext < 1000) {
3313 >                                if (E_vl_ext < 1000 && lowlight ==1) {
3314                                  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 ;}
3315                                  dgp =low_light_corr*dgp;
3316 <                                dgp =age_corr_factor*dgp;
3317 <                printf("%f\n", dgp);
3316 >
3317 >                     if (calcfast == 2) {
3318 >                    
3319 >                         lum_backg_cos=(E_vl_ext-dir_ill)/3.1415927;
3320 >                         ugr = get_ugr(p, lum_backg_cos, igs, posindex_2);
3321 >                         printf("%f %f \n", dgp,ugr);
3322 >                     }else{      
3323 >                         printf("%f\n", dgp);
3324 >                }
3325          }
3326 + }
3327  
3328 + }else{
3329 + /* only multiimagemode                        
3330  
3331 + */
3332 +
3333 +
3334 +                       for (j = 0; j < num_images; j++) {
3335 +
3336 +                                
3337 + /* loop over temporal images */
3338 +
3339 + pict_read(pcoeff,temp_image_name[j] );
3340 +
3341 + /*printf ("num_img:%i x-size:%i xpos:%i y-size: %i ypos %i bigimg-xsize %i %i ",num_images,pict_get_xsize(pcoeff),x_temp_img[j],pict_get_ysize(pcoeff),y_temp_img[j],pict_get_xsize(p),pict_get_ysize(p));
3342 + */
3343 +
3344 +
3345 + /* loop over scans
3346 +    empty of */
3347 + for (jj=0;jj<=num_scans;jj++){
3348 +
3349 + /*printf("scale %i %i %i %f ",j,jj,num_scans,scale_image_scans[j*(num_scans+1)+jj]); */
3350 +
3351 +
3352 + /* copy luminance value into big image and remove glare sources*/
3353 +        for (x = 0; x < pict_get_xsize(pcoeff); x++) {
3354 +                for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3355 +                        xmap=x_temp_img[j]+x;
3356 +                        ymap=y_temp_img[j]+y;
3357 +                        if (xmap <0) { xmap=0;}
3358 +                        if (ymap <0) { ymap=0;}
3359 +                        
3360 +                        pict_get_color(p, xmap, ymap)[RED] = scale_image_scans[j*(num_scans+1)+jj]*pict_get_color(pcoeff, x, y)[RED];
3361 +                        pict_get_color(p, xmap, ymap)[GRN] = scale_image_scans[j*(num_scans+1)+jj]*pict_get_color(pcoeff, x, y)[GRN];
3362 +                        pict_get_color(p, xmap, ymap)[BLU] = scale_image_scans[j*(num_scans+1)+jj]*pict_get_color(pcoeff, x, y)[BLU];
3363 +                        pict_get_gsn(p, xmap, ymap) = 0;
3364 +                        pict_get_pgs(p, xmap, ymap) = 0;
3365 + }}
3366 +
3367 +
3368 +
3369 +
3370 +
3371 + actual_igs =0;
3372 +
3373 + /* first glare source scan: find primary glare sources */
3374 +        for (x = 0; x < pict_get_xsize(pcoeff); x++) {
3375 +                lastpixelwas_gs=0;
3376 + /*              lastpixelwas_peak=0; */
3377 +                for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3378 +                        xmap=x_temp_img[j]+x;
3379 +                        ymap=y_temp_img[j]+y;
3380 +                        if (pict_get_hangle(p, xmap, ymap, p->view.vdir, p->view.vup, &ang)) {
3381 +                                if (pict_is_validpixel(p, xmap, ymap)) {
3382 +                                        act_lum = luminance(pict_get_color(p, xmap, ymap));
3383 +                                        if (act_lum> lum_source) {
3384 +                                                if (act_lum >lum_total_max) {
3385 +                                                                             lum_total_max=act_lum;
3386 +                                                                               }
3387 +                      
3388 +                                                if (lastpixelwas_gs==0 || search_pix <= 1.0 ) {
3389 +                                                actual_igs = find_near_pgs(p, xmap, ymap, max_angle, 0, igs, 1);
3390 +  }
3391 +                                                if (actual_igs == 0) {
3392 +                                                        igs = igs + 1;
3393 +                                                        pict_new_gli(p);
3394 +                                                        pict_get_Eglare(p,igs) = 0.0;
3395 + /*  not necessary here                                  pict_get_lum_min(p, igs) = HUGE_VAL;
3396 +                                                        pict_get_Eglare(p,igs) = 0.0;
3397 +                                                        pict_get_Dglare(p,igs) = 0.0;
3398 +                                                        pict_get_z1_gsn(p,igs) = 0;
3399 +                                                        pict_get_z2_gsn(p,igs) = 0; */
3400 +                                                        actual_igs = igs;
3401 +                                                        
3402 +                                                }
3403 +                                                pict_get_gsn(p, xmap, ymap) = actual_igs;
3404 +                                                pict_get_pgs(p, xmap, ymap) = 1;
3405 +                                                add_pixel_to_gs(p, xmap, ymap, actual_igs);
3406 +                                                lastpixelwas_gs=actual_igs;
3407 +                                                
3408 +                                                
3409 +                                                
3410 +                                                delta_E=DOT(p->view.vdir, pict_get_cached_dir(p, xmap, ymap))* pict_get_omega(p, xmap, ymap)* luminance(pict_get_color(p, xmap, ymap));
3411 +                                                pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
3412 +
3413 +                                              
3414 +                                                
3415 +
3416 +                                        } else {
3417 +                                                pict_get_pgs(p, xmap, ymap) = 0;
3418 +                                                lastpixelwas_gs=0;
3419 +                                        }
3420 +                                }
3421 +                        }
3422 +                }
3423 +             }
3424 +
3425 +
3426 + /*                              here should be peak extraction  */
3427 + i_max=igs;
3428 +                r_split = max_angle / 2.0;
3429 +                for (i = 0; i <= i_max; i++) {
3430 +
3431 +                        if (pict_get_npix(p, i) > 0) {
3432 +                                l_max = pict_get_lum_max(p, i);
3433 +                                i_splitstart = igs + 1;
3434 +                                if (l_max >= limit) {
3435 +                                        for (x = 0; x < pict_get_xsize(pcoeff); x++)
3436 +                                                for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3437 +                                                xmap=x_temp_img[j]+x;
3438 +                                                ymap=y_temp_img[j]+y;
3439 +                                                
3440 +                                                
3441 +                                                        if (pict_get_hangle
3442 +                                                                (p, xmap, ymap, p->view.vdir, p->view.vup, &ang))
3443 +                                                        {
3444 +                                                                if (pict_is_validpixel(p, xmap, ymap)
3445 +                                                                        && luminance(pict_get_color(p, xmap, ymap))
3446 +                                                                        >= limit
3447 +                                                                        && pict_get_gsn(p, xmap, ymap) == i) {
3448 +                                                                        if (i_splitstart == (igs + 1)) {
3449 +                                                                                igs = igs + 1;
3450 +                                                                                pict_new_gli(p);
3451 +                                                                                pict_get_z1_gsn(p,igs) = 0;
3452 +                                                                                pict_get_z2_gsn(p,igs) = 0;
3453 +
3454 +                                                                                pict_get_Eglare(p,igs) = 0.0;
3455 +                                                                                pict_get_Dglare(p,igs) = 0.0;
3456 +                                                                                pict_get_lum_min(p, igs) =
3457 +                                                                                        99999999999999.999;
3458 +                                                                                i_split = igs;
3459 +                                                                        } else {
3460 +                                                                                i_split =
3461 +                                                                                        find_split(p, xmap, ymap, r_split,
3462 +                                                                                                           i_splitstart, igs);
3463 +                                                                        }
3464 +                                                                        if (i_split == 0) {
3465 +                                                                                igs = igs + 1;
3466 +                                                                                pict_new_gli(p);
3467 +                                                                                pict_get_z1_gsn(p,igs) = 0;
3468 +                                                                                pict_get_z2_gsn(p,igs) = 0;
3469 +
3470 +                                                                                pict_get_Eglare(p,igs) = 0.0;
3471 +                                                                                pict_get_Dglare(p,igs) = 0.0;
3472 +                                                                                pict_get_lum_min(p, igs) =
3473 +                                                                                        99999999999999.999;
3474 +                                                                                i_split = igs;
3475 +                                                                        }
3476 +                                                                        split_pixel_from_gs(p, xmap, ymap, i_split, uniform_gs, u_r, u_g , u_b);
3477 +
3478 +                                                                }
3479 +                                                        }
3480 +                                                }
3481 +
3482 +                                }
3483 +                                change = 1;
3484 +                                while (change == 1) {
3485 +                                        change = 0;
3486 +                                        for (x = 0; x < pict_get_xsize(pcoeff); x++)
3487 +                                                for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3488 +                                                xmap=x_temp_img[j]+x;
3489 +                                                ymap=y_temp_img[j]+y;
3490 +                                                        before_igs = pict_get_gsn(p, xmap, ymap);
3491 +                                                        if (before_igs >= i_splitstart) {
3492 +                                                                if (pict_get_hangle
3493 +                                                                        (p, xmap, ymap, p->view.vdir, p->view.vup,
3494 +                                                                         &ang)) {
3495 +                                                                        if (pict_is_validpixel(p, xmap, ymap)
3496 +                                                                                && before_igs > 0) {
3497 +                                                                                actual_igs =
3498 +                                                                                        find_near_pgs(p, xmap, ymap,
3499 +                                                                                                                  max_angle,
3500 +                                                                                                                  before_igs, igs,
3501 +                                                                                                                  i_splitstart);
3502 +                                                                                if (!(actual_igs == before_igs)) {
3503 +                                                                                        change = 1;
3504 +                                                                                }
3505 +                                                                                if (before_igs > 0) {
3506 +                                                                                        actual_igs =
3507 +                                                                                                pict_get_gsn(p, xmap, ymap);
3508 + /*                                                                                     icol =
3509 +                                                                                                setglcolor(p, x, y,
3510 +                                                                                                                   actual_igs, uniform_gs, u_r, u_g , u_b);*/
3511 +                                                                                }
3512 +
3513 +                                                                        }
3514 +                                                                }
3515 +                                                        }
3516 +
3517 +                                                }
3518 +                                }
3519 +
3520 +
3521 +                        }
3522 +                }
3523 +
3524 + /*                              end peak extraction  */
3525 +
3526 +
3527 + /* calculation of direct vertical illuminance for th multi-image-mode */
3528 +
3529 +
3530 +        for (x = 0; x < pict_get_xsize(pcoeff); x++)
3531 +                for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3532 +                        xmap=x_temp_img[j]+x;
3533 +                        ymap=y_temp_img[j]+y;
3534 +                        if (pict_get_hangle(p, xmap, ymap, p->view.vdir, p->view.vup, &ang)) {
3535 +                                if (pict_is_validpixel(p, xmap, ymap)) {
3536 +                                        if (pict_get_gsn(p, xmap, ymap) > 0) {
3537 +                                                actual_igs = pict_get_gsn(p, xmap, ymap);
3538 +                                                delta_E=DOT(p->view.vdir, pict_get_cached_dir(p, xmap, ymap))* pict_get_omega(p, xmap, ymap)* luminance(pict_get_color(p, xmap, ymap));
3539 +                                                pict_get_Eglare(p,actual_igs ) = pict_get_Eglare(p,actual_igs ) + delta_E;
3540 +                                        }
3541 +                                }
3542 +                        }
3543 +                }
3544 +
3545 +
3546 +
3547 +
3548 +
3549 +
3550 +
3551 +
3552 +                        i = 0;
3553 +                        for (x = 0; x <= igs; x++) {
3554 +                                if (pict_get_npix(p, x) > 0) {
3555 +                                        i = i + 1;
3556 +                                        }}
3557 + no_glaresources=i;                      
3558 +
3559 + /*
3560 +
3561 + sprintf(file_out, "%s%i%s","ray2/img_",j,".hdr");
3562 + pict_write(p, file_out);
3563 + */
3564 + printf("%i ",no_glaresources);
3565 +                        i = 0;
3566 +                        for (x = 0; x <= igs; x++) {
3567 +                                if (pict_get_npix(p, x) > 0) {
3568 +                                        i = i + 1;
3569 +                                        pict_get_sigma(p, pict_get_av_posx(p, x), pict_get_av_posy(p, x), p->view.vdir, p->view.vup, &sigma);
3570 +                                                                  
3571 +                                        x2=pict_get_av_posx(p, x);
3572 +                                        y2=pict_get_av_posy(p, x);
3573 +                                        printf("%f %.10f %f %f %f %f %f ",      pict_get_av_lum(p, x), pict_get_av_omega(p, x), get_posindex(p, pict_get_av_posx(p, x), pict_get_av_posy(p, x),
3574 +                                                                posindex_2),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) );
3575 +                                }
3576 + pict_get_npix(p, x)=0;
3577 + pict_get_av_lum(p, x)=0;
3578 + pict_get_av_posy(p, x)=0;
3579 + pict_get_av_posx(p, x)=0;
3580 + pict_get_av_omega(p, x)=0;
3581 +                        }
3582 + printf("\n");
3583 +
3584 +
3585 + /* empty big image and remove glare sources*/
3586 +        for (x = 0; x < pict_get_xsize(pcoeff); x++) {
3587 +                for (y = 0; y < pict_get_ysize(pcoeff); y++) {
3588 +                        xmap=x_temp_img[j]+x;
3589 +                        ymap=y_temp_img[j]+y;
3590 +                        pict_get_color(p, xmap, ymap)[RED] = 0;
3591 +                        pict_get_color(p, xmap, ymap)[GRN] = 0;
3592 +                        pict_get_color(p, xmap, ymap)[BLU] = 0;
3593 +                        pict_get_gsn(p, xmap, ymap) = 0;
3594 +                        pict_get_pgs(p, xmap, ymap) = 0;
3595 + }}
3596 + igs=0;
3597 +
3598 +
3599 + }
3600 +
3601 +
3602 + }
3603 +
3604 + }
3605 +
3606 + /* end multi-image-mode */
3607 +
3608   /*printf("lowlight=%f\n", low_light_corr); */
3609  
3610  
# Line 2268 | Line 3632 | has to be re-written from scratch....
3632  
3633   /*output  */
3634          if (checkfile == 1) {
3635 +                        
3636 +        
3637                  pict_write(p, file_out);
3638          }
3639  

Diff Legend

Removed lines
+ Added lines
< Changed lines (old)
> Changed lines (new)