MNE-CPP 0.1.9
A Framework for Electrophysiology
Loading...
Searching...
No Matches
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
71using namespace DISPLIB;
72using namespace UTILSLIB;
73using namespace FIFFLIB;
74using namespace Eigen;
75using namespace RTPROCESSINGLIB;
76
77//=============================================================================================================
78// DEFINE MEMBER METHODS
79//=============================================================================================================
80
81EVENTSLIB::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
119
120//=============================================================================================================
121//virtual functions
122int 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
133int RtFiffRawViewModel::columnCount(const QModelIndex & /*parent*/) const
134{
135 return 3;
136}
137
138//=============================================================================================================
139
140QVariant 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
206QVariant 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
238void 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
298void 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
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
375void 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
417void 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
574fiff_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
586fiff_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
598fiff_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
610void 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
631void 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
663void 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
686void RtFiffRawViewModel::setScaling(const QMap< qint32,float >& p_qMapChScaling)
687{
688 beginResetModel();
689 m_qMapChScaling = p_qMapChScaling;
690 endResetModel();
691}
692
693//=============================================================================================================
694
695void 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
816void 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
895void 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
925{
926 m_colBackground = color;
927}
928
929//=============================================================================================================
930
931void 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
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
1011void 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
1038void 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
1083void 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
1108void 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
1118void 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
1188void 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
1304void 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 m_EventManager.addEvent(iSample);
1387}
1388
1389//=============================================================================================================
1390
1391std::unique_ptr<std::vector<EVENTSLIB::Event> > RtFiffRawViewModel::getEventsToDisplay(int iBegin, int iEnd) const
1392{
1393 return m_EventManager.getEventsBetween(iBegin, iEnd);
1394}
Declaration of the Sphara class.
RTPROCESINGSHARED_EXPORT Eigen::MatrixXd makeSpharaProjector(const Eigen::MatrixXd &matBaseFct, const Eigen::VectorXi &vecIndices, int iOperatorDim, int iNBaseFct, int iSkip=0)
RTPROCESINGSHARED_EXPORT Eigen::MatrixXd filterData(const Eigen::MatrixXd &matData, int type, double dCenterfreq, double dBandwidth, double dTransition, double dSFreq, int iOrder=1024, int designMethod=FilterKernel::m_designMethods.indexOf(FilterParameter("Cosine")), const Eigen::RowVectorXi &vecPicks=Eigen::RowVectorXi(), bool bUseThreads=true, bool bKeepOverhead=false)
DetectTrigger declarations.
RTPROCESINGSHARED_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)
Declaration of the RtFiffRawViewModel Class.
FiffInfo class declaration.
#define FIFFV_REF_MEG_CH
int k
Definition fiff_tag.cpp:324
Definitions for describing the objects in a FIFF file.
IOUtils class declaration.
MNEMath class declaration.
const QMap< qint32, float > & getScaling() const
void setBackgroundColor(const QColor &color)
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)
std::unique_ptr< std::vector< EVENTSLIB::Event > > getEventsToDisplay(int iBegin, int iEnd) const
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 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
void triggerInfoChanged(const QMap< double, QColor > &colorMap, bool active, QString triggerCh, double threshold)
void setFilter(QList< RTPROCESSINGLIB::FilterKernel > filterData)
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
Event addEvent(int sample)
std::unique_ptr< std::vector< Event > > getEventsBetween(int sampleStart, int sampleEnd) const
CTF software compensation data.
FiffNamedMatrix::SDPtr data
FIFF measurement file information.
Definition fiff_info.h:85
static Eigen::RowVectorXi pick_channels(const QStringList &ch_names, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList)
The FilterKernel class provides methods to create/design a FIR filter kernel.
static QVector< FilterParameter > m_designMethods
static QVector< FilterParameter > m_filterTypes
static bool read_eigen_matrix(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &out, const QString &path)
Definition ioutils.h:481
static double log2(const T d)
Definition mnemath.h:607