cgv
marching_cubes.h
1 #pragma once
2 
3 #include <vector>
4 #include <deque>
5 #include <cgv/utils/progression.h>
6 #include <cgv/math/fvec.h>
7 #include <cgv/math/mfunc.h>
8 #include <cgv/media/axis_aligned_box.h>
9 #include <cgv/media/mesh/streaming_mesh.h>
10 
11 #include <cgv/media/lib_begin.h>
12 
13 namespace cgv {
14  namespace media {
15  namespace mesh {
16 
17 extern CGV_API int get_nr_cube_triangles(int idx);
18 extern CGV_API void put_cube_triangle(int idx, int t, int& vi, int& vj, int& vk);
19 
20 
22 template <typename T>
23 struct slice_info
24 {
25  unsigned int resx, resy;
26  std::vector<T> values;
27  std::vector<bool> flags;
28  std::vector<int> indices;
30  slice_info(unsigned int _resx, unsigned int _resy) : resx(_resx), resy(_resy) {
31  unsigned int n = resx*resy;
32  values.resize(n);
33  flags.resize(n);
34  indices.resize(4 * n);
35  }
37  void init() {
38  int n = resx*resy;
39  for (int i = 0; i < 4 * n; ++i)
40  indices[i] = -1;
41  }
43  bool flag(int x, int y) const { return flags[y*resx + x]; }
45  int get_bit_code(int x, int y, int step = 1) const {
46  int i = y*resx + x;
47  return (flags[i] ? 1 : 0) + (flags[i + step] ? 2 : 0) +
48  (flags[i + step*(resx + 1)] ? 4 : 0) + (flags[i + step*resx] ? 8 : 0);
49  }
51  const T& value(int x, int y) const { return values[y*resx + x]; }
52  T& value(int x, int y) { return values[y*resx + x]; }
54  void set_value(int x, int y, T value, T iso_value) {
55  int i = y*resx + x;
56  values[i] = value;
57  flags[i] = value > iso_value;
58  }
60  const int& index(int x, int y, int e) const { return indices[4 * (y*resx + x) + e]; }
61  int& index(int x, int y, int e) { return indices[4 * (y*resx + x) + e]; }
63  const int& snap_index(int x, int y) const { return indices[4 * (y*resx + x) + 3]; }
64  int& snap_index(int x, int y) { return indices[4 * (y*resx + x) + 3]; }
65 };
66 
68 template <typename X, typename T>
69 class marching_cubes_base : public streaming_mesh<X>
70 {
71 public:
72  typedef streaming_mesh<X> base_type;
74  typedef cgv::math::fvec<X,3> pnt_type;
76  typedef cgv::math::fvec<X,3> vec_type;
77 private:
78  pnt_type p;
79  vec_type d;
80  T iso_value;
81 protected:
82  X epsilon;
83  X grid_epsilon;
84 public:
86  marching_cubes_base(streaming_mesh_callback_handler* _smcbh,
87  const X& _grid_epsilon = 0.01f,
88  const X& _epsilon = 1e-6f) : epsilon(_epsilon), grid_epsilon(_grid_epsilon)
89  {
90  base_type::set_callback_handler(_smcbh);
91  }
93  void construct_vertex(slice_info<T> *info_ptr_1, int i_1, int j_1, int e,
94  slice_info<T> *info_ptr_2, int i_2, int j_2)
95  {
96  // read values at edge ends
97  T v_1 = info_ptr_1->value(i_1, j_1);
98  T v_2 = info_ptr_2->value(i_2, j_2);
99  // from values compute affin location
100  X f = (fabs(v_2 - v_1) > epsilon) ? (X)(iso_value - v_1) / (v_2 - v_1) : (X) 0.5;
101  // check whether to snap to edge start
102  int vi = base_type::get_nr_vertices();
103  pnt_type q = p;
104  if (f < grid_epsilon) {
105  int vj = info_ptr_1->snap_index(i_1, j_1);
106  if (vj != -1) {
107  info_ptr_1->index(i_1, j_1, e) = vj;
108  return;
109  }
110  info_ptr_1->snap_index(i_1, j_1) = vi;
111  q(e) -= d(e);
112  }
113  else if (1 - f < grid_epsilon) {
114  int vj = info_ptr_2->snap_index(i_2, j_2);
115  if (vj != -1) {
116  info_ptr_1->index(i_1, j_1, e) = vj;
117  return;
118  }
119  info_ptr_2->snap_index(i_2, j_2) = vi;
120  }
121  else
122  q(e) -= (1 - f)*d(e);
123 
124  info_ptr_1->index(i_1, j_1, e) = vi;
125  this->new_vertex(q);
126  }
127 
129  template <typename Eval, typename Valid>
130  void extract_impl(const T& _iso_value,
131  const axis_aligned_box<X, 3>& box,
132  unsigned int resx, unsigned int resy, unsigned int resz,
133  const Eval& eval, const Valid& valid, bool show_progress = false)
134  {
135  // prepare private members
136  p = box.get_min_pnt();
137  d = box.get_extent();
138  d(0) /= (resx - 1); d(1) /= (resy - 1); d(2) /= (resz - 1);
139  iso_value = _iso_value;
140 
141  // prepare progression
142  cgv::utils::progression prog;
143  if (show_progress) prog.init("extraction", resz, 10);
144 
145  // construct two slice infos
146  slice_info<T> slice_info_1(resx, resy), slice_info_2(resx, resy);
147  slice_info<T> *slice_info_ptrs[2] = { &slice_info_1, &slice_info_2 };
148 
149  // iterate through all slices
150  unsigned int nr_vertices[3] = { 0, 0, 0 };
151  unsigned int i, j, k, n;
152  for (k = 0; k < resz; ++k, p(2) += d(2)) {
153  n = (int)base_type::get_nr_vertices();
154  // evaluate function on next slice and construct slice interior vertices
155  slice_info<T> *info_ptr = slice_info_ptrs[k & 1];
156  info_ptr->init();
157  for (j = 0, p(1) = box.get_min_pnt()(1); j < resy; ++j, p(1) += d(1))
158  for (i = 0, p(0) = box.get_min_pnt()(0); i < resx; ++i, p(0) += d(0)) {
159  T v = eval(i, j, k, p);
160  info_ptr->set_value(i, j, v, iso_value);
161  if (valid(v)) {
162  if (i > 0 && info_ptr->flag(i - 1, j) != info_ptr->flag(i, j) && valid(info_ptr->value(i - 1, j)))
163  construct_vertex(info_ptr, i - 1, j, 0, info_ptr, i, j);
164  if (j > 0 && info_ptr->flag(i, j - 1) != info_ptr->flag(i, j) && valid(info_ptr->value(i, j - 1)))
165  construct_vertex(info_ptr, i, j - 1, 1, info_ptr, i, j);
166  }
167  }
168  // show progression
169  if (show_progress)
170  prog.step();
171  // if this is the first considered slice, construct the next one
172  if (k != 0) {
173  // get info of previous slice
174  slice_info<T> *prev_info_ptr = slice_info_ptrs[1 - (k & 1)];
175  // construct vertices on edges between previous and new slice
176  for (j = 0, p(1) = box.get_min_pnt()(1); j < resy; ++j, p(1) += d(1))
177  for (i = 0, p(0) = box.get_min_pnt()(0); i < resx; ++i, p(0) += d(0))
178  if (prev_info_ptr->flag(i, j) != info_ptr->flag(i, j) && valid(prev_info_ptr->value(i, j)) && valid(info_ptr->value(i, j)))
179  construct_vertex(prev_info_ptr, i, j, 2, info_ptr, i, j);
180 
181  // construct triangles
182  for (j = 0; j < resy - 1; ++j) {
183  for (i = 0; i < resx - 1; ++i) {
184  // compute the bit index for the current cube
185  int idx = prev_info_ptr->get_bit_code(i, j) +
186  16 * info_ptr->get_bit_code(i, j);
187  // skip empty cubes
188  if (idx == 0 || idx == 255)
189  continue;
190  // set edge vertices
191  int vis[12] = {
192  prev_info_ptr->index(i, j, 0),
193  prev_info_ptr->index(i + 1, j, 1),
194  prev_info_ptr->index(i, j + 1, 0),
195  prev_info_ptr->index(i, j, 1),
196  info_ptr->index(i, j, 0),
197  info_ptr->index(i + 1, j, 1),
198  info_ptr->index(i, j + 1, 0),
199  info_ptr->index(i, j, 1),
200  prev_info_ptr->index(i, j, 2),
201  prev_info_ptr->index(i + 1, j, 2),
202  prev_info_ptr->index(i, j + 1, 2),
203  prev_info_ptr->index(i + 1, j + 1, 2)
204  };
205  // lookup triangles and construct them
206  int n = get_nr_cube_triangles(idx);
207  for (int t = 0; t < n; ++t) {
208  int vi, vj, vk;
209  put_cube_triangle(idx, t, vi, vj, vk);
210  vi = vis[vi];
211  vj = vis[vj];
212  vk = vis[vk];
213  if (vi == -1 || vj == -1 || vk == -1)
214  continue;
215  if ((vi != vj) && (vi != vk) && (vj != vk))
216  base_type::new_triangle(vk, vj, vi);
217  }
218  }
219  }
220  }
221  n = (int)base_type::get_nr_vertices() - n;
222  nr_vertices[k % 3] = n;
223  n = nr_vertices[(k + 2) % 3];
224  if (n > 0)
225  base_type::drop_vertices(n);
226  }
227  }
228 };
229 
230 template <typename T>
231 struct always_valid
232 {
233  bool operator () (const T& v) const { return true; }
234 };
235 
236 #include <limits>
237 
238 template <typename T>
239 struct not_max_valid
240 {
241  bool operator () (const T& v) const { return v != std::numeric_limits<T>::max(); }
242 };
243 
244 template <typename T>
245 struct not_nan_valid
246 {
247  bool operator () (const T& v) const { return !std::isnan(v); }
248 };
249 
251 template <typename X, typename T>
252 class marching_cubes : public marching_cubes_base<X, T>
253 {
254 public:
255  // GCC does not examine the base class scope for templates,
256  // so they must be explicitely redefined
257  typedef typename marching_cubes_base<X, T>::base_type base_type;
259  typedef typename marching_cubes_base<X, T>::pnt_type pnt_type;
261  typedef typename marching_cubes_base<X, T>::vec_type vec_type;
262 
263  const cgv::math::v3_func<X, T>& func;
265  marching_cubes(const cgv::math::v3_func<X, T>& _func,
266  streaming_mesh_callback_handler* _smcbh,
267  const X& _grid_epsilon = 0.01f,
268  const X& _epsilon = 1e-6f) : marching_cubes_base<X, T>(_smcbh, _grid_epsilon, _epsilon), func(_func)
269  {
270  }
271  T operator () (unsigned i, unsigned j, unsigned k, const pnt_type& p) const {
272  return func.evaluate(p.to_vec());
273  }
274  void extract(const T& _iso_value,
275  const axis_aligned_box<X, 3>& box,
276  unsigned int resx, unsigned int resy, unsigned int resz,
277  bool show_progress = false)
278  {
279  always_valid<T> valid;
280  this->extract_impl(_iso_value, box, resx, resy, resz, *this, valid, show_progress);
281  }
282 };
283 
284  }
285  }
286 }
287 
288 #include <cgv/config/lib_end.h>
cgv::media::mesh::streaming_mesh< X >::set_callback_handler
void set_callback_handler(streaming_mesh_callback_handler *_smcbh)
set a new callback handler
Definition: streaming_mesh.h:47
cgv::media::mesh::marching_cubes_base::construct_vertex
void construct_vertex(slice_info< T > *info_ptr_1, int i_1, int j_1, int e, slice_info< T > *info_ptr_2, int i_2, int j_2)
construct a new vertex on an edge
Definition: marching_cubes.h:93
cgv::math::fvec< X, 3 >
cgv::media::mesh::marching_cubes
class used to perform the marching cubes algorithm
Definition: marching_cubes.h:253
cgv::math::v3_func
Definition: mfunc.h:63
cgv::media::mesh::streaming_mesh_callback_handler
Definition: streaming_mesh.h:13
cgv::media::mesh::slice_info::set_value
void set_value(int x, int y, T value, T iso_value)
set value and flag at once
Definition: marching_cubes.h:54
cgv::media::mesh::marching_cubes_base::vec_type
cgv::math::fvec< X, 3 > vec_type
vectors must have three components
Definition: marching_cubes.h:76
cgv::media::mesh::streaming_mesh< X >::new_vertex
unsigned int new_vertex(const pnt_type &p)
add a new vertex with the given location and call the callback of the callback handler
Definition: streaming_mesh.h:80
cgv::media::axis_aligned_box::get_extent
fvec_type get_extent() const
return a vector with the extents in the different dimensions
Definition: axis_aligned_box.h:80
cgv::media::mesh::marching_cubes::marching_cubes
marching_cubes(const cgv::math::v3_func< X, T > &_func, streaming_mesh_callback_handler *_smcbh, const X &_grid_epsilon=0.01f, const X &_epsilon=1e-6f)
construct marching cubes object
Definition: marching_cubes.h:265
cgv::media::axis_aligned_box
Definition: axis_aligned_box.h:15
cgv::media::mesh::marching_cubes_base::extract_impl
void extract_impl(const T &_iso_value, const axis_aligned_box< X, 3 > &box, unsigned int resx, unsigned int resy, unsigned int resz, const Eval &eval, const Valid &valid, bool show_progress=false)
extract iso surface and send triangles to marching cubes handler
Definition: marching_cubes.h:130
cgv::media::mesh::marching_cubes_base
class used to perform the marching cubes algorithm
Definition: marching_cubes.h:70
cgv::media::mesh::marching_cubes_base::marching_cubes_base
marching_cubes_base(streaming_mesh_callback_handler *_smcbh, const X &_grid_epsilon=0.01f, const X &_epsilon=1e-6f)
construct marching cubes object
Definition: marching_cubes.h:86
cgv::media::mesh::marching_cubes::vec_type
marching_cubes_base< X, T >::vec_type vec_type
vectors must have three components
Definition: marching_cubes.h:261
cgv::media::mesh::marching_cubes::pnt_type
marching_cubes_base< X, T >::pnt_type pnt_type
points must have three components
Definition: marching_cubes.h:259
cgv::media::axis_aligned_box::get_min_pnt
const fpnt_type & get_min_pnt() const
return a const reference to corner 0
Definition: axis_aligned_box.h:60
cgv::media::mesh::streaming_mesh< X >::new_triangle
void new_triangle(unsigned int vi, unsigned int vj, unsigned int vk)
construct a new triangle by calling the new polygon method of the callback handler
Definition: streaming_mesh.h:89
cgv::media::mesh::streaming_mesh
class used to perform the marching cubes algorithm
Definition: streaming_mesh.h:25
cgv::utils::progression
Definition: progression.h:14
cgv::media::mesh::slice_info
Definition: marching_cubes.h:24
cgv::media::mesh::streaming_mesh< X >::get_nr_vertices
unsigned int get_nr_vertices() const
return the number of vertices
Definition: streaming_mesh.h:53
cgv::math::mfunc::evaluate
virtual T evaluate(const pnt_type &p) const =0
interface for evaluation of the multivariate function
cgv::media::mesh::streaming_mesh< X >::drop_vertices
void drop_vertices(unsigned int n)
drop n vertices from the front of the deque
Definition: streaming_mesh.h:67
cgv
the cgv namespace
Definition: vr_calib.cxx:9
cgv::media::mesh::marching_cubes_base::pnt_type
cgv::math::fvec< X, 3 > pnt_type
points must have three components
Definition: marching_cubes.h:74