| 1 | #ifndef lint | 
| 2 | static const char RCSid[] = "$Id: rfluxmtx.c,v 2.9 2014/07/25 16:58:20 greg Exp $"; | 
| 3 | #endif | 
| 4 | /* | 
| 5 | * Calculate flux transfer matrix or matrices using rcontrib | 
| 6 | */ | 
| 7 |  | 
| 8 | #include "copyright.h" | 
| 9 |  | 
| 10 | #include <ctype.h> | 
| 11 | #include <stdlib.h> | 
| 12 | #include "rtio.h" | 
| 13 | #include "rtmath.h" | 
| 14 | #include "bsdf.h" | 
| 15 | #include "bsdf_m.h" | 
| 16 | #include "random.h" | 
| 17 | #include "triangulate.h" | 
| 18 | #include "platform.h" | 
| 19 |  | 
| 20 | #ifdef getc_unlocked            /* avoid horrendous overhead of flockfile */ | 
| 21 | #undef getc | 
| 22 | #define getc    getc_unlocked | 
| 23 | #endif | 
| 24 |  | 
| 25 | #ifdef _WIN32 | 
| 26 | #define SPECIALS        " \t\"$*?" | 
| 27 | #define QUOTCHAR        '"' | 
| 28 | #else | 
| 29 | #define SPECIALS        " \t\n'\"()${}*?[];|&" | 
| 30 | #define QUOTCHAR        '\'' | 
| 31 | #define ALTQUOT         '"' | 
| 32 | #endif | 
| 33 |  | 
| 34 | #define MAXRCARG        512 | 
| 35 |  | 
| 36 | char            *progname;              /* global argv[0] */ | 
| 37 |  | 
| 38 | int             verbose = 0;            /* verbose mode? */ | 
| 39 |  | 
| 40 | char            *rcarg[MAXRCARG+1] = {"rcontrib", "-fo+"}; | 
| 41 | int             nrcargs = 2; | 
| 42 |  | 
| 43 | const char      overflowerr[] = "%s: too many arguments!\n"; | 
| 44 |  | 
| 45 | #define CHECKARGC(n)    if (nrcargs >= MAXRCARG-(n)) \ | 
| 46 | { fprintf(stderr, overflowerr, progname); exit(1); } | 
| 47 |  | 
| 48 | int             sampcnt = 0;            /* sample count (0==unset) */ | 
| 49 |  | 
| 50 | char            *reinhfn = "reinhartb.cal"; | 
| 51 | char            *shirchiufn = "disk2square.cal"; | 
| 52 | char            *kfullfn = "klems_full.cal"; | 
| 53 | char            *khalffn = "klems_half.cal"; | 
| 54 | char            *kquarterfn = "klems_quarter.cal"; | 
| 55 |  | 
| 56 | /* string indicating parameters */ | 
| 57 | const char      PARAMSTART[] = "@rfluxmtx"; | 
| 58 |  | 
| 59 | /* surface type IDs */ | 
| 60 | #define ST_NONE         0 | 
| 61 | #define ST_POLY         1 | 
| 62 | #define ST_RING         2 | 
| 63 | #define ST_SOURCE       3 | 
| 64 |  | 
| 65 | typedef struct surf_s { | 
| 66 | struct surf_s   *next;          /* next surface in list */ | 
| 67 | void            *priv;          /* private data (malloc'ed) */ | 
| 68 | char            sname[32];      /* surface name */ | 
| 69 | FVECT           snrm;           /* surface normal */ | 
| 70 | double          area;           /* surface area / proj. solid angle */ | 
| 71 | short           styp;           /* surface type */ | 
| 72 | short           nfargs;         /* number of real arguments */ | 
| 73 | double          farg[1];        /* real values (extends struct) */ | 
| 74 | } SURF;                         /* surface structure */ | 
| 75 |  | 
| 76 | typedef struct { | 
| 77 | FVECT   uva[2];                 /* tangent axes */ | 
| 78 | int     ntris;                  /* number of triangles */ | 
| 79 | struct ptri { | 
| 80 | float   afrac;                  /* fraction of total area */ | 
| 81 | short   vndx[3];                /* vertex indices */ | 
| 82 | }       tri[1];                 /* triangle array (extends struct) */ | 
| 83 | } POLYTRIS;                     /* triangulated polygon */ | 
| 84 |  | 
| 85 | typedef struct param_s { | 
| 86 | char            hemis[32];      /* hemispherical sampling spec. */ | 
| 87 | int             hsiz;           /* hemisphere basis size */ | 
| 88 | int             nsurfs;         /* number of surfaces */ | 
| 89 | SURF            *slist;         /* list of surfaces */ | 
| 90 | FVECT           vup;            /* up vector (zero if unset) */ | 
| 91 | FVECT           nrm;            /* average normal direction */ | 
| 92 | FVECT           udir, vdir;     /* v-up tangent axes */ | 
| 93 | char            *outfn;         /* output file name (receiver) */ | 
| 94 | int             (*sample_basis)(struct param_s *p, int, FILE *); | 
| 95 | } PARAMS;                       /* sender/receiver parameters */ | 
| 96 |  | 
| 97 | PARAMS          curparams; | 
| 98 | char            curmod[128]; | 
| 99 | char            newparams[1024]; | 
| 100 |  | 
| 101 | typedef int     SURFSAMP(FVECT, SURF *, double); | 
| 102 |  | 
| 103 | static SURFSAMP ssamp_bad, ssamp_poly, ssamp_ring; | 
| 104 |  | 
| 105 | SURFSAMP        *orig_in_surf[4] = { | 
| 106 | ssamp_bad, ssamp_poly, ssamp_ring, ssamp_bad | 
| 107 | }; | 
| 108 |  | 
| 109 | /* Clear parameter set */ | 
| 110 | static void | 
| 111 | clear_params(PARAMS *p, int reset_only) | 
| 112 | { | 
| 113 | while (p->slist != NULL) { | 
| 114 | SURF    *sdel = p->slist; | 
| 115 | p->slist = sdel->next; | 
| 116 | if (sdel->priv != NULL) | 
| 117 | free(sdel->priv); | 
| 118 | free(sdel); | 
| 119 | } | 
| 120 | if (reset_only) { | 
| 121 | p->nsurfs = 0; | 
| 122 | memset(p->nrm, 0, sizeof(FVECT)); | 
| 123 | memset(p->vup, 0, sizeof(FVECT)); | 
| 124 | p->outfn = NULL; | 
| 125 | return; | 
| 126 | } | 
| 127 | memset(p, 0, sizeof(PARAMS)); | 
| 128 | } | 
| 129 |  | 
| 130 | /* Get surface type from name */ | 
| 131 | static int | 
| 132 | surf_type(const char *otype) | 
| 133 | { | 
| 134 | if (!strcmp(otype, "polygon")) | 
| 135 | return(ST_POLY); | 
| 136 | if (!strcmp(otype, "ring")) | 
| 137 | return(ST_RING); | 
| 138 | if (!strcmp(otype, "source")) | 
| 139 | return(ST_SOURCE); | 
| 140 | return(ST_NONE); | 
| 141 | } | 
| 142 |  | 
| 143 | /* Add arguments to oconv command */ | 
| 144 | static char * | 
| 145 | oconv_command(int ac, char *av[]) | 
| 146 | { | 
| 147 | static char     oconvbuf[2048] = "!oconv -f"; | 
| 148 | char            *cp = oconvbuf + 9; | 
| 149 |  | 
| 150 | while (ac-- > 0) { | 
| 151 | if (cp >= oconvbuf+(sizeof(oconvbuf)-32)) { | 
| 152 | fputs(progname, stderr); | 
| 153 | fputs(": too many file arguments!\n", stderr); | 
| 154 | exit(1); | 
| 155 | } | 
| 156 | *cp++ = ' '; | 
| 157 | strcpy(cp, *av++); | 
| 158 | while (*cp) cp++; | 
| 159 | } | 
| 160 | *cp = '\0'; | 
| 161 | return(oconvbuf); | 
| 162 | } | 
| 163 |  | 
| 164 | /* Check if any of the characters in str2 are found in str1 */ | 
| 165 | static int | 
| 166 | matchany(const char *str1, const char *str2) | 
| 167 | { | 
| 168 | while (*str1) { | 
| 169 | const char      *cp = str2; | 
| 170 | while (*cp) | 
| 171 | if (*cp++ == *str1) | 
| 172 | return(*str1); | 
| 173 | ++str1; | 
| 174 | } | 
| 175 | return(0); | 
| 176 | } | 
| 177 |  | 
| 178 |  | 
| 179 | /* Convert a set of arguments into a command line for pipe() or system() */ | 
| 180 | static char * | 
| 181 | convert_commandline(char *cmd, const int len, char *av[]) | 
| 182 | { | 
| 183 | int     match; | 
| 184 | char    *cp; | 
| 185 |  | 
| 186 | for (cp = cmd; *av != NULL; av++) { | 
| 187 | const int       n = strlen(*av); | 
| 188 | if (cp+n >= cmd+(len-3)) { | 
| 189 | fputs(progname, stderr); | 
| 190 | return(NULL); | 
| 191 | } | 
| 192 | if ((match = matchany(*av, SPECIALS))) { | 
| 193 | const int       quote = | 
| 194 | #ifdef ALTQUOT | 
| 195 | (match == QUOTCHAR) ? ALTQUOT : | 
| 196 | #endif | 
| 197 | QUOTCHAR; | 
| 198 | *cp++ = quote; | 
| 199 | strcpy(cp, *av); | 
| 200 | cp += n; | 
| 201 | *cp++ = quote; | 
| 202 | } else { | 
| 203 | strcpy(cp, *av); | 
| 204 | cp += n; | 
| 205 | } | 
| 206 | *cp++ = ' '; | 
| 207 | } | 
| 208 | if (cp <= cmd) | 
| 209 | return(NULL); | 
| 210 | *--cp = '\0'; | 
| 211 | return(cmd); | 
| 212 | } | 
| 213 |  | 
| 214 | /* Open a pipe to/from a command given as an argument list */ | 
| 215 | static FILE * | 
| 216 | popen_arglist(char *av[], char *mode) | 
| 217 | { | 
| 218 | char    cmd[10240]; | 
| 219 |  | 
| 220 | if (!convert_commandline(cmd, sizeof(cmd), av)) { | 
| 221 | fputs(progname, stderr); | 
| 222 | fputs(": command line too long in popen_arglist()\n", stderr); | 
| 223 | return(NULL); | 
| 224 | } | 
| 225 | if (verbose) | 
| 226 | fprintf(stderr, "%s: opening pipe %s: %s\n", | 
| 227 | progname, (*mode=='w') ? "to" : "from", cmd); | 
| 228 | return(popen(cmd, mode)); | 
| 229 | } | 
| 230 |  | 
| 231 | #ifdef _WIN32 | 
| 232 | /* Execute system command (Windows version) */ | 
| 233 | static int | 
| 234 | my_exec(char *av[]) | 
| 235 | { | 
| 236 | char    cmd[10240]; | 
| 237 |  | 
| 238 | if (!convert_commandline(cmd, sizeof(cmd), av)) { | 
| 239 | fputs(progname, stderr); | 
| 240 | fputs(": command line too long in my_exec()\n", stderr); | 
| 241 | return(1); | 
| 242 | } | 
| 243 | if (verbose) | 
| 244 | fprintf(stderr, "%s: running: %s\n", progname, cmd); | 
| 245 | return(system(cmd)); | 
| 246 | } | 
| 247 | #else | 
| 248 | /* Execute system command in our stead (Unix version) */ | 
| 249 | static int | 
| 250 | my_exec(char *av[]) | 
| 251 | { | 
| 252 | char    *compath; | 
| 253 |  | 
| 254 | if ((compath = getpath((char *)av[0], getenv("PATH"), X_OK)) == NULL) { | 
| 255 | fprintf(stderr, "%s: cannot locate %s\n", progname, av[0]); | 
| 256 | return(1); | 
| 257 | } | 
| 258 | if (verbose) { | 
| 259 | char    cmd[4096]; | 
| 260 | if (!convert_commandline(cmd, sizeof(cmd), av)) | 
| 261 | strcpy(cmd, "COMMAND TOO LONG TO SHOW"); | 
| 262 | fprintf(stderr, "%s: running: %s\n", progname, cmd); | 
| 263 | } | 
| 264 | execv(compath, av);     /* successful call never returns */ | 
| 265 | perror(compath); | 
| 266 | return(1); | 
| 267 | } | 
| 268 | #endif | 
| 269 |  | 
| 270 | /* Get normalized direction vector from string specification */ | 
| 271 | static int | 
| 272 | get_direction(FVECT dv, const char *s) | 
| 273 | { | 
| 274 | int     sign = 1; | 
| 275 | int     axis = 0; | 
| 276 |  | 
| 277 | memset(dv, 0, sizeof(FVECT)); | 
| 278 | nextchar: | 
| 279 | switch (*s) { | 
| 280 | case '+': | 
| 281 | ++s; | 
| 282 | goto nextchar; | 
| 283 | case '-': | 
| 284 | sign = -sign; | 
| 285 | ++s; | 
| 286 | goto nextchar; | 
| 287 | case 'z': | 
| 288 | case 'Z': | 
| 289 | ++axis; | 
| 290 | case 'y': | 
| 291 | case 'Y': | 
| 292 | ++axis; | 
| 293 | case 'x': | 
| 294 | case 'X': | 
| 295 | dv[axis] = sign; | 
| 296 | return(!s[1] | isspace(s[1])); | 
| 297 | default: | 
| 298 | break; | 
| 299 | } | 
| 300 | #ifdef SMLFLT | 
| 301 | if (sscanf(s, "%f,%f,%f", &dv[0], &dv[1], &dv[2]) != 3) | 
| 302 | #else | 
| 303 | if (sscanf(s, "%lf,%lf,%lf", &dv[0], &dv[1], &dv[2]) != 3) | 
| 304 | #endif | 
| 305 | return(0); | 
| 306 | dv[0] *= (RREAL)sign; | 
| 307 | return(normalize(dv) > 0); | 
| 308 | } | 
| 309 |  | 
| 310 | /* Parse program parameters (directives) */ | 
| 311 | static int | 
| 312 | parse_params(PARAMS *p, char *pargs) | 
| 313 | { | 
| 314 | char    *cp = pargs; | 
| 315 | int     nparams = 0; | 
| 316 | int     i; | 
| 317 |  | 
| 318 | for ( ; ; ) | 
| 319 | switch (*cp++) { | 
| 320 | case 'h': | 
| 321 | if (*cp++ != '=') | 
| 322 | break; | 
| 323 | p->hsiz = 0; | 
| 324 | i = 0; | 
| 325 | while (*cp && !isspace(*cp)) { | 
| 326 | if (isdigit(*cp)) | 
| 327 | p->hsiz = 10*p->hsiz + *cp - '0'; | 
| 328 | p->hemis[i++] = *cp++; | 
| 329 | } | 
| 330 | if (!i) | 
| 331 | break; | 
| 332 | p->hemis[i] = '\0'; | 
| 333 | p->hsiz += !p->hsiz; | 
| 334 | ++nparams; | 
| 335 | continue; | 
| 336 | case 'u': | 
| 337 | if (*cp++ != '=') | 
| 338 | break; | 
| 339 | if (!get_direction(p->vup, cp)) | 
| 340 | break; | 
| 341 | ++nparams; | 
| 342 | continue; | 
| 343 | case 'o': | 
| 344 | if (*cp++ != '=') | 
| 345 | break; | 
| 346 | i = 0; | 
| 347 | while (*cp && !isspace(*cp++)) | 
| 348 | i++; | 
| 349 | if (!i) | 
| 350 | break; | 
| 351 | *--cp = '\0'; | 
| 352 | p->outfn = savqstr(cp-i); | 
| 353 | *cp++ = ' '; | 
| 354 | ++nparams; | 
| 355 | continue; | 
| 356 | case ' ': | 
| 357 | case '\t': | 
| 358 | case '\r': | 
| 359 | case '\n': | 
| 360 | continue; | 
| 361 | case '\0': | 
| 362 | return(nparams); | 
| 363 | default: | 
| 364 | break; | 
| 365 | } | 
| 366 | fprintf(stderr, "%s: bad parameter string '%s'\n", progname, pargs); | 
| 367 | exit(1); | 
| 368 | return(-1);     /* pro forma return */ | 
| 369 | } | 
| 370 |  | 
| 371 | /* Add receiver arguments (directives) corresponding to the current modifier */ | 
| 372 | static void | 
| 373 | finish_receiver(void) | 
| 374 | { | 
| 375 | char    sbuf[256]; | 
| 376 | int     uniform = 0; | 
| 377 | char    *calfn = NULL; | 
| 378 | char    *params = NULL; | 
| 379 | char    *binv = NULL; | 
| 380 | char    *binf = NULL; | 
| 381 | char    *nbins = NULL; | 
| 382 |  | 
| 383 | if (!curmod[0]) { | 
| 384 | fputs(progname, stderr); | 
| 385 | fputs(": missing receiver surface!\n", stderr); | 
| 386 | exit(1); | 
| 387 | } | 
| 388 | if (curparams.outfn != NULL) {  /* add output file spec. */ | 
| 389 | CHECKARGC(2); | 
| 390 | rcarg[nrcargs++] = "-o"; | 
| 391 | rcarg[nrcargs++] = curparams.outfn; | 
| 392 | } | 
| 393 | /* check arguments */ | 
| 394 | if (!curparams.hemis[0]) { | 
| 395 | fputs(progname, stderr); | 
| 396 | fputs(": missing hemisphere sampling type!\n", stderr); | 
| 397 | exit(1); | 
| 398 | } | 
| 399 | if (normalize(curparams.nrm) == 0) { | 
| 400 | fputs(progname, stderr); | 
| 401 | fputs(": undefined normal for hemisphere sampling\n", stderr); | 
| 402 | exit(1); | 
| 403 | } | 
| 404 | if (normalize(curparams.vup) == 0) { | 
| 405 | if (fabs(curparams.nrm[2]) < .7) | 
| 406 | curparams.vup[2] = 1; | 
| 407 | else | 
| 408 | curparams.vup[1] = 1; | 
| 409 | } | 
| 410 | /* determine sample type/bin */ | 
| 411 | if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1') { | 
| 412 | binv = "0";             /* uniform sampling -- one bin */ | 
| 413 | uniform = 1; | 
| 414 | } else if (tolower(curparams.hemis[0]) == 's' && | 
| 415 | tolower(curparams.hemis[1]) == 'c') { | 
| 416 | /* assign parameters */ | 
| 417 | if (curparams.hsiz <= 1) { | 
| 418 | fputs(progname, stderr); | 
| 419 | fputs(": missing size for Shirley-Chiu sampling!\n", stderr); | 
| 420 | exit(1); | 
| 421 | } | 
| 422 | calfn = shirchiufn; shirchiufn = NULL; | 
| 423 | sprintf(sbuf, "SCdim=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g", | 
| 424 | curparams.hsiz, | 
| 425 | curparams.nrm[0], curparams.nrm[1], curparams.nrm[2], | 
| 426 | curparams.vup[0], curparams.vup[1], curparams.vup[2]); | 
| 427 | params = savqstr(sbuf); | 
| 428 | binv = "scbin"; | 
| 429 | nbins = "SCdim*SCdim"; | 
| 430 | } else if ((tolower(curparams.hemis[0]) == 'r') | | 
| 431 | (tolower(curparams.hemis[0]) == 't')) { | 
| 432 | calfn = reinhfn; reinhfn = NULL; | 
| 433 | sprintf(sbuf, "MF=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g", | 
| 434 | curparams.hsiz, | 
| 435 | curparams.nrm[0], curparams.nrm[1], curparams.nrm[2], | 
| 436 | curparams.vup[0], curparams.vup[1], curparams.vup[2]); | 
| 437 | params = savqstr(sbuf); | 
| 438 | binv = "rbin"; | 
| 439 | nbins = "Nrbins"; | 
| 440 | } else if (tolower(curparams.hemis[0]) == 'k' && | 
| 441 | !curparams.hemis[1] | | 
| 442 | (tolower(curparams.hemis[1]) == 'f') | | 
| 443 | (curparams.hemis[1] == '1')) { | 
| 444 | calfn = kfullfn; kfullfn = NULL; | 
| 445 | binf = "kbin"; | 
| 446 | nbins = "Nkbins"; | 
| 447 | } else if (tolower(curparams.hemis[0]) == 'k' && | 
| 448 | (tolower(curparams.hemis[1]) == 'h') | | 
| 449 | (curparams.hemis[1] == '2')) { | 
| 450 | calfn = khalffn; khalffn = NULL; | 
| 451 | binf = "khbin"; | 
| 452 | nbins = "Nkhbins"; | 
| 453 | } else if (tolower(curparams.hemis[0]) == 'k' && | 
| 454 | (tolower(curparams.hemis[1]) == 'q') | | 
| 455 | (curparams.hemis[1] == '4')) { | 
| 456 | calfn = kquarterfn; kquarterfn = NULL; | 
| 457 | binf = "kqbin"; | 
| 458 | nbins = "Nkqbins"; | 
| 459 | } else { | 
| 460 | fprintf(stderr, "%s: unrecognized hemisphere sampling: h=%s\n", | 
| 461 | progname, curparams.hemis); | 
| 462 | exit(1); | 
| 463 | } | 
| 464 | if (!uniform & (curparams.slist->styp == ST_SOURCE)) { | 
| 465 | SURF    *sp; | 
| 466 | for (sp = curparams.slist; sp != NULL; sp = sp->next) | 
| 467 | if (fabs(sp->area - PI) > 1e-3) { | 
| 468 | fprintf(stderr, "%s: source '%s' must be 180-degrees\n", | 
| 469 | progname, sp->sname); | 
| 470 | exit(1); | 
| 471 | } | 
| 472 | } | 
| 473 | if (calfn != NULL) {            /* add cal file if needed */ | 
| 474 | CHECKARGC(2); | 
| 475 | rcarg[nrcargs++] = "-f"; | 
| 476 | rcarg[nrcargs++] = calfn; | 
| 477 | } | 
| 478 | if (params != NULL) {           /* parameters _after_ cal file */ | 
| 479 | CHECKARGC(2); | 
| 480 | rcarg[nrcargs++] = "-p"; | 
| 481 | rcarg[nrcargs++] = params; | 
| 482 | } | 
| 483 | if (nbins != NULL) {            /* add #bins if set */ | 
| 484 | CHECKARGC(2); | 
| 485 | rcarg[nrcargs++] = "-bn"; | 
| 486 | rcarg[nrcargs++] = nbins; | 
| 487 | } | 
| 488 | if (binv != NULL) { | 
| 489 | CHECKARGC(2);           /* assign bin variable */ | 
| 490 | rcarg[nrcargs++] = "-b"; | 
| 491 | rcarg[nrcargs++] = binv; | 
| 492 | } else if (binf != NULL) { | 
| 493 | CHECKARGC(2);           /* assign bin function */ | 
| 494 | rcarg[nrcargs++] = "-b"; | 
| 495 | sprintf(sbuf, "%s(%g,%g,%g,%g,%g,%g)", binf, | 
| 496 | curparams.nrm[0], curparams.nrm[1], curparams.nrm[2], | 
| 497 | curparams.vup[0], curparams.vup[1], curparams.vup[2]); | 
| 498 | rcarg[nrcargs++] = savqstr(sbuf); | 
| 499 | } | 
| 500 | CHECKARGC(2);                           /* modifier argument goes last */ | 
| 501 | rcarg[nrcargs++] = "-m"; | 
| 502 | rcarg[nrcargs++] = savqstr(curmod); | 
| 503 | } | 
| 504 |  | 
| 505 | /* Make randomly oriented tangent plane axes for given normal direction */ | 
| 506 | static void | 
| 507 | make_axes(FVECT uva[2], const FVECT nrm) | 
| 508 | { | 
| 509 | int     i; | 
| 510 |  | 
| 511 | uva[1][0] = 0.5 - frandom(); | 
| 512 | uva[1][1] = 0.5 - frandom(); | 
| 513 | uva[1][2] = 0.5 - frandom(); | 
| 514 | for (i = 3; i--; ) | 
| 515 | if ((-0.6 < nrm[i]) & (nrm[i] < 0.6)) | 
| 516 | break; | 
| 517 | if (i < 0) { | 
| 518 | fputs(progname, stderr); | 
| 519 | fputs(": bad surface normal in make_axes!\n", stderr); | 
| 520 | exit(1); | 
| 521 | } | 
| 522 | uva[1][i] = 1.0; | 
| 523 | VCROSS(uva[0], uva[1], nrm); | 
| 524 | normalize(uva[0]); | 
| 525 | VCROSS(uva[1], nrm, uva[0]); | 
| 526 | } | 
| 527 |  | 
| 528 | /* Illegal sender surfaces end up here */ | 
| 529 | static int | 
| 530 | ssamp_bad(FVECT orig, SURF *sp, double x) | 
| 531 | { | 
| 532 | fprintf(stderr, "%s: illegal sender surface '%s'\n", | 
| 533 | progname, sp->sname); | 
| 534 | return(0); | 
| 535 | } | 
| 536 |  | 
| 537 | /* Generate origin on ring surface from uniform random variable */ | 
| 538 | static int | 
| 539 | ssamp_ring(FVECT orig, SURF *sp, double x) | 
| 540 | { | 
| 541 | FVECT   *uva = (FVECT *)sp->priv; | 
| 542 | double  samp2[2]; | 
| 543 | double  uv[2]; | 
| 544 | int     i; | 
| 545 |  | 
| 546 | if (uva == NULL) {              /* need tangent axes */ | 
| 547 | uva = (FVECT *)malloc(sizeof(FVECT)*2); | 
| 548 | if (uva == NULL) { | 
| 549 | fputs(progname, stderr); | 
| 550 | fputs(": out of memory in ssamp_ring!\n", stderr); | 
| 551 | return(0); | 
| 552 | } | 
| 553 | make_axes(uva, sp->snrm); | 
| 554 | sp->priv = (void *)uva; | 
| 555 | } | 
| 556 | SDmultiSamp(samp2, 2, x); | 
| 557 | samp2[0] = sqrt(samp2[0]*sp->area*(1./PI) + sp->farg[6]*sp->farg[6]); | 
| 558 | samp2[1] *= 2.*PI; | 
| 559 | uv[0] = samp2[0]*tcos(samp2[1]); | 
| 560 | uv[1] = samp2[0]*tsin(samp2[1]); | 
| 561 | for (i = 3; i--; ) | 
| 562 | orig[i] = sp->farg[i] + uv[0]*uva[0][i] + uv[1]*uva[1][i]; | 
| 563 | return(1); | 
| 564 | } | 
| 565 |  | 
| 566 | /* Add triangle to polygon's list (call-back function) */ | 
| 567 | static int | 
| 568 | add_triangle(const Vert2_list *tp, int a, int b, int c) | 
| 569 | { | 
| 570 | POLYTRIS        *ptp = (POLYTRIS *)tp->p; | 
| 571 | struct ptri     *trip = ptp->tri + ptp->ntris++; | 
| 572 |  | 
| 573 | trip->vndx[0] = a; | 
| 574 | trip->vndx[1] = b; | 
| 575 | trip->vndx[2] = c; | 
| 576 | return(1); | 
| 577 | } | 
| 578 |  | 
| 579 | /* Generate origin on polygon surface from uniform random variable */ | 
| 580 | static int | 
| 581 | ssamp_poly(FVECT orig, SURF *sp, double x) | 
| 582 | { | 
| 583 | POLYTRIS        *ptp = (POLYTRIS *)sp->priv; | 
| 584 | double          samp2[2]; | 
| 585 | double          *v0, *v1, *v2; | 
| 586 | int             i; | 
| 587 |  | 
| 588 | if (ptp == NULL) {              /* need to triangulate */ | 
| 589 | ptp = (POLYTRIS *)malloc(sizeof(POLYTRIS) + | 
| 590 | sizeof(struct ptri)*(sp->nfargs/3 - 3)); | 
| 591 | if (ptp == NULL) | 
| 592 | goto memerr; | 
| 593 | if (sp->nfargs == 3) {  /* simple case */ | 
| 594 | ptp->ntris = 1; | 
| 595 | ptp->tri[0].vndx[0] = 0; | 
| 596 | ptp->tri[0].vndx[1] = 1; | 
| 597 | ptp->tri[0].vndx[2] = 2; | 
| 598 | ptp->tri[0].afrac = 1; | 
| 599 | } else { | 
| 600 | Vert2_list      *v2l = polyAlloc(sp->nfargs/3); | 
| 601 | if (v2l == NULL) | 
| 602 | goto memerr; | 
| 603 | make_axes(ptp->uva, sp->snrm); | 
| 604 | for (i = v2l->nv; i--; ) { | 
| 605 | v2l->v[i].mX = DOT(sp->farg+3*i, ptp->uva[0]); | 
| 606 | v2l->v[i].mY = DOT(sp->farg+3*i, ptp->uva[1]); | 
| 607 | } | 
| 608 | ptp->ntris = 0; | 
| 609 | v2l->p = (void *)ptp; | 
| 610 | if (!polyTriangulate(v2l, add_triangle)) { | 
| 611 | fprintf(stderr, | 
| 612 | "%s: cannot triangulate polygon '%s'\n", | 
| 613 | progname, sp->sname); | 
| 614 | return(0); | 
| 615 | } | 
| 616 | for (i = ptp->ntris; i--; ) { | 
| 617 | int     a = ptp->tri[i].vndx[0]; | 
| 618 | int     b = ptp->tri[i].vndx[1]; | 
| 619 | int     c = ptp->tri[i].vndx[2]; | 
| 620 | ptp->tri[i].afrac = | 
| 621 | (v2l->v[a].mX*v2l->v[b].mY - | 
| 622 | v2l->v[b].mX*v2l->v[a].mY + | 
| 623 | v2l->v[b].mX*v2l->v[c].mY - | 
| 624 | v2l->v[c].mX*v2l->v[b].mY + | 
| 625 | v2l->v[c].mX*v2l->v[a].mY - | 
| 626 | v2l->v[a].mX*v2l->v[c].mY) / | 
| 627 | (2.*sp->area); | 
| 628 | } | 
| 629 | polyFree(v2l); | 
| 630 | } | 
| 631 | sp->priv = (void *)ptp; | 
| 632 | } | 
| 633 | /* pick triangle by partial area */ | 
| 634 | for (i = 0; i < ptp->ntris && x > ptp->tri[i].afrac; i++) | 
| 635 | x -= ptp->tri[i].afrac; | 
| 636 | SDmultiSamp(samp2, 2, x/ptp->tri[i].afrac); | 
| 637 | samp2[0] *= samp2[1] = sqrt(samp2[1]); | 
| 638 | samp2[1] = 1. - samp2[1]; | 
| 639 | v0 = sp->farg + 3*ptp->tri[i].vndx[0]; | 
| 640 | v1 = sp->farg + 3*ptp->tri[i].vndx[1]; | 
| 641 | v2 = sp->farg + 3*ptp->tri[i].vndx[2]; | 
| 642 | for (i = 3; i--; ) | 
| 643 | orig[i] = v0[i] + samp2[0]*(v1[i] - v0[i]) | 
| 644 | + samp2[1]*(v2[i] - v0[i]) ; | 
| 645 | return(1); | 
| 646 | memerr: | 
| 647 | fputs(progname, stderr); | 
| 648 | fputs(": out of memory in ssamp_poly!\n", stderr); | 
| 649 | return(0); | 
| 650 | } | 
| 651 |  | 
| 652 | /* Compute sample origin based on projected areas of sender subsurfaces */ | 
| 653 | static int | 
| 654 | sample_origin(PARAMS *p, FVECT orig, const FVECT rdir, double x) | 
| 655 | { | 
| 656 | static double   *projsa; | 
| 657 | static int      nall; | 
| 658 | double          tarea = 0; | 
| 659 | int             i; | 
| 660 | SURF            *sp; | 
| 661 | /* special case for lone surface */ | 
| 662 | if (p->nsurfs == 1) { | 
| 663 | sp = p->slist; | 
| 664 | if (DOT(sp->snrm, rdir) >= -FTINY) { | 
| 665 | fprintf(stderr, | 
| 666 | "%s: internal - sample behind sender '%s'\n", | 
| 667 | progname, sp->sname); | 
| 668 | return(0); | 
| 669 | } | 
| 670 | return((*orig_in_surf[sp->styp])(orig, sp, x)); | 
| 671 | } | 
| 672 | if (p->nsurfs > nall) {         /* (re)allocate surface area cache */ | 
| 673 | if (projsa) free(projsa); | 
| 674 | projsa = (double *)malloc(sizeof(double)*p->nsurfs); | 
| 675 | if (!projsa) return(0); | 
| 676 | nall = p->nsurfs; | 
| 677 | } | 
| 678 | /* compute projected areas */ | 
| 679 | for (i = 0, sp = p->slist; sp != NULL; i++, sp = sp->next) { | 
| 680 | projsa[i] = -DOT(sp->snrm, rdir) * sp->area; | 
| 681 | tarea += projsa[i] *= (double)(projsa[i] > FTINY); | 
| 682 | } | 
| 683 | if (tarea <= FTINY) {           /* wrong side of sender? */ | 
| 684 | fputs(progname, stderr); | 
| 685 | fputs(": internal - sample behind all sender elements!\n", | 
| 686 | stderr); | 
| 687 | return(0); | 
| 688 | } | 
| 689 | tarea *= x;                     /* get surface from list */ | 
| 690 | for (i = 0, sp = p->slist; tarea > projsa[i]; sp = sp->next) | 
| 691 | tarea -= projsa[i++]; | 
| 692 | return((*orig_in_surf[sp->styp])(orig, sp, tarea/projsa[i])); | 
| 693 | } | 
| 694 |  | 
| 695 | /* Uniform sample generator */ | 
| 696 | static int | 
| 697 | sample_uniform(PARAMS *p, int b, FILE *fp) | 
| 698 | { | 
| 699 | int     n = sampcnt; | 
| 700 | double  samp3[3]; | 
| 701 | double  duvw[3]; | 
| 702 | FVECT   orig_dir[2]; | 
| 703 | int     i; | 
| 704 |  | 
| 705 | if (fp == NULL)                 /* just requesting number of bins? */ | 
| 706 | return(1); | 
| 707 |  | 
| 708 | while (n--) {                   /* stratified hemisphere sampling */ | 
| 709 | SDmultiSamp(samp3, 3, (n+frandom())/sampcnt); | 
| 710 | SDsquare2disk(duvw, samp3[1], samp3[2]); | 
| 711 | duvw[2] = -sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]); | 
| 712 | for (i = 3; i--; ) | 
| 713 | orig_dir[1][i] = duvw[0]*p->udir[i] + | 
| 714 | duvw[1]*p->vdir[i] + | 
| 715 | duvw[2]*p->nrm[i] ; | 
| 716 | if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0])) | 
| 717 | return(0); | 
| 718 | if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2) | 
| 719 | return(0); | 
| 720 | } | 
| 721 | return(1); | 
| 722 | } | 
| 723 |  | 
| 724 | /* Shirly-Chiu sample generator */ | 
| 725 | static int | 
| 726 | sample_shirchiu(PARAMS *p, int b, FILE *fp) | 
| 727 | { | 
| 728 | int     n = sampcnt; | 
| 729 | double  samp3[3]; | 
| 730 | double  duvw[3]; | 
| 731 | FVECT   orig_dir[2]; | 
| 732 | int     i; | 
| 733 |  | 
| 734 | if (fp == NULL)                 /* just requesting number of bins? */ | 
| 735 | return(p->hsiz*p->hsiz); | 
| 736 |  | 
| 737 | while (n--) {                   /* stratified sampling */ | 
| 738 | SDmultiSamp(samp3, 3, (n+frandom())/sampcnt); | 
| 739 | SDsquare2disk(duvw, (b/p->hsiz + samp3[1])/curparams.hsiz, | 
| 740 | (b%p->hsiz + samp3[2])/curparams.hsiz); | 
| 741 | duvw[2] = sqrt(1. - duvw[0]*duvw[0] - duvw[1]*duvw[1]); | 
| 742 | for (i = 3; i--; ) | 
| 743 | orig_dir[1][i] = -duvw[0]*p->udir[i] - | 
| 744 | duvw[1]*p->vdir[i] - | 
| 745 | duvw[2]*p->nrm[i] ; | 
| 746 | if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0])) | 
| 747 | return(0); | 
| 748 | if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2) | 
| 749 | return(0); | 
| 750 | } | 
| 751 | return(1); | 
| 752 | } | 
| 753 |  | 
| 754 | /* Reinhart/Tregenza sample generator */ | 
| 755 | static int | 
| 756 | sample_reinhart(PARAMS *p, int b, FILE *fp) | 
| 757 | { | 
| 758 | #define T_NALT  7 | 
| 759 | static const int        tnaz[T_NALT] = {30, 30, 24, 24, 18, 12, 6}; | 
| 760 | const int               RowMax = T_NALT*p->hsiz + 1; | 
| 761 | const double            RAH = (.5*PI)/(RowMax-.5); | 
| 762 | #define rnaz(r)                 (r >= RowMax-1 ? 1 : p->hsiz*tnaz[r/p->hsiz]) | 
| 763 | int                     n = sampcnt; | 
| 764 | int                     row, col; | 
| 765 | double                  samp3[3]; | 
| 766 | double                  alt, azi; | 
| 767 | double                  duvw[3]; | 
| 768 | FVECT                   orig_dir[2]; | 
| 769 | int                     i; | 
| 770 |  | 
| 771 | if (fp == NULL) {               /* just requesting number of bins? */ | 
| 772 | n = 0; | 
| 773 | for (row = RowMax; row--; ) n += rnaz(row); | 
| 774 | return(n); | 
| 775 | } | 
| 776 | row = 0;                        /* identify row & column */ | 
| 777 | col = b; | 
| 778 | while (col >= rnaz(row)) { | 
| 779 | col -= rnaz(row); | 
| 780 | ++row; | 
| 781 | } | 
| 782 | while (n--) {                   /* stratified sampling */ | 
| 783 | SDmultiSamp(samp3, 3, (n+frandom())/sampcnt); | 
| 784 | alt = (row+samp3[1])*RAH; | 
| 785 | azi = (2.*PI)*(col+samp3[2]-.5)/rnaz(row); | 
| 786 | duvw[2] = cos(alt);     /* measured from horizon */ | 
| 787 | duvw[0] = tcos(azi)*duvw[2]; | 
| 788 | duvw[1] = tsin(azi)*duvw[2]; | 
| 789 | duvw[2] = sqrt(1. - duvw[2]*duvw[2]); | 
| 790 | for (i = 3; i--; ) | 
| 791 | orig_dir[1][i] = -duvw[0]*p->udir[i] - | 
| 792 | duvw[1]*p->vdir[i] - | 
| 793 | duvw[2]*p->nrm[i] ; | 
| 794 | if (!sample_origin(p, orig_dir[0], orig_dir[1], samp3[0])) | 
| 795 | return(0); | 
| 796 | if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2) | 
| 797 | return(0); | 
| 798 | } | 
| 799 | return(1); | 
| 800 | #undef rnaz | 
| 801 | #undef T_NALT | 
| 802 | } | 
| 803 |  | 
| 804 | /* Klems sample generator */ | 
| 805 | static int | 
| 806 | sample_klems(PARAMS *p, int b, FILE *fp) | 
| 807 | { | 
| 808 | static const char       bname[4][20] = { | 
| 809 | "LBNL/Klems Full", | 
| 810 | "LBNL/Klems Half", | 
| 811 | "INTERNAL ERROR", | 
| 812 | "LBNL/Klems Quarter" | 
| 813 | }; | 
| 814 | static ANGLE_BASIS      *kbasis[4]; | 
| 815 | const int               bi = p->hemis[1] - '1'; | 
| 816 | int                     n = sampcnt; | 
| 817 | double                  samp2[2]; | 
| 818 | double                  duvw[3]; | 
| 819 | FVECT                   orig_dir[2]; | 
| 820 | int                     i; | 
| 821 |  | 
| 822 | if (!kbasis[bi]) {              /* need to get basis, first */ | 
| 823 | for (i = 4; i--; ) | 
| 824 | if (!strcasecmp(abase_list[i].name, bname[bi])) { | 
| 825 | kbasis[bi] = &abase_list[i]; | 
| 826 | break; | 
| 827 | } | 
| 828 | if (i < 0) { | 
| 829 | fprintf(stderr, "%s: unknown hemisphere basis '%s'\n", | 
| 830 | progname, bname[bi]); | 
| 831 | return(0); | 
| 832 | } | 
| 833 | } | 
| 834 | if (fp == NULL)                 /* just requesting number of bins? */ | 
| 835 | return(kbasis[bi]->nangles); | 
| 836 |  | 
| 837 | while (n--) {                   /* stratified sampling */ | 
| 838 | SDmultiSamp(samp2, 2, (n+frandom())/sampcnt); | 
| 839 | if (!bo_getvec(duvw, b+samp2[1], kbasis[bi])) | 
| 840 | return(0); | 
| 841 | for (i = 3; i--; ) | 
| 842 | orig_dir[1][i] = duvw[0]*p->udir[i] + | 
| 843 | duvw[1]*p->vdir[i] + | 
| 844 | duvw[2]*p->nrm[i] ; | 
| 845 | if (!sample_origin(p, orig_dir[0], orig_dir[1], samp2[0])) | 
| 846 | return(0); | 
| 847 | if (fwrite(orig_dir, sizeof(FVECT), 2, fp) != 2) | 
| 848 | return(0); | 
| 849 | } | 
| 850 | return(1); | 
| 851 | } | 
| 852 |  | 
| 853 | /* Prepare hemisphere basis sampler that will send rays to rcontrib */ | 
| 854 | static int | 
| 855 | prepare_sampler(void) | 
| 856 | { | 
| 857 | if (curparams.slist == NULL) {  /* missing sample surface! */ | 
| 858 | fputs(progname, stderr); | 
| 859 | fputs(": no sender surface!\n", stderr); | 
| 860 | return(-1); | 
| 861 | } | 
| 862 | if (curparams.outfn != NULL)    /* misplaced output file spec. */ | 
| 863 | fprintf(stderr, "%s: warning - ignoring output file in sender ('%s')\n", | 
| 864 | progname, curparams.outfn); | 
| 865 | /* check/set basis hemisphere */ | 
| 866 | if (!curparams.hemis[0]) { | 
| 867 | fputs(progname, stderr); | 
| 868 | fputs(": missing sender sampling type!\n", stderr); | 
| 869 | return(-1); | 
| 870 | } | 
| 871 | if (normalize(curparams.nrm) == 0) { | 
| 872 | fputs(progname, stderr); | 
| 873 | fputs(": undefined normal for sender sampling\n", stderr); | 
| 874 | return(-1); | 
| 875 | } | 
| 876 | if (normalize(curparams.vup) == 0) { | 
| 877 | if (fabs(curparams.nrm[2]) < .7) | 
| 878 | curparams.vup[2] = 1; | 
| 879 | else | 
| 880 | curparams.vup[1] = 1; | 
| 881 | } | 
| 882 | VCROSS(curparams.udir, curparams.vup, curparams.nrm); | 
| 883 | if (normalize(curparams.udir) == 0) { | 
| 884 | fputs(progname, stderr); | 
| 885 | fputs(": up vector coincides with sender normal\n", stderr); | 
| 886 | return(-1); | 
| 887 | } | 
| 888 | VCROSS(curparams.vdir, curparams.nrm, curparams.udir); | 
| 889 | if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1') | 
| 890 | curparams.sample_basis = sample_uniform; | 
| 891 | else if (tolower(curparams.hemis[0]) == 's' && | 
| 892 | tolower(curparams.hemis[1]) == 'c') | 
| 893 | curparams.sample_basis = sample_shirchiu; | 
| 894 | else if ((tolower(curparams.hemis[0]) == 'r') | | 
| 895 | (tolower(curparams.hemis[0]) == 't')) | 
| 896 | curparams.sample_basis = sample_reinhart; | 
| 897 | else if (tolower(curparams.hemis[0]) == 'k') { | 
| 898 | switch (curparams.hemis[1]) { | 
| 899 | case '1': | 
| 900 | case '2': | 
| 901 | case '4': | 
| 902 | break; | 
| 903 | case 'f': | 
| 904 | case 'F': | 
| 905 | case '\0': | 
| 906 | curparams.hemis[1] = '1'; | 
| 907 | break; | 
| 908 | case 'h': | 
| 909 | case 'H': | 
| 910 | curparams.hemis[1] = '2'; | 
| 911 | break; | 
| 912 | case 'q': | 
| 913 | case 'Q': | 
| 914 | curparams.hemis[1] = '4'; | 
| 915 | break; | 
| 916 | default: | 
| 917 | goto unrecognized; | 
| 918 | } | 
| 919 | curparams.hemis[2] = '\0'; | 
| 920 | curparams.sample_basis = sample_klems; | 
| 921 | } else | 
| 922 | goto unrecognized; | 
| 923 | /* return number of bins */ | 
| 924 | return((*curparams.sample_basis)(&curparams,0,NULL)); | 
| 925 | unrecognized: | 
| 926 | fprintf(stderr, "%s: unrecognized sender sampling: h=%s\n", | 
| 927 | progname, curparams.hemis); | 
| 928 | return(-1); | 
| 929 | } | 
| 930 |  | 
| 931 | /* Compute normal and area for polygon */ | 
| 932 | static int | 
| 933 | finish_polygon(SURF *p) | 
| 934 | { | 
| 935 | const int       nv = p->nfargs / 3; | 
| 936 | FVECT           e1, e2, vc; | 
| 937 | int             i; | 
| 938 |  | 
| 939 | memset(p->snrm, 0, sizeof(FVECT)); | 
| 940 | VSUB(e1, p->farg+3, p->farg); | 
| 941 | for (i = 2; i < nv; i++) { | 
| 942 | VSUB(e2, p->farg+3*i, p->farg); | 
| 943 | VCROSS(vc, e1, e2); | 
| 944 | p->snrm[0] += vc[0]; | 
| 945 | p->snrm[1] += vc[1]; | 
| 946 | p->snrm[2] += vc[2]; | 
| 947 | VCOPY(e1, e2); | 
| 948 | } | 
| 949 | p->area = normalize(p->snrm)*0.5; | 
| 950 | return(p->area > FTINY); | 
| 951 | } | 
| 952 |  | 
| 953 | /* Add a surface to our current parameters */ | 
| 954 | static void | 
| 955 | add_surface(int st, const char *oname, FILE *fp) | 
| 956 | { | 
| 957 | SURF    *snew; | 
| 958 | int     n; | 
| 959 | /* get floating-point arguments */ | 
| 960 | if (!fscanf(fp, "%d", &n)) return; | 
| 961 | while (n-- > 0) fscanf(fp, "%*s"); | 
| 962 | if (!fscanf(fp, "%d", &n)) return; | 
| 963 | while (n-- > 0) fscanf(fp, "%*d"); | 
| 964 | if (!fscanf(fp, "%d", &n) || n <= 0) return; | 
| 965 | snew = (SURF *)malloc(sizeof(SURF) + sizeof(double)*(n-1)); | 
| 966 | if (snew == NULL) { | 
| 967 | fputs(progname, stderr); | 
| 968 | fputs(": out of memory!\n", stderr); | 
| 969 | exit(1); | 
| 970 | } | 
| 971 | strncpy(snew->sname, oname, sizeof(snew->sname)-1); | 
| 972 | snew->sname[sizeof(snew->sname)-1] = '\0'; | 
| 973 | snew->styp = st; | 
| 974 | snew->priv = NULL; | 
| 975 | snew->nfargs = n; | 
| 976 | for (n = 0; n < snew->nfargs; n++) | 
| 977 | if (fscanf(fp, "%lf", &snew->farg[n]) != 1) { | 
| 978 | fprintf(stderr, "%s: error reading arguments for '%s'\n", | 
| 979 | progname, oname); | 
| 980 | exit(1); | 
| 981 | } | 
| 982 | switch (st) { | 
| 983 | case ST_RING: | 
| 984 | if (snew->nfargs != 8) | 
| 985 | goto badcount; | 
| 986 | VCOPY(snew->snrm, snew->farg+3); | 
| 987 | if (normalize(snew->snrm) == 0) | 
| 988 | goto badnorm; | 
| 989 | if (snew->farg[7] < snew->farg[6]) { | 
| 990 | double  t = snew->farg[7]; | 
| 991 | snew->farg[7] = snew->farg[6]; | 
| 992 | snew->farg[6] = t; | 
| 993 | } | 
| 994 | snew->area = PI*(snew->farg[7]*snew->farg[7] - | 
| 995 | snew->farg[6]*snew->farg[6]); | 
| 996 | break; | 
| 997 | case ST_POLY: | 
| 998 | if (snew->nfargs < 9 || snew->nfargs % 3) | 
| 999 | goto badcount; | 
| 1000 | finish_polygon(snew); | 
| 1001 | break; | 
| 1002 | case ST_SOURCE: | 
| 1003 | if (snew->nfargs != 4) | 
| 1004 | goto badcount; | 
| 1005 | for (n = 3; n--; )      /* need to reverse "normal" */ | 
| 1006 | snew->snrm[n] = -snew->farg[n]; | 
| 1007 | if (normalize(snew->snrm) == 0) | 
| 1008 | goto badnorm; | 
| 1009 | snew->area = sin((PI/180./2.)*snew->farg[3]); | 
| 1010 | snew->area *= PI*snew->area; | 
| 1011 | break; | 
| 1012 | } | 
| 1013 | if (snew->area <= FTINY) { | 
| 1014 | fprintf(stderr, "%s: warning - zero area for surface '%s'\n", | 
| 1015 | progname, oname); | 
| 1016 | free(snew); | 
| 1017 | return; | 
| 1018 | } | 
| 1019 | VSUM(curparams.nrm, curparams.nrm, snew->snrm, snew->area); | 
| 1020 | snew->next = curparams.slist; | 
| 1021 | curparams.slist = snew; | 
| 1022 | curparams.nsurfs++; | 
| 1023 | return; | 
| 1024 | badcount: | 
| 1025 | fprintf(stderr, "%s: bad argument count for surface element '%s'\n", | 
| 1026 | progname, oname); | 
| 1027 | exit(1); | 
| 1028 | badnorm: | 
| 1029 | fprintf(stderr, "%s: bad orientation for surface element '%s'\n", | 
| 1030 | progname, oname); | 
| 1031 | exit(1); | 
| 1032 | } | 
| 1033 |  | 
| 1034 | /* Parse a receiver object (look for modifiers to add to rcontrib command) */ | 
| 1035 | static int | 
| 1036 | add_recv_object(FILE *fp) | 
| 1037 | { | 
| 1038 | int             st; | 
| 1039 | char            thismod[128], otype[32], oname[128]; | 
| 1040 | int             n; | 
| 1041 |  | 
| 1042 | if (fscanf(fp, "%s %s %s", thismod, otype, oname) != 3) | 
| 1043 | return(0);              /* must have hit EOF! */ | 
| 1044 | if (!strcmp(otype, "alias")) { | 
| 1045 | fscanf(fp, "%*s");      /* skip alias */ | 
| 1046 | return(0); | 
| 1047 | } | 
| 1048 | /* is it a new receiver? */ | 
| 1049 | if ((st = surf_type(otype)) != ST_NONE) { | 
| 1050 | if (curparams.slist != NULL && (st == ST_SOURCE) ^ | 
| 1051 | (curparams.slist->styp == ST_SOURCE)) { | 
| 1052 | fputs(progname, stderr); | 
| 1053 | fputs(": cannot mix source/non-source receivers!\n", stderr); | 
| 1054 | return(-1); | 
| 1055 | } | 
| 1056 | if (strcmp(thismod, curmod)) { | 
| 1057 | if (curmod[0]) {        /* output last receiver? */ | 
| 1058 | finish_receiver(); | 
| 1059 | clear_params(&curparams, 1); | 
| 1060 | } | 
| 1061 | parse_params(&curparams, newparams); | 
| 1062 | newparams[0] = '\0'; | 
| 1063 | strcpy(curmod, thismod); | 
| 1064 | } | 
| 1065 | add_surface(st, oname, fp);     /* read & store surface */ | 
| 1066 | return(1); | 
| 1067 | } | 
| 1068 | /* else skip arguments */ | 
| 1069 | if (!fscanf(fp, "%d", &n)) return(0); | 
| 1070 | while (n-- > 0) fscanf(fp, "%*s"); | 
| 1071 | if (!fscanf(fp, "%d", &n)) return(0); | 
| 1072 | while (n-- > 0) fscanf(fp, "%*d"); | 
| 1073 | if (!fscanf(fp, "%d", &n)) return(0); | 
| 1074 | while (n-- > 0) fscanf(fp, "%*f"); | 
| 1075 | return(0); | 
| 1076 | } | 
| 1077 |  | 
| 1078 | /* Parse a sender object */ | 
| 1079 | static int | 
| 1080 | add_send_object(FILE *fp) | 
| 1081 | { | 
| 1082 | int             st; | 
| 1083 | char            otype[32], oname[128]; | 
| 1084 | int             n; | 
| 1085 |  | 
| 1086 | if (fscanf(fp, "%*s %s %s", otype, oname) != 2) | 
| 1087 | return(0);              /* must have hit EOF! */ | 
| 1088 | if (!strcmp(otype, "alias")) { | 
| 1089 | fscanf(fp, "%*s");      /* skip alias */ | 
| 1090 | return(0); | 
| 1091 | } | 
| 1092 | /* is it a new surface? */ | 
| 1093 | if ((st = surf_type(otype)) != ST_NONE) { | 
| 1094 | if (st == ST_SOURCE) { | 
| 1095 | fputs(progname, stderr); | 
| 1096 | fputs(": cannot use source as a sender!\n", stderr); | 
| 1097 | return(-1); | 
| 1098 | } | 
| 1099 | parse_params(&curparams, newparams); | 
| 1100 | newparams[0] = '\0'; | 
| 1101 | add_surface(st, oname, fp);     /* read & store surface */ | 
| 1102 | return(0); | 
| 1103 | } | 
| 1104 | /* else skip arguments */ | 
| 1105 | if (!fscanf(fp, "%d", &n)) return(0); | 
| 1106 | while (n-- > 0) fscanf(fp, "%*s"); | 
| 1107 | if (!fscanf(fp, "%d", &n)) return(0); | 
| 1108 | while (n-- > 0) fscanf(fp, "%*d"); | 
| 1109 | if (!fscanf(fp, "%d", &n)) return(0); | 
| 1110 | while (n-- > 0) fscanf(fp, "%*f"); | 
| 1111 | return(0); | 
| 1112 | } | 
| 1113 |  | 
| 1114 | /* Load a Radiance scene using the given callback function for objects */ | 
| 1115 | static int | 
| 1116 | load_scene(const char *inspec, int (*ocb)(FILE *)) | 
| 1117 | { | 
| 1118 | int     rv = 0; | 
| 1119 | char    inpbuf[1024]; | 
| 1120 | FILE    *fp; | 
| 1121 | int     c; | 
| 1122 |  | 
| 1123 | if (*inspec == '!') | 
| 1124 | fp = popen(inspec+1, "r"); | 
| 1125 | else | 
| 1126 | fp = fopen(inspec, "r"); | 
| 1127 | if (fp == NULL) { | 
| 1128 | fprintf(stderr, "%s: cannot load '%s'\n", progname, inspec); | 
| 1129 | return(-1); | 
| 1130 | } | 
| 1131 | while ((c = getc(fp)) != EOF) { /* load receiver data */ | 
| 1132 | if (isspace(c))         /* skip leading white space */ | 
| 1133 | continue; | 
| 1134 | if (c == '!') {         /* read from a new command */ | 
| 1135 | inpbuf[0] = c; | 
| 1136 | if (fgetline(inpbuf+1, sizeof(inpbuf)-1, fp) != NULL) { | 
| 1137 | if ((c = load_scene(inpbuf, ocb)) < 0) | 
| 1138 | return(c); | 
| 1139 | rv += c; | 
| 1140 | } | 
| 1141 | continue; | 
| 1142 | } | 
| 1143 | if (c == '#') {         /* parameters/comment */ | 
| 1144 | if ((c = getc(fp)) == EOF || ungetc(c,fp) == EOF) | 
| 1145 | break; | 
| 1146 | if (!isspace(c) && fscanf(fp, "%s", inpbuf) == 1 && | 
| 1147 | !strcmp(inpbuf, PARAMSTART)) { | 
| 1148 | if (fgets(inpbuf, sizeof(inpbuf), fp) != NULL) | 
| 1149 | strcat(newparams, inpbuf); | 
| 1150 | continue; | 
| 1151 | } | 
| 1152 | while ((c = getc(fp)) != EOF && c != '\n') | 
| 1153 | ;       /* else skipping comment */ | 
| 1154 | continue; | 
| 1155 | } | 
| 1156 | ungetc(c, fp);          /* else check object for receiver */ | 
| 1157 | c = (*ocb)(fp); | 
| 1158 | if (c < 0) | 
| 1159 | return(c); | 
| 1160 | rv += c; | 
| 1161 | } | 
| 1162 | /* close our input stream */ | 
| 1163 | c = (*inspec == '!') ? pclose(fp) : fclose(fp); | 
| 1164 | if (c != 0) { | 
| 1165 | fprintf(stderr, "%s: error loading '%s'\n", progname, inspec); | 
| 1166 | return(-1); | 
| 1167 | } | 
| 1168 | return(rv); | 
| 1169 | } | 
| 1170 |  | 
| 1171 | /* Get command arguments and run program */ | 
| 1172 | int | 
| 1173 | main(int argc, char *argv[]) | 
| 1174 | { | 
| 1175 | char    fmtopt[6] = "-faa";     /* default output is ASCII */ | 
| 1176 | char    *xrs=NULL, *yrs=NULL, *ldopt=NULL; | 
| 1177 | int     wantIrradiance = 0; | 
| 1178 | char    *sendfn; | 
| 1179 | char    sampcntbuf[32], nsbinbuf[32]; | 
| 1180 | FILE    *rcfp; | 
| 1181 | int     nsbins; | 
| 1182 | int     a, i; | 
| 1183 | /* screen rcontrib options */ | 
| 1184 | progname = argv[0]; | 
| 1185 | for (a = 1; a < argc-2 && argv[a][0] == '-'; a++) { | 
| 1186 | int     na = 1;         /* !! Keep consistent !! */ | 
| 1187 | switch (argv[a][1]) { | 
| 1188 | case 'v':               /* verbose mode */ | 
| 1189 | verbose = !verbose; | 
| 1190 | na = 0; | 
| 1191 | continue; | 
| 1192 | case 'f':               /* special case for -fo, -ff, etc. */ | 
| 1193 | switch (argv[a][2]) { | 
| 1194 | case '\0':              /* cal file */ | 
| 1195 | goto userr; | 
| 1196 | case 'o':               /* force output */ | 
| 1197 | goto userr; | 
| 1198 | case 'a':               /* output format */ | 
| 1199 | case 'f': | 
| 1200 | case 'd': | 
| 1201 | case 'c': | 
| 1202 | if (!(fmtopt[3] = argv[a][3])) | 
| 1203 | fmtopt[3] = argv[a][2]; | 
| 1204 | fmtopt[2] = argv[a][2]; | 
| 1205 | na = 0; | 
| 1206 | continue;       /* will pass later */ | 
| 1207 | default: | 
| 1208 | goto userr; | 
| 1209 | } | 
| 1210 | break; | 
| 1211 | case 'x':               /* x-resolution */ | 
| 1212 | xrs = argv[++a]; | 
| 1213 | na = 0; | 
| 1214 | continue; | 
| 1215 | case 'y':               /* y-resolution */ | 
| 1216 | yrs = argv[++a]; | 
| 1217 | na = 0; | 
| 1218 | continue; | 
| 1219 | case 'c':               /* number of samples */ | 
| 1220 | sampcnt = atoi(argv[++a]); | 
| 1221 | if (sampcnt <= 0) | 
| 1222 | goto userr; | 
| 1223 | na = 0;         /* we re-add this later */ | 
| 1224 | continue; | 
| 1225 | case 'I':               /* only for pass-through mode */ | 
| 1226 | wantIrradiance = 1; | 
| 1227 | na = 0; | 
| 1228 | continue; | 
| 1229 | case 'V':               /* options without arguments */ | 
| 1230 | case 'w': | 
| 1231 | case 'u': | 
| 1232 | case 'i': | 
| 1233 | case 'h': | 
| 1234 | case 'r': | 
| 1235 | break; | 
| 1236 | case 'n':               /* options with 1 argument */ | 
| 1237 | case 's': | 
| 1238 | case 'o': | 
| 1239 | na = 2; | 
| 1240 | break; | 
| 1241 | case 'b':               /* special case */ | 
| 1242 | if (argv[a][2] != 'v') goto userr; | 
| 1243 | break; | 
| 1244 | case 'l':               /* special case */ | 
| 1245 | if (argv[a][2] == 'd') { | 
| 1246 | ldopt = argv[a]; | 
| 1247 | na = 0; | 
| 1248 | continue; | 
| 1249 | } | 
| 1250 | na = 2; | 
| 1251 | break; | 
| 1252 | case 'd':               /* special case */ | 
| 1253 | if (argv[a][2] != 'v') na = 2; | 
| 1254 | break; | 
| 1255 | case 'a':               /* special case */ | 
| 1256 | na = (argv[a][2] == 'v') ? 4 : 2; | 
| 1257 | break; | 
| 1258 | case 'm':               /* special case */ | 
| 1259 | if (!argv[a][2]) goto userr; | 
| 1260 | na = (argv[a][2] == 'e') | (argv[a][2] == 'a') ? 4 : 2; | 
| 1261 | break; | 
| 1262 | case '\0':              /* pass-through mode */ | 
| 1263 | goto done_opts; | 
| 1264 | default:                /* anything else is verbotten */ | 
| 1265 | goto userr; | 
| 1266 | } | 
| 1267 | if (na <= 0) continue; | 
| 1268 | CHECKARGC(na);          /* pass on option */ | 
| 1269 | rcarg[nrcargs++] = argv[a]; | 
| 1270 | while (--na)            /* + arguments if any */ | 
| 1271 | rcarg[nrcargs++] = argv[++a]; | 
| 1272 | } | 
| 1273 | done_opts: | 
| 1274 | if (a > argc-2) | 
| 1275 | goto userr;             /* check at end of options */ | 
| 1276 | sendfn = argv[a++];             /* assign sender & receiver inputs */ | 
| 1277 | if (sendfn[0] == '-') {         /* user wants pass-through mode? */ | 
| 1278 | if (sendfn[1]) goto userr; | 
| 1279 | sendfn = NULL; | 
| 1280 | if (wantIrradiance) { | 
| 1281 | CHECKARGC(1); | 
| 1282 | rcarg[nrcargs++] = "-I"; | 
| 1283 | } | 
| 1284 | if (xrs) { | 
| 1285 | CHECKARGC(2); | 
| 1286 | rcarg[nrcargs++] = "-x"; | 
| 1287 | rcarg[nrcargs++] = xrs; | 
| 1288 | } | 
| 1289 | if (yrs) { | 
| 1290 | CHECKARGC(2); | 
| 1291 | rcarg[nrcargs++] = "-y"; | 
| 1292 | rcarg[nrcargs++] = yrs; | 
| 1293 | } | 
| 1294 | if (ldopt) { | 
| 1295 | CHECKARGC(1); | 
| 1296 | rcarg[nrcargs++] = ldopt; | 
| 1297 | } | 
| 1298 | if (sampcnt <= 0) sampcnt = 1; | 
| 1299 | } else {                        /* else in sampling mode */ | 
| 1300 | if (wantIrradiance) { | 
| 1301 | fputs(progname, stderr); | 
| 1302 | fputs(": -I supported for pass-through only\n", stderr); | 
| 1303 | return(1); | 
| 1304 | } | 
| 1305 | fmtopt[2] = (sizeof(RREAL)==sizeof(double)) ? 'd' : 'f'; | 
| 1306 | if (sampcnt <= 0) sampcnt = 10000; | 
| 1307 | } | 
| 1308 | sprintf(sampcntbuf, "%d", sampcnt); | 
| 1309 | CHECKARGC(3);                   /* add our format & sample count */ | 
| 1310 | rcarg[nrcargs++] = fmtopt; | 
| 1311 | rcarg[nrcargs++] = "-c"; | 
| 1312 | rcarg[nrcargs++] = sampcntbuf; | 
| 1313 | /* add receiver arguments to rcontrib */ | 
| 1314 | if (load_scene(argv[a], add_recv_object) < 0) | 
| 1315 | return(1); | 
| 1316 | finish_receiver(); | 
| 1317 | if (sendfn == NULL) {           /* pass-through mode? */ | 
| 1318 | CHECKARGC(1);           /* add octree */ | 
| 1319 | rcarg[nrcargs++] = oconv_command(argc-a, argv+a); | 
| 1320 | rcarg[nrcargs] = NULL; | 
| 1321 | return(my_exec(rcarg)); /* rcontrib does everything */ | 
| 1322 | } | 
| 1323 | clear_params(&curparams, 0);    /* else load sender surface & params */ | 
| 1324 | if (load_scene(sendfn, add_send_object) < 0) | 
| 1325 | return(1); | 
| 1326 | if ((nsbins = prepare_sampler()) <= 0) | 
| 1327 | return(1); | 
| 1328 | CHECKARGC(3);                   /* add row count and octree */ | 
| 1329 | rcarg[nrcargs++] = "-y"; | 
| 1330 | sprintf(nsbinbuf, "%d", nsbins); | 
| 1331 | rcarg[nrcargs++] = nsbinbuf; | 
| 1332 | rcarg[nrcargs++] = oconv_command(argc-a, argv+a); | 
| 1333 | rcarg[nrcargs] = NULL; | 
| 1334 | /* open pipe to rcontrib process */ | 
| 1335 | if ((rcfp = popen_arglist(rcarg, "w")) == NULL) | 
| 1336 | return(1); | 
| 1337 | SET_FILE_BINARY(rcfp); | 
| 1338 | #ifdef getc_unlocked | 
| 1339 | flockfile(rcfp); | 
| 1340 | #endif | 
| 1341 | if (verbose) { | 
| 1342 | fprintf(stderr, "%s: sampling %d directions", progname, nsbins); | 
| 1343 | if (curparams.nsurfs > 1) | 
| 1344 | fprintf(stderr, " (%d elements)\n", curparams.nsurfs); | 
| 1345 | else | 
| 1346 | fputc('\n', stderr); | 
| 1347 | } | 
| 1348 | for (i = 0; i < nsbins; i++)    /* send rcontrib ray samples */ | 
| 1349 | if (!(*curparams.sample_basis)(&curparams, i, rcfp)) | 
| 1350 | return(1); | 
| 1351 | return(pclose(rcfp) == 0);      /* all finished! */ | 
| 1352 | userr: | 
| 1353 | if (a < argc-2) | 
| 1354 | fprintf(stderr, "%s: unsupported option '%s'", progname, argv[a]); | 
| 1355 | fprintf(stderr, "Usage: %s [-v][rcontrib options] sender.rad receiver.rad [system.rad ..]\n", | 
| 1356 | progname); | 
| 1357 | return(1); | 
| 1358 | } |