OpenMEEG
geometry_reader.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 #include <vector>
43 #include <interface.h>
44 #include <GeometryExceptions.H>
45 #include <IOUtils.H>
46 #include <PropertiesSpecialized.h>
47 #include <geometry_io.h>
48 
49 namespace OpenMEEG {
50 
53  public:
55 
57  void read_geom(const std::string&);
58 
60  void read_cond(const std::string&);
61 
62  private:
64 
66  bool is_relative_path(const std::string& name);
67  #if WIN32
68  static const char PathSeparator[];
69  #else
70  static const char PathSeparator = '/';
71  #endif
72  };
73  #if WIN32
74  const char GeometryReader::PathSeparator[] = "/\\";
75  #endif
76 
77  bool GeometryReader::is_relative_path(const std::string& name) {
78  #if WIN32
79  const char c0 = name[0];
80  if ( c0 == '/' || c0 == '\\' ) {
81  return false;
82  }
83  const char c1 = name[1];
84  const char c2 = name[2];
85  return !( std::isalpha(c0) && c1 == ':' && ( c2 == '/' || c2 == '\\' ) );
86  #else
87  return ( name[0] != PathSeparator );
88  #endif
89  }
90 
91  void GeometryReader::read_geom(const std::string& geometry) {
92  // Read the head file description and load the information into the data structures.
93 
94  // check file data/README.rst for complete details about the .geom format.
95  // The syntax of the head description is a header ("# Domain Description 1.1") followed
96  // by two or three sections:
97  //
98  // - The first section is optional (for backward compatibility).
99  // Starting with the keyword "MeshFile", it follows the path to the VTK/vtp file containing the meshes.
100  // OR
101  // Starting with the keyword "Meshes", it follows the number of meshes, and the paths to the meshes (with the keyword "Mesh").
102  //
103  // - the second section is introduced by the keyword "Interfaces" followed by a number
104  // (the number of interfaces) and optionnally the keyword "Mesh" (for backward compatibility).
105  // The section contains the list of interfaces preceded by keyword "Interface".
106  //
107  // - the third section is introduced by the keyword "Domains" and the number of domains
108  // (everything else on the line containing the keyword is currently ignored). The section
109  // contains domains descriptions, one per line. Each domain consist of:
110  //
111  // o a domain name.
112  // o a list of IDs (signed numbers or signed names): the sign ('+'by default) of the ID depicts
113  // whether the interior or the exterior of the interface should be used to select the domain.
114  //
115  // Any line starting with # is considered a comment and is silently ignored.
116 
117  std::ifstream ifs(geometry.c_str());
118 
119  if ( !ifs.is_open() ) {
120  throw OpenMEEG::OpenError(geometry);
121  }
122 
123  // Get the version of the geometry file format.
124 
125  unsigned version[2];
126  ifs >> io_utils::match("# Domain Description ") >> version[0] >> io_utils::match(".") >> version[1];
127 
128  if ( ifs.fail() ) {
129  throw OpenMEEG::WrongFileFormat(geometry);
130  }
131 
133  if (version[0]==1) {
134  if (version[1]==0) {
136  std::cerr << "(DEPRECATED) Please consider updating your geometry file to the new format 1.1 (see data/README.rst): "
137  << geometry << std::endl;
138  }
139  if (version[1]==1)
141  }
142 
144  std::cerr << "Domain Description version not available !" << std::endl;
145  throw OpenMEEG::WrongFileFormat(geometry);
146  }
147 
148  // extract the absolut path of geometry file
149  const std::string::size_type pos = geometry.find_last_of(this->PathSeparator);
150  const std::string path = (pos == std::string::npos) ? "" : geometry.substr(0, pos+1);
151 
152  // Process meshes. -----------------------------------------------------------------------------------
154  // Read the mesh section of the description file.
155  // Try to load the meshfile (VTK::vtp file)
156  // or try to load the meshes
157  bool Is_MeshFile, Is_Meshes;
158  ifs >> io_utils::skip_comments("#") >> io_utils::match_optional("MeshFile", Is_MeshFile);
159  ifs >> io_utils::skip_comments("#") >> io_utils::match_optional("Meshes", Is_Meshes);
160  if ( Is_MeshFile ) {
161  std::string name;
162  ifs >> io_utils::skip_comments("#") >> io_utils::filename(name, '"', false);
163  const std::string& full_name = (is_relative_path(name))?path+name:name;
164  geom.load_vtp(full_name);
165  } else if ( Is_Meshes ) {
166  unsigned nb_meshes;
167  ifs >> nb_meshes;
168  std::vector<std::string> meshname(nb_meshes); // names
169  std::vector<std::string> filename(nb_meshes);
170  std::vector<std::string> fullname(nb_meshes);
171  Meshes meshes(nb_meshes);
172  for ( unsigned i = 0; i < nb_meshes; ++i ) {
173  bool unnamed;
174  ifs >> io_utils::skip_comments("#") >> io_utils::match_optional("Mesh:", unnamed);
175  if ( unnamed ) {
176  ifs >> io_utils::filename(filename[i], '"', false);
177  std::stringstream defaultname;
178  defaultname << i+1;
179  meshname[i] = defaultname.str();
180  } else {
181  ifs >> io_utils::match("Mesh")
182  >> io_utils::token(meshname[i], ':')
183  >> io_utils::filename(filename[i], '"', false);
184  }
185  fullname[i] = (is_relative_path(filename[i]))?path+filename[i]:filename[i];
186  // Load the mesh.
187  meshes[i].load(fullname[i], false);
188  meshes[i].name() = meshname[i];
189  }
190  // Now properly load the meshes into the geometry (not dupplicated vertices)
191  geom.import_meshes(meshes);
192  }
193  }
194 
195  // Process interfaces. -----------------------------------------------------------------------------------
196  unsigned nb_interfaces;
197  bool trash; // backward compatibility, catch "Mesh" optionnally.
198  ifs >> io_utils::skip_comments('#')
199  >> io_utils::match("Interfaces") >> nb_interfaces >> io_utils::match_optional("Mesh", trash);
200 
201  if ( ifs.fail() ) {
202  throw OpenMEEG::WrongFileFormat(geometry);
203  }
204 
205  // load the interfaces
206  std::string id; // id of mesh/interface/domain
207  Interfaces interfaces;
208  // if meshes are not already loaded
209  if ( geom.nb_meshes() == 0 ) { // ---------------------------------------
210  geom.meshes_.reserve(nb_interfaces);
211  std::vector<std::string> interfacename(nb_interfaces);
212  std::vector<std::string> filename(nb_interfaces);
213  std::vector<std::string> fullname(nb_interfaces);
214 
215  // First read the total number of vertices
216 
217  unsigned nb_vertices = 0;
218  for ( unsigned i = 0; i < nb_interfaces; ++i ) {
219  bool unnamed;
220  ifs >> io_utils::skip_comments("#") >> io_utils::match_optional("Interface:", unnamed);
221  if ( unnamed ) {
222  ifs >> io_utils::filename(filename[i], '"', false);
223  std::stringstream defaultname;
224  defaultname << i+1;
225  interfacename[i] = defaultname.str();
226  } else if (geom.version_id==Geometry::VERSION10) { // backward compatibility
227  std::stringstream defaultname;
228  defaultname << i+1;
229  interfacename[i] = defaultname.str();
230  ifs >> io_utils::filename(filename[i], '"', false);
231  } else {
232  ifs >> io_utils::match("Interface")
233  >> io_utils::token(interfacename[i], ':')
234  >> io_utils::filename(filename[i], '"', false);
235  }
236  Mesh m;
237  fullname[i] = (is_relative_path(filename[i]))?path+filename[i]:filename[i];
238  nb_vertices += m.load(fullname[i], false, false);
239  }
240  geom.vertices_.reserve(nb_vertices);
241  // Second really load the meshes
242  for ( unsigned i = 0; i < nb_interfaces; ++i ) {
243  geom.meshes_.push_back(Mesh(geom.vertices_, interfacename[i]));
244  geom.meshes_[i].load(fullname[i], false);
245  interfaces.push_back( Interface(interfacename[i]) );
246  interfaces[i].push_back(OrientedMesh(geom.meshes_[i], true)); // one mesh per interface, (well oriented)
247  }
248  } else { // -----------------------
249  std::string interfacename;
250  for ( unsigned i = 0; i < nb_interfaces; ++i ) {
251  bool unnamed;
252  std::string line; // extract a line and parse it
253  ifs >> io_utils::skip_comments("#");
254  std::getline(ifs, line);
255  std::istringstream iss(line);
256  iss >> io_utils::match_optional("Interface:", unnamed);
257  if ( unnamed ) {
258  std::stringstream defaultname;
259  defaultname << i+1;
260  interfacename = defaultname.str();
261  } else {
262  iss >> io_utils::match("Interface")
263  >> io_utils::token(interfacename, ':');
264  }
265  interfaces.push_back( interfacename );
266  while ( iss >> id ) {
267  bool oriented = true; // does the id starts with a '-' or a '+' ?
268  if ( ( id[0] == '-' ) || ( id[0] == '+' ) ) {
269  oriented = ( id[0] == '+' );
270  id = id.substr(1, id.size());
271  }
272  interfaces[i].push_back(OrientedMesh(geom.mesh(id), oriented));
273  }
274  }
275  }
276 
277  // Process domains. -----------------------------------------------------------------------------------
278 
279  unsigned num_domains;
280  ifs >> io_utils::skip_comments('#') >> io_utils::match("Domains") >> num_domains;
281 
282  if ( ifs.fail() ) {
283  throw OpenMEEG::WrongFileFormat(geometry);
284  }
285 
286  geom.domains_.resize(num_domains);
287  for ( Domains::iterator dit = geom.domain_begin(); dit != geom.domain_end(); ++dit) {
288  std::string line;
289  if (geom.version_id==Geometry::VERSION10) { // backward compatibility
290  ifs >> io_utils::skip_comments('#') >> io_utils::match("Domain") >> dit->name();
291  } else {
292  ifs >> io_utils::skip_comments('#') >> io_utils::match("Domain") >> io_utils::token(dit->name(), ':');
293  }
294  getline(ifs, line);
295  std::istringstream iss(line);
296  while ( iss >> id ) {
297  bool found = false;
298  bool inside = false; // does the id starts with a '-' or a '+' ?
299  if ( ( id[0] == '-' ) || ( id[0] == '+' ) ) {
300  inside = ( id[0] == '-' );
301  id = id.substr(1, id.size());
302  } else if ( id == "shared" ) {
303  std::cerr << "(DEPRECATED) Keyword shared is useless. Please consider updating your geometry file to the new format 1.1 (see data/README.rst): " << geometry << std::endl;
304  break;
305  }
306  for ( Interfaces::iterator iit = interfaces.begin(); iit != interfaces.end() ; ++iit) {
307  if ( iit->name() == id ) {
308  found = true;
309  if ( !iit->check() ) { // check and correct global orientation
310  std::cerr << "Interface \"" << iit->name() << "\" is not closed !" << std::endl;
311  std::cerr << "Please correct a mesh orientation when defining the interface in the geometry file." << std::endl;
312  exit(1);
313  }
314  dit->push_back(HalfSpace(*iit, inside));
315  }
316  }
317  if ( !found ) {
318  throw OpenMEEG::NonExistingDomain<std::string>(dit->name(), id);
319  }
320  }
321  }
322 
323  // Search for the outermost domain and set boolean OUTERMOST on the domain in the vector domains.
324  // An outermost domain is (here) defined as the only domain outside represented by only one interface.
325 
326  Domains::iterator dit_out;
327  for ( Domains::iterator dit = geom.domain_begin(); dit != geom.domain_end(); ++dit) {
328  bool outer = true;
329  for ( Domain::iterator hit = dit->begin(); hit != dit->end(); ++hit)
330  outer = outer && !(hit->inside());
331 
332  if (outer) {
333  dit_out = dit;
334  dit_out->outermost() = true;
335  for (Domain::iterator hit = dit_out->begin(); hit != dit_out->end(); ++hit) {
336  hit->interface().set_to_outermost();
337  }
338  break;
339  }
340  }
341 
342  // Determine if the geometry is nested or not
343  // The geometry is considered non nested if (at least) one domain is defined as being outside two or more interfaces
344  // OR
345  // if 2 interfaces are composed by a same mesh oriented once correctly once wrongly.
346 
347  bool nested = true;
348  for ( Domains::const_iterator dit = geom.domain_begin(); dit != geom.domain_end(); ++dit) {
349  unsigned out_interface = 0;
350  if ( dit != dit_out ) {
351  for ( Domain::const_iterator hit = dit->begin(); hit != dit->end(); ++hit) {
352  if ( !hit->inside() ) {
353  out_interface++;
354  }
355  }
356  }
357  if ( out_interface >= 2 ) {
358  nested = false;
359  break;
360  }
361  }
362 
363  if ( nested ) {
364  for ( Geometry::const_iterator mit = geom.begin(); mit != geom.end(); ++mit) {
365  unsigned m_oriented = 0;
366  for ( Domains::const_iterator dit = geom.domain_begin(); dit != geom.domain_end(); ++dit) {
367  for ( Domain::const_iterator hit = dit->begin(); hit != dit->end(); ++hit) {
368  for ( Interface::const_iterator iit = hit->first.begin(); iit != hit->first.end(); ++iit) {
369  if ( iit->mesh() == *mit ) {
370  m_oriented += (iit->orientation());
371  }
372  }
373  }
374  }
375  if ( m_oriented == 0 ) {
376  nested = false; // TODO unless a mesh is defined but unused ...
377  break;
378  }
379  }
380  }
381  geom.is_nested_ = nested;
382 
383  if ( ifs.fail() ) {
384  throw OpenMEEG::WrongFileFormat(geometry);
385  }
386 
387  // Close the input file. -----------------------------------------------------------------------------------
388  ifs.close();
389  }
390 
391  void GeometryReader::read_cond(const std::string& condFileName) {
392 
393  typedef Utils::Properties::Named<std::string, Conductivity<double> > HeadProperties;
394  HeadProperties properties(condFileName.c_str());
395 
396  // Store the internal conductivity of the external boundary of domain i
397  // and store the external conductivity of the internal boundary of domain i
398  for ( Domains::iterator dit = geom.domain_begin(); dit != geom.domain_end(); ++dit) {
399  try {
400  const Conductivity<double>& cond = properties.find(dit->name());
401  dit->sigma() = cond.sigma();
402  } catch( const Utils::Properties::UnknownProperty<HeadProperties::Id>& e) {
403  throw OpenMEEG::BadDomain(dit->name());
404  }
405  }
406  }
407 }
A class to read geometry and cond file.
static const char PathSeparator
bool is_relative_path(const std::string &name)
void read_geom(const std::string &)
read a geometry file
void read_cond(const std::string &)
read a cond file
Geometry contains the electrophysiological model Here are stored the vertices, meshes and domains.
Definition: geometry.h:61
iterator end()
Definition: geometry.h:78
void import_meshes(const Meshes &m)
imports meshes from a list of meshes
Domains::iterator domain_begin()
Definition: geometry.h:84
iterator begin()
Iterators.
Definition: geometry.h:76
const Mesh & mesh(const std::string &id) const
void load_vtp(const std::string &filename)
Definition: geometry.h:126
Domains::iterator domain_end()
Definition: geometry.h:86
Vertices vertices_
Definition: geometry.h:146
size_t nb_meshes() const
Definition: geometry.h:106
Meshes::const_iterator const_iterator
Definition: geometry.h:73
Domains domains_
Definition: geometry.h:148
VersionId version_id
Members.
Definition: geometry.h:145
a HalfSpace is a pair of Interface and boolean A simple domain (HalfSpace) is given by an interface (...
Definition: domain.h:58
Interface class An interface is a closed-shape composed of oriented meshes (here pointer to meshes)
Definition: interface.h:72
Mesh class.
Definition: mesh.h:85
unsigned load(const std::string &filename, const bool &verbose=true, const bool &read_all=true)
Read mesh from file.
An Oriented Mesh is a mesh associated with a boolean stating if it is well oriented.
Definition: interface.h:54
std::vector< Mesh > Meshes
A vector of Mesh is called Meshes.
Definition: mesh.h:309
std::vector< Interface > Interfaces
A vector of Interface is called Interfaces.
Definition: interface.h:127