MNE-CPP  0.1.9
A Framework for Electrophysiology
geometryinfo.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //=============================================================================================================
37 // INCLUDES
38 //=============================================================================================================
39 
40 #include "geometryinfo.h"
41 
42 #include <fiff/fiff_info.h>
43 
44 //=============================================================================================================
45 // INCLUDES
46 //=============================================================================================================
47 
48 #include <cmath>
49 #include <fstream>
50 #include <set>
51 
52 //=============================================================================================================
53 // QT INCLUDES
54 //=============================================================================================================
55 
56 #include <QtConcurrent/QtConcurrent>
57 
58 //=============================================================================================================
59 // EIGEN INCLUDES
60 //=============================================================================================================
61 
62 //=============================================================================================================
63 // USED NAMESPACES
64 //=============================================================================================================
65 
66 using namespace DISP3DLIB;
67 using namespace Eigen;
68 using namespace FIFFLIB;
69 
70 //=============================================================================================================
71 // DEFINE GLOBAL METHODS
72 //=============================================================================================================
73 
74 //=============================================================================================================
75 // DEFINE MEMBER METHODS
76 //=============================================================================================================
77 
78 QSharedPointer<MatrixXd> GeometryInfo::scdc(const MatrixX3f &matVertices,
79  const QVector<QVector<int> > &vecNeighborVertices,
80  QVector<int> &vecVertSubset,
81  double dCancelDist)
82 {
83  // create matrix and check for empty subset:
84  qint32 iCols = vecVertSubset.size();
85  if(vecVertSubset.empty()) {
86  // caller passed an empty subset, need to fill in all vertex IDs
87  qDebug() << "[WARNING] SCDC received empty subset, calculating full distance table, make sure you have enough memory !";
88  vecVertSubset.reserve(matVertices.rows());
89  for(qint32 id = 0; id < matVertices.rows(); ++id) {
90  vecVertSubset.push_back(id);
91  }
92  iCols = matVertices.rows();
93  }
94 
95  // convention: first dimension in distance table is "from", second dimension "to"
96  QSharedPointer<MatrixXd> returnMat = QSharedPointer<MatrixXd>::create(matVertices.rows(), iCols);
97 
98  // distribute calculation on cores
99  int iCores = QThread::idealThreadCount();
100  if (iCores <= 0) {
101  // assume that we have at least two available cores
102  iCores = 2;
103  }
104 
105  // start threads with their respective parts of the final subset
106  qint32 iSubArraySize = int(double(vecVertSubset.size()) / double(iCores));
107  QVector<QFuture<void> > vecThreads(iCores);
108  qint32 iBegin = 0;
109  qint32 iEnd = iSubArraySize;
110 
111  for (int i = 0; i < vecThreads.size(); ++i) {
112  //last round
113  if(i == vecThreads.size()-1)
114  {
115  vecThreads[i] = QtConcurrent::run(std::bind(iterativeDijkstra,
116  returnMat,
117  std::cref(matVertices),
118  std::cref(vecNeighborVertices),
119  std::cref(vecVertSubset),
120  iBegin,
121  vecVertSubset.size(),
122  dCancelDist));
123  break;
124  }
125  else
126  {
127  vecThreads[i] = QtConcurrent::run(std::bind(iterativeDijkstra,
128  returnMat,
129  std::cref(matVertices),
130  std::cref(vecNeighborVertices),
131  std::cref(vecVertSubset),
132  iBegin,
133  iEnd,
134  dCancelDist));
135  iBegin += iSubArraySize;
136  iEnd += iSubArraySize;
137  }
138  }
139 
140  // wait for all other threads to finish
141  for (QFuture<void>& f : vecThreads) {
142  f.waitForFinished();
143  }
144 
145  return returnMat;
146 }
147 
148 //=============================================================================================================
149 
150 QVector<int> GeometryInfo::projectSensors(const MatrixX3f &matVertices,
151  const QVector<Vector3f> &vecSensorPositions)
152 {
153  QVector<int> vecOutputArray;
154 
155  qint32 iCores = QThread::idealThreadCount();
156  if (iCores <= 0)
157  {
158  // assume that we have at least two available cores
159  iCores = 2;
160  }
161 
162  const qint32 iSubArraySize = int(double(vecSensorPositions.size()) / double(iCores));
163 
164  //small input size no threads needed
165  if(iSubArraySize <= 1)
166  {
167  vecOutputArray.append(nearestNeighbor(matVertices,
168  vecSensorPositions.constBegin(),
169  vecSensorPositions.constEnd()));
170  return vecOutputArray;
171  }
172 
173  // split input array + thread start
174  QVector<QFuture<QVector<int> > > vecThreads(iCores);
175  qint32 iBeginOffset = 0;
176  qint32 iEndOffset = iBeginOffset + iSubArraySize;
177  for(qint32 i = 0; i < vecThreads.size(); ++i)
178  {
179  //last round
180  if(i == vecThreads.size()-1)
181  {
182  vecThreads[i] = QtConcurrent::run(nearestNeighbor,
183  matVertices,
184  vecSensorPositions.constBegin() + iBeginOffset,
185  vecSensorPositions.constEnd());
186  break;
187  }
188  else
189  {
190  vecThreads[i] = QtConcurrent::run(nearestNeighbor,
191  matVertices,
192  vecSensorPositions.constBegin() + iBeginOffset,
193  vecSensorPositions.constBegin() + iEndOffset);
194  iBeginOffset = iEndOffset;
195  iEndOffset += iSubArraySize;
196  }
197  }
198 
199  //wait for threads to finish
200  for (QFuture<QVector<int> >& f : vecThreads) {
201  f.waitForFinished();
202  }
203 
204  //move sub arrays back into output
205  for(qint32 i = 0; i < vecThreads.size(); ++i)
206  {
207  vecOutputArray.append(std::move(vecThreads[i].result()));
208  }
209 
210  return vecOutputArray;
211 }
212 
213 //=============================================================================================================
214 
215 QVector<int> GeometryInfo::nearestNeighbor(const MatrixX3f &matVertices,
216  QVector<Vector3f>::const_iterator itSensorBegin,
217  QVector<Vector3f>::const_iterator itSensorEnd)
218 {
219  //lin search sensor positions
220  QVector<int> vecMappedSensors;
221  vecMappedSensors.reserve(std::distance(itSensorBegin, itSensorEnd));
222 int u =0;
223  for(auto sensor = itSensorBegin; sensor != itSensorEnd; ++sensor)
224  {
225  qint32 iChampionId;
226  double iChampDist = std::numeric_limits<double>::max();
227  for(qint32 i = 0; i < matVertices.rows(); ++i)
228  {
229  //calculate 3d euclidian distance
230  double dDist = sqrt(squared(matVertices(i, 0) - (*sensor)[0]) // x-cord
231  + squared(matVertices(i, 1) - (*sensor)[1]) // y-cord
232  + squared(matVertices(i, 2) - (*sensor)[2])); // z-cord
233  if(dDist < iChampDist)
234  {
235  iChampionId = i;
236  iChampDist = dDist;
237  }
238  }
239  // qDebug() << "[GeometryInfo::nearestNeighbor] u" << u++;
240  vecMappedSensors.push_back(iChampionId);
241  }
242 
243  return vecMappedSensors;
244 }
245 
246 //=============================================================================================================
247 
248 void GeometryInfo::iterativeDijkstra(QSharedPointer<MatrixXd> matOutputDistMatrix,
249  const MatrixX3f &matVertices,
250  const QVector<QVector<int> > &vecNeighborVertices,
251  const QVector<int> &vecVertSubset,
252  qint32 iBegin,
253  qint32 iEnd,
254  double dCancelDistance) {
255  // initialization
256  const QVector<QVector<int> > &vecAdjacency = vecNeighborVertices;
257  qint32 n = vecAdjacency.size();
258  QVector<double> vecMinDists(n);
259  std::set< std::pair< double, qint32> > vertexQ;
260  const double INF = FLOAT_INFINITY;
261 
262  // outer loop, iterated for each vertex of 'vertSubset' between 'begin' and 'end'
263  for (qint32 i = iBegin; i < iEnd; ++i) {
264  // init phase of dijkstra: set source node for current iteration and reset data fields
265  if(i ==123)
266  qDebug() << "blabla";
267  qint32 iRoot = vecVertSubset.at(i);
268  vertexQ.clear();
269  vecMinDists.fill(INF);
270  vecMinDists[iRoot] = 0.0;
271  vertexQ.insert(std::make_pair(vecMinDists[iRoot], iRoot));
272 
273  // dijkstra main loop
274  while (vertexQ.empty() == false) {
275  // remove next vertex from queue
276  const double dDist = vertexQ.begin()->first;
277  const qint32 u = vertexQ.begin()->second;
278  vertexQ.erase(vertexQ.begin());
279 
280  // check if we are still below cancel distance
281  if (dDist <= dCancelDistance) {
282  // visit each neighbour of u
283  const QVector<int>& vecNeighbours = vecAdjacency[u];
284 
285  for (qint32 ne = 0; ne < vecNeighbours.length(); ++ne) {
286  qint32 v = vecNeighbours[ne];
287 
288  // distance from source (i.e. root) to v, using u as its predecessor
289  // calculate inline since designated function was magnitudes slower (even when declared as inline)
290  const double dDistX = matVertices(u, 0) - matVertices(v, 0);
291  const double dDistY = matVertices(u, 1) - matVertices(v, 1);
292  const double dDistZ = matVertices(u, 2) - matVertices(v, 2);
293  const double dDistWithU = dDist + sqrt(dDistX * dDistX + dDistY * dDistY + dDistZ * dDistZ);
294 
295  if (dDistWithU < vecMinDists[v]) {
296  // this is a combination of insert and decreaseKey
297  vertexQ.erase(std::make_pair(vecMinDists[v], v));
298  vecMinDists[v] = dDistWithU;
299  vertexQ.insert(std::make_pair(vecMinDists[v], v));
300  }
301  }
302  }
303  }
304 
305  // save results for current root in matrix
306  for (qint32 m = 0; m < vecMinDists.size(); ++m) {
307  matOutputDistMatrix->coeffRef(m , i) = vecMinDists[m];
308  }
309  }
310 }
311 
312 //=============================================================================================================
313 
314 QVector<int> GeometryInfo::filterBadChannels(QSharedPointer<Eigen::MatrixXd> matDistanceTable,
315  const FIFFLIB::FiffInfo& fiffInfo,
316  qint32 iSensorType) {
317  // use pointer to avoid copying of FiffChInfo objects
318  QVector<int> vecBadColumns;
319  QVector<const FiffChInfo*> vecSensors;
320  for(const FiffChInfo& s : fiffInfo.chs){
321  //Only take EEG with V as unit or MEG magnetometers with T as unit
322  if(s.kind == iSensorType && (s.unit == FIFF_UNIT_T || s.unit == FIFF_UNIT_V)){
323  vecSensors.push_back(&s);
324  }
325  }
326 
327  // inefficient: going through all bad sensors, i.e. also the ones which are of different type than the passed one
328  for(const QString& b : fiffInfo.bads){
329  for(int col = 0; col < vecSensors.size(); ++col){
330  if(vecSensors[col]->ch_name == b){
331  // found index of our bad channel, set whole column to infinity
332  vecBadColumns.push_back(col);
333  for(int row = 0; row < matDistanceTable->rows(); ++row){
334  matDistanceTable->coeffRef(row, col) = FLOAT_INFINITY;
335  }
336  break;
337  }
338  }
339  }
340  return vecBadColumns;
341 }
Channel info descriptor.
Definition: fiff_ch_info.h:74
FIFF measurement file information.
Definition: fiff_info.h:84
static QVector< int > nearestNeighbor(const Eigen::MatrixX3f &matVertices, QVector< Eigen::Vector3f >::const_iterator itSensorBegin, QVector< Eigen::Vector3f >::const_iterator itSensorEnd)
nearestNeighbor Calculates the nearest vertex of an MNEmatVertices for each position between the two ...
GeometryInfo class declaration.
static QSharedPointer< Eigen::MatrixXd > scdc(const Eigen::MatrixX3f &matVertices, const QVector< QVector< int > > &vecNeighborVertices, QVector< int > &pVecVertSubset, double dCancelDist=FLOAT_INFINITY)
scdc Calculates surface constrained distances on a mesh.
static QVector< int > projectSensors(const Eigen::MatrixX3f &matVertices, const QVector< Eigen::Vector3f > &vecSensorPositions)
Calculates the nearest neighbor (euclidian distance) vertex to each sensor.
static void iterativeDijkstra(QSharedPointer< Eigen::MatrixXd > matOutputDistMatrix, const Eigen::MatrixX3f &matVertices, const QVector< QVector< int > > &vecNeighborVertices, const QVector< int > &vecVertSubset, qint32 iBegin, qint32 iEnd, double dCancelDistance)
iterativeDijkstra Calculates shortest distances on the mesh that is held by the MNEmatVertices for ea...
FiffInfo class declaration.
static QVector< int > filterBadChannels(QSharedPointer< Eigen::MatrixXd > matDistanceTable, const FIFFLIB::FiffInfo &fiffInfo, qint32 iSensorType)
filterBadChannels Filters bad channels from distance table
QList< FiffChInfo > chs