OpenMEEG
geometry_io.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 <map>
43 #include <sstream>
44 
45 #ifdef USE_VTK
46 #include <vtkPolyData.h>
47 #include <vtkPoints.h>
48 #include <vtkXMLPolyDataReader.h>
49 #include <vtkXMLPolyDataWriter.h>
50 #include <vtkCellArray.h>
51 #include <vtkDoubleArray.h>
52 #include <vtkUnsignedIntArray.h>
53 #include <vtkPointData.h>
54 #include <vtkCellData.h>
55 #include <vtkDataArray.h>
56 #include <vtkAbstractArray.h>
57 #include <vtkStringArray.h>
58 #include <vtkSmartPointer.h>
59 #endif
60 
61 namespace OpenMEEG {
62 
64 
65  void Geometry::load_vtp(const std::string& filename, Matrix& data, const bool READ_DATA) {
66  #ifdef USE_VTK
67  vtkSmartPointer<vtkXMLPolyDataReader> reader = vtkSmartPointer<vtkXMLPolyDataReader>::New();
68  reader->SetFileName(filename.c_str()); // Specify file name of vtp data file to read
69  reader->Update();
70 
71  vtkSmartPointer<vtkPolyData> vtkMesh = reader->GetOutput();
72 
73  int trash;
74  vtkSmartPointer<vtkUnsignedIntArray> v_indices = vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetPointData()->GetArray("Indices", trash));
75  vtkSmartPointer<vtkUnsignedIntArray> c_indices = vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetCellData()->GetArray("Indices", trash));
76 
77  const unsigned npts = vtkMesh->GetNumberOfPoints();
78  vertices_.reserve(npts); // alocate memory for the vertices
79 
80  for (unsigned i = 0; i < npts; ++i) {
81  if (trash == -1) { // no indices provided
82  vertices_.push_back(Vertex(vtkMesh->GetPoint(i)[0], vtkMesh->GetPoint(i)[1], vtkMesh->GetPoint(i)[2]));
83  } else {
84  vertices_.push_back(Vertex(vtkMesh->GetPoint(i)[0], vtkMesh->GetPoint(i)[1], vtkMesh->GetPoint(i)[2], v_indices->GetValue(i)));
85  }
86  }
87 
88  // find the number of different meshes reading the string array associated with the cells
89  std::vector<std::string> meshes_name;
90 
91  // ids/mesh name
92  vtkSmartPointer<vtkStringArray> cell_id = vtkStringArray::SafeDownCast(vtkMesh->GetCellData()->GetAbstractArray("Names"));
93 
94  const unsigned ntrgs = vtkMesh->GetNumberOfCells();
95  for (unsigned i = 0; i < ntrgs; ++i) {
96  if (std::find(meshes_name.begin(), meshes_name.end(), cell_id->GetValue(i)) == meshes_name.end()) {
97  meshes_name.push_back(cell_id->GetValue(i));
98  }
99  }
100 
101  // create the meshes
102  for (std::vector<std::string>::const_iterator vit = meshes_name.begin(); vit != meshes_name.end(); ++vit) {
103  meshes_.push_back(Mesh(vertices_, *vit));
104  }
105 
106  // insert the triangle and mesh vertices address into the right mesh
107  vtkSmartPointer<vtkIdList> l;
108 
109  for (iterator mit = begin(); mit != end(); ++mit) {
110  for (unsigned i = 0; i < ntrgs; ++i) {
111  // get the mesh which has this name
112  if (cell_id->GetValue(i) == mit->name()) {
113  if (vtkMesh->GetCellType(i) == VTK_TRIANGLE) {
114  l = vtkMesh->GetCell(i)->GetPointIds();
115  if (trash != -1) { // no indices provided
116  mit->push_back(Triangle(vertices_[l->GetId(0)],
117  vertices_[l->GetId(1)],
118  vertices_[l->GetId(2)],
119  c_indices->GetValue(i)));
120  } else {
121  mit->push_back(Triangle(vertices_[l->GetId(0)],
122  vertices_[l->GetId(1)],
123  vertices_[l->GetId(2)]));
124  }
125  } else {
126  std::cerr << "This is not a triangulation" << std::endl;
127  exit(1);
128  }
129  }
130  }
131 
132  mit->build_mesh_vertices();
133  mit->update();
134  if (READ_DATA) {
135  const unsigned nc = vtkMesh->GetPointData()->GetNumberOfArrays()-1;
136  const unsigned nl = vtkMesh->GetNumberOfPoints()+vtkMesh->GetNumberOfCells();
137  if (nc != 0) {
138  data = Matrix(nl, nc);
139  for (unsigned ic = 0; ic < nc; ++ic) {
140  std::ostringstream sic;
141  sic << ic;
142  for (unsigned iil = 0; iil < vtkMesh->GetPointData()->GetArray(("Potentials-"+sic.str()).c_str(), trash)->GetSize(); ++iil) {
143  data(vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetPointData()->GetArray("Indices", trash))->GetValue(iil), ic) =
144  vtkDoubleArray::SafeDownCast(vtkMesh->GetPointData()->GetArray(("Potentials-"+sic.str()).c_str(), trash))->GetValue(iil);
145  }
146  for (unsigned iil = 0; iil < vtkMesh->GetCellData()->GetArray(("Currents-"+sic.str()).c_str(), trash)->GetSize(); ++iil) {
147  data(vtkUnsignedIntArray::SafeDownCast(vtkMesh->GetCellData()->GetArray("Indices", trash))->GetValue(iil), ic) =
148  vtkDoubleArray::SafeDownCast(vtkMesh->GetCellData()->GetArray(("Currents-"+sic.str()).c_str(), trash))->GetValue(iil);
149  }
150  }
151  }
152  }
153  }
154  #else
155  std::cerr << "Error: please specify USE_VTK to cmake" << std::endl; // TODO in Legacy format ? // Exceptions
156  #endif
157  }
158 
160 
161  void Geometry::write_vtp(const std::string& filename, const Matrix& data) const {
162 
163  #ifdef USE_VTK
164  vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
165  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New(); // vertices
166  vtkSmartPointer<vtkDoubleArray> normals = vtkSmartPointer<vtkDoubleArray>::New(); // normals
167  vtkSmartPointer<vtkStringArray> cell_id = vtkSmartPointer<vtkStringArray>::New(); // ids/mesh name
168  vtkSmartPointer<vtkUnsignedIntArray> cell_indices = vtkSmartPointer<vtkUnsignedIntArray>::New(); // indices
169  vtkSmartPointer<vtkUnsignedIntArray> point_indices = vtkSmartPointer<vtkUnsignedIntArray>::New(); // indices
170  std::vector<vtkSmartPointer<vtkDoubleArray> > potentials(data.ncol()); // potential on vertices
171  std::vector<vtkSmartPointer<vtkDoubleArray> > currents(data.ncol()); // current on triangles
172 
173  normals->SetNumberOfComponents(3); // 3d normals (ie x, y, z)
174  normals->SetName("Normals");
175  cell_id->SetName("Names");
176  cell_indices->SetName("Indices");
177  point_indices->SetName("Indices");
178 
179  if (data.nlin() != 0) {
180  // Check the data corresponds to the geometry
181  bool HAS_OUTERMOST = false; // data has last p values ?
182  if (data.nlin() == size()) {
183  HAS_OUTERMOST = true;
184  } else {
185  om_error(data.nlin() == size()-outermost_interface().nb_triangles());
186  }
187  for (unsigned j = 0; j < data.ncol(); ++j) {
188  std::stringstream sdip;
189  sdip << j;
190  potentials[j] = vtkSmartPointer<vtkDoubleArray>::New();
191  currents[j] = vtkSmartPointer<vtkDoubleArray>::New();
192  potentials[j]->SetName(("Potentials-"+sdip.str()).c_str());
193  currents[j]->SetName(("Currents-"+sdip.str()).c_str());
194 
195  for (Vertices::const_iterator vit = vertex_begin(); vit != vertex_end(); ++vit) {
196  potentials[j]->InsertNextValue(data(vit->index(), j));
197  }
198 
199  for (Meshes::const_iterator mit = meshes_.begin(); mit != meshes_.end(); ++mit) {
200  for (Mesh::const_iterator tit = mit->begin(); tit != mit->end(); ++tit) {
201  currents[j]->InsertNextValue(((mit->outermost() && !HAS_OUTERMOST) ? 0.0 : data(tit->index(), j)));
202  }
203  }
204  }
205  }
206 
207  std::map<const Vertex*, unsigned> map;
208 
209  unsigned i = 0;
210  for (Vertices::const_iterator vit = vertex_begin(); vit != vertex_end(); ++vit, ++i) {
211  points->InsertNextPoint((*vit)(0), (*vit)(1), (*vit)(2));
212  point_indices->InsertNextValue(vit->index());
213  map[&*vit] = i;
214  }
215  // Add the vertices to polydata
216  polydata->SetPoints(points);
217 
218  vtkSmartPointer<vtkCellArray> polys = vtkSmartPointer<vtkCellArray>::New(); // the triangles
219 
220  for (Meshes::const_iterator mit = meshes_.begin(); mit != meshes_.end(); ++mit) {
221  for ( Mesh::const_iterator tit = mit->begin(); tit != mit->end(); ++tit) {
222  vtkIdType triangle[3] = { map[&(tit->s1())], map[&(tit->s2())], map[&(tit->s3())] };
223  polys->InsertNextCell(3, triangle);
224  cell_id->InsertNextValue(mit->name());
225  cell_indices->InsertNextValue(tit->index());
226  }
227  }
228 
229  // Add the triangles to polydata
230  polydata->SetPolys(polys);
231  // Add the array of Names to polydata
232  polydata->GetCellData()->AddArray(cell_id);
233  // Add the array of Indices to polydata cells
234  polydata->GetCellData()->AddArray(cell_indices);
235  // Add the array of Indices to polydata points
236  polydata->GetPointData()->AddArray(point_indices);
237  // Add optional potentials and currents
238  if (data.nlin() != 0) {
239  for (unsigned j = 0; j < data.ncol(); ++j) {
240  polydata->GetPointData()->AddArray(potentials[j]);
241  polydata->GetCellData()->AddArray(currents[j]);
242  }
243  }
244 
245  vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
246  writer->SetFileName(filename.c_str());
247 
248  #if VTK_MAJOR_VERSION <= 5
249  writer->SetInput(polydata);
250  #else
251  writer->SetInputData(polydata);
252  #endif
253  writer->Write();
254 
255  #else
256  std::cerr << "Error: please specify USE_VTK to cmake" << std::endl; // TODO write in Legacy format ? // Exceptions
257  #endif
258  }
259 }
OpenMEEG::Geometry::iterator
Meshes::iterator iterator
Default iterator of a Geometry is an Iterator on the meshes.
Definition: geometry.h:72
OpenMEEG::Triangle
Triangle.
Definition: triangle.h:54
OpenMEEG::Geometry::meshes_
Meshes meshes_
Definition: geometry.h:147
OpenMEEG::Vertex
Vertex.
Definition: vertex.h:52
OpenMEEG::Geometry::size
size_t size() const
the total number of vertices + triangles
Definition: geometry.h:102
OpenMEEG::Mesh
Mesh class.
Definition: mesh.h:85
OpenMEEG::Geometry::vertices_
Vertices vertices_
Definition: geometry.h:146
OpenMEEG::LinOpInfo::nlin
size_t nlin() const
Definition: linop.h:77
OpenMEEG
Definition: analytics.h:45
OpenMEEG::Geometry::begin
iterator begin()
Iterators.
Definition: geometry.h:76
OpenMEEG::Geometry::vertex_end
Vertices::iterator vertex_end()
Definition: geometry.h:82
OpenMEEG::LinOpInfo::ncol
virtual size_t ncol() const
Definition: linop.h:80
OpenMEEG::Matrix
Matrix class.
Definition: matrix.h:61
OpenMEEG::Geometry::load_vtp
void load_vtp(const std::string &filename)
Definition: geometry.h:126
OpenMEEG::Geometry::nb_triangles
size_t nb_triangles() const
Definition: geometry.h:104
OpenMEEG::Geometry::vertex_begin
Vertices::iterator vertex_begin()
Definition: geometry.h:80
OpenMEEG::Geometry::write_vtp
void write_vtp(const std::string &filename, const Matrix &data=Matrix()) const
write a VTK\vtp file
Definition: geometry_io.h:161
OpenMEEG::Geometry::end
iterator end()
Definition: geometry.h:78
OpenMEEG::Geometry::outermost_interface
const Interface & outermost_interface() const
returns the outermost interface (only valid for nested geometries).