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

Comparing ray/src/cv/pabopto2xml.c (file contents):
Revision 2.7 by greg, Sun Sep 2 15:33:15 2012 UTC vs.
Revision 2.10 by greg, Tue Sep 18 02:50:13 2012 UTC

# Line 16 | Line 16 | static const char RCSid[] = "$Id$";
16   #include <math.h>
17   #include "bsdf.h"
18  
19 + #define DEBUG           1
20 +
21   #ifndef GRIDRES
22 < #define GRIDRES         200             /* max. grid resolution per side */
22 > #define GRIDRES         200             /* grid resolution per side */
23   #endif
24  
25   #define RSCA            2.7             /* radius scaling factor (empirical) */
# Line 51 | Line 53 | typedef struct s_rbfnode {
53          struct s_rbfnode        *next;          /* next in global RBF list */
54          MIGRATION               *ejl;           /* edge list for this vertex */
55          FVECT                   invec;          /* incident vector direction */
56 +        double                  vtotal;         /* volume for normalization */
57          int                     nrbf;           /* number of RBFs */
58          RBFVAL                  rbfa[1];        /* RBF array (extends struct) */
59 < } RBFLIST;                      /* RBF representation of DSF @ 1 incidence */
59 > } RBFNODE;                      /* RBF representation of DSF @ 1 incidence */
60  
61                                  /* our loaded grid for this incident angle */
62 < static double   theta_in_deg, phi_in_deg;
63 < static GRIDVAL  dsf_grid[GRIDRES][GRIDRES];
62 > static double           theta_in_deg, phi_in_deg;
63 > static GRIDVAL          dsf_grid[GRIDRES][GRIDRES];
64  
65 +                                /* all incident angles in-plane so far? */
66 + static int              single_plane_incident = -1;
67 +
68 +                                /* input/output orientations */
69 + static int              input_orient = 0;
70 + static int              output_orient = 0;
71 +
72                                  /* processed incident DSF measurements */
73 < static RBFLIST          *dsf_list = NULL;
73 > static RBFNODE          *dsf_list = NULL;
74  
75 <                                /* edge (linking) matrices */
75 >                                /* RBF-linking matrices (edges) */
76   static MIGRATION        *mig_list = NULL;
77  
78 < #define mtxval(m,i,j)   (m)->mtx[(i)*(m)->rbfv[1]->nrbf+(j)]
79 < #define nextedge(rbf,m) (m)->enxt[(rbf)==(m)->rbfv[1]]
78 >                                /* migration edges drawn in raster fashion */
79 > static MIGRATION        *mig_grid[GRIDRES][GRIDRES];
80  
81 + #define mtx_nrows(m)    ((m)->rbfv[0]->nrbf)
82 + #define mtx_ncols(m)    ((m)->rbfv[1]->nrbf)
83 + #define mtx_ndx(m,i,j)  ((i)*mtx_ncols(m) + (j))
84 + #define is_src(rbf,m)   ((rbf) == (m)->rbfv[0])
85 + #define is_dest(rbf,m)  ((rbf) == (m)->rbfv[1])
86 + #define nextedge(rbf,m) (m)->enxt[is_dest(rbf,m)]
87 + #define opp_rbf(rbf,m)  (m)->rbfv[is_src(rbf,m)]
88 +
89 + #define round(v)        (int)((v) + .5 - ((v) < -.5))
90 +
91 + /* Compute volume associated with Gaussian lobe */
92 + static double
93 + rbf_volume(const RBFVAL *rbfp)
94 + {
95 +        double  rad = R2ANG(rbfp->crad);
96 +
97 +        return((2.*M_PI) * rbfp->peak * rad*rad);
98 + }
99 +
100   /* Compute outgoing vector from grid position */
101   static void
102 < vec_from_pos(FVECT vec, int xpos, int ypos)
102 > ovec_from_pos(FVECT vec, int xpos, int ypos)
103   {
104          double  uv[2];
105          double  r2;
# Line 81 | Line 110 | vec_from_pos(FVECT vec, int xpos, int ypos)
110          vec[0] = vec[1] = sqrt(2. - r2);
111          vec[0] *= uv[0];
112          vec[1] *= uv[1];
113 <        vec[2] = 1. - r2;
113 >        vec[2] = output_orient*(1. - r2);
114   }
115  
116 < /* Compute grid position from normalized outgoing vector */
116 > /* Compute grid position from normalized input/output vector */
117   static void
118   pos_from_vec(int pos[2], const FVECT vec)
119   {
120          double  sq[2];          /* uniform hemispherical projection */
121 <        double  norm = 1./sqrt(1. + vec[2]);
121 >        double  norm = 1./sqrt(1. + fabs(vec[2]));
122  
123          SDdisk2square(sq, vec[0]*norm, vec[1]*norm);
124  
# Line 99 | Line 128 | pos_from_vec(int pos[2], const FVECT vec)
128  
129   /* Evaluate RBF for DSF at the given normalized outgoing direction */
130   static double
131 < eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
131 > eval_rbfrep(const RBFNODE *rp, const FVECT outvec)
132   {
133          double          res = .0;
134          const RBFVAL    *rbfp;
# Line 109 | Line 138 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
138  
139          rbfp = rp->rbfa;
140          for (n = rp->nrbf; n--; rbfp++) {
141 <                vec_from_pos(odir, rbfp->gx, rbfp->gy);
141 >                ovec_from_pos(odir, rbfp->gx, rbfp->gy);
142                  sig2 = R2ANG(rbfp->crad);
143                  sig2 = (DOT(odir,outvec) - 1.) / (sig2*sig2);
144                  if (sig2 > -19.)
# Line 118 | Line 147 | eval_rbfrep(const RBFLIST *rp, const FVECT outvec)
147          return(res);
148   }
149  
150 + /* Insert a new directional scattering function in our global list */
151 + static void
152 + insert_dsf(RBFNODE *newrbf)
153 + {
154 +        RBFNODE         *rbf, *rbf_last;
155 +
156 +                                        /* keep in ascending theta order */
157 +        for (rbf_last = NULL, rbf = dsf_list;
158 +                        single_plane_incident & (rbf != NULL);
159 +                                        rbf_last = rbf, rbf = rbf->next)
160 +                if (input_orient*rbf->invec[2] < input_orient*newrbf->invec[2])
161 +                        break;
162 +        if (rbf_last == NULL) {
163 +                newrbf->next = dsf_list;
164 +                dsf_list = newrbf;
165 +                return;
166 +        }
167 +        newrbf->next = rbf;
168 +        rbf_last->next = newrbf;
169 + }
170 +
171   /* Count up filled nodes and build RBF representation from current grid */
172 < static RBFLIST *
172 > static RBFNODE *
173   make_rbfrep(void)
174   {
175          int     niter = 16;
176          double  lastVar, thisVar = 100.;
177          int     nn;
178 <        RBFLIST *newnode;
178 >        RBFNODE *newnode;
179          int     i, j;
180          
181          nn = 0;                 /* count selected bins */
# Line 133 | Line 183 | make_rbfrep(void)
183              for (j = 0; j < GRIDRES; j++)
184                  nn += dsf_grid[i][j].nval;
185                                  /* allocate RBF array */
186 <        newnode = (RBFLIST *)malloc(sizeof(RBFLIST) + sizeof(RBFVAL)*(nn-1));
186 >        newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
187          if (newnode == NULL) {
188 <                fputs("Out of memory in make_rbfrep\n", stderr);
188 >                fputs("Out of memory in make_rbfrep()\n", stderr);
189                  exit(1);
190          }
191          newnode->next = NULL;
# Line 143 | Line 193 | make_rbfrep(void)
193          newnode->invec[2] = sin(M_PI/180.*theta_in_deg);
194          newnode->invec[0] = cos(M_PI/180.*phi_in_deg)*newnode->invec[2];
195          newnode->invec[1] = sin(M_PI/180.*phi_in_deg)*newnode->invec[2];
196 <        newnode->invec[2] = sqrt(1. - newnode->invec[2]*newnode->invec[2]);
196 >        newnode->invec[2] = input_orient*sqrt(1. - newnode->invec[2]*newnode->invec[2]);
197 >        newnode->vtotal = 0;
198          newnode->nrbf = nn;
199          nn = 0;                 /* fill RBF array */
200          for (i = 0; i < GRIDRES; i++)
# Line 164 | Line 215 | make_rbfrep(void)
215                          if (dsf_grid[i][j].nval) {
216                                  FVECT   odir;
217                                  double  corr;
218 <                                vec_from_pos(odir, i, j);
218 >                                ovec_from_pos(odir, i, j);
219                                  newnode->rbfa[nn++].peak *= corr =
220                                          dsf_grid[i][j].vsum /
221                                                  eval_rbfrep(newnode, odir);
# Line 173 | Line 224 | make_rbfrep(void)
224                          }
225                  lastVar = thisVar;
226                  thisVar = dsum2/(double)nn;
227 <                /*
227 > #ifdef DEBUG
228                  fprintf(stderr, "Avg., RMS error: %.1f%%  %.1f%%\n",
229                                          100.*dsum/(double)nn,
230                                          100.*sqrt(thisVar));
231 <                */
231 > #endif
232          } while (--niter > 0 && lastVar-thisVar > 0.02*lastVar);
233  
234 <        newnode->next = dsf_list;
235 <        return(dsf_list = newnode);
234 >        nn = 0;                 /* compute sum for normalization */
235 >        while (nn < newnode->nrbf)
236 >                newnode->vtotal += rbf_volume(&newnode->rbfa[nn++]);
237 >
238 >        insert_dsf(newnode);
239 >        return(newnode);
240   }
241  
242   /* Load a set of measurements corresponding to a particular incident angle */
243   static int
244 < load_bsdf_meas(const char *fname)
244 > load_pabopto_meas(const char *fname)
245   {
246          FILE    *fp = fopen(fname, "r");
247          int     inp_is_DSF = -1;
248 <        double  theta_out, phi_out, val;
248 >        double  new_phi, theta_out, phi_out, val;
249          char    buf[2048];
250          int     n, c;
251          
# Line 200 | Line 255 | load_bsdf_meas(const char *fname)
255                  return(0);
256          }
257          memset(dsf_grid, 0, sizeof(dsf_grid));
258 + #ifdef DEBUG
259 +        fprintf(stderr, "Loading measurement file '%s'...\n", fname);
260 + #endif
261                                  /* read header information */
262          while ((c = getc(fp)) == '#' || c == EOF) {
263                  if (fgets(buf, sizeof(buf), fp) == NULL) {
# Line 218 | Line 276 | load_bsdf_meas(const char *fname)
276                  }
277                  if (sscanf(buf, "intheta %lf", &theta_in_deg) == 1)
278                          continue;
279 <                if (sscanf(buf, "inphi %lf", &phi_in_deg) == 1)
279 >                if (sscanf(buf, "inphi %lf", &new_phi) == 1)
280                          continue;
281                  if (sscanf(buf, "incident_angle %lf %lf",
282 <                                &theta_in_deg, &phi_in_deg) == 2)
282 >                                &theta_in_deg, &new_phi) == 2)
283                          continue;
284          }
285          if (inp_is_DSF < 0) {
# Line 230 | Line 288 | load_bsdf_meas(const char *fname)
288                  fclose(fp);
289                  return(0);
290          }
291 <        ungetc(c, fp);          /* read actual data */
291 >        if (!input_orient)              /* check input orientation */
292 >                input_orient = 1 - 2*(theta_in_deg > 90.);
293 >        else if (input_orient > 0 ^ theta_in_deg < 90.) {
294 >                fputs("Cannot handle input angles on both sides of surface\n",
295 >                                stderr);
296 >                exit(1);
297 >        }
298 >        if (single_plane_incident > 0)  /* check if still in plane */
299 >                single_plane_incident = (round(new_phi) == round(phi_in_deg));
300 >        else if (single_plane_incident < 0)
301 >                single_plane_incident = 1;
302 >        phi_in_deg = new_phi;
303 >        ungetc(c, fp);                  /* read actual data */
304          while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
305                  FVECT   ovec;
306                  int     pos[2];
307  
308 +                if (!output_orient)     /* check output orientation */
309 +                        output_orient = 1 - 2*(theta_out > 90.);
310 +                else if (output_orient > 0 ^ theta_out < 90.) {
311 +                        fputs("Cannot handle output angles on both sides of surface\n",
312 +                                        stderr);
313 +                        exit(1);
314 +                }
315                  ovec[2] = sin(M_PI/180.*theta_out);
316                  ovec[0] = cos(M_PI/180.*phi_out) * ovec[2];
317                  ovec[1] = sin(M_PI/180.*phi_out) * ovec[2];
# Line 276 | Line 353 | compute_radii(void)
353                  j = (i&1) ? jn : GRIDRES-1-jn;
354                  if (dsf_grid[i][j].nval)        /* find empty grid pos. */
355                          continue;
356 <                vec_from_pos(ovec0, i, j);
356 >                ovec_from_pos(ovec0, i, j);
357                  inear = jnear = -1;             /* find nearest non-empty */
358                  lastang2 = M_PI*M_PI;
359                  for (ii = i-r; ii <= i+r; ii++) {
# Line 287 | Line 364 | compute_radii(void)
364                          if (jj >= GRIDRES) break;
365                          if (!dsf_grid[ii][jj].nval)
366                                  continue;
367 <                        vec_from_pos(ovec1, ii, jj);
367 >                        ovec_from_pos(ovec1, ii, jj);
368                          ang2 = 2. - 2.*DOT(ovec0,ovec1);
369                          if (ang2 >= lastang2)
370                                  continue;
# Line 304 | Line 381 | compute_radii(void)
381                  if (r > dsf_grid[inear][jnear].crad)
382                          dsf_grid[inear][jnear].crad = r;
383                                                  /* next search radius */
384 <                r = ang2*(2.*GRIDRES/M_PI) + 1;
384 >                r = ang2*(2.*GRIDRES/M_PI) + 3;
385              }
386                                                  /* blur radii over hemisphere */
387          memset(fill_grid, 0, sizeof(fill_grid));
# Line 348 | Line 425 | cull_values(void)
425                          continue;
426                  if (!dsf_grid[i][j].crad)
427                          continue;               /* shouldn't happen */
428 <                vec_from_pos(ovec0, i, j);
428 >                ovec_from_pos(ovec0, i, j);
429                  maxang = 2.*R2ANG(dsf_grid[i][j].crad);
430                  if (maxang > ovec0[2])          /* clamp near horizon */
431                          maxang = ovec0[2];
# Line 364 | Line 441 | cull_values(void)
441                                  continue;
442                          if ((ii == i) & (jj == j))
443                                  continue;       /* don't get self-absorbed */
444 <                        vec_from_pos(ovec1, ii, jj);
444 >                        ovec_from_pos(ovec1, ii, jj);
445                          if (2. - 2.*DOT(ovec0,ovec1) >= maxang2)
446                                  continue;
447                                                  /* absorb sum */
# Line 385 | Line 462 | cull_values(void)
462                  }
463   }
464  
465 + /* Compute (and allocate) migration price matrix for optimization */
466 + static float *
467 + price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
468 + {
469 +        float   *pmtx = (float *)malloc(sizeof(float) *
470 +                                        from_rbf->nrbf * to_rbf->nrbf);
471 +        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
472 +        int     i, j;
473  
474 +        if ((pmtx == NULL) | (vto == NULL)) {
475 +                fputs("Out of memory in migration_costs()\n", stderr);
476 +                exit(1);
477 +        }
478 +        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
479 +                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
480 +
481 +        for (i = from_rbf->nrbf; i--; ) {
482 +            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
483 +            FVECT               vfrom;
484 +            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
485 +            for (j = to_rbf->nrbf; j--; )
486 +                pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
487 +                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
488 +        }
489 +        free(vto);
490 +        return(pmtx);
491 + }
492 +
493 + /* Comparison routine needed for sorting price row */
494 + static const float      *price_arr;
495 + static int
496 + msrt_cmp(const void *p1, const void *p2)
497 + {
498 +        float   c1 = price_arr[*(const int *)p1];
499 +        float   c2 = price_arr[*(const int *)p2];
500 +
501 +        if (c1 > c2) return(1);
502 +        if (c1 < c2) return(-1);
503 +        return(0);
504 + }
505 +
506 + /* Compute minimum (optimistic) cost for moving the given source material */
507 + static double
508 + min_cost(double amt2move, const double *avail, const float *price, int n)
509 + {
510 +        static int      *price_sort = NULL;
511 +        static int      n_alloc = 0;
512 +        double          total_cost = 0;
513 +        int             i;
514 +
515 +        if (amt2move <= FTINY)                  /* pre-emptive check */
516 +                return(0.);
517 +        if (n > n_alloc) {                      /* (re)allocate sort array */
518 +                if (n_alloc) free(price_sort);
519 +                price_sort = (int *)malloc(sizeof(int)*n);
520 +                if (price_sort == NULL) {
521 +                        fputs("Out of memory in min_cost()\n", stderr);
522 +                        exit(1);
523 +                }
524 +                n_alloc = n;
525 +        }
526 +        for (i = n; i--; )
527 +                price_sort[i] = i;
528 +        price_arr = price;
529 +        qsort(price_sort, n, sizeof(int), &msrt_cmp);
530 +                                                /* move cheapest first */
531 +        for (i = 0; i < n && amt2move > FTINY; i++) {
532 +                int     d = price_sort[i];
533 +                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
534 +
535 +                total_cost += amt * price[d];
536 +                amt2move -= amt;
537 +        }
538 +        return(total_cost);
539 + }
540 +
541 + /* Take a step in migration by choosing optimal bucket to transfer */
542 + static double
543 + migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
544 + {
545 +        static double   *src_cost = NULL;
546 +        int             n_alloc = 0;
547 +        const double    maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
548 +        double          amt = 0;
549 +        struct {
550 +                int     s, d;   /* source and destination */
551 +                double  price;  /* price estimate per amount moved */
552 +                double  amt;    /* amount we can move */
553 +        } cur, best;
554 +        int             i;
555 +
556 +        if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
557 +                if (n_alloc)
558 +                        free(src_cost);
559 +                src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
560 +                if (src_cost == NULL) {
561 +                        fputs("Out of memory in migration_step()\n", stderr);
562 +                        exit(1);
563 +                }
564 +                n_alloc = mtx_nrows(mig);
565 +        }
566 +        for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
567 +                src_cost[i] = min_cost(src_rem[i], dst_rem,
568 +                                        pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
569 +
570 +                                                /* find best source & dest. */
571 +        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
572 +        for (cur.s = mtx_nrows(mig); cur.s--; ) {
573 +            const float *price = pmtx + cur.s*mtx_ncols(mig);
574 +            double      cost_others = 0;
575 +            if (src_rem[cur.s] <= FTINY)
576 +                    continue;
577 +            cur.d = -1;                         /* examine cheapest dest. */
578 +            for (i = mtx_ncols(mig); i--; )
579 +                if (dst_rem[i] > FTINY &&
580 +                                (cur.d < 0 || price[i] < price[cur.d]))
581 +                        cur.d = i;
582 +            if (cur.d < 0)
583 +                    return(.0);
584 +            if ((cur.price = price[cur.d]) >= best.price)
585 +                    continue;                   /* no point checking further */
586 +            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
587 +                                src_rem[cur.s] : dst_rem[cur.d];
588 +            if (cur.amt > maxamt) cur.amt = maxamt;
589 +            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
590 +            for (i = mtx_nrows(mig); i--; ) {
591 +                if (i == cur.s) continue;
592 +                cost_others += min_cost(src_rem[i], dst_rem, price, mtx_ncols(mig))
593 +                                        - src_cost[i];
594 +            }
595 +            dst_rem[cur.d] += cur.amt;          /* undo trial move */
596 +            cur.price += cost_others/cur.amt;   /* adjust effective price */
597 +            if (cur.price < best.price)         /* are we better than best? */
598 +                    best = cur;
599 +        }
600 +        if ((best.s < 0) | (best.d < 0))
601 +                return(.0);
602 +                                                /* make the actual move */
603 +        mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
604 +        src_rem[best.s] -= best.amt;
605 +        dst_rem[best.d] -= best.amt;
606 +        return(best.amt);
607 + }
608 +
609 + /* Compute (and insert) migration along directed edge */
610 + static MIGRATION *
611 + make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
612 + {
613 +        const double    end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
614 +        float           *pmtx = price_routes(from_rbf, to_rbf);
615 +        MIGRATION       *newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
616 +                                                        sizeof(float) *
617 +                                        (from_rbf->nrbf*to_rbf->nrbf - 1));
618 +        double          *src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
619 +        double          *dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
620 +        double          total_rem = 1.;
621 +        int             i;
622 +
623 +        if ((newmig == NULL) | (src_rem == NULL) | (dst_rem == NULL)) {
624 +                fputs("Out of memory in make_migration()\n", stderr);
625 +                exit(1);
626 +        }
627 + #ifdef DEBUG
628 +        {
629 +                double  theta, phi;
630 +                theta = acos(from_rbf->invec[2])*(180./M_PI);
631 +                phi = atan2(from_rbf->invec[1],from_rbf->invec[0])*(180./M_PI);
632 +                fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
633 +                                round(theta), round(phi));
634 +                theta = acos(to_rbf->invec[2])*(180./M_PI);
635 +                phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
636 +                fprintf(stderr, "(%d,%d)\n", round(theta), round(phi));
637 +        }
638 + #endif
639 +        newmig->next = NULL;
640 +        newmig->rbfv[0] = from_rbf;
641 +        newmig->rbfv[1] = to_rbf;
642 +        newmig->enxt[0] = newmig->enxt[1] = NULL;
643 +        memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
644 +                                                /* starting quantities */
645 +        for (i = from_rbf->nrbf; i--; )
646 +                src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
647 +        for (i = to_rbf->nrbf; i--; )
648 +                dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
649 +                                                /* move a bit at a time */
650 +        while (total_rem > end_thresh)
651 +                total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
652 +
653 +        free(pmtx);                             /* free working arrays */
654 +        free(src_rem);
655 +        free(dst_rem);
656 +        for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
657 +            float       nf = rbf_volume(&from_rbf->rbfa[i]);
658 +            int         j;
659 +            if (nf <= FTINY) continue;
660 +            nf = from_rbf->vtotal / nf;
661 +            for (j = to_rbf->nrbf; j--; )
662 +                newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
663 +        }
664 +                                                /* insert in edge lists */
665 +        newmig->enxt[0] = from_rbf->ejl;
666 +        from_rbf->ejl = newmig;
667 +        newmig->enxt[1] = to_rbf->ejl;
668 +        to_rbf->ejl = newmig;
669 +        newmig->next = mig_list;
670 +        return(mig_list = newmig);
671 + }
672 +
673 + /* Get triangle surface orientation (unnormalized) */
674 + static void
675 + tri_orient(FVECT vres, const FVECT v1, const FVECT v2, const FVECT v3)
676 + {
677 +        FVECT   v2minus1, v3minus2;
678 +
679 +        VSUB(v2minus1, v2, v1);
680 +        VSUB(v3minus2, v3, v2);
681 +        VCROSS(vres, v2minus1, v3minus2);
682 + }
683 +
684 + /* Determine if vertex order is reversed (inward normal) */
685 + static int
686 + is_rev_tri(const FVECT v1, const FVECT v2, const FVECT v3)
687 + {
688 +        FVECT   tor;
689 +
690 +        tri_orient(tor, v1, v2, v3);
691 +
692 +        return(DOT(tor, v2) < 0.);
693 + }
694 +
695 + /* Find vertices completing triangles on either side of the given edge */
696 + static int
697 + get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
698 + {
699 +        const MIGRATION *ej, *ej2;
700 +        RBFNODE         *tv;
701 +
702 +        rbfv[0] = rbfv[1] = NULL;
703 +        for (ej = mig->rbfv[0]->ejl; ej != NULL;
704 +                                ej = nextedge(mig->rbfv[0],ej)) {
705 +                if (ej == mig)
706 +                        continue;
707 +                tv = opp_rbf(mig->rbfv[0],ej);
708 +                for (ej2 = tv->ejl; ej2 != NULL; ej2 = nextedge(tv,ej2))
709 +                        if (opp_rbf(tv,ej2) == mig->rbfv[1]) {
710 +                                rbfv[is_rev_tri(mig->rbfv[0]->invec,
711 +                                                mig->rbfv[1]->invec,
712 +                                                tv->invec)] = tv;
713 +                                break;
714 +                        }
715 +        }
716 +        return((rbfv[0] != NULL) + (rbfv[1] != NULL));
717 + }
718 +
719 + /* Find context hull vertex to complete triangle (oriented call) */
720 + static RBFNODE *
721 + find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
722 + {
723 +        FVECT   vmid, vor;
724 +        RBFNODE *rbf, *rbfbest = NULL;
725 +        double  dprod2, bestdprod2 = 0.5;
726 +
727 +        VADD(vmid, rbf0->invec, rbf1->invec);
728 +        if (normalize(vmid) == 0)
729 +                return(NULL);
730 +                                                /* XXX exhaustive search */
731 +        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
732 +                if ((rbf == rbf0) | (rbf == rbf1))
733 +                        continue;
734 +                tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
735 +                dprod2 = DOT(vor, vmid);
736 +                if (dprod2 <= FTINY)
737 +                        continue;               /* wrong orientation */
738 +                dprod2 *= dprod2 / DOT(vor,vor);
739 +                if (dprod2 > bestdprod2) {      /* more convex than prev? */
740 +                        rbfbest = rbf;
741 +                        bestdprod2 = dprod2;
742 +                }
743 +        }
744 +        return(rbf);
745 + }
746 +
747 + /* Create new migration edge and grow mesh recursively around it */
748 + static void
749 + mesh_from_edge(RBFNODE *rbf0, RBFNODE *rbf1)
750 + {
751 +        MIGRATION       *newej;
752 +        RBFNODE         *tvert[2];
753 +
754 +        if (rbf0 < rbf1)                        /* avoid migration loops */
755 +                newej = make_migration(rbf0, rbf1);
756 +        else
757 +                newej = make_migration(rbf1, rbf0);
758 +                                                /* triangle on either side? */
759 +        get_triangles(tvert, newej);
760 +        if (tvert[0] == NULL) {                 /* recurse on new right edge */
761 +                tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
762 +                if (tvert[0] != NULL) {
763 +                        mesh_from_edge(rbf0, tvert[0]);
764 +                        mesh_from_edge(rbf1, tvert[0]);
765 +                }
766 +        }
767 +        if (tvert[1] == NULL) {                 /* recurse on new left edge */
768 +                tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
769 +                if (tvert[1] != NULL) {
770 +                        mesh_from_edge(rbf0, tvert[1]);
771 +                        mesh_from_edge(rbf1, tvert[1]);
772 +                }
773 +        }
774 + }
775 +
776 + /* Draw edge list into mig_grid array */
777 + static void
778 + draw_edges()
779 + {
780 +        int             nnull = 0, ntot = 0;
781 +        MIGRATION       *ej;
782 +        int             p0[2], p1[2];
783 +
784 +        /* memset(mig_grid, 0, sizeof(mig_grid)); */
785 +        for (ej = mig_list; ej != NULL; ej = ej->next) {
786 +                ++ntot;
787 +                pos_from_vec(p0, ej->rbfv[0]->invec);
788 +                pos_from_vec(p1, ej->rbfv[1]->invec);
789 +                if ((p0[0] == p1[0]) & (p0[1] == p1[1])) {
790 +                        ++nnull;
791 +                        mig_grid[p0[0]][p0[1]] = ej;
792 +                        continue;
793 +                }
794 +                if (abs(p1[0]-p0[0]) > abs(p1[1]-p0[1])) {
795 +                        const int       xstep = 2*(p1[0] > p0[0]) - 1;
796 +                        const double    ystep = (double)((p1[1]-p0[1])*xstep) /
797 +                                                        (double)(p1[0]-p0[0]);
798 +                        int             x;
799 +                        double          y;
800 +                        for (x = p0[0], y = p0[1]+.5; x != p1[0];
801 +                                                x += xstep, y += ystep)
802 +                                mig_grid[x][(int)y] = ej;
803 +                        mig_grid[x][(int)y] = ej;
804 +                } else {
805 +                        const int       ystep = 2*(p1[1] > p0[1]) - 1;
806 +                        const double    xstep = (double)((p1[0]-p0[0])*ystep) /
807 +                                                        (double)(p1[1]-p0[1]);
808 +                        int             y;
809 +                        double          x;
810 +                        for (y = p0[1], x = p0[0]+.5; y != p1[1];
811 +                                                y += ystep, x += xstep)
812 +                                mig_grid[(int)x][y] = ej;
813 +                        mig_grid[(int)x][y] = ej;
814 +                }
815 +        }
816 +        if (nnull)
817 +                fprintf(stderr, "Warning: %d of %d edges are null\n",
818 +                                nnull, ntot);
819 + }
820 +        
821 + /* Build our triangle mesh from recorded RBFs */
822 + static void
823 + build_mesh()
824 + {
825 +        double          best2 = M_PI*M_PI;
826 +        RBFNODE         *rbf, *rbf_near = NULL;
827 +                                                /* check if isotropic */
828 +        if (single_plane_incident) {
829 +                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
830 +                        if (rbf->next != NULL)
831 +                                make_migration(rbf, rbf->next);
832 +                return;
833 +        }
834 +                                                /* find RBF nearest to head */
835 +        if (dsf_list == NULL)
836 +                return;
837 +        for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
838 +                double  dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
839 +                if (dist2 < best2) {
840 +                        rbf_near = rbf;
841 +                        best2 = dist2;
842 +                }
843 +        }
844 +        if (rbf_near == NULL) {
845 +                fputs("Cannot find nearest point for first edge\n", stderr);
846 +                exit(1);
847 +        }
848 +                                                /* build mesh from this edge */
849 +        mesh_from_edge(dsf_list, rbf_near);
850 +                                                /* draw edge list into grid */
851 +        draw_edges();
852 + }
853 +
854 + /* Identify enclosing triangle for this position (flood fill raster check) */
855 + static int
856 + identify_tri(MIGRATION *miga[3], unsigned char vmap[GRIDRES][(GRIDRES+7)/8],
857 +                        int px, int py)
858 + {
859 +        const int       btest = 1<<(py&07);
860 +
861 +        if (vmap[px][py>>3] & btest)            /* already visited here? */
862 +                return(1);
863 +                                                /* else mark it */
864 +        vmap[px][py>>3] |= btest;
865 +
866 +        if (mig_grid[px][py] != NULL) {         /* are we on an edge? */
867 +                int     i;
868 +                for (i = 0; i < 3; i++) {
869 +                        if (miga[i] == mig_grid[px][py])
870 +                                return(1);
871 +                        if (miga[i] != NULL)
872 +                                continue;
873 +                        while (i > 0 && miga[i-1] > mig_grid[px][py]) {
874 +                                miga[i] = miga[i-1];
875 +                                --i;            /* order vertices by pointer */
876 +                        }
877 +                        miga[i] = mig_grid[px][py];
878 +                        return(1);
879 +                }
880 +                return(0);                      /* outside triangle! */
881 +        }
882 +                                                /* check neighbors (flood) */
883 +        if (px > 0 && !identify_tri(miga, vmap, px-1, py))
884 +                return(0);
885 +        if (px < GRIDRES-1 && !identify_tri(miga, vmap, px+1, py))
886 +                return(0);
887 +        if (py > 0 && !identify_tri(miga, vmap, px, py-1))
888 +                return(0);
889 +        if (py < GRIDRES-1 && !identify_tri(miga, vmap, px, py+1))
890 +                return(0);
891 +        return(1);                              /* this neighborhood done */
892 + }
893 +
894 + /* Find edge(s) for interpolating the given incident vector */
895 + static int
896 + get_interp(MIGRATION *miga[3], const FVECT invec)
897 + {
898 +        miga[0] = miga[1] = miga[2] = NULL;
899 +        if (single_plane_incident) {            /* isotropic BSDF? */
900 +                RBFNODE *rbf;                   /* find edge we're on */
901 +                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
902 +                        if (input_orient*rbf->invec[2] < input_orient*invec[2])
903 +                                break;
904 +                        if (rbf->next != NULL &&
905 +                                        input_orient*rbf->next->invec[2] <
906 +                                                        input_orient*invec[2]) {
907 +                                for (miga[0] = rbf->ejl; miga[0] != NULL;
908 +                                                miga[0] = nextedge(rbf,miga[0]))
909 +                                        if (opp_rbf(rbf,miga[0]) == rbf->next)
910 +                                                return(1);
911 +                                break;
912 +                        }
913 +                }
914 +                return(0);                      /* outside range! */
915 +        }
916 +        {                                       /* else use triagnle mesh */
917 +                unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
918 +                int             pstart[2];
919 +
920 +                pos_from_vec(pstart, invec);
921 +                memset(floodmap, 0, sizeof(floodmap));
922 +                                                /* call flooding function */
923 +                if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
924 +                        return(0);              /* outside mesh */
925 +                if ((miga[0] == NULL) | (miga[2] == NULL))
926 +                        return(0);              /* should never happen */
927 +                if (miga[1] == NULL)
928 +                        return(1);              /* on edge */
929 +                return(3);                      /* else in triangle */
930 +        }
931 + }
932 +
933 + /* Advect and allocate new RBF along edge */
934 + static RBFNODE *
935 + e_advect_rbf(const MIGRATION *mig, const FVECT invec)
936 + {
937 +        RBFNODE         *rbf;
938 +        int             n, i, j;
939 +        double          t, full_dist;
940 +                                                /* get relative position */
941 +        t = acos(DOT(invec, mig->rbfv[0]->invec));
942 +        if (t < M_PI/GRIDRES) {                 /* near first DSF */
943 +                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[0]->nrbf-1);
944 +                rbf = (RBFNODE *)malloc(n);
945 +                if (rbf == NULL)
946 +                        goto memerr;
947 +                memcpy(rbf, mig->rbfv[0], n);   /* just duplicate */
948 +                return(rbf);
949 +        }
950 +        full_dist = acos(DOT(mig->rbfv[0]->invec, mig->rbfv[1]->invec));
951 +        if (t > full_dist-M_PI/GRIDRES) {       /* near second DSF */
952 +                n = sizeof(RBFNODE) + sizeof(RBFVAL)*(mig->rbfv[1]->nrbf-1);
953 +                rbf = (RBFNODE *)malloc(n);
954 +                if (rbf == NULL)
955 +                        goto memerr;
956 +                memcpy(rbf, mig->rbfv[1], n);   /* just duplicate */
957 +                return(rbf);
958 +        }
959 +        t /= full_dist;
960 +        n = 0;                                  /* count migrating particles */
961 +        for (i = 0; i < mtx_nrows(mig); i++)
962 +            for (j = 0; j < mtx_ncols(mig); j++)
963 +                n += (mig->mtx[mtx_ndx(mig,i,j)] > FTINY);
964 +        
965 +        rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
966 +        if (rbf == NULL)
967 +                goto memerr;
968 +        rbf->next = NULL; rbf->ejl = NULL;
969 +        VCOPY(rbf->invec, invec);
970 +        rbf->nrbf = n;
971 +        rbf->vtotal = 1.-t + t*mig->rbfv[1]->vtotal/mig->rbfv[0]->vtotal;
972 +        n = 0;                                  /* advect RBF lobes */
973 +        for (i = 0; i < mtx_nrows(mig); i++) {
974 +            const RBFVAL        *rbf0i = &mig->rbfv[0]->rbfa[i];
975 +            const float         peak0 = rbf0i->peak;
976 +            const double        rad0 = R2ANG(rbf0i->crad);
977 +            FVECT               v0;
978 +            float               mv;
979 +            ovec_from_pos(v0, rbf0i->gx, rbf0i->gy);
980 +            for (j = 0; j < mtx_ncols(mig); j++)
981 +                if ((mv = mig->mtx[mtx_ndx(mig,i,j)]) > FTINY) {
982 +                        const RBFVAL    *rbf1j = &mig->rbfv[1]->rbfa[j];
983 +                        double          rad1 = R2ANG(rbf1j->crad);
984 +                        FVECT           v;
985 +                        int             pos[2];
986 +                        rbf->rbfa[n].peak = peak0 * mv * rbf->vtotal;
987 +                        rbf->rbfa[n].crad = ANG2R(sqrt(rad0*rad0*(1.-t) +
988 +                                                        rad1*rad1*t));
989 +                        ovec_from_pos(v, rbf1j->gx, rbf1j->gy);
990 +                        geodesic(v, v0, v, t, GEOD_REL);
991 +                        pos_from_vec(pos, v);
992 +                        rbf->rbfa[n].gx = pos[0];
993 +                        rbf->rbfa[n].gy = pos[1];
994 +                        ++n;
995 +                }
996 +        }
997 +        rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
998 +        return(rbf);
999 + memerr:
1000 +        fputs("Out of memory in e_advect_rbf()\n", stderr);
1001 +        exit(1);
1002 +        return(NULL);   /* pro forma return */
1003 + }
1004 +
1005 + /* Partially advect between recorded incident angles and allocate new RBF */
1006 + static RBFNODE *
1007 + advect_rbf(const FVECT invec)
1008 + {
1009 +        MIGRATION       *miga[3];
1010 +        RBFNODE         *rbf;
1011 +        int             n, i, j;
1012 +        double          s, t;
1013 +
1014 +        if (!get_interp(miga, invec))           /* can't interpolate? */
1015 +                return(NULL);
1016 +        if (miga[1] == NULL)                    /* along edge? */
1017 +                return(e_advect_rbf(miga[0], invec));
1018 +                                                /* figure out position */
1019 +        
1020 +        rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1021 +        if (rbf == NULL) {
1022 +                fputs("Out of memory in advect_rbf()\n", stderr);
1023 +                exit(1);
1024 +        }
1025 +        rbf->next = NULL; rbf->ejl = NULL;
1026 +        VCOPY(rbf->invec, invec);
1027 +        rbf->vtotal = 0;
1028 +        rbf->nrbf = n;
1029 +        
1030 +        return(rbf);
1031 + }
1032 +
1033   #if 1
1034   /* Test main produces a Radiance model from the given input file */
1035   int
# Line 401 | Line 1045 | main(int argc, char *argv[])
1045                  fprintf(stderr, "Usage: %s input.dat > output.rad\n", argv[0]);
1046                  return(1);
1047          }
1048 <        if (!load_bsdf_meas(argv[1]))
1048 >        if (!load_pabopto_meas(argv[1]))
1049                  return(1);
1050  
1051          compute_radii();
# Line 414 | Line 1058 | main(int argc, char *argv[])
1058          for (i = 0; i < GRIDRES; i++)
1059              for (j = 0; j < GRIDRES; j++)
1060                  if (dsf_grid[i][j].vsum > .0f) {
1061 <                        vec_from_pos(dir, i, j);
1061 >                        ovec_from_pos(dir, i, j);
1062                          bsdf = dsf_grid[i][j].vsum / dir[2];
1063                          if (dsf_grid[i][j].nval) {
1064                                  printf("pink cone c%04d\n0\n0\n8\n", ++n);
# Line 425 | Line 1069 | main(int argc, char *argv[])
1069                                          dir[2]*(bsdf+.005));
1070                                  puts("\t.003\t0\n");
1071                          } else {
1072 <                                vec_from_pos(dir, i, j);
1072 >                                ovec_from_pos(dir, i, j);
1073                                  printf("yellow sphere s%04d\n0\n0\n", ++n);
1074                                  printf("4 %.6g %.6g %.6g .0015\n\n",
1075                                          dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);
# Line 443 | Line 1087 | main(int argc, char *argv[])
1087          }
1088          for (i = 0; i < GRIDRES; i++)
1089              for (j = 0; j < GRIDRES; j++) {
1090 <                vec_from_pos(dir, i, j);
1090 >                ovec_from_pos(dir, i, j);
1091                  bsdf = eval_rbfrep(dsf_list, dir) / dir[2];
1092                  fprintf(pfp, "%.8e %.8e %.8e\n",
1093                                  dir[0]*bsdf, dir[1]*bsdf, dir[2]*bsdf);

Diff Legend

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