Biomechanical Joint Model  Author: Anderson Maciel

---

marching_cubes.cpp

Go to the documentation of this file.

/*@#---------------------------------------------------------------------------
 PROJECT            : PHD THESIS
 MODULE             : PHYSICALLY-BASED DEFORMATION

 FILE               : marching_cubes.C
 CREATION DATE      : Thu Nov 11 11:30:03 MET 1999
 CREATION AUTHOR(S) : aubel 
 ------------------------------------------------------------------------------
 KEYWORD(S)         : marching cubes, voxelisation, volume
 DEFINED TYPE(S)    : 
 RELATED TYPE(S)    : 

 FUNCTIONALITY      : 
 ------------------------------------------------------------------------------
 LANGAGE            : Ansi C++
 SYSTEM             : UNIX V 6.x / SGI
 GRAPHIC LIBRARY    : SGI  
 ------------------------------------------------------------------------------
 PROJECT DIRECTOR   : D. THALMANN
 LAB                : EPFL- LIG
 --------------------------------------------------------------------------#@*/
#include "marching_cubes.h"
#include "hash_map.h"
#include "primitive.h"
#include <math.h>
#include <vector>  
#include <list>

#include <time.h>

//namespace MarchingCubes {

// Freely adapted from basic code written by Paul Bourke 
// http://www.swin.edu.au/astronomy/pbourke/modelling/polygonise/

// LL = lower left corner of BBox, UR = upper right corner...
idx3
marchingCubes::findSeedVoxel()
{
   Ray ray;
   float inc_x = (UR[0]-LL[0])/64.0;
   float inc_y = (UR[1]-LL[1])/64.0;
   VEC M = 0.5*(UR+LL);
    
   ray.O[2] = LL[2]; // correct in Z
   ray.D[0] = ray.D[1] = 0; ray.D[2] = UR[2]-LL[2];    // ray along +Z
   ray.segment = true;

   while ( inc_x < (UR[0]-LL[0])/2 )
   {
      // ray-casting starts from center because of increased chances of hitting the surface
      for (float i=-inc_x; i<=inc_x; i += inc_x)
      {
         ray.O[0] = i + M[0];

         if (ray.O[0] <= LL[0] || ray.O[0] >= UR[0]) break;
         for (float j=-inc_y; j<=inc_y; j += inc_y)
         {
            ray.O[1] = j + M[1];

            bool found=false;
            REAL tmin = infinity;
            for (int k=0; k<primlist.size(); k++)
            {          
               if ( primlist[k]->intersectedBy(ray) && ray.t < tmin) 
               {
                  tmin = ray.t;
                  found = true;
               }
            }
            if (found)
            { 
               // construct voxel-cube that encompasses that (intersection) point
               // this'll be our starting cube
               ray.P = ray.O + tmin * ray.D;
               cout << "intersection found for ray " << ray.P << endl; 
               VEC D = (ray.P - LL)/res;
               // printf("Starting voxel index=[%d,%d,%d]\n",(int)D[0], (int)D[1], (int)D[2]); getchar();
               return idx3((int)D[0], (int)D[1], (int)D[2]);   
            }
            else cout << "warning: intersection not found for ray " << ray.O << endl; 
         }
      }
      inc_x *= 2;
      inc_y *= 2;
   }
       
   cout << "\nFatal error: could not find any intersection by ray-casting for marching cubes\n";
   return idx3();
}

/*
static void
printCube(GridCell &cube)
{
   for (int i=0; i<cube.val[0].size(); i++)
      for (int p=0; p<8; p++)
      {
         cout << cube.p[p] << " has density " << cube.val[p][i];
         cout << endl;
      }

}
*/

// given a cube whose vertices have associated density values, create
// a triangulation (at most 5 triangles)
int
marchingCubes::polygonise(GridCell &cube)
{

#include "tables.cpp"

// printCube(cube);

   // for each vertex of cube the corresponding bit in cubeindex
   //  is set to one iff it is outside of ALL primitives
   int cubeindex=0;
   for (int p=0, bit=1; p<8; p++, bit<<=1)
   {
      bool outside = true;
      for (int k=0; k<primlist.size() && outside; k++)
         if (cube.val[p][k] >= 0) outside = false;

      if ( outside ) cubeindex |= bit;
   }

   if (edgeTable[cubeindex] == 0) 
      return faceTable[cubeindex];


   int index[12];
   static VEC V(3), N(3);

   // Find the vertices (and their normals) where the surface intersects the cube
   // Always check first intersection has not been calculated before

   // check 12 edges
   int i;
   for (int bit2=1, i=0; i<12; bit2 <<= 1, i++)
   {
      if (edgeTable[cubeindex] & bit2)
      {
         // index of each vertex of current edge
         int idx1 = cubeEdges[i][0];
         int idx2 = cubeEdges[i][1];

         Edge &edge = edges[cube.indices[idx1]+cube.indices[idx2]];
         // has intersection already been computed ?
         if (edge.done == false)
         {      
            for (int k=0; k<primlist.size(); k++)
            {
               // is primitive intersecting this edge?
               if (cube.val[idx1][k]*cube.val[idx2][k] <= 0)
               {                  
                  primlist[k]->findIntersection(cube.p[idx1], cube.p[idx2], 
                                                cube.val[idx1][k], cube.val[idx2][k], V, N);

                  // now check if the point found is outside all other primitives
                  bool outside=true;
                  for (int j=0; outside && j<primlist.size(); j++)
                  {
                     // skip the current primitive
                     if (j == k) continue;

                     // skip primitives that do not intersect this edge
                     if (cube.val[idx1][j]*cube.val[idx2][j] > 0) continue;
                     
                     outside = primlist[j]->outside(V);
                  }
                  
                  // look no further?
                  if (outside == true) break;
               }
            }
            vertlist.push_back(V);
            normlist.push_back(N);
            edge.done = true;
            index[i] = edge.vtx_idx = vertlist.size()-1;            
         }          
         else index[i] = edge.vtx_idx;
      }
   }

   // Create the triangles
   for (i=0;triTable[cubeindex][i]!=-1;i+=3) 
   {
      Triangle tri;
      for (int k=0; k<3; k++)
      {
         int idx = triTable[cubeindex][i+k];
         tri.indices[k] = index[idx];
      }
      trilist.push_back(tri);
   }

   return(faceTable[cubeindex]);   
}


