MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_bem_surface.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mne_bem_surface.h"
42 #include <fstream>
43 
44 //=============================================================================================================
45 // USED NAMESPACES
46 //=============================================================================================================
47 
48 using namespace MNELIB;
49 using namespace Eigen;
50 using namespace FIFFLIB;
51 
52 //=============================================================================================================
53 // DEFINE MEMBER METHODS
54 //=============================================================================================================
55 
57 : id(-1)
58 , np(-1)
59 , ntri(-1)
60 , coord_frame(-1)
61 , sigma(-1)
62 , rr(MatrixX3f::Zero(0,3))
63 , nn(MatrixX3f::Zero(0,3))
64 , tris(MatrixX3i::Zero(0,3))
65 , tri_cent(MatrixX3d::Zero(0,3))
66 , tri_nn(MatrixX3d::Zero(0,3))
67 , tri_area(VectorXd::Zero(0))
68 {
69 }
70 
71 //=============================================================================================================
72 
74 : id(p_MNEBemSurface.id)
75 , np(p_MNEBemSurface.np)
76 , ntri(p_MNEBemSurface.ntri)
77 , coord_frame(p_MNEBemSurface.coord_frame)
78 , sigma(p_MNEBemSurface.sigma)
79 , rr(p_MNEBemSurface.rr)
80 , nn(p_MNEBemSurface.nn)
81 , tris(p_MNEBemSurface.tris)
82 , tri_cent(p_MNEBemSurface.tri_cent)
83 , tri_nn(p_MNEBemSurface.tri_nn)
84 , tri_area(p_MNEBemSurface.tri_area)
85 , neighbor_tri(p_MNEBemSurface.neighbor_tri)
86 , neighbor_vert(p_MNEBemSurface.neighbor_vert)
87 {
88  //*m_pGeometryData = *p_MNEBemSurface.m_pGeometryData;
89 }
90 
91 //=============================================================================================================
92 
94 {
95 }
96 
97 //=============================================================================================================
98 
100 {
101  id = -1;
102  np = -1;
103  ntri = -1;
104  coord_frame = -1;
105  sigma = -1;
106  rr = MatrixX3f::Zero(0,3);
107  nn = MatrixX3f::Zero(0,3);
108  tris = MatrixX3i::Zero(0,3);
109  tri_cent = MatrixX3d::Zero(0,3);
110  tri_nn = MatrixX3d::Zero(0,3);
111  tri_area = VectorXd::Zero(0);
112 }
113 
114 //=============================================================================================================
115 
117 {
118  //
119  // Main triangulation
120  //
121  printf("\tCompleting triangulation info...");
122  this->tri_cent = MatrixX3d::Zero(this->ntri,3);
123  this->tri_nn = MatrixX3d::Zero(this->ntri,3);
124  this->tri_area = VectorXd::Zero(this->ntri);
125 
126  Matrix3d r;
127  Vector3d a, b;
128  int k = 0;
129  float size = 0;
130  for (qint32 i = 0; i < this->ntri; ++i)
131  {
132  for ( qint32 j = 0; j < 3; ++j)
133  {
134  k = this->tris(i, j);
135 
136  r(j,0) = this->rr(k, 0);
137  r(j,1) = this->rr(k, 1);
138  r(j,2) = this->rr(k, 2);
139 
140  this->tri_cent(i, 0) += this->rr(k, 0);
141  this->tri_cent(i, 1) += this->rr(k, 1);
142  this->tri_cent(i, 2) += this->rr(k, 2);
143  }
144  this->tri_cent.row(i) /= 3.0f;
145 
146  //cross product {cross((r2-r1),(r3-r1))}
147  a = r.row(1) - r.row(0 );
148  b = r.row(2) - r.row(0);
149  this->tri_nn(i,0) = a(1)*b(2)-a(2)*b(1);
150  this->tri_nn(i,1) = a(2)*b(0)-a(0)*b(2);
151  this->tri_nn(i,2) = a(0)*b(1)-a(1)*b(0);
152 
153  //area
154  size = this->tri_nn.row(i)*this->tri_nn.row(i).transpose();
155  size = std::pow(size, 0.5f );
156 
157  this->tri_area(i) = size/2.0f;
158  this->tri_nn.row(i) /= size;
159  }
160 
161  std::fstream doc("./Output/tri_area.dat", std::ofstream::out | std::ofstream::trunc);
162  if(doc) // if succesfully opened
163  {
164  // instructions
165  doc << this->tri_area << "\n";
166  doc.close();
167  }
168 
169  printf("Adding additional geometry info\n");
171 
172  printf("[done]\n");
173 
174  return true;
175 }
176 
177 //=============================================================================================================
178 
180 {
181  int k,c,p,q;
182  bool found;
183 
184  //Create neighboring triangle vector
185  neighbor_tri = QVector<QVector<int> >(this->tris.rows());
186 
187  //Create neighbor_tri information
188  for (p = 0; p < this->tris.rows(); p++) {
189  for (k = 0; k < 3; k++) {
190  this->neighbor_tri[this->tris(p,k)].append(p);
191  }
192  }
193 
194  //Create the neighboring vertices vector
195  neighbor_vert = QVector<QVector<int> >(this->np);
196 
197  for (k = 0; k < this->np; k++) {
198  for (p = 0; p < this->neighbor_tri[k].size(); p++) {
199  //Fit in the other vertices of the neighboring triangle
200  for (c = 0; c < 3; c++) {
201  int vert = this->tris(this->neighbor_tri[k][p], c);
202 
203  if (vert != k) {
204  found = false;
205 
206  for (q = 0; q < this->neighbor_vert[k].size(); q++) {
207  if (this->neighbor_vert[k][q] == vert) {
208  found = true;
209  break;
210  }
211  }
212 
213  if(!found) {
214  this->neighbor_vert[k].append(vert);
215  }
216  }
217  }
218  }
219  }
220 
221  return true;
222 }
223 
224 //=============================================================================================================
225 
227 {
228 
229  //
230  // Accumulate the vertex normals
231  //
232 
233 // this->nn.resize(this->np,3);
234 
235  for (qint32 p = 0; p < this->ntri; ++p) //check each triangle
236  {
237  for (qint32 j = 0; j < 3 ; ++j)
238  {
239  int nodenr;
240  nodenr = this->tris(p,j); //find the corners(nodes) of the triangles
241  this->nn(nodenr,0) += this->tri_nn(p,0); //add the triangle normal to the nodenormal
242  this->nn(nodenr,1) += this->tri_nn(p,1);
243  this->nn(nodenr,2) += this->tri_nn(p,2);
244  }
245  }
246 
247  // normalize
248  for (qint32 p = 0; p < this->np; ++p)
249  {
250  float size = 0;
251  size = this->nn.row(p)*this->nn.row(p).transpose();
252  size = std::pow(size, 0.5f );
253  this->nn.row(p) /= size;
254  }
255 
256 return true;
257 }
258 
259 //=============================================================================================================
260 
262 {
263  if(this->id <=0)
264  this->id=FIFFV_MNE_SURF_UNKNOWN;
265  if(this->sigma>0.0)
266  p_pStream->write_float(FIFF_BEM_SIGMA, &this->sigma);
267  p_pStream->write_int(FIFF_BEM_SURF_ID, &this->id);
268  p_pStream->write_int(FIFF_MNE_COORD_FRAME, &this->coord_frame);
269  p_pStream->write_int(FIFF_BEM_SURF_NNODE, &this->np);
270  p_pStream->write_int(FIFF_BEM_SURF_NTRI, &this->ntri);
271  p_pStream->write_float_matrix(FIFF_BEM_SURF_NODES, this->rr);
272  if (this->ntri > 0)
273  p_pStream->write_int_matrix(FIFF_BEM_SURF_TRIANGLES, this->tris.array() + 1);
274  p_pStream->write_float_matrix(FIFF_BEM_SURF_NORMALS, this->nn);
275 }
276 
277 //=============================================================================================================
278 
279 QString MNEBemSurface::id_name(int id)
280 {
281  switch(id) {
282  case FIFFV_BEM_SURF_ID_BRAIN: return "Brain";
283  case FIFFV_BEM_SURF_ID_SKULL: return "Skull";
284  case FIFFV_BEM_SURF_ID_HEAD: return "Head";
285  case FIFFV_BEM_SURF_ID_UNKNOWN: return "Unknown";
286  default: return "Unknown";
287  }
288 }
FIFFLIB::fiff_int_t id
Eigen::MatrixX3f nn
#define FIFF_BEM_SURF_TRIANGLES
Definition: fiff_file.h:739
FIFFLIB::fiff_float_t sigma
#define FIFF_BEM_SURF_NNODE
Definition: fiff_file.h:736
BEM surface provides geometry information.
fiff_long_t write_float_matrix(fiff_int_t kind, const Eigen::MatrixXf &mat)
Eigen::VectorXd tri_area
Eigen::MatrixX3f rr
#define FIFF_BEM_SURF_ID
Definition: fiff_file.h:734
fiff_long_t write_int_matrix(fiff_int_t kind, const Eigen::MatrixXi &mat)
FIFFLIB::fiff_int_t np
fiff_long_t write_float(fiff_int_t kind, const float *data, fiff_int_t nel=1)
Eigen::MatrixX3d tri_nn
fiff_long_t write_int(fiff_int_t kind, const fiff_int_t *data, fiff_int_t nel=1, fiff_int_t next=FIFFV_NEXT_SEQ)
void writeToStream(FIFFLIB::FiffStream *p_pStream)
FIFF File I/O routines.
Definition: fiff_stream.h:104
Eigen::MatrixX3d tri_cent
MNEBemSurface class declaration.
QVector< QVector< int > > neighbor_tri
#define FIFF_BEM_SURF_NTRI
Definition: fiff_file.h:737
QVector< QVector< int > > neighbor_vert
#define FIFF_BEM_SURF_NORMALS
Definition: fiff_file.h:740
Eigen::MatrixX3i tris
FIFFLIB::fiff_int_t ntri
static QString id_name(int id)
#define FIFF_BEM_SURF_NODES
Definition: fiff_file.h:738
#define FIFF_BEM_SIGMA
Definition: fiff_file.h:747
#define FIFF_MNE_COORD_FRAME
FIFFLIB::fiff_int_t coord_frame