Back to index

radiance  4R0+20100331
findglare.c
Go to the documentation of this file.
00001 #ifndef lint
00002 static const char    RCSid[] = "$Id: findglare.c,v 2.13 2009/02/07 05:40:47 greg Exp $";
00003 #endif
00004 /*
00005  * Find glare sources in a scene or image.
00006  *
00007  *     Greg Ward     March 1991
00008  */
00009 
00010 #include "glare.h"
00011 
00012 #define FEQ(a,b)     ((a)-(b)<=FTINY&&(b)-(a)<=FTINY)
00013 #define VEQ(v1,v2)   (FEQ((v1)[0],(v2)[0])&&FEQ((v1)[1],(v2)[1]) \
00014                             &&FEQ((v1)[2],(v2)[2]))
00015 
00016 char   *rtargv[64] = {"rtrace", "-h-", "-ov", "-fff", "-ld-", "-i-", "-I-"};
00017 int    rtargc = 7;
00018 
00019 VIEW   ourview = STDVIEW;          /* our view */
00020 VIEW   pictview = STDVIEW;         /* picture view */
00021 VIEW   leftview, rightview;        /* leftmost and rightmost views */
00022 
00023 char   *picture = NULL;            /* picture file name */
00024 char   *octree = NULL;                    /* octree file name */
00025 
00026 int    verbose = 0;                /* verbose reporting */
00027 char   *progname;                  /* global argv[0] */
00028 
00029 double threshold = 0.;                    /* glare threshold */
00030 
00031 int    sampdens = SAMPDENS;        /* sample density */
00032 ANGLE  glarang[180] = {AEND};             /* glare calculation angles */
00033 int    nglarangs = 0;
00034 double maxtheta;                   /* maximum angle (in radians) */
00035 int    hsize;                      /* horizontal size */
00036 
00037 struct illum  *indirect;           /* array of indirect illuminances */
00038 
00039 long   npixinvw;                   /* number of pixels in view */
00040 long   npixmiss;                   /* number of pixels missed */
00041 
00042 static int angcmp(const void *ap1, const void *ap2);
00043 static void init(void);
00044 static void cleanup(void);
00045 static void printsources(void);
00046 static void printillum(void);
00047 
00048 
00049 
00050 int
00051 main(
00052        int    argc,
00053        char   *argv[]
00054 )
00055 {
00056        int    combine = 1;
00057        int    gotview = 0;
00058        int    rval, i;
00059        char   *err;
00060 
00061        progname = argv[0];
00062                                    /* process options */
00063        for (i = 1; i < argc && argv[i][0] == '-'; i++) {
00064                                           /* expand arguments */
00065               while ((rval = expandarg(&argc, &argv, i)) > 0)
00066                      ;
00067               if (rval < 0) {
00068                      fprintf(stderr, "%s: cannot expand '%s'\n",
00069                                    argv[0], argv[i]);
00070                      exit(1);
00071               }
00072               rval = getviewopt(&ourview, argc-i, argv+i);
00073               if (rval >= 0) {
00074                      i += rval;
00075                      gotview++;
00076                      continue;
00077               }
00078               switch (argv[i][1]) {
00079               case 't':
00080                      threshold = atof(argv[++i]);
00081                      break;
00082               case 'r':
00083                      sampdens = atoi(argv[++i])/2;
00084                      break;
00085               case 'v':
00086                      if (argv[i][2] == '\0') {
00087                             verbose++;
00088                             break;
00089                      }
00090                      if (argv[i][2] != 'f')
00091                             goto userr;
00092                      rval = viewfile(argv[++i], &ourview, NULL);
00093                      if (rval < 0) {
00094                             fprintf(stderr,
00095                             "%s: cannot open view file \"%s\"\n",
00096                                           progname, argv[i]);
00097                             exit(1);
00098                      } else if (rval == 0) {
00099                             fprintf(stderr,
00100                                    "%s: bad view file \"%s\"\n",
00101                                           progname, argv[i]);
00102                             exit(1);
00103                      } else
00104                             gotview++;
00105                      break;
00106               case 'g':
00107                      if (argv[i][2] != 'a')
00108                             goto userr;
00109                      if (setscan(glarang, argv[++i]) < 0) {
00110                             fprintf(stderr, "%s: bad angle spec \"%s\"\n",
00111                                    progname, argv[i]);
00112                             exit(1);
00113                      }
00114                      break;
00115               case 'p':
00116                      picture = argv[++i];
00117                      break;
00118               case 'c':
00119                      combine = !combine;
00120                      break;
00121               case 'd':
00122                      rtargv[rtargc++] = argv[i];
00123                      if (argv[i][2] != 'v')
00124                             rtargv[rtargc++] = argv[++i];
00125                      break;
00126               case 'l':
00127                      if (argv[i][2] == 'd')
00128                             break;
00129                      /* FALL THROUGH */
00130               case 's':
00131               case 'P':
00132                      rtargv[rtargc++] = argv[i];
00133                      rtargv[rtargc++] = argv[++i];
00134                      break;
00135               case 'w':
00136                      rtargv[rtargc++] = argv[i];
00137                      break;
00138               case 'a':
00139                      rtargv[rtargc++] = argv[i];
00140                      if (argv[i][2] == 'v') {
00141                             rtargv[rtargc++] = argv[++i];
00142                             rtargv[rtargc++] = argv[++i];
00143                      }
00144                      rtargv[rtargc++] = argv[++i];
00145                      break;
00146               default:
00147                      goto userr;
00148               }
00149        }
00150                                           /* get octree */
00151        if (i < argc-1)
00152               goto userr;
00153        if (i == argc-1) {
00154               rtargv[rtargc++] = octree = argv[i];
00155               rtargv[rtargc] = NULL;
00156        }
00157                                           /* get view */
00158        if (picture != NULL) {
00159               rval = viewfile(picture, &pictview, NULL);
00160               if (rval < 0) {
00161                      fprintf(stderr, "%s: cannot open picture file \"%s\"\n",
00162                                    progname, picture);
00163                      exit(1);
00164               } else if (rval == 0) {
00165                      fprintf(stderr,
00166                             "%s: cannot get view from picture \"%s\"\n",
00167                                    progname, picture);
00168                      exit(1);
00169               }
00170               if (pictview.type == VT_PAR) {
00171                      fprintf(stderr, "%s: %s: cannot use parallel view\n",
00172                                    progname, picture);
00173                      exit(1);
00174               }
00175               if ((err = setview(&pictview)) != NULL) {
00176                      fprintf(stderr, "%s: %s\n", picture, err);
00177                      exit(1);
00178               }
00179        }
00180        if (!gotview) {
00181               if (picture == NULL) {
00182                      fprintf(stderr, "%s: must have view or picture\n",
00183                                    progname);
00184                      exit(1);
00185               }
00186               ourview = pictview;
00187        } else if (picture != NULL && !VEQ(ourview.vp, pictview.vp)) {
00188               fprintf(stderr, "%s: picture must have same viewpoint\n",
00189                             progname);
00190               exit(1);
00191        }
00192        ourview.type = VT_HEM;
00193        ourview.horiz = ourview.vert = 180.0;
00194        ourview.hoff = ourview.voff = 0.0;
00195        fvsum(ourview.vdir, ourview.vdir, ourview.vup,
00196                      -DOT(ourview.vdir,ourview.vup));
00197        if ((err = setview(&ourview)) != NULL) {
00198               fprintf(stderr, "%s: %s\n", progname, err);
00199               exit(1);
00200        }
00201        if (octree == NULL && picture == NULL) {
00202               fprintf(stderr,
00203               "%s: must specify at least one of picture or octree\n",
00204                             progname);
00205               exit(1);
00206        }
00207        init();                                   /* initialize program */
00208        if (threshold <= FTINY)
00209               comp_thresh();                     /* compute glare threshold */
00210        analyze();                         /* analyze view */
00211        if (combine)
00212               absorb_specks();            /* eliminate tiny sources */
00213        cleanup();                         /* tidy up */
00214                                           /* print header */
00215        newheader("RADIANCE", stdout);
00216        printargs(argc, argv, stdout);
00217        fputs(VIEWSTR, stdout);
00218        fprintview(&ourview, stdout);
00219        printf("\n");
00220        fputformat("ascii", stdout);
00221        printf("\n");
00222        printsources();                           /* print glare sources */
00223        printillum();                      /* print illuminances */
00224        exit(0);
00225 userr:
00226        fprintf(stderr,
00227 "Usage: %s [view options][-ga angles][-p picture][[rtrace options] octree]\n",
00228                      progname);
00229        exit(1);
00230 }
00231 
00232 
00233 static int
00234 angcmp(              /* compare two angles */
00235        const void    *ap1,
00236        const void    *ap2
00237 )
00238 {
00239        register int  a1, a2;
00240 
00241        a1 = *(ANGLE *)ap1;
00242        a2 = *(ANGLE *)ap2;
00243        if (a1 == a2) {
00244               fprintf(stderr, "%s: duplicate glare angle (%d)\n",
00245                             progname, a1);
00246               exit(1);
00247        }
00248        return(a1-a2);
00249 }
00250 
00251 
00252 static void
00253 init(void)                         /* initialize global variables */
00254 {
00255        double d;
00256        register int  i;
00257 
00258        if (verbose)
00259               fprintf(stderr, "%s: initializing data structures...\n",
00260                             progname);
00261                                           /* set direction vectors */
00262        for (i = 0; glarang[i] != AEND; i++)
00263               ;
00264        qsort(glarang, i, sizeof(ANGLE), angcmp);
00265        if (i > 0 && (glarang[0] <= 0 || glarang[i-1] > 180)) {
00266               fprintf(stderr, "%s: glare angles must be between 1 and 180\n",
00267                             progname);
00268               exit(1);
00269        }
00270        nglarangs = i;
00271        /* nglardirs = 2*nglarangs + 1; */
00272        /* vsize = sampdens - 1; */
00273        if (nglarangs > 0)
00274               maxtheta = (PI/180.)*glarang[nglarangs-1];
00275        else
00276               maxtheta = 0.0;
00277        hsize = hlim(0) + sampdens - 1;
00278        if (hsize > (int)(PI*sampdens))
00279               hsize = PI*sampdens;
00280        indirect = (struct illum *)calloc(nglardirs, sizeof(struct illum));
00281        if (indirect == NULL)
00282               memerr("indirect illuminances");
00283        npixinvw = npixmiss = 0L;
00284        leftview = ourview;
00285        rightview = ourview;
00286        spinvector(leftview.vdir, ourview.vdir, ourview.vup, maxtheta);
00287        spinvector(rightview.vdir, ourview.vdir, ourview.vup, -maxtheta);
00288        setview(&leftview);
00289        setview(&rightview);
00290        indirect[nglarangs].lcos =
00291        indirect[nglarangs].rcos = cos(maxtheta);
00292        indirect[nglarangs].rsin =
00293        -(indirect[nglarangs].lsin = sin(maxtheta));
00294        indirect[nglarangs].theta = 0.0;
00295        for (i = 0; i < nglarangs; i++) {
00296               d = (glarang[nglarangs-1] - glarang[i])*(PI/180.);
00297               indirect[nglarangs-i-1].lcos =
00298               indirect[nglarangs+i+1].rcos = cos(d);
00299               indirect[nglarangs+i+1].rsin =
00300               -(indirect[nglarangs-i-1].lsin = sin(d));
00301               d = (glarang[nglarangs-1] + glarang[i])*(PI/180.);
00302               indirect[nglarangs-i-1].rcos =
00303               indirect[nglarangs+i+1].lcos = cos(d);
00304               indirect[nglarangs-i-1].rsin =
00305               -(indirect[nglarangs+i+1].lsin = sin(d));
00306               indirect[nglarangs-i-1].theta = (PI/180.)*glarang[i];
00307               indirect[nglarangs+i+1].theta = -(PI/180.)*glarang[i];
00308        }
00309                                           /* open picture */
00310        if (picture != NULL) {
00311               if (verbose)
00312                      fprintf(stderr, "%s: opening picture file \"%s\"...\n",
00313                                    progname, picture);
00314               open_pict(picture);
00315        }
00316                                           /* start rtrace */
00317        if (octree != NULL) {
00318               if (verbose) {
00319                      fprintf(stderr,
00320                             "%s: starting luminance calculation...\n\t",
00321                                    progname);
00322                      printargs(rtargc, rtargv, stderr);
00323               }
00324               fork_rtrace(rtargv);
00325        }
00326 }
00327 
00328 
00329 static void
00330 cleanup(void)                      /* close files, wait for children */
00331 {
00332        if (verbose)
00333               fprintf(stderr, "%s: cleaning up...        \n", progname);
00334        if (picture != NULL)
00335               close_pict();
00336        if (octree != NULL)
00337               done_rtrace();
00338        if (npixinvw < 100*npixmiss)
00339               fprintf(stderr, "%s: warning -- missing %d%% of samples\n",
00340                             progname, (int)(100L*npixmiss/npixinvw));
00341 }
00342 
00343 
00344 extern int
00345 compdir(                    /* compute direction for x,y */
00346        FVECT  vd,
00347        int    x,
00348        int    y
00349 )
00350 {
00351        int    hl;
00352        FVECT  org;                 /* dummy variable */
00353 
00354        hl = hlim(y);
00355        if (x <= -hl) {                    /* left region */
00356               if (x <= -hl-sampdens)
00357                      return(-1);
00358               return(viewray(org, vd, &leftview,
00359                             (double)(x+hl)/(2*sampdens)+.5,
00360                             (double)y/(2*sampdens)+.5));
00361        }
00362        if (x >= hl) {                     /* right region */
00363               if (x >= hl+sampdens)
00364                      return(-1);
00365               return(viewray(org, vd, &rightview,
00366                             (double)(x-hl)/(2*sampdens)+.5,
00367                             (double)y/(2*sampdens)+.5));
00368        }
00369                                    /* central region */
00370        if (viewray(org, vd, &ourview, .5, (double)y/(2*sampdens)+.5) < 0)
00371               return(-1);
00372        spinvector(vd, vd, ourview.vup, h_theta(x,y));
00373        return(0);
00374 }
00375 
00376 
00377 extern double
00378 pixsize(             /* return the solid angle of pixel at (x,y) */
00379        int    x,
00380        int    y
00381 )
00382 {
00383        register int  hl, xo;
00384        double disc;
00385 
00386        hl = hlim(y);
00387        if (x < -hl)
00388               xo = x+hl;
00389        else if (x > hl)
00390               xo = x-hl;
00391        else
00392               xo = 0;
00393        disc = 1. - (double)((long)xo*xo + (long)y*y)/((long)sampdens*sampdens);
00394        if (disc <= FTINY*FTINY)
00395               return(0.);
00396        return(1./(sampdens*sampdens*sqrt(disc)));
00397 }
00398 
00399 
00400 extern void
00401 memerr(                     /* malloc failure */
00402        char   *s
00403 )
00404 {
00405        fprintf(stderr, "%s: out of memory for %s\n", progname, s);
00406        exit(1);
00407 }
00408 
00409 
00410 static void
00411 printsources(void)                 /* print out glare sources */
00412 {
00413        register struct source      *sp;
00414 
00415        printf("BEGIN glare source\n");
00416        for (sp = donelist; sp != NULL; sp = sp->next)
00417               printf("\t%f %f %f\t%f\t%f\n",
00418                             sp->dir[0], sp->dir[1], sp->dir[2],
00419                             sp->dom, sp->brt);
00420        printf("END glare source\n");
00421 }
00422 
00423 
00424 static void
00425 printillum(void)                   /* print out indirect illuminances */
00426 {
00427        register int  i;
00428 
00429        printf("BEGIN indirect illuminance\n");
00430        for (i = 0; i < nglardirs; i++)
00431               if (indirect[i].n > FTINY)
00432                      printf("\t%.0f\t%f\n", (180.0/PI)*indirect[i].theta,
00433                                    PI * indirect[i].sum / indirect[i].n);
00434        printf("END indirect illuminance\n");
00435 }