MNE-CPP  beta 1.0
rtsss.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "rtsss.h"
42 #include "rtsssalgo.h"
44 
45 //*************************************************************************************************************
46 //=============================================================================================================
47 // QT INCLUDES
48 //=============================================================================================================
49 
50 #include <QtCore/QtPlugin>
51 #include <QDebug>
52 #include <QFuture>
53 #include <QtConcurrent/QtConcurrentMap>
54 #include <QSettings>
55 
56 RtSssAlgo rsss;
57 
58 //*************************************************************************************************************
59 //=============================================================================================================
60 // USED NAMESPACES
61 //=============================================================================================================
62 
63 using namespace RtSssPlugin;
64 using namespace FIFFLIB;
65 using namespace MNEX;
66 using namespace XMEASLIB;
67 
68 
69 //*************************************************************************************************************
70 //=============================================================================================================
71 // DEFINE MEMBER METHODS
72 //=============================================================================================================
73 
74 RtSss::RtSss()
75 : m_bIsRunning(false)
76 , m_bReceiveData(false)
77 , m_bProcessData(false)
78 {
79 }
80 
81 
82 //*************************************************************************************************************
83 
84 RtSss::~RtSss()
85 {
86  if(this->isRunning())
87  stop();
88 }
89 
90 
91 //*************************************************************************************************************
92 
93 QSharedPointer<IPlugin> RtSss::clone() const
94 {
95  QSharedPointer<RtSss> pRtSssClone(new RtSss);
96  return pRtSssClone;
97 }
98 
99 
100 //*************************************************************************************************************
101 //=============================================================================================================
102 // Creating required display instances and set configurations
103 //=============================================================================================================
104 
105 void RtSss::init()
106 {
107  qDebug() << "*********** Initialization ************";
108 
109  //Delete Buffer - will be initailzed with first incoming data
110  if(!m_pRtSssBuffer.isNull())
111  m_pRtSssBuffer = CircularMatrixBuffer<double>::SPtr();
112 
113  // Input
114  m_pRTMSAInput = PluginInputData<NewRealTimeMultiSampleArray>::create(this, "RtSssIn", "RtSss input data");
115  connect(m_pRTMSAInput.data(), &PluginInputConnector::notify, this, &RtSss::update, Qt::DirectConnection);
116  m_inputConnectors.append(m_pRTMSAInput);
117 
118  // Output
119  m_pRTMSAOutput = PluginOutputData<NewRealTimeMultiSampleArray>::create(this, "RtSssOut", "RtSss output data");
120  m_outputConnectors.append(m_pRTMSAOutput);
121 
122 
123  //
124  // Load Settings
125  //
126 // QSettings settings;
127 // yourValue = settings.value(QString("Plugin/%1/yourValue").arg(this->getName()), 400).toInt();
128 }
129 
130 
131 //*************************************************************************************************************
132 
133 void RtSss::unload()
134 {
135  //
136  // Store Settings
137  //
138 // QSettings settings;
139 // settings.setValue(QString("Plugin/%1/yourValue").arg(this->getName()), yourValue);
140 
141 }
142 
143 
144 //*************************************************************************************************************
145 
146 bool RtSss::start()
147 {
148  qDebug() << "*********** Start ************";
149 
150  QThread::start();
151  return true;
152 }
153 
154 
155 //*************************************************************************************************************
156 
157 bool RtSss::stop()
158 {
159  qDebug() << "*********** Stop *************";
160 
161  m_bIsRunning = false;
162 
163  // Stop threads
164  QThread::terminate();
165  QThread::wait();
166 
167  if(m_pRtSssBuffer)
168  m_pRtSssBuffer->clear();
169 
170  m_bReceiveData = false;
171 
172  return true;
173 }
174 
175 
176 //*************************************************************************************************************
177 
178 IPlugin::PluginType RtSss::getType() const
179 {
180  return _IAlgorithm;
181 }
182 
183 
184 //*************************************************************************************************************
185 
186 QString RtSss::getName() const
187 {
188  return "RtSss Toolbox";
189 }
190 
191 
192 //*************************************************************************************************************
193 
194 QWidget* RtSss::setupWidget()
195 {
196  qDebug() << "*********** SetUpWidget ************";
197 
198  RtSssSetupWidget* widget = new RtSssSetupWidget(this); //widget is later distroyed by CentralWidget - so it has to be created everytime new
199 
200  connect(widget, &RtSssSetupWidget::signalNewLinRR, this, &RtSss::setLinRR);
201  connect(widget, &RtSssSetupWidget::signalNewLoutRR, this, &RtSss::setLoutRR);
202  connect(widget, &RtSssSetupWidget::signalNewLin, this, &RtSss::setLin);
203  connect(widget, &RtSssSetupWidget::signalNewLout, this, &RtSss::setLout);
204 
205  LinRR = widget->getLinRR();
206  LoutRR = widget->getLoutRR();
207  Lin = widget->getLin();
208  Lout = widget->getLout();
209 
210  return widget;
211 }
212 
213 
214 void RtSss::setLinRR(int val)
215 {
216 // qDebug() <<" new LinRR: " << val;
217  LinRR = val;
218 }
219 
220 void RtSss::setLoutRR(int val)
221 {
222 // qDebug() <<" new LoutRR: " << val;
223  LoutRR = val;
224 }
225 
226 void RtSss::setLin(int val)
227 {
228 // qDebug() <<" new Lin: " << val;
229  Lin = val;
230 }
231 
232 void RtSss::setLout(int val)
233 {
234 // qDebug() <<" new Lout: " << val;
235  Lout = val;
236 }
237 
238 //*************************************************************************************************************
239 
240 void RtSss::update(XMEASLIB::NewMeasurement::SPtr pMeasurement)
241 {
242 // qDebug() << "*********** Update ************";
243 
244  QSharedPointer<NewRealTimeMultiSampleArray> pRTMSA = pMeasurement.dynamicCast<NewRealTimeMultiSampleArray>();
245 
246  if(pRTMSA && m_bReceiveData)
247  {
248  //Check if buffer initialized
249  if(!m_pRtSssBuffer)
250  m_pRtSssBuffer = CircularMatrixBuffer<double>::SPtr(new CircularMatrixBuffer<double>(64, pRTMSA->getNumChannels(), pRTMSA->getMultiSampleArray()[0].cols()));
251 
252  //Fiff information
253  if(!m_pFiffInfo)
254  m_pFiffInfo = pRTMSA->info();
255 
256  if(m_bProcessData)
257  {
258  MatrixXd in_mat;
259  for(unsigned char i = 0; i < pRTMSA->getMultiArraySize(); ++i)
260  {
261  in_mat = pRTMSA->getMultiSampleArray()[i];
262  m_pRtSssBuffer->push(&in_mat);
263  }
264  }
265  }
266 }
267 
268 
269 //*************************************************************************************************************
270 
271 
272 MatrixXd rt_sss(const MatrixXd &p_mat)
273 {
274  MatrixXd out;
275  out = rsss.getSSSRR(p_mat);
276  return out;
277 }
278 
279 
280 //*************************************************************************************************************
281 
283 {
284 // QList<MatrixXd> lineqn;
285  MatrixXd lineqn;
286 
287  m_bIsRunning = true;
288 
289  // start receiving data
290  //
291  m_bReceiveData = true;
292 
293  // Read Fiff Info
294  //
295  while(!m_pFiffInfo)
296  msleep(10);// Wait for fiff Info
297 
298  // Initialize output
299  m_pRTMSAOutput->data()->initFromFiffInfo(m_pFiffInfo);
300  m_pRTMSAOutput->data()->setMultiArraySize(100);
301  // m_pRTMSAOutput->data()->setSamplingRate(m_pFiffInfo->sfreq);
302  m_pRTMSAOutput->data()->setVisibility(true);
303 
304  // Set MEG channel infomation to rtSSS
305  rsss.setMEGInfo(m_pFiffInfo);
306 
307  // Get channel information
308  qint32 nmegchan = rsss.getNumMEGChan();
309  qint32 nmegchanused = rsss.getNumMEGChanUsed();
310  qDebug() << "number of meg channels(run): " << nmegchan;
311  VectorXi badch = rsss.getBadChan();
312 
313  // Load and set the number of spherical harmonics expansion
314  qDebug() << "LinRR (run): " << LinRR << ", LoutRR (run): " << LoutRR <<", Lin (run): " << Lin <<", Lout (run): " << Lout;
315  QList<int> expOrder;
316  expOrder << LinRR << LoutRR << Lin << Lout;
317  rsss.setSSSParameter(expOrder);
318 
319 // // Find out if coils are all gradiometers, all magnetometers, or both.
320 // // When both gradiometers and magnetometers are used,
321 // // MagScale facor of 100 must be appiled to magnetomters.
322 // float MagScale;
323 // if ((0 < CoilGrad.sum()) && (CoilGrad.sum() < NumCoil)) MagScale = 100;
324 // else MagScale = 1;
325 
326 // VectorXd CoilScale;
327 // CoilScale.setOnes(NumCoil);
328 // for(int i=0; i<NumCoil; i++)
329 // {
330 // if (CoilGrad(i) == 0) CoilScale(i) = MagScale;
332 // }
333 
334  // Find a starting MEG channel index fiff
335  // When the MEG recording of the first channel is saved starting from the first row in the signal matrix, the startID_MEGch will be 0.
336  // When the MEG recording of the first channel is saved starting not from the first row and being followed by others like EEG,
337  // the startID_MEGch will be the first encounter of MEG channel other than zero.
338  //
339  qint32 startID_MEGch = 0;
340  for(qint32 i = 0; i < m_pFiffInfo->nchan; ++i)
341  if(m_pFiffInfo->chs[i].kind == FIFFV_MEG_CH)
342  {
343  startID_MEGch = i;
344  break;
345  }
346 // qDebug() << "strat id: " << startID_MEGch;
347 
348  // Build linear equation
349  qDebug() << "building SSS linear equation .....";
350  lineqn = rsss.buildLinearEqn();
351 
352  qDebug() << "..finished !!";
353 
354  // start processing data
355  m_bProcessData = true;
356  qint32 HEADMOV_COR_cnt = 15 ;
357  qint32 cnt=0;
358  qDebug() << "rtSSS started.....";
359 
360  bool m_bIsHeadMov = true;
361 
362  while(m_bIsRunning)
363  {
364 // if (m_bIsHeadMov)
365 // {
366  lineqn = rsss.buildLinearEqn();
367  qDebug() << "rebuilt SSS linear equation .....";
368 // m_bIsHeadMov = false;
369 // }
370 // else
371 // {
372 // m_bIsHeadMov = true;
373 // }
374 
375  qint16 nrows = m_pRtSssBuffer->rows();
376 
377  if(nrows > 0) // check if init
378  {
379  // * Dispatch the inputs * //
380  MatrixXd in_mat = m_pRtSssBuffer->pop();
381 // qDebug() << "size of in_mat (run): " << in_mat.rows() << " x " << in_mat.cols();
382 
383  // Remove bad channel signals
384  MatrixXd in_mat_used(nmegchanused, in_mat.cols());
385 // qDebug() << "size of in_mat_used (run): " << in_mat_used.rows() << " x " << in_mat_used.cols();
386  for(qint32 i = 0, k = 0; i < nmegchan; ++i)
387  if (badch(i) == 0)
388  {
389  in_mat_used.row(k) = in_mat.row(i);
390  k++;
391  }
392 // in_mat_used = in_mat.block(0,0,nmegchanused,in_mat.cols());
393 
394  // Implement Concurrent mapreduced for parallel processing
395  // divide the in_mat_used into 2 or 4 matrices, which renders 50ms or 25ms data
396  QList<MatrixXd> list_in_mat;
397  qint32 nSubSample = 20;
398  qint32 nThread = in_mat.cols() / nSubSample;
399 
400  for(qint32 ith = 0; ith < nThread; ++ ith)
401  list_in_mat.append(in_mat_used.block(0,ith*nSubSample,nmegchanused,nSubSample));
402 
403  QFuture<MatrixXd> res = QtConcurrent::mapped(list_in_mat, rt_sss);
404  res.waitForFinished();
405 
406  for(qint32 ith = 0; ith < nThread; ++ ith)
407  in_mat_used.block(0,ith*nSubSample,nmegchanused,nSubSample)= res.resultAt(ith);
408 
409  // Replace raw signal by SSS signal
410  for(qint32 i = 0, k = 0; i < nmegchan; ++i)
411  if (badch(i) == 0)
412  {
413  in_mat.row(i) = in_mat_used.row(k);
414  k++;
415  }
416 
417  // Output to display
418  for(qint32 i = 0; i <in_mat.cols(); ++i)
419  m_pRTMSAOutput->data()->setValue(1e7 * in_mat.col(i));
420 // m_pRTMSAOutput->data()->setValue(1e-16 * in_mat.col(i));
421 
422  cnt++;
423  qDebug() << cnt;
424  }
425  }
426 
427  m_bProcessData = false;
428  m_bReceiveData = false;
429  qDebug() << "rtSSS stopped.";
430 }
431 
432 
433 
434 //
435 // A possible way of passing parameters to QtConcurrent::mapped
436 
437 //QFuture<MatrixXd> thumbnails = QtConcurrent::mapped(list_in_mat, rsss.getSSSRR);
438 
439 
440 //class obj
441 //{
442 // MatrixXd data;
443 // int lineq1;
444 // int lineq2;
445 // int lineq3;
446 // int lineq4;
447 
448 // MatrixXd getSSSRR()
449 // {
450 
451 // lineq1;
452 // lineq2;
453 // lineq3;
454 // lineq4;
455 // MatrixXd res;
456 // return res;
457 // }
458 //};
459 
460 
461 //QList<obj> qListObj;
462 
464 //obj obj0;
465 //obj0.data = list_in_mat(0);
466 //obj0.lineq1 = 0; obj0.lineq2 = 0; obj0.lineq3 = 0; obj0.lineq4 = 0;
467 
468 //obj obj1;
469 //obj1.data = list_in_mat(1);
470 //obj1.lineq1 = 0; obj1.lineq2 = 0; obj1.lineq3 = 1; obj1.lineq4 = 0;
471 
472 //qListObj << obj0;
473 //qListObj << obj1;
475 
476 //QFuture<MatrixXd> thumbnails = QtConcurrent::mapped(qListObj, &obj::getSSSRR);
477 
478 
virtual void run()
Definition: rtsss.cpp:282
OutputConnectorList m_outputConnectors
Definition: IPlugin.h:222
Contains the declaration of the RtSssAlgo class.
The RtSss class provides a rtsss algorithm structure.
Definition: rtsss.h:102
QSharedPointer< CircularMatrixBuffer > SPtr
The circular matrix buffer.
Contains the declaration of the RtSssSetupWidget class.
void push(const Matrix< _Tp, Dynamic, Dynamic > *pMatrix)
The RtSssSetupWidget class provides the RtSss configuration window.
Definition: arrow.h:75
Matrix< _Tp, Dynamic, Dynamic > pop()
Definition: fiff.h:98
QSharedPointer< NewMeasurement > SPtr
Contains the declaration of the RTSSS class.
static QSharedPointer< PluginOutputData< T > > create(IPlugin *parent, const QString &name, const QString &descr)
InputConnectorList m_inputConnectors
Definition: IPlugin.h:221
static QSharedPointer< PluginInputData< T > > create(IPlugin *parent, const QString &name, const QString &descr)
The RealTimeMultiSampleArrayNew class is the base class of every RealTimeMultiSampleArrayNew Measurem...