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

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

3D图形编程

开发平台：

Visual C++

VolCarve.cc：源码内容
  1. //############################################################
  2. //
  3. // VolCarve.cc
  4. //
  5. // Matt Ginzton, Kari Pulli
  6. // Mon Jul 13 11:53:22 PDT 1998
  7. //
  8. // Perform a volumetric space carving on RigidScans.
  9. // Use octrees.
  10. // For occlusion testing, treat octree nodes as spheres.
  11. //
  12. //############################################################
  14. #include <map.h>
  15. #include <stack.h>
  16. #include <deque.h>
  17. #ifdef WIN32
  18. # include <float.h>
  19. # define isnanf _isnan
  20. #endif
  21. #ifdef sgi
  22. # include <ieeefp.h>
  23. #endif
  24. #include "VolCarve.h"
  25. #include "MCEdgeTable.h"
  28. #define VOLCARVE_MAIN 0
  30. // face indices:
  31. // 0, 2, 4 = dimension (x, y, z); +0 negative, +1 positive
  32. // so 0 is left, 1 is right
  33. //    2 is down, 3 is up
  34. //    4 is in,   5 is out
  35. static const int face_negX = 0;
  36. //static const int face_posX = 1;
  37. static const int face_negY = 2;
  38. //static const int face_posY = 3;
  39. static const int face_negZ = 4;
  40. //static const int face_posZ = 5;
  44. vector<Pnt3> vertex_list;
  46. /*
  47. #ifdef WIN32
  48.   typedef __int64 CUBEKEY;
  49. #else
  50.   typedef long long CUBEKEY;
  51. #endif
  52. typedef map<CUBEKEY, int> VERTEXMAP;
  53. VERTEXMAP g_vertexMap;
  55. Pnt3 fblo, fbhi; // full box lower and upper corners
  56. // use the two upper two bits to encode which edge of a cube corner
  57. // the vertex belongs to
  58. static const int encode_bits  = ((8 * sizeof(CUBEKEY))-2) / 3;
  59. static const int encode_shift = (1 << encode_bits);
  60. static float encode_scale;
  63. static inline CUBEKEY
  64. encode(const Pnt3 &p)
  65. {
  66.   CUBEKEY x = (p[0]-fblo[0]) * encode_scale;
  67.   CUBEKEY y = (p[1]-fblo[1]) * encode_scale;
  68.   CUBEKEY z = (p[2]-fblo[2]) * encode_scale;
  69.   return (((x << encode_bits) + y) << encode_bits) + z;
  70. }
  71. */
  74. class CTL_memhandler {
  75. private:
  76.   int blockSize;
  77.   stack<CubeTreeLeafData*> free;
  78.   int count;
  80. public:
  82.   CTL_memhandler(void) : blockSize(0) {}
  83.   ~CTL_memhandler(void) 
  84.     {
  85.       //SHOW(free.size());
  86.       //SHOW(count);
  87.       while (!free.empty()) {
  88. delete[] free.top();
  89. free.pop();
  90.       }
  91.     }
  93.   void set_blocksize(int bs) 
  94.     {
  95.       if (bs != blockSize) {
  96. while (!free.empty()) {
  97.   delete[] free.top();
  98.   free.pop();
  99. }
  100. blockSize = bs;
  101. count     = 0;
  102.       }
  103.     }
  105.   CubeTreeLeafData *operator()(void)
  106.     {
  107.       CubeTreeLeafData *ret;
  108.       if (!free.empty()) { 
  109. ret = free.top(); 
  110. free.pop(); 
  111. return ret; 
  112.       } else {
  113. count++;
  114. return new CubeTreeLeafData[blockSize];
  115.       }
  116.     }
  118.   void store(CubeTreeLeafData *p)
  119.     {
  120.       free.push(p);
  121.     }
  122. };
  124. CTL_memhandler mem_handler;
  127. // first, implementation of methods used by all CubeTrees (CubeTreeBase),
  128. // followed by specifics for octree nodes referenced by pointers
  129. // (CubeTree, which doesn't need to have a complete set of children)
  130. // and octree nodes that store all their children in a flat array
  131. // (CubeTreeLeaf, which always has a complete set of children).
  134. int   CubeTreeBase::leafDepth = CubeTree::defaultLeafDepth;
  135. float CubeTreeBase::leafSize  = 0.0;
  136. int   CubeTreeBase::nVoxels   = 0;
  138. CubeTreeBase::CubeTreeBase (Pnt3 c, float s, CubeTree* _par, int _iParIdx)
  139.   : mCtr(c), mSize(s), parent(_par), parIdx(_iParIdx), type(INDETERMINATE)
  140. { }
  143. CubeTreeBase::~CubeTreeBase()
  144. { }
  146. void
  147. CubeTreeBase::set_leaf_depth (int ld)
  148. {
  149.   if (ld >= 0)  leafDepth = ld;
  150.   else          leafDepth = defaultLeafDepth;
  151. }
  154. CubeTreeBase*
  155. CubeTreeBase::find_neighbor (int face)
  156. {
  157.   if (!parent)
  158.     return NULL;
  160.   int iChild;
  161.   bool sibling = is_neighbor_cell(face, iChild);
  162.   if (sibling) {
  163.     CubeTreeBase* theSibling = parent->child[iChild];
  164.     if (!theSibling) return NULL;
  165.     return theSibling;
  166.   }
  168.   // not a sibling: push child and go up
  169.   stack<int> path;
  170.   path.push (iChild);
  171.   return parent->recursive_find (face, path, false);
  172. }
  177. bool
  178. CubeTreeBase::is_neighbor_cell (int face, int& iChild)
  179. {
  180.   // given a direction, return the child index of the neighbor
  181.   // that lies in that direction.  If that neighbor is a sibling
  182.   // (iChild is valid for the same parent as us), return true.
  183.   // Otherwise, return false.
  185.   int axis = face / 2;
  186.   int sign = face % 2;
  188.   // neighbor will be in same position as we are, except
  189.   // the bit representing the axis we're moving along will be 
  190.   // flipped
  191.   iChild = parIdx ^ (1 << axis);
  193.   // if we have the bit for that axis set and the direction of
  194.   // movement is negative, it's a sibling, but if the direction 
  195.   // of movement is positive, it's not a sibling.  
  196.   // Vice-versa if our child index does
  197.   // not have the bit for the movement axis set.
  198.   int haveSign = (0 < (parIdx & (1 << axis)));
  199.   return (sign != haveSign);
  200. }
  203. //
  204. // Nonmember function to determine the type of a node with respect to
  205. // the views, following the precedence order:
  206. //   OUTSIDE
  207. //   BOUNDARY
  208. //   SILHOUETTE
  209. //   INSIDE
  210. //   NOT_IN_FRUSTUM
  211. //   INDETERMINATE
  212. //
  213. OccSt
  214. view_query_status (vector<RigidScan*> &views,
  215.    vector<RigidScan*> &remaining_views,
  216.    const Pnt3 &ctr, float size,
  217.    OccSt type = INDETERMINATE)
  218. {
  219.   vector<RigidScan*>::iterator vit; 
  220.   for (vit = views.begin(); vit != views.end(); vit++) {
  221.     OccSt t = (*vit)->carve_cube(ctr, size);
  222.     // OccSt types listed in descending precedence order
  223.     if (int(type) > int(t))
  224.       type = t;
  226.     if (t != NOT_IN_FRUSTUM) {
  227.       remaining_views.push_back(*vit);
  228.     }
  230.     if (t == OUTSIDE) {
  231.       // can't get any better than this
  232.       return t;
  233.     }
  234.   }
  236.   return type;
  237. }
  240. float
  241. view_query_distance (vector<RigidScan*> &views, const Pnt3& ctr,
  242.      float& distance)
  243. {
  244.   float totalConf = 0;
  245.   OccSt occ;
  246.   Pnt3 closest;
  247.   distance = 0;
  249.   //
  250.   // HOW SHOULD THE DISTANCES REALLY BE COMBINED????
  251.   //
  253.   float thr  = 15.0;
  254.   float ramp = 30.0;
  256.   vector<RigidScan*>::iterator vit; 
  257.   for (vit = views.begin(); vit != views.end(); vit++) {
  258.     float conf = (*vit)->closest_along_line_of_sight (ctr, closest, occ);
  259.     //float conf = (*vit)->closest_point_on_mesh(ctr, closest, occ);
  261.     float thisDist = dist(ctr, closest);
  263.     if (occ == INSIDE) {
  264.       if (thisDist > thr) {
  265. float d = thisDist - thr;
  266. if (d > ramp) conf = 0;
  267. else          conf *= (ramp - d) / ramp;
  268.       }
  269.     }
  270.     totalConf += conf;
  271.     if (occ == INSIDE)
  272.       thisDist *= -1;
  273.     distance += conf * thisDist;
  274.   }
  276.   if (totalConf) distance /= totalConf;
  277.   else           distance = -1000;
  279.   return totalConf;
  280. }
  285. ///////////////////////////////////////////////////////////////////////////
  286. //
  287. // Pointerized octree nodes: CubeTree
  288. //
  289. ///////////////////////////////////////////////////////////////////////////
  291. void 
  292. CubeTree::delete_children (void)
  293. { 
  294.   for (int i = 0; i < 8; i++) {
  295.     delete child[i]; 
  296.     child[i] = NULL;
  297.   }
  298. }
  301. CubeTreeBase*
  302. CubeTree::recursive_find(int face, stack<int>& path, 
  303.  bool bReachedTop)
  304. {
  305.   if (path.empty()) {
  306.     /*
  307.     CubeTreeLeaf* leaf = dynamic_cast<CubeTreeLeaf*> (this);
  308.     if (!leaf) {
  309.       assert (0);
  310.       return this;
  311.     }
  312.     return leaf;
  313.     */
  314.     return this;
  315.   }
  317.   int iChild;
  318.   if (bReachedTop) {
  319.     // then we want to descend
  320.     iChild = path.top(); path.pop();
  321.     CubeTreeBase* theChild = child[iChild];
  322.     if (!theChild)
  323.       return this; // child does not exist
  324.     return theChild->recursive_find (face, path, true);
  325.   } else {
  326.     // we still need to traverse horizontally...
  327.     if (!parent)
  328.       return NULL;
  329.     bool sibling = is_neighbor_cell (face, iChild);
  330.     if (sibling) {
  331.       // if we have a sibling, we can traverse 
  332.       // (and we've reached the top!)
  333.       CubeTreeBase* sibling = parent->child[iChild];
  334.       if (!sibling) return NULL;
  335.       return sibling->recursive_find (face, path, true);
  336.     } else {
  337.       // otherwise we have to go farther up.
  338.       path.push (iChild);
  339.       return parent->recursive_find (face, path, false);
  340.     }
  341.   }
  342. }
  344. // find the neighbors of child i into negative directions
  345. // in x, y, and z, and tell them that data along that
  346. // face is not needed any more
  347. void
  348. CubeTree::release_space(int i)
  349. {
  350.   if (i&1) {
  351.     child[i-1]->keep_walls(0,1,1);
  352.   } else {
  353.     CubeTreeBase *ctb = find_neighbor(face_negX);
  354.     if (ctb) {
  355.       ctb = (dynamic_cast<CubeTree*>(ctb))->child[i+1];
  356.       if (ctb) ctb->keep_walls(0,1,1);
  357.     }
  358.   }
  359.   if (i&2) {
  360.     child[i-2]->keep_walls(1,0,1);
  361.   } else {
  362.     CubeTreeBase *ctb = find_neighbor(face_negY);
  363.     if (ctb) {
  364.       ctb = (dynamic_cast<CubeTree*>(ctb))->child[i+2];
  365.       if (ctb) ctb->keep_walls(1,0,1);
  366.     }
  367.   }
  368.   if (i&4) {
  369.     child[i-4]->keep_walls(1,1,0);
  370.   } else {
  371.     CubeTreeBase *ctb = find_neighbor(face_negZ);
  372.     if (ctb) {
  373.       ctb = (dynamic_cast<CubeTree*>(ctb))->child[i+4];
  374.       if (ctb) ctb->keep_walls(1,1,0);
  375.     }
  376.   }
  377. }
  380. CubeTree::CubeTree (Pnt3 c, float s, CubeTree* _par, int _i)
  381.   : CubeTreeBase (c, s, _par, _i)
  382. {
  383.   for (int i = 0; i < 8; i++)
  384.     child[i] = NULL; 
  385. }
  388. CubeTree::~CubeTree (void)
  389. {
  390.   delete_children();
  391. }
  394. // recursively descend up to given level
  395. // then if the neighboring cube is outside,
  396. // create a face
  397. void
  398. CubeTree::extract_faces(vector<int> &tri_inds)
  399. {
  400.   // treat NOT_IN_FRUSTUM as OUTSIDE for now
  401.   // for INSIDE should perhaps test that it is truly
  402.   // all the way inside
  403.   //if (!boundary())
  404.   if (!go_on())
  405.     return;
  407.   // descend
  408.   for (int i = 0; i < 8; i++) {
  409.     assert(child[i]);
  410.     child[i]->extract_faces(tri_inds);
  411.   }
  412. }
  415. void
  416. CubeTree::keep_walls(bool x, bool y, bool z)
  417. {
  418.   int l=0;
  419.   for (int i=0; i<2; i++) {
  420.     for (int j=0; j<2; j++) {
  421.       for (int k=0; k<2; k++, l++) {
  422. if (child[l] == NULL) continue;
  423. if (k&&x || j&&y || i&&z) {
  424.   child[l]->keep_walls(x&&k, y&&j, z&&i);
  425. } else {
  426.   delete child[l];
  427.   child[l] = NULL;
  428. }
  429.       }
  430.     }
  431.   }    
  432. }
  435. void 
  436. CubeTree::carve_help(vector<RigidScan*> &views, int levels,
  437.      vector<int> &tri_inds,
  438.      float perc_min, float perc_max)
  439. {
  440.   // When checking status inflate size by the size of a leaf
  441.   // so if the surface goes close to the  CubeTreeLeaf block boundary, 
  442.   // its neighbors exist in the tree. 
  443.   vector<RigidScan*> remaining_views;
  444.   type = ::view_query_status (views, remaining_views,
  445.       mCtr, mSize+leafSize, type);
  447.   SHOW(levels);
  449.   //// only boundary nodes can have children
  450.   //if (!boundary()) return;
  451.   if (!go_on()) return;
  452.   if (levels == 0) return;
  453.   --levels;
  455.   float childSize = mSize*.5;
  456.   float stepSize  = mSize*.25;
  458.   bool allChildrenOut = true;
  459.   bool allChildrenIn  = true;
  461.   Pnt3 childCtr;
  463.   float perc_frac = (perc_max-perc_min) / 8.0;
  465.   for (int i=0; i<8; i++) {
  466.     childCtr.set(mCtr[0] + ((i&1) ? stepSize : -stepSize),
  467.  mCtr[1] + ((i&2) ? stepSize : -stepSize),
  468.  mCtr[2] + ((i&4) ? stepSize : -stepSize));
  470.     if (levels > leafDepth) {
  471.       CubeTree* ct = new CubeTree (childCtr, childSize, this, i);
  472.       child[i] = ct;
  473.       float _perc_min = perc_min + i*perc_frac;
  474.       ct->carve_help(remaining_views, levels, tri_inds,
  475.      _perc_min, _perc_min + perc_frac);
  476.       cout << "Carved " << _perc_min + perc_frac << " percent" << endl;
  477.     } else {
  478.       CubeTreeLeaf* ctl = new CubeTreeLeaf(childCtr, childSize,
  479.    this, i, views);
  480.       child[i] = ctl;
  481.       ctl->extract_faces(tri_inds);
  482.     }
  484.     /*    
  485.     // release unneeded parts of the tree
  486.     if (i&1) child[i-1]->keep_walls(0,1,1);
  487.     else {
  488.       CubeTreeBase *ctb = find_neighbor(face_negX);
  489.       if (ctb) {
  490. CubeTree *ct = dynamic_cast<CubeTree*>(ctb);
  491. if (ct->child[i+1]) ct->child[i+1]->keep_walls(0,1,1);
  492.       }
  493.     }
  494.     if (i&2) child[i-2]->keep_walls(1,0,1);
  495.     if (i&4) child[i-4]->keep_walls(1,1,0);
  496.     */
  497.     release_space(i);
  499.     OccSt t = child[i]->type;
  500.     if (t != OUTSIDE) allChildrenOut = false;
  501.     if (t != INSIDE)  allChildrenIn  = false;
  502.   }
  504.   if (allChildrenOut) {
  505.     delete_children();
  506.     type = OUTSIDE;
  507.   } 
  508.   if (allChildrenIn) {
  509.     delete_children();
  510.     type = INSIDE;
  511.   }
  512. }
  514. //#define VIS_TRIOCTREE
  516. #ifdef VIS_TRIOCTREE
  517. #include "CyraScan.h"
  518. #include "TriOctree.h"
  519. #endif
  521. // first split, then
  522. // extract the boundary between outsides and the rest
  523. void 
  524. CubeTree::carve(vector<RigidScan*> &views, int levels,
  525. vector<Pnt3> &coords, vector<int> &tri_inds)
  526. {
  527. #ifdef VIS_TRIOCTREE
  528.   CyraScan *cs = dynamic_cast<CyraScan*> (views[0]);
  529.   if (cs) cs->TriOctreeMesh(coords, tri_inds);
  530.   return;
  531. #endif
  533.   // Figure out the size of the leaves.
  534.   leafSize = mSize; 
  535.   for (int j = 0; j<levels; j++) leafSize *= .5;
  536.   nVoxels  = 1<<3*leafDepth;
  537.   mem_handler.set_blocksize(nVoxels);
  539.   tri_inds.clear();
  540.   vertex_list.clear();
  541.   /*
  542.   fbhi = fblo = mCtr;
  543.   Pnt3 step(mSize*.5,mSize*.5,mSize*.5);
  544.   fbhi += step;
  545.   fblo -= step;
  546.   encode_scale = .99999 * encode_shift / (fbhi[0] - fblo[0]);
  547.   */
  549.   carve_help(views, levels, tri_inds, 0, 100);
  551.   delete_children();
  553.   coords = vertex_list;
  554.   SHOW(vertex_list.size());
  555.   SHOW(tri_inds.size());
  556.   vertex_list.clear();
  558.   for (int i=0; i<coords.size(); i++) {
  559.     if (isnanf(coords[i][0]) ||
  560. isnanf(coords[i][1]) ||
  561. isnanf(coords[i][2])) {
  562.       SHOW(i); SHOW(coords[i]);
  563.       break;
  564.     }
  565.   }
  566. }
  568. //#define CARVE_TEST
  570. static int leaf_count = 0;
  572. CubeTreeLeaf::CubeTreeLeaf (const Pnt3 &c, float s,
  573.     CubeTree* _par, int _i,
  574.     vector<RigidScan*> &views)
  575.   : CubeTreeBase (c, s, _par, _i), data(NULL)
  576. {
  577.   // When checking status inflate size by the size of a leaf
  578.   // so if the surface goes close to the  CubeTreeLeaf block boundary, 
  579.   // it's neighbors exist in the tree. 
  580.   vector<RigidScan*> remaining_views;
  581.   type = ::view_query_status (views, remaining_views,
  582.       c, s+leafSize);
  583.   //type = ::view_query_status (views, c, 1.25*leafSize);
  585.   //if (!boundary())
  586.   if (!go_on())
  587.     return;
  589.   leaf_count++;
  590.   SHOW(leaf_count);
  592. #ifdef CARVE_TEST
  593.   return;
  594. #endif
  596.   nLeavesPerSide = 1 << leafDepth;
  597.   mask = nLeavesPerSide - 1;
  599.   //data = new CubeTreeLeafData[nVoxels];
  600.   data = mem_handler();
  602.   float tmp = .5 * (nLeavesPerSide - 1) * leafSize;
  603.   Pnt3 base(mCtr[0]-tmp, mCtr[1]-tmp, mCtr[2]-tmp);
  605.   Pnt3 lCtr;
  606.   int depth2 = 2*leafDepth;
  607.   int x, y, z;
  609.   // evaluate the function at all the leaf nodes
  610.   // also interpolate the edges here
  611.   int ofsx = 1;
  612.   int ofsy = 1<<leafDepth;
  613.   int ofsz = 1<<depth2;
  614.   // length of the step to go to the other end along this dimension
  615.   int lenx = ofsy-ofsx;
  616.   int leny = ofsz-ofsy;
  617.   int lenz = (1<<leafDepth+depth2)-ofsz;
  618.   CubeTreeLeafData *vp;
  620.   CubeTreeLeaf* nborX = find_neighbor (face_negX);
  621.   if (nborX && !nborX->data) nborX = NULL;
  622.   CubeTreeLeaf* nborY = find_neighbor (face_negY);
  623.   if (nborY && !nborY->data) nborY = NULL;
  624.   CubeTreeLeaf* nborZ = find_neighbor (face_negZ);
  625.   if (nborZ && !nborZ->data) nborZ = NULL;
  627.   int i = 0;
  628.   for (z=0; z<nLeavesPerSide; z++) {
  629.     for (y=0; y<nLeavesPerSide; y++) {
  630.       for (x=0; x<nLeavesPerSide; x++, i++) {
  632. float &dist = data[i].distance;
  634. lCtr.set(base[0] + leafSize*x, 
  635.  base[1] + leafSize*y, 
  636.  base[2] + leafSize*z);
  637. data[i].confidence = view_query_distance(views, lCtr,dist);
  639. // if zero confidence, assign arbitrarily dist = 1.0
  640. if (data[i].confidence < 1e-5) dist = 1.0;
  642. assert(dist > -1e22 && dist < 1e22);
  643. // interpolate edges
  644. if (x) {
  645.   vp = &data[i-ofsx];
  646. } else {
  647.   if (nborX)  vp = &nborX->data[i+lenx];
  648.   else        vp = NULL;
  649. }
  650. if (vp && dist * vp->distance <= 0.0) {
  651.   // signs differ, interpolate
  652.   vp->eixi = vertex_list.size();
  653.   float d = leafSize * dist / (dist - vp->distance);
  654.   assert(!isnanf(d));
  655.   vertex_list.push_back(Pnt3(lCtr[0]-d, lCtr[1], lCtr[2]));
  656. }
  658. if (y) {
  659.   vp = &data[i-ofsy];
  660. } else {
  661.   if (nborY)  vp = &nborY->data[i+leny];
  662.   else        vp = NULL;
  663. }
  664. if (vp && dist * vp->distance <= 0.0) {
  665.   // signs differ, interpolate
  666.   vp->eiyi = vertex_list.size();
  667.   float d = leafSize * dist / (dist - vp->distance);
  668.   assert(!isnanf(d));
  669.   vertex_list.push_back(Pnt3(lCtr[0], lCtr[1]-d, lCtr[2]));
  670. }
  672. if (z) {
  673.   vp = &data[i-ofsz];
  674. } else {
  675.   if (nborZ)  vp = &nborZ->data[i+lenz];
  676.   else        vp = NULL;
  677. }
  678. if (vp && dist * vp->distance <= 0.0) {
  679.   // signs differ, interpolate
  680.   vp->eizi = vertex_list.size();
  681.   float d = leafSize * dist / (dist - vp->distance);
  682.   assert(!isnanf(d));
  683.   vertex_list.push_back(Pnt3(lCtr[0], lCtr[1], lCtr[2]-d));
  684. }
  685.       }
  686.     }
  687.   }
  688. }
  691. CubeTreeLeaf::~CubeTreeLeaf (void)
  692. {
  693.   //delete[] data;
  694.   if (data) mem_handler.store(data);
  695. }
  698. static void
  699. cube_bdry(const Pnt3 &ctr, float size, 
  700.   vector<Pnt3> &p, vector<int> &i)
  701. {
  702.   float w = .5 * size;
  703.   p.push_back(ctr + Pnt3(-w,-w,-w));
  704.   p.push_back(ctr + Pnt3(-w,-w, w));
  705.   p.push_back(ctr + Pnt3(-w, w,-w));
  706.   p.push_back(ctr + Pnt3(-w, w, w));
  707.   p.push_back(ctr + Pnt3( w,-w,-w));
  708.   p.push_back(ctr + Pnt3( w,-w, w));
  709.   p.push_back(ctr + Pnt3( w, w,-w));
  710.   p.push_back(ctr + Pnt3( w, w, w));
  712.   static int triinds[] = {
  713.     0,1,3,2,0,3, 0,4,5,0,5,1, 2,6,4,2,4,0,
  714.     3,5,7,3,1,5, 2,7,6,2,3,7, 7,4,6,7,5,4
  715.   };
  717.   int   s = vertex_list.size();
  718.   for (int j=0; j<36; j++) 
  719.     i.push_back(s+triinds[j]);
  720. }
  722. static EdgeTable TheEdgeTable;
  724. void
  725. CubeTreeLeaf::extract_faces(vector<int> &tri_inds)
  726. {
  727.   //if (!boundary()) return;
  728.   if (!go_on()) return;
  730. #ifdef CARVE_TEST
  731.   cube_bdry(mCtr, mSize, vertex_list, tri_inds);
  732.   return;
  733. #endif
  735.   int depth2 = 2*leafDepth;
  736.   int x, y, z;
  738.   int ofsx = 1;
  739.   int ofsy = 1<<leafDepth;
  740.   int ofsz = 1<<depth2;
  741.   // length of the step to go to the other end along this dimension
  742.   int lenx = ofsy-ofsx;
  743.   int leny = ofsz-ofsy;
  744.   int lenz = nVoxels-ofsz;
  745.   for (int i = 0; i < nVoxels; i++) {
  746.     // set CubeTreeLeafData static center and size
  747.     // nVoxels has depth bits each for x, y, z
  748.     // find center of corresponding voxel i
  749.     x = (i)              & mask;
  750.     y = (i >> leafDepth) & mask;
  751.     z = (i >> depth2)    & mask;
  753.     //if (x&&y&&z) continue;
  755.     // get the eight pointers to vertices
  756.     // create also the index to the edge table
  757.     CubeTreeLeafData *vp[8];
  758.     int edgeTableIndex = 0;
  759.     CubeTreeLeaf *nborX = NULL, *nborY = NULL, *nborZ = NULL;
  760.     CubeTreeLeaf *nborXY = NULL, *nborXZ = NULL, *nborYZ = NULL;
  762.     // these are the default vertices
  763.     // the ones with a negative index to data will be replaced
  764.     int j;
  765.     vp[0] = ((j=i-ofsx-ofsy-ofsz) < 0 ? NULL : &data[j]);
  766.     vp[1] = ((j=i     -ofsy-ofsz) < 0 ? NULL : &data[j]);
  767.     vp[2] = ((j=i          -ofsz) < 0 ? NULL : &data[j]);
  768.     vp[3] = ((j=i-ofsx     -ofsz) < 0 ? NULL : &data[j]);
  769.     vp[4] = ((j=i-ofsx-ofsy     ) < 0 ? NULL : &data[j]);
  770.     vp[5] = ((j=i     -ofsy     ) < 0 ? NULL : &data[j]);
  771.     vp[6] = &data[i];
  772.     vp[7] = ((j=i-ofsx          ) < 0 ? NULL : &data[j]);
  774.     if (x == 0) {
  775.       nborX = find_neighbor (face_negX);
  776.       if (nborX && nborX->data) {
  777. vp[0] = ((j=i+lenx-ofsy-ofsz) < 0 ? NULL : &nborX->data[j]);
  778. vp[3] = ((j=i+lenx     -ofsz) < 0 ? NULL : &nborX->data[j]);
  779. vp[4] = ((j=i+lenx-ofsy     ) < 0 ? NULL : &nborX->data[j]);
  780. vp[7] = &nborX->data[i+lenx];
  781.       } else {
  782. continue; // no neighbor, skip
  783.       }
  784.     }
  786.     if (y == 0) {
  787.       nborY = find_neighbor (face_negY);
  788.       if (nborY && nborY->data) {
  789. vp[0] = ((j=i+leny-ofsx-ofsz) < 0 ? NULL : &nborY->data[j]);
  790. vp[1] = ((j=i+leny     -ofsz) < 0 ? NULL : &nborY->data[j]);
  791. vp[4] = &nborY->data[i+leny-ofsx];
  792. vp[5] = &nborY->data[i+leny];
  793.       } else {
  794. continue; // no neighbor, skip
  795.       }
  796.     }
  797.     
  798.     if (z == 0) {
  799.       nborZ = find_neighbor (face_negZ);
  800.       if (nborZ && nborZ->data) {
  801. vp[0] = &nborZ->data[i+lenz-ofsx-ofsy];
  802. vp[1] = &nborZ->data[i+lenz     -ofsy];
  803. vp[2] = &nborZ->data[i+lenz          ];
  804. vp[3] = &nborZ->data[i+lenz-ofsx     ];
  805.       } else {
  806. continue; // no neighbor, skip
  807.       }
  808.     }
  810.     if (nborX && nborY) { // then out in both x and y
  811.       nborXY = nborY->find_neighbor (face_negX);
  812.       if (nborXY && nborXY->data) {
  813. vp[0] = ((j=i+lenx+leny-ofsz) < 0 ? NULL : &nborXY->data[j]);
  814. vp[4] = &nborXY->data[i+lenx+leny];
  815.       } else {
  816. continue; // no neighbor, skip
  817.       }
  818.     }
  820.     if (nborX && nborZ) {
  821.       nborXZ = nborX->find_neighbor (face_negZ);
  822.       if (nborXZ && nborXZ->data) {
  823. vp[0] = &nborXZ->data[i+lenx+lenz-ofsy];
  824. vp[3] = &nborXZ->data[i+lenx+lenz];
  825.       } else {
  826. continue; // no neighbor, skip
  827.       }
  828.     }
  830.     if (nborY && nborZ) {
  831.       nborYZ = nborZ->find_neighbor (face_negY);
  832.       if (nborYZ && nborYZ->data) {
  833. vp[0] = &nborYZ->data[i+leny+lenz-ofsx];
  834. vp[1] = &nborYZ->data[i+leny+lenz];
  835.       } else {
  836. continue; // no neighbor, skip
  837.       }
  838.     }
  839.     
  840.     if (nborXY && nborZ) {
  841.       CubeTreeLeaf* nborXYZ = nborXY->find_neighbor (face_negZ);
  842.       if (nborXYZ && nborXYZ->data) {
  843. vp[0] = &nborXYZ->data[i+lenx+leny+lenz];
  844.       } else {
  845. continue; // no neighbor, skip
  846.       }
  847.     }
  849.     for (j=0; j<8; j++) 
  850.       assert(vp[j]->distance > -1e22 && vp[j]->distance < 1e22);
  851.     
  852.     if (vp[0]->distance > 0.0) edgeTableIndex |=   1;
  853.     if (vp[1]->distance > 0.0) edgeTableIndex |=   2;
  854.     if (vp[2]->distance > 0.0) edgeTableIndex |=   4;
  855.     if (vp[3]->distance > 0.0) edgeTableIndex |=   8;
  856.     if (vp[4]->distance > 0.0) edgeTableIndex |=  16;
  857.     if (vp[5]->distance > 0.0) edgeTableIndex |=  32;
  858.     if (vp[6]->distance > 0.0) edgeTableIndex |=  64;
  859.     if (vp[7]->distance > 0.0) edgeTableIndex |= 128;
  861.     // for all the 12 edges, if the table tells that there should
  862.     // be an edge, copy a pointer to the corresponding precalculated vertex
  863.     bool *edge = TheEdgeTable[edgeTableIndex].edge;
  864.     int evtx[12];
  865.     if (edge[0])  evtx[0]  = vp[0]->eixi;
  866.     if (edge[1])  evtx[1]  = vp[1]->eiyi;
  867.     if (edge[2])  evtx[2]  = vp[3]->eixi;
  868.     if (edge[3])  evtx[3]  = vp[0]->eiyi;
  869.     if (edge[4])  evtx[4]  = vp[4]->eixi;
  870.     if (edge[5])  evtx[5]  = vp[5]->eiyi;
  871.     if (edge[6])  evtx[6]  = vp[7]->eixi;
  872.     if (edge[7])  evtx[7]  = vp[4]->eiyi;
  873.     if (edge[8])  evtx[8]  = vp[0]->eizi;
  874.     if (edge[9])  evtx[9]  = vp[1]->eizi;
  875.     if (edge[10]) evtx[10] = vp[3]->eizi;
  876.     if (edge[11]) evtx[11] = vp[2]->eizi;
  878.     int nTrisCreated = TheEdgeTable[edgeTableIndex].Ntriangles;
  879.     Triple *triList  = TheEdgeTable[edgeTableIndex].TriangleList;
  880.     for (j=0; j<nTrisCreated; j++) {
  881. #if 0
  882.       SHOW (evtx[triList[j].A]);
  883.       SHOW (evtx[triList[j].B]);
  884.       SHOW (evtx[triList[j].C]);
  885. #endif
  886.       assert(evtx[triList[j].A]>-1);
  887.       assert(evtx[triList[j].B]>-1);
  888.       assert(evtx[triList[j].C]>-1);
  889.       tri_inds.push_back(evtx[triList[j].A]);
  890.       tri_inds.push_back(evtx[triList[j].B]);
  891.       tri_inds.push_back(evtx[triList[j].C]);
  892.     }
  893.   }
  894. }
  897. void
  898. CubeTreeLeaf::keep_walls(bool , bool , bool )
  899. {
  900. }
  903. OccSt
  904. CubeTreeLeaf::neighbor_status (int face, int child)
  905. {
  906.   // get address, and child status, of neighbor of given child
  907.   bool sibling = is_neighbor_cell (face, child);
  908.   if (sibling)
  909.     return type;
  911.   // not a child of this block; have to go up tree
  912.   CubeTreeBase* ctb = CubeTreeBase::find_neighbor (face);
  913.   CubeTreeLeaf* ctl = dynamic_cast<CubeTreeLeaf*> (ctb);
  914.   if (!ctb)
  915.     return OUTSIDE;
  917.   return ctl->type;
  918. }
  921. CubeTreeBase*
  922. CubeTreeLeaf::recursive_find(int face, 
  923.      stack<int>& path, bool bReachedTop)
  924. {
  925.   return this;
  926. }
  929. CubeTreeLeaf*
  930. CubeTreeLeaf::find_neighbor (int face)
  931. {
  932.   return dynamic_cast<CubeTreeLeaf*> (CubeTreeBase::find_neighbor (face));
  933. }
  936. bool
  937. CubeTreeLeaf::is_neighbor_cell (int face, int& child)
  938. {
  939.   // child is both input and output parameter...
  940.   bool contained = true;
  942.   int axis = face / 2;
  943.   int sign = face % 2;
  945.   // need to extract relevant bits, munge them, and replace them into child:
  946.   // extract
  947.   int shift = axis * leafDepth;
  948.   int coord = (child & (mask << shift)) >> shift;
  950.   // munge
  951.   coord += (sign ? +1 : -1);
  952.   if (coord >= nLeavesPerSide) {
  953.     contained = false;
  954.     coord -= nLeavesPerSide;
  955.   }
  956.   else if (coord < 0) {
  957.     contained = false;
  958.     coord += nLeavesPerSide;
  959.   }
  961.   // replace
  962.   child = (child & ~(mask << shift)) | (coord << shift);
  964.   return contained;
  965. }