BALL  1.4.1
contourSurface.h
Go to the documentation of this file.
00001 // -*- Mode: C++; tab-width: 2; -*-
00002 // vi: set ts=2:
00003 //
00004 
00005 #ifndef BALL_DATATYPE_CONTOURSURFACE_H
00006 #define BALL_DATATYPE_CONTOURSURFACE_H
00007 
00008 #ifndef BALL_COMMON_H
00009 # include <BALL/common.h>
00010 #endif
00011 
00012 #ifndef BALL_DATATYPE_REGULARDATA3D_H
00013 # include <BALL/DATATYPE/regularData3D.h>
00014 #endif
00015 
00016 #ifndef BALL_MATHS_SURFACE_H
00017 # include <BALL/MATHS/surface.h>
00018 #endif
00019 
00020 #ifndef BALL_DATATYPE_HASHMAP_H
00021 # include <BALL/DATATYPE/hashMap.h>
00022 #endif
00023 
00024 #include <vector>
00025 #include <math.h>
00026 
00027 #ifdef BALL_HAS_HASH_MAP
00028 namespace BALL_MAP_NAMESPACE
00029 {
00030   template<>
00031   struct hash<std::pair<BALL::Position, BALL::Position> >
00032   {
00033     size_t operator () (const std::pair<BALL::Position, BALL::Position>& p) const
00034     {return (size_t)(p.first + p.second);}
00035   };
00036 }
00037 #endif
00038 
00039 namespace BALL
00040 {
00041   template<>
00042   BALL_EXPORT HashIndex Hash(const std::pair<Position, Position>& p);
00043 
00044   // 
00045   typedef Index FacetArray[256][12];
00046 
00047   // function defined in contourSurface.C to precompute some tables.
00048   BALL_EXPORT extern const FacetArray& getContourSurfaceFacetData(double threshold);
00049 
00056   template <typename T>  
00057   class TContourSurface 
00058     : public Surface
00059   {
00060     public:
00061 
00065     
00068     typedef std::pair<Position, Position> KeyType;
00069 
00073     typedef Vector3 PointType;
00074     
00078     typedef std::vector<std::pair<PointType, std::pair<Position, Position> > > VectorType;
00080 
00084 
00086     TContourSurface();
00087 
00089     TContourSurface(T threshold);
00090 
00092     TContourSurface(const TContourSurface& surface);
00093 
00095     TContourSurface(const TRegularData3D<T>& data, T threshold = 0.0);
00096 
00098     virtual ~TContourSurface();
00100 
00104     
00106     const TContourSurface& operator = (const TContourSurface<T>& surface);
00107 
00109     const TContourSurface<T>& operator << (const TRegularData3D<T>& data);
00110 
00112     virtual void clear();
00114 
00118 
00120     bool operator == (const TContourSurface<T>& surface) const;
00121 
00123 
00124     
00125     protected:
00126 
00132     class Cube
00133     { 
00134       public:
00135 
00136       Cube(const TRegularData3D<T>& grid)
00137         : grid_(&grid),
00138           current_position_(0),
00139           ptr_(0),
00140           spacing_(grid.getSpacing().x, grid.getSpacing().y, grid.getSpacing().z)
00141       {
00142         // Retrieve the number of points in the grid along the x- and y-axes.
00143         Size nx = grid.getSize().x;
00144         Size ny = grid.getSize().y;
00145 
00146         // Compute the offsets in the grid for the eight 
00147         // corners of the cube (in absolute grid indices).
00148         grid_offset_[0] = nx * ny;
00149         grid_offset_[1] = 0;
00150         grid_offset_[2] = 1;
00151         grid_offset_[3] = 1 + nx * ny;
00152         grid_offset_[4] = nx * ny + nx;
00153         grid_offset_[5] = nx;
00154         grid_offset_[6] = nx + 1;
00155         grid_offset_[7] = nx * ny + nx + 1;
00156       }
00157 
00158       void setTo(Position p) 
00159       {
00160         current_position_ = p;
00161 
00162         ptr_ = &(const_cast<TRegularData3D<T>*>(grid_)->getData(current_position_));
00163         for (Position i = 0; i < 8; i++)
00164         {
00165           values[i] = *(ptr_ + grid_offset_[i]);
00166         }
00167       }
00168 
00169       inline Vector3 getOrigin() const
00170       {
00171         return grid_->getCoordinates(current_position_);
00172       }
00173 
00174       inline const Vector3& getSpacing() const
00175       {
00176         return spacing_;
00177       }
00178 
00179       inline Vector3 getCoordinates(Position index) const
00180       {
00181         return grid_->getCoordinates(index);
00182       }
00183 
00185       inline Position getIndex(Position corner) const
00186       {
00187         return current_position_ + grid_offset_[corner];
00188       }
00189 
00190       void shift() 
00191       {
00192         // Shift the cube by one along the x-axis.
00193         current_position_++;
00194         ptr_++;
00195         
00196         // Keep the four old values for x = 1.
00197         values[0] = values[3];
00198         values[1] = values[2];
00199         values[4] = values[7];
00200         values[5] = values[6];
00201         
00202         // Retrieve the four new values.
00203         values[3] = *(ptr_ + grid_offset_[3]);
00204         values[2] = *(ptr_ + grid_offset_[2]);
00205         values[7] = *(ptr_ + grid_offset_[7]);
00206         values[6] = *(ptr_ + grid_offset_[6]);
00207       }
00208 
00210       Position computeTopology(double threshold)  
00211       {
00212         static const Position topology_modifier[8] = {1, 2, 4, 8, 16, 32, 64, 128};
00213 
00214         // The topology is a bitvector constructed by ORing the
00215         // bits corresponding to each of the inside/outside status of
00216         // of each individual corner.
00217         Position topology = 0;
00218         for (Position i = 0; i < 8; i++)
00219         {
00220           topology |= ((values[i] > threshold) ? topology_modifier[i] : 0);
00221         }
00222         
00223         return topology;
00224       }
00225 
00226       // The values at the eight corners of the cube.
00227       double                    values[8];
00228 
00229       protected:
00230 
00231       // A pointer to the grid.
00232       const TRegularData3D<T>*  grid_;
00233 
00234       // The current position of the cube as an absolute index into the grid.
00235       Position                  current_position_;
00236 
00237       // The gridd offsets: what to add to current_index_ to get to the correct
00238       // grid coordinate.
00239       Position                  grid_offset_[8];
00240 
00241       // A pointer into (nasty hack!) the grid itself. For speed's sake.
00242       const T*                  ptr_;
00243     
00244       // The spacing of the grid.
00245       Vector3 spacing_;
00246     };
00247 
00248 
00250     void addTriangles_(Cube& cube, const FacetArray& facet_data);
00251 
00253     void computeTriangles(Size topology, const TRegularData3D<T>& data);
00254 
00256     T threshold_;
00257 
00258     //
00259     HashMap<std::pair<Position, Position>, Position> cut_hash_map_;
00260   };
00261 
00263   typedef TContourSurface<float> ContourSurface;
00264 
00265   template <typename T>
00266   TContourSurface<T>::TContourSurface()
00267     : threshold_(0.0)
00268   {
00269   }
00270 
00271   template <typename T>
00272   TContourSurface<T>::TContourSurface(T threshold)
00273     : threshold_(threshold)
00274   {
00275   }
00276    
00277   template <typename T>
00278   TContourSurface<T>::TContourSurface(const TRegularData3D<T>& data, T threshold)
00279     : threshold_(threshold)
00280   {
00281     this->operator << (data);
00282   }
00283    
00284   template <typename T>
00285   TContourSurface<T>::~TContourSurface()  
00286   {
00287   }
00288 
00289   template <typename T>
00290   TContourSurface<T>::TContourSurface(const TContourSurface<T>& from)
00291     : threshold_(from.threshold_)
00292   {
00293   }
00294 
00295   template <typename T>
00296   void TContourSurface<T>::clear()
00297   {
00298     Surface::clear();
00299     cut_hash_map_.clear();
00300   }
00301 
00302   template <typename T>
00303   const TContourSurface<T>& TContourSurface<T>::operator = (const TContourSurface<T>& data)
00304   {
00305     // Avoid self-assignment
00306     if (&data != this)
00307     {
00308       threshold_ = data.threshold_;
00309     }
00310 
00311     return *this;
00312   }
00313 
00314   template <typename T>
00315   bool TContourSurface<T>:: operator == (const TContourSurface<T>& data) const
00316   {
00317     return ((threshold_    == data.threshold_)
00318              && Surface::operator == (data.data_));
00319   }
00320 
00321   template <typename T>
00322   const TContourSurface<T>& TContourSurface<T>::operator << (const TRegularData3D<T>& data)
00323   {
00324     // Clear the old stuff:
00325     clear();
00326     
00327     
00328     // Marching cube algorithm: construct a contour surface from
00329     // a volume data set.
00330 
00331     // Get the dimensions of the volume data set.
00332     Vector3 origin = data.getOrigin();
00333     Size number_of_cells_x = (Size)data.getSize().x;
00334     Size number_of_cells_y = (Size)data.getSize().y;
00335     Size number_of_cells_z = (Size)data.getSize().z;
00336 
00337     // Precompute the facet data. This depends on the threshold!
00338     const FacetArray& facet_data = getContourSurfaceFacetData(threshold_);
00339 
00340     // We start in the left-front-bottom-most corner of the grid.
00341     Position current_index = 0;
00342     Cube cube(data);
00343     for (Position curr_cell_z = 0; curr_cell_z < (number_of_cells_z - 1); curr_cell_z++)
00344     { 
00345       // Determine the start position in the current XY plane.
00346       current_index = curr_cell_z * number_of_cells_y * number_of_cells_x;
00347 
00348       // Walk along the y-axis....
00349       for (Position curr_cell_y = 0; curr_cell_y < (number_of_cells_y - 1); curr_cell_y++)
00350       {
00351         // Retrieve the cube from the current grid position (the first position along
00352         // along the x-axis).
00353         cube.setTo(current_index);
00354 
00355         // Walk along the x-axis....
00356         for (Position curr_cell_x = 0; (curr_cell_x < (number_of_cells_x - 2)); )
00357         {
00358           // Compute topology, triangles, and add those triangles to the surface.
00359           addTriangles_(cube, facet_data);
00360             
00361           // Done. cube.shift() will now shift the cube
00362           // along the x-axis and efficently retrieve the four new values.
00363           curr_cell_x++;
00364           cube.shift();
00365         }
00366 
00367         // Add the triangles from the last cube position.
00368         addTriangles_(cube, facet_data);
00369   
00370         // Shift the cube by one along the y-axis.
00371         current_index += number_of_cells_x;
00372       }
00373     }
00374 
00375     // Normalize the vertex normals.
00376     for (Position i = 0; i < normal.size(); i++)
00377     {
00378       try
00379       {
00380         normal[i].normalize();
00381       }
00382       catch (...)
00383       {
00384       }
00385     }
00386 
00387     // Return this (stream operator, for command chaining...)
00388     return *this;
00389   }
00390 
00391   template <typename T>
00392   void TContourSurface<T>::addTriangles_
00393     (typename TContourSurface<T>::Cube& cube, const FacetArray& facet_data)
00394   { 
00395     // Some static variables we need below -- since we will
00396     // call this rather often, we would rather want to avoid
00397     // to many ctor/dtor calls.
00398     static Triangle t;
00399     static std::pair<Position, Position> key;
00400     static std::pair<Position, Position> indices;
00401 
00402     // The indices of the corners of a cube's twelve edges.
00403     static const Position edge_indices[12][2] 
00404       = {{1, 0}, {1, 2}, {2, 3}, {0, 3}, {5, 4}, {5, 6},
00405          {6, 7}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
00406         };
00407 
00408     // The index (into Vector3) of the axis along which the
00409     // current edged runs (0: x, 1: y, 2: z).
00410     static const Position edge_axis[12] 
00411       = {2, 0, 2, 0, 2, 0, 2, 0, 1, 1, 1, 1};
00412 
00413     // Retrieve some basic grid properties.
00414     const Vector3& spacing = cube.getSpacing();
00415 
00416     // The indices (into Surface::vertex) of the triangle
00417     // under construction.
00418     TVector3<Position> triangle_vertices;
00419 
00420     // A counter for the number of vertices already in triangle_vertices
00421     Size vertex_counter = 0;
00422     
00423     // Compute the cube's topology
00424     Position topology = cube.computeTopology(threshold_);
00425     if (topology == 0)
00426     {
00427       return;
00428     }
00429 
00430     // Iterate over all 12 edges and determine whether
00431     // there's a cut. 
00432     for (Position i = 0; i < 12; i++) 
00433     { 
00434       // facet_data_ defines whether there's a cut for
00435       // a given topology and a given edge.
00436       Index facet_index = facet_data[topology][i];
00437 
00438       // There's a cut only for values larger than -1
00439       if (facet_index != -1)
00440       { 
00441         // There is a cut -- determine its position along the edge.
00442 
00443         // The axis: x = 0, y = 1, z = 2 -- used as in index into Vector3.
00444         Position edge = edge_axis[facet_index];
00445 
00446         indices.first = edge_indices[facet_index][0];
00447         indices.second = edge_indices[facet_index][1];
00448         key.first = cube.getIndex(indices.first);
00449         key.second = cube.getIndex(indices.second);
00450 
00451         // Check whether we computed this cut already.
00452         if (!cut_hash_map_.has(key))
00453         {
00454           // Compute the position of the cut.
00455           
00456           // Get the position of the d1 point
00457           Vector3 pos = cube.getCoordinates(key.first);
00458 
00459           // Compute the position of the cut along the edge.
00460           const double& d1 = cube.values[indices.first];
00461           const double& d2 = cube.values[indices.second];
00462           pos[edge] += ((double)threshold_ - d1) / (d2 - d1) * spacing[edge];
00463           
00464           // Store it as a triangle vertex.
00465           triangle_vertices[vertex_counter++] = vertex.size();
00466 
00467           // Store the index of the vertex in the hash map under the
00468           // indices of its grid points.
00469           cut_hash_map_.insert(std::pair<std::pair<Position, Position>, Position>(key, (Size)vertex.size()));
00470 
00471           // Create a vertex and a normal (the normal reamins empty for now).
00472           vertex.push_back(pos);
00473           static Vector3 null_normal(0.0, 0.0, 0.0);
00474           normal.push_back(null_normal);
00475         }
00476         else
00477         {
00478           // This one we know already! Retrieve it from the hash map.
00479           triangle_vertices[vertex_counter++] = cut_hash_map_[key];
00480         }
00481         
00482         // For every three vertices, create a new triangle.
00483         if (vertex_counter == 3)
00484         {
00485           // Create a new triangle.
00486           t.v1 = triangle_vertices.x;
00487           t.v2 = triangle_vertices.y;
00488           t.v3 = triangle_vertices.z;
00489           triangle.push_back(t);
00490 
00491           // We can start with the next one.
00492           vertex_counter = 0;
00493 
00494           // Compute the normals: add the triangle
00495           // normals to each of the triangle vertices.
00496           // We will average them out to the correct normals later.
00497           Vector3 h1(vertex[t.v1] - vertex[t.v2]);
00498           Vector3 h2(vertex[t.v3] - vertex[t.v2]);
00499           Vector3 current_normal(h1.y * h2.z - h1.z * h2.y,
00500                                  h1.z * h2.x - h2.z * h1.x,
00501                                  h1.x * h2.y - h1.y * h2.x);
00502           normal[t.v1] += current_normal;
00503           normal[t.v2] += current_normal;
00504           normal[t.v3] += current_normal;
00505         }
00506       }
00507     }
00508   }
00509 }
00510 #endif
00511 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Defines