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.13 by greg, Fri Sep 21 05:17:22 2012 UTC vs.
Revision 2.15 by greg, Sun Sep 23 16:45:20 2012 UTC

# Line 164 | Line 164 | static void
164   insert_dsf(RBFNODE *newrbf)
165   {
166          RBFNODE         *rbf, *rbf_last;
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 561 | 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 623 | 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) *
634 <                                        (from_rbf->nrbf*to_rbf->nrbf - 1));
635 <        double          *src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
636 <        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);
649 <                fprintf(stderr, "Building path from (theta,phi) (%d,%d) to ",
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)", 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 667 | Line 691 | make_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
691          while (total_rem > end_thresh) {
692                  total_rem -= migration_step(newmig, src_rem, dst_rem, pmtx);
693   #ifdef DEBUG
694 <                fputc('.', stderr);
695 < /*XXX*/break;
694 >                /* fputc('.', stderr); */
695 >                fprintf(stderr, "\n%.9f remaining...", total_rem);
696   #endif
697          }
698   #ifdef DEBUG
# Line 750 | Line 774 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
774          const MIGRATION *ej;
775          RBFNODE         *vother[2];
776          int             im_rev;
777 <                                        /* find shared edge in mesh */
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) {
# Line 764 | Line 788 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
788                          break;
789                  }
790          }
791 <        if (!get_triangles(vother, ej))
791 >        if (!get_triangles(vother, ej))         /* triangle on same side? */
792                  return(0);
793          return(vother[im_rev] != NULL);
794   }
# Line 773 | Line 797 | overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, c
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, area2, bestarea2 = FHUGE, 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);
789 <                if (dprod2 <= FTINY)
813 >                tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
814 >                if (DOT(vp, vmid) <= FTINY)
815                          continue;               /* wrong orientation */
816 <                area2 = DOT(vor, vor);
817 <                dprod2 *= dprod2 / area2;
818 <                if (dprod2 > bestdprod2 +
819 <                                FTINY*(1 - 2*(area2 < bestarea2)) &&
820 <                                !overlaps_tri(rbf0, rbf1, rbf)) {
821 <                        rbfbest = rbf;
822 <                        bestdprod2 = dprod2;
823 <                        bestarea2 = area2;
824 <                }
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(rbfbest);
830   }
# Line 807 | Line 835 | mesh_from_edge(MIGRATION *edge)
835   {
836          MIGRATION       *ej0, *ej1;
837          RBFNODE         *tvert[2];
838 +        
839 +        if (edge == NULL)
840 +                return;
841                                                  /* triangle on either side? */
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                          if (tvert[0] > edge->rbfv[0])
847 <                                ej0 = make_migration(edge->rbfv[0], tvert[0]);
847 >                                ej0 = make_migration(edge->rbfv[0], tvert[0], 1);
848                          else
849 <                                ej0 = make_migration(tvert[0], edge->rbfv[0]);
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]);
851 >                                ej1 = make_migration(edge->rbfv[1], tvert[0], 1);
852                          else
853 <                                ej1 = make_migration(tvert[0], edge->rbfv[1]);
853 >                                ej1 = make_migration(tvert[0], edge->rbfv[1], 1);
854                          mesh_from_edge(ej0);
855                          mesh_from_edge(ej1);
856                  }
857 <        }
827 <        if (tvert[1] == NULL) {                 /* grow mesh on left */
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                          if (tvert[1] > edge->rbfv[0])
861 <                                ej0 = make_migration(edge->rbfv[0], tvert[1]);
861 >                                ej0 = make_migration(edge->rbfv[0], tvert[1], 1);
862                          else
863 <                                ej0 = make_migration(tvert[1], edge->rbfv[0]);
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]);
865 >                                ej1 = make_migration(edge->rbfv[1], tvert[1], 1);
866                          else
867 <                                ej1 = make_migration(tvert[1], edge->rbfv[1]);
867 >                                ej1 = make_migration(tvert[1], edge->rbfv[1], 1);
868                          mesh_from_edge(ej0);
869                          mesh_from_edge(ej1);
870                  }
# Line 844 | Line 874 | mesh_from_edge(MIGRATION *edge)
874   #ifdef DEBUG
875   #include "random.h"
876   #include "bmpfile.h"
877 < /* Hash pointer to byte value */
877 > /* Hash pointer to byte value (must return 0 for NULL) */
878   static int
879   byte_hash(const void *p)
880   {
# Line 866 | Line 896 | write_edge_image(const char *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].r = random() & 0xff;
900 <                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 */
# Line 945 | Line 978 | build_mesh()
978          if (single_plane_incident) {
979                  for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
980                          if (rbf0->next != NULL)
981 <                                make_migration(rbf0, rbf0->next);
981 >                                make_migration(rbf0, rbf0->next, 1);
982                  return;
983          }
984                                                  /* start w/ shortest edge */
# Line 965 | Line 998 | build_mesh()
998          }
999                                                  /* build mesh from this edge */
1000          if (shrt_edj[0] < shrt_edj[1])
1001 <                mesh_from_edge(make_migration(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]));
1003 >                mesh_from_edge(make_migration(shrt_edj[1], shrt_edj[0], 0));
1004                                                  /* draw edge list into grid */
1005          draw_edges();
1006   }
# Line 1008 | 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 1033 | 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 1043 | 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 1122 | Line 1266 | memerr:
1266          return(NULL);   /* pro forma return */
1267   }
1268  
1125 – /* Insert vertex in ordered list */
1126 – static void
1127 – insert_vert(RBFNODE **vlist, RBFNODE *v)
1128 – {
1129 –        int     i, j;
1130 –
1131 –        for (i = 0; vlist[i] != NULL; i++) {
1132 –                if (v == vlist[i])
1133 –                        return;
1134 –                if (v < vlist[i])
1135 –                        break;
1136 –        }
1137 –        for (j = i; vlist[j] != NULL; j++)
1138 –                ;
1139 –        while (j > i) {
1140 –                vlist[j] = vlist[j-1];
1141 –                --j;
1142 –        }
1143 –        vlist[i] = v;
1144 – }
1145 –
1146 – /* Sort triangle edges in standard order */
1147 – static void
1148 – order_triangle(MIGRATION *miga[3])
1149 – {
1150 –        RBFNODE         *vert[4];
1151 –        MIGRATION       *ord[3];
1152 –        int             i;
1153 –                                                /* order vertices, first */
1154 –        memset(vert, 0, sizeof(vert));
1155 –        for (i = 0; i < 3; i++) {
1156 –                insert_vert(vert, miga[i]->rbfv[0]);
1157 –                insert_vert(vert, miga[i]->rbfv[1]);
1158 –        }
1159 –                                                /* identify edge 0 */
1160 –        for (i = 0; i < 3; i++)
1161 –                if (miga[i]->rbfv[0] == vert[0] &&
1162 –                                miga[i]->rbfv[1] == vert[1]) {
1163 –                        ord[0] = miga[i];
1164 –                        break;
1165 –                }
1166 –                                                /* identify edge 1 */
1167 –        for (i = 0; i < 3; i++)
1168 –                if (miga[i]->rbfv[0] == vert[1] &&
1169 –                                miga[i]->rbfv[1] == vert[2]) {
1170 –                        ord[1] = miga[i];
1171 –                        break;
1172 –                }
1173 –                                                /* identify edge 2 */
1174 –        for (i = 0; i < 3; i++)
1175 –                if (miga[i]->rbfv[0] == vert[0] &&
1176 –                                miga[i]->rbfv[1] == vert[2]) {
1177 –                        ord[2] = miga[i];
1178 –                        break;
1179 –                }
1180 –        miga[0] = ord[0]; miga[1] = ord[1]; miga[2] = ord[2];
1181 – }
1182 –
1269   /* Partially advect between recorded incident angles and allocate new RBF */
1270   static RBFNODE *
1271   advect_rbf(const FVECT invec)
# Line 1195 | Line 1281 | advect_rbf(const FVECT invec)
1281                  return(NULL);
1282          if (miga[1] == NULL)                    /* advect along edge? */
1283                  return(e_advect_rbf(miga[0], invec));
1198 –                                                /* put in standard order */
1199 –        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] |

Diff Legend

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