BALL  1.5.0
contourSurface.h
Go to the documentation of this file.
1 // -*- Mode: C++; tab-width: 2; -*-
2 // vi: set ts=2:
3 //
4 
5 #ifndef BALL_DATATYPE_CONTOURSURFACE_H
6 #define BALL_DATATYPE_CONTOURSURFACE_H
7 
8 #ifndef BALL_COMMON_H
9 # include <BALL/common.h>
10 #endif
11 
12 #ifndef BALL_DATATYPE_REGULARDATA3D_H
14 #endif
15 
16 #ifndef BALL_MATHS_SURFACE_H
17 # include <BALL/MATHS/surface.h>
18 #endif
19 
20 #ifndef BALL_DATATYPE_HASHMAP_H
21 # include <BALL/DATATYPE/hashMap.h>
22 #endif
23 
24 #include <vector>
25 
26 namespace BALL
27 {
28  template<>
29  BALL_EXPORT HashIndex Hash(const std::pair<Position, Position>& p);
30 
31  //
32  typedef Index FacetArray[256][12];
33 
34  // function defined in contourSurface.C to precompute some tables.
35  BALL_EXPORT extern const FacetArray& getContourSurfaceFacetData(double threshold);
36 
43  template <typename T>
45  : public Surface
46  {
47  public:
48 
52 
55  typedef std::pair<Position, Position> KeyType;
56 
60  typedef Vector3 PointType;
61 
65  typedef std::vector<std::pair<PointType, std::pair<Position, Position> > > VectorType;
67 
71 
74 
76  TContourSurface(T threshold);
77 
79  TContourSurface(const TContourSurface& surface);
80 
82  TContourSurface(const TRegularData3D<T>& data, T threshold = 0.0);
83 
85  virtual ~TContourSurface();
87 
91 
93  const TContourSurface& operator = (const TContourSurface<T>& surface);
94 
97 
99  virtual void clear();
101 
105 
107  bool operator == (const TContourSurface<T>& surface) const;
108 
110 
111 
112  protected:
113 
119  class Cube
120  {
121  public:
122 
123  Cube(const TRegularData3D<T>& grid)
124  : grid_(&grid),
126  ptr_(0),
127  spacing_(grid.getSpacing().x, grid.getSpacing().y, grid.getSpacing().z)
128  {
129  // Retrieve the number of points in the grid along the x- and y-axes.
130  Size nx = grid.getSize().x;
131  Size ny = grid.getSize().y;
132 
133  // Compute the offsets in the grid for the eight
134  // corners of the cube (in absolute grid indices).
135  grid_offset_[0] = nx * ny;
136  grid_offset_[1] = 0;
137  grid_offset_[2] = 1;
138  grid_offset_[3] = 1 + nx * ny;
139  grid_offset_[4] = nx * ny + nx;
140  grid_offset_[5] = nx;
141  grid_offset_[6] = nx + 1;
142  grid_offset_[7] = nx * ny + nx + 1;
143  }
144 
145  void setTo(Position p)
146  {
147  current_position_ = p;
148 
149  ptr_ = &(const_cast<TRegularData3D<T>*>(grid_)->getData(current_position_));
150  for (Position i = 0; i < 8; i++)
151  {
152  values[i] = *(ptr_ + grid_offset_[i]);
153  }
154  }
155 
156  inline Vector3 getOrigin() const
157  {
159  }
160 
161  inline const Vector3& getSpacing() const
162  {
163  return spacing_;
164  }
165 
166  inline Vector3 getCoordinates(Position index) const
167  {
168  return grid_->getCoordinates(index);
169  }
170 
172  inline Position getIndex(Position corner) const
173  {
174  return current_position_ + grid_offset_[corner];
175  }
176 
177  void shift()
178  {
179  // Shift the cube by one along the x-axis.
181  ptr_++;
182 
183  // Keep the four old values for x = 1.
184  values[0] = values[3];
185  values[1] = values[2];
186  values[4] = values[7];
187  values[5] = values[6];
188 
189  // Retrieve the four new values.
190  values[3] = *(ptr_ + grid_offset_[3]);
191  values[2] = *(ptr_ + grid_offset_[2]);
192  values[7] = *(ptr_ + grid_offset_[7]);
193  values[6] = *(ptr_ + grid_offset_[6]);
194  }
195 
197  Position computeTopology(double threshold)
198  {
199  static const Position topology_modifier[8] = {1, 2, 4, 8, 16, 32, 64, 128};
200 
201  // The topology is a bitvector constructed by ORing the
202  // bits corresponding to each of the inside/outside status of
203  // of each individual corner.
204  Position topology = 0;
205  for (Position i = 0; i < 8; i++)
206  {
207  topology |= ((values[i] > threshold) ? topology_modifier[i] : 0);
208  }
209 
210  return topology;
211  }
212 
213  // The values at the eight corners of the cube.
214  double values[8];
215 
216  protected:
217 
218  // A pointer to the grid.
220 
221  // The current position of the cube as an absolute index into the grid.
223 
224  // The gridd offsets: what to add to current_index_ to get to the correct
225  // grid coordinate.
227 
228  // A pointer into (nasty hack!) the grid itself. For speed's sake.
229  const T* ptr_;
230 
231  // The spacing of the grid.
233  };
234 
235 
237  void addTriangles_(Cube& cube, const FacetArray& facet_data);
238 
240  void computeTriangles(Size topology, const TRegularData3D<T>& data);
241 
244 
245  //
247  };
248 
251 
252  template <typename T>
254  : threshold_(0.0)
255  {
256  }
257 
258  template <typename T>
260  : threshold_(threshold)
261  {
262  }
263 
264  template <typename T>
266  : threshold_(threshold)
267  {
268  this->operator << (data);
269  }
270 
271  template <typename T>
273  {
274  }
275 
276  template <typename T>
278  : threshold_(from.threshold_)
279  {
280  }
281 
282  template <typename T>
284  {
285  Surface::clear();
286  cut_hash_map_.clear();
287  }
288 
289  template <typename T>
291  {
292  // Avoid self-assignment
293  if (&data != this)
294  {
295  threshold_ = data.threshold_;
296  }
297 
298  return *this;
299  }
300 
301  template <typename T>
303  {
304  return ((threshold_ == data.threshold_)
305  && Surface::operator == (data.data_));
306  }
307 
308  template <typename T>
310  {
311  // Clear the old stuff:
312  clear();
313 
314 
315  // Marching cube algorithm: construct a contour surface from
316  // a volume data set.
317 
318  // Get the dimensions of the volume data set.
319  Vector3 origin = data.getOrigin();
320  Size number_of_cells_x = (Size)data.getSize().x;
321  Size number_of_cells_y = (Size)data.getSize().y;
322  Size number_of_cells_z = (Size)data.getSize().z;
323 
324  // Precompute the facet data. This depends on the threshold!
325  const FacetArray& facet_data = getContourSurfaceFacetData(threshold_);
326 
327  // We start in the left-front-bottom-most corner of the grid.
328  Position current_index = 0;
329  Cube cube(data);
330  for (Position curr_cell_z = 0; curr_cell_z < (number_of_cells_z - 1); curr_cell_z++)
331  {
332  // Determine the start position in the current XY plane.
333  current_index = curr_cell_z * number_of_cells_y * number_of_cells_x;
334 
335  // Walk along the y-axis....
336  for (Position curr_cell_y = 0; curr_cell_y < (number_of_cells_y - 1); curr_cell_y++)
337  {
338  // Retrieve the cube from the current grid position (the first position along
339  // along the x-axis).
340  cube.setTo(current_index);
341 
342  // Walk along the x-axis....
343  for (Position curr_cell_x = 0; (curr_cell_x < (number_of_cells_x - 2)); )
344  {
345  // Compute topology, triangles, and add those triangles to the surface.
346  addTriangles_(cube, facet_data);
347 
348  // Done. cube.shift() will now shift the cube
349  // along the x-axis and efficently retrieve the four new values.
350  curr_cell_x++;
351  cube.shift();
352  }
353 
354  // Add the triangles from the last cube position.
355  addTriangles_(cube, facet_data);
356 
357  // Shift the cube by one along the y-axis.
358  current_index += number_of_cells_x;
359  }
360  }
361 
362  // Normalize the vertex normals.
363  for (Position i = 0; i < normal.size(); i++)
364  {
365  try
366  {
367  normal[i].normalize();
368  }
369  catch (...)
370  {
371  }
372  }
373 
374  // Return this (stream operator, for command chaining...)
375  return *this;
376  }
377 
378  template <typename T>
380  (typename TContourSurface<T>::Cube& cube, const FacetArray& facet_data)
381  {
382  // Some static variables we need below -- since we will
383  // call this rather often, we would rather want to avoid
384  // to many ctor/dtor calls.
385  static Triangle t;
386  static std::pair<Position, Position> key;
387  static std::pair<Position, Position> indices;
388 
389  // The indices of the corners of a cube's twelve edges.
390  static const Position edge_indices[12][2]
391  = {{1, 0}, {1, 2}, {2, 3}, {0, 3}, {5, 4}, {5, 6},
392  {6, 7}, {4, 7}, {0, 4}, {1, 5}, {2, 6}, {3, 7}
393  };
394 
395  // The index (into Vector3) of the axis along which the
396  // current edged runs (0: x, 1: y, 2: z).
397  static const Position edge_axis[12]
398  = {2, 0, 2, 0, 2, 0, 2, 0, 1, 1, 1, 1};
399 
400  // Retrieve some basic grid properties.
401  const Vector3& spacing = cube.getSpacing();
402 
403  // The indices (into Surface::vertex) of the triangle
404  // under construction.
405  TVector3<Position> triangle_vertices;
406 
407  // A counter for the number of vertices already in triangle_vertices
408  Size vertex_counter = 0;
409 
410  // Compute the cube's topology
411  Position topology = cube.computeTopology(threshold_);
412  if (topology == 0)
413  {
414  return;
415  }
416 
417  // Iterate over all 12 edges and determine whether
418  // there's a cut.
419  for (Position i = 0; i < 12; i++)
420  {
421  // facet_data_ defines whether there's a cut for
422  // a given topology and a given edge.
423  Index facet_index = facet_data[topology][i];
424 
425  // There's a cut only for values larger than -1
426  if (facet_index != -1)
427  {
428  // There is a cut -- determine its position along the edge.
429 
430  // The axis: x = 0, y = 1, z = 2 -- used as in index into Vector3.
431  Position edge = edge_axis[facet_index];
432 
433  indices.first = edge_indices[facet_index][0];
434  indices.second = edge_indices[facet_index][1];
435  key.first = cube.getIndex(indices.first);
436  key.second = cube.getIndex(indices.second);
437 
438  // Check whether we computed this cut already.
439  if (!cut_hash_map_.has(key))
440  {
441  // Compute the position of the cut.
442 
443  // Get the position of the d1 point
444  Vector3 pos = cube.getCoordinates(key.first);
445 
446  // Compute the position of the cut along the edge.
447  const double& d1 = cube.values[indices.first];
448  const double& d2 = cube.values[indices.second];
449  pos[edge] += ((double)threshold_ - d1) / (d2 - d1) * spacing[edge];
450 
451  // Store it as a triangle vertex.
452  triangle_vertices[vertex_counter++] = vertex.size();
453 
454  // Store the index of the vertex in the hash map under the
455  // indices of its grid points.
456  cut_hash_map_.insert(std::pair<std::pair<Position, Position>, Position>(key, (Size)vertex.size()));
457 
458  // Create a vertex and a normal (the normal reamins empty for now).
459  vertex.push_back(pos);
460  static Vector3 null_normal(0.0, 0.0, 0.0);
461  normal.push_back(null_normal);
462  }
463  else
464  {
465  // This one we know already! Retrieve it from the hash map.
466  triangle_vertices[vertex_counter++] = cut_hash_map_[key];
467  }
468 
469  // For every three vertices, create a new triangle.
470  if (vertex_counter == 3)
471  {
472  // Create a new triangle.
473  t.v1 = triangle_vertices.x;
474  t.v2 = triangle_vertices.y;
475  t.v3 = triangle_vertices.z;
476  triangle.push_back(t);
477 
478  // We can start with the next one.
479  vertex_counter = 0;
480 
481  // Compute the normals: add the triangle
482  // normals to each of the triangle vertices.
483  // We will average them out to the correct normals later.
484  Vector3 h1(vertex[t.v1] - vertex[t.v2]);
485  Vector3 h2(vertex[t.v3] - vertex[t.v2]);
486  Vector3 current_normal(h1.y * h2.z - h1.z * h2.y,
487  h1.z * h2.x - h2.z * h1.x,
488  h1.x * h2.y - h1.y * h2.x);
489  normal[t.v1] += current_normal;
490  normal[t.v2] += current_normal;
491  normal[t.v3] += current_normal;
492  }
493  }
494  }
495  }
496 }
497 #endif
498 
BALL::TVector3::x
T x
Definition: vector3.h:489
BALL::TVector3::y
T y
Definition: vector3.h:493
BALL::TContourSurface::cut_hash_map_
HashMap< std::pair< Position, Position >, Position > cut_hash_map_
Definition: contourSurface.h:246
BALL_INDEX_TYPE
BALL::TContourSurface::Cube::getSpacing
const Vector3 & getSpacing() const
Definition: contourSurface.h:161
BALL_EXPORT
#define BALL_EXPORT
Definition: COMMON/global.h:50
BALL::TContourSurface::Cube::grid_
const TRegularData3D< T > * grid_
Definition: contourSurface.h:219
BALL::TContourSurface::operator==
bool operator==(const TContourSurface< T > &surface) const
Equality operator.
Definition: contourSurface.h:302
BALL::getContourSurfaceFacetData
const BALL_EXPORT FacetArray & getContourSurfaceFacetData(double threshold)
hashMap.h
BALL::TSurface< float >
BALL::Size
BALL_SIZE_TYPE Size
Definition: COMMON/global.h:114
BALL::HashIndex
BALL_SIZE_TYPE HashIndex
Definition: COMMON/global.h:131
BALL::HashMap
HashMap class based on the STL map (containing serveral convenience functions)
Definition: hashMap.h:73
BALL::TContourSurface::clear
virtual void clear()
Clear method.
Definition: contourSurface.h:283
BALL::TRegularData3D< T >
BALL::TVector3< float >
BALL::TContourSurface::Cube::ptr_
const T * ptr_
Definition: contourSurface.h:229
BALL_SIZE_TYPE
double
BALL::TContourSurface::Cube::Cube
Cube(const TRegularData3D< T > &grid)
Definition: contourSurface.h:123
BALL::TContourSurface::Cube::grid_offset_
Position grid_offset_[8]
Definition: contourSurface.h:226
BALL
Definition: constants.h:12
BALL::TContourSurface::TContourSurface
TContourSurface()
Default constructor.
Definition: contourSurface.h:253
BALL::TContourSurface::operator=
const TContourSurface & operator=(const TContourSurface< T > &surface)
Assignment operator.
Definition: contourSurface.h:290
BALL::TContourSurface::Cube::shift
void shift()
Definition: contourSurface.h:177
BALL::TRegularData3D::getSize
const IndexType & getSize() const
Definition: regularData3D.h:306
BALL::TContourSurface::VectorType
std::vector< std::pair< PointType, std::pair< Position, Position > > > VectorType
Definition: contourSurface.h:65
BALL::TContourSurface::KeyType
std::pair< Position, Position > KeyType
Definition: contourSurface.h:55
BALL::FacetArray
Index FacetArray[256][12]
Definition: contourSurface.h:32
BALL::TContourSurface::Cube::getOrigin
Vector3 getOrigin() const
Definition: contourSurface.h:156
BALL::TContourSurface::addTriangles_
void addTriangles_(Cube &cube, const FacetArray &facet_data)
Definition: contourSurface.h:380
BALL::TContourSurface
Definition: contourSurface.h:44
BALL::TContourSurface::Cube::getCoordinates
Vector3 getCoordinates(Position index) const
Definition: contourSurface.h:166
BALL::TContourSurface::computeTriangles
void computeTriangles(Size topology, const TRegularData3D< T > &data)
BALL::TSurface::Triangle
Definition: surface.h:39
BALL::Hash
HashIndex Hash(const T &key)
Definition: hash.h:47
BALL::TRegularData3D::getOrigin
const CoordinateType & getOrigin() const
Definition: regularData3D.h:312
BALL::TContourSurface::Cube::setTo
void setTo(Position p)
Definition: contourSurface.h:145
BALL::TContourSurface::~TContourSurface
virtual ~TContourSurface()
Destructor.
Definition: contourSurface.h:272
BALL::TSurface::Triangle::v2
Index v2
Definition: surface.h:43
BALL::TSurface::Triangle::v1
Index v1
Definition: surface.h:42
BALL::TSurface::Triangle::v3
Index v3
Definition: surface.h:44
BALL::TContourSurface::PointType
Vector3 PointType
Definition: contourSurface.h:60
BALL::TContourSurface::operator<<
const TContourSurface< T > & operator<<(const TRegularData3D< T > &data)
Create a contour surface from a given data set.
Definition: contourSurface.h:309
BALL::TContourSurface::Cube::getIndex
Position getIndex(Position corner) const
Return the absolute grid position for a given corner.
Definition: contourSurface.h:172
common.h
BALL::TVector3::z
T z
Definition: vector3.h:497
regularData3D.h
BALL::TContourSurface::Cube
Definition: contourSurface.h:119
BALL::TContourSurface::Cube::spacing_
Vector3 spacing_
Definition: contourSurface.h:232
BALL::TSurface< float >::clear
void clear()
Definition: surface.h:220
BALL::ContourSurface
TContourSurface< float > ContourSurface
Default type.
Definition: contourSurface.h:250
BALL::TRegularData3D::getCoordinates
CoordinateType getCoordinates(const IndexType &index) const
surface.h
BALL::TContourSurface::Cube::values
double values[8]
Definition: contourSurface.h:214
BALL::TContourSurface::Cube::current_position_
Position current_position_
Definition: contourSurface.h:222
BALL::TContourSurface::threshold_
T threshold_
The threshold separating inside and outside.
Definition: contourSurface.h:243
BALL::TContourSurface::Cube::computeTopology
Position computeTopology(double threshold)
Compute the topology code for the current cube.
Definition: contourSurface.h:197