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.12 by greg, Thu Sep 20 01:23:36 2012 UTC vs.
Revision 2.14 by greg, Sat Sep 22 23:10:24 2012 UTC

# Line 89 | Line 89 | static MIGRATION       *mig_grid[GRIDRES][GRIDRES];
89   #define round(v)        (int)((v) + .5 - ((v) < -.5))
90  
91   char                    *progname;
92 <                                /* percentage to cull (<0 to turn off) */
92 >
93 > #ifdef DEBUG                    /* percentage to cull (<0 to turn off) */
94 > int                     pctcull = -1;
95 > #else
96   int                     pctcull = 90;
97 <                                /* sampling order */
97 > #endif
98 >                                /* sampling order (set by data density) */
99   int                     samp_order = 0;
100  
101   /* Compute volume associated with Gaussian lobe */
# Line 160 | Line 164 | static void
164   insert_dsf(RBFNODE *newrbf)
165   {
166          RBFNODE         *rbf, *rbf_last;
167 <
167 >                                        /* check for redundant meas. */
168 >        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
169 >                if (DOT(rbf->invec, newrbf->invec) >= 1.-FTINY) {
170 >                        fputs("Duplicate incident measurement (ignored)\n", stderr);
171 >                        free(newrbf);
172 >                        return;
173 >                }
174                                          /* keep in ascending theta order */
175          for (rbf_last = NULL, rbf = dsf_list;
176                          single_plane_incident & (rbf != NULL);
# Line 219 | Line 229 | make_rbfrep(void)
229                  }
230                                  /* iterate to improve interpolation accuracy */
231          do {
232 <                double  dsum = .0, dsum2 = .0;
232 >                double  dsum = 0, dsum2 = 0;
233                  nn = 0;
234                  for (i = 0; i < GRIDRES; i++)
235                      for (j = 0; j < GRIDRES; j++)
# Line 248 | Line 258 | make_rbfrep(void)
258  
259          insert_dsf(newnode);
260                                  /* adjust sampling resolution */
261 <        samp_order = log(2./R2ANG(minrad))/log(2.) + .5;
261 >        samp_order = log(2./R2ANG(minrad))/M_LN2 + .5;
262  
263          return(newnode);
264   }
# Line 558 | Line 568 | migration_step(MIGRATION *mig, double *src_rem, double
568   {
569          static double   *src_cost = NULL;
570          int             n_alloc = 0;
571 <        const double    maxamt = 2./(mtx_nrows(mig)*mtx_ncols(mig));
571 >        const double    maxamt = .1; /* 2./(mtx_nrows(mig)*mtx_ncols(mig)); */
572          double          amt = 0;
573          struct {
574                  int     s, d;   /* source and destination */
# Line 620 | Line 630 | migration_step(MIGRATION *mig, double *src_rem, double
630          return(best.amt);
631   }
632  
633 + #ifdef DEBUG
634 + static char *
635 + thetaphi(const FVECT v)
636 + {
637 +        static char     buf[128];
638 +        double          theta, phi;
639 +
640 +        theta = 180./M_PI*acos(v[2]);
641 +        phi = 180./M_PI*atan2(v[1],v[0]);
642 +        sprintf(buf, "(%.0f,%.0f)", theta, phi);
643 +
644 +        return(buf);
645 + }
646 + #endif
647 +
648   /* Compute (and insert) migration along directed edge */
649   static MIGRATION *
650 < make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
650 > make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf, int creat_only)
651   {
652          const double    end_thresh = 0.02/(from_rbf->nrbf*to_rbf->nrbf);
653 <        float           *pmtx = price_routes(from_rbf, to_rbf);
654 <        MIGRATION       *newmig = (MIGRATION *)malloc(sizeof(MIGRATION) +
655 <                                                        sizeof(float) *
631 <                                        (from_rbf->nrbf*to_rbf->nrbf - 1));
632 <        double          *src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
633 <        double          *dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
653 >        float           *pmtx;
654 >        MIGRATION       *newmig;
655 >        double          *src_rem, *dst_rem;
656          double          total_rem = 1.;
657          int             i;
658 <
658 >                                                /* check if exists already */
659 >        for (newmig = from_rbf->ejl; newmig != NULL;
660 >                        newmig = nextedge(from_rbf,newmig))
661 >                if (newmig->rbfv[1] == to_rbf)
662 >                        return(creat_only ? (MIGRATION *)NULL : newmig);
663 >                                                /* else allocate */
664 >        pmtx = price_routes(from_rbf, to_rbf);
665 >        newmig = (MIGRATION *)malloc(sizeof(MIGRATION) + sizeof(float) *
666 >                                        (from_rbf->nrbf*to_rbf->nrbf - 1));
667 >        src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
668 >        dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
669          if ((newmig == NULL) | (src_rem == NULL) | (dst_rem == NULL)) {
670                  fputs("Out of memory in make_migration()\n", stderr);
671                  exit(1);
672          }
673   #ifdef DEBUG
674          {
675 <                double  theta, phi;
676 <                theta = acos(from_rbf->invec[2])*(180./M_PI);
677 <                phi = atan2(from_rbf->invec[1],from_rbf->invec[0])*(180./M_PI);
646 <                fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
647 <                                round(theta), round(phi));
648 <                theta = acos(to_rbf->invec[2])*(180./M_PI);
649 <                phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
650 <                fprintf(stderr, "(%d,%d)\n", round(theta), round(phi));
675 >                fprintf(stderr, "Building path from (theta,phi) %s ",
676 >                                thetaphi(from_rbf->invec));
677 >                fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
678          }
679   #endif
680          newmig->next = NULL;
# Line 661 | Line 688 | make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
688          for (i = to_rbf->nrbf; i--; )
689                  dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
690                                                  /* move a bit at a time */
691 <        while (total_rem > end_thresh)
691 >        while (total_rem > end_thresh) {
692                  total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
693 + #ifdef DEBUG
694 +                /* fputc('.', stderr); */
695 +                fprintf(stderr, "\n%.9f remaining...", total_rem);
696 + #endif
697 +        }
698 + #ifdef DEBUG
699 +        fputs("done.\n", stderr);
700 + #endif
701  
702          free(pmtx);                             /* free working arrays */
703          free(src_rem);
# Line 714 | Line 749 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
749          RBFNODE         *tv;
750  
751          rbfv[0] = rbfv[1] = NULL;
752 +        if (mig == NULL)
753 +                return(0);
754          for (ej = mig->rbfv[0]->ejl; ej != NULL;
755                                  ej = nextedge(mig->rbfv[0],ej)) {
756                  if (ej == mig)
# Line 730 | Line 767 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
767          return((rbfv[0] != NULL) + (rbfv[1] != NULL));
768   }
769  
770 + /* Check if prospective vertex would create overlapping triangle */
771 + static int
772 + overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
773 + {
774 +        const MIGRATION *ej;
775 +        RBFNODE         *vother[2];
776 +        int             im_rev;
777 +                                        /* find shared edge in mesh */
778 +        for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
779 +                const RBFNODE   *tv = opp_rbf(pv,ej);
780 +                if (tv == bv0) {
781 +                        im_rev = is_rev_tri(ej->rbfv[0]->invec,
782 +                                        ej->rbfv[1]->invec, bv1->invec);
783 +                        break;
784 +                }
785 +                if (tv == bv1) {
786 +                        im_rev = is_rev_tri(ej->rbfv[0]->invec,
787 +                                        ej->rbfv[1]->invec, bv0->invec);
788 +                        break;
789 +                }
790 +        }
791 +        if (!get_triangles(vother, ej))
792 +                return(0);
793 +        return(vother[im_rev] != NULL);
794 + }
795 +
796   /* Find context hull vertex to complete triangle (oriented call) */
797   static RBFNODE *
798   find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
799   {
800 <        FVECT   vmid, vor;
800 >        FVECT   vmid, vejn, vp;
801          RBFNODE *rbf, *rbfbest = NULL;
802 <        double  dprod2, bestdprod2 = 0.5;
802 >        double  dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
803  
804 +        VSUB(vejn, rbf1->invec, rbf0->invec);
805          VADD(vmid, rbf0->invec, rbf1->invec);
806 <        if (normalize(vmid) == 0)
806 >        if (normalize(vejn) == 0 || normalize(vmid) == 0)
807                  return(NULL);
808                                                  /* XXX exhaustive search */
809          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
810                  if ((rbf == rbf0) | (rbf == rbf1))
811                          continue;
812 <                tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
813 <                dprod2 = DOT(vor, vmid);
750 <                if (dprod2 <= FTINY)
812 >                tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
813 >                if (DOT(vp, vmid) <= FTINY)
814                          continue;               /* wrong orientation */
815 <                dprod2 *= dprod2 / DOT(vor,vor);
816 <                if (dprod2 > bestdprod2) {      /* more convex than prev? */
817 <                        rbfbest = rbf;
818 <                        bestdprod2 = dprod2;
819 <                }
815 >                area2 = DOT(vp,vp);
816 >                VSUB(vp, rbf->invec, rbf0->invec);
817 >                dprod = -DOT(vp, vejn);
818 >                VSUM(vp, vp, vejn, dprod);
819 >                dprod = DOT(vp, vmid) / VLEN(vp);
820 >                if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
821 >                        continue;               /* found better already */
822 >                if (overlaps_tri(rbf0, rbf1, rbf))
823 >                        continue;               /* overlaps another triangle */
824 >                rbfbest = rbf;
825 >                bestdprod = dprod;              /* new one to beat */
826 >                bestarea2 = area2;
827          }
828 <        return(rbf);
828 >        return(rbfbest);
829   }
830  
831   /* Create new migration edge and grow mesh recursively around it */
832   static void
833 < mesh_from_edge(RBFNODE *rbf0, RBFNODE *rbf1)
833 > mesh_from_edge(MIGRATION *edge)
834   {
835 <        MIGRATION       *newej;
835 >        MIGRATION       *ej0, *ej1;
836          RBFNODE         *tvert[2];
837 <
838 <        if (rbf0 < rbf1)                        /* avoid migration loops */
839 <                newej = make_migration(rbf0, rbf1);
770 <        else
771 <                newej = make_migration(rbf1, rbf0);
837 >        
838 >        if (edge == NULL)
839 >                return;
840                                                  /* triangle on either side? */
841 <        get_triangles(tvert, newej);
842 <        if (tvert[0] == NULL) {                 /* recurse on new right edge */
843 <                tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
841 >        get_triangles(tvert, edge);
842 >        if (tvert[0] == NULL) {                 /* grow mesh on right */
843 >                tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
844                  if (tvert[0] != NULL) {
845 <                        mesh_from_edge(rbf0, tvert[0]);
846 <                        mesh_from_edge(rbf1, tvert[0]);
845 >                        if (tvert[0] > edge->rbfv[0])
846 >                                ej0 = make_migration(edge->rbfv[0], tvert[0], 1);
847 >                        else
848 >                                ej0 = make_migration(tvert[0], edge->rbfv[0], 1);
849 >                        if (tvert[0] > edge->rbfv[1])
850 >                                ej1 = make_migration(edge->rbfv[1], tvert[0], 1);
851 >                        else
852 >                                ej1 = make_migration(tvert[0], edge->rbfv[1], 1);
853 >                        mesh_from_edge(ej0);
854 >                        mesh_from_edge(ej1);
855                  }
856 <        }
857 <        if (tvert[1] == NULL) {                 /* recurse on new left edge */
782 <                tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
856 >        } else if (tvert[1] == NULL) {          /* grow mesh on left */
857 >                tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
858                  if (tvert[1] != NULL) {
859 <                        mesh_from_edge(rbf0, tvert[1]);
860 <                        mesh_from_edge(rbf1, tvert[1]);
859 >                        if (tvert[1] > edge->rbfv[0])
860 >                                ej0 = make_migration(edge->rbfv[0], tvert[1], 1);
861 >                        else
862 >                                ej0 = make_migration(tvert[1], edge->rbfv[0], 1);
863 >                        if (tvert[1] > edge->rbfv[1])
864 >                                ej1 = make_migration(edge->rbfv[1], tvert[1], 1);
865 >                        else
866 >                                ej1 = make_migration(tvert[1], edge->rbfv[1], 1);
867 >                        mesh_from_edge(ej0);
868 >                        mesh_from_edge(ej1);
869                  }
870          }
871   }
872  
873 + #ifdef DEBUG
874 + #include "random.h"
875 + #include "bmpfile.h"
876 + /* Hash pointer to byte value */
877 + static int
878 + byte_hash(const void *p)
879 + {
880 +        size_t  h = (size_t)p;
881 +        h ^= (size_t)p >> 8;
882 +        h ^= (size_t)p >> 16;
883 +        h ^= (size_t)p >> 24;
884 +        return(h & 0xff);
885 + }
886 + /* Write out BMP image showing edges */
887 + static void
888 + write_edge_image(const char *fname)
889 + {
890 +        BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
891 +        BMPWriter       *wtr;
892 +        int             i, j;
893 +
894 +        fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
895 +        hdr->compr = BI_RLE8;
896 +        for (i = 256; --i; ) {                  /* assign random color map */
897 +                hdr->palette[i].r = random() & 0xff;
898 +                hdr->palette[i].r = random() & 0xff;
899 +                hdr->palette[i].r = random() & 0xff;
900 +        }
901 +        hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
902 +                                                /* open output */
903 +        wtr = BMPopenOutputFile(fname, hdr);
904 +        if (wtr == NULL) {
905 +                free(hdr);
906 +                return;
907 +        }
908 +        for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
909 +                for (j = 0; j < GRIDRES; j++)
910 +                        wtr->scanline[j] = byte_hash(mig_grid[i][j]);
911 +                if (BMPwriteScanline(wtr) != BIR_OK)
912 +                        break;
913 +        }
914 +        BMPcloseOutput(wtr);                    /* close & clean up */
915 + }
916 + #endif
917 +
918   /* Draw edge list into mig_grid array */
919   static void
920   draw_edges()
# Line 830 | Line 958 | draw_edges()
958          if (nnull)
959                  fprintf(stderr, "Warning: %d of %d edges are null\n",
960                                  nnull, ntot);
961 + #ifdef DEBUG
962 +        write_edge_image("bsdf_edges.bmp");
963 + #endif
964   }
965          
966   /* Build our triangle mesh from recorded RBFs */
# Line 837 | Line 968 | static void
968   build_mesh()
969   {
970          double          best2 = M_PI*M_PI;
971 <        RBFNODE         *rbf, *rbf_near = NULL;
971 >        RBFNODE         *shrt_edj[2];
972 >        RBFNODE         *rbf0, *rbf1;
973                                                  /* check if isotropic */
974          if (single_plane_incident) {
975 <                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
976 <                        if (rbf->next != NULL)
977 <                                make_migration(rbf, rbf->next);
975 >                for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
976 >                        if (rbf0->next != NULL)
977 >                                make_migration(rbf0, rbf0->next, 1);
978                  return;
979          }
980 <                                                /* find RBF nearest to head */
981 <        if (dsf_list == NULL)
982 <                return;
983 <        for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
984 <                double  dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
980 >                                                /* start w/ shortest edge */
981 >        shrt_edj[0] = shrt_edj[1] = NULL;
982 >        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
983 >            for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
984 >                double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
985                  if (dist2 < best2) {
986 <                        rbf_near = rbf;
986 >                        shrt_edj[0] = rbf0;
987 >                        shrt_edj[1] = rbf1;
988                          best2 = dist2;
989                  }
990          }
991 <        if (rbf_near == NULL) {
992 <                fputs("Cannot find nearest point for first edge\n", stderr);
991 >        if (shrt_edj[0] == NULL) {
992 >                fputs("Cannot find shortest edge\n", stderr);
993                  exit(1);
994          }
995                                                  /* build mesh from this edge */
996 <        mesh_from_edge(dsf_list, rbf_near);
996 >        if (shrt_edj[0] < shrt_edj[1])
997 >                mesh_from_edge(make_migration(shrt_edj[0], shrt_edj[1], 0));
998 >        else
999 >                mesh_from_edge(make_migration(shrt_edj[1], shrt_edj[0], 0));
1000                                                  /* draw edge list into grid */
1001          draw_edges();
1002   }
# Line 1086 | Line 1222 | advect_rbf(const FVECT invec)
1222  
1223          if (!get_interp(miga, invec))           /* can't interpolate? */
1224                  return(NULL);
1225 <        if (miga[1] == NULL)                    /* along edge? */
1225 >        if (miga[1] == NULL)                    /* advect along edge? */
1226                  return(e_advect_rbf(miga[0], invec));
1227                                                  /* put in standard order */
1228          order_triangle(miga);
# Line 1186 | Line 1322 | interp_isotropic()
1322          int             ix, ox, oy;
1323          FVECT           ivec, ovec;
1324          double          bsdf;
1325 <
1325 > #if DEBUG
1326 >        fprintf(stderr, "Writing isotropic order %d ", samp_order);
1327 >        if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1328 >        else fputs("raw data\n", stderr);
1329 > #endif
1330          if (pctcull >= 0) {                     /* begin output */
1331                  sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1332                                  pctcull, samp_order);
# Line 1240 | Line 1380 | interp_anisotropic()
1380          int             ix, iy, ox, oy;
1381          FVECT           ivec, ovec;
1382          double          bsdf;
1383 <
1383 > #if DEBUG
1384 >        fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1385 >        if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1386 >        else fputs("raw data\n", stderr);
1387 > #endif
1388          if (pctcull >= 0) {                     /* begin output */
1389                  sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1390                                  pctcull, samp_order);

Diff Legend

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