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.15 by greg, Sun Sep 23 16:45:20 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))         /* triangle on same side? */
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 +        /* Find triangle with minimum rotation from perpendicular */
810          for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
811                  if ((rbf == rbf0) | (rbf == rbf1))
812                          continue;
813 <                tri_orient(vor, rbf0->invec, rbf1->invec, rbf->invec);
814 <                dprod2 = DOT(vor, vmid);
750 <                if (dprod2 <= FTINY)
813 >                tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
814 >                if (DOT(vp, vmid) <= FTINY)
815                          continue;               /* wrong orientation */
816 <                dprod2 *= dprod2 / DOT(vor,vor);
817 <                if (dprod2 > bestdprod2) {      /* more convex than prev? */
818 <                        rbfbest = rbf;
819 <                        bestdprod2 = dprod2;
820 <                }
816 >                area2 = .25*DOT(vp,vp);
817 >                VSUB(vp, rbf->invec, rbf0->invec);
818 >                dprod = -DOT(vp, vejn);
819 >                VSUM(vp, vp, vejn, dprod);      /* above guarantees non-zero */
820 >                dprod = DOT(vp, vmid) / VLEN(vp);
821 >                if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
822 >                        continue;               /* found better already */
823 >                if (overlaps_tri(rbf0, rbf1, rbf))
824 >                        continue;               /* overlaps another triangle */
825 >                rbfbest = rbf;
826 >                bestdprod = dprod;              /* new one to beat */
827 >                bestarea2 = area2;
828          }
829 <        return(rbf);
829 >        return(rbfbest);
830   }
831  
832   /* Create new migration edge and grow mesh recursively around it */
833   static void
834 < mesh_from_edge(RBFNODE *rbf0, RBFNODE *rbf1)
834 > mesh_from_edge(MIGRATION *edge)
835   {
836 <        MIGRATION       *newej;
836 >        MIGRATION       *ej0, *ej1;
837          RBFNODE         *tvert[2];
838 <
839 <        if (rbf0 < rbf1)                        /* avoid migration loops */
840 <                newej = make_migration(rbf0, rbf1);
770 <        else
771 <                newej = make_migration(rbf1, rbf0);
838 >        
839 >        if (edge == NULL)
840 >                return;
841                                                  /* triangle on either side? */
842 <        get_triangles(tvert, newej);
843 <        if (tvert[0] == NULL) {                 /* recurse on new right edge */
844 <                tvert[0] = find_chull_vert(newej->rbfv[0], newej->rbfv[1]);
842 >        get_triangles(tvert, edge);
843 >        if (tvert[0] == NULL) {                 /* grow mesh on right */
844 >                tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
845                  if (tvert[0] != NULL) {
846 <                        mesh_from_edge(rbf0, tvert[0]);
847 <                        mesh_from_edge(rbf1, tvert[0]);
846 >                        if (tvert[0] > edge->rbfv[0])
847 >                                ej0 = make_migration(edge->rbfv[0], tvert[0], 1);
848 >                        else
849 >                                ej0 = make_migration(tvert[0], edge->rbfv[0], 1);
850 >                        if (tvert[0] > edge->rbfv[1])
851 >                                ej1 = make_migration(edge->rbfv[1], tvert[0], 1);
852 >                        else
853 >                                ej1 = make_migration(tvert[0], edge->rbfv[1], 1);
854 >                        mesh_from_edge(ej0);
855 >                        mesh_from_edge(ej1);
856                  }
857 <        }
858 <        if (tvert[1] == NULL) {                 /* recurse on new left edge */
782 <                tvert[1] = find_chull_vert(newej->rbfv[1], newej->rbfv[0]);
857 >        } else if (tvert[1] == NULL) {          /* grow mesh on left */
858 >                tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
859                  if (tvert[1] != NULL) {
860 <                        mesh_from_edge(rbf0, tvert[1]);
861 <                        mesh_from_edge(rbf1, tvert[1]);
860 >                        if (tvert[1] > edge->rbfv[0])
861 >                                ej0 = make_migration(edge->rbfv[0], tvert[1], 1);
862 >                        else
863 >                                ej0 = make_migration(tvert[1], edge->rbfv[0], 1);
864 >                        if (tvert[1] > edge->rbfv[1])
865 >                                ej1 = make_migration(edge->rbfv[1], tvert[1], 1);
866 >                        else
867 >                                ej1 = make_migration(tvert[1], edge->rbfv[1], 1);
868 >                        mesh_from_edge(ej0);
869 >                        mesh_from_edge(ej1);
870                  }
871          }
872   }
873  
874 + #ifdef DEBUG
875 + #include "random.h"
876 + #include "bmpfile.h"
877 + /* Hash pointer to byte value (must return 0 for NULL) */
878 + static int
879 + byte_hash(const void *p)
880 + {
881 +        size_t  h = (size_t)p;
882 +        h ^= (size_t)p >> 8;
883 +        h ^= (size_t)p >> 16;
884 +        h ^= (size_t)p >> 24;
885 +        return(h & 0xff);
886 + }
887 + /* Write out BMP image showing edges */
888 + static void
889 + write_edge_image(const char *fname)
890 + {
891 +        BMPHeader       *hdr = BMPmappedHeader(GRIDRES, GRIDRES, 0, 256);
892 +        BMPWriter       *wtr;
893 +        int             i, j;
894 +
895 +        fprintf(stderr, "Writing incident mesh drawing to '%s'\n", fname);
896 +        hdr->compr = BI_RLE8;
897 +        for (i = 256; --i; ) {                  /* assign random color map */
898 +                hdr->palette[i].r = random() & 0xff;
899 +                hdr->palette[i].g = random() & 0xff;
900 +                hdr->palette[i].b = random() & 0xff;
901 +                                                /* reject dark colors */
902 +                i += (hdr->palette[i].r + hdr->palette[i].g +
903 +                                                hdr->palette[i].b < 128);
904 +        }
905 +        hdr->palette[0].r = hdr->palette[0].g = hdr->palette[0].b = 0;
906 +                                                /* open output */
907 +        wtr = BMPopenOutputFile(fname, hdr);
908 +        if (wtr == NULL) {
909 +                free(hdr);
910 +                return;
911 +        }
912 +        for (i = 0; i < GRIDRES; i++) {         /* write scanlines */
913 +                for (j = 0; j < GRIDRES; j++)
914 +                        wtr->scanline[j] = byte_hash(mig_grid[i][j]);
915 +                if (BMPwriteScanline(wtr) != BIR_OK)
916 +                        break;
917 +        }
918 +        BMPcloseOutput(wtr);                    /* close & clean up */
919 + }
920 + #endif
921 +
922   /* Draw edge list into mig_grid array */
923   static void
924   draw_edges()
# Line 830 | Line 962 | draw_edges()
962          if (nnull)
963                  fprintf(stderr, "Warning: %d of %d edges are null\n",
964                                  nnull, ntot);
965 + #ifdef DEBUG
966 +        write_edge_image("bsdf_edges.bmp");
967 + #endif
968   }
969          
970   /* Build our triangle mesh from recorded RBFs */
# Line 837 | Line 972 | static void
972   build_mesh()
973   {
974          double          best2 = M_PI*M_PI;
975 <        RBFNODE         *rbf, *rbf_near = NULL;
975 >        RBFNODE         *shrt_edj[2];
976 >        RBFNODE         *rbf0, *rbf1;
977                                                  /* check if isotropic */
978          if (single_plane_incident) {
979 <                for (rbf = dsf_list; rbf != NULL; rbf = rbf->next)
980 <                        if (rbf->next != NULL)
981 <                                make_migration(rbf, rbf->next);
979 >                for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
980 >                        if (rbf0->next != NULL)
981 >                                make_migration(rbf0, rbf0->next, 1);
982                  return;
983          }
984 <                                                /* find RBF nearest to head */
985 <        if (dsf_list == NULL)
986 <                return;
987 <        for (rbf = dsf_list->next; rbf != NULL; rbf = rbf->next) {
988 <                double  dist2 = 2. - 2.*DOT(dsf_list->invec,rbf->invec);
984 >                                                /* start w/ shortest edge */
985 >        shrt_edj[0] = shrt_edj[1] = NULL;
986 >        for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
987 >            for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
988 >                double  dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
989                  if (dist2 < best2) {
990 <                        rbf_near = rbf;
990 >                        shrt_edj[0] = rbf0;
991 >                        shrt_edj[1] = rbf1;
992                          best2 = dist2;
993                  }
994          }
995 <        if (rbf_near == NULL) {
996 <                fputs("Cannot find nearest point for first edge\n", stderr);
995 >        if (shrt_edj[0] == NULL) {
996 >                fputs("Cannot find shortest edge\n", stderr);
997                  exit(1);
998          }
999                                                  /* build mesh from this edge */
1000 <        mesh_from_edge(dsf_list, rbf_near);
1000 >        if (shrt_edj[0] < shrt_edj[1])
1001 >                mesh_from_edge(make_migration(shrt_edj[0], shrt_edj[1], 0));
1002 >        else
1003 >                mesh_from_edge(make_migration(shrt_edj[1], shrt_edj[0], 0));
1004                                                  /* draw edge list into grid */
1005          draw_edges();
1006   }
# Line 901 | Line 1041 | identify_tri(MIGRATION *miga[3], unsigned char vmap[GR
1041          return(1);                              /* this neighborhood done */
1042   }
1043  
1044 + /* Insert vertex in ordered list */
1045 + static void
1046 + insert_vert(RBFNODE **vlist, RBFNODE *v)
1047 + {
1048 +        int     i, j;
1049 +
1050 +        for (i = 0; vlist[i] != NULL; i++) {
1051 +                if (v == vlist[i])
1052 +                        return;
1053 +                if (v < vlist[i])
1054 +                        break;
1055 +        }
1056 +        for (j = i; vlist[j] != NULL; j++)
1057 +                ;
1058 +        while (j > i) {
1059 +                vlist[j] = vlist[j-1];
1060 +                --j;
1061 +        }
1062 +        vlist[i] = v;
1063 + }
1064 +
1065 + /* Sort triangle edges in standard order */
1066 + static int
1067 + order_triangle(MIGRATION *miga[3])
1068 + {
1069 +        RBFNODE         *vert[7];
1070 +        MIGRATION       *ord[3];
1071 +        int             i;
1072 +                                                /* order vertices, first */
1073 +        memset(vert, 0, sizeof(vert));
1074 +        for (i = 3; i--; ) {
1075 +                if (miga[i] == NULL)
1076 +                        return(0);
1077 +                insert_vert(vert, miga[i]->rbfv[0]);
1078 +                insert_vert(vert, miga[i]->rbfv[1]);
1079 +        }
1080 +                                                /* should be just 3 vertices */
1081 +        if ((vert[3] == NULL) | (vert[4] != NULL))
1082 +                return(0);
1083 +                                                /* identify edge 0 */
1084 +        for (i = 3; i--; )
1085 +                if (miga[i]->rbfv[0] == vert[0] &&
1086 +                                miga[i]->rbfv[1] == vert[1]) {
1087 +                        ord[0] = miga[i];
1088 +                        break;
1089 +                }
1090 +        if (i < 0)
1091 +                return(0);
1092 +                                                /* identify edge 1 */
1093 +        for (i = 3; i--; )
1094 +                if (miga[i]->rbfv[0] == vert[1] &&
1095 +                                miga[i]->rbfv[1] == vert[2]) {
1096 +                        ord[1] = miga[i];
1097 +                        break;
1098 +                }
1099 +        if (i < 0)
1100 +                return(0);
1101 +                                                /* identify edge 2 */
1102 +        for (i = 3; i--; )
1103 +                if (miga[i]->rbfv[0] == vert[0] &&
1104 +                                miga[i]->rbfv[1] == vert[2]) {
1105 +                        ord[2] = miga[i];
1106 +                        break;
1107 +                }
1108 +        if (i < 0)
1109 +                return(0);
1110 +                                                /* reassign order */
1111 +        miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1112 +        return(1);
1113 + }
1114 +
1115   /* Find edge(s) for interpolating the given incident vector */
1116   static int
1117   get_interp(MIGRATION *miga[3], const FVECT invec)
# Line 926 | Line 1137 | get_interp(MIGRATION *miga[3], const FVECT invec)
1137          {                                       /* else use triangle mesh */
1138                  unsigned char   floodmap[GRIDRES][(GRIDRES+7)/8];
1139                  int             pstart[2];
1140 +                RBFNODE         *vother;
1141 +                MIGRATION       *ej;
1142 +                int             i;
1143  
1144                  pos_from_vec(pstart, invec);
1145                  memset(floodmap, 0, sizeof(floodmap));
# Line 936 | Line 1150 | get_interp(MIGRATION *miga[3], const FVECT invec)
1150                          return(0);              /* should never happen */
1151                  if (miga[1] == NULL)
1152                          return(1);              /* on edge */
1153 <                return(3);                      /* else in triangle */
1153 >                                                /* verify triangle */
1154 >                if (!order_triangle(miga)) {
1155 > #ifdef DEBUG
1156 >                        fputs("Munged triangle in get_interp()\n", stderr);
1157 > #endif
1158 >                        vother = NULL;          /* find triangle from edge */
1159 >                        for (i = 3; i--; ) {
1160 >                            RBFNODE     *tpair[2];
1161 >                            if (get_triangles(tpair, miga[i]) &&
1162 >                                        (vother = tpair[ is_rev_tri(
1163 >                                                        miga[i]->rbfv[0]->invec,
1164 >                                                        miga[i]->rbfv[1]->invec,
1165 >                                                        invec) ]) != NULL)
1166 >                                        break;
1167 >                        }
1168 >                        if (vother == NULL) {   /* couldn't find 3rd vertex */
1169 > #ifdef DEBUG
1170 >                                fputs("No triangle in get_interp()\n", stderr);
1171 > #endif
1172 >                                return(0);
1173 >                        }
1174 >                                                /* reassign other two edges */
1175 >                        for (ej = vother->ejl; ej != NULL;
1176 >                                                ej = nextedge(vother,ej)) {
1177 >                                RBFNODE *vorig = opp_rbf(vother,ej);
1178 >                                if (vorig == miga[i]->rbfv[0])
1179 >                                        miga[(i+1)%3] = ej;
1180 >                                else if (vorig == miga[i]->rbfv[1])
1181 >                                        miga[(i+2)%3] = ej;
1182 >                        }
1183 >                        if (!order_triangle(miga)) {
1184 > #ifdef DEBUG
1185 >                                fputs("Bad triangle in get_interp()\n", stderr);
1186 > #endif
1187 >                                return(0);
1188 >                        }
1189 >                }
1190 >                return(3);                      /* return in standard order */
1191          }
1192   }
1193  
# Line 1015 | Line 1266 | memerr:
1266          return(NULL);   /* pro forma return */
1267   }
1268  
1018 – /* Insert vertex in ordered list */
1019 – static void
1020 – insert_vert(RBFNODE **vlist, RBFNODE *v)
1021 – {
1022 –        int     i, j;
1023 –
1024 –        for (i = 0; vlist[i] != NULL; i++) {
1025 –                if (v == vlist[i])
1026 –                        return;
1027 –                if (v < vlist[i])
1028 –                        break;
1029 –        }
1030 –        for (j = i; vlist[j] != NULL; j++)
1031 –                ;
1032 –        while (j > i) {
1033 –                vlist[j] = vlist[j-1];
1034 –                --j;
1035 –        }
1036 –        vlist[i] = v;
1037 – }
1038 –
1039 – /* Sort triangle edges in standard order */
1040 – static void
1041 – order_triangle(MIGRATION *miga[3])
1042 – {
1043 –        RBFNODE         *vert[4];
1044 –        MIGRATION       *ord[3];
1045 –        int             i;
1046 –                                                /* order vertices, first */
1047 –        memset(vert, 0, sizeof(vert));
1048 –        for (i = 0; i < 3; i++) {
1049 –                insert_vert(vert, miga[i]->rbfv[0]);
1050 –                insert_vert(vert, miga[i]->rbfv[1]);
1051 –        }
1052 –                                                /* identify edge 0 */
1053 –        for (i = 0; i < 3; i++)
1054 –                if (miga[i]->rbfv[0] == vert[0] &&
1055 –                                miga[i]->rbfv[1] == vert[1]) {
1056 –                        ord[0] = miga[i];
1057 –                        break;
1058 –                }
1059 –                                                /* identify edge 1 */
1060 –        for (i = 0; i < 3; i++)
1061 –                if (miga[i]->rbfv[0] == vert[1] &&
1062 –                                miga[i]->rbfv[1] == vert[2]) {
1063 –                        ord[1] = miga[i];
1064 –                        break;
1065 –                }
1066 –                                                /* identify edge 2 */
1067 –        for (i = 0; i < 3; i++)
1068 –                if (miga[i]->rbfv[0] == vert[0] &&
1069 –                                miga[i]->rbfv[1] == vert[2]) {
1070 –                        ord[2] = miga[i];
1071 –                        break;
1072 –                }
1073 –        miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1074 – }
1075 –
1269   /* Partially advect between recorded incident angles and allocate new RBF */
1270   static RBFNODE *
1271   advect_rbf(const FVECT invec)
# Line 1086 | Line 1279 | advect_rbf(const FVECT invec)
1279  
1280          if (!get_interp(miga, invec))           /* can't interpolate? */
1281                  return(NULL);
1282 <        if (miga[1] == NULL)                    /* along edge? */
1282 >        if (miga[1] == NULL)                    /* advect along edge? */
1283                  return(e_advect_rbf(miga[0], invec));
1091 –                                                /* put in standard order */
1092 –        order_triangle(miga);
1284   #ifdef DEBUG
1285          if (miga[0]->rbfv[0] != miga[2]->rbfv[0] |
1286                          miga[0]->rbfv[1] != miga[1]->rbfv[0] |
# Line 1186 | Line 1377 | interp_isotropic()
1377          int             ix, ox, oy;
1378          FVECT           ivec, ovec;
1379          double          bsdf;
1380 <
1380 > #if DEBUG
1381 >        fprintf(stderr, "Writing isotropic order %d ", samp_order);
1382 >        if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1383 >        else fputs("raw data\n", stderr);
1384 > #endif
1385          if (pctcull >= 0) {                     /* begin output */
1386                  sprintf(cmd, "rttree_reduce -h -a -fd -r 3 -t %d -g %d",
1387                                  pctcull, samp_order);
# Line 1240 | Line 1435 | interp_anisotropic()
1435          int             ix, iy, ox, oy;
1436          FVECT           ivec, ovec;
1437          double          bsdf;
1438 <
1438 > #if DEBUG
1439 >        fprintf(stderr, "Writing anisotropic order %d ", samp_order);
1440 >        if (pctcull >= 0) fprintf(stderr, "data with %d%% culling\n", pctcull);
1441 >        else fputs("raw data\n", stderr);
1442 > #endif
1443          if (pctcull >= 0) {                     /* begin output */
1444                  sprintf(cmd, "rttree_reduce -h -a -fd -r 4 -t %d -g %d",
1445                                  pctcull, samp_order);

Diff Legend

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