MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_project_to_surface.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //=============================================================================================================
37 // INCLUDES
38 //=============================================================================================================
39 
40 #include "mne_project_to_surface.h"
41 
42 //=============================================================================================================
43 // INCLUDES
44 //=============================================================================================================
45 
46 #include <mne/mne_bem_surface.h>
47 #include <mne/mne_surface.h>
48 
49 //=============================================================================================================
50 // QT INCLUDES
51 //=============================================================================================================
52 
53 //=============================================================================================================
54 // EIGEN INCLUDES
55 //=============================================================================================================
56 
57 #include <Eigen/Geometry>
58 
59 //=============================================================================================================
60 // USED NAMESPACES
61 //=============================================================================================================
62 
63 using namespace MNELIB;
64 using namespace Eigen;
65 
66 //=============================================================================================================
67 // DEFINE GLOBAL METHODS
68 //=============================================================================================================
69 
70 //=============================================================================================================
71 // DEFINE MEMBER METHODS
72 //=============================================================================================================
73 
75 : r1(MatrixX3f::Zero(1,3))
76 , r12(MatrixX3f::Zero(1,3))
77 , r13(MatrixX3f::Zero(1,3))
78 , nn(MatrixX3f::Zero(1,3))
79 , a(VectorXf::Zero(1))
80 , b(VectorXf::Zero(1))
81 , c(VectorXf::Zero(1))
82 , det(VectorXf::Zero(1))
83 {
84 }
85 
86 //=============================================================================================================
87 
89 : r1(MatrixX3f::Zero(p_MNEBemSurf.ntri,3))
90 , r12(MatrixX3f::Zero(p_MNEBemSurf.ntri,3))
91 , r13(MatrixX3f::Zero(p_MNEBemSurf.ntri,3))
92 , nn(MatrixX3f::Zero(p_MNEBemSurf.ntri,3))
93 , a(VectorXf::Zero(p_MNEBemSurf.ntri))
94 , b(VectorXf::Zero(p_MNEBemSurf.ntri))
95 , c(VectorXf::Zero(p_MNEBemSurf.ntri))
96 , det(VectorXf::Zero(p_MNEBemSurf.ntri))
97 {
98  for (int i = 0; i < p_MNEBemSurf.ntri; ++i)
99  {
100  r1.row(i) = p_MNEBemSurf.rr.row(p_MNEBemSurf.tris(i,0));
101  r12.row(i) = p_MNEBemSurf.rr.row(p_MNEBemSurf.tris(i,1)) - r1.row(i);
102  r13.row(i) = p_MNEBemSurf.rr.row(p_MNEBemSurf.tris(i,2)) - r1.row(i);
103  a(i) = r12.row(i) * r12.row(i).transpose();
104  b(i) = r13.row(i) * r13.row(i).transpose();
105  c(i) = r12.row(i) * r13.row(i).transpose();
106  }
107 
108  if (!(p_MNEBemSurf.tri_nn.isZero(0)))
109  {
110  nn = p_MNEBemSurf.tri_nn.cast<float>();
111  }
112  else
113  {
114  for (int i = 0; i < p_MNEBemSurf.ntri; ++i)
115  {
116  nn.row(i) = r12.row(i).transpose().cross(r13.row(i).transpose()).transpose();
117  }
118  }
119  det = (a.array()*b.array() - c.array()*c.array()).matrix();
120 }
121 
122 //=============================================================================================================
123 
125 : r1(MatrixX3f::Zero(p_MNESurf.ntri,3))
126 , r12(MatrixX3f::Zero(p_MNESurf.ntri,3))
127 , r13(MatrixX3f::Zero(p_MNESurf.ntri,3))
128 , nn(MatrixX3f::Zero(p_MNESurf.ntri,3))
129 , a(VectorXf::Zero(p_MNESurf.ntri))
130 , b(VectorXf::Zero(p_MNESurf.ntri))
131 , c(VectorXf::Zero(p_MNESurf.ntri))
132 , det(VectorXf::Zero(p_MNESurf.ntri))
133 {
134  for (int i = 0; i < p_MNESurf.ntri; ++i)
135  {
136  r1.row(i) = p_MNESurf.rr.row(p_MNESurf.tris(i,0));
137  r12.row(i) = p_MNESurf.rr.row(p_MNESurf.tris(i,1)) - r1.row(i);
138  r13.row(i) = p_MNESurf.rr.row(p_MNESurf.tris(i,2)) - r1.row(i);
139  nn.row(i) = r12.row(i).transpose().cross(r13.row(i).transpose()).transpose();
140  a(i) = r12.row(i) * r12.row(i).transpose();
141  b(i) = r13.row(i) * r13.row(i).transpose();
142  c(i) = r12.row(i) * r13.row(i).transpose();
143  }
144 
145  det = (a.array()*b.array() - c.array()*c.array()).matrix();
146 }
147 
148 //=============================================================================================================
149 
150 bool MNEProjectToSurface::mne_find_closest_on_surface(const MatrixXf &r, const int np, MatrixXf &rTri,
151  VectorXi &nearest, VectorXf &dist)
152 {
153  // resize output
154  nearest.resize(np);
155  dist.resize(np);
156  rTri.resize(np,3);
157 
158  if (this->r1.isZero(0))
159  {
160  qDebug() << "No surface loaded to make the projection./n";
161  return false;
162  }
163  int bestTri = -1;
164  float bestDist = -1;
165  Vector3f rTriK;
166  for (int k = 0; k < np; ++k)
167  {
168  /*
169  * To do: decide_search_restriction for the use in an iterative closest point to plane algorithm
170  * For now it's OK to go through all triangles.
171  */
172  if (!this->mne_project_to_surface(r.row(k).transpose(), rTriK, bestTri, bestDist))
173  {
174  qDebug() << "The projection of point number " << k << " didn't work./n";
175  return false;
176  }
177  rTri.row(k) = rTriK.transpose();
178  nearest[k] = bestTri;
179  dist[k] = bestDist;
180  }
181  return true;
182 }
183 
184 //=============================================================================================================
185 
186 bool MNEProjectToSurface::mne_project_to_surface(const Vector3f &r, Vector3f &rTri, int &bestTri, float &bestDist)
187 {
188  float p = 0, q = 0, p0 = 0, q0 = 0, dist0 = 0;
189  bestDist = 0.0f;
190  bestTri = -1;
191  for (int tri = 0; tri < a .size(); ++tri)
192  {
193  if (!this->nearest_triangle_point(r, tri, p0, q0, dist0))
194  {
195  qDebug() << "The projection on triangle " << tri << " didn't work./n";
196  return false;
197  }
198 
199  if ((bestTri < 0) || (std::fabs(dist0) < std::fabs(bestDist)))
200  {
201  bestDist = dist0;
202  p = p0;
203  q = q0;
204  bestTri = tri;
205  }
206  }
207 
208  if (bestTri >= 0)
209  {
210  if (!this->project_to_triangle(rTri, p, q, bestTri))
211  {
212  qDebug() << "The coordinate transform to cartesian system didn't work./n";
213  return false;
214  }
215  return true;
216  }
217 
218  qDebug() << "No best Triangle found./n";
219  return false;
220 }
221 
222 //=============================================================================================================
223 
224 bool MNEProjectToSurface::nearest_triangle_point(const Vector3f &r, const int tri, float &p, float &q, float &dist)
225 {
226  //Calculate some helpers
227  Vector3f rr = r - this->r1.row(tri).transpose(); //Vector from triangle corner #1 to r
228  float v1 = this->r12.row(tri)*rr;
229  float v2 = this->r13.row(tri)*rr;
230 
231  //Calculate the orthogonal projection of the point r on the plane
232  dist = this->nn.row(tri)*rr;
233  p = (this->b(tri)*v1 - this->c(tri)*v2)/det(tri);
234  q = (this->a(tri)*v2 - this->c(tri)*v1)/det(tri);
235 
236  //If the point projects into the triangle we are done
237  if (p >= 0.0 && p <= 1.0 && q >= 0.0 && q <= 1.0 && (p+q) <= 1.0)
238  {
239  return true;
240  }
241 
242  /*
243  * Tough: must investigate the sides
244  * We might do something intelligent here. However, for now it is ok
245  * to do it in the hard way
246  */
247  float p0, q0, t0, dist0, best, bestp, bestq;
248 
249  /*
250  * Side 1 -> 2
251  */
252  p0 = p + (q * this->c(tri)) / this->a(tri);
253  // Place the point in the corner if it is not on the side
254  if (p0 < 0.0)
255  {
256  p0 = 0.0;
257  }
258  else if (p0 > 1.0)
259  {
260  p0 = 1.0;
261  }
262  q0 = 0;
263  // Distance
264  dist0 = sqrt((p-p0)*(p-p0)*this->a(tri) +
265  (q-q0)*(q-q0)*this->b(tri) +
266  2*(p-p0)*(q-q0)*this->c(tri) +
267  dist*dist);
268 
269  best = dist0;
270  bestp = p0;
271  bestq = q0;
272  /*
273  * Side 2 -> 3
274  */
275  t0 = ((a(tri)-c(tri))*(-p) + (b(tri)-c(tri))*q)/(a(tri)+b(tri)-2*c(tri));
276  // Place the point in the corner if it is not on the side
277  if (t0 < 0.0)
278  {
279  t0 = 0.0;
280  }
281  else if (t0 > 1.0)
282  {
283  t0 = 1.0;
284  }
285  p0 = 1.0 - t0;
286  q0 = t0;
287  // Distance
288  dist0 = sqrt((p-p0)*(p-p0)*this->a(tri) +
289  (q-q0)*(q-q0)*this->b(tri) +
290  2*(p-p0)*(q-q0)*this->c(tri) +
291  dist*dist);
292  if (dist0 < best)
293  {
294  best = dist0;
295  bestp = p0;
296  bestq = q0;
297  }
298  /*
299  * Side 1 -> 3
300  */
301  p0 = 0.0;
302  q0 = q + (p * c(tri))/b(tri);
303  // Place the point in the corner if it is not on the side
304  if (q0 < 0.0)
305  {
306  q0 = 0.0;
307 
308  }
309  else if (q0 > 1.0)
310  {
311  q0 = 1.0;
312  }
313  // Distance
314  dist0 = sqrt((p-p0)*(p-p0)*this->a(tri) +
315  (q-q0)*(q-q0)*this->b(tri) +
316  2*(p-p0)*(q-q0)*this->c(tri) +
317  dist*dist);
318  if (dist0 < best)
319  {
320  best = dist0;
321  bestp = p0;
322  bestq = q0;
323  }
324  dist = best;
325  p = bestp;
326  q = bestq;
327  return true;
328 }
329 
330 //=============================================================================================================
331 
332 bool MNEProjectToSurface::project_to_triangle(Vector3f &rTri, const float p, const float q, const int tri)
333 {
334  rTri = this->r1.row(tri) + p*this->r12.row(tri) + q*this->r13.row(tri);
335  return true;
336 }
BEM Surface.
Definition: mne_surface.h:82
BEM surface provides geometry information.
FIFFLIB::fiff_int_t ntri
Definition: mne_surface.h:137
Contains the declaration of the MNESurface class.
Eigen::MatrixX3f rr
Eigen::MatrixX3d tri_nn
MNEProjectToSurface class declaration.
MNEBemSurface class declaration.
Eigen::MatrixX3i tris
FIFFLIB::fiff_int_t ntri
bool mne_find_closest_on_surface(const Eigen::MatrixXf &r, const int np, Eigen::MatrixXf &rTri, Eigen::VectorXi &nearest, Eigen::VectorXf &dist)
mne_find_closest_on_surface