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.2 by greg, Sun Feb 9 22:19:30 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"
# Line 87 | Line 88 | add2zone(IDF_PARAMETER *param, const char *zname)
88          return(zptr);
89   }
90  
91 < /* Determine if a parameter is a surface in the indicated zone */
92 < static int
93 < in_zone(IDF_PARAMETER *param, const char *zname)
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(0);
101 <                                                /* get zone field */
102 <        fptr = idf_getfield(param, surf_type[i].zone_fld);
103 <        if (fptr == NULL)
104 <                return(0);
105 <                                                /* check for match */
106 <        if (strcmp(fptr->arg, zname))
106 <                return(0);
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(surf_type[i].vert_fld);
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(const char *zname, IDF_PARAMETER *param, FILE *ofp)
113 > rad_surface(IDF_PARAMETER *param, FILE *ofp)
114   {
115 <        const char      *sname = idf_getfield(param, NAME_FLD)->arg;
116 <        IDF_FIELD       *fptr = idf_getfield(param, in_zone(param, zname));
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->arg)) < 3) {
119 >        if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
120                  fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
121                  return(0);
122          }
# Line 125 | Line 125 | rad_surface(const char *zname, IDF_PARAMETER *param, F
125          while (nvert--) {
126                  for (i = 3; i--; ) {
127                          fptr = fptr->next;
128 <                        if (fptr == NULL) {
128 >                        if (fptr == NULL || !isflt(fptr->val)) {
129                                  fprintf(stderr,
130 <                                "%s: missing vertex fields in surface '%s'\n",
131 <                                                progname, sname);
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->arg, ofp);
135 >                        fputs(fptr->val, ofp);
136                  }
137                  fputc('\n', ofp);
138          }
139          return(!ferror(ofp));
140   }
141  
142 < /* Strat rcontrib process */
142 > /* Start rcontrib process */
143   static int
144   start_rcontrib(SUBPROC *pd, ZONE *zp)
145   {
# Line 169 | Line 169 | start_rcontrib(SUBPROC *pd, ZONE *zp)
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) {
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 <                av[i++] = "-m";
179 <                av[i++] = fptr->arg;
180 <                if (!rad_surface(zp->zname, pptr, ofp))
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);
# Line 197 | Line 197 | start_rcontrib(SUBPROC *pd, ZONE *zp)
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 *sp)
356 > compute_uvfs(SUBPROC *pd, ZONE *zp)
357   {
358 <        
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 */
# Line 291 | Line 505 | main(int argc, char *argv[])
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 */

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines