MNE-CPP  0.1.9
A Framework for Electrophysiology
rtfiffrawviewmodel.cpp
Go to the documentation of this file.
1 //=============================================================================================================
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include "rtfiffrawviewmodel.h"
43 
44 #include <fiff/fiff_types.h>
45 #include <fiff/fiff_info.h>
46 
47 #include <utils/mnemath.h>
48 #include <utils/ioutils.h>
49 
50 #include <rtprocessing/sphara.h>
52 
53 //=============================================================================================================
54 // QT INCLUDES
55 //=============================================================================================================
56 
57 #include <QBrush>
58 #include <QCoreApplication>
59 #include <QtConcurrent>
60 #include <QFuture>
61 #include <QDebug>
62 
63 //=============================================================================================================
64 // EIGEN INCLUDES
65 //=============================================================================================================
66 
67 //=============================================================================================================
68 // USED NAMESPACES
69 //=============================================================================================================
70 
71 using namespace DISPLIB;
72 using namespace UTILSLIB;
73 using namespace FIFFLIB;
74 using namespace Eigen;
75 using namespace RTPROCESSINGLIB;
76 
77 //=============================================================================================================
78 // DEFINE MEMBER METHODS
79 //=============================================================================================================
80 
81 EVENTSLIB::EventManager RtFiffRawViewModel::m_EventManager;
82 
84 : QAbstractTableModel(parent)
85 , m_bProjActivated(false)
86 , m_bCompActivated(false)
87 , m_bSpharaActivated(false)
88 , m_bIsFreezed(false)
89 , m_bDrawFilterFront(true)
90 , m_bPerformFiltering(false)
91 , m_bTriggerDetectionActive(false)
92 , m_fSps(1024.0f)
93 , m_dTriggerThreshold(0.01)
94 , m_iT(10)
95 , m_iDownsampling(10)
96 , m_iMaxSamples(1024)
97 , m_iCurrentSample(0)
98 , m_iCurrentStartingSample(0)
99 , m_iCurrentSampleFreeze(0)
100 , m_iMaxFilterLength(128)
101 , m_iCurrentBlockSize(1024)
102 , m_iResidual(0)
103 , m_iCurrentTriggerChIndex(0)
104 , m_iDistanceTimerSpacer(1000)
105 , m_iDetectedTriggers(0)
106 , m_sFilterChannelType("MEG")
107 , m_pFiffInfo(FiffInfo::SPtr::create())
108 , m_colBackground(Qt::white)
109 {
110  m_EventManager.initSharedMemory(EVENTSLIB::SharedMemoryMode::READWRITE);
111 }
112 
113 //=============================================================================================================
114 
116 {
117  m_EventManager.stopSharedMemory();
118 }
119 
120 //=============================================================================================================
121 //virtual functions
122 int RtFiffRawViewModel::rowCount(const QModelIndex & /*parent*/) const
123 {
124  if(!m_pFiffInfo->chs.isEmpty()) {
125  return m_pFiffInfo->chs.size();
126  } else {
127  return 0;
128  }
129 }
130 
131 //=============================================================================================================
132 
133 int RtFiffRawViewModel::columnCount(const QModelIndex & /*parent*/) const
134 {
135  return 3;
136 }
137 
138 //=============================================================================================================
139 
140 QVariant RtFiffRawViewModel::data(const QModelIndex &index, int role) const
141 {
142  if(role != Qt::DisplayRole && role != Qt::BackgroundRole) {
143  return QVariant();
144  }
145 
146  if (role == Qt::BackgroundRole) {
147  return QVariant(QBrush(m_colBackground));
148  }
149 
150  if (index.isValid()) {
151  qint32 row = m_qMapIdxRowSelection.value(index.row(),0);
152 
153  //******** first column (chname) ********
154  if(index.column() == 0 && role == Qt::DisplayRole)
155  return QVariant(m_pFiffInfo->ch_names[row]);
156 
157  //******** second column (data plot) ********
158  if(index.column() == 1) {
159  QVariant v;
160  RowVectorPair rowVectorPair;
161 
162  switch(role) {
163  case Qt::DisplayRole: {
164  if(m_bIsFreezed) {
165  // data freeze
166  if(!m_filterKernel.isEmpty() && m_bPerformFiltering) {
167  rowVectorPair.first = m_matDataFilteredFreeze.data() + row*m_matDataFilteredFreeze.cols();
168  rowVectorPair.second = m_matDataFilteredFreeze.cols();
169  v.setValue(rowVectorPair);
170  } else {
171  rowVectorPair.first = m_matDataRawFreeze.data() + row*m_matDataRawFreeze.cols();
172  rowVectorPair.second = m_matDataRawFreeze.cols();
173  v.setValue(rowVectorPair);
174  }
175  }
176  else {
177  // data stream
178  if(!m_filterKernel.isEmpty() && m_bPerformFiltering) {
179  rowVectorPair.first = m_matDataFiltered.data() + row*m_matDataFiltered.cols();
180  rowVectorPair.second = m_matDataFiltered.cols();
181  v.setValue(rowVectorPair);
182  } else {
183  rowVectorPair.first = m_matDataRaw.data() + row*m_matDataRaw.cols();
184  rowVectorPair.second = m_matDataRaw.cols();
185  v.setValue(rowVectorPair);
186  }
187  }
188 
189  return v;
190  }
191  } // end role switch
192  } // end column check
193 
194  //******** third column (bad channel) ********
195  if(index.column() == 2 && role == Qt::DisplayRole) {
196  return QVariant(m_pFiffInfo->bads.contains(m_pFiffInfo->ch_names[row]));
197  } // end column check
198 
199  } // end index.valid() check
200 
201  return QVariant();
202 }
203 
204 //=============================================================================================================
205 
206 QVariant RtFiffRawViewModel::headerData(int section, Qt::Orientation orientation, int role) const
207 {
208  if(role != Qt::DisplayRole && role != Qt::TextAlignmentRole)
209  return QVariant();
210 
211  if(orientation == Qt::Horizontal) {
212  switch(section) {
213  case 0: //chname column
214  return QVariant();
215  case 1: //data plot column
216  switch(role) {
217  case Qt::DisplayRole:
218  return QVariant("data plot");
219  case Qt::TextAlignmentRole:
220  return QVariant(Qt::AlignLeft);
221  }
222  return QVariant("data plot");
223  }
224  }
225  else if(orientation == Qt::Vertical) {
226  QModelIndex chname = createIndex(section,0);
227  switch(role) {
228  case Qt::DisplayRole:
229  return QVariant(data(chname).toString());
230  }
231  }
232 
233  return QVariant();
234 }
235 
236 //=============================================================================================================
237 
238 void RtFiffRawViewModel::initSphara()
239 {
240  //Load SPHARA matrices for babymeg and vectorview
241  IOUtils::read_eigen_matrix(m_matSpharaVVGradLoaded, QCoreApplication::applicationDirPath() + QString("/resources/mne_scan/plugins/noisereduction/SPHARA/Vectorview_SPHARA_InvEuclidean_Grad.txt"));
242  IOUtils::read_eigen_matrix(m_matSpharaVVMagLoaded, QCoreApplication::applicationDirPath() + QString("/resources/mne_scan/plugins/noisereduction/SPHARA/Vectorview_SPHARA_InvEuclidean_Mag.txt"));
243 
244  IOUtils::read_eigen_matrix(m_matSpharaBabyMEGInnerLoaded, QCoreApplication::applicationDirPath() + QString("/resources/mne_scan/plugins/noisereduction/SPHARA/BabyMEG_SPHARA_InvEuclidean_Inner.txt"));
245  IOUtils::read_eigen_matrix(m_matSpharaBabyMEGOuterLoaded, QCoreApplication::applicationDirPath() + QString("/resources/mne_scan/plugins/noisereduction/SPHARA/BabyMEG_SPHARA_InvEuclidean_Outer.txt"));
246 
247  IOUtils::read_eigen_matrix(m_matSpharaEEGLoaded, QCoreApplication::applicationDirPath() + QString("/resources/mne_scan/plugins/noisereduction/SPHARA/Current_SPHARA_EEG.txt"));
248 
249  //Generate indices used to create the SPHARA operators for VectorView
250  m_vecIndicesFirstVV.resize(0);
251  m_vecIndicesSecondVV.resize(0);
252 
253  for(int r = 0; r < m_pFiffInfo->chs.size(); ++r) {
254  //Find GRADIOMETERS
255  if(m_pFiffInfo->chs.at(r).chpos.coil_type == 3012) {
256  m_vecIndicesFirstVV.conservativeResize(m_vecIndicesFirstVV.rows()+1);
257  m_vecIndicesFirstVV(m_vecIndicesFirstVV.rows()-1) = r;
258  }
259 
260  //Find Magnetometers
261  if(m_pFiffInfo->chs.at(r).chpos.coil_type == 3024) {
262  m_vecIndicesSecondVV.conservativeResize(m_vecIndicesSecondVV.rows()+1);
263  m_vecIndicesSecondVV(m_vecIndicesSecondVV.rows()-1) = r;
264  }
265  }
266 
267  //Generate indices used to create the SPHARA operators for babyMEG
268  m_vecIndicesFirstBabyMEG.resize(0);
269  for(int r = 0; r < m_pFiffInfo->chs.size(); ++r) {
270  //Find INNER LAYER
271  if(m_pFiffInfo->chs.at(r).chpos.coil_type == 7002) {
272  m_vecIndicesFirstBabyMEG.conservativeResize(m_vecIndicesFirstBabyMEG.rows()+1);
273  m_vecIndicesFirstBabyMEG(m_vecIndicesFirstBabyMEG.rows()-1) = r;
274  }
275 
276  //TODO: Find outer layer
277  }
278 
279  //Generate indices used to create the SPHARA operators for EEG layouts
280  m_vecIndicesFirstEEG.resize(0);
281  for(int r = 0; r < m_pFiffInfo->chs.size(); ++r) {
282  //Find EEG
283  if(m_pFiffInfo->chs.at(r).kind == FIFFV_EEG_CH) {
284  m_vecIndicesFirstEEG.conservativeResize(m_vecIndicesFirstEEG.rows()+1);
285  m_vecIndicesFirstEEG(m_vecIndicesFirstEEG.rows()-1) = r;
286  }
287  }
288 
289  //Create Sphara operator for the first time
290  updateSpharaOptions("BabyMEG", 270, 105);
291 
292  qDebug()<<"RtFiffRawViewModel::initSphara - Read VectorView mag matrix "<<m_matSpharaVVMagLoaded.rows()<<m_matSpharaVVMagLoaded.cols()<<"and grad matrix"<<m_matSpharaVVGradLoaded.rows()<<m_matSpharaVVGradLoaded.cols();
293  qDebug()<<"RtFiffRawViewModel::initSphara - Read BabyMEG inner layer matrix "<<m_matSpharaBabyMEGInnerLoaded.rows()<<m_matSpharaBabyMEGInnerLoaded.cols()<<"and outer layer matrix"<<m_matSpharaBabyMEGOuterLoaded.rows()<<m_matSpharaBabyMEGOuterLoaded.cols();
294 }
295 
296 //=============================================================================================================
297 
298 void RtFiffRawViewModel::setFiffInfo(QSharedPointer<FIFFLIB::FiffInfo> &p_pFiffInfo)
299 {
300  if(p_pFiffInfo) {
301  RowVectorXi sel;// = RowVectorXi(0,0);
302  QStringList emptyExclude;
303 
304  if(p_pFiffInfo->bads.size() > 0) {
305  sel = FiffInfoBase::pick_channels(p_pFiffInfo->ch_names, p_pFiffInfo->bads, emptyExclude);
306  }
307 
308  m_vecBadIdcs = sel;
309 
310  m_pFiffInfo = p_pFiffInfo;
311 
312  resetSelection();
313 
314  //Resize data matrix without touching the stored values
315  m_matDataRaw.conservativeResize(m_pFiffInfo->chs.size(), m_iMaxSamples);
316  m_matDataRaw.setZero();
317 
318  m_matDataFiltered.conservativeResize(m_pFiffInfo->chs.size(), m_iMaxSamples);
319  m_matDataFiltered.setZero();
320 
321  m_vecLastBlockFirstValuesFiltered.conservativeResize(m_pFiffInfo->chs.size());
322  m_vecLastBlockFirstValuesFiltered.setZero();
323 
324  m_vecLastBlockFirstValuesRaw.conservativeResize(m_pFiffInfo->chs.size());
325  m_vecLastBlockFirstValuesRaw.setZero();
326 
327  m_matOverlap.conservativeResize(m_pFiffInfo->chs.size(), m_iMaxFilterLength);
328 
329  m_matSparseProjMult = SparseMatrix<double>(m_pFiffInfo->chs.size(),m_pFiffInfo->chs.size());
330  m_matSparseCompMult = SparseMatrix<double>(m_pFiffInfo->chs.size(),m_pFiffInfo->chs.size());
331  m_matSparseSpharaMult = SparseMatrix<double>(m_pFiffInfo->chs.size(),m_pFiffInfo->chs.size());
332  m_matSparseProjCompMult = SparseMatrix<double>(m_pFiffInfo->chs.size(),m_pFiffInfo->chs.size());
333 
334  m_matSparseProjMult.setIdentity();
335  m_matSparseCompMult.setIdentity();
336  m_matSparseSpharaMult.setIdentity();
337  m_matSparseProjCompMult.setIdentity();
338 
339  //Create the initial Compensator projector
341 
342  //Initialize filter channel names
343  int visibleInit = 20;
344  QStringList filterChannels;
345 
346  if(visibleInit > m_pFiffInfo->chs.size()) {
347  while(visibleInit>m_pFiffInfo->chs.size()) {
348  visibleInit--;
349  }
350  }
351 
352  for(qint32 b = 0; b < visibleInit; ++b) {
353  filterChannels.append(m_pFiffInfo->ch_names.at(b));
354  }
355 
356  createFilterChannelList(filterChannels);
357 
358 // //Look for trigger channels and initialise detected trigger map
359 // for(int i = 0; i<m_pFiffInfo->chs.size(); ++i) {
360 // if(m_pFiffInfo->chs[i].kind == FIFFV_STIM_CH/* && m_pFiffInfo->chs[i].ch_name == "STI 001"*/)
361 // m_lTriggerChannelIndices.append(i);
362 // }
363 
364  //Init the sphara operators
365  initSphara();
366  } else {
367  m_vecBadIdcs = RowVectorXi(0,0);
368  m_matProj = MatrixXd(0,0);
369  m_matComp = MatrixXd(0,0);
370  }
371 }
372 
373 //=============================================================================================================
374 
375 void RtFiffRawViewModel::setSamplingInfo(float sps, int T, bool bSetZero)
376 {
377  beginResetModel();
378 
379  m_iT = T;
380 
381  m_iMaxSamples = (qint32) ceil(sps * T);
382 
383  //Resize data matrix without touching the stored values
384  m_matDataRaw.conservativeResize(m_pFiffInfo->chs.size(), m_iMaxSamples);
385  m_matDataFiltered.conservativeResize(m_pFiffInfo->chs.size(), m_iMaxSamples);
386  m_vecLastBlockFirstValuesRaw.conservativeResize(m_pFiffInfo->chs.size());
387  m_vecLastBlockFirstValuesFiltered.conservativeResize(m_pFiffInfo->chs.size());
388 
389  if(bSetZero) {
390  m_matDataRaw.setZero();
391  m_matDataFiltered.setZero();
392  m_vecLastBlockFirstValuesRaw.setZero();
393  m_vecLastBlockFirstValuesFiltered.setZero();
394  }
395 
396  if(m_iCurrentSample>m_iMaxSamples) {
397  m_iCurrentStartingSample += m_iCurrentSample;
398  m_iCurrentSample = 0;
399  }
400 
401  endResetModel();
402 }
403 
404 //=============================================================================================================
405 
407 {
408  if(!m_filterKernel.isEmpty() && m_bPerformFiltering) {
409  return m_matDataFiltered.block(0, m_iCurrentSample-m_iCurrentBlockSize, m_matDataFiltered.rows(), m_iCurrentBlockSize);
410  }
411 
412  return m_matDataRaw.block(0, m_iCurrentSample-m_iCurrentBlockSize, m_matDataRaw.rows(), m_iCurrentBlockSize);
413 }
414 
415 //=============================================================================================================
416 
417 void RtFiffRawViewModel::addData(const QList<MatrixXd> &data)
418 {
419  //SSP
420  bool doProj = m_bProjActivated && m_matDataRaw.cols() > 0 && m_matDataRaw.rows() == m_matProj.cols() ? true : false;
421 
422  //Compensator
423  bool doComp = m_bCompActivated && m_matDataRaw.cols() > 0 && m_matDataRaw.rows() == m_matComp.cols() ? true : false;
424 
425  //SPHARA
426  bool doSphara = m_bSpharaActivated && m_matSparseSpharaMult.cols() > 0 && m_matDataRaw.rows() == m_matSparseSpharaMult.cols() ? true : false;
427 
428  //Copy new data into the global data matrix
429  for(qint32 b = 0; b < data.size(); ++b) {
430  int nCol = data.at(b).cols();
431  int nRow = data.at(b).rows();
432 
433  if(nRow != m_matDataRaw.rows()) {
434  qDebug()<<"incoming data does not match internal data row size. Returning...";
435  return;
436  }
437 
438  //Reset m_iCurrentSample and start filling the data matrix from the beginning again. Also add residual amount of data to the end of the matrix.
439  if(m_iCurrentSample+nCol > m_matDataRaw.cols()) {
440  m_iResidual = nCol - ((m_iCurrentSample+nCol) % m_matDataRaw.cols());
441 
442  if(m_iResidual == nCol) {
443  m_iResidual = 0;
444  }
445 
446 // std::cout<<"incoming data exceeds internal data cols by: "<<(m_iCurrentSample+nCol) % m_matDataRaw.cols()<<std::endl;
447 // std::cout<<"m_iCurrentSample+nCol: "<<m_iCurrentSample+nCol<<std::endl;
448 // std::cout<<"m_matDataRaw.cols(): "<<m_matDataRaw.cols()<<std::endl;
449 // std::cout<<"nCol-m_iResidual: "<<nCol-m_iResidual<<std::endl<<std::endl;
450 
451  if(doComp) {
452  if(doProj) {
453  //Comp + Proj
454  m_matDataRaw.block(0, m_iCurrentSample, nRow, m_iResidual) = m_matSparseProjCompMult * data.at(b).block(0,0,nRow,m_iResidual);
455  } else {
456  //Comp
457  m_matDataRaw.block(0, m_iCurrentSample, nRow, m_iResidual) = m_matSparseCompMult * data.at(b).block(0,0,nRow,m_iResidual);
458  }
459  } else {
460  if(doProj)
461  {
462  //Proj
463  m_matDataRaw.block(0, m_iCurrentSample, nRow, m_iResidual) = m_matSparseProjMult * data.at(b).block(0,0,nRow,m_iResidual);
464  } else {
465  //None - Raw
466  m_matDataRaw.block(0, m_iCurrentSample, nRow, m_iResidual) = data.at(b).block(0,0,nRow,m_iResidual);
467  }
468  }
469 
470  m_iCurrentStartingSample += m_iCurrentSample;
471  m_iCurrentStartingSample += m_iResidual;
472 
473  m_iCurrentSample = 0;
474 
475  if(!m_bIsFreezed) {
476  m_vecLastBlockFirstValuesFiltered = m_matDataFiltered.col(0);
477  m_vecLastBlockFirstValuesRaw = m_matDataRaw.col(0);
478  }
479 
480  //Store old detected triggers
481  m_qMapDetectedTriggerOld = m_qMapDetectedTrigger;
482 
483  //Clear detected triggers
484  if(m_bTriggerDetectionActive) {
485  QMutableMapIterator<int,QList<QPair<int,double> > > i(m_qMapDetectedTrigger);
486  while (i.hasNext()) {
487  i.next();
488  i.value().clear();
489  }
490  }
491  } else {
492  m_iResidual = 0;
493  }
494 
495  //std::cout<<"incoming data is ok"<<std::endl;
496 
497  if(doComp) {
498  if(doProj) {
499  //Comp + Proj
500  m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol) = m_matSparseProjCompMult * data.at(b);
501  } else {
502  //Comp
503  m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol) = m_matSparseCompMult * data.at(b);
504  }
505  } else {
506  if(doProj) {
507  //Proj
508  m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol) = m_matSparseProjMult * data.at(b);
509  } else {
510  //None - Raw
511  m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol) = data.at(b);
512  }
513  }
514 
515  //Filter if neccessary else set filtered data matrix to zero
516  if(!m_filterKernel.isEmpty() && m_bPerformFiltering) {
517  filterDataBlock(m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol), m_iCurrentSample);
518 
519  //Perform SPHARA on filtered data after actual filtering - SPHARA should be applied on the best possible data
520  if(doSphara) {
521  if(m_iCurrentSample-m_iMaxFilterLength/2 >= 0) {
522  m_matDataFiltered.block(0, m_iCurrentSample-m_iMaxFilterLength/2, nRow, nCol) = m_matSparseSpharaMult * m_matDataFiltered.block(0, m_iCurrentSample-m_iMaxFilterLength/2, nRow, nCol);
523  }
524  else {
525  if(m_iCurrentSample-m_iMaxFilterLength/2 < 0) {
526  m_matDataFiltered.block(0, 0, nRow, nCol) = m_matSparseSpharaMult * m_matDataFiltered.block(0, 0, nRow, nCol);
527  int iResidual = m_iResidual+m_iMaxFilterLength/2;
528  m_matDataFiltered.block(0, m_matDataFiltered.cols()-iResidual, nRow, iResidual) = m_matSparseSpharaMult * m_matDataFiltered.block(0, m_matDataFiltered.cols()-iResidual, nRow, iResidual);
529  }
530  }
531  }
532  } else {
533  m_matDataFiltered.block(0, m_iCurrentSample, nRow, nCol).setZero();// = m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol);
534 
535  //Perform SPHARA on raw data data
536  if(doSphara) {
537  m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol) = m_matSparseSpharaMult * m_matDataRaw.block(0, m_iCurrentSample, nRow, nCol);
538  }
539  }
540 
541  m_iCurrentSample += nCol;
542  m_iCurrentBlockSize = nCol;
543 
544  //detect the trigger flanks in the trigger channels
545  if(m_bTriggerDetectionActive) {
546  int iOldDetectedTriggers = m_qMapDetectedTrigger[m_iCurrentTriggerChIndex].size();
547 
548  QList<QPair<int,double> > qMapDetectedTrigger = RTPROCESSINGLIB::detectTriggerFlanksMax(data.at(b), m_iCurrentTriggerChIndex, m_iCurrentSample-nCol, m_dTriggerThreshold, true, 500);
549  //QList<QPair<int,double> > qMapDetectedTrigger = RTPROCESSINGLIB::detectTriggerFlanksGrad(data.at(b), m_iCurrentTriggerChIndex, m_iCurrentSample-nCol, m_dTriggerThreshold, false, "Rising");
550 
551  //Append results to already found triggers
552  m_qMapDetectedTrigger[m_iCurrentTriggerChIndex].append(qMapDetectedTrigger);
553 
554  //Compute newly counted triggers
555  int newTriggers = m_qMapDetectedTrigger[m_iCurrentTriggerChIndex].size() - iOldDetectedTriggers;
556 
557  if(newTriggers!=0) {
558  m_iDetectedTriggers += newTriggers;
559  emit triggerDetected(m_iDetectedTriggers, m_qMapDetectedTrigger);
560  }
561  }
562  }
563 
564  //Update data content
565  QModelIndex topLeft = this->index(0,1);
566  QModelIndex bottomRight = this->index(m_pFiffInfo->ch_names.size()-1,1);
567  QVector<int> roles; roles << Qt::DisplayRole;
568 
569  emit dataChanged(topLeft, bottomRight, roles);
570 }
571 
572 //=============================================================================================================
573 
574 fiff_int_t RtFiffRawViewModel::getKind(qint32 row) const
575 {
576  if(row < m_qMapIdxRowSelection.size()) {
577  qint32 chRow = m_qMapIdxRowSelection[row];
578  return m_pFiffInfo->chs.at(chRow).kind;
579  }
580 
581  return 0;
582 }
583 
584 //=============================================================================================================
585 
586 fiff_int_t RtFiffRawViewModel::getUnit(qint32 row) const
587 {
588  if(row < m_qMapIdxRowSelection.size()) {
589  qint32 chRow = m_qMapIdxRowSelection[row];
590  return m_pFiffInfo->chs.at(chRow).unit;
591  }
592 
593  return FIFF_UNIT_NONE;
594 }
595 
596 //=============================================================================================================
597 
598 fiff_int_t RtFiffRawViewModel::getCoil(qint32 row) const
599 {
600  if(row < m_qMapIdxRowSelection.size()) {
601  qint32 chRow = m_qMapIdxRowSelection[row];
602  return m_pFiffInfo->chs.at(chRow).chpos.coil_type;
603  }
604 
605  return FIFFV_COIL_NONE;
606 }
607 
608 //=============================================================================================================
609 
610 void RtFiffRawViewModel::selectRows(const QList<qint32> &selection)
611 {
612  beginResetModel();
613 
614  m_qMapIdxRowSelection.clear();
615 
616  qint32 count = 0;
617  for(qint32 i = 0; i < selection.size(); ++i) {
618  if(selection[i] < m_pFiffInfo->chs.size()) {
619  m_qMapIdxRowSelection.insert(count,selection[i]);
620  ++count;
621  }
622  }
623 
624  emit newSelection(selection);
625 
626  endResetModel();
627 }
628 
629 //=============================================================================================================
630 
631 void RtFiffRawViewModel::hideRows(const QList<qint32> &selection)
632 {
633  beginResetModel();
634 
635  for(qint32 i = 0; i < selection.size(); ++i) {
636  if(m_qMapIdxRowSelection.contains(selection.at(i))) {
637  m_qMapIdxRowSelection.remove(selection.at(i));
638  }
639  }
640 
641  emit newSelection(selection);
642 
643  endResetModel();
644 }
645 
646 //=============================================================================================================
647 
649 {
650  beginResetModel();
651 
652  m_qMapIdxRowSelection.clear();
653 
654  for(qint32 i = 0; i < m_pFiffInfo->chs.size(); ++i) {
655  m_qMapIdxRowSelection.insert(i,i);
656  }
657 
658  endResetModel();
659 }
660 
661 //=============================================================================================================
662 
663 void RtFiffRawViewModel::toggleFreeze(const QModelIndex &)
664 {
665  m_bIsFreezed = !m_bIsFreezed;
666 
667  if(m_bIsFreezed) {
668  m_matDataRawFreeze = m_matDataRaw;
669  m_matDataFilteredFreeze = m_matDataFiltered;
670  m_qMapDetectedTriggerFreeze = m_qMapDetectedTrigger;
671  m_qMapDetectedTriggerOldFreeze = m_qMapDetectedTriggerOld;
672 
673  m_iCurrentSampleFreeze = m_iCurrentSample;
674  }
675 
676  //Update data content
677  QModelIndex topLeft = this->index(0,1);
678  QModelIndex bottomRight = this->index(m_pFiffInfo->chs.size()-1,1);
679  QVector<int> roles; roles << Qt::DisplayRole;
680 
681  emit dataChanged(topLeft, bottomRight, roles);
682 }
683 
684 //=============================================================================================================
685 
686 void RtFiffRawViewModel::setScaling(const QMap< qint32,float >& p_qMapChScaling)
687 {
688  beginResetModel();
689  m_qMapChScaling = p_qMapChScaling;
690  endResetModel();
691 }
692 
693 //=============================================================================================================
694 
695 void RtFiffRawViewModel::updateProjection(const QList<FIFFLIB::FiffProj>& projs)
696 {
697  // Update the SSP projector
698  if(m_pFiffInfo) {
699  //If a minimum of one projector is active set m_bProjActivated to true so that this model applies the ssp to the incoming data
700  m_bProjActivated = false;
701  m_pFiffInfo->projs = projs;
702 
703  for(qint32 i = 0; i < this->m_pFiffInfo->projs.size(); ++i) {
704  if(this->m_pFiffInfo->projs[i].active) {
705  m_bProjActivated = true;
706  break;
707  }
708  }
709 
710  this->m_pFiffInfo->make_projector(m_matProj);
711 
712  qDebug() << "RtFiffRawViewModel::updateProjection - New projection calculated.";
713 
714  //set columns of matrix to zero depending on bad channels indexes
715  for(qint32 j = 0; j < m_vecBadIdcs.cols(); ++j) {
716  m_matProj.col(m_vecBadIdcs[j]).setZero();
717  }
718 
719 // std::cout << "Bads\n" << m_vecBadIdcs << std::endl;
720 // std::cout << "Proj\n";
721 // std::cout << m_matProj.block(0,0,10,10) << std::endl;
722 
723  //
724  // Make proj sparse
725  //
726  qint32 nchan = this->m_pFiffInfo->nchan;
727  qint32 i, k;
728 
729  typedef Eigen::Triplet<double> T;
730  std::vector<T> tripletList;
731  tripletList.reserve(nchan);
732 
733  tripletList.clear();
734  tripletList.reserve(m_matProj.rows()*m_matProj.cols());
735  for(i = 0; i < m_matProj.rows(); ++i) {
736  for(k = 0; k < m_matProj.cols(); ++k) {
737  if(m_matProj(i,k) != 0) {
738  tripletList.push_back(T(i, k, m_matProj(i,k)));
739  }
740  }
741  }
742 
743  m_matSparseProjMult = SparseMatrix<double>(m_matProj.rows(),m_matProj.cols());
744  if(tripletList.size() > 0) {
745  m_matSparseProjMult.setFromTriplets(tripletList.begin(), tripletList.end());
746  }
747 
748  //Create full multiplication matrix
749  m_matSparseProjCompMult = m_matSparseProjMult * m_matSparseCompMult;
750  }
751 }
752 
753 //=============================================================================================================
754 
756 {
757  // Update the compensator
758  if(m_pFiffInfo) {
759  if(to == 0) {
760  m_bCompActivated = false;
761  } else {
762  m_bCompActivated = true;
763  }
764 
765 // qDebug()<<"to"<<to;
766 // qDebug()<<"from"<<from;
767 // qDebug()<<"m_bCompActivated"<<m_bCompActivated;
768 
769  FiffCtfComp newComp;
770  this->m_pFiffInfo->make_compensator(0, to, newComp);//Do this always from 0 since we always read new raw data, we never actually perform a multiplication on already existing data
771 
772  //We do not need to call this->m_pFiffInfo->set_current_comp(to);
773  //Because we will set the compensators to the coil in the same FiffInfo which is already used to write to file.
774  //Note that the data is written in raw form not in compensated form.
775  m_matComp = newComp.data->data;
776 
777  //
778  // Make proj sparse
779  //
780  qint32 nchan = this->m_pFiffInfo->nchan;
781  qint32 i, k;
782 
783  typedef Eigen::Triplet<double> T;
784  std::vector<T> tripletList;
785  tripletList.reserve(nchan);
786 
787  tripletList.clear();
788  tripletList.reserve(m_matComp.rows()*m_matComp.cols());
789  for(i = 0; i < m_matComp.rows(); ++i) {
790  for(k = 0; k < m_matComp.cols(); ++k) {
791  if(m_matComp(i,k) != 0) {
792  tripletList.push_back(T(i, k, m_matComp(i,k)));
793  }
794  }
795  }
796 
797  m_matSparseCompMult = SparseMatrix<double>(m_matComp.rows(),m_matComp.cols());
798  if(tripletList.size() > 0) {
799  m_matSparseCompMult.setFromTriplets(tripletList.begin(), tripletList.end());
800  }
801 
802  //Create full multiplication matrix
803  m_matSparseProjCompMult = m_matSparseProjMult * m_matSparseCompMult;
804  }
805 }
806 
807 //=============================================================================================================
808 
810 {
811  m_bSpharaActivated = state;
812 }
813 
814 //=============================================================================================================
815 
816 void RtFiffRawViewModel::updateSpharaOptions(const QString& sSytemType, int nBaseFctsFirst, int nBaseFctsSecond)
817 {
818  if(m_pFiffInfo) {
819  qDebug()<<"RtFiffRawViewModel::updateSpharaOptions - Creating SPHARA operator for"<<sSytemType;
820 
821  MatrixXd matSpharaMultFirst = MatrixXd::Identity(m_pFiffInfo->chs.size(), m_pFiffInfo->chs.size());
822  MatrixXd matSpharaMultSecond = MatrixXd::Identity(m_pFiffInfo->chs.size(), m_pFiffInfo->chs.size());
823 
824  if(sSytemType == "VectorView" && m_matSpharaVVGradLoaded.size() != 0 && m_matSpharaVVMagLoaded.size() != 0) {
825  matSpharaMultFirst = RTPROCESSINGLIB::makeSpharaProjector(m_matSpharaVVGradLoaded, m_vecIndicesFirstVV, m_pFiffInfo->nchan, nBaseFctsFirst, 1); //GRADIOMETERS
826  matSpharaMultSecond = RTPROCESSINGLIB::makeSpharaProjector(m_matSpharaVVMagLoaded, m_vecIndicesSecondVV, m_pFiffInfo->nchan, nBaseFctsSecond, 0); //Magnetometers
827  }
828 
829  if(sSytemType == "BabyMEG" && m_matSpharaBabyMEGInnerLoaded.size() != 0) {
830  matSpharaMultFirst = RTPROCESSINGLIB::makeSpharaProjector(m_matSpharaBabyMEGInnerLoaded, m_vecIndicesFirstBabyMEG, m_pFiffInfo->nchan, nBaseFctsFirst, 0); //InnerLayer
831  }
832 
833  if(sSytemType == "EEG" && m_matSpharaEEGLoaded.size() != 0) {
834  matSpharaMultFirst = RTPROCESSINGLIB::makeSpharaProjector(m_matSpharaEEGLoaded, m_vecIndicesFirstEEG, m_pFiffInfo->nchan, nBaseFctsFirst, 0); //InnerLayer
835  }
836 
837  //Write final operator matrices to file
838 // IOUtils::write_eigen_matrix(matSpharaMultFirst, QString(QCoreApplication::applicationDirPath() + "/resources/mne_scan/plugins/noisereduction/SPHARA/matSpharaMultFirst.txt"));
839 // IOUtils::write_eigen_matrix(matSpharaMultSecond, QString(QCoreApplication::applicationDirPath() + "/resources/mne_scan/plugins/noisereduction/SPHARA/matSpharaMultSecond.txt"));
840 // IOUtils::write_eigen_matrix(m_matSpharaEEGLoaded, QString(QCoreApplication::applicationDirPath() + "/resources/mne_scan/plugins/noisereduction/SPHARA/m_matSpharaEEGLoaded.txt"));
841 
842  //
843  // Make operators sparse
844  //
845  qint32 nchan = this->m_pFiffInfo->nchan;
846  qint32 i, k;
847 
848  typedef Eigen::Triplet<double> T;
849  std::vector<T> tripletList;
850  tripletList.reserve(nchan);
851 
852  //First operator
853  tripletList.clear();
854  tripletList.reserve(matSpharaMultFirst.rows()*matSpharaMultFirst.cols());
855  for(i = 0; i < matSpharaMultFirst.rows(); ++i) {
856  for(k = 0; k < matSpharaMultFirst.cols(); ++k) {
857  if(matSpharaMultFirst(i,k) != 0) {
858  tripletList.push_back(T(i, k, matSpharaMultFirst(i,k)));
859  }
860  }
861  }
862 
863  Eigen::SparseMatrix<double> matSparseSpharaMultFirst = SparseMatrix<double>(m_pFiffInfo->chs.size(),m_pFiffInfo->chs.size());
864 
865  matSparseSpharaMultFirst = SparseMatrix<double>(matSpharaMultFirst.rows(),matSpharaMultFirst.cols());
866  if(tripletList.size() > 0) {
867  matSparseSpharaMultFirst.setFromTriplets(tripletList.begin(), tripletList.end());
868  }
869 
870  //Second operator
871  tripletList.clear();
872  tripletList.reserve(matSpharaMultSecond.rows()*matSpharaMultSecond.cols());
873 
874  for(i = 0; i < matSpharaMultSecond.rows(); ++i) {
875  for(k = 0; k < matSpharaMultSecond.cols(); ++k) {
876  if(matSpharaMultSecond(i,k) != 0) {
877  tripletList.push_back(T(i, k, matSpharaMultSecond(i,k)));
878  }
879  }
880  }
881 
882  Eigen::SparseMatrix<double>matSparseSpharaMultSecond = SparseMatrix<double>(m_pFiffInfo->chs.size(),m_pFiffInfo->chs.size());
883 
884  if(tripletList.size() > 0) {
885  matSparseSpharaMultSecond.setFromTriplets(tripletList.begin(), tripletList.end());
886  }
887 
888  //Create full multiplication matrix
889  m_matSparseSpharaMult = matSparseSpharaMultFirst * matSparseSpharaMultSecond;
890  }
891 }
892 
893 //=============================================================================================================
894 
895 void RtFiffRawViewModel::setFilter(QList<FilterKernel> filterData)
896 {
897  m_filterKernel = filterData;
898 
899  m_iMaxFilterLength = 1;
900  for(int i=0; i<filterData.size(); ++i) {
901  if(m_iMaxFilterLength<filterData.at(i).getFilterOrder()) {
902  m_iMaxFilterLength = filterData.at(i).getFilterOrder();
903  }
904  }
905 
906  m_matOverlap.conservativeResize(m_pFiffInfo->chs.size(), m_iMaxFilterLength);
907  m_matOverlap.setZero();
908 
909  m_bDrawFilterFront = false;
910 
911  //Filter all visible data channels at once
912  //filterDataBlock();
913 }
914 
915 //=============================================================================================================
916 
918 {
919  m_bPerformFiltering = state;
920 }
921 
922 //=============================================================================================================
923 
924 void RtFiffRawViewModel::setBackgroundColor(const QColor& color)
925 {
926  m_colBackground = color;
927 }
928 
929 //=============================================================================================================
930 
931 void RtFiffRawViewModel::setFilterChannelType(const QString &channelType)
932 {
933  m_sFilterChannelType = channelType;
934  m_filterChannelList = m_visibleChannelList;
935 
936  //This version is for when all channels of a type are to be filtered (not only the visible ones).
937  //Create channel filter list independent from channelNames
938  m_filterChannelList.clear();
939 
940  for(int i = 0; i<m_pFiffInfo->chs.size(); ++i) {
941  if((m_pFiffInfo->chs.at(i).kind == FIFFV_MEG_CH || m_pFiffInfo->chs.at(i).kind == FIFFV_EEG_CH ||
942  m_pFiffInfo->chs.at(i).kind == FIFFV_EOG_CH || m_pFiffInfo->chs.at(i).kind == FIFFV_ECG_CH ||
943  m_pFiffInfo->chs.at(i).kind == FIFFV_EMG_CH)/* && !m_pFiffInfo->bads.contains(m_pFiffInfo->chs.at(i).ch_name)*/) {
944  if(m_sFilterChannelType == "All") {
945  m_filterChannelList << m_pFiffInfo->chs.at(i).ch_name;
946  } else if(m_pFiffInfo->chs.at(i).ch_name.contains(m_sFilterChannelType)) {
947  m_filterChannelList << m_pFiffInfo->chs.at(i).ch_name;
948  }
949  }
950  }
951 
952 // if(channelType != "All") {
953 // QMutableListIterator<QString> i(m_filterChannelList);
954 // while(i.hasNext()) {
955 // QString val = i.next();
956 // if(!val.contains(channelType, Qt::CaseInsensitive)) {
957 // i.remove();
958 // }
959 // }
960 // }
961 
962 // m_bDrawFilterFront = false;
963 
964  //Filter all visible data channels at once
965  //filterDataBlock();
966 }
967 
968 //=============================================================================================================
969 
970 void RtFiffRawViewModel::createFilterChannelList(QStringList channelNames)
971 {
972  m_filterChannelList.clear();
973  m_visibleChannelList = channelNames;
974 
975 // //Create channel fiter list based on channelNames
976 // for(int i = 0; i<m_pFiffInfo->chs.size(); ++i) {
977 // if((m_pFiffInfo->chs.at(i).kind == FIFFV_MEG_CH || m_pFiffInfo->chs.at(i).kind == FIFFV_EEG_CH ||
978 // m_pFiffInfo->chs.at(i).kind == FIFFV_EOG_CH || m_pFiffInfo->chs.at(i).kind == FIFFV_ECG_CH ||
979 // m_pFiffInfo->chs.at(i).kind == FIFFV_EMG_CH) && !m_pFiffInfo->bads.contains(m_pFiffInfo->chs.at(i).ch_name)) {
980 // if(m_sFilterChannelType == "All" && channelNames.contains(m_pFiffInfo->chs.at(i).ch_name))
981 // m_filterChannelList << m_pFiffInfo->chs.at(i).ch_name;
982 // else if(m_pFiffInfo->chs.at(i).ch_name.contains(m_sFilterChannelType) && channelNames.contains(m_pFiffInfo->chs.at(i).ch_name))
983 // m_filterChannelList << m_pFiffInfo->chs.at(i).ch_name;
984 // }
985 // }
986 
987  //Create channel filter list independent from channelNames
988  for(int i = 0; i < m_pFiffInfo->chs.size(); ++i) {
989  if((m_pFiffInfo->chs.at(i).kind == FIFFV_MEG_CH || m_pFiffInfo->chs.at(i).kind == FIFFV_EEG_CH ||
990  m_pFiffInfo->chs.at(i).kind == FIFFV_EOG_CH || m_pFiffInfo->chs.at(i).kind == FIFFV_ECG_CH ||
991  m_pFiffInfo->chs.at(i).kind == FIFFV_EMG_CH)/* && !m_pFiffInfo->bads.contains(m_pFiffInfo->chs.at(i).ch_name)*/) {
992  if(m_sFilterChannelType == "All") {
993  m_filterChannelList << m_pFiffInfo->chs.at(i).ch_name;
994  } else if(m_pFiffInfo->chs.at(i).ch_name.contains(m_sFilterChannelType)) {
995  m_filterChannelList << m_pFiffInfo->chs.at(i).ch_name;
996  }
997  }
998  }
999 
1000 // m_bDrawFilterFront = false;
1001 
1002 // for(int i = 0; i<m_filterChannelList.size(); ++i)
1003 // std::cout<<m_filterChannelList.at(i).toStdString()<<std::endl;
1004 
1005  //Filter all visible data channels at once
1006  //filterDataBlock();
1007 }
1008 
1009 //=============================================================================================================
1010 
1011 void RtFiffRawViewModel::markChBad(QModelIndex ch, bool status)
1012 {
1013  QList<FiffChInfo> chInfolist = m_pFiffInfo->chs;
1014 
1015  if(status) {
1016  if(!m_pFiffInfo->bads.contains(chInfolist[ch.row()].ch_name))
1017  m_pFiffInfo->bads.append(chInfolist[ch.row()].ch_name);
1018  qDebug() << "RawModel:" << chInfolist[ch.row()].ch_name << "marked as bad.";
1019  } else if(m_pFiffInfo->bads.contains(chInfolist[ch.row()].ch_name)) {
1020  int index = m_pFiffInfo->bads.indexOf(chInfolist[ch.row()].ch_name);
1021  m_pFiffInfo->bads.removeAt(index);
1022  qDebug() << "RawModel:" << chInfolist[ch.row()].ch_name << "marked as good.";
1023  }
1024 
1025  //Redefine channels which are to be filtered
1026  QStringList channelNames;
1027  createFilterChannelList(channelNames);
1028 
1029  //Update indeices of bad channels (this vector is needed when creating new ssp operators)
1030  QStringList emptyExclude;
1031  m_vecBadIdcs = FiffInfoBase::pick_channels(m_pFiffInfo->ch_names, m_pFiffInfo->bads, emptyExclude);
1032 
1033  emit dataChanged(ch,ch);
1034 }
1035 
1036 //=============================================================================================================
1037 
1038 void RtFiffRawViewModel::triggerInfoChanged(const QMap<double, QColor>& colorMap, bool active, QString triggerCh, double threshold)
1039 {
1040  m_qMapTriggerColor = colorMap;
1041  m_bTriggerDetectionActive = active;
1042  m_dTriggerThreshold = threshold;
1043 
1044  //Find channel index and initialise detected trigger map if channel name changed
1045  if(m_sCurrentTriggerCh != triggerCh) {
1046  m_sCurrentTriggerCh = triggerCh;
1047 
1048  QList<QPair<int,double> > temp;
1049  m_qMapDetectedTrigger.clear();
1050 
1051  for(int i = 0; i < m_pFiffInfo->chs.size(); ++i) {
1052  if(m_pFiffInfo->chs[i].ch_name == m_sCurrentTriggerCh) {
1053  m_iCurrentTriggerChIndex = i;
1054  m_qMapDetectedTrigger.insert(i, temp);
1055  break;
1056  }
1057  }
1058  }
1059 
1060  m_sCurrentTriggerCh = triggerCh;
1061 }
1062 
1063 //=============================================================================================================
1064 
1066 {
1067  if(value <= 0) {
1068  m_iDistanceTimerSpacer = 1000;
1069  } else {
1070  m_iDistanceTimerSpacer = value;
1071  }
1072 }
1073 
1074 //=============================================================================================================
1075 
1077 {
1078  m_iDetectedTriggers = 0;
1079 }
1080 
1081 //=============================================================================================================
1082 
1083 void RtFiffRawViewModel::markChBad(QModelIndexList chlist, bool status)
1084 {
1085  QList<FiffChInfo> chInfolist = m_pFiffInfo->chs;
1086 
1087  for(int i = 0; i < chlist.size(); ++i) {
1088  if(status) {
1089  if(!m_pFiffInfo->bads.contains(chInfolist[chlist[i].row()].ch_name))
1090  m_pFiffInfo->bads.append(chInfolist[chlist[i].row()].ch_name);
1091  } else {
1092  if(m_pFiffInfo->bads.contains(chInfolist[chlist[i].row()].ch_name)) {
1093  int index = m_pFiffInfo->bads.indexOf(chInfolist[chlist[i].row()].ch_name);
1094  m_pFiffInfo->bads.removeAt(index);
1095  }
1096  }
1097 
1098  emit dataChanged(chlist[i],chlist[i]);
1099  }
1100 
1101  //Update indeices of bad channels (this vector is needed when creating new ssp operators)
1102  QStringList emptyExclude;
1103  m_vecBadIdcs = FiffInfoBase::pick_channels(m_pFiffInfo->ch_names, m_pFiffInfo->bads, emptyExclude);
1104 }
1105 
1106 //=============================================================================================================
1107 
1108 void RtFiffRawViewModel::doFilterPerChannelRTMSA(QPair<QList<FilterKernel>,QPair<int,RowVectorXd> > &channelDataTime)
1109 {
1110  for(int i = 0; i < channelDataTime.first.size(); ++i) {
1111  //channelDataTime.second.second = channelDataTime.first.at(i).applyConvFilter(channelDataTime.second.second, true);
1112  channelDataTime.first[i].applyFftFilter(channelDataTime.second.second, true); //FFT Convolution for rt is not suitable. FFT make the signal filtering non causal.
1113  }
1114 }
1115 
1116 //=============================================================================================================
1117 
1118 void RtFiffRawViewModel::filterDataBlock()
1119 {
1120  //std::cout<<"START RtFiffRawViewModel::filterDataBlock"<<std::endl;
1121 
1122  if(m_filterKernel.isEmpty() || !m_bPerformFiltering) {
1123  return;
1124  }
1125 
1126  //Create temporary filters with higher fft length because we are going to filter all available data at once for one time
1127  QList<FilterKernel> tempFilterList;
1128 
1129  int fftLength = m_matDataRaw.row(0).cols() + 4 * m_iMaxFilterLength;
1130  int exp = ceil(MNEMath::log2(fftLength));
1131  fftLength = pow(2, exp) < 512 ? 512 : pow(2, exp);
1132 
1133  for(int i = 0; i<m_filterKernel.size(); ++i) {
1134  FilterKernel tempFilter(m_filterKernel.at(i).getName(),
1135  FilterKernel::m_filterTypes.indexOf(m_filterKernel.at(i).getFilterType()),
1136  m_filterKernel.at(i).getFilterOrder(),
1137  m_filterKernel.at(i).getCenterFrequency(),
1138  m_filterKernel.at(i).getBandwidth(),
1139  m_filterKernel.at(i).getParksWidth(),
1140  m_filterKernel.at(i).getSamplingFrequency(),
1141  FilterKernel::m_designMethods.indexOf(m_filterKernel.at(i).getDesignMethod()));
1142 
1143  tempFilterList.append(tempFilter);
1144  }
1145 
1146  //Generate QList structure which can be handled by the QConcurrent framework
1147  QList<QPair<QList<FilterKernel>,QPair<int,RowVectorXd> > > timeData;
1148  QList<int> notFilterChannelIndex;
1149 
1150  //Also append mirrored data in front and back to get rid of edge effects
1151  for(qint32 i=0; i<m_matDataRaw.rows(); ++i) {
1152  if(m_filterChannelList.contains(m_pFiffInfo->chs.at(i).ch_name)) {
1153  RowVectorXd datTemp(m_matDataRaw.row(i).cols() + 2 * m_iMaxFilterLength);
1154  datTemp << m_matDataRaw.row(i).head(m_iMaxFilterLength).reverse(), m_matDataRaw.row(i), m_matDataRaw.row(i).tail(m_iMaxFilterLength).reverse();
1155  timeData.append(QPair<QList<FilterKernel>,QPair<int,RowVectorXd> >(tempFilterList,QPair<int,RowVectorXd>(i,datTemp)));
1156  } else {
1157  notFilterChannelIndex.append(i);
1158  }
1159  }
1160 
1161  //Do the concurrent filtering
1162  if(!timeData.isEmpty()) {
1163  QFuture<void> future = QtConcurrent::map(timeData,
1164  doFilterPerChannelRTMSA);
1165 
1166  future.waitForFinished();
1167 
1168  for(int r = 0; r < timeData.size(); ++r) {
1169  m_matDataFiltered.row(timeData.at(r).second.first) = timeData.at(r).second.second.segment(m_iMaxFilterLength+m_iMaxFilterLength/2, m_matDataRaw.cols());
1170  m_matOverlap.row(timeData.at(r).second.first) = timeData.at(r).second.second.tail(m_iMaxFilterLength);
1171  }
1172  }
1173 
1174  //Fill filtered data with raw data if the channel was not filtered
1175  for(int i = 0; i < notFilterChannelIndex.size(); ++i) {
1176  m_matDataFiltered.row(notFilterChannelIndex.at(i)) = m_matDataRaw.row(notFilterChannelIndex.at(i));
1177  }
1178 
1179  if(!m_bIsFreezed) {
1180  m_vecLastBlockFirstValuesFiltered = m_matDataFiltered.col(0);
1181  }
1182 
1183  //std::cout<<"END RtFiffRawViewModel::filterDataBlock"<<std::endl;
1184 }
1185 
1186 //=============================================================================================================
1187 
1188 void RtFiffRawViewModel::filterDataBlock(const MatrixXd &data, int iDataIndex)
1189 {
1190  //std::cout<<"START RtFiffRawViewModel::filterDataBlock"<<std::endl;
1191 
1192  if(iDataIndex >= m_matDataFiltered.cols() || data.cols() < m_iMaxFilterLength) {
1193  return;
1194  }
1195 
1196  //Generate QList structure which can be handled by the QConcurrent framework
1197  QList<QPair<QList<FilterKernel>,QPair<int,RowVectorXd> > > timeData;
1198  QList<int> notFilterChannelIndex;
1199 
1200  for(qint32 i = 0; i < data.rows(); ++i) {
1201  if(m_filterChannelList.contains(m_pFiffInfo->chs.at(i).ch_name)) {
1202  timeData.append(QPair<QList<FilterKernel>,QPair<int,RowVectorXd> >(m_filterKernel,QPair<int,RowVectorXd>(i,data.row(i))));
1203  } else {
1204  notFilterChannelIndex.append(i);
1205  }
1206  }
1207 
1208  //Do the concurrent filtering
1209  if(!timeData.isEmpty()) {
1210  QFuture<void> future = QtConcurrent::map(timeData,
1211  doFilterPerChannelRTMSA);
1212 
1213  future.waitForFinished();
1214 
1215  //Do the overlap add method and store in m_matDataFiltered
1216  int iFilterDelay = m_iMaxFilterLength/2;
1217  int iFilteredNumberCols = timeData.at(0).second.second.cols();
1218 
1219  for(int r = 0; r<timeData.size(); ++r) {
1220  if(iDataIndex+2*data.cols() > m_matDataRaw.cols()) {
1221  //Handle last data block
1222  //std::cout<<"Handle last data block"<<std::endl;
1223 
1224  if(m_bDrawFilterFront) {
1225  //Get the currently filtered data. This data has a delay of filterLength/2 in front and back.
1226  RowVectorXd tempData = timeData.at(r).second.second;
1227 
1228  //Perform the actual overlap add by adding the last filterlength data to the newly filtered one
1229  tempData.head(m_iMaxFilterLength) += m_matOverlap.row(timeData.at(r).second.first);
1230 
1231  //Write the newly calulated filtered data to the filter data matrix. Keep in mind that the current block also effect last part of the last block (begin at dataIndex-iFilterDelay).
1232  int start = iDataIndex-iFilterDelay < 0 ? 0 : iDataIndex-iFilterDelay;
1233  m_matDataFiltered.row(timeData.at(r).second.first).segment(start,iFilteredNumberCols-m_iMaxFilterLength) = tempData.head(iFilteredNumberCols-m_iMaxFilterLength);
1234  } else {
1235  //Perform this else case everytime the filter was changed. Do not begin to plot from dataIndex-iFilterDelay because the impsulse response and m_matOverlap do not match with the new filter anymore.
1236  m_matDataFiltered.row(timeData.at(r).second.first).segment(iDataIndex-iFilterDelay,m_iMaxFilterLength) = timeData.at(r).second.second.segment(m_iMaxFilterLength,m_iMaxFilterLength);
1237  m_matDataFiltered.row(timeData.at(r).second.first).segment(iDataIndex+iFilterDelay,iFilteredNumberCols-2*m_iMaxFilterLength) = timeData.at(r).second.second.segment(m_iMaxFilterLength,iFilteredNumberCols-2*m_iMaxFilterLength);
1238  }
1239 
1240  //Refresh the m_matOverlap with the new calculated filtered data.
1241  m_matOverlap.row(timeData.at(r).second.first) = timeData.at(r).second.second.tail(m_iMaxFilterLength);
1242  } else if(iDataIndex == 0) {
1243  //Handle first data block
1244  //std::cout<<"Handle first data block"<<std::endl;
1245 
1246  if(m_bDrawFilterFront) {
1247  //Get the currently filtered data. This data has a delay of filterLength/2 in front and back.
1248  RowVectorXd tempData = timeData.at(r).second.second;
1249 
1250  //Add newly calculate data to the tail of the current filter data matrix
1251  m_matDataFiltered.row(timeData.at(r).second.first).segment(m_matDataFiltered.cols()-iFilterDelay-m_iResidual, iFilterDelay) = tempData.head(iFilterDelay) + m_matOverlap.row(timeData.at(r).second.first).head(iFilterDelay);
1252 
1253  //Perform the actual overlap add by adding the last filterlength data to the newly filtered one
1254  tempData.head(m_iMaxFilterLength) += m_matOverlap.row(timeData.at(r).second.first);
1255  m_matDataFiltered.row(timeData.at(r).second.first).head(iFilteredNumberCols-m_iMaxFilterLength-iFilterDelay) = tempData.segment(iFilterDelay,iFilteredNumberCols-m_iMaxFilterLength-iFilterDelay);
1256 
1257  //Copy residual data from the front to the back. The residual is != 0 if the chosen block size cannot be evenly fit into the matrix size
1258  m_matDataFiltered.row(timeData.at(r).second.first).tail(m_iResidual) = m_matDataFiltered.row(timeData.at(r).second.first).head(m_iResidual);
1259  } else {
1260  //Perform this else case everytime the filter was changed. Do not begin to plot from dataIndex-iFilterDelay because the impsulse response and m_matOverlap do not match with the new filter anymore.
1261  m_matDataFiltered.row(timeData.at(r).second.first).head(m_iMaxFilterLength) = timeData.at(r).second.second.segment(m_iMaxFilterLength,m_iMaxFilterLength);
1262  m_matDataFiltered.row(timeData.at(r).second.first).segment(iFilterDelay,iFilteredNumberCols-2*m_iMaxFilterLength) = timeData.at(r).second.second.segment(m_iMaxFilterLength,iFilteredNumberCols-2*m_iMaxFilterLength);
1263  }
1264 
1265  //Refresh the m_matOverlap with the new calculated filtered data.
1266  m_matOverlap.row(timeData.at(r).second.first) = timeData.at(r).second.second.tail(m_iMaxFilterLength);
1267  } else {
1268  //Handle middle data blocks
1269  //std::cout<<"Handle middle data block"<<std::endl;
1270 
1271  if(m_bDrawFilterFront) {
1272  //Get the currently filtered data. This data has a delay of filterLength/2 in front and back.
1273  RowVectorXd tempData = timeData.at(r).second.second;
1274 
1275  //Perform the actual overlap add by adding the last filterlength data to the newly filtered one
1276  tempData.head(m_iMaxFilterLength) += m_matOverlap.row(timeData.at(r).second.first);
1277 
1278  //Write the newly calulated filtered data to the filter data matrix. Keep in mind that the current block also effect last part of the last block (begin at dataIndex-iFilterDelay).
1279  m_matDataFiltered.row(timeData.at(r).second.first).segment(iDataIndex-iFilterDelay,iFilteredNumberCols-m_iMaxFilterLength) = tempData.head(iFilteredNumberCols-m_iMaxFilterLength);
1280  } else {
1281  //Perform this else case everytime the filter was changed. Do not begin to plot from dataIndex-iFilterDelay because the impsulse response and m_matOverlap do not match with the new filter anymore.
1282  m_matDataFiltered.row(timeData.at(r).second.first).segment(iDataIndex-iFilterDelay,m_iMaxFilterLength).setZero();// = timeData.at(r).second.second.segment(m_iMaxFilterLength,m_iMaxFilterLength);
1283  m_matDataFiltered.row(timeData.at(r).second.first).segment(iDataIndex+iFilterDelay,iFilteredNumberCols-2*m_iMaxFilterLength) = timeData.at(r).second.second.segment(m_iMaxFilterLength,iFilteredNumberCols-2*m_iMaxFilterLength);
1284  }
1285 
1286  //Refresh the m_matOverlap with the new calculated filtered data.
1287  m_matOverlap.row(timeData.at(r).second.first) = timeData.at(r).second.second.tail(m_iMaxFilterLength);
1288  }
1289  }
1290  }
1291 
1292  m_bDrawFilterFront = true;
1293 
1294  //Fill filtered data with raw data if the channel was not filtered
1295  for(int i = 0; i < notFilterChannelIndex.size(); ++i) {
1296  m_matDataFiltered.row(notFilterChannelIndex.at(i)).segment(iDataIndex,data.row(notFilterChannelIndex.at(i)).cols()) = data.row(notFilterChannelIndex.at(i));
1297  }
1298 
1299  //std::cout<<"END RtFiffRawViewModel::filterDataBlock"<<std::endl;
1300 }
1301 
1302 //=============================================================================================================
1303 
1304 void RtFiffRawViewModel::clearModel()
1305 {
1306  beginResetModel();
1307 
1308  m_matDataRaw.setZero();
1309  m_matDataFiltered.setZero();
1310  m_matDataRawFreeze.setZero();
1311  m_matDataFilteredFreeze.setZero();
1312  m_vecLastBlockFirstValuesFiltered.setZero();
1313  m_vecLastBlockFirstValuesRaw.setZero();
1314  m_matOverlap.setZero();
1315 
1316  endResetModel();
1317 }
1318 
1319 //=============================================================================================================
1320 
1322 {
1323  double dMaxValue;
1324  qint32 kind = getKind(row);
1325 
1326  switch(kind) {
1327  case FIFFV_MEG_CH: {
1328  dMaxValue = 1e-11f;
1329  qint32 unit = getUnit(row);
1330  if(unit == FIFF_UNIT_T_M) { //gradiometers
1331  dMaxValue = 1e-10f;
1332  if(getScaling().contains(FIFF_UNIT_T_M))
1333  dMaxValue = getScaling()[FIFF_UNIT_T_M];
1334  }
1335  else if(unit == FIFF_UNIT_T) //magnetometers
1336  {
1337  dMaxValue = 1e-11f;
1338  if(getScaling().contains(FIFF_UNIT_T))
1339  dMaxValue = getScaling()[FIFF_UNIT_T];
1340  }
1341  break;
1342  }
1343 
1344  case FIFFV_REF_MEG_CH: {
1345  dMaxValue = 1e-11f;
1346  if( getScaling().contains(FIFF_UNIT_T))
1347  dMaxValue = getScaling()[FIFF_UNIT_T];
1348  break;
1349  }
1350  case FIFFV_EEG_CH: {
1351  dMaxValue = 1e-4f;
1352  if( getScaling().contains(FIFFV_EEG_CH))
1353  dMaxValue = getScaling()[FIFFV_EEG_CH];
1354  break;
1355  }
1356  case FIFFV_EOG_CH: {
1357  dMaxValue = 1e-3f;
1358  if( getScaling().contains(FIFFV_EOG_CH))
1359  dMaxValue = getScaling()[FIFFV_EOG_CH];
1360  break;
1361  }
1362  case FIFFV_STIM_CH: {
1363  dMaxValue = 5;
1364  if( getScaling().contains(FIFFV_STIM_CH))
1365  dMaxValue = getScaling()[FIFFV_STIM_CH];
1366  break;
1367  }
1368  case FIFFV_MISC_CH: {
1369  dMaxValue = 1e-3f;
1370  if( getScaling().contains(FIFFV_MISC_CH))
1371  dMaxValue = getScaling()[FIFFV_MISC_CH];
1372  break;
1373  }
1374  default :
1375  dMaxValue = 1e-9f;
1376  break;
1377  }
1378 
1379  return dMaxValue;
1380 }
1381 
1382 //=============================================================================================================
1383 
1385 {
1386  auto pGroups = m_EventManager.getAllGroups();
1387  for(auto g : *pGroups){
1388  qDebug() << "Group: " << g.name.c_str() << "- Id: " << g.id;
1389  }
1390 
1391  m_EventManager.addEvent(iSample);
1392 
1393  auto pEvents = m_EventManager.getAllEvents();
1394 
1395  for(auto e : *pEvents){
1396  qDebug() << "Event> Sample: " << e.sample << "- Id: " << e.id;
1397  }
1398 }
1399 
1400 //=============================================================================================================
1401 
1402 std::unique_ptr<std::vector<EVENTSLIB::Event> > RtFiffRawViewModel::getEventsToDisplay(int iBegin, int iEnd) const
1403 {
1404  return m_EventManager.getEventsBetween(iBegin, iEnd);
1405 }
Old fiff_type declarations - replace them.
virtual int columnCount(const QModelIndex &parent=QModelIndex()) const
IOUtils class declaration.
void triggerDetected(int numberDetectedTriggers, const QMap< int, QList< QPair< int, double > > > &mapDetectedTriggers)
virtual QVariant data(const QModelIndex &index, int role=Qt::DisplayRole) const
FIFF measurement file information.
Definition: fiff_info.h:84
void toggleFreeze(const QModelIndex &index)
The FilterKernel class provides methods to create/design a FIR filter kernel.
Definition: filterkernel.h:132
void markChBad(QModelIndex ch, bool status)
DetectTrigger declarations.
void triggerInfoChanged(const QMap< double, QColor > &colorMap, bool active, QString triggerCh, double threshold)
QSharedPointer< RtFiffRawViewModel > SPtr
FIFFLIB::fiff_int_t getKind(qint32 row) const
Event addEvent(int sample)
void selectRows(const QList< qint32 > &selection)
CTF software compensation data.
Definition: fiff_ctf_comp.h:73
virtual int rowCount(const QModelIndex &parent=QModelIndex()) const
void setFiffInfo(QSharedPointer< FIFFLIB::FiffInfo > &p_pFiffInfo)
void setBackgroundColor(const QColor &color)
void updateProjection(const QList< FIFFLIB::FiffProj > &projs)
FIFFLIB::fiff_int_t getCoil(qint32 row) const
void setFilterChannelType(const QString &channelType)
FiffNamedMatrix::SDPtr data
Declaration of the RtFiffRawViewModel Class.
RTPROCESINGSHARED_EXPORT QList< QPair< int, double > > detectTriggerFlanksMax(const Eigen::MatrixXd &data, int iTriggerChannelIdx, int iOffsetIndex, double dThreshold, bool bRemoveOffset, int iBurstLengthSamp=100)
MNEMath class declaration.
void addData(const QList< Eigen::MatrixXd > &data)
void setSamplingInfo(float sps, int T, bool bSetZero=false)
void newSelection(const QList< qint32 > &selection)
std::unique_ptr< std::vector< Event > > getAllEvents() const
virtual QVariant headerData(int section, Qt::Orientation orientation, int role=Qt::DisplayRole) const
std::unique_ptr< std::vector< EVENTSLIB::Event > > getEventsToDisplay(int iBegin, int iEnd) const
void setScaling(const QMap< qint32, float > &p_qMapChScaling)
std::unique_ptr< std::vector< EventGroup > > getAllGroups() const
RtFiffRawViewModel(QObject *parent=0)
std::unique_ptr< std::vector< Event > > getEventsBetween(int sampleStart, int sampleEnd) const
void setFilter(QList< RTPROCESSINGLIB::FilterKernel > filterData)
void hideRows(const QList< qint32 > &selection)
double getMaxValueFromRawViewModel(int row) const
FIFFLIB::fiff_int_t getUnit(qint32 row) const
Declaration of the Sphara class.
#define FIFFV_REF_MEG_CH
#define FIFFV_COIL_NONE
void updateSpharaOptions(const QString &sSytemType, int nBaseFctsFirst, int nBaseFctsSecond)
FiffInfo class declaration.
const QMap< qint32, float > & getScaling() const
void createFilterChannelList(QStringList channelNames)
RTPROCESINGSHARED_EXPORT Eigen::MatrixXd makeSpharaProjector(const Eigen::MatrixXd &matBaseFct, const Eigen::VectorXi &vecIndices, int iOperatorDim, int iNBaseFct, int iSkip=0)