开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 3D图形编程 > 查看源码
TriMeshUtils.cc
资源名称：scanalyze-1.0.3_source_code.rar　[点击查看]
上传用户：kellyonhid
上传日期：2013-10-12
资源大小：932k
文件大小：43k
源码类别：

3D图形编程

开发平台：

Visual C++

TriMeshUtils.cc：源码内容
  1. //############################################################
  2. // 
  3. // TriMeshUtils.cc
  4. //
  5. // Utility functions for triangle meshes that are represented
  6. // as list of vertices and list of indices that define the
  7. // mesh connectivity.
  8. //
  9. //############################################################
  11. #include <stdio.h>
  12. #include <stdlib.h>
  13. #include <iostream.h>
  14. #include <vector.h>
  15. #include <hash_map.h>
  16. #include <hash_set.h>
  18. #include "TriMeshUtils.h"
  19. #include "Median.h"
  20. #include "ply++.h"
  21. #include "plvGlobals.h"
  22. #include "Progress.h"
  23. #include "Timer.h"
  27. // calculate vertex normals by averaging from triangle
  28. // normals (obtained from cross products)
  29. // possibly weighted with triangle areas
  30. void
  31. getVertexNormals(const vector<Pnt3> &vtx,
  32.  const vector<int>  &tri,
  33.  bool                strips,
  34.  vector<short>      &nrm_s,
  35.  int useArea)
  36. {
  37.   vector<Pnt3> nrm (vtx.size(), Pnt3());
  39.   int ntris = tri.size();
  40.   if (!useArea) {
  41.     if (strips) {
  42.       bool flip = false;
  43.       vector<int>::const_iterator i;
  44.       for (i = tri.begin(); i < tri.end(); i++) {
  45. if (*(i+2) == -1) {
  46.   // end of strip
  47.   i += 2;
  48.   flip = false;
  49.   continue;
  50. }
  52. int b, c;
  53. if (flip) {
  54.   b = *(i+2); c = *(i+1); flip = false;
  55. } else {
  56.   b = *(i+1); c = *(i+2); flip = true;
  57. }
  59. Pnt3 norm = normal(vtx[*i], vtx[b], vtx[c]);
  60. nrm[*i] += norm;
  61. nrm[b]  += norm;
  62. nrm[c]  += norm;
  63.       }
  64.     } else {
  65.       for (int i = 0; i < ntris; i += 3) {
  66. Pnt3 norm = normal(vtx[tri[i+0]], 
  67.    vtx[tri[i+1]], 
  68.    vtx[tri[i+2]]);
  69. nrm[tri[i+0]] += norm;
  70. nrm[tri[i+1]] += norm;
  71. nrm[tri[i+2]] += norm;
  72.       }
  73.     }
  75.   } else { // do use area norm
  77.     if (strips) {
  78.       bool flip = false;
  79.       for (int i = 0; i < ntris; i++) {
  80. if (tri[i+2] == -1) {
  81.   // end of strip
  82.   i += 2;
  83.   flip = false;
  84.   continue;
  85. }
  87. int a, b, c;
  88. if (flip) {
  89.   a = tri[i]; b = tri[i+2]; c = tri[i+1]; flip = false;
  90. } else {
  91.   a = tri[i]; b = tri[i+1]; c = tri[i+2]; flip = true;
  92. }
  94. //float area = .5 * cross(vtx[a], vtx[b], vtx[c]).norm();
  95. //Pnt3 areaNorm = normal(vtx[a], vtx[b], vtx[c]);
  96. //areaNorm *= area;
  97. Pnt3 areaNorm = cross(vtx[a], vtx[b], vtx[c]);
  98. nrm[a] += areaNorm;
  99. nrm[b] += areaNorm;
  100. nrm[c] += areaNorm;
  101.       }
  102.     } else {
  103.       for (int i = 0; i < ntris; i += 3) {
  105. /*
  106. float area = .5 * cross(vtx[tri[i+0]], vtx[tri[i+1]],
  107. vtx[tri[i+2]]).norm();
  108. Pnt3 areaNorm = normal(vtx[tri[i+0]], 
  109.        vtx[tri[i+1]], 
  110.        vtx[tri[i+2]]);
  111. areaNorm *= area;
  112. */
  113. Pnt3 areaNorm = cross(vtx[tri[i+0]], 
  114.       vtx[tri[i+1]], 
  115.       vtx[tri[i+2]]);
  116. nrm[tri[i+0]] += areaNorm;
  117. nrm[tri[i+1]] += areaNorm;
  118. nrm[tri[i+2]] += areaNorm;
  119.       }
  120.     }
  121.   }
  123.   nrm_s.clear();
  124.   nrm_s.reserve (3 * vtx.size());
  125.   for (Pnt3* n = nrm.begin(); n != nrm.end(); n++) {
  126.     n->set_norm(32767.0);
  127.     nrm_s.push_back (short(n->operator[](0)));
  128.     nrm_s.push_back (short(n->operator[](1)));
  129.     nrm_s.push_back (short(n->operator[](2)));
  130.   }
  131. }
  134. float
  135. median_edge_length(vector<Pnt3> &vtx,
  136.    vector<int>  &tri,
  137.    bool          strips,
  138.    int           percentile)
  139. {
  140.   int n = tri.size();
  141.   if (n < 3) return 0.0;
  143.   Median<float> med(percentile, n);
  144.     
  145.   if (strips) {
  146.     int *end = &tri[n];
  147.     // reuse the "inner edge" length in the strip
  148.     float inner_edge = dist2(vtx[0], vtx[1]);
  149.     for (int *i=&tri[2]; i<end; i++) {
  150.       if (*i == -1) {
  151. // skip over the start part of the next strip
  152. int cnt = 0;
  153. while (cnt < 3) {
  154.   i++;
  155.   if (i >= end) goto done;
  156.   if (*i == -1) cnt = 0;
  157.   else          cnt++;
  158. }
  159. inner_edge = dist2(vtx[*(i-1)], vtx[*(i-2)]);
  160.       }
  161.       med += inner_edge;
  162.       inner_edge = dist2(vtx[*(i-1)], vtx[*i]);
  163.       med += inner_edge;
  164.       med += dist2(vtx[*(i-2)], vtx[*i]);
  165.     }
  166.   } else {
  167.     // edge lengths
  168.     for (int i=0; i<n; i+=3) {
  169.       med += dist2(vtx[tri[i+0]], vtx[tri[i+1]]);
  170.       med += dist2(vtx[tri[i+2]], vtx[tri[i+1]]);
  171.       med += dist2(vtx[tri[i+0]], vtx[tri[i+2]]);
  172.     } 
  173.   }
  174.  done:
  175.   return sqrtf(med.find());
  176. }
  178. // find the median edge length
  179. // if a triangle has an edge that's longer than
  180. // factor times the median (or percentile), remove the triangle
  181. void 
  182. remove_stepedges(vector<Pnt3> &vtx,
  183.  vector<int>  &tri,
  184.  int           factor, 
  185.  int           percentile,
  186.  bool          strips)
  187. {
  188.   // calculate the threshold
  189.   float thr = median_edge_length(vtx, tri, strips, percentile);
  190.   thr *= thr * float(factor * factor);
  192.   int n = tri.size();
  193.   if (strips) {
  194.     vector<int> ntri;
  195.     ntri.reserve(n);
  196.     bool strip_on = false;
  197.     int cnt  = 0;
  198.     for (int i=0; i<n; i++) {
  199.       if (strip_on) {
  200. // strip is on
  201. bool end_found = (tri[i] == -1);
  202. bool end_here  = false;
  203. if (!end_found) {
  204.   // check for too long edges
  205.   if (cnt == 1) 
  206.     end_here = (dist2(vtx[tri[i]], vtx[tri[i-1]]) > thr);
  207.   else
  208.     end_here = (dist2(vtx[tri[i]], vtx[tri[i-1]]) > thr ||
  209. dist2(vtx[tri[i]], vtx[tri[i-2]]) > thr);
  210. }
  211. if (end_found || end_here) {
  212.   // found an end, clean up
  213.   if (cnt < 3) {
  214.     ntri.pop_back();
  215.     if (cnt == 2) ntri.pop_back();
  216.     if (ntri.size()) assert(ntri.back() == -1);
  217.   } else {
  218.     ntri.push_back(-1);
  219.   }
  220.   strip_on = false;
  221. }
  222. if (end_here) {
  223.   if (cnt % 2) {
  224.     // cnt odd, check if can add a single triangle
  225.     if (i > n-4) break; // almost done, break...
  226.     if (tri[i+1] != -1 &&
  227. tri[i+2] != -1 &&
  228. dist2(vtx[tri[i  ]], vtx[tri[i+1]]) < thr &&
  229. dist2(vtx[tri[i  ]], vtx[tri[i+2]]) < thr &&
  230. dist2(vtx[tri[i+1]], vtx[tri[i+2]]) < thr) {
  231.       ntri.push_back(tri[i]);
  232.       ntri.push_back(tri[i+1]);
  233.       ntri.push_back(tri[i+2]);
  234.       ntri.push_back(-1);
  235.     }
  236.   } else {
  237.     // cnt even, restart a strip
  238.     strip_on = true;
  239.     cnt = 1;
  240.     ntri.push_back(tri[i]);
  241.   }
  242. }
  243. if (!end_found && !end_here) {
  244.   // continue this strip
  245.   ntri.push_back(tri[i]);
  246.   cnt++;
  247. }
  248.       } else {
  249. // strip is not on
  250. if (tri[i] != -1) {
  251.   ntri.push_back(tri[i]);
  252.   strip_on = true;
  253.   cnt = 1;
  254. }
  255.       }
  256.       if (!strip_on) cnt = 0;
  257.     }
  258.     tri = ntri;
  260.   } else {
  261.     // could make faster by not recalculating edges!
  262.     for (int i=n-3; i>=0; i-=3) {
  263.       if (dist2(vtx[tri[i+0]], vtx[tri[i+1]]) > thr ||
  264.   dist2(vtx[tri[i+2]], vtx[tri[i+1]]) > thr ||
  265.   dist2(vtx[tri[i+0]], vtx[tri[i+2]]) > thr) {
  266. // remove the triangle
  267. tri[i+2] = tri.back(); tri.pop_back();
  268. tri[i+1] = tri.back(); tri.pop_back();
  269. tri[i+0] = tri.back(); tri.pop_back();
  270.       }
  271.     }
  272.   }
  273. }
  276. // check whether some vertices in vtx are not being used in tri
  277. // if so, remove them and also adjust the tri indices
  278. void
  279. remove_unused_vtxs(vector<Pnt3> &vtx,
  280.    vector<int>  &tri)
  281. {
  282.   // prepare vertex index map
  283.   vector<int> vtx_ind(vtx.size(), -1);
  284.   
  285.   // march through the triangles
  286.   // and mark the vertices that are actually used
  287.   int n = tri.size();
  288.   for (int i=0; i<n; i+=3) {
  289.     vtx_ind[tri[i+0]] = tri[i+0];
  290.     vtx_ind[tri[i+1]] = tri[i+1];
  291.     vtx_ind[tri[i+2]] = tri[i+2];
  292.   }
  294.   // remove the vertices that were not marked,
  295.   // also keep tab on how the indices change
  296.   int cnt = 0;
  297.   n = vtx.size();
  298.   for (i=0; i<n; i++) {
  299.     if (vtx_ind[i] != -1) {
  300.       vtx_ind[i] = cnt;
  301.       vtx[cnt] = vtx[i];
  302.       cnt++;
  303.     }
  304.   }
  305.   vtx.erase(&vtx[cnt], vtx.end());
  306.   // march through triangles and correct the indices
  307.   n = tri.size();
  308.   for (i=0; i<n; i+=3) {
  309.     tri[i+0] = vtx_ind[tri[i+0]];
  310.     tri[i+1] = vtx_ind[tri[i+1]];
  311.     tri[i+2] = vtx_ind[tri[i+2]];
  312.     assert(tri[i+0] >=0 && tri[i+0] < vtx.size());
  313.     assert(tri[i+1] >=0 && tri[i+1] < vtx.size());
  314.     assert(tri[i+2] >=0 && tri[i+2] < vtx.size());
  315.   }
  316. }
  319. // count the number of triangles in a tstrip
  320. int
  321. count_tris(const vector<int> &strips)
  322. {
  323.   int cnt = 0;
  324.   int i = 0;
  325.   while (i < strips.size()) {
  326.     // skip over the 2 first inds
  327.     i+=2;
  328.     while (strips[i++] != -1) {
  329.       cnt++;
  330.     }
  331.   }
  332.   return cnt;
  333. }
  337. struct _edge {
  338.   int a,b;
  339.   _edge(void) : a(0),b(0) {}
  340.   _edge(int _a, int _b)
  341.     {
  342.       if (_a < _b) { a=_a; b=_b; }
  343.       else         { b=_a; a=_b; }
  344.     }
  345.   void set(int _a, int _b)
  346.     {
  347.       if (_a < _b) { a=_a; b=_b; }
  348.       else         { b=_a; a=_b; }
  349.     }
  350.   bool operator==(const _edge &e) const
  351.     {
  352.       return a == e.a && b == e.b;
  353.     }
  354. };
  356. struct equal_edge
  357. {
  358.   bool operator()(const _edge &e1, const _edge &e2) const
  359.     {
  360.       return e1.a == e2.a && e1.b == e2.b;
  361.     }
  362. };
  364. struct hash_edge
  365. {
  366.   size_t operator()(const _edge &e) const
  367.     {
  368.       return e.a * 1000000 + e.b;
  369.     }
  370. };
  373. static void
  374. remove_this_triangle(vector<int> &ts, int i)
  375. {
  376.   for (int j=0; j<3; j++) {
  377.     int curr = ts[i+j]; // current neighbor
  378.     if (curr >= 0) {
  379.       if      (ts[curr]   == -2) continue; // nbor has been removed
  380.       // mark the correct edge of the neighbor to 
  381.       // not have a neighbor any more
  382.       if      (ts[curr]   == i) ts[curr]   = -1;
  383.       else if (ts[curr+1] == i) ts[curr+1] = -1;
  384.       else if (ts[curr+2] == i) ts[curr+2] = -1;
  385.     }
  386.     ts[i+j] = -2; // remove
  387.   }
  388. }
  391. static int Verbose = 1;
  392. static int Warnings = 1;
  394. static void
  395. old_tris_to_strips(const vector<int> &tris, 
  396.    vector<int>& tstripinds,
  397.    const char* name)
  398. {
  399.   int nTris = tris.size();
  400.   tstripinds.clear();
  401.   tstripinds.reserve(1.1*nTris/3);
  403.   // build a multimap: map edges to incident triangles
  404.   typedef hash_multimap<_edge,int,hash_edge,equal_edge> hmm;
  405.   hmm edges;
  406.   hmm::iterator edge_it;
  407.   int i,j;
  408.   _edge e;
  410.   cout << "Creating triangle strips from raw triangle data..." 
  411.        << endl;
  412.   if (Verbose)
  413.     printf ("input: %d vertices (%d triangles)n", nTris, nTris/3);
  415.   // show progress graphically
  416.   Progress progress (nTris * 3, "Create %s T-Strips", name);
  418.   for (i = 0; i < nTris; i+=3) {  // for all tris
  419.     if (i%300 == 0)
  420.       progress.update (i);
  422.     for (j = 0; j < 3; j++) {      // for all edges
  423.       // create an edge
  424.       e.set(tris[i+j], tris[i+(j+1)%3]);
  425.       // store the edge and starting vertex of the triangle
  426.       edges.insert(pair<const _edge,int>(e,i));
  427.     }
  428.   }
  429.   
  431.   // march through triangles
  432.   // store for each edge the neighboring triangle 
  433.   //  -1 if neighbor doesn't exist, 
  434.   //  later -2 if the triangle is removed/processed
  435.   vector<int> tri_nbors(nTris);
  436.   for (i = 0; i < nTris; i+=3) {
  437.     if (i%300 == 0)
  438.       progress.update (i + nTris);
  440.     for (j = 0; j < 3; j++) {
  441.       // find all instances of current edge
  442.       e.set(tris[i+j], tris[i+(j+1)%3]);
  443.       pair<hmm::iterator, hmm::iterator> p = edges.equal_range(e);
  445.       edge_it = p.first;
  446.       // shouldn't have more than two triangles sharing this edge:
  447.       p.first++;
  448.       if (p.first != p.second) {
  449. p.first++;
  450. if (p.first != p.second) {
  451.   if (Warnings) {
  452.     fprintf(stderr, "More than two pointers to same "
  453.     "edge (%d,%d)n", e.a, e.b);
  454.     while (Verbose && edge_it != p.second) {
  455.       int tri = (*edge_it).second;
  456.       fprintf(stderr, 
  457.       "t: offending tri: #%d (%d,%d,%d)n",
  458.       tri, tris[tri+0],
  459.       tris[tri+1], tris[tri+2]);
  460.       edge_it++;
  461.     }
  462.   }
  463.   // non-manifold edge, no neighbors across this
  464.   tri_nbors[i+j] = -1;
  465.   continue;
  466. }
  467.       }
  468.       // there's only 1 or 2 triangles for this edge
  469.       if ((*edge_it).second != i) { 
  470. // if the first match is not this triangle,
  471. // the neighbor across this edge is the current triangle
  472. tri_nbors[i+j] = (*edge_it).second;
  473.       } else {
  474. // the first of range is the current one
  475. edge_it++;
  476. if (edge_it == p.second) {
  477.   // that was the only edge, so this edge is not shared.
  478.   tri_nbors[i+j] = -1;
  479. } else {
  480.   // found another edge
  481.   // so the neighbor across this edge is now in .second
  482.   assert((*edge_it).first == e);
  483.   assert((*edge_it).second != i);
  484.   tri_nbors[i+j] = (*edge_it).second;
  485. }
  486.       }
  487.     }
  488.   }
  489.   
  490.   int nStrips = 0; //for stat-keeping
  492.   // start removing triangles from the tri_nbors
  493.   // and build the t-strips
  494.   int nNbors = tri_nbors.size();
  495.   for (i = 0; i < nNbors; i += 3) {
  496.     if (i%300 == 0)
  497.       progress.update (2* nTris + i);
  498.       
  499.     if (tri_nbors[i] == -2) continue; // has been removed
  500.     // try finding a neighbor
  501.     int nbor = -1;
  502.     for (j = 0; j < 3; j++) {
  503.       if (tri_nbors[i+j] >= 0) {
  504. nbor = j;
  505. break;
  506.       }
  507.     }
  508.     if (nbor == -1) {
  509.       // single triangle
  510.       tstripinds.push_back(tris[i+0]);
  511.       tstripinds.push_back(tris[i+1]);
  512.       tstripinds.push_back(tris[i+2]);
  513.     } else {
  514.       // initialize the tstrip
  515.       tstripinds.push_back(tris[i+(nbor+2)%3]);
  516.       tstripinds.push_back(tris[i+(nbor  )  ]);
  517.       tstripinds.push_back(tris[i+(nbor+4)%3]);
  518.       
  519.       int prev = i;
  520.       int dir = 1;
  521.       int curr = tri_nbors[prev+nbor];
  522.       while (curr != -1) {
  523. // find from which nbor we came to curr
  524. if      (tri_nbors[curr]   == prev) nbor = 0;
  525. else if (tri_nbors[curr+1] == prev) nbor = 1;
  526. else if (tri_nbors[curr+2] == prev) nbor = 2;
  527. else printf ("Warning: really lost!n");
  529. remove_this_triangle(tri_nbors, prev);
  530. tstripinds.push_back(tris[curr+(nbor+2)%3]);
  531. // find next nbor
  532. nbor = (nbor+1+dir)%3;
  533. dir ^= 1;
  534. prev = curr;
  535. curr = tri_nbors[prev+nbor];
  536.       }
  537.       remove_this_triangle(tri_nbors, prev);
  538.     }
  540.     ++nStrips;   //just for stat-keeping
  541.     tstripinds.push_back(-1);
  542.   }
  544.   if (Verbose)
  545.     printf ("output: %ld vertices (%d strips); avg. length %.1fn",
  546.   tstripinds.size()-nStrips, nStrips,
  547.   ((float)tstripinds.size()/nStrips) - 3);
  549.   printf ("Done.n");
  550. }
  559. typedef unsigned* adjacentfacelist;
  561. adjacentfacelist*
  562. TriMesh_FindAdjacentFaces (int numvertices,
  563.    const vector<int>& faces,
  564.    int*& numadjacentfaces)
  565. {
  566.   cout << "  Computing vtx->face mappings..." << flush;
  567.   int numfaces = faces.size();
  568.   assert (numvertices != 0);
  569.   assert (numfaces != 0);
  571.   // Step I - compute numadjacentfaces
  572.   numadjacentfaces = new int[numvertices];
  573.   memset(numadjacentfaces, 0, numvertices*sizeof(int));
  575.   for (int i = 0; i < numfaces; i++) {
  576.     numadjacentfaces[faces[i]]++;
  577.   }
  579.   // allocate one chunk of memory for all adjacent face lists
  580.   // total number of adjacent faces is numfaces
  581.   unsigned* adjacentfacedata = new unsigned [numfaces];
  582.   // this pointer will be incremented as needed but adjacentfaces[0]
  583.   // will always point to the beginning so it can later be freed
  584.   
  585.   // Step II - compute the actual vertex->tri lists...
  586.   adjacentfacelist* adjacentfaces = new adjacentfacelist[numvertices];
  587.   for (i = 0; i < numvertices; i++) {
  588.     adjacentfaces[i] = adjacentfacedata;
  589.     adjacentfacedata += numadjacentfaces[i];
  591.     //for (int j=0; j<numadjacentfaces[i]; j++)
  592.     //  adjacentfaces[i][j] = numfaces;
  593.   }
  595.   assert (adjacentfacedata == adjacentfaces[0] + numfaces);
  596.   for (unsigned* afdp = adjacentfaces[0]; afdp < adjacentfacedata; afdp++)
  597.     *afdp = numfaces;
  598.   
  599.   for (i = 0; i < numfaces; i++) {
  600.     unsigned *p = adjacentfaces[faces[i]];
  601.     while (*p != numfaces)
  602.       p++;
  603.     // snap to nearest multiple of 3, the start of tri data for that face
  604.     *p = (i/3) * 3;
  605.   }
  607.   return adjacentfaces;
  608. }
  613. #define FOR_EACH_VERTEX_OF_FACE(i,j) 
  614.   for (int jtmp = 0, j = faces[i]; 
  615.        jtmp < 3; 
  616.        jtmp++, j = faces[i + jtmp])
  618. #define FOR_EACH_ADJACENT_FACE(i,j) 
  619.   for (int jtmp=0, j = adjacentfaces[i][0]; 
  620.        jtmp < numadjacentfaces[i]; 
  621.        jtmp++, j = adjacentfaces[i][jtmp])
  623. static bool* done = NULL;
  624. static unsigned* stripsp = NULL;
  625. static int* numadjacentfaces = NULL;
  626. static adjacentfacelist* adjacentfaces = NULL;
  627. static const int* faces = NULL;
  628. static int nstrips = 0;
  629. static int nEvilTriangles;
  633. // Figure out the next triangle we're headed for...
  634. static inline int
  635. Tstrip_Next_Tri(unsigned tri, unsigned v1, unsigned v2)
  636. {
  637.   FOR_EACH_ADJACENT_FACE(v1, f1) {
  638.     if ((f1 == tri) || done[f1/3])
  639.       continue;
  640.     FOR_EACH_ADJACENT_FACE(v2, f2) {
  641.       if ((f2 == tri) || done[f2/3])
  642. continue;
  643.       if (f1 == f2)
  644. return f1;
  645.     }
  646.   }
  647.   
  648.   return -1;
  649. }
  651. // Build a whole strip of triangles, as long as possible...
  652. static void Tstrip_Crawl(unsigned v1, unsigned v2, unsigned v3,
  653.  unsigned next)
  654. {
  655.   // Insert the first tri...
  656.   *stripsp++ = v1;
  657.   *stripsp++ = v2;
  658.   *stripsp++ = v3;
  659.   
  660.   unsigned vlast1 = v3;
  661.   unsigned vlast2 = v2;
  663.   bool shouldbeflipped = true;
  665.   // Main loop...
  666.   do {
  667.     // Find the next vertex
  668.     int vnext;
  669.     FOR_EACH_VERTEX_OF_FACE(next,vnext_tmp) {
  670.       if ((vnext_tmp == vlast1) || (vnext_tmp == vlast2))
  671. continue;
  672.       vnext = vnext_tmp;
  673.       break;
  674.     }
  676.     bool thisflipped = true;
  677.     if ((faces[next+0] == vlast2) &&
  678. (faces[next+1] == vlast1) &&
  679. (faces[next+2] == vnext))
  680.       thisflipped = false;
  681.     if ((faces[next+2] == vlast2) &&
  682. (faces[next+0] == vlast1) &&
  683. (faces[next+1] == vnext))
  684.       thisflipped = false;
  685.     if ((faces[next+1] == vlast2) &&
  686. (faces[next+2] == vlast1) &&
  687. (faces[next+0] == vnext))
  688.       thisflipped = false;
  689.     
  690.     if (thisflipped != shouldbeflipped) {
  691.       if (nEvilTriangles-- > 0) {
  692. cerr << "Tstrip generation: inconsistent triangle orientation, "
  693.      << "tri " << next << endl;
  694.       }
  695.       goto bail;
  696.     }
  697.     
  698.     
  699.     // Record it
  701.     *stripsp++ = vnext;
  702.     vlast2 = vlast1;
  703.     vlast1 = vnext;
  704.     done[next/3] = true;
  705.     shouldbeflipped = !shouldbeflipped;
  706.     
  707.     // Try to find the next tri
  708.   } while ((next = Tstrip_Next_Tri(next, vlast1, vlast2)) != -1);
  709.   
  710.  bail:
  711.   // OK, done.  Mark end of strip
  712.   *stripsp++ = -1;
  713.   ++nstrips;
  714. }
  716. // Begin a tstrip, starting with triangle tri
  717. // tri is ordinal, not index (counts by 1)
  718. static void Tstrip_Bootstrap(unsigned tri)
  719. {
  720.   done[tri] = true;
  721.   
  722.   // Find two vertices with which to start.
  723.   // We do only a bit of lookahead, starting with vertices that will
  724.   // let us form a strip of length at least 2...
  725.   
  726.   tri *= 3;
  727.   unsigned vert1 = faces[tri];
  728.   unsigned vert2 = faces[tri+1];
  729.   unsigned vert3 = faces[tri+2];
  730.   
  731.   // Try vertices 1 and 2...
  732.   int nextface = Tstrip_Next_Tri(tri, vert1, vert2);
  733.   if (nextface != -1) {
  734.     Tstrip_Crawl(vert3, vert1, vert2, nextface);
  735.     return;
  736.   }
  737.   
  738.   // Try vertices 2 and 3...
  739.   nextface = Tstrip_Next_Tri(tri, vert2, vert3);
  740.   if (nextface != -1) {
  741.     Tstrip_Crawl(vert1, vert2, vert3, nextface);
  742.     return;
  743.   }
  744.   
  745.   // Try vertices 3 and 1...
  746.   nextface = Tstrip_Next_Tri(tri, vert3, vert1);
  747.   if (nextface != -1) {
  748.     Tstrip_Crawl(vert2, vert3, vert1, nextface);
  749.     return;
  750.   }
  751.   
  752.   // OK, nothing we can do. Do a single-triangle-long tstrip.
  753.   *stripsp++ = vert1;
  754.   *stripsp++ = vert2;
  755.   *stripsp++ = vert3;
  756.   *stripsp++ = -1;
  757.   ++nstrips;
  758. }
  761. static unsigned *Build_Tstrips(int numvertices,
  762.        const vector<int>& tris,
  763.        unsigned*& endstrips,
  764.        int& outNstrips)
  765. {
  766.   adjacentfaces = TriMesh_FindAdjacentFaces
  767.     (numvertices, tris, numadjacentfaces);
  769.   cout << " stripping... " << flush;
  770.   int numfaces = tris.size() / 3;
  771.   
  772.   // Allocate more than enough memory
  773.   unsigned* strips = new unsigned[4*numfaces+1];
  774.   stripsp = strips;
  775.   nEvilTriangles = 3;
  776.   
  777.   // Allocate array to record what triangles we've already done
  778.   done = new bool[numfaces];
  779.   memset(done, 0, numfaces*sizeof(bool));
  780.   faces = &tris[0];
  781.   nstrips = 0;
  782.   
  783.   // Build the tstrips
  784.   for (int i = 0; i < numfaces; i++) {
  785.     if (!done[i])
  786.       Tstrip_Bootstrap (i);
  787.   }
  788.   endstrips = stripsp;
  789.   outNstrips = nstrips;
  791.   if (nEvilTriangles < 0) {
  792.     cerr << "And there were " << -nEvilTriangles
  793.  << " more evil triangles for which no warnings were printed."
  794.  << endl << endl;
  795.   }
  796.   
  797.   // cleanup global arrays
  798.   delete [] done; done = NULL;
  799.   delete [] numadjacentfaces; numadjacentfaces = NULL;
  800.   delete [] adjacentfaces[0]; // ptr to one chunk of data for all
  801.   delete [] adjacentfaces; adjacentfaces = NULL;
  802.   stripsp = NULL;
  803.   faces = NULL;
  805.   cout << " done." << endl;
  806.   return strips;
  807. }
  810. void
  811. tris_to_strips(int numvertices,
  812.        const vector<int>& tris,
  813.        vector<int>& tstripinds,
  814.        const char* name)
  815. {
  816. #if 0
  817.   old_tris_to_strips (tris, tstripinds, name);
  818. #else
  819.   cout << "Tstripping " << tris.size()/3 << " triangles ("
  820.        << tris.size() << " vertices)..." << endl;
  822.   tstripinds.clear();
  823.   unsigned *strips, *end;
  824.   int nStrips = 0;
  826.   if (tris.size()) {
  827.     strips = Build_Tstrips (numvertices, tris, end, nStrips);
  829.     // this is a hair faster...
  830.     //tstripinds.reserve (end-strips);
  831.     //memcpy (&tstripinds[0], strips, 4*(end-strips));
  832.     
  833.     // but this is legal :)
  834.     //#ifdef __STL_MEMBER_TEMPLATES
  835. #if 0
  836.     tstripinds.assign (strips, end);
  837. #else
  838.     // as long as your compiler supports member templates :(
  839.     tstripinds.reserve (end - strips);
  840.     while (strips < end)
  841.       tstripinds.push_back (*strips++);
  842. #endif
  843.   }
  845.   cout << "Tstrip results: " << nstrips << " strips ("
  846.        << tstripinds.size() - nStrips << " vertices, avg. length "
  847.        << ((float)tstripinds.size()/nStrips) - 3 << ")." << endl;
  848. #endif
  849. }
  855. void 
  856. strips_to_tris(const vector<int> &tstrips, 
  857.        vector<int>& tris, int nTris)
  858. {
  859.   if (!tstrips.size())
  860.      return;
  862.   assert(tstrips.back() == -1);
  863.   tris.clear();
  864.   if (!nTris)
  865.     nTris = tstrips.size() * 3;  // estimate
  866.   tris.reserve (nTris);
  868.   vector<int>::const_iterator vert;
  869.   for (vert = tstrips.begin(); vert != tstrips.end(); vert++) {
  870.     while (*vert == -1) { // handle 0-length strips
  871.       ++vert;
  872.       if (vert == tstrips.end()) break;
  873.     }
  874.     if (vert == tstrips.end()) break;
  876.     vert += 2;   //we'll look backwards at these
  878.     int dir = 0;
  879.     while (*vert != -1)  {
  880.       tris.push_back(vert[-2 + dir]);
  881.       tris.push_back(vert[-1 - dir]);
  882.       tris.push_back(vert[0]);
  883.       vert++;
  884.       dir ^= 1;
  885.     }
  886.   }
  887. }
  890. void
  891. flip_tris(vector<int> &inds, bool stripped)
  892. {
  893.   int n = inds.size();
  894.   
  895.   if (stripped) {
  896.     int i=0;
  897.     while (i<n) {
  898.       i++; // skip the first one
  899.       // after that, swap all complete pairs until -1
  900.       while (1) {
  901. int tmp = inds[i++];
  902. // after break, inds[i-1] == -1
  903. if (tmp == -1) break;
  904. if (inds[i] == -1) { i++; break; }
  905. inds[i-1] = inds[i];
  906. inds[i++] = tmp;
  907.       }
  908.     }
  909.   } else {
  910.     for (int i = 0; i < n; i += 3)
  911.       swap(inds[i], inds[i+1]);
  912.   }
  913. }
  917. #define BEFORE(mod, var) 
  918.   ((mod == 0) ? tris[var + 2] : tris[var - 1])
  920. #define AFTER(mod, var) 
  921.   ((mod == 2) ? tris[var - 2] : tris[var + 1])
  923. // sort a ring of neighbors
  925. static void
  926. sort_nbors(int *nbors, int n, const  vector<int> &tris)
  927. {
  928.   for (int i=0; i<n; i++) {
  929.     int vb = BEFORE(nbors[i]%3, nbors[i]);
  930.     for (int j=i+2; j<n; j++) {
  931.       if (AFTER(nbors[j]%3, nbors[j]) == vb) {
  932. // swap
  933. int tmp    = nbors[i+1];
  934. nbors[i+1] = nbors[j];
  935. nbors[j]   = tmp;
  936. break;
  937.       }
  938.     }
  939.   }
  940. }
  942. static void
  943. nbors_tris(int    n_vtx,
  944.    const  vector<int> &tris,
  945.    int*  &n_nbors, // number of neighbors for vtx i
  946.    int** &nbors,   // the actual neighbors
  947.    bool   sorted = false)
  948. {
  949.   // Step I - compute numadjacentfaces
  950.   n_nbors = new int[n_vtx];
  951.   memset(n_nbors, 0, n_vtx*sizeof(int));
  953.   int n_tri_inds = tris.size();
  954.   for (int i=0; i<n_tri_inds; i++) {
  955.     n_nbors[tris[i]]++;
  956.   }
  958.   // allocate one chunk of memory for all neighbor triangle lists
  959.   // total number of neighboring triangles is n_tri_inds
  960.   int *data = new int[n_tri_inds];
  961.   // this pointer will be incremented as needed but nbors[0]
  962.   // will always point to the beginning so it can later be freed
  963.   
  964.   // Step II - compute the actual vertex->tri lists...
  965.   nbors = new int*[n_vtx];
  966.   for (i=0; i<n_vtx; i++) {
  967.     nbors[i] = data;
  968.     data += n_nbors[i];
  969.   }
  971.   assert (data == nbors[0] + n_tri_inds);
  972.   // mark as empty
  973.   int *p;
  974.   for (p = nbors[0]; p<data; p++) {
  975.     *p = -1;
  976.   }
  978.   for (i=0; i<n_tri_inds; i++) {
  979.     int *p = nbors[tris[i]];
  980.     while (*p != -1) p++;
  981.     *p = i;
  982.   }
  984.   if (sorted) {
  985.     for (i=0; i<n_vtx; i++) {
  986.       sort_nbors(nbors[i], n_nbors[i], tris);
  987.     }
  988.   }
  989. }
  993. // assume bdry has the right size and has been initialized with
  994. // zeroes
  995. void
  996. mark_boundary_verts(vector<char> &bdry,
  997.     const vector<int> &tris)
  998. {
  999.   cout << "Marking boundary vertices ... " << flush;
  1000.   int nv = bdry.size();
  1001.   int *n_nbors, **nbors;
  1002.   nbors_tris(nv, tris, n_nbors, nbors, true);
  1004.   // for each vertex, try to find a full loop 
  1005.   // around it, if can't its boundary
  1006.   int va, vb; // vertex after, vertex before
  1007.   for (int i=0; i<nv; i++) {
  1008.     if (n_nbors[i] <= 1) {
  1009.       bdry[i] = 1;
  1010.       continue;
  1011.     }
  1012.     int jend = n_nbors[i] - 1;
  1013.     for (int j=0; j<jend; j++) {
  1014.       // which vtx comes before this (corner j)?
  1015.       int k = nbors[i][j];
  1016.       vb = BEFORE(k % 3, k);
  1017.       // which vtx comes after this (corner j+1)?
  1018.       k = nbors[i][j+1];
  1019.       va = AFTER(k % 3, k);
  1020.       if (va != vb) {
  1021. // not a continuous neighbor chain
  1022. bdry[i] = 1;
  1023. break;
  1024.       }
  1025.     }
  1026.     if (bdry[i] == 0) {
  1027.       // check the last possible link (from end to start)
  1028.       // which vtx comes before this (corner jend)?
  1029.       int k = nbors[i][jend];
  1030.       vb = BEFORE(k % 3, k);
  1031.       // which vtx comes after this (corner jstart)?
  1032.       k = nbors[i][0];
  1033.       va = AFTER(k % 3, k);
  1034.       if (va != vb) {
  1035. // not a continuous neighbor chain
  1036. bdry[i] = 1;
  1037.       }
  1038.     }      
  1039.   }
  1041.   // cleanup
  1042.   delete[] n_nbors;
  1043.   delete[] nbors[0]; // ptr to one chunk of data for all
  1044.   delete[] nbors;
  1045.   cout << "done" << endl;
  1046. }
  1049. typedef vector< vector<int> > vecvec;
  1050. typedef vector<int>::const_iterator vicit;
  1051. static void
  1052. build_nbor_map(vecvec &nlist, const vector<int> &tris)
  1053. {
  1054.   // build an edge set
  1055.   typedef hash_set<_edge,hash_edge,equal_edge> hs;
  1056.   hs edges;
  1057.   int n = tris.size();
  1058.   for (int i=0; i<n; i+=3) {
  1059.     edges.insert(_edge(tris[i+0], tris[i+1]));
  1060.     edges.insert(_edge(tris[i+1], tris[i+2]));
  1061.     edges.insert(_edge(tris[i+2], tris[i+0]));
  1062.   }
  1064.   // for every edge, create an entry to neighborhood map
  1065.   hs::const_iterator edge_it;
  1066.   for (edge_it = edges.begin(); edge_it != edges.end(); edge_it++) {
  1067.     nlist[(*edge_it).a].push_back((*edge_it).b);
  1068.     nlist[(*edge_it).b].push_back((*edge_it).a);
  1069.   }
  1070. }
  1072. #if 0
  1074. // assume bdry has the right size and has been initialized with
  1075. // zeroes
  1076. void
  1077. distance_from_boundary(vector<float> &distances,
  1078.        const vector<Pnt3> &pnts,
  1079.        const vector<int> &tris)
  1080. {
  1081.   TIMER(dist);
  1083.   int n = pnts.size();
  1085.   // find boundary verts
  1086.   vector<char> bdry(n, 0);
  1087.   mark_boundary_verts(bdry, tris);
  1089.   // initialize distance vector
  1090.   distances.clear();
  1091.   distances.insert(distances.begin(), n, 1e33);
  1092.   int nbdry = 0;
  1093.   for (int i=0; i<n; i++) {
  1094.     if (bdry[i]) { distances[i] = 0.0; nbdry++; }
  1095.   }
  1097.   // build neighborhood map
  1098.   vecvec nmap(n);
  1099.   build_nbor_map(nmap, tris);
  1101.   // find the neighbors of bdry vertices
  1102.   hash_set<int,hash<int>,equal_to<int> > prevset, workset, nextset;
  1103.   hash_set<int,hash<int>,equal_to<int> >::const_iterator hcit;
  1104.   for (i=0; i<n; i++) {
  1105.     if (bdry[i]) {
  1106.       vicit end = nmap[i].end();
  1107.       for (vicit it = nmap[i].begin(); it != end; it++) {
  1108. if (!bdry[*it]) {
  1109.   workset.insert(*it);
  1110. }
  1111.       }
  1112.       prevset.insert(i);
  1113.     }
  1114.   }
  1116.   // first round: just go through the mesh in the
  1117.   // topological order in increasing distance from the
  1118.   // mesh boundary
  1119.   int cnt = nbdry;
  1120.   cout << pnts.size() << "..." << flush;
  1121.   while (!workset.empty()) {
  1122.     
  1123.     cout << pnts.size() - cnt << "..." << flush;
  1124.     cnt += workset.size();
  1126.     for (hcit = workset.begin(); hcit != workset.end(); hcit++) {
  1127.       i = *hcit;
  1128.       vicit end = nmap[i].end();
  1129.       for (vicit it = nmap[i].begin(); it != end; it++) {
  1130. // calculate new distance
  1131. float dij = dist(pnts[i], pnts[*it]);
  1132. float d   = dij + distances[*it];
  1133. if (d < distances[i]) distances[i] = d;
  1134. // should neighbor go to the nextset?
  1135. if (workset.find(*it) != workset.end()) continue;
  1136. if (prevset.find(*it) != prevset.end()) continue;
  1137. nextset.insert(*it);
  1138. d = dij + distances[i];
  1139. if (d < distances[*it]) distances[*it] = d;
  1140.       }
  1141.     }
  1143.     prevset = workset;
  1144.     workset = nextset;
  1145.     nextset.clear();
  1146.   }
  1148.   cout << endl << "Restart" << endl;
  1150.   // second round: put all the nonboundary vertices into
  1151.   // the workset, all the 
  1153.   workset.clear();
  1154.   for (i=0; i<n; i++) {
  1155.     if (!bdry[i]) {
  1156.       workset.insert(i);
  1157.     }
  1158.   }
  1160.   while (!workset.empty()) {
  1161.     
  1162.     cout << workset.size() << "..." << flush;
  1164.     nextset.clear();
  1165.     // calculate distances
  1166.     for (hcit = workset.begin(); hcit != workset.end(); hcit++) {
  1167.       i = *hcit;
  1168.       vicit end = nmap[i].end();
  1169.       for (vicit it = nmap[i].begin(); it != end; it++) {
  1170. float d = dist(pnts[i], pnts[*it]) + distances[*it];
  1171. if (d < distances[i]) distances[i] = d;
  1172.       }
  1173.     }
  1175.     // see if neighbors should be in set
  1176.     for (hcit = workset.begin(); hcit != workset.end(); hcit++) {
  1177.       i = *hcit;
  1178.       vicit end = nmap[i].end();
  1179.       for (vicit it = nmap[i].begin(); it != end; it++) {
  1180. float d = dist(pnts[i], pnts[*it]) + distances[i];
  1181. if (d < distances[*it]) {
  1182.   distances[*it] = d;
  1183.   nextset.insert(*it);
  1184. }
  1185.       }
  1186.     }
  1187.     workset = nextset;
  1188.   }
  1189.   cout << endl;
  1190. }
  1191. #else
  1193. // assume bdry has the right size and has been initialized with
  1194. // zeroes
  1195. void
  1196. distance_from_boundary(vector<float> &distances,
  1197.        const vector<Pnt3> &pnts,
  1198.        const vector<int> &tris)
  1199. {
  1200.   TIMER(dist);
  1202.   int n = pnts.size();
  1204.   // find boundary verts
  1205.   vector<char> bdry(n, (char)0);
  1207.   int *n_nbors, **nbors;
  1208.   nbors_tris(n, tris, n_nbors, nbors, true);
  1210.   // for each vertex, try to find a full loop 
  1211.   // around it, if can't its boundary
  1212.   int va, vb; // vertex after, vertex before
  1213.   for (int i=0; i<n; i++) {
  1214.     if (n_nbors[i] <= 1) {
  1215.       bdry[i] = 1;
  1216.       continue;
  1217.     }
  1218.     int jend = n_nbors[i] - 1;
  1219.     for (int j=0; j<jend; j++) {
  1220.       // which vtx comes before this (corner j)?
  1221.       int k = nbors[i][j];
  1222.       vb = BEFORE(k % 3, k);
  1223.       // which vtx comes after this (corner j+1)?
  1224.       k = nbors[i][j+1];
  1225.       va = AFTER(k % 3, k);
  1226.       if (va != vb) {
  1227. // not a continuous neighbor chain
  1228. bdry[i] = 1;
  1229. break;
  1230.       }
  1231.     }
  1232.     if (bdry[i] == 0) {
  1233.       // check the last possible link (from end to start)
  1234.       // which vtx comes before this (corner jend)?
  1235.       int k = nbors[i][jend];
  1236.       vb = BEFORE(k % 3, k);
  1237.       // which vtx comes after this (corner jstart)?
  1238.       k = nbors[i][0];
  1239.       va = AFTER(k % 3, k);
  1240.       if (va != vb) {
  1241. // not a continuous neighbor chain
  1242. bdry[i] = 1;
  1243.       }
  1244.     }      
  1245.   }
  1247.   cout << "done boundary" << endl;
  1249.   // initialize distance vector
  1250.   distances.clear();
  1251.   distances.insert(distances.begin(), n, 1e33);
  1252.   int nbdry = 0;
  1253.   for (i=0; i<n; i++) {
  1254.     if (bdry[i]) { distances[i] = 0.0; nbdry++; }
  1255.   }
  1257.   // find the neighbors of bdry vertices
  1258.   hash_set<int,hash<int>,equal_to<int> > prevset, workset, nextset;
  1259.   hash_set<int,hash<int>,equal_to<int> >::const_iterator hcit;
  1260.   for (i=0; i<n; i++) {
  1261.     if (bdry[i]) {
  1262.       for (int j=0; j<n_nbors[i]; j++) {
  1263. int k = AFTER(nbors[i][j]%3, nbors[i][j]);
  1264. if (!bdry[k]) workset.insert(k);
  1265. k = BEFORE(nbors[i][j]%3, nbors[i][j]);
  1266. if (!bdry[k]) workset.insert(k);
  1267.       }
  1268.       prevset.insert(i);
  1269.     }
  1270.   }
  1272.   // first round: just go through the mesh in the
  1273.   // topological order in increasing distance from the
  1274.   // mesh boundary
  1275.   int cnt = nbdry;
  1276.   cout << pnts.size() << "..." << flush;
  1277.   while (!workset.empty()) {
  1278.     
  1279.     cout << pnts.size() - cnt << "..." << flush;
  1280.     cnt += workset.size();
  1282.     for (hcit = workset.begin(); hcit != workset.end(); hcit++) {
  1283.       i = *hcit;
  1284.       for (int j=0; j<n_nbors[i]; j++) {
  1285. int k = AFTER(nbors[i][j]%3, nbors[i][j]);
  1286. // calculate new distance
  1287. float dij = dist(pnts[i], pnts[k]);
  1288. float d   = dij + distances[k];
  1289. if (d < distances[i]) distances[i] = d;
  1290. // should neighbor go to the nextset?
  1291. if (workset.find(k) != workset.end()) continue;
  1292. if (prevset.find(k) != prevset.end()) continue;
  1293. nextset.insert(k);
  1294. d = dij + distances[i];
  1295. if (d < distances[k]) distances[k] = d;
  1296.       }
  1297.     }
  1299.     prevset = workset;
  1300.     workset = nextset;
  1301.     nextset.clear();
  1302.   }
  1304.   cout << endl << "Restart" << endl;
  1306.   // second round: put all the nonboundary vertices into
  1307.   // the workset, all the 
  1309.   workset.clear();
  1310.   for (i=0; i<n; i++) {
  1311.     if (!bdry[i]) {
  1312.       workset.insert(i);
  1313.     }
  1314.   }
  1316.   while (!workset.empty()) {
  1317.     
  1318.     cout << workset.size() << "..." << flush;
  1320.     nextset.clear();
  1321.     // calculate distances
  1322.     for (hcit = workset.begin(); hcit != workset.end(); hcit++) {
  1323.       i = *hcit;
  1324.       for (int j=0; j<n_nbors[i]; j++) {
  1325. int k = AFTER(nbors[i][j]%3, nbors[i][j]);
  1326. float d = dist(pnts[i], pnts[k]) + distances[k];
  1327. if (d < distances[i]) distances[i] = d;
  1328.       }
  1329.     }
  1331.     // see if neighbors should be in set
  1332.     for (hcit = workset.begin(); hcit != workset.end(); hcit++) {
  1333.       i = *hcit;
  1334.       for (int j=0; j<n_nbors[i]; j++) {
  1335. int k = AFTER(nbors[i][j]%3, nbors[i][j]);
  1336. float d = dist(pnts[i], pnts[k]) + distances[i];
  1337. if (d < distances[k]) {
  1338.   distances[k] = d;
  1339.   nextset.insert(k);
  1340. }
  1341.       }
  1342.     }
  1343.     workset = nextset;
  1344.   }
  1345.   cout << endl;
  1346.   // cleanup
  1347.   delete[] n_nbors;
  1348.   delete[] nbors[0]; // ptr to one chunk of data for all
  1349.   delete[] nbors;
  1350. }
  1351. #endif
  1354. void
  1355. quadric_simplify(const vector<Pnt3> &vtx_in,
  1356.  const vector<int>  &tri_in,
  1357.  vector<Pnt3> &vtx_out,
  1358.  vector<int>  &tri_out,
  1359.  int           goal,
  1360.  int           optLevel,
  1361.  float         errLevel,
  1362.  float         boundWeight)
  1363. {
  1364.    char inName[L_tmpnam], outName[L_tmpnam];
  1365.    FILE *f;
  1366.    int i;
  1367.    char ins[200];
  1369.    vtx_out.clear();
  1370.    tri_out.clear();
  1372.    // Write out the original mesh data to a temporary input file
  1374.    if (!tmpnam(inName))  {
  1375.       cerr << "ERROR: Unable to get temporary filename for QSlim input" << endl;
  1376.       return;
  1377.    }
  1378.    if (!(f = fopen(inName, "w")))  {
  1379.       cerr << "ERROR: Unable to open file " << inName << " for QSlim input"
  1380.            << endl;
  1381.       return;
  1382.    }
  1383.    for (i=0; i<vtx_in.size(); ++i) 
  1384.       fprintf(f, "v %f %f %fn", vtx_in[i][0], vtx_in[i][1], vtx_in[i][2]);
  1385.    for (i=0; i<tri_in.size(); i+=3) 
  1386.       fprintf(f, "f %d %d %dn", tri_in[i]+1, tri_in[i+1]+1, tri_in[i+2]+1);
  1387.    fclose(f);
  1389.    // Get a filename for the temporary output file
  1390.    
  1391.    if (!tmpnam(outName))  {
  1392.       cerr << "ERROR: Unable to get temporary filename for QSlim output"
  1393.            << endl;
  1394.       remove(inName);
  1395.       return;
  1396.    }
  1398.    // Run QSlim as an external program to simplify the mesh
  1400.    char qslimCmd[200];
  1401.    sprintf(qslimCmd, "qslim -t %d -B %f -o %s %s", goal, boundWeight,
  1402.                                                    outName, inName);
  1403.    if (system(qslimCmd))  {
  1404.       cerr << "nERROR: Unable to execute call to external QSlim program:n"
  1405.            << qslimCmd << "nPlease verify that qslim is properly installed"
  1406.            << " and accessible on your system." << endl;
  1407.       remove(inName);
  1408.       remove(outName);
  1409.       return;
  1410.    }
  1412.    // Read in the simplified mesh from the resulting output file
  1414.    float v1, v2, v3;
  1415.    int f1, f2, f3;
  1417.    if (!(f = fopen(outName, "r"))) {
  1418.       cerr << "ERROR: Unable to open file " << outName << " from QSlim"
  1419.            << endl;
  1420.       remove(inName);
  1421.       return;
  1422.    }
  1423.    while (fscanf(f, "%s", ins) == 1)  {
  1424.       if (!strcmp(ins, "v")) {
  1425.          fscanf(f, "%f %f %f", &v1, &v2, &v3);
  1426.          vtx_out.push_back(Pnt3(v1, v2, v3));
  1427.       }
  1428.       else if (!strcmp(ins, "f")) {
  1429.          fscanf(f, "%d %d %d", &f1, &f2, &f3);
  1430.          tri_out.push_back(f1-1);
  1431.          tri_out.push_back(f2-1);
  1432.          tri_out.push_back(f3-1);
  1433.       }
  1434.    }
  1435.    fclose(f);
  1437.    remove(inName);
  1438.    remove(outName);
  1439. }
  1443. /////////////////////////////////////
  1444. // BEGIN write_ply_file() HELPER CODE
  1446. struct PlyVertex {
  1447.   float x, y, z;
  1448.   float nx, ny, nz;
  1449.   uchar diff_r, diff_g, diff_b;
  1450.   float tex_u, tex_v;
  1451.   float intensity;
  1452.   float std_dev;
  1453.   float confidence;
  1454. };
  1456. const int MAX_FACE_VERTS = 100;
  1458. struct PlyFace {
  1459.   uchar nverts;
  1460.   int verts[MAX_FACE_VERTS];
  1461. };
  1464. struct CoordConfVert {
  1465.    float v[3];
  1466.    float conf;
  1467. };
  1469. struct TriLocal {
  1470.    uchar nverts;
  1471.    int index[3];
  1472. };
  1475. static PlyProperty vert_prop_x =  
  1476.    {"x", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1477. static PlyProperty vert_prop_y =  
  1478.   {"y", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1479. static PlyProperty vert_prop_z =  
  1480.   {"z", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1481. static PlyProperty vert_prop_nx =  
  1482.    {"nx", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1483. static PlyProperty vert_prop_ny =  
  1484.   {"ny", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1485. static PlyProperty vert_prop_nz =  
  1486.   {"nz", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1487. static PlyProperty vert_prop_intens =  
  1488.   {"intensity", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1489. static PlyProperty vert_prop_conf =
  1490.   {"confidence", PLY_FLOAT, PLY_FLOAT, 0, 0, PLY_START_TYPE,
  1491.    PLY_START_TYPE, 0};
  1494. static PlyProperty vert_prop_diff_r =  
  1495.   {"diffuse_red", PLY_UCHAR, PLY_UCHAR, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1496. static PlyProperty vert_prop_diff_g =  
  1497.   {"diffuse_green", PLY_UCHAR, PLY_UCHAR, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1498. static PlyProperty vert_prop_diff_b =  
  1499.   {"diffuse_blue", PLY_UCHAR, PLY_UCHAR, 0, 0, PLY_START_TYPE, PLY_START_TYPE, 0};
  1501. static PlyProperty face_props[] = { 
  1502.   {"vertex_indices", PLY_INT, PLY_INT, 0, true, PLY_UCHAR, PLY_UCHAR, 0},
  1503. };
  1505. struct tstrip_info {
  1506.   int nverts;
  1507.   const int* vertData;
  1508. };
  1510. static PlyProperty tstrips_props[] = {
  1511.   {"vertex_indices", PLY_INT, PLY_INT, 4, true, PLY_INT, PLY_INT, 0},
  1512. };
  1514. static void
  1515. set_offsets(void)
  1516. {
  1517.   static int first = 1;
  1518.   if (first) {
  1519.     first = 0;
  1520.     vert_prop_x.offset          = offsetof(PlyVertex,x);
  1521.     vert_prop_y.offset          = offsetof(PlyVertex,y);
  1522.     vert_prop_z.offset          = offsetof(PlyVertex,z);
  1523.     vert_prop_nx.offset         = offsetof(PlyVertex,nx);
  1524.     vert_prop_ny.offset         = offsetof(PlyVertex,ny);
  1525.     vert_prop_nz.offset         = offsetof(PlyVertex,nz);
  1526.     vert_prop_intens.offset     = offsetof(PlyVertex,intensity);
  1527.     vert_prop_conf.offset       = offsetof(PlyVertex,confidence);
  1528.     vert_prop_diff_r.offset     = offsetof(PlyVertex,diff_r);
  1529.     vert_prop_diff_g.offset     = offsetof(PlyVertex,diff_g);
  1530.     vert_prop_diff_b.offset     = offsetof(PlyVertex,diff_b);
  1532.     tstrips_props[0].offset     = offsetof(tstrip_info,vertData);
  1533.   }
  1534. }
  1536. ////////////////////////////
  1537. // write_ply_file() : called by all exported wrappers to handle the different
  1538. // variants of ply files.
  1540. static void
  1541. write_ply_file(const char *fname,
  1542.        const vector<Pnt3> &vtx,
  1543.        const vector<int> &tris, bool strips,
  1544.        const vector<Pnt3> &nrm,
  1545.        const vector<uchar> &intensity,
  1546.        const vector<float> &confidence,
  1547.        bool hasNrm, bool hasIntensity, bool hasConfidence)
  1548. {
  1549.   int i;
  1550.   char *elem_names[] = {"vertex", strips ? "tristrips" : "face"};
  1552.   bool hasColor = false;
  1553.   bool hasAlpha = false;
  1554.   if (hasIntensity) {
  1555.     if (intensity.size() == 3 * vtx.size()) {
  1556.       hasColor = true;
  1557.       hasIntensity = false;
  1558.     } else if (intensity.size() == 4 * vtx.size()) {
  1559.       hasColor = true;
  1560.       hasAlpha = true;
  1561.       hasIntensity = false;
  1562.     }
  1563.   }
  1565.   PlyFile ply;
  1566.   if (!ply.open_for_writing(fname, 2, elem_names,PLY_BINARY_BE))
  1567.     return;
  1569.   vector<PlyProperty> vert_props;
  1570.   set_offsets();
  1572.   vert_props.push_back(vert_prop_x);
  1573.   vert_props.push_back(vert_prop_y);
  1574.   vert_props.push_back(vert_prop_z);
  1576.   if (hasNrm) {
  1577.     vert_props.push_back(vert_prop_nx);
  1578.     vert_props.push_back(vert_prop_ny);
  1579.     vert_props.push_back(vert_prop_nz);
  1580.   }
  1582.   if (hasIntensity) {
  1583.     vert_props.push_back(vert_prop_intens);
  1584.   }
  1585.     
  1586.   if (hasConfidence) {
  1587.     vert_props.push_back(vert_prop_conf);
  1588.   }
  1590.   if (hasColor) {
  1591.     vert_props.push_back(vert_prop_diff_r);
  1592.     vert_props.push_back(vert_prop_diff_g);
  1593.     vert_props.push_back(vert_prop_diff_b);
  1594.   }
  1595.     
  1596.   // count offset
  1597.   face_props[0].offset = offsetof(PlyFace, verts);
  1598.   face_props[0].count_offset = offsetof(PlyFace, nverts);  
  1599.     
  1600.   ply.describe_element ("vertex", vtx.size(), 
  1601. vert_props.size(), &vert_props[0]);
  1602.   if (strips)
  1603.     ply.describe_element ("tristrips", 1, 1, tstrips_props);
  1604.   else
  1605.     ply.describe_element ("face", tris.size()/3, 1, face_props);
  1606.   ply.header_complete();
  1607.     
  1608.     // set up and write the vertex elements
  1609.   PlyVertex plyVert;
  1611.   ply.put_element_setup ("vertex");
  1613.   for (i = 0; i < vtx.size(); i++) {
  1614.     plyVert.x = vtx[i][0];
  1615.     plyVert.y = vtx[i][1];
  1616.     plyVert.z = vtx[i][2];
  1617.     if (hasIntensity)
  1618.       plyVert.intensity = intensity[i] / 255.0;
  1620.     if (hasColor) {
  1621.       int ci = (hasAlpha ? 4 : 3) * i;
  1622.       plyVert.diff_r = intensity[ci  ];
  1623.       plyVert.diff_g = intensity[ci+1];
  1624.       plyVert.diff_b = intensity[ci+2];
  1625.     }
  1627.     if (hasNrm) {
  1628.       plyVert.nx = nrm[i][0];
  1629.       plyVert.ny = nrm[i][1];
  1630.       plyVert.nz = nrm[i][2];
  1631.     }
  1632.     if (hasConfidence) {
  1633.       plyVert.confidence = confidence[i];
  1634.     }
  1635.     ply.put_element((void *) &plyVert);
  1636.   }
  1638.   if (strips) {
  1639.     ply.put_element_setup ("tristrips");
  1641.     tstrip_info ts_info = { tris.size(), tris.begin() };
  1642.     ply.put_element (&ts_info);
  1644.   } else {
  1646.     ply.put_element_setup ("face");
  1647.     
  1648.     PlyFace plyFace;
  1649.     for (i = 0; i < tris.size(); i+=3) {
  1650.       plyFace.nverts = 3;
  1651.       plyFace.verts[0] = tris[i+0];
  1652.       plyFace.verts[1] = tris[i+1];
  1653.       plyFace.verts[2] = tris[i+2];
  1654.       
  1655.       ply.put_element_static_strg((void *) &plyFace);
  1656.     }
  1657.   }
  1658. }
  1660. // exported wrappers for write_ply_file
  1662. void
  1663. write_ply_file(const char *fname,
  1664.        const vector<Pnt3> &vtx,
  1665.        const vector<int> &tris, bool strips,
  1666.        const vector<Pnt3> &nrm,
  1667.        const vector<uchar> &intensity,
  1668.        const vector<float> &confidence)
  1669. {
  1670.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence, 
  1671.  true, true, true);
  1672. }
  1674. void
  1675. write_ply_file(const char *fname,
  1676.        const vector<Pnt3> &vtx,
  1677.        const vector<int> &tris, bool strips,
  1678.        const vector<Pnt3> &nrm,
  1679.        const vector<uchar> &intensity)
  1680. {
  1681.   vector<float> confidence;
  1682.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence, 
  1683.  true, true, false);
  1684. }
  1686. void
  1687. write_ply_file(const char *fname,
  1688.        const vector<Pnt3> &vtx,
  1689.        const vector<int> &tris, bool strips,
  1690.        const vector<Pnt3> &nrm)
  1691. {
  1692.   vector<uchar> intensity;
  1693.   vector<float> confidence;
  1694.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence, 
  1695.  true, false, false);
  1696. }
  1698. void
  1699. write_ply_file(const char *fname,
  1700.        const vector<Pnt3> &vtx,
  1701.        const vector<int> &tris, bool strips)
  1702. {
  1703.   vector<Pnt3> nrm;
  1704.   vector<uchar> intensity;
  1705.   vector<float> confidence;
  1706.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence,
  1707.  false, false, false);
  1708. }
  1710. void
  1711. write_ply_file(const char *fname,
  1712.        const vector<Pnt3> &vtx,
  1713.        const vector<int> &tris, bool strips,
  1714.        const vector<uchar> &intensity)
  1715. {
  1716.   vector<Pnt3> nrm;
  1717.   vector<float> confidence;
  1718.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence,
  1719.  false, true, false);
  1720. }
  1722. void
  1723. write_ply_file(const char *fname,
  1724.        const vector<Pnt3> &vtx,
  1725.        const vector<int> &tris, bool strips,
  1726.        const vector<float> &confidence)
  1727. {
  1728.   vector<Pnt3> nrm;
  1729.   vector<uchar> intensity;
  1730.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence,
  1731.  false, false, true);
  1732. }
  1734. void
  1735. write_ply_file(const char *fname,
  1736.        const vector<Pnt3> &vtx,
  1737.        const vector<int> &tris, bool strips,
  1738.        const vector<uchar> &intensity,
  1739.        const vector<float> &confidence)
  1740. {
  1741.   vector<Pnt3> nrm;
  1742.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence,
  1743.  false, true, true);
  1744. }
  1746. void
  1747. write_ply_file(const char *fname,
  1748.        const vector<Pnt3> &vtx,
  1749.        const vector<int> &tris, bool strips,
  1750.        const vector<Pnt3> &nrm,
  1751.        const vector<float> &confidence)
  1752. {
  1753.   vector<uchar> intensity;
  1754.   write_ply_file(fname, vtx, tris, strips, nrm, intensity, confidence, 
  1755.  true, false, true);
  1756. }
  1759. void
  1760. pushNormalAsShorts (vector<short>& nrms, Pnt3 n)
  1761. {
  1762.   n *= 32767;
  1763.   nrms.push_back (n[0]);
  1764.   nrms.push_back (n[1]);
  1765.   nrms.push_back (n[2]);
  1766. }
  1769. void
  1770. pushNormalAsPnt3 (vector<Pnt3>& nrms, short* n, int i)
  1771. {
  1772.   i *= 3;
  1773.   Pnt3 nrm (n[i]/32767.0, n[i+1]/32767.0, n[i+2]/32767.0);
  1774.   nrms.push_back (nrm);
  1775. }