void
marchingCubes::computeVoxelFieldValues(GridCell &cell, const idx3 &id)
                        
{
   static unsigned int tab[8][3] = 
   {
      {0,0,0},{1,0,0},{1,0,1},{0,0,1},
      {0,1,0},{1,1,0},{1,1,1},{0,1,1}
   };

   //   cout<<"Voxel coord:"<<endl;
   // test for all 8 vertices of cube
   for (int p=0; p<8; p++)
   {
      int x = tab[p][0];
      int y = tab[p][1];
      int z = tab[p][2];
     
      idx3 id_vtx(id.i+x, id.j+y, id.k+z);
      Voxel &voxel = surf[id_vtx];
    
      // compute, if necessary, pos and field value
      if (voxel.field_values.empty())
      {
         for (int i=0; i<3; i++) 
            voxel.X[i] = LL[i] + id_vtx[i]*res;

         for (int k=0; k<primlist.size(); k++)
            voxel.field_values.push_back(primlist[k]->getDistance(voxel.X));
      }         
#if 0
      std::cout<<"["<<id_vtx<<"]="<<voxel.X<<std::endl;
#endif
      cell.p[p] = voxel.X;
      cell.val[p] = voxel.field_values;
      cell.indices[p] = id_vtx;
   }
}


// voxlist = voxels to test; this list grows continuously until all voxels belonging to 
// the surface have been tested.
// vertlist = resulting list of vertices
// trilist = resulting list of triangles
void
marchingCubes::constructSurface(list<idx3> &voxlist)
{

   VEC D = (UR-LL)/res;
   //   cout<<"Size of the bbox cube:"<<int(D[0])<<","<<int(D[1])<<","<<int(D[2])<<endl;

   while (voxlist.empty() == false)
   {
      idx3 id = voxlist.front();
      //      std::cout<<"Current voxel Id="<<id<<std::endl; getchar();
      if (surf[id].done == false)
      {  
         GridCell cell;
         computeVoxelFieldValues(cell, id);
         surf[id].done = true;
   
         // add new vertices
         if (int faces = polygonise(cell) )
         {
            // possibly add new voxel
            if (faces & 1)
            {
               idx3 id_adj(id); id_adj.k--;
               if (id_adj.k>=0 && surf[id_adj].done == false) voxlist.push_back(id_adj);
            }
            // possibly add new voxel
            if (faces & 2)
            {
               idx3 id_adj(id); id_adj.k++;
               if (id_adj.k<=(int)D[2] && surf[id_adj].done == false) voxlist.push_back(id_adj);
            }

            // possibly add new voxel
            if (faces & 4)
            {
               idx3 id_adj(id); id_adj.j--;
               if (id_adj.j>=0 && surf[id_adj].done == false) voxlist.push_back(id_adj);
            }

            // possibly add new voxel
            if (faces & 8)
            {
               idx3 id_adj(id); id_adj.j++;
               if (id_adj.j<=(int)D[1] && surf[id_adj].done == false) voxlist.push_back(id_adj);
            }   
 
            // possibly add new voxel
            if (faces & 16)
            {
               idx3 id_adj(id); id_adj.i--;
               if (id_adj.i>=0 && surf[id_adj].done == false) voxlist.push_back(id_adj);
            }

            // possibly add new voxel
            if (faces & 32)
            {
               idx3 id_adj(id); id_adj.i++;
               if (id_adj.i<=(int)D[0] && surf[id_adj].done == false) voxlist.push_back(id_adj);
            }   
         }
      }
      voxlist.pop_front();
   }
}


marchingCubes::marchingCubes( vector<Primitive*> primitives, REAL resolution)/*
   : vertlist(), normlist(), trilist(), surf(), edges(),
     LL(dimspace), UR(dimspace), res(resolution), primlist(primitives)*/
{  
   clock_t start=clock(), finish;
   
   // find bounding box
   BBox bb;
   for (int k=0; k<primlist.size(); k++) bb.extendBy( primlist[k]->getBoundingBox() );
   bb.scale(1.05);
   LL = bb.getMin();
   UR = bb.getMax();  
#if 0
   std::cout<<"LL "<<LL<<std::endl;
   std::cout<<"UR "<<UR<<std::endl;
   getchar();
#endif
     
   if (VEC(UR-LL).norm2() < epsilon){ 
     cout << " LL = " << LL << " UR = " << UR << endl;
     cout << "Error: found a bounding box with dimensions equal to 0 for marching cubes\n";
   }
   else
   {
      // find seed cube  by ray-casting
      idx3 id_start = findSeedVoxel();
      list<idx3> voxlist(1, id_start);
      constructSurface(voxlist);

      finish = clock();
      cout << "voxelisation took " << 
         (double)(finish-start)/(double)CLOCKS_PER_SEC << " seconds to execute\n";
      cout << "NVerts  :"<<vertlist.size()<<endl;
      cout << "NFacets :"<<trilist.size()<<endl; 
   }
}
//}

---