cgv
dual_contouring.h
1 #pragma once
2 
3 #include <vector>
4 #include <deque>
5 #include <cgv/utils/progression.h>
6 #include <cgv/math/qem.h>
7 #include <cgv/math/mfunc.h>
8 #include <cgv/media/axis_aligned_box.h>
9 #include "streaming_mesh.h"
10 
11 namespace cgv {
12  namespace media {
13  namespace mesh {
14 
16 template <typename T>
17 struct cell_info {
18  cgv::math::qem<T> Q;
19  cgv::math::fvec<T,3> center;
20  int count;
21 };
22 
23 
25 template <typename T>
26 struct dc_slice_info
27 {
28  unsigned int resx, resy;
29  std::vector<T> values;
30  std::vector<bool> flags;
31  std::vector<int> indices;
32  typedef cgv::math::qem<T> qem_type;
33  typedef cgv::math::fvec<T,3> pnt_type;
34  std::vector<cell_info<T> > cell_infos;
36  dc_slice_info(unsigned int _resx, unsigned int _resy) : resx(_resx), resy(_resy) {
37  unsigned int n = resx*resy;
38  values.resize(n);
39  flags.resize(n);
40  indices.resize(n);
41  cell_infos.resize(n);
42  }
44  void init() {
45  int n = resx*resy;
46  for (int i=0; i<n; ++i) {
47  indices[i] = -1;
48  cell_infos[i].Q = qem_type();
49  cell_infos[i].count = 0;
50  cell_infos[i].center = pnt_type();
51  }
52  }
54  bool flag(int x, int y) const { return flags[y*resx+x]; }
56  const T& value(int x, int y) const { return values[y*resx+x]; }
57  T& value(int x, int y) { return values[y*resx+x]; }
59  void set_value(int x, int y, T value, T iso_value) {
60  int i = y*resx+x;
61  values[i] = value;
62  flags[i] = value>iso_value;
63  }
65  const int& index(int x, int y) const { return indices[y*resx+x]; }
66  int& index(int x, int y) { return indices[y*resx+x]; }
68  const qem_type& get_qem(int x, int y) const { return cell_infos[y*resx+x].Q; }
69  qem_type& ref_qem(int x, int y) { return cell_infos[y*resx+x].Q; }
71  int count(int x, int y) const { return cell_infos[y*resx+x].count; }
72  int& count(int x, int y) { return cell_infos[y*resx+x].count; }
74  const pnt_type& center(int x, int y) const { return cell_infos[y*resx+x].center; }
75  pnt_type& center(int x, int y) { return cell_infos[y*resx+x].center; }
77  const cell_info<T>& info(int x, int y) const { return cell_infos[y*resx+x]; }
78  cell_info<T>& info(int x, int y) { return cell_infos[y*resx+x]; }
79 };
80 
82 template <typename X, typename T>
83 class dual_contouring : public streaming_mesh<X>
84 {
85 public:
86  typedef streaming_mesh<X> base_type;
88  typedef cgv::math::fvec<X,3> pnt_type;
90  typedef cgv::math::fvec<X,3> vec_type;
92  typedef cgv::math::qem<X> qem_type;
94  typedef cell_info<X> cell_info_type;
95 private:
96  pnt_type p, minp;
97  unsigned int resx, resy, resz;
98  vec_type d;
99  T iso_value;
100 protected:
101  const cgv::math::v3_func<X,T>& func;
102  X epsilon;
103  unsigned int max_nr_iters;
104  X consistency_threshold;
105 public:
107  dual_contouring(const cgv::math::v3_func<X,T>& _func,
108  streaming_mesh_callback_handler* _smcbh,
109  const X& _consistency_threshold = 0.01f, unsigned int _max_nr_iters = 10,
110  const X& _epsilon = 1e-6f) :
111  func(_func), max_nr_iters(_max_nr_iters), consistency_threshold(_consistency_threshold), epsilon(_epsilon)
112  {
113  base_type::set_callback_handler(_smcbh);
114  }
116  void compute_cell_vertex(dc_slice_info<T> *info_ptr, int i, int j)
117  {
118  if (info_ptr->count(i,j) == 0) {
119  info_ptr->index(i,j) = -1;
120  return;
121  }
122  pnt_type p_ref = info_ptr->center(i,j) / X(info_ptr->count(i,j));
123  cgv::math::vec<X> min_pnt = info_ptr->get_qem(i, j).minarg(p_ref.to_vec(), X(0.1), d.length());
124  pnt_type q(min_pnt.size(), min_pnt);
125  info_ptr->index(i,j) = this->new_vertex(q);
126  }
128  void generate_quad(unsigned int vi, unsigned int vj, unsigned int vk, unsigned int vl, bool reorient)
129  {
130  if (reorient)
131  base_type::new_quad(vi,vl,vk,vj);
132  else
133  base_type::new_quad(vi,vj,vk,vl);
134  }
136  void compute_edge_point(const T& _v_1, const T& _v_2, int e, const pnt_type& p, const vec_type& d,
137  pnt_type& q, vec_type& n)
138  {
139  X de = d(e);
140  T v_1 = _v_1;
141  T v_2 = _v_2;
142  pnt_type p_end = p;
143  q = p;
144  unsigned int nr_iters = 0;
145  bool finished;
146  X ref = consistency_threshold*fabs(v_2-v_1);
147  do {
148  finished = false;
149  // compute point on edge
150  X alpha;
151  if (fabs(v_2-v_1) > epsilon) {
152  finished = true;
153  alpha = (X)(iso_value-v_1)/(v_2-v_1);
154  }
155  else
156  alpha = (X) 0.5;
157 
158  q(e) = p_end(e) - (1-alpha)*de;
159  // compute normal at point
160  cgv::math::vec<X> nml_vec = func.evaluate_gradient(q.to_vec());
161  n = vec_type(nml_vec.size(), nml_vec);
162  n.normalize();
163 
164  // stop if maximum number of iterations reached
165  if (++nr_iters >= max_nr_iters)
166  break;
167  // check if gradient is in accordance with value difference
168  if (finished) {
169  if (alpha < 0 || alpha > 1) {
170  std::cout << "invalid alpha = " << alpha << std::endl;
171  }
172  X f_e = (v_2 - v_1) / de;
173  if (fabs(f_e - n(e)) > consistency_threshold*f_e) {
174  X v = func.evaluate(q.to_vec());
175  if (fabs(v-iso_value) > ref) {
176  if ((v > iso_value) == (v_2 > iso_value)) {
177  de *= alpha;
178  v_2 = v;
179  p_end(e) = q(e);
180  }
181  else {
182  v_1 = v;
183  de *= 1-alpha;
184  }
185  finished = false;
186  }
187  }
188  }
189  } while (!finished);
190  if (nr_iters > 15)
191  std::cout << "not converged after " << nr_iters << " iterations" << std::endl;
192  // if gradient does not define normal
193  if (n.length() < 1e-6) {
194  // define normal from edge direction
195  n = vec_type(0, 0, 0);
196  n(e) = T((v_1 < v_2) ? 1 : 0);
197  }
198  else
199  n.normalize();
200  }
202  void process_edge_plane(const T& v_1, const T& v_2, int e,
203  cell_info_type* C1, cell_info_type* C2, cell_info_type* C3, cell_info_type* C4)
204  {
205  pnt_type q;
206  vec_type n;
207  compute_edge_point(v_1, v_2, e, p, d, q, n);
208  // construct qem
209  qem_type Q(q.to_vec(),n.to_vec());
210  // add qem to incident cells
211  if (C1) { if (C1->count == 0) { C1->Q = Q; C1->center = q; } else { C1->Q += Q; C1->center += q; } ++C1->count; }
212  if (C2) { if (C2->count == 0) { C2->Q = Q; C2->center = q; } else { C2->Q += Q; C2->center += q; } ++C2->count; }
213  if (C3) { if (C3->count == 0) { C3->Q = Q; C3->center = q; } else { C3->Q += Q; C3->center += q; } ++C3->count; }
214  if (C4) { if (C4->count == 0) { C4->Q = Q; C4->center = q; } else { C4->Q += Q; C4->center += q; } ++C4->count; }
215  }
217  void process_slice(dc_slice_info<T> *prev_info_ptr, dc_slice_info<T> *info_ptr)
218  {
219  unsigned int i,j;
220  info_ptr->init();
221  for (j = 0, p(1) = minp(1); j < resy; ++j, p(1) += d(1))
222  for (i = 0, p(0) = minp(0); i < resx; ++i, p(0)+=d(0)) {
223  // eval function on slice
224  info_ptr->set_value(i,j,func.evaluate(p.to_vec()),iso_value);
225  // process slice internal edges
226  if (i > 0 && info_ptr->flag(i-1,j) != info_ptr->flag(i,j))
227  process_edge_plane(info_ptr->value(i-1,j),
228  info_ptr->value(i,j), 0,
229  &info_ptr->info(i-1,j),
230  j > 0 ? &info_ptr->info(i-1,j-1) : 0,
231  prev_info_ptr ? &prev_info_ptr->info(i-1,j) : 0,
232  (j > 0 && prev_info_ptr) ? &prev_info_ptr->info(i-1,j-1) : 0);
233  if (j > 0 && info_ptr->flag(i,j-1) != info_ptr->flag(i,j))
234  process_edge_plane(info_ptr->value(i,j-1),
235  info_ptr->value(i,j), 1,
236  &info_ptr->info(i,j-1),
237  i > 0 ? &info_ptr->info(i-1,j-1) : 0,
238  prev_info_ptr ? &prev_info_ptr->info(i,j-1) : 0,
239  (i > 0 && prev_info_ptr) ? &prev_info_ptr->info(i-1,j-1) : 0);
240  }
241  }
243  void process_slab(dc_slice_info<T> *info_ptr_1, dc_slice_info<T> *info_ptr_2)
244  {
245  unsigned int i,j;
246  for (j = 0, p(1) = minp(1); j < resy; ++j, p(1) += d(1))
247  for (i = 0, p(0) = minp(0); i < resx; ++i, p(0)+=d(0)) {
248  // process slab edges
249  if (info_ptr_1->flag(i,j) != info_ptr_2->flag(i,j))
250  process_edge_plane(info_ptr_1->value(i,j),
251  info_ptr_2->value(i,j), 2,
252  &info_ptr_1->info(i,j),
253  j > 0 ? &info_ptr_1->info(i,j-1) : 0,
254  i > 0 ? &info_ptr_1->info(i-1,j) : 0,
255  (i > 0 && j > 0) ? &info_ptr_1->info(i-1,j-1) : 0);
256  // compute cell vertices
257  if (i>0 && j>0)
258  compute_cell_vertex(info_ptr_1, i-1,j-1);
259  }
260  // generate the quads of inner edges inside the slab
261  for (j = 1, p(1) = minp(1)+d(1); j < resy-1; ++j, p(1) += d(1))
262  for (i = 1, p(0) = minp(0)+d(0); i < resx-1; ++i, p(0)+=d(0))
263  if (info_ptr_1->flag(i,j) != info_ptr_2->flag(i,j))
264  generate_quad(info_ptr_1->index(i,j),
265  info_ptr_1->index(i-1,j),
266  info_ptr_1->index(i-1,j-1),
267  info_ptr_1->index(i,j-1),
268  info_ptr_1->flag(i,j));
269  }
271  void generate_slice_quads(dc_slice_info<T> *info_ptr_1, dc_slice_info<T> *info_ptr_2)
272  {
273  unsigned int i,j;
274  for (j = 1, p(1) = minp(1)+d(1); j < resy-1; ++j, p(1) += d(1))
275  for (i = 1, p(0) = minp(0)+d(0); i < resx-1; ++i, p(0)+=d(0)) {
276  if (info_ptr_2->flag(i-1,j) != info_ptr_2->flag(i,j))
277  generate_quad(info_ptr_1->index(i-1,j),
278  info_ptr_2->index(i-1,j),
279  info_ptr_2->index(i-1,j-1),
280  info_ptr_1->index(i-1,j-1),
281  info_ptr_2->flag(i-1,j));
282  if (info_ptr_2->flag(i,j-1) != info_ptr_2->flag(i,j))
283  generate_quad(info_ptr_1->index(i,j-1),
284  info_ptr_1->index(i-1,j-1),
285  info_ptr_2->index(i-1,j-1),
286  info_ptr_2->index(i,j-1),
287  info_ptr_2->flag(i,j-1));
288  }
289  }
291  void extract(const T& _iso_value,
292  const axis_aligned_box<X,3>& box,
293  unsigned int _resx, unsigned int _resy, unsigned int _resz,
294  bool show_progress = false)
295  {
296  // prepare private members
297  resx = _resx; resy = _resy; resz = _resz;
298  minp = p = box.get_min_pnt();
299  d = box.get_extent();
300  d(0) /= (resx-1); d(1) /= (resy-1); d(2) /= (resz-1);
301  iso_value = _iso_value;
302 
303  // prepare progression
304  cgv::utils::progression prog;
305  if (show_progress) prog.init("extraction", resz, 10);
306 
307  // construct three slice infos
308  dc_slice_info<T> slice_info_1(resx,resy), slice_info_2(resx,resy), slice_info_3(resx,resy);
309  dc_slice_info<T> *slice_info_ptrs[3] = { &slice_info_1, &slice_info_2, &slice_info_3 };
310 
311  // iterate through all slices
312  unsigned int nr_vertices[4] = { 0, 0, 0, 0 };
313  unsigned int k, n;
314 
315  process_slice(0, slice_info_ptrs[0]);
316  p(2) += d(2);
317  process_slice(slice_info_ptrs[0], slice_info_ptrs[1]);
318  process_slab(slice_info_ptrs[0], slice_info_ptrs[1]);
319  p(2) += d(2);
320  // show progression
321  if (show_progress) {
322  prog.step();
323  prog.step();
324  }
325  for (k=2; k<resz; ++k, p(2) += d(2)) {
326  n = base_type::get_nr_vertices();
327  // evaluate function on next slice and construct slice interior vertices
328  dc_slice_info<T> *info_ptr_0 = slice_info_ptrs[(k-2)%3];
329  dc_slice_info<T> *info_ptr_1 = slice_info_ptrs[(k-1)%3];
330  dc_slice_info<T> *info_ptr_2 = slice_info_ptrs[k%3];
331  process_slice(info_ptr_1, info_ptr_2);
332  process_slab(info_ptr_1, info_ptr_2);
333  generate_slice_quads(info_ptr_0, info_ptr_1);
334 
335  n = base_type::get_nr_vertices()-n;
336  nr_vertices[k%4] = n;
337  n = nr_vertices[(k+3)%4];
338  if (n > 0)
339  base_type::drop_vertices(n);
340  // show progression
341  if (show_progress)
342  prog.step();
343  }
344  }
345 };
346  }
347  }
348 }
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::math::qem::minarg
vec< T > minarg(const vec< T > &p_ref, T relative_epsilon, T max_distance=-1, T epsilon=1e-10) const
Definition: qem.h:98
cgv::math::fvec::length
T length() const
length of the vector L2-Norm
Definition: fvec.h:219
cgv::media::mesh::cell_info
information stored per cell
Definition: dual_contouring.h:17
cgv::media::mesh::dc_slice_info::set_value
void set_value(int x, int y, T value, T iso_value)
set value and flag at once
Definition: dual_contouring.h:59
cgv::media::mesh::dual_contouring::process_slice
void process_slice(dc_slice_info< T > *prev_info_ptr, dc_slice_info< T > *info_ptr)
process a slice
Definition: dual_contouring.h:217
cgv::media::mesh::dual_contouring::cell_info_type
cell_info< X > cell_info_type
qem type must have dimension three
Definition: dual_contouring.h:94
cgv::math::fvec< T, 3 >
cgv::media::mesh::dual_contouring::generate_quad
void generate_quad(unsigned int vi, unsigned int vj, unsigned int vk, unsigned int vl, bool reorient)
construct a quadrilateral
Definition: dual_contouring.h:128
cgv::media::mesh::dual_contouring::vec_type
cgv::math::fvec< X, 3 > vec_type
vectors must have three components
Definition: dual_contouring.h:90
cgv::math::vec
A column vector class.
Definition: fvec.h:13
cgv::math::vec::size
unsigned size() const
number of elements
Definition: vec.h:97
cgv::math::v3_func
Definition: mfunc.h:63
cgv::media::mesh::streaming_mesh_callback_handler
Definition: streaming_mesh.h:13
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::mesh::dual_contouring::pnt_type
cgv::math::fvec< X, 3 > pnt_type
points must have three components
Definition: dual_contouring.h:88
cgv::media::mesh::dual_contouring
class used to perform the marching cubes algorithm
Definition: dual_contouring.h:84
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::axis_aligned_box
Definition: axis_aligned_box.h:15
cgv::media::mesh::dual_contouring::qem_type
cgv::math::qem< X > qem_type
qem type must have dimension three
Definition: dual_contouring.h:92
cgv::math::fvec::to_vec
vec< T > to_vec() const
conversion to vector type
Definition: fvec.h:411
cgv::math::qem
dimension independent implementation of quadric error metrics
Definition: qem.h:14
cgv::media::mesh::dc_slice_info
Definition: dual_contouring.h:27
cgv::media::mesh::dual_contouring::process_edge_plane
void process_edge_plane(const T &v_1, const T &v_2, int e, cell_info_type *C1, cell_info_type *C2, cell_info_type *C3, cell_info_type *C4)
construct plane through edge and add it to the incident qems
Definition: dual_contouring.h:202
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::dual_contouring::extract
void extract(const T &_iso_value, const axis_aligned_box< X, 3 > &box, unsigned int _resx, unsigned int _resy, unsigned int _resz, bool show_progress=false)
extract iso surface and send quads to dual contouring handler
Definition: dual_contouring.h:291
cgv::media::mesh::streaming_mesh
class used to perform the marching cubes algorithm
Definition: streaming_mesh.h:25
cgv::media::mesh::streaming_mesh< X >::new_quad
void new_quad(unsigned int vi, unsigned int vj, unsigned int vk, unsigned int vl)
construct a new quad by calling the new polygon method of the callback handler
Definition: streaming_mesh.h:99
cgv::utils::progression
Definition: progression.h:14
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::math::mfunc::evaluate_gradient
virtual vec_type evaluate_gradient(const pnt_type &p) const
Definition: mfunc.h:31
cgv::media::mesh::dual_contouring::compute_cell_vertex
void compute_cell_vertex(dc_slice_info< T > *info_ptr, int i, int j)
construct a new vertex on an edge
Definition: dual_contouring.h:116
cgv::media::mesh::dual_contouring::dual_contouring
dual_contouring(const cgv::math::v3_func< X, T > &_func, streaming_mesh_callback_handler *_smcbh, const X &_consistency_threshold=0.01f, unsigned int _max_nr_iters=10, const X &_epsilon=1e-6f)
construct dual contouring object
Definition: dual_contouring.h:107