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.18 by greg, Wed Oct 17 19:01:47 2012 UTC vs.
Revision 2.19 by greg, Thu Oct 18 04:28:20 2012 UTC

# Line 72 | Line 72 | static GRIDVAL         dsf_grid[GRIDRES][GRIDRES];
72                                  /* all incident angles in-plane so far? */
73   static int              single_plane_incident = -1;
74  
75 +                                /* represented incident quadrants */
76 + #define INP_QUAD1       1               /* 0-90 degree quadrant */
77 + #define INP_QUAD2       2               /* 90-180 degree quadrant */
78 + #define INP_QUAD3       4               /* 180-270 degree quadrant */
79 + #define INP_QUAD4       8               /* 270-360 degree quadrant */
80 +
81 + static int              inp_coverage = 0;
82 +
83 +                                /* symmetry operations */
84 + #define MIRROR_X        1               /* mirror(ed) x-coordinate */
85 + #define MIRROR_Y        2               /* mirror(ed) y-coordinate */
86 +
87                                  /* input/output orientations */
88   static int              input_orient = 0;
89   static int              output_orient = 0;
# Line 108 | Line 120 | int                    nchild = 0;
120                                  /* sampling order (set by data density) */
121   int                     samp_order = 0;
122  
123 +                                /* get phi value in degrees, [0,360) range */
124 + #define get_phi360(v)   ((180./M_PI)*atan2((v)[1],(v)[0]) + 180.)
125 +
126 + /* Apply symmetry to the given vector based on distribution */
127 + static int
128 + use_symmetry(FVECT vec)
129 + {
130 +        double  phi = get_phi360(vec);
131 +
132 +        switch (inp_coverage) {
133 +        case INP_QUAD1|INP_QUAD2|INP_QUAD3|INP_QUAD4:
134 +                break;
135 +        case INP_QUAD1|INP_QUAD2:
136 +                if ((-FTINY > phi) | (phi > 180.+FTINY))
137 +                        goto mir_y;
138 +                break;
139 +        case INP_QUAD2|INP_QUAD3:
140 +                if ((90.-FTINY > phi) | (phi > 270.+FTINY))
141 +                        goto mir_x;
142 +                break;
143 +        case INP_QUAD3|INP_QUAD4:
144 +                if ((180.-FTINY > phi) | (phi > 360.+FTINY))
145 +                        goto mir_y;
146 +                break;
147 +        case INP_QUAD4|INP_QUAD1:
148 +                if ((270.-FTINY > phi) & (phi > 90.+FTINY))
149 +                        goto mir_x;
150 +                break;
151 +        case INP_QUAD1:
152 +                if ((-FTINY > phi) | (phi > 90.+FTINY))
153 +                        switch ((int)(phi*(1./90.))) {
154 +                        case 1: goto mir_x;
155 +                        case 2: goto mir_xy;
156 +                        case 3: goto mir_y;
157 +                        }
158 +                break;
159 +        case INP_QUAD2:
160 +                if ((90.-FTINY > phi) | (phi > 180.+FTINY))
161 +                        switch ((int)(phi*(1./90.))) {
162 +                        case 0: goto mir_x;
163 +                        case 2: goto mir_y;
164 +                        case 3: goto mir_xy;
165 +                        }
166 +                break;
167 +        case INP_QUAD3:
168 +                if ((180.-FTINY > phi) | (phi > 270.+FTINY))
169 +                        switch ((int)(phi*(1./90.))) {
170 +                        case 0: goto mir_xy;
171 +                        case 1: goto mir_y;
172 +                        case 3: goto mir_x;
173 +                        }
174 +                break;
175 +        case INP_QUAD4:
176 +                if ((270.-FTINY > phi) | (phi > 360.+FTINY))
177 +                        switch ((int)(phi*(1./90.))) {
178 +                        case 0: goto mir_y;
179 +                        case 1: goto mir_xy;
180 +                        case 2: goto mir_x;
181 +                        }
182 +                break;
183 +        default:
184 +                fprintf(stderr, "%s: Illegal input coverage (%d)\n",
185 +                                        progname, inp_coverage);
186 +                exit(1);
187 +        }
188 +        return(0);              /* in range */
189 + mir_x:
190 +        vec[0] = -vec[0];
191 +        return(MIRROR_X);
192 + mir_y:
193 +        vec[1] = -vec[1];
194 +        return(MIRROR_Y);
195 + mir_xy:
196 +        vec[0] = -vec[0];
197 +        vec[1] = -vec[1];
198 +        return(MIRROR_X|MIRROR_Y);
199 + }
200 +
201 + /* Reverse symmetry based on what was done before */
202 + static void
203 + rev_symmetry(FVECT vec, int sym)
204 + {
205 +        if (sym & MIRROR_X)
206 +                vec[0] = -vec[0];
207 +        if (sym & MIRROR_Y)
208 +                vec[1] = -vec[1];
209 + }
210 +
211 + /* Reverse symmetry for an RBF distribution */
212 + static void
213 + rev_rbf_symmetry(RBFNODE *rbf, int sym)
214 + {
215 +        int     n;
216 +
217 +        rev_symmetry(rbf->invec, sym);
218 +        if (sym & MIRROR_X)
219 +                for (n = rbf->nrbf; n-- > 0; )
220 +                        rbf->rbfa[n].gx = GRIDRES-1 - rbf->rbfa[n].gx;
221 +        if (sym & MIRROR_Y)
222 +                for (n = rbf->nrbf; n-- > 0; )
223 +                        rbf->rbfa[n].gy = GRIDRES-1 - rbf->rbfa[n].gy;
224 + }
225 +
226   /* Compute volume associated with Gaussian lobe */
227   static double
228   rbf_volume(const RBFVAL *rbfp)
# Line 177 | Line 292 | insert_dsf(RBFNODE *newrbf)
292                                          /* check for redundant meas. */
293          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
294                  if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
295 <                        fputs("Duplicate incident measurement (ignored)\n", stderr);
295 >                        fprintf(stderr,
296 >                                "%s: Duplicate incident measurement (ignored)\n",
297 >                                        progname);
298                          free(newrbf);
299                          return;
300                  }
# Line 214 | Line 331 | make_rbfrep(void)
331                                  /* allocate RBF array */
332          newnode = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(nn-1));
333          if (newnode == NULL) {
334 <                fputs("Out of memory in make_rbfrep()\n", stderr);
334 >                fprintf(stderr, "%s: Out of memory in make_rbfrep()\n", progname);
335                  exit(1);
336          }
337          newnode->next = NULL;
# Line 331 | Line 448 | load_pabopto_meas(const char *fname)
448                                  stderr);
449                  exit(1);
450          }
451 <        if (single_plane_incident > 0)  /* check if still in plane */
451 >        if (single_plane_incident > 0)  /* check input coverage */
452                  single_plane_incident = (round(new_phi) == round(phi_in_deg));
453          else if (single_plane_incident < 0)
454                  single_plane_incident = 1;
455          phi_in_deg = new_phi;
456 +        new_phi += 360.*(new_phi < -FTINY);
457 +        if ((1. < new_phi) & (new_phi < 89.))
458 +                inp_coverage |= INP_QUAD1;
459 +        else if ((91. < new_phi) & (new_phi < 179.))
460 +                inp_coverage |= INP_QUAD2;
461 +        else if ((181. < new_phi) & (new_phi < 269.))
462 +                inp_coverage |= INP_QUAD3;
463 +        else if ((271. < new_phi) & (new_phi < 359.))
464 +                inp_coverage |= INP_QUAD4;
465          ungetc(c, fp);                  /* read actual data */
466          while (fscanf(fp, "%lf %lf %lf\n", &theta_out, &phi_out, &val) == 3) {
467                  FVECT   ovec;
# Line 409 | Line 535 | compute_radii(void)
535                      }
536                  }
537                  if (inear < 0) {
538 <                        fputs("Could not find non-empty neighbor!\n", stderr);
538 >                        fprintf(stderr,
539 >                                "%s: Could not find non-empty neighbor!\n",
540 >                                        progname);
541                          exit(1);
542                  }
543                  ang2 = sqrt(lastang2);
# Line 508 | Line 636 | price_routes(const RBFNODE *from_rbf, const RBFNODE *t
636          int     i, j;
637  
638          if ((pmtx == NULL) | (vto == NULL)) {
639 <                fputs("Out of memory in migration_costs()\n", stderr);
639 >                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
640 >                                progname);
641                  exit(1);
642          }
643          for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
# Line 554 | Line 683 | min_cost(double amt2move, const double *avail, const f
683                  if (n_alloc) free(price_sort);
684                  price_sort = (int *)malloc(sizeof(int)*n);
685                  if (price_sort == NULL) {
686 <                        fputs("Out of memory in min_cost()\n", stderr);
686 >                        fprintf(stderr, "%s: Out of memory in min_cost()\n",
687 >                                        progname);
688                          exit(1);
689                  }
690                  n_alloc = n;
# Line 579 | Line 709 | static double
709   migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
710   {
711          const double    maxamt = .1;
712 +        const double    minamt = maxamt*.0001;
713          static double   *src_cost = NULL;
714          static int      n_alloc = 0;
715          struct {
# Line 593 | Line 724 | migration_step(MIGRATION *mig, double *src_rem, double
724                          free(src_cost);
725                  src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
726                  if (src_cost == NULL) {
727 <                        fputs("Out of memory in migration_step()\n", stderr);
727 >                        fprintf(stderr, "%s: Out of memory in migration_step()\n",
728 >                                        progname);
729                          exit(1);
730                  }
731                  n_alloc = mtx_nrows(mig);
# Line 607 | Line 739 | migration_step(MIGRATION *mig, double *src_rem, double
739          for (cur.s = mtx_nrows(mig); cur.s--; ) {
740              const float *price = pmtx + cur.s*mtx_ncols(mig);
741              double      cost_others = 0;
742 <            if (src_rem[cur.s] <= FTINY)
742 >            if (src_rem[cur.s] < minamt)
743                      continue;
744              cur.d = -1;                         /* examine cheapest dest. */
745              for (i = mtx_ncols(mig); i--; )
746 <                if (dst_rem[i] > FTINY &&
746 >                if (dst_rem[i] > minamt &&
747                                  (cur.d < 0 || price[i] < price[cur.d]))
748                          cur.d = i;
749              if (cur.d < 0)
# Line 758 | Line 890 | static MIGRATION *
890   create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
891   {
892          const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
893 <        const double    rel_thresh = 0.0001;
893 >        const double    check_thresh = 0.01;
894 >        const double    rel_thresh = 5e-6;
895          float           *pmtx;
896          MIGRATION       *newmig;
897          double          *src_rem, *dst_rem;
# Line 777 | Line 910 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
910          src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
911          dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
912          if ((src_rem == NULL) | (dst_rem == NULL)) {
913 <                fputs("Out of memory in create_migration()\n", stderr);
913 >                fprintf(stderr, "%s: Out of memory in create_migration()\n",
914 >                                progname);
915                  exit(1);
916          }
917   #ifdef DEBUG
# Line 800 | Line 934 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
934                          /* fputc('.', stderr); */
935                          fprintf(stderr, "%.9f remaining...\r", total_rem);
936   #endif
937 <        } while ((total_rem > end_thresh) & (move_amt > rel_thresh*total_rem));
937 >        } while (total_rem > end_thresh && (total_rem > check_thresh) |
938 >                                        (move_amt > rel_thresh*total_rem));
939   #ifdef DEBUG
940          if (!nchild) fputs("\ndone.\n", stderr);
941 +        else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
942   #endif
943          for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
944              float       nf = rbf_volume(&from_rbf->rbfa[i]);
# Line 1082 | Line 1218 | build_mesh()
1218                  await_children(nchild);
1219                  return;
1220          }
1221 <                                                /* start w/ shortest edge */
1086 <        shrt_edj[0] = shrt_edj[1] = NULL;
1221 >        shrt_edj[0] = shrt_edj[1] = NULL;       /* start w/ shortest edge */
1222          for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
1223              for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
1224                  double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
# Line 1094 | Line 1229 | build_mesh()
1229                  }
1230          }
1231          if (shrt_edj[0] == NULL) {
1232 <                fputs("Cannot find shortest edge\n", stderr);
1232 >                fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
1233                  exit(1);
1234          }
1235                                                  /* build mesh from this edge */
# Line 1215 | Line 1350 | order_triangle(MIGRATION *miga[3])
1350          return(1);
1351   }
1352  
1353 < /* Find edge(s) for interpolating the given incident vector */
1353 > /* Find edge(s) for interpolating the given vector, applying symmetry */
1354   static int
1355 < get_interp(MIGRATION *miga[3], const FVECT invec)
1355 > get_interp(MIGRATION *miga[3], FVECT invec)
1356   {
1357          miga[0] = miga[1] = miga[2] = NULL;
1358          if (single_plane_incident) {            /* isotropic BSDF? */
# Line 1231 | Line 1366 | get_interp(MIGRATION *miga[3], const FVECT invec)
1366                                  for (miga[0] = rbf->ejl; miga[0] != NULL;
1367                                                  miga[0] = nextedge(rbf,miga[0]))
1368                                          if (opp_rbf(rbf,miga[0]) == rbf->next)
1369 <                                                return(1);
1369 >                                                return(0);
1370                                  break;
1371                          }
1372                  }
1373 <                return(0);                      /* outside range! */
1373 >                return(-1);                     /* outside range! */
1374          }
1375          {                                       /* else use triangle mesh */
1376 +                const int       sym = use_symmetry(invec);
1377                  unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
1378                  int             pstart[2];
1379                  RBFNODE         *vother;
# Line 1248 | Line 1384 | get_interp(MIGRATION *miga[3], const FVECT invec)
1384                  memset(floodmap, 0, sizeof(floodmap));
1385                                                  /* call flooding function */
1386                  if (!identify_tri(miga, floodmap, pstart[0], pstart[1]))
1387 <                        return(0);              /* outside mesh */
1387 >                        return(-1);             /* outside mesh */
1388                  if ((miga[0] == NULL) | (miga[2] == NULL))
1389 <                        return(0);              /* should never happen */
1389 >                        return(-1);             /* should never happen */
1390                  if (miga[1] == NULL)
1391 <                        return(1);              /* on edge */
1391 >                        return(sym);            /* on edge */
1392                                                  /* verify triangle */
1393                  if (!order_triangle(miga)) {
1394   #ifdef DEBUG
# Line 1272 | Line 1408 | get_interp(MIGRATION *miga[3], const FVECT invec)
1408   #ifdef DEBUG
1409                                  fputs("No triangle in get_interp()\n", stderr);
1410   #endif
1411 <                                return(0);
1411 >                                return(-1);
1412                          }
1413                                                  /* reassign other two edges */
1414                          for (ej = vother->ejl; ej != NULL;
# Line 1287 | Line 1423 | get_interp(MIGRATION *miga[3], const FVECT invec)
1423   #ifdef DEBUG
1424                                  fputs("Bad triangle in get_interp()\n", stderr);
1425   #endif
1426 <                                return(0);
1426 >                                return(-1);
1427                          }
1428                  }
1429 <                return(3);                      /* return in standard order */
1429 >                return(sym);                    /* return in standard order */
1430          }
1431   }
1432  
# Line 1364 | Line 1500 | e_advect_rbf(const MIGRATION *mig, const FVECT invec)
1500          rbf->vtotal *= mig->rbfv[0]->vtotal;    /* turn ratio into actual */
1501          return(rbf);
1502   memerr:
1503 <        fputs("Out of memory in e_advect_rbf()\n", stderr);
1503 >        fprintf(stderr, "%s: Out of memory in e_advect_rbf()\n", progname);
1504          exit(1);
1505          return(NULL);   /* pro forma return */
1506   }
# Line 1373 | Line 1509 | memerr:
1509   static RBFNODE *
1510   advect_rbf(const FVECT invec)
1511   {
1512 +        FVECT           sivec;
1513          MIGRATION       *miga[3];
1514          RBFNODE         *rbf;
1515 +        int             sym;
1516          float           mbfact, mcfact;
1517          int             n, i, j, k;
1518          FVECT           v0, v1, v2;
1519          double          s, t;
1520  
1521 <        if (!get_interp(miga, invec))           /* can't interpolate? */
1521 >        VCOPY(sivec, invec);                    /* find triangle/edge */
1522 >        sym = get_interp(miga, sivec);
1523 >        if (sym < 0)                            /* can't interpolate? */
1524                  return(NULL);
1525 <        if (miga[1] == NULL)                    /* advect along edge? */
1526 <                return(e_advect_rbf(miga[0], invec));
1525 >        if (miga[1] == NULL) {                  /* advect along edge? */
1526 >                rbf = e_advect_rbf(miga[0], sivec);
1527 >                rev_rbf_symmetry(rbf, sym);
1528 >                return(rbf);
1529 >        }
1530   #ifdef DEBUG
1531          if (miga[0]->rbfv[0] != miga[2]->rbfv[0] |
1532                          miga[0]->rbfv[1] != miga[1]->rbfv[0] |
1533                          miga[1]->rbfv[1] != miga[2]->rbfv[1]) {
1534 <                fputs("Triangle vertex screw-up!\n", stderr);
1534 >                fprintf(stderr, "%s: Triangle vertex screw-up!\n", progname);
1535                  exit(1);
1536          }
1537   #endif
# Line 1397 | Line 1540 | advect_rbf(const FVECT invec)
1540          normalize(v0);
1541          fcross(v2, miga[1]->rbfv[0]->invec, miga[1]->rbfv[1]->invec);
1542          normalize(v2);
1543 <        fcross(v1, invec, miga[1]->rbfv[1]->invec);
1543 >        fcross(v1, sivec, miga[1]->rbfv[1]->invec);
1544          normalize(v1);
1545          s = acos(DOT(v0,v1)) / acos(DOT(v0,v2));
1546          geodesic(v1, miga[0]->rbfv[0]->invec, miga[0]->rbfv[1]->invec,
1547                          s, GEOD_REL);
1548 <        t = acos(DOT(v1,invec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
1548 >        t = acos(DOT(v1,sivec)) / acos(DOT(v1,miga[1]->rbfv[1]->invec));
1549          n = 0;                                  /* count migrating particles */
1550          for (i = 0; i < mtx_nrows(miga[0]); i++)
1551              for (j = 0; j < mtx_ncols(miga[0]); j++)
# Line 1417 | Line 1560 | advect_rbf(const FVECT invec)
1560   #endif
1561          rbf = (RBFNODE *)malloc(sizeof(RBFNODE) + sizeof(RBFVAL)*(n-1));
1562          if (rbf == NULL) {
1563 <                fputs("Out of memory in advect_rbf()\n", stderr);
1563 >                fprintf(stderr, "%s: Out of memory in advect_rbf()\n", progname);
1564                  exit(1);
1565          }
1566          rbf->next = NULL; rbf->ejl = NULL;
1567 <        VCOPY(rbf->invec, invec);
1567 >        VCOPY(rbf->invec, sivec);
1568          rbf->nrbf = n;
1569          n = 0;                                  /* compute RBF lobes */
1570          mbfact = s * miga[0]->rbfv[1]->vtotal/miga[0]->rbfv[0]->vtotal *
# Line 1467 | Line 1610 | advect_rbf(const FVECT invec)
1610              }
1611          }
1612          rbf->vtotal = miga[0]->rbfv[0]->vtotal * (mbfact + mcfact);
1613 +        rev_rbf_symmetry(rbf, sym);
1614          return(rbf);
1615   }
1616  

Diff Legend

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