OpenMEEG
mesh.h
Go to the documentation of this file.
1 /*
2 Project Name : OpenMEEG
3 
4 © INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
5 GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
6 Emmanuel OLIVI
7 Maureen.Clerc.AT.inria.fr, keriven.AT.certis.enpc.fr,
8 kybic.AT.fel.cvut.cz, papadop.AT.inria.fr)
9 
10 The OpenMEEG software is a C++ package for solving the forward/inverse
11 problems of electroencephalography and magnetoencephalography.
12 
13 This software is governed by the CeCILL-B license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL-B
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's authors, the holders of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL-B license and that you accept its terms.
38 */
39 
40 #pragma once
41 
42 // for IO:s
43 #include <iostream>
44 #include <fstream>
45 
46 #include <vector>
47 #include <set>
48 #include <map>
49 #include <stack>
50 #include <string>
51 
52 #include <om_common.h>
53 #include <triangle.h>
54 #include <IOUtils.H>
55 #include <om_utils.h>
56 #include <sparse_matrix.h>
57 
58 #ifdef USE_VTK
59 #include <vtkPolyData.h>
60 #include <vtkPoints.h>
61 #include <vtkPolyDataReader.h>
62 #include <vtkXMLPolyDataReader.h>
63 #include <vtkDataReader.h>
64 #include <vtkCellArray.h>
65 #include <vtkCharArray.h>
66 #include <vtkPointData.h>
67 #include <vtkDataArray.h>
68 #endif
69 
70 #ifdef USE_GIFTI
71 extern "C" {
72  #include <gifti_io.h>
73 }
74 #endif
75 
76 namespace OpenMEEG {
77 
78  enum Filetype { VTK, TRI, BND, MESH, OFF, GIFTI };
79 
85  class OPENMEEG_EXPORT Mesh: public Triangles {
86  public:
87 
88  typedef std::vector<Triangle*> VectPTriangle;
89  typedef std::vector<Vertex*> VectPVertex;
90  typedef VectPVertex::iterator vertex_iterator;
91  typedef VectPVertex::const_iterator const_vertex_iterator;
92  typedef VectPVertex::const_reverse_iterator const_vertex_reverse_iterator;
93 
94  // Constructors:
96 
97  Mesh(): Triangles(), name_(""), all_vertices_(0), outermost_(false), allocate_(false), current_barrier_(false), isolated_(false) { }
98 
102 
103  Mesh(const unsigned& nv,const unsigned& nt): name_(""), outermost_(false), allocate_(true), current_barrier_(false), isolated_(false) {
104  all_vertices_ = new Vertices;
105  all_vertices_->reserve(nv); // allocates space for the vertices
106  reserve(nt);
107  }
108 
110 
111  Mesh(const Mesh& m): Triangles(), current_barrier_(false), isolated_(false) { *this = m; }
112 
114 
115  Mesh(Vertices& av,const std::string name = ""): name_(name), all_vertices_(&av), outermost_(false), allocate_(false), current_barrier_(false), isolated_(false) {
116  set_vertices_.insert(all_vertices_->begin(), all_vertices_->end());
117  }
118 
120 
121  Mesh(std::string filename,const bool verbose=true,const std::string n=""): name_(n), outermost_(false), allocate_(true), current_barrier_(false), isolated_(false) {
122  unsigned nb_v = load(filename, false, false);
123  all_vertices_ = new Vertices(nb_v); // allocates space for the vertices
124  load(filename, verbose);
125  }
126 
128 
129  ~Mesh() { destroy(); }
130 
131  // Iterators on vertices
132 
133  vertex_iterator vertex_begin() { return vertices_.begin(); }
134  vertex_iterator vertex_end() { return vertices_.end(); }
135 
136  const size_t vertex_size() const { return vertices_.size(); } // Just for old OpenMP implementations.
137 
138  const_vertex_iterator vertex_begin() const { return vertices_.begin(); }
139  const_vertex_iterator vertex_end() const { return vertices_.end(); }
140 
141  const_vertex_reverse_iterator vertex_rbegin() const { return vertices_.rbegin(); }
142  const_vertex_reverse_iterator vertex_rend() const { return vertices_.rend(); }
143 
144  std::string & name() { return name_; }
145  const std::string & name() const { return name_; }
146 
147  const VectPVertex & vertices() const { return vertices_; }
148  const size_t nb_vertices() const { return vertices_.size(); }
149  const size_t nb_triangles() const { return size(); }
150 
151  Vertices all_vertices() const { return *all_vertices_; }
152  const size_t nb_all_vertices() const { return all_vertices_->size(); }
153 
155 
156  void add_vertex(const Vertex& v);
157 
161 
162  void info(const bool verbous = false) const;
163  bool has_self_intersection() const;
164  bool intersection(const Mesh&) const;
165  bool has_correct_orientation() const;
166  void build_mesh_vertices();
167  void generate_indices();
168  void update();
169  void merge(const Mesh&, const Mesh&);
170 
172 
173  void flip_triangles() {
174  for (iterator tit = begin(); tit != end(); ++tit)
175  tit->flip();
176  }
177 
178  void correct_local_orientation();
179  void correct_global_orientation();
180  double compute_solid_angle(const Vect3& p) const;
181  const VectPTriangle& get_triangles_for_vertex(const Vertex& V) const;
182  VectPTriangle adjacent_triangles(const Triangle&) const;
183  Normal normal(const Vertex& v) const;
184  void laplacian(SymMatrix &A) const;
185 
186  bool& outermost() { return outermost_; }
187  const bool& outermost() const { return outermost_; }
188 
193 
194  void smooth(const double& smoothing_intensity, const unsigned& niter);
195 
197 
198  void gradient_norm2(SymMatrix &A) const;
199 
200  // for IO:s --------------------------------------------------------------------
205 
206  unsigned load(const std::string& filename,const bool& verbose=true,const bool& read_all=true);
207  unsigned load_tri(std::istream& , const bool& read_all = true);
208  unsigned load_tri(const std::string&, const bool& read_all = true);
209  unsigned load_bnd(std::istream& , const bool& read_all = true);
210  unsigned load_bnd(const std::string&, const bool& read_all = true);
211  unsigned load_off(std::istream& , const bool& read_all = true);
212  unsigned load_off(const std::string&, const bool& read_all = true);
213  unsigned load_mesh(std::istream& , const bool& read_all = true);
214  unsigned load_mesh(const std::string&, const bool& read_all = true);
215 
216  #ifdef USE_VTK
217  unsigned load_vtk(std::istream& , const bool& read_all = true);
218  unsigned load_vtk(const std::string&, const bool& read_all = true);
219  unsigned get_data_from_vtk_reader(vtkPolyDataReader* vtkMesh, const bool& read_all);
220  #else
221  template <typename T>
222  unsigned load_vtk(T, const bool& read_all = true) {
223  std::cerr << "You have to compile OpenMEEG with VTK to read VTK/VTP files. (specify USE_VTK to cmake)" << std::endl;
224  exit(1);
225  }
226  #endif
227 
228  #ifdef USE_GIFTI
229  unsigned load_gifti(const std::string&, const bool& read_all = true);
230  void save_gifti(const std::string&) const;
231  #else
232  template <typename T>
233  unsigned load_gifti(T, const bool&) {
234  std::cerr << "You have to compile OpenMEEG with GIFTI to read GIFTI files" << std::endl;
235  exit(1);
236  }
237  template <typename T>
238  void save_gifti(T) const {
239  std::cerr << "You have to compile OpenMEEG with GIFTI to read GIFTI files" << std::endl;
240  exit(1);
241  }
242  #endif
243 
244 
247 
248  void save(const std::string& filename) const ;
249  void save_vtk(const std::string&) const;
250  void save_bnd(const std::string&) const;
251  void save_tri(const std::string&) const;
252  void save_off(const std::string&) const;
253  void save_mesh(const std::string&) const;
254 
255  // IO:s ----------------------------------------------------------------------------
256 
257  Mesh& operator=(const Mesh& m) {
258  if ( this != &m )
259  copy(m);
260  return *this;
261  }
262 
263  friend std::istream& operator>>(std::istream& is, Mesh& m);
264 
265  private:
266 
268 
269  typedef std::map<std::pair<const Vertex *, const Vertex *>, int> EdgeMap;
270 
271  void destroy();
272  void copy(const Mesh&);
273 
274  // regarding mesh orientation
275 
276  const EdgeMap compute_edge_map() const;
277  void orient_adjacent_triangles(std::stack<Triangle*>& t_stack,std::map<Triangle*,bool>& tri_reoriented);
278  bool triangle_intersection(const Triangle&,const Triangle&) const;
279 
281 
282  Vect3 P1gradient(const Vect3& p0,const Vect3& p1,const Vect3& p2) const { return p1^p2/(p0*(p1^p2)); }
283 
285 
286  double P0gradient_norm2(const Triangle& t1,const Triangle& t2) const {
287  return std::pow(t1.normal()*t2.normal(),2)/(t1.center()-t2.center()).norm2();
288  }
289 
290  std::string name_;
291  std::map<const Vertex*,VectPTriangle> links_;
293  VectPVertex vertices_;
294  bool outermost_;
295  bool allocate_;
296  std::set<Vertex> set_vertices_;
299  bool isolated_;
300  public:
301  const bool& current_barrier() const { return current_barrier_; }
302  bool& current_barrier() { return current_barrier_; }
303  const bool& isolated() const { return isolated_; }
304  bool& isolated() { return isolated_; }
305  };
306 
308 
309  typedef std::vector<Mesh> Meshes;
310 }
Normal & normal()
Definition: triangle.h:117
const VectPVertex & vertices() const
Definition: mesh.h:147
const_vertex_reverse_iterator vertex_rend() const
Definition: mesh.h:142
std::vector< Vertex > Vertices
Definition: vertex.h:73
const size_t nb_triangles() const
Definition: mesh.h:149
Vertices * all_vertices_
Pointer to all the vertices.
Definition: mesh.h:292
std::string name_
Name of the mesh.
Definition: mesh.h:290
bool current_barrier_
handle multiple 0 conductivity domains
Definition: mesh.h:298
const std::string & name() const
Definition: mesh.h:145
std::istream & operator>>(std::istream &is, Conductivity< REP > &m)
std::set< Vertex > set_vertices_
Definition: mesh.h:296
vertex_iterator vertex_begin()
Definition: mesh.h:133
std::vector< Triangle > Triangles
Definition: triangle.h:176
VectPVertex::const_iterator const_vertex_iterator
Definition: mesh.h:91
Mesh(const unsigned &nv, const unsigned &nt)
constructor from scratch (add vertices/triangles one by one)
Definition: mesh.h:103
const_vertex_iterator vertex_begin() const
Definition: mesh.h:138
Vertices all_vertices() const
Definition: mesh.h:151
const bool & isolated() const
Definition: mesh.h:303
bool & current_barrier()
Definition: mesh.h:302
VectPVertex::const_reverse_iterator const_vertex_reverse_iterator
Definition: mesh.h:92
Mesh(const Mesh &m)
constructor from another mesh
Definition: mesh.h:111
Mesh class.
Definition: mesh.h:85
std::vector< Mesh > Meshes
A vector of Mesh is called Meshes.
Definition: mesh.h:309
std::vector< Vertex * > VectPVertex
Definition: mesh.h:89
std::string & name()
Definition: mesh.h:144
vertex_iterator vertex_end()
Definition: mesh.h:134
Mesh()
default constructor
Definition: mesh.h:97
bool & outermost()
Definition: mesh.h:186
double P0gradient_norm2(const Triangle &t1, const Triangle &t2) const
P0gradient_norm2 : aux function to compute the square norm of the surfacic gradient.
Definition: mesh.h:286
std::map< const Vertex *, VectPTriangle > links_
links[&v] are the triangles that contain vertex v.
Definition: mesh.h:291
Triangle.
Definition: triangle.h:54
Vertex.
Definition: vertex.h:52
bool allocate_
Are the vertices allocate within the mesh or shared ?
Definition: mesh.h:295
unsigned load_vtk(T, const bool &read_all=true)
Definition: mesh.h:222
const size_t nb_vertices() const
Definition: mesh.h:148
const size_t vertex_size() const
Definition: mesh.h:136
bool & isolated()
Definition: mesh.h:304
bool isolated_
Definition: mesh.h:299
const_vertex_iterator vertex_end() const
Definition: mesh.h:139
std::map< std::pair< const Vertex *, const Vertex * >, int > EdgeMap
map the edges with an unsigned
Definition: mesh.h:269
const_vertex_reverse_iterator vertex_rbegin() const
Definition: mesh.h:141
VectPVertex::iterator vertex_iterator
Definition: mesh.h:90
Vect3 P1gradient(const Vect3 &p0, const Vect3 &p1, const Vect3 &p2) const
P1gradient : aux function to compute the surfacic gradient.
Definition: mesh.h:282
Mesh(std::string filename, const bool verbose=true, const std::string n="")
constructor loading directly a mesh file named
Definition: mesh.h:121
Mesh(Vertices &av, const std::string name="")
constructor using an outisde storage for vertices
Definition: mesh.h:115
Vect3.
Definition: vect3.h:62
VectPVertex vertices_
Vector of pointers to the mesh vertices.
Definition: mesh.h:293
unsigned load_gifti(T, const bool &)
Definition: mesh.h:233
Vect3 center() const
Definition: triangle.h:153
Mesh & operator=(const Mesh &m)
Definition: mesh.h:257
void save_gifti(T) const
Definition: mesh.h:238
const bool & outermost() const
Returns True if it is an outermost mesh.
Definition: mesh.h:187
~Mesh()
Destructor.
Definition: mesh.h:129
bool outermost_
Is it an outermost mesh ? (i.e does it touch the Air domain)
Definition: mesh.h:294
Filetype
Definition: mesh.h:78
std::vector< Triangle * > VectPTriangle
Definition: mesh.h:88
const size_t nb_all_vertices() const
Definition: mesh.h:152
void flip_triangles()
Flip all triangles.
Definition: mesh.h:173
const bool & current_barrier() const
Definition: mesh.h:301