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.13 by greg, Fri Sep 21 05:17:22 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;
163 –
167                                          /* keep in ascending theta order */
168          for (rbf_last = NULL, rbf = dsf_list;
169                          single_plane_incident & (rbf != NULL);
# Line 219 | Line 222 | make_rbfrep(void)
222                  }
223                                  /* iterate to improve interpolation accuracy */
224          do {
225 <                double  dsum = .0, dsum2 = .0;
225 >                double  dsum = 0, dsum2 = 0;
226                  nn = 0;
227                  for (i = 0; i < GRIDRES; i++)
228                      for (j = 0; j < GRIDRES; j++)
# Line 248 | Line 251 | make_rbfrep(void)
251  
252          insert_dsf(newnode);
253                                  /* adjust sampling resolution */
254 <        samp_order = log(2./R2ANG(minrad))/log(2.) + .5;
254 >        samp_order = log(2./R2ANG(minrad))/M_LN2 + .5;
255  
256          return(newnode);
257   }
# Line 647 | Line 650 | make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
650                                  round(theta), round(phi));
651                  theta = acos(to_rbf->invec[2])*(180./M_PI);
652                  phi = atan2(to_rbf->invec[1],to_rbf->invec[0])*(180./M_PI);
653 <                fprintf(stderr, "(%d,%d)\n", round(theta), round(phi));
653 >                fprintf(stderr, "(%d,%d)", round(theta), round(phi));
654          }
655   #endif
656          newmig->next = NULL;
# Line 661 | Line 664 | make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
664          for (i = to_rbf->nrbf; i--; )
665                  dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
666                                                  /* move a bit at a time */
667 <        while (total_rem > end_thresh)
667 >        while (total_rem > end_thresh) {
668                  total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
669 + #ifdef DEBUG
670 +                fputc('.', stderr);
671 + /*XXX*/break;
672 + #endif
673 +        }
674 + #ifdef DEBUG
675 +        fputs("done.\n", stderr);
676 + #endif
677  
678          free(pmtx);                             /* free working arrays */
679          free(src_rem);
# Line 714 | Line 725 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
725          RBFNODE         *tv;
726  
727          rbfv[0] = rbfv[1] = NULL;
728 +        if (mig == NULL)
729 +                return(0);
730          for (ej = mig->rbfv[0]->ejl; ej != NULL;
731                                  ej = nextedge(mig->rbfv[0],ej)) {
732                  if (ej == mig)
# Line 730 | Line 743 | get_triangles(RBFNODE *rbfv[2], const MIGRATION *mig)
743          return((rbfv[0] != NULL) + (rbfv[1] != NULL));
744   }
745  
746 + /* Check if prospective vertex would create overlapping triangle */
747 + static int
748 + overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
749 + {
750 +        const MIGRATION *ej;
751 +        RBFNODE         *vother[2];
752 +        int             im_rev;
753 +                                        /* find shared edge in mesh */
754 +        for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
755 +                const RBFNODE   *tv = opp_rbf(pv,ej);
756 +                if (tv == bv0) {
757 +                        im_rev = is_rev_tri(ej->rbfv[0]->invec,
758 +                                        ej->rbfv[1]->invec, bv1->invec);
759 +                        break;
760 +                }
761 +                if (tv == bv1) {
762 +                        im_rev = is_rev_tri(ej->rbfv[0]->invec,
763 +                                        ej->rbfv[1]->invec, bv0->invec);
764 +                        break;
765 +                }
766 +        }
767 +        if (!get_triangles(vother, ej))
768 +                return(0);
769 +        return(vother[im_rev] != NULL);
770 + }
771 +
772   /* Find context hull vertex to complete triangle (oriented call) */
773   static RBFNODE *
774   find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
775   {
776          FVECT   vmid, vor;
777          RBFNODE *rbf, *rbfbest = NULL;
778 <        double  dprod2, bestdprod2 = 0.5;
778 >        double  dprod2, area2, bestarea2 = FHUGE, bestdprod2 = 0.5;
779  
780          VADD(vmid, rbf0->invec, rbf1->invec);
781          if (normalize(vmid) == 0)
# Line 749 | Line 788 | find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rb
788                  dprod2 = DOT(vor, vmid);
789                  if (dprod2 <= FTINY)
790                          continue;               /* wrong orientation */
791 <                dprod2 *= dprod2 / DOT(vor,vor);
792 <                if (dprod2 > bestdprod2) {      /* more convex than prev? */
791 >                area2 = DOT(vor, vor);
792 >                dprod2 *= dprod2 / area2;
793 >                if (dprod2 > bestdprod2 +
794 >                                FTINY*(1 - 2*(area2 < bestarea2)) &&
795 >                                !overlaps_tri(rbf0, rbf1, rbf)) {
796                          rbfbest = rbf;
797                          bestdprod2 = dprod2;
798 +                        bestarea2 = area2;
799                  }
800          }
801 <        return(rbf);
801 >        return(rbfbest);
802   }
803  
804   /* Create new migration edge and grow mesh recursively around it */
805   static void
806 < mesh_from_edge(RBFNODE *rbf0, RBFNODE *rbf1)
806 > mesh_from_edge(MIGRATION *edge)
807   {
808 <        MIGRATION       *newej;
808 >        MIGRATION       *ej0, *ej1;
809          RBFNODE         *tvert[2];
767 –
768 –        if (rbf0 < rbf1)                        /* avoid migration loops */
769 –                newej = make_migration(rbf0, rbf1);
770 –        else
771 –                newej = make_migration(rbf1, rbf0);
810                                                  /* triangle on either side? */
811 <        get_triangles(tvert, newej);
812 <        if (tvert[0] == NULL) {                 /* recurse on new right edge */
813 <                tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
811 >        get_triangles(tvert, edge);
812 >        if (tvert[0] == NULL) {                 /* grow mesh on right */
813 >                tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
814                  if (tvert[0] != NULL) {
815 <                        mesh_from_edge(rbf0, tvert[0]);
816 <                        mesh_from_edge(rbf1, tvert[0]);
815 >                        if (tvert[0] > edge->rbfv[0])
816 >                                ej0 = make_migration(edge->rbfv[0], tvert[0]);
817 >                        else
818 >                                ej0 = make_migration(tvert[0], edge->rbfv[0]);
819 >                        if (tvert[0] > edge->rbfv[1])
820 >                                ej1 = make_migration(edge->rbfv[1], tvert[0]);
821 >                        else
822 >                                ej1 = make_migration(tvert[0], edge->rbfv[1]);
823 >                        mesh_from_edge(ej0);
824 >                        mesh_from_edge(ej1);
825                  }
826          }
827 <        if (tvert[1] == NULL) {                 /* recurse on new left edge */
828 <                tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
827 >        if (tvert[1] == NULL) {                 /* grow mesh on left */
828 >                tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
829                  if (tvert[1] != NULL) {
830 <                        mesh_from_edge(rbf0, tvert[1]);
831 <                        mesh_from_edge(rbf1, tvert[1]);
830 >                        if (tvert[1] > edge->rbfv[0])
831 >                                ej0 = make_migration(edge->rbfv[0], tvert[1]);
832 >                        else
833 >                                ej0 = make_migration(tvert[1], edge->rbfv[0]);
834 >                        if (tvert[1] > edge->rbfv[1])
835 >                                ej1 = make_migration(edge->rbfv[1], tvert[1]);
836 >                        else
837 >                                ej1 = make_migration(tvert[1], edge->rbfv[1]);
838 >                        mesh_from_edge(ej0);
839 >                        mesh_from_edge(ej1);
840                  }
841          }
842   }
843  
844 + #ifdef DEBUG
845 + #include "random.h"
846 + #include "bmpfile.h"
847 + /* Hash pointer to byte value */
848 + static int
849 + byte_hash(const void *p)
850 + {
851 +        size_t  h = (size_t)p;
852 +        h ^= (size_t)p >> 8;
853 +        h ^= (size_t)p >> 16;
854 +        h ^= (size_t)p >> 24;
855 +        return(h & 0xff);
856 + }
857 + /* Write out BMP image showing edges */
858 + static void
859 + write_edge_image(const char *fname)
860 + {
861 +        BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
862 +        BMPWriter       *wtr;
863 +        int             i, j;
864 +
865 +        fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
866 +        hdr->compr = BI_RLE8;
867 +        for (i = 256; --i; ) {                  /* assign random color map */
868 +                hdr->palette[i].r = random() & 0xff;
869 +                hdr->palette[i].r = random() & 0xff;
870 +                hdr->palette[i].r = random() & 0xff;
871 +        }
872 +        hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
873 +                                                /* open output */
874 +        wtr = BMPopenOutputFile(fname, hdr);
875 +        if (wtr == NULL) {
876 +                free(hdr);
877 +                return;
878 +        }
879 +        for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
880 +                for (j = 0; j < GRIDRES; j++)
881 +                        wtr->scanline[j] = byte_hash(mig_grid[i][j]);
882 +                if (BMPwriteScanline(wtr) != BIR_OK)
883 +                        break;
884 +        }
885 +        BMPcloseOutput(wtr);                    /* close & clean up */
886 + }
887 + #endif
888 +
889   /* Draw edge list into mig_grid array */
890   static void
891   draw_edges()
# Line 830 | Line 929 | draw_edges()
929          if (nnull)
930                  fprintf(stderr, "Warning: %d of %d edges are null\n",
931                                  nnull, ntot);
932 + #ifdef DEBUG
933 +        write_edge_image("bsdf_edges.bmp");
934 + #endif
935   }
936          
937   /* Build our triangle mesh from recorded RBFs */
# Line 837 | Line 939 | static void
939   build_mesh()
940   {
941          double          best2 = M_PI*M_PI;
942 <        RBFNODE         *rbf, *rbf_near = NULL;
942 >        RBFNODE         *shrt_edj[2];
943 >        RBFNODE         *rbf0, *rbf1;
944                                                  /* check if isotropic */
945          if (single_plane_incident) {
946 <                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
947 <                        if (rbf->next != NULL)
948 <                                make_migration(rbf, rbf->next);
946 >                for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
947 >                        if (rbf0->next != NULL)
948 >                                make_migration(rbf0, rbf0->next);
949                  return;
950          }
951 <                                                /* find RBF nearest to head */
952 <        if (dsf_list == NULL)
953 <                return;
954 <        for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
955 <                double  dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
951 >                                                /* start w/ shortest edge */
952 >        shrt_edj[0] = shrt_edj[1] = NULL;
953 >        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
954 >            for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
955 >                double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
956                  if (dist2 < best2) {
957 <                        rbf_near = rbf;
957 >                        shrt_edj[0] = rbf0;
958 >                        shrt_edj[1] = rbf1;
959                          best2 = dist2;
960                  }
961          }
962 <        if (rbf_near == NULL) {
963 <                fputs("Cannot find nearest point for first edge\n", stderr);
962 >        if (shrt_edj[0] == NULL) {
963 >                fputs("Cannot find shortest edge\n", stderr);
964                  exit(1);
965          }
966                                                  /* build mesh from this edge */
967 <        mesh_from_edge(dsf_list, rbf_near);
967 >        if (shrt_edj[0] < shrt_edj[1])
968 >                mesh_from_edge(make_migration(shrt_edj[0], shrt_edj[1]));
969 >        else
970 >                mesh_from_edge(make_migration(shrt_edj[1], shrt_edj[0]));
971                                                  /* draw edge list into grid */
972          draw_edges();
973   }
# Line 1086 | Line 1193 | advect_rbf(const FVECT invec)
1193  
1194          if (!get_interp(miga, invec))           /* can't interpolate? */
1195                  return(NULL);
1196 <        if (miga[1] == NULL)                    /* along edge? */
1196 >        if (miga[1] == NULL)                    /* advect along edge? */
1197                  return(e_advect_rbf(miga[0], invec));
1198                                                  /* put in standard order */
1199          order_triangle(miga);
# Line 1186 | Line 1293 | interp_isotropic()
1293          int             ix, ox, oy;
1294          FVECT           ivec, ovec;
1295          double          bsdf;
1296 <
1296 > #if DEBUG
1297 >        fprintf(stderr, "Writing isotropic order %d ", samp_order);
1298 >        if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1299 >        else fputs("raw data\n", stderr);
1300 > #endif
1301          if (pctcull >= 0) {                     /* begin output */
1302                  sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1303                                  pctcull, samp_order);
# Line 1240 | Line 1351 | interp_anisotropic()
1351          int             ix, iy, ox, oy;
1352          FVECT           ivec, ovec;
1353          double          bsdf;
1354 <
1354 > #if DEBUG
1355 >        fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1356 >        if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1357 >        else fputs("raw data\n", stderr);
1358 > #endif
1359          if (pctcull >= 0) {                     /* begin output */
1360                  sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1361                                  pctcull, samp_order);

Diff Legend

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