--- ray/src/util/rfluxmtx.c 2014/07/22 02:12:48 2.3 +++ ray/src/util/rfluxmtx.c 2014/07/25 16:58:20 2.9 @@ -1,5 +1,5 @@ #ifndef lint -static const char RCSid[] = "$Id: rfluxmtx.c,v 2.3 2014/07/22 02:12:48 greg Exp $"; +static const char RCSid[] = "$Id: rfluxmtx.c,v 2.9 2014/07/25 16:58:20 greg Exp $"; #endif /* * Calculate flux transfer matrix or matrices using rcontrib @@ -53,7 +53,8 @@ char *kfullfn = "klems_full.cal"; char *khalffn = "klems_half.cal"; char *kquarterfn = "klems_quarter.cal"; -#define PARAMSTART "@rfluxmtx" /* string indicating parameters */ + /* string indicating parameters */ +const char PARAMSTART[] = "@rfluxmtx"; /* surface type IDs */ #define ST_NONE 0 @@ -66,7 +67,7 @@ typedef struct surf_s { void *priv; /* private data (malloc'ed) */ char sname[32]; /* surface name */ FVECT snrm; /* surface normal */ - double area; /* surface area (or solid angle) */ + double area; /* surface area / proj. solid angle */ short styp; /* surface type */ short nfargs; /* number of real arguments */ double farg[1]; /* real values (extends struct) */ @@ -95,6 +96,7 @@ typedef struct param_s { PARAMS curparams; char curmod[128]; +char newparams[1024]; typedef int SURFSAMP(FVECT, SURF *, double); @@ -108,7 +110,6 @@ SURFSAMP *orig_in_surf[4] = { static void clear_params(PARAMS *p, int reset_only) { - curmod[0] = '\0'; while (p->slist != NULL) { SURF *sdel = p->slist; p->slist = sdel->next; @@ -123,7 +124,7 @@ clear_params(PARAMS *p, int reset_only) p->outfn = NULL; return; } - memset(p, 0, sizeof(curparams)); + memset(p, 0, sizeof(PARAMS)); } /* Get surface type from name */ @@ -308,7 +309,7 @@ nextchar: /* Parse program parameters (directives) */ static int -parse_params(char *pargs) +parse_params(PARAMS *p, char *pargs) { char *cp = pargs; int nparams = 0; @@ -319,24 +320,23 @@ parse_params(char *pargs) case 'h': if (*cp++ != '=') break; - curparams.hsiz = 0; + p->hsiz = 0; i = 0; while (*cp && !isspace(*cp)) { if (isdigit(*cp)) - curparams.hsiz = 10*curparams.hsiz + - *cp - '0'; - curparams.hemis[i++] = *cp++; + p->hsiz = 10*p->hsiz + *cp - '0'; + p->hemis[i++] = *cp++; } if (!i) break; - curparams.hemis[i] = '\0'; - curparams.hsiz += !curparams.hsiz; + p->hemis[i] = '\0'; + p->hsiz += !p->hsiz; ++nparams; continue; case 'u': if (*cp++ != '=') break; - if (!get_direction(curparams.vup, cp)) + if (!get_direction(p->vup, cp)) break; ++nparams; continue; @@ -349,7 +349,7 @@ parse_params(char *pargs) if (!i) break; *--cp = '\0'; - curparams.outfn = savqstr(cp-i); + p->outfn = savqstr(cp-i); *cp++ = ' '; ++nparams; continue; @@ -372,12 +372,13 @@ parse_params(char *pargs) static void finish_receiver(void) { + char sbuf[256]; + int uniform = 0; char *calfn = NULL; char *params = NULL; char *binv = NULL; char *binf = NULL; char *nbins = NULL; - char sbuf[256]; if (!curmod[0]) { fputs(progname, stderr); @@ -389,7 +390,7 @@ finish_receiver(void) rcarg[nrcargs++] = "-o"; rcarg[nrcargs++] = curparams.outfn; } - /* add bin specification */ + /* check arguments */ if (!curparams.hemis[0]) { fputs(progname, stderr); fputs(": missing hemisphere sampling type!\n", stderr); @@ -405,8 +406,10 @@ finish_receiver(void) curparams.vup[2] = 1; else curparams.vup[1] = 1; + /* determine sample type/bin */ if (tolower(curparams.hemis[0]) == 'u' | curparams.hemis[0] == '1') { binv = "0"; /* uniform sampling -- one bin */ + uniform = 1; } else if (tolower(curparams.hemis[0]) == 's' && tolower(curparams.hemis[1]) == 'c') { /* assign parameters */ @@ -416,7 +419,7 @@ finish_receiver(void) exit(1); } calfn = shirchiufn; shirchiufn = NULL; - sprintf(sbuf, "SCdim=%d,Nx=%g,Ny=%g,Nz=%g,Ux=%g,Uy=%g,Uz=%g", + sprintf(sbuf, "SCdim=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g", curparams.hsiz, curparams.nrm[0], curparams.nrm[1], curparams.nrm[2], curparams.vup[0], curparams.vup[1], curparams.vup[2]); @@ -426,7 +429,7 @@ finish_receiver(void) } else if ((tolower(curparams.hemis[0]) == 'r') | (tolower(curparams.hemis[0]) == 't')) { calfn = reinhfn; reinhfn = NULL; - sprintf(sbuf, "MF=%d,Nx=%g,Ny=%g,Nz=%g,Ux=%g,Uy=%g,Uz=%g", + sprintf(sbuf, "MF=%d,rNx=%g,rNy=%g,rNz=%g,Ux=%g,Uy=%g,Uz=%g", curparams.hsiz, curparams.nrm[0], curparams.nrm[1], curparams.nrm[2], curparams.vup[0], curparams.vup[1], curparams.vup[2]); @@ -457,6 +460,15 @@ finish_receiver(void) progname, curparams.hemis); exit(1); } + if (!uniform & (curparams.slist->styp == ST_SOURCE)) { + SURF *sp; + for (sp = curparams.slist; sp != NULL; sp = sp->next) + if (fabs(sp->area - PI) > 1e-3) { + fprintf(stderr, "%s: source '%s' must be 180-degrees\n", + progname, sp->sname); + exit(1); + } + } if (calfn != NULL) { /* add cal file if needed */ CHECKARGC(2); rcarg[nrcargs++] = "-f"; @@ -541,7 +553,7 @@ ssamp_ring(FVECT orig, SURF *sp, double x) sp->priv = (void *)uva; } SDmultiSamp(samp2, 2, x); - samp2[0] = sp->farg[6] + sqrt(samp2[0]*sp->area*(1./PI)); + samp2[0] = sqrt(samp2[0]*sp->area*(1./PI) + sp->farg[6]*sp->farg[6]); samp2[1] *= 2.*PI; uv[0] = samp2[0]*tcos(samp2[1]); uv[1] = samp2[0]*tsin(samp2[1]); @@ -594,8 +606,12 @@ ssamp_poly(FVECT orig, SURF *sp, double x) } ptp->ntris = 0; v2l->p = (void *)ptp; - if (!polyTriangulate(v2l, add_triangle)) + if (!polyTriangulate(v2l, add_triangle)) { + fprintf(stderr, + "%s: cannot triangulate polygon '%s'\n", + progname, sp->sname); return(0); + } for (i = ptp->ntris; i--; ) { int a = ptp->tri[i].vndx[0]; int b = ptp->tri[i].vndx[1]; @@ -644,8 +660,12 @@ sample_origin(PARAMS *p, FVECT orig, const FVECT rdir, /* special case for lone surface */ if (p->nsurfs == 1) { sp = p->slist; - if (DOT(sp->snrm, rdir) >= -FTINY) - return(0); /* behind surface! */ + if (DOT(sp->snrm, rdir) >= -FTINY) { + fprintf(stderr, + "%s: internal - sample behind sender '%s'\n", + progname, sp->sname); + return(0); + } return((*orig_in_surf[sp->styp])(orig, sp, x)); } if (p->nsurfs > nall) { /* (re)allocate surface area cache */ @@ -659,8 +679,12 @@ sample_origin(PARAMS *p, FVECT orig, const FVECT rdir, projsa[i] = -DOT(sp->snrm, rdir) * sp->area; tarea += projsa[i] *= (double)(projsa[i] > FTINY); } - if (tarea <= FTINY) /* wrong side of sender? */ + if (tarea <= FTINY) { /* wrong side of sender? */ + fputs(progname, stderr); + fputs(": internal - sample behind all sender elements!\n", + stderr); return(0); + } tarea *= x; /* get surface from list */ for (i = 0, sp = p->slist; tarea > projsa[i]; sp = sp->next) tarea -= projsa[i++]; @@ -733,7 +757,7 @@ sample_reinhart(PARAMS *p, int b, FILE *fp) #define T_NALT 7 static const int tnaz[T_NALT] = {30, 30, 24, 24, 18, 12, 6}; const int RowMax = T_NALT*p->hsiz + 1; - const double RAH = (.25*PI)/(RowMax-.5); + const double RAH = (.5*PI)/(RowMax-.5); #define rnaz(r) (r >= RowMax-1 ? 1 : p->hsiz*tnaz[r/p->hsiz]) int n = sampcnt; int row, col; @@ -758,7 +782,7 @@ sample_reinhart(PARAMS *p, int b, FILE *fp) SDmultiSamp(samp3, 3, (n+frandom())/sampcnt); alt = (row+samp3[1])*RAH; azi = (2.*PI)*(col+samp3[2]-.5)/rnaz(row); - duvw[2] = tcos(alt); /* measured from horizon */ + duvw[2] = cos(alt); /* measured from horizon */ duvw[0] = tcos(azi)*duvw[2]; duvw[1] = tsin(azi)*duvw[2]; duvw[2] = sqrt(1. - duvw[2]*duvw[2]); @@ -976,10 +1000,12 @@ add_surface(int st, const char *oname, FILE *fp) case ST_SOURCE: if (snew->nfargs != 4) goto badcount; - VCOPY(snew->snrm, snew->farg); + for (n = 3; n--; ) /* need to reverse "normal" */ + snew->snrm[n] = -snew->farg[n]; if (normalize(snew->snrm) == 0) goto badnorm; - snew->area = 2.*PI*(1. - cos((PI/180./2.)*snew->farg[3])); + snew->area = sin((PI/180./2.)*snew->farg[3]); + snew->area *= PI*snew->area; break; } if (snew->area <= FTINY) { @@ -994,11 +1020,11 @@ add_surface(int st, const char *oname, FILE *fp) curparams.nsurfs++; return; badcount: - fprintf(stderr, "%s: bad argument count for surface '%s'\n", + fprintf(stderr, "%s: bad argument count for surface element '%s'\n", progname, oname); exit(1); badnorm: - fprintf(stderr, "%s: bad orientation for surface '%s'\n", + fprintf(stderr, "%s: bad orientation for surface element '%s'\n", progname, oname); exit(1); } @@ -1030,13 +1056,15 @@ add_recv_object(FILE *fp) finish_receiver(); clear_params(&curparams, 1); } + parse_params(&curparams, newparams); + newparams[0] = '\0'; strcpy(curmod, thismod); } add_surface(st, oname, fp); /* read & store surface */ return(1); } /* else skip arguments */ - if (!fscanf(fp, "%d", &n)) return; + if (!fscanf(fp, "%d", &n)) return(0); while (n-- > 0) fscanf(fp, "%*s"); if (!fscanf(fp, "%d", &n)) return; while (n-- > 0) fscanf(fp, "%*d"); @@ -1066,11 +1094,13 @@ add_send_object(FILE *fp) fputs(": cannot use source as a sender!\n", stderr); return(-1); } + parse_params(&curparams, newparams); + newparams[0] = '\0'; add_surface(st, oname, fp); /* read & store surface */ return(0); } /* else skip arguments */ - if (!fscanf(fp, "%d", &n)) return; + if (!fscanf(fp, "%d", &n)) return(0); while (n-- > 0) fscanf(fp, "%*s"); if (!fscanf(fp, "%d", &n)) return; while (n-- > 0) fscanf(fp, "%*d"); @@ -1109,10 +1139,12 @@ load_scene(const char *inspec, int (*ocb)(FILE *)) continue; } if (c == '#') { /* parameters/comment */ - if (fscanf(fp, "%s", inpbuf) == 1 && + if ((c = getc(fp)) == EOF || ungetc(c,fp) == EOF) + break; + if (!isspace(c) && fscanf(fp, "%s", inpbuf) == 1 && !strcmp(inpbuf, PARAMSTART)) { if (fgets(inpbuf, sizeof(inpbuf), fp) != NULL) - parse_params(inpbuf); + strcat(newparams, inpbuf); continue; } while ((c = getc(fp)) != EOF && c != '\n'); @@ -1139,6 +1171,8 @@ int main(int argc, char *argv[]) { char fmtopt[6] = "-faa"; /* default output is ASCII */ + char *xrs=NULL, *yrs=NULL, *ldopt=NULL; + int wantIrradiance = 0; char *sendfn; char sampcntbuf[32], nsbinbuf[32]; FILE *rcfp; @@ -1163,26 +1197,39 @@ main(int argc, char *argv[]) case 'f': case 'd': case 'c': - if (!(fmtopt[4] = argv[a][3])) - fmtopt[4] = argv[a][2]; - fmtopt[3] = argv[a][2]; + if (!(fmtopt[3] = argv[a][3])) + fmtopt[3] = argv[a][2]; + fmtopt[2] = argv[a][2]; na = 0; continue; /* will pass later */ default: goto userr; } break; + case 'x': /* x-resolution */ + xrs = argv[++a]; + na = 0; + continue; + case 'y': /* y-resolution */ + yrs = argv[++a]; + na = 0; + continue; case 'c': /* number of samples */ - sampcnt = atoi(argv[a+1]); + sampcnt = atoi(argv[++a]); if (sampcnt <= 0) goto userr; na = 0; /* we re-add this later */ continue; + case 'I': /* only for pass-through mode */ + wantIrradiance = 1; + na = 0; + continue; case 'V': /* options without arguments */ case 'w': case 'u': case 'i': case 'h': + case 'r': break; case 'n': /* options with 1 argument */ case 's': @@ -1193,7 +1240,11 @@ main(int argc, char *argv[]) if (argv[a][2] != 'v') goto userr; break; case 'l': /* special case */ - if (argv[a][2] == 'd') goto userr; + if (argv[a][2] == 'd') { + ldopt = argv[a]; + na = 0; + continue; + } na = 2; break; case 'd': /* special case */ @@ -1224,9 +1275,32 @@ done_opts: if (sendfn[0] == '-') { /* user wants pass-through mode? */ if (sendfn[1]) goto userr; sendfn = NULL; + if (wantIrradiance) { + CHECKARGC(1); + rcarg[nrcargs++] = "-I"; + } + if (xrs) { + CHECKARGC(2); + rcarg[nrcargs++] = "-x"; + rcarg[nrcargs++] = xrs; + } + if (yrs) { + CHECKARGC(2); + rcarg[nrcargs++] = "-y"; + rcarg[nrcargs++] = yrs; + } + if (ldopt) { + CHECKARGC(1); + rcarg[nrcargs++] = ldopt; + } if (sampcnt <= 0) sampcnt = 1; - } else { /* else FVECT determines input format */ - fmtopt[3] = (sizeof(RREAL)==sizeof(double)) ? 'd' : 'f'; + } else { /* else in sampling mode */ + if (wantIrradiance) { + fputs(progname, stderr); + fputs(": -I supported for pass-through only\n", stderr); + return(1); + } + fmtopt[2] = (sizeof(RREAL)==sizeof(double)) ? 'd' : 'f'; if (sampcnt <= 0) sampcnt = 10000; } sprintf(sampcntbuf, "%d", sampcnt); @@ -1265,7 +1339,7 @@ done_opts: if (verbose) { fprintf(stderr, "%s: sampling %d directions", progname, nsbins); if (curparams.nsurfs > 1) - fprintf(stderr, " (%d surfaces)\n", curparams.nsurfs); + fprintf(stderr, " (%d elements)\n", curparams.nsurfs); else fputc('\n', stderr); } @@ -1276,7 +1350,7 @@ done_opts: userr: if (a < argc-2) fprintf(stderr, "%s: unsupported option '%s'", progname, argv[a]); - fprintf(stderr, "Usage: %s [rcontrib options] sender.rad receiver.rad [system.rad ..]\n", + fprintf(stderr, "Usage: %s [-v][rcontrib options] sender.rad receiver.rad [system.rad ..]\n", progname); return(1); }