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