BALL
1.4.1
|
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