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

Comparing ray/src/util/eplus_adduvf.c (file contents):
Revision 2.1 by greg, Sun Feb 9 05:49:21 2014 UTC vs.
Revision 2.3 by greg, Mon Feb 10 04:51:26 2014 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7   *      G.Ward for LBNL
8   */
9  
10 – #include <stdio.h>
10   #include <stdlib.h>
11 < #include <string.h>
11 > #include "rtio.h"
12 > #include "rtmath.h"
13 > #include "random.h"
14   #include "eplus_idf.h"
15   #include "triangulate.h"
16 + #include "rtprocess.h"
17  
18   #ifndef NSAMPLES
19   #define NSAMPLES        100000                  /* number of samples to use */
20   #endif
21  
22 < #define UVF_PNAME       "ZoneProperty:UserViewFactor:bySurfaceName"
22 > char            *progname;                      /* global argv[0] */
23  
24 + char            temp_octree[128];                       /* temporary octree */
25 +
26 + const char      UVF_PNAME[] =
27 +                        "ZoneProperty:UserViewFactor:bySurfaceName";
28 +
29   const char      ADD_HEADER[] =
30 <                "!+++ User View Factors computed by Radiance +++!\n";
30 >                        "!+++ User View Factors computed by Radiance +++!\n";
31  
32   #define NAME_FLD        1                       /* name field always first? */
33  
# Line 35 | Line 42 | const SURF_PTYPE       surf_type[] = {
42                  {NULL}
43          };
44  
45 < struct s_zone {
45 > typedef struct s_zone {
46          const char      *zname;                 /* zone name */
47          struct s_zone   *next;                  /* next zone in list */
48 +        int             nsurf;                  /* surface count */
49          IDF_PARAMETER   *pfirst;                /* first matching parameter */
50          IDF_PARAMETER   *plast;                 /* last matching parameter */
51 < } *zone_list = NULL;
51 > } ZONE;                 /* a list of collected zone surfaces */
52  
53 + ZONE            *zone_list = NULL;      /* our list of zones */
54 +
55   IDF_LOADED      *our_idf = NULL;        /* loaded/modified IDF */
56  
57   /* Create a new zone and push to top of our list */
58 < static struct s_zone *
58 > static ZONE *
59   new_zone(const char *zname, IDF_PARAMETER *param)
60   {
61 <        struct s_zone   *znew = (struct s_zone *)malloc(sizeof(struct s_zone));
61 >        ZONE    *znew = (ZONE *)malloc(sizeof(ZONE));
62  
63          if (znew == NULL)
64                  return(NULL);
65          znew->zname = zname;                    /* assumes static */
66          znew->next = zone_list;
67          znew->pfirst = znew->plast = param;
68 +        znew->nsurf = 1;
69          return(zone_list = znew);
70   }
71  
72   /* Add the detailed surface (polygon) to the named zone */
73 < static struct s_zone *
73 > static ZONE *
74   add2zone(IDF_PARAMETER *param, const char *zname)
75   {
76 <        struct s_zone   *zptr;
76 >        ZONE    *zptr;
77  
78          for (zptr = zone_list; zptr != NULL; zptr = zptr->next)
79                  if (!strcmp(zptr->zname, zname))
# Line 73 | Line 84 | add2zone(IDF_PARAMETER *param, const char *zname)
84          if (!idf_movparam(our_idf, param, zptr->plast))
85                  return(NULL);
86          zptr->plast = param;
87 +        zptr->nsurf++;
88          return(zptr);
89   }
90  
91 + /* Return field for vertices in the given parameter */
92 + static IDF_FIELD *
93 + get_vlist(IDF_PARAMETER *param, const char *zname)
94 + {
95 +        int             i = 0;
96 +        IDF_FIELD       *fptr;
97 +                                                /* check for surface type */
98 +        while (strcmp(surf_type[i].pname, param->pname))
99 +                if (surf_type[++i].pname == NULL)
100 +                        return(NULL);
101 +
102 +        if (zname != NULL) {                    /* matches specified zone? */
103 +                fptr = idf_getfield(param, surf_type[i].zone_fld);
104 +                if (fptr == NULL || strcmp(fptr->val, zname))
105 +                        return(NULL);
106 +        }
107 +                                                /* return field for #verts */
108 +        return(idf_getfield(param, surf_type[i].vert_fld));
109 + }
110 +
111 + /* Convert surface to Radiance with modifier based on unique name */
112 + static int
113 + rad_surface(IDF_PARAMETER *param, FILE *ofp)
114 + {
115 +        const char      *sname = idf_getfield(param,NAME_FLD)->val;
116 +        IDF_FIELD       *fptr = get_vlist(param, NULL);
117 +        int             nvert, i;
118 +
119 +        if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
120 +                fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
121 +                return(0);
122 +        }
123 +        fprintf(ofp, "\nvoid glow '%s'\n0\n0\n4 1 1 1 0\n", sname);
124 +        fprintf(ofp, "\n'%s' polygon 's_%s'\n0\n0\n%d\n", sname, sname, 3*nvert);
125 +        while (nvert--) {
126 +                for (i = 3; i--; ) {
127 +                        fptr = fptr->next;
128 +                        if (fptr == NULL || !isflt(fptr->val)) {
129 +                                fprintf(stderr,
130 +                                        "%s: missing/bad vertex for %s '%s'\n",
131 +                                                progname, param->pname, sname);
132 +                                return(0);
133 +                        }
134 +                        fputc('\t', ofp);
135 +                        fputs(fptr->val, ofp);
136 +                }
137 +                fputc('\n', ofp);
138 +        }
139 +        return(!ferror(ofp));
140 + }
141 +
142 + /* Start rcontrib process */
143 + static int
144 + start_rcontrib(SUBPROC *pd, ZONE *zp)
145 + {
146 + #define BASE_AC         3
147 +        static char     *base_av[BASE_AC] = {
148 +                                "rcontrib", "-ff", "-h"
149 +                        };
150 +        char            cbuf[300];
151 +        char            **av;
152 +        FILE            *ofp;
153 +        IDF_PARAMETER   *pptr;
154 +        int             i, n;
155 +                                                /* start oconv command */
156 +        sprintf(cbuf, "oconv - > '%s'", temp_octree);
157 +        if ((ofp = popen(cbuf, "w")) == NULL) {
158 +                fputs(progname, stderr);
159 +                fputs(": cannot open oconv process\n", stderr);
160 +                return(0);
161 +        }
162 +                                                /* allocate argument list */
163 +        av = (char **)malloc(sizeof(char *)*(BASE_AC+4+2*zp->nsurf));
164 +        if (av == NULL)
165 +                return(0);
166 +        for (i = 0; i < BASE_AC; i++)
167 +                av[i] = base_av[i];
168 +        sprintf(cbuf, "%d", NSAMPLES);
169 +        av[i++] = "-c";
170 +        av[i++] = cbuf;                         /* add modifier arguments */
171 +        for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
172 +                IDF_FIELD       *fptr = idf_getfield(pptr,NAME_FLD);
173 +                if (fptr == NULL || !fptr->val[0]) {
174 +                        fputs(progname, stderr);
175 +                        fputs(": missing name for surface parameter\n", stderr);
176 +                        return(0);
177 +                }
178 +                if (!rad_surface(pptr, ofp))    /* add surface to octree */
179 +                        return(0);
180 +                av[i++] = "-m";
181 +                av[i++] = fptr->val;
182 +        }
183 +        if (pclose(ofp) != 0) {                 /* finish oconv */
184 +                fputs(progname, stderr);
185 +                fputs(": error running oconv process\n", stderr);
186 +                return(0);
187 +        }
188 +        av[i++] = temp_octree;                  /* add final octree argument */
189 +        av[i] = NULL;
190 +        if (!open_process(pd, av)) {            /* start process */
191 +                fputs(progname, stderr);
192 +                fputs(": cannot start rcontrib process\n", stderr);
193 +                return(0);
194 +        }
195 +        free(av);                               /* all done -- clean up */
196 +        return(1);
197 + #undef BASE_AC
198 + }
199 +
200 + typedef struct {
201 +        FVECT           sdir[3];        /* XYZ unit sampling vectors */
202 +        double          poff;           /* Z-offset for plane of polygon */
203 +        double          area_left;      /* area left to sample */
204 +        int             samp_left;      /* remaining samples */
205 +        int             wd;             /* output file descriptor */
206 + } POLYSAMP;             /* structure for polygon sampling */
207 +
208 + /* Initialize polygon sampling */
209 + static Vert2_list *
210 + init_poly(POLYSAMP *ps, IDF_FIELD *f0, int nv)
211 + {
212 +        IDF_FIELD       *fptr = f0;
213 +        int             i, j;
214 +        FVECT           *vl3, e1, e2, vc;
215 +        Vert2_list      *vl2 = polyAlloc(nv);
216 +
217 +        if (vl2 == NULL)
218 +                return(NULL);
219 +        vl2->p = ps;
220 +                                        /* get 3-D vertices */
221 +        vl3 = (FVECT *)malloc(sizeof(FVECT)*nv);
222 +        if (vl3 == NULL)
223 +                return(NULL);
224 +        for (i = 0; i < nv; i++)
225 +                for (j = 0; j < 3; j++) {
226 +                        if (fptr == NULL) {
227 +                                fputs(progname, stderr);
228 +                                fputs(": bad vertex in init_poly()\n", stderr);
229 +                                return(NULL);
230 +                        }
231 +                        vl3[i][j] = atof(fptr->val);
232 +                }
233 +                                        /* compute area and normal */
234 +        ps->sdir[2][0] = ps->sdir[2][1] = ps->sdir[2][2] = 0;
235 +        VSUB(e1, vl3[1], vl3[0]);
236 +        for (i = 2; i < nv; i++) {
237 +                VSUB(e2, vl3[i], vl3[0]);
238 +                fcross(vc, e1, e2);
239 +                ps->sdir[2][0] += vc[0];
240 +                ps->sdir[2][1] += vc[1];
241 +                ps->sdir[2][2] += vc[2];
242 +                VCOPY(e1, e2);
243 +        }
244 +        ps->area_left = .5 * normalize(ps->sdir[2]);
245 +        if (ps->area_left == .0) {
246 +                fputs(progname, stderr);
247 +                fputs(": degenerate polygon in init_poly()\n", stderr);
248 +                return(0);
249 +        }
250 +                                        /* create X & Y axes */
251 +        VCOPY(ps->sdir[0], e1);
252 +        normalize(ps->sdir[0]);
253 +        fcross(ps->sdir[1], ps->sdir[2], ps->sdir[0]);
254 +                                        /* compute plane offset */
255 +        ps->poff = DOT(vl3[0], ps->sdir[2]) + FTINY;
256 +                                        /* assign 2-D vertices */
257 +        for (i = 0; i < nv; i++) {
258 +                vl2->v[i].mX = DOT(vl3[i], ps->sdir[0]);
259 +                vl2->v[i].mY = DOT(vl3[i], ps->sdir[1]);
260 +        }
261 +        free(vl3);                      /* it's ready! */
262 +        return(vl2);
263 + }
264 +
265 + /* Generate samples on 2-D triangle */
266 + static int
267 + sample_triangle(const Vert2_list *vl2, int a, int b, int c)
268 + {
269 +        POLYSAMP        *ps = (POLYSAMP *)vl2->p;
270 +        float           *samp;
271 +        FVECT           orig;
272 +        FVECT           ab, ac;
273 +        double          area;
274 +        int             i, j, ns;
275 +                                        /* compute sampling axes */
276 +        for (i = 3; i--; ) {
277 +                orig[i] = vl2->v[a].mX*ps->sdir[0][i] +
278 +                                vl2->v[a].mY*ps->sdir[1][i] +
279 +                                ps->poff*ps->sdir[2][i];
280 +                ab[i] = (vl2->v[b].mX - vl2->v[a].mX)*ps->sdir[0][i] +
281 +                                (vl2->v[b].mY - vl2->v[a].mY)*ps->sdir[1][i];
282 +                ac[i] = (vl2->v[c].mX - vl2->v[a].mX)*ps->sdir[0][i] +
283 +                                (vl2->v[c].mY - vl2->v[a].mY)*ps->sdir[1][i];
284 +        }
285 +                                        /* compute number of samples to take */
286 +        area = .5*(vl2->v[a].mX*vl2->v[b].mY - vl2->v[b].mX*vl2->v[a].mY +
287 +                        vl2->v[b].mX*vl2->v[c].mY - vl2->v[c].mX*vl2->v[b].mY +
288 +                        vl2->v[c].mX*vl2->v[a].mY - vl2->v[a].mX*vl2->v[c].mY);
289 +        if (area < .0) {
290 +                fputs(progname, stderr);
291 +                fputs(": negative triangle area in sample_triangle()\n", stderr);
292 +                return(0);
293 +        }
294 +        if (area >= ps->area_left) {
295 +                ns = ps->samp_left;
296 +                ps->area_left = 0;
297 +        } else {
298 +                ns = (ps->samp_left*area/ps->area_left + .5);
299 +                ps->samp_left -= ns;
300 +                ps->area_left -= area;
301 +        }
302 +        if (ns <= 0)                    /* XXX should be error? */
303 +                return(1);
304 +                                        /* buffer sample rays */
305 +        samp = (float *)malloc(sizeof(float)*6*ns);
306 +        if (samp == NULL)
307 +                return(0);
308 +        for (i = ns; i--; ) {
309 +                double  sv[4];
310 +                multisamp(sv, 4, (i+frandom())/(double)ns);
311 +                sv[0] *= sv[1] = sqrt(sv[1]);
312 +                sv[1] = 1. - sv[1];
313 +                for (j = 3; j--; )
314 +                        samp[ns*6 + j] = orig[j] + sv[0]*ab[j] + sv[1]*ac[j];
315 +                sv[3] = sqrt(sv[3]);
316 +                sv[4] *= 2.*PI;
317 +                samp[ns*6 + 3] = tcos(sv[4]) * sv[3];
318 +                samp[ns*6 + 4] = tsin(sv[4]) * sv[3];
319 +                samp[ns*6 + 5] = sqrt(1. - sv[3]*sv[3]);
320 +        }
321 +                                        /* send to our process */
322 +        writebuf(ps->wd, (char *)samp, sizeof(float)*6*ns);
323 +        free(samp);                     /* that's it! */
324 +        return(1);
325 + }
326 +
327 + /* Sample the given surface */
328 + static int
329 + sample_surface(IDF_PARAMETER *param, int wd)
330 + {
331 +        IDF_FIELD       *fptr = get_vlist(param, NULL);
332 +        POLYSAMP        psamp;
333 +        int             nv;
334 +        Vert2_list      *vlist2;
335 +                                        /* set up our polygon sampler */
336 +        if (fptr == NULL || (nv = atoi(fptr->val)) < 3) {
337 +                fprintf(stderr, "%s: error in %s '%s'\n",
338 +                                progname, param->pname,
339 +                                idf_getfield(param,NAME_FLD)->val);
340 +                return(0);
341 +        }
342 +        vlist2 = init_poly(&psamp, fptr->next, nv);
343 +        if (vlist2 == NULL)
344 +                return(0);
345 +        psamp.samp_left = NSAMPLES;     /* assign samples & destination */
346 +        psamp.wd = wd;
347 +                                        /* sample each subtriangle */
348 +        if (!polyTriangulate(vlist2, &sample_triangle))
349 +                return(0);
350 +        polyFree(vlist2);               /* clean up and return */
351 +        return(1);
352 + }
353 +
354 + /* Compute User View Factors using open rcontrib process */
355 + static int
356 + compute_uvfs(SUBPROC *pd, ZONE *zp)
357 + {
358 +        IDF_PARAMETER   *pptr, *pout, *pptr1;
359 +        float           *uvfa;
360 +        char            uvfbuf[24];
361 +        int             n, m;
362 +                                                /* create output parameter */
363 +        pout = idf_newparam(our_idf, UVF_PNAME,
364 +                        "    ! computed by Radiance\n        ", zp->plast);
365 +        if (pout == NULL) {
366 +                fputs(progname, stderr);
367 +                fputs(": cannot create new IDF parameter\n", stderr);
368 +                return(0);
369 +        }
370 +        if (!idf_addfield(pout, zp->zname,
371 +                        "    !- Zone Name\n        ")) {
372 +                fputs(progname, stderr);
373 +                fputs(": cannot add zone name field\n", stderr);
374 +                return(0);
375 +        }
376 +                                                /* allocate read buffer */
377 +        uvfa = (float *)malloc(sizeof(float)*3*zp->nsurf);
378 +        if (uvfa == NULL)
379 +                return(0);
380 +                                                /* UVFs from each surface */
381 +        for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
382 +                                                /* send samples to rcontrib */
383 +                if (!sample_surface(pptr, pd->w))
384 +                        return(0);
385 +                                                /* read results */
386 +                if (readbuf(pd->r, (char *)uvfa, sizeof(float)*3*zp->nsurf) !=
387 +                                sizeof(float)*3*zp->nsurf) {
388 +                        fputs(progname, stderr);
389 +                        fputs(": read error from rcontrib process\n", stderr);
390 +                        return(0);
391 +                }
392 +                                                /* append UVF fields */
393 +                for (m = 0, pptr1 = zp->pfirst;
394 +                                m < zp->nsurf; m++, pptr1 = pptr1->dnext) {
395 +                        if (pptr1 == pptr) {
396 +                                if (uvfa[3*m + 1] > .001)
397 +                                        fprintf(stderr,
398 +                "%s: warning - non-zero self-VF (%.3e) for surface '%s'\n",
399 +                                                progname, uvfa[3*m + 1],
400 +                                                idf_getfield(pptr,NAME_FLD)->val);
401 +                                continue;       /* don't record self-factor */
402 +                        }
403 +                        sprintf(uvfbuf, "%.6f", uvfa[3*m + 1]);
404 +                        if (!idf_addfield(pout,
405 +                                        idf_getfield(pptr,NAME_FLD)->val, NULL) ||
406 +                                !idf_addfield(pout,
407 +                                        idf_getfield(pptr1,NAME_FLD)->val, NULL) ||
408 +                                !idf_addfield(pout, uvfbuf,
409 +                                                (n || m < zp->nsurf-2) ?
410 +                                                        "\n        " : "\n\n")) {
411 +                                fputs(progname, stderr);
412 +                                fputs(": error adding UVF fields\n", stderr);
413 +                                return(0);
414 +                        }
415 +                }
416 +        }
417 +        free(uvfa);                             /* clean up and return */
418 +        return(1);
419 + }
420 +
421   /* Compute zone User View Factors */
422   static int
423   compute_zones(void)
424   {
425 <        fputs("compute_zones() unimplemented\n", stderr);
426 <        return(0);
425 >        ZONE    *zptr;
426 >                                                /* temporary octree name */
427 >        if (temp_filename(temp_octree, sizeof(temp_octree), TEMPLATE) == NULL) {
428 >                fputs(progname, stderr);
429 >                fputs(": cannot create temporary octree\n", stderr);
430 >                return(0);
431 >        }
432 >                                                /* compute each zone */
433 >        for (zptr = zone_list; zptr != NULL; zptr = zptr->next) {
434 >                SUBPROC rcproc;
435 >                                                /* start rcontrib process */
436 >                if (!start_rcontrib(&rcproc, zptr))
437 >                        return(0);
438 >                                                /* compute+add view factors */
439 >                if (!compute_uvfs(&rcproc, zptr))
440 >                        return(0);
441 >                if (close_process(&rcproc) != 0) {
442 >                        fputs(progname, stderr);
443 >                        fputs(": bad return status from rcontrib\n", stderr);
444 >                        return(0);
445 >                }
446 >        }
447 >        unlink(temp_octree);                    /* remove octree file */
448 >        return(1);
449   }
450  
451   /* Load IDF and compute User View Factors */
452   int
453   main(int argc, char *argv[])
454   {
455 +        int             incl_comments = 1;
456          char            *origIDF, *revIDF;
457          IDF_PARAMETER   *pptr;
458          int             i;
459  
460 +        progname = argv[0];
461 +        if (argc > 2 && !strcmp(argv[1], "-c")) {
462 +                incl_comments = -1;             /* output header only */
463 +                ++argv; --argc;
464 +        }
465          if ((argc < 2) | (argc > 3)) {
466                  fputs("Usage: ", stderr);
467 <                fputs(argv[0], stderr);
468 <                fputs(" Model.idf [Revised.idf]\n", stderr);
467 >                fputs(progname, stderr);
468 >                fputs(" [-c] Model.idf [Revised.idf]\n", stderr);
469                  return(1);
470          }
471          origIDF = argv[1];
# Line 103 | Line 473 | main(int argc, char *argv[])
473                                                  /* load Input Data File */
474          our_idf = idf_load(origIDF);
475          if (our_idf == NULL) {
476 <                fputs(argv[0], stderr);
476 >                fputs(progname, stderr);
477                  fputs(": cannot load IDF '", stderr);
478                  fputs(origIDF, stderr);
479                  fputs("'\n", stderr);
# Line 112 | Line 482 | main(int argc, char *argv[])
482                                                  /* remove existing UVFs */
483          if ((pptr = idf_getparam(our_idf, UVF_PNAME)) != NULL) {
484                  IDF_PARAMETER   *pnext;
485 <                fputs(argv[0], stderr);
485 >                fputs(progname, stderr);
486                  fputs(": removing previous User View Factors\n", stderr);
487                  do {
488                          pnext = pptr->pnext;
# Line 131 | Line 501 | main(int argc, char *argv[])
501                          IDF_FIELD       *fptr = idf_getfield(pptr,
502                                                          surf_type[i].zone_fld);
503                          if (fptr == NULL) {
504 <                                fputs("Warning: missing zone field!\n", stderr);
504 >                                fputs(progname, stderr);
505 >                                fputs(": warning - missing zone field\n", stderr);
506                                  continue;
507                          }
508 <                        if (add2zone(pptr, fptr->arg) == NULL)
508 >                        if (add2zone(pptr, fptr->val) == NULL)
509                                  return(1);
510                  }
511                                                  /* run rcontrib on each zone */
512          if (!compute_zones())
513                  return(1);
514                                                  /* write out modified IDF */
515 <        if (!idf_write(our_idf, revIDF, 1)) {
516 <                fputs(argv[0], stderr);
515 >        if (!idf_write(our_idf, revIDF, incl_comments)) {
516 >                fputs(progname, stderr);
517                  fputs(": error writing IDF '", stderr);
518                  fputs(revIDF, stderr);
519                  fputs("'\n", stderr);

Diff Legend

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