A.I.M.S algorithms


geodesic_mesh.h
Go to the documentation of this file.
1 //Copyright (C) 2008 Danil Kirsanov, MIT License
2 #ifndef GEODESIC_MESH_20071231
3 #define GEODESIC_MESH_20071231
4 
5 #include <cstddef>
6 #include <vector>
7 #include <cmath>
8 #include <iostream>
9 #include <algorithm>
10 #include <fstream>
11 
12 #include "geodesic_mesh_elements.h"
13 #include "geodesic_memory.h"
15 
16 namespace geodesic{
17 
19 {
20  unsigned source;
22 };
23 
24 class Mesh
25 {
26 public:
27  Mesh()
28  {};
29 
30  ~Mesh(){};
31 
32  template<class Points, class Faces>
33  void initialize_mesh_data(unsigned num_vertices,
34  Points& p,
35  unsigned num_faces,
36  Faces& tri, const float* curvature, int mode,int strain); //build mesh from regular point-triangle representation
37 
38  template<class Points, class Faces>
39  void initialize_mesh_data(Points& p, Faces& tri, const float* curvature, int mode,int strain); //build mesh from regular point-triangle representation
40 
41  void update_weight(const float *curvature,int mode,double strain, double sigmo); //build internal structure of the mesh
42 
43  std::vector<Vertex>& vertices(){return m_vertices;};
44  std::vector<Edge>& edges(){return m_edges;};
45  std::vector<Face>& faces(){return m_faces;};
46 
47  unsigned closest_vertices(SurfacePoint* p,
48  std::vector<vertex_pointer>* storage = NULL); //list vertices closest to the point
49 
50 private:
51 
52  void build_adjacencies(const float *curvature,int mode,int strain); //build internal structure of the mesh
53  bool verify(); //verifies connectivity of the mesh and prints some debug info
54 
55  typedef void* void_pointer;
56  void_pointer allocate_pointers(unsigned n)
57  {
58  return m_pointer_allocator.allocate(n);
59  }
60 
61  std::vector<Vertex> m_vertices;
62  std::vector<Edge> m_edges;
63  std::vector<Face> m_faces;
64 
65  SimlpeMemoryAllocator<void_pointer> m_pointer_allocator; //fast memory allocating for Face/Vertex/Edge cross-references
66 };
67 
69  std::vector<vertex_pointer>* storage)
70 {
71  assert(p->type() != UNDEFINED_POINT);
72 
73  if(p->type() == VERTEX)
74  {
75  if(storage)
76  {
77  storage->push_back(static_cast<vertex_pointer>(p->base_element()));
78  }
79  return 1;
80  }
81  else if(p->type() == FACE)
82  {
83  if(storage)
84  {
86  storage->push_back(*vp);
87  storage->push_back(*(vp+1));
88  storage->push_back(*(vp+2));
89  }
90  return 2;
91  }
92  else if(p->type() == EDGE) //for edge include all 4 adjacent vertices
93  {
94  edge_pointer edge = static_cast<edge_pointer>(p->base_element());
95 
96  if(storage)
97  {
98  storage->push_back(edge->adjacent_vertices()[0]);
99  storage->push_back(edge->adjacent_vertices()[1]);
100 
101  for(unsigned i = 0; i < edge->adjacent_faces().size(); ++i)
102  {
103  face_pointer face = edge->adjacent_faces()[i];
104  storage->push_back(face->opposite_vertex(edge));
105  }
106  }
107  return 2 + edge->adjacent_faces().size();
108  }
109 
110  assert(0);
111  return 0;
112 }
113 
114 template<class Points, class Faces>
115 void Mesh::initialize_mesh_data(Points& p, Faces& tri, const float* curvature, int mode, int strain) //build mesh from regular point-triangle representation
116 {
117  assert(p.size() % 3 == 0);
118  unsigned const num_vertices = p.size() / 3;
119  assert(tri.size() % 3 == 0);
120  unsigned const num_faces = tri.size() / 3;
121 
122  initialize_mesh_data(num_vertices, p, num_faces, tri, curvature, mode, strain);
123 }
124 
125 template<class Points, class Faces>
126 void Mesh::initialize_mesh_data(unsigned num_vertices,
127  Points& p,
128  unsigned num_faces,
129  Faces& tri, const float* curvature, int mode, int strain)
130 {
131  //printf("initialize_mesh_data\n");
132 
133  unsigned const approximate_number_of_internal_pointers = (num_vertices + num_faces)*4;
134 
135  //printf("approximate_number_of_internal_pointers = \n");
136 
137  unsigned const max_number_of_pointer_blocks = 100;
138  m_pointer_allocator.reset(approximate_number_of_internal_pointers,
139  max_number_of_pointer_blocks);
140 
141  //printf("copy coordinates to vertices\n");
142 
143  m_vertices.resize(num_vertices);
144  for(unsigned i=0; i<num_vertices; ++i) //copy coordinates to vertices
145  {
146  Vertex& v = m_vertices[i];
147  v.id() = i;
148 
149  unsigned shift = 3*i;
150  v.x() = p[shift];
151  v.y() = p[shift + 1];
152  v.z() = p[shift + 2];
153  }
154 
155  //printf("copy adjacent vertices to polygons/faces \n");
156 
157  m_faces.resize(num_faces);
158  for(unsigned i=0; i<num_faces; ++i) //copy adjacent vertices to polygons/faces
159  {
160  //std::cout << i << std::endl;
161  Face& f = m_faces[i];
162  f.id() = i;
163  f.adjacent_vertices().set_allocation(allocate_pointers(3),3); //allocate three units of memory
164 
165  unsigned shift = 3*i;
166  unsigned vertex_index;
167  for(unsigned j=0; j<3; ++j)
168  {
169  //std::cout << shift + j;
170  vertex_index = tri[shift + j];
171  assert(vertex_index < num_vertices);
172  //if (vertex_index < num_vertices)
173  f.adjacent_vertices()[j] = &m_vertices[vertex_index];
174  }
175  }
176 
177  //printf("build the structure of the mesh \n");
178 
179  build_adjacencies(curvature,mode,strain); //build the structure of the mesh
180 }
181 
182 inline void Mesh::build_adjacencies(const float *curvature, int mode, int strain)
183 {
184  // Vertex->adjacent Faces
185  std::vector<unsigned> count(m_vertices.size()); //count adjacent vertices
186  for(unsigned i=0; i<m_faces.size(); ++i)
187  {
188  Face& f = m_faces[i];
189  for(unsigned j=0; j<3; ++j)
190  {
191  unsigned vertex_id = f.adjacent_vertices()[j]->id();
192 // assert(vertex_id < m_vertices.size());
193  count[vertex_id]++;
194  }
195  }
196 
197  for(unsigned i=0; i<m_vertices.size(); ++i) //reserve space
198  {
199  Vertex& v = m_vertices[i];
200  unsigned num_adjacent_faces = count[i];
201 
202  v.adjacent_faces().set_allocation(allocate_pointers(num_adjacent_faces), //allocate three units of memory
203  num_adjacent_faces);
204  }
205 
206  std::fill(count.begin(), count.end(), 0);
207  for(unsigned i=0; i<m_faces.size(); ++i)
208  {
209  Face& f = m_faces[i];
210  for(unsigned j=0; j<3; ++j)
211  {
212  vertex_pointer v = f.adjacent_vertices()[j];
213  v->adjacent_faces()[count[v->id()]++] = &f;
214  }
215  }
216 
217  //find all edges
218  //i.e. find all half-edges, sort and combine them into edges
219  std::vector<HalfEdge> half_edges(m_faces.size()*3);
220  unsigned k = 0;
221  for(unsigned i=0; i<m_faces.size(); ++i)
222  {
223  Face& f = m_faces[i];
224  for(unsigned j=0; j<3; ++j)
225  {
226  half_edges[k].face_id = i;
227  unsigned vertex_id_1 = f.adjacent_vertices()[j]->id();
228  unsigned vertex_id_2 = f.adjacent_vertices()[(j+1) % 3]->id();
229  half_edges[k].vertex_0 = std::min(vertex_id_1, vertex_id_2);
230  half_edges[k].vertex_1 = std::max(vertex_id_1, vertex_id_2);
231 
232  k++;
233  }
234  }
235  std::sort(half_edges.begin(), half_edges.end());
236 
237  unsigned number_of_edges = 1;
238  for(unsigned i=1; i<half_edges.size(); ++i)
239  {
240  if(half_edges[i] != half_edges[i-1])
241  {
242  ++number_of_edges;
243  }
244  else
245  {
246  if(i<half_edges.size()-1) //sanity check: there should be at most two equal half-edges
247  { //if it fails, most likely the input data are messed up
248 // assert(half_edges[i] != half_edges[i+1]);
249  }
250  }
251  }
252 
253  // Edges->adjacent Vertices and Faces
254  m_edges.resize(number_of_edges);
255  unsigned edge_id = 0;
256  for(unsigned i=0; i<half_edges.size();)
257  {
258  Edge& e = m_edges[edge_id];
259  e.id() = edge_id++;
260 
261  e.adjacent_vertices().set_allocation(allocate_pointers(2),2); //allocate two units of memory
262 
263  e.adjacent_vertices()[0] = &m_vertices[half_edges[i].vertex_0];
264  e.adjacent_vertices()[1] = &m_vertices[half_edges[i].vertex_1];
265 
266  //ARN
267  vertex_pointer v1 = e.adjacent_vertices()[0];
268  vertex_pointer v2 = e.adjacent_vertices()[1];
269 
270  //cout << i << "v1 = " << v1->id() << " curv = " << curvature[v1->id()];;
271  //cout << " v2 = " << v2->id() << " curv = " << curvature[v2->id()] << endl;
272 
273 // if (curvature!= NULL)
274 // e.length() = (v1->distance(v2)*(fabs(curvature[v2->id()] - curvature[v1->id()])) ) ;
275 // else
276 // e.length() = e.adjacent_vertices()[0]->distance(e.adjacent_vertices()[1]);
277 //
278  double a1,a2;
279 
280  if (curvature!= NULL && mode != 0)
281  {
282  // Sulcal
283  if (mode == 1)
284  {
285  a1 = pow ((1.0)/(1.0 + exp(-strain*curvature[v1->id()])), 2);
286  a2 = pow ((1.0)/(1.0 + exp(-strain*curvature[v2->id()])), 2);
287  }
288 
289  // Gyral
290  if (mode == 2)
291  {
292  a1 = pow ((1.0)/(1.0 + exp(strain*curvature[v1->id()])), 2);
293  a2 = pow ((1.0)/(1.0 + exp(strain*curvature[v2->id()])), 2);
294  }
295 
296  e.length() = (v1->distance(v2) * (a1 + a2));
297  }
298  else
299  e.length() = e.adjacent_vertices()[0]->distance(e.adjacent_vertices()[1]);
300 
301  //assert(e.length() > 1e-100);
302  //algorithm works well with non-degenerate meshes only
303 
304  if(i != half_edges.size()-1 && half_edges[i] == half_edges[i+1]) //double edge
305  {
306  e.adjacent_faces().set_allocation(allocate_pointers(2),2);
307  e.adjacent_faces()[0] = &m_faces[half_edges[i].face_id];
308  e.adjacent_faces()[1] = &m_faces[half_edges[i+1].face_id];
309  i += 2;
310  }
311  else //single edge
312  {
313  e.adjacent_faces().set_allocation(allocate_pointers(1),1); //one adjucent faces
314  e.adjacent_faces()[0] = &m_faces[half_edges[i].face_id];
315  i += 1;
316  }
317  }
318 
319  // Vertices->adjacent Edges
320  std::fill(count.begin(), count.end(), 0);
321  for(unsigned i=0; i<m_edges.size(); ++i)
322  {
323  Edge& e = m_edges[i];
324  assert(e.adjacent_vertices().size()==2);
325  count[e.adjacent_vertices()[0]->id()]++;
326  count[e.adjacent_vertices()[1]->id()]++;
327  }
328  for(unsigned i=0; i<m_vertices.size(); ++i)
329  {
330  m_vertices[i].adjacent_edges().set_allocation(allocate_pointers(count[i]),
331  count[i]);
332  }
333  std::fill(count.begin(), count.end(), 0);
334  for(unsigned i=0; i<m_edges.size(); ++i)
335  {
336  Edge& e = m_edges[i];
337  for(unsigned j=0; j<2; ++j)
338  {
339  vertex_pointer v = e.adjacent_vertices()[j];
340  v->adjacent_edges()[count[v->id()]++] = &e;
341  }
342  }
343 
344  // Faces->adjacent Edges
345  for(unsigned i=0; i<m_faces.size(); ++i)
346  {
347  m_faces[i].adjacent_edges().set_allocation(allocate_pointers(3),3);
348  }
349 
350  count.resize(m_faces.size());
351  std::fill(count.begin(), count.end(), 0);
352  for(unsigned i=0; i<m_edges.size(); ++i)
353  {
354  Edge& e = m_edges[i];
355  for(unsigned j=0; j<e.adjacent_faces().size(); ++j)
356  {
357  face_pointer f = e.adjacent_faces()[j];
358  assert(count[f->id()]<3);
359  f->adjacent_edges()[count[f->id()]++] = &e;
360  }
361  }
362 
363  //compute angles for the faces
364  for(unsigned i=0; i<m_faces.size(); ++i)
365  {
366  Face& f = m_faces[i];
367  double abc[3];
368  double sum = 0;
369  for(unsigned j=0; j<3; ++j) //compute angle adjacent to the vertex j
370  {
371  for(unsigned k=0; k<3; ++k)
372  {
373  vertex_pointer v = f.adjacent_vertices()[(j + k)%3];
374  abc[k] = f.opposite_edge(v)->length();
375  }
376 
377  double angle = angle_from_edges(abc[0], abc[1], abc[2]);
378 
379  //ARN
380  //assert(angle>1e-5); //algorithm works well with non-degenerate meshes only
381 
382  f.corner_angles()[j] = angle;
383  sum += angle;
384  }
385  //assert(std::abs(sum - M_PI) < 1e-5); //algorithm works well with non-degenerate meshes only
386  }
387 
388  //define m_turn_around_flag for vertices
389  std::vector<double> total_vertex_angle(m_vertices.size());
390  for(unsigned i=0; i<m_faces.size(); ++i)
391  {
392  Face& f = m_faces[i];
393  for(unsigned j=0; j<3; ++j)
394  {
396  total_vertex_angle[v->id()] += f.corner_angles()[j];
397  }
398  }
399 
400  for(unsigned i=0; i<m_vertices.size(); ++i)
401  {
402  Vertex& v = m_vertices[i];
403  v.saddle_or_boundary() = (total_vertex_angle[v.id()] > 2.0*M_PI - 1e-5);
404  }
405 
406  for(unsigned i=0; i<m_edges.size(); ++i)
407  {
408  Edge& e = m_edges[i];
409  if(e.is_boundary())
410  {
411  e.adjacent_vertices()[0]->saddle_or_boundary() = true;
412  e.adjacent_vertices()[1]->saddle_or_boundary() = true;
413  }
414  }
415 
416  //assert(verify());
417 }
418 
419 inline void Mesh::update_weight(const float *curvature, int mode,double strain,double sigmo)
420 {
421  // Vertex->adjacent Faces
422  std::vector<unsigned> count(m_vertices.size()); //count adjacent vertices
423  for(unsigned i=0; i<m_faces.size(); ++i)
424  {
425  Face& f = m_faces[i];
426  for(unsigned j=0; j<3; ++j)
427  {
428  unsigned vertex_id = f.adjacent_vertices()[j]->id();
429  assert(vertex_id < m_vertices.size());
430  count[vertex_id]++;
431  }
432  }
433 
434  for(unsigned i=0; i<m_vertices.size(); ++i) //reserve space
435  {
436  Vertex& v = m_vertices[i];
437  unsigned num_adjacent_faces = count[i];
438 
439  v.adjacent_faces().set_allocation(allocate_pointers(num_adjacent_faces), //allocate three units of memory
440  num_adjacent_faces);
441  }
442 
443  std::fill(count.begin(), count.end(), 0);
444  for(unsigned i=0; i<m_faces.size(); ++i)
445  {
446  Face& f = m_faces[i];
447  for(unsigned j=0; j<3; ++j)
448  {
450  v->adjacent_faces()[count[v->id()]++] = &f;
451  }
452  }
453 
454  //find all edges
455  //i.e. find all half-edges, sort and combine them into edges
456  std::vector<HalfEdge> half_edges(m_faces.size()*3);
457  unsigned k = 0;
458  for(unsigned i=0; i<m_faces.size(); ++i)
459  {
460  Face& f = m_faces[i];
461  for(unsigned j=0; j<3; ++j)
462  {
463  half_edges[k].face_id = i;
464  unsigned vertex_id_1 = f.adjacent_vertices()[j]->id();
465  unsigned vertex_id_2 = f.adjacent_vertices()[(j+1) % 3]->id();
466  half_edges[k].vertex_0 = std::min(vertex_id_1, vertex_id_2);
467  half_edges[k].vertex_1 = std::max(vertex_id_1, vertex_id_2);
468 
469  k++;
470  }
471  }
472  std::sort(half_edges.begin(), half_edges.end());
473 
474  unsigned number_of_edges = 1;
475  for(unsigned i=1; i<half_edges.size(); ++i)
476  {
477  if(half_edges[i] != half_edges[i-1])
478  {
479  ++number_of_edges;
480  }
481  else
482  {
483  if(i<half_edges.size()-1) //sanity check: there should be at most two equal half-edges
484  { //if it fails, most likely the input data are messed up
485  assert(half_edges[i] != half_edges[i+1]);
486  }
487  }
488  }
489 
490  // Edges->adjacent Vertices and Faces
491  m_edges.resize(number_of_edges);
492  unsigned edge_id = 0;
493  for(unsigned i=0; i<half_edges.size();)
494  {
495  Edge& e = m_edges[edge_id];
496  e.id() = edge_id++;
497 
498  e.adjacent_vertices().set_allocation(allocate_pointers(2),2); //allocate two units of memory
499 
500  e.adjacent_vertices()[0] = &m_vertices[half_edges[i].vertex_0];
501  e.adjacent_vertices()[1] = &m_vertices[half_edges[i].vertex_1];
502 
503  //ARN
504  vertex_pointer v1 = e.adjacent_vertices()[0];
505  vertex_pointer v2 = e.adjacent_vertices()[1];
506 
507  //cout << i << "v1 = " << v1->id() << " curv = " << curvature[v1->id()];;
508  //cout << " v2 = " << v2->id() << " curv = " << curvature[v2->id()] << endl;
509 
510 // if (curvature!= NULL)
511 // e.length() = (v1->distance(v2)*(fabs(curvature[v2->id()] - curvature[v1->id()])) ) ;
512 // else
513 // e.length() = e.adjacent_vertices()[0]->distance(e.adjacent_vertices()[1]);
514 //
515  double a1,a2;
516 
517  if (curvature!= NULL && mode != 0)
518  {
519  // Sulcal
520  if (mode == 1)
521  {
522  a1 = pow ((double)(1.0)/(1.0 + exp(-strain*curvature[v1->id()])), (double)sigmo);
523  a2 = pow ((double)(1.0)/(1.0 + exp(-strain*curvature[v2->id()])), (double)sigmo);
524  }
525 
526  // Gyral
527  if (mode == 2)
528  {
529  a1 = pow ((double)(1.0)/(1.0 + exp(strain*curvature[v1->id()])), (double)sigmo);
530  a2 = pow ((double)(1.0)/(1.0 + exp(strain*curvature[v2->id()])), (double)sigmo);
531  }
532 
533  e.length() = (v1->distance(v2) * (a1 + a2));
534  }
535  else
536  e.length() = e.adjacent_vertices()[0]->distance(e.adjacent_vertices()[1]);
537 
538  //assert(e.length() > 1e-100); //algorithm works well with non-degenerate meshes only
539 
540  if(i != half_edges.size()-1 && half_edges[i] == half_edges[i+1]) //double edge
541  {
542  e.adjacent_faces().set_allocation(allocate_pointers(2),2);
543  e.adjacent_faces()[0] = &m_faces[half_edges[i].face_id];
544  e.adjacent_faces()[1] = &m_faces[half_edges[i+1].face_id];
545  i += 2;
546  }
547  else //single edge
548  {
549  e.adjacent_faces().set_allocation(allocate_pointers(1),1); //one adjucent faces
550  e.adjacent_faces()[0] = &m_faces[half_edges[i].face_id];
551  i += 1;
552  }
553  }
554 
555  // Vertices->adjacent Edges
556  std::fill(count.begin(), count.end(), 0);
557  for(unsigned i=0; i<m_edges.size(); ++i)
558  {
559  Edge& e = m_edges[i];
560  assert(e.adjacent_vertices().size()==2);
561  count[e.adjacent_vertices()[0]->id()]++;
562  count[e.adjacent_vertices()[1]->id()]++;
563  }
564  for(unsigned i=0; i<m_vertices.size(); ++i)
565  {
566  m_vertices[i].adjacent_edges().set_allocation(allocate_pointers(count[i]),
567  count[i]);
568  }
569  std::fill(count.begin(), count.end(), 0);
570  for(unsigned i=0; i<m_edges.size(); ++i)
571  {
572  Edge& e = m_edges[i];
573  for(unsigned j=0; j<2; ++j)
574  {
576  v->adjacent_edges()[count[v->id()]++] = &e;
577  }
578  }
579 
580  // Faces->adjacent Edges
581  for(unsigned i=0; i<m_faces.size(); ++i)
582  {
583  m_faces[i].adjacent_edges().set_allocation(allocate_pointers(3),3);
584  }
585 
586  count.resize(m_faces.size());
587  std::fill(count.begin(), count.end(), 0);
588  for(unsigned i=0; i<m_edges.size(); ++i)
589  {
590  Edge& e = m_edges[i];
591  for(unsigned j=0; j<e.adjacent_faces().size(); ++j)
592  {
593  face_pointer f = e.adjacent_faces()[j];
594  assert(count[f->id()]<3);
595  f->adjacent_edges()[count[f->id()]++] = &e;
596  }
597  }
598 
599  //compute angles for the faces
600  for(unsigned i=0; i<m_faces.size(); ++i)
601  {
602  Face& f = m_faces[i];
603  double abc[3];
604  double sum = 0;
605  for(unsigned j=0; j<3; ++j) //compute angle adjacent to the vertex j
606  {
607  for(unsigned k=0; k<3; ++k)
608  {
609  vertex_pointer v = f.adjacent_vertices()[(j + k)%3];
610  abc[k] = f.opposite_edge(v)->length();
611  }
612 
613  double angle = angle_from_edges(abc[0], abc[1], abc[2]);
614 
615  //ARN
616  //assert(angle>1e-5); //algorithm works well with non-degenerate meshes only
617 
618  f.corner_angles()[j] = angle;
619  sum += angle;
620  }
621  //assert(std::abs(sum - M_PI) < 1e-5); //algorithm works well with non-degenerate meshes only
622  }
623 
624  //define m_turn_around_flag for vertices
625  std::vector<double> total_vertex_angle(m_vertices.size());
626  for(unsigned i=0; i<m_faces.size(); ++i)
627  {
628  Face& f = m_faces[i];
629  for(unsigned j=0; j<3; ++j)
630  {
632  total_vertex_angle[v->id()] += f.corner_angles()[j];
633  }
634  }
635 
636  for(unsigned i=0; i<m_vertices.size(); ++i)
637  {
638  Vertex& v = m_vertices[i];
639  v.saddle_or_boundary() = (total_vertex_angle[v.id()] > 2.0*M_PI - 1e-5);
640  }
641 
642  for(unsigned i=0; i<m_edges.size(); ++i)
643  {
644  Edge& e = m_edges[i];
645  if(e.is_boundary())
646  {
647  e.adjacent_vertices()[0]->saddle_or_boundary() = true;
648  e.adjacent_vertices()[1]->saddle_or_boundary() = true;
649  }
650  }
651 
652  //assert(verify());
653  //assert(verify());
654 }
655 inline bool Mesh::verify() //verifies connectivity of the mesh and prints some debug info
656 {
657  std::cout << std::endl;
658  // make sure that all vertices are mentioned at least once.
659  // though the loose vertex is not a bug, it most likely indicates that something is wrong with the mesh
660  std::vector<bool> map(m_vertices.size(), false);
661  for(unsigned i=0; i<m_edges.size(); ++i)
662  {
663  edge_pointer e = &m_edges[i];
664  map[e->adjacent_vertices()[0]->id()] = true;
665  map[e->adjacent_vertices()[1]->id()] = true;
666  }
667  assert(std::find(map.begin(), map.end(), false) == map.end());
668 
669  //make sure that the mesh is connected trough its edges
670  //if mesh has more than one connected component, it is most likely a bug
671  std::vector<face_pointer> stack(1,&m_faces[0]);
672  stack.reserve(m_faces.size());
673 
674  map.resize(m_faces.size());
675  std::fill(map.begin(), map.end(), false);
676  map[0] = true;
677 
678  while(!stack.empty())
679  {
680  face_pointer f = stack.back();
681  stack.pop_back();
682 
683  for(unsigned i=0; i<3; ++i)
684  {
685  edge_pointer e = f->adjacent_edges()[i];
686  face_pointer f_adjacent = e->opposite_face(f);
687  if(f_adjacent && !map[f_adjacent->id()])
688  {
689  map[f_adjacent->id()] = true;
690  stack.push_back(f_adjacent);
691  }
692  }
693  }
694  assert(std::find(map.begin(), map.end(), false) == map.end());
695 
696  //print some mesh statistics that can be useful in debugging
697  std::cout << "mesh has " << m_vertices.size()
698  << " vertices, " << m_faces.size()
699  << " faces, " << m_edges.size()
700  << " edges\n";
701 
702  unsigned total_boundary_edges = 0;
703  double longest_edge = 0;
704  double shortest_edge = 1e100;
705  for(unsigned i=0; i<m_edges.size(); ++i)
706  {
707  Edge& e = m_edges[i];
708  total_boundary_edges += e.is_boundary() ? 1 : 0;
709  longest_edge = std::max(longest_edge, e.length());
710  shortest_edge = std::min(shortest_edge, e.length());
711  }
712  std::cout << total_boundary_edges << " edges are boundary edges\n";
713  std::cout << "shortest/longest edges are "
714  << shortest_edge << "/"
715  << longest_edge << " = "
716  << shortest_edge/longest_edge
717  << std::endl;
718 
719  double minx = 1e100;
720  double maxx = -1e100;
721  double miny = 1e100;
722  double maxy = -1e100;
723  double minz = 1e100;
724  double maxz = -1e100;
725  for(unsigned i=0; i<m_vertices.size(); ++i)
726  {
727  Vertex& v = m_vertices[i];
728  minx = std::min(minx, v.x());
729  maxx = std::max(maxx, v.x());
730  miny = std::min(miny, v.y());
731  maxy = std::max(maxy, v.y());
732  minz = std::min(minz, v.z());
733  maxz = std::max(maxz, v.z());
734  }
735  std::cout << "enclosing XYZ box:"
736  <<" X[" << minx << "," << maxx << "]"
737  <<" Y[" << miny << "," << maxy << "]"
738  <<" Z[" << minz << "," << maxz << "]"
739  << std::endl;
740 
741  double dx = maxx - minx;
742  double dy = maxy - miny;
743  double dz = maxz - minz;
744  std::cout << "approximate diameter of the mesh is "
745  << sqrt(dx*dx + dy*dy + dz*dz)
746  << std::endl;
747 
748  double min_angle = 1e100;
749  double max_angle = -1e100;
750  for(unsigned i=0; i<m_faces.size(); ++i)
751  {
752  Face& f = m_faces[i];
753  for(unsigned j=0; j<3; ++j)
754  {
755  double angle = f.corner_angles()[j];
756  min_angle = std::min(min_angle, angle);
757  max_angle = std::max(max_angle, angle);
758  }
759  }
760  std::cout << "min/max face angles are "
761  << min_angle/M_PI*180.0 << "/"
762  << max_angle/M_PI*180.0
763  << " degrees\n";
764 
765  std::cout << std::endl;
766  return true;
767 }
768 
770  double* data,
771  Mesh* mesh)
772 {
773  point->set(data);
774  unsigned type = (unsigned) data[3];
775  unsigned id = (unsigned) data[4];
776 
777 
778  if(type == 0) //vertex
779  {
780  point->base_element() = &mesh->vertices()[id];
781  }
782  else if(type == 1) //edge
783  {
784  point->base_element() = &mesh->edges()[id];
785  }
786  else //face
787  {
788  point->base_element() = &mesh->faces()[id];
789  }
790 }
791 
793  double* data,
794  long mesh_id)
795 {
796  data[0] = point->x();
797  data[1] = point->y();
798  data[2] = point->z();
799  data[4] = point->base_element()->id();
800 
801  if(point->type() == VERTEX) //vertex
802  {
803  data[3] = 0;
804  }
805  else if(point->type() == EDGE) //edge
806  {
807  data[3] = 1;
808  }
809  else //face
810  {
811  data[3] = 2;
812  }
813 }
814 
815 } //geodesic
816 
817 #endif
float min(float x, float y)
Definition: thickness.h:105
void reset(unsigned block_size, unsigned max_number_of_blocks)
float max(float x, float y)
Definition: thickness.h:97
pointer allocate(unsigned const n)
double angle_from_edges(double const a, double const b, double const c)
vertex_pointer_vector & adjacent_vertices()
edge_pointer_vector & adjacent_edges()
void set(double new_x, double new_y, double new_z)
void set_allocation(DataPointer begin, unsigned size)
void initialize_mesh_data(unsigned num_vertices, Points &p, unsigned num_faces, Faces &tri, const float *curvature, int mode, int strain)
std::vector< Edge > & edges()
Definition: geodesic_mesh.h:44
Vertex * vertex_pointer
void update_weight(const float *curvature, int mode, double strain, double sigmo)
void fill_surface_point_double(geodesic::SurfacePoint *point, double *data, long mesh_id)
std::vector< Vertex > & vertices()
Definition: geodesic_mesh.h:43
std::vector< Face > & faces()
Definition: geodesic_mesh.h:45
vertex_pointer opposite_vertex(edge_pointer e)
edge_pointer opposite_edge(vertex_pointer v)
face_pointer_vector & adjacent_faces()
double distance(double *v)
void fill_surface_point_structure(geodesic::SurfacePoint *point, double *data, Mesh *mesh)
unsigned closest_vertices(SurfacePoint *p, std::vector< vertex_pointer > *storage=NULL)
Definition: geodesic_mesh.h:68