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  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  size_t nb_vertices() const { return vertices_.size(); }
149  size_t nb_triangles() const { return size(); }
150 
151  Vertices all_vertices() const { return *all_vertices_; }
152  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;
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 
180  double compute_solid_angle(const Vect3& p) 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_;
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 }
Mesh class.
Definition: mesh.h:85
unsigned load_mesh(const std::string &, const bool &read_all=true)
void update()
recompute triangles normals, area, and links
std::vector< Triangle * > VectPTriangle
Definition: mesh.h:88
bool current_barrier_
handle multiple 0 conductivity domains
Definition: mesh.h:298
void laplacian(SymMatrix &A) const
compute mesh laplacian
unsigned load_vtk(T, const bool &read_all=true)
Definition: mesh.h:222
const_vertex_reverse_iterator vertex_rend() const
Definition: mesh.h:142
bool has_correct_orientation() const
check the local orientation of the mesh triangles
std::string & name()
Definition: mesh.h:144
const bool & outermost() const
Returns True if it is an outermost mesh.
Definition: mesh.h:187
void correct_local_orientation()
correct the local orientation of the mesh triangles
Mesh(const unsigned &nv, const unsigned &nt)
constructor from scratch (add vertices/triangles one by one)
Definition: mesh.h:103
VectPVertex vertices_
Vector of pointers to the mesh vertices.
Definition: mesh.h:293
void save_off(const std::string &) const
VectPVertex::iterator vertex_iterator
Definition: mesh.h:90
bool allocate_
Are the vertices allocate within the mesh or shared ?
Definition: mesh.h:295
unsigned load_tri(const std::string &, const bool &read_all=true)
size_t vertex_size() const
Definition: mesh.h:136
bool & isolated()
Definition: mesh.h:304
void save_bnd(const std::string &) const
size_t nb_vertices() const
Definition: mesh.h:148
unsigned load_bnd(std::istream &, const bool &read_all=true)
size_t nb_all_vertices() const
Definition: mesh.h:152
unsigned load(const std::string &filename, const bool &verbose=true, const bool &read_all=true)
Read mesh from file.
Vect3 P1gradient(const Vect3 &p0, const Vect3 &p1, const Vect3 &p2) const
P1gradient : aux function to compute the surfacic gradient.
Definition: mesh.h:282
bool intersection(const Mesh &) const
check if the mesh intersects another mesh
unsigned load_off(const std::string &, const bool &read_all=true)
VectPVertex::const_iterator const_vertex_iterator
Definition: mesh.h:91
void smooth(const double &smoothing_intensity, const unsigned &niter)
Smooth Mesh.
bool outermost_
Is it an outermost mesh ? (i.e does it touch the Air domain)
Definition: mesh.h:294
unsigned load_mesh(std::istream &, const bool &read_all=true)
bool triangle_intersection(const Triangle &, const Triangle &) const
size_t nb_triangles() const
Definition: mesh.h:149
Normal normal(const Vertex &v) const
get the Normal at vertex
unsigned load_off(std::istream &, const bool &read_all=true)
void correct_global_orientation()
correct the global orientation (if there is one)
void merge(const Mesh &, const Mesh &)
properly merge two meshes into one
void save_gifti(T) const
Definition: mesh.h:238
friend std::istream & operator>>(std::istream &is, Mesh &m)
insert a triangle into the mesh
vertex_iterator vertex_end()
Definition: mesh.h:134
void save_mesh(const std::string &) const
bool isolated_
Definition: mesh.h:299
std::map< const Vertex *, VectPTriangle > links_
links[&v] are the triangles that contain vertex v.
Definition: mesh.h:291
~Mesh()
Destructor.
Definition: mesh.h:129
bool has_self_intersection() const
check if the mesh self-intersects
Mesh(std::string filename, const bool verbose=true, const std::string n="")
constructor loading directly a mesh file named
Definition: mesh.h:121
void generate_indices()
generate indices (if allocate)
const std::string & name() const
Definition: mesh.h:145
Mesh(const Mesh &m)
constructor from another mesh
Definition: mesh.h:111
unsigned load_tri(std::istream &, const bool &read_all=true)
const EdgeMap compute_edge_map() const
const bool & current_barrier() const
Definition: mesh.h:301
double compute_solid_angle(const Vect3 &p) const
Given a point p, it computes the solid angle.
bool & outermost()
Definition: mesh.h:186
void save(const std::string &filename) const
Save mesh to file.
unsigned load_gifti(T, const bool &)
Definition: mesh.h:233
void build_mesh_vertices()
construct the list of the mesh vertices out of its triangles
void add_vertex(const Vertex &v)
properly add vertex to the list.
const VectPVertex & vertices() const
Definition: mesh.h:147
const_vertex_reverse_iterator vertex_rbegin() const
Definition: mesh.h:141
void info(const bool verbous=false) const
Print info Print to std::cout some info about the mesh.
std::set< Vertex > set_vertices_
Definition: mesh.h:296
const VectPTriangle & get_triangles_for_vertex(const Vertex &V) const
get the triangles associated with vertex V
std::string name_
Name of the mesh.
Definition: mesh.h:290
VectPVertex::const_reverse_iterator const_vertex_reverse_iterator
Definition: mesh.h:92
bool & current_barrier()
Definition: mesh.h:302
const bool & isolated() const
Definition: mesh.h:303
vertex_iterator vertex_begin()
Definition: mesh.h:133
Mesh()
default constructor
Definition: mesh.h:97
unsigned load_bnd(const std::string &, const bool &read_all=true)
void copy(const Mesh &)
std::vector< Vertex * > VectPVertex
Definition: mesh.h:89
VectPTriangle adjacent_triangles(const Triangle &) const
get the adjacent triangles
void save_tri(const std::string &) const
void gradient_norm2(SymMatrix &A) const
Compute the square norm of the surfacic gradient.
Mesh(Vertices &av, const std::string name="")
constructor using an outisde storage for vertices
Definition: mesh.h:115
Vertices all_vertices() const
Definition: mesh.h:151
const_vertex_iterator vertex_begin() const
Definition: mesh.h:138
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
void orient_adjacent_triangles(std::stack< Triangle * > &t_stack, std::map< Triangle *, bool > &tri_reoriented)
std::map< std::pair< const Vertex *, const Vertex * >, int > EdgeMap
map the edges with an unsigned
Definition: mesh.h:269
void save_vtk(const std::string &) const
const_vertex_iterator vertex_end() const
Definition: mesh.h:139
Mesh & operator=(const Mesh &m)
Definition: mesh.h:257
void flip_triangles()
Flip all triangles.
Definition: mesh.h:173
Vertices * all_vertices_
Pointer to all the vertices.
Definition: mesh.h:292
Triangle.
Definition: triangle.h:54
Vect3 center() const
Definition: triangle.h:153
Normal & normal()
Definition: triangle.h:117
Vect3.
Definition: vect3.h:62
Vertex.
Definition: vertex.h:52
std::vector< Vertex > Vertices
Definition: vertex.h:73
std::vector< Triangle > Triangles
Definition: triangle.h:176
std::vector< Mesh > Meshes
A vector of Mesh is called Meshes.
Definition: mesh.h:309
Filetype
Definition: mesh.h:78
@ BND
Definition: mesh.h:78
@ TRI
Definition: mesh.h:78
@ OFF
Definition: mesh.h:78
@ GIFTI
Definition: mesh.h:78
@ MESH
Definition: mesh.h:78
@ VTK
Definition: mesh.h:78