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
22#endif
23
24#include <vector>
25
26namespace 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
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
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
146 {
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 {
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
Definition: constants.h:13
Index FacetArray[256][12]
BALL_SIZE_TYPE Size
TContourSurface< float > ContourSurface
Default type.
HashIndex Hash(const T &key)
Definition: hash.h:47
BALL_EXPORT const FacetArray & getContourSurfaceFacetData(double threshold)
T threshold_
The threshold separating inside and outside.
std::vector< std::pair< PointType, std::pair< Position, Position > > > VectorType
void computeTriangles(Size topology, const TRegularData3D< T > &data)
HashMap< std::pair< Position, Position >, Position > cut_hash_map_
void addTriangles_(Cube &cube, const FacetArray &facet_data)
const TContourSurface< T > & operator<<(const TRegularData3D< T > &data)
Create a contour surface from a given data set.
const TContourSurface & operator=(const TContourSurface< T > &surface)
Assignment operator.
std::pair< Position, Position > KeyType
bool operator==(const TContourSurface< T > &surface) const
Equality operator.
TContourSurface()
Default constructor.
virtual void clear()
Clear method.
virtual ~TContourSurface()
Destructor.
const TRegularData3D< T > * grid_
const Vector3 & getSpacing() const
Cube(const TRegularData3D< T > &grid)
Position getIndex(Position corner) const
Return the absolute grid position for a given corner.
Position computeTopology(double threshold)
Compute the topology code for the current cube.
Vector3 getCoordinates(Position index) const
HashMap class based on the STL map (containing serveral convenience functions)
Definition: hashMap.h:74
const IndexType & getSize() const
const CoordinateType & getOrigin() const
CoordinateType getCoordinates(const IndexType &index) const
#define BALL_EXPORT
Definition: COMMON/global.h:50