MNE-CPP  beta 1.0
rthpi.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 #include <Eigen/Dense>
41 #include <iostream>
42 #include "rthpi.h"
44 #include <math.h>
45 #include <string>
46 #include <vector>
47 #include <algorithm>
48 #include <fstream> // std::ifstream
49 
50 //*************************************************************************************************************
51 //=============================================================================================================
52 // QT INCLUDES
53 //=============================================================================================================
54 
55 #include <QtCore/QtPlugin>
56 #include <QDebug>
57 
58 
59 //*************************************************************************************************************
60 //=============================================================================================================
61 // USED NAMESPACES
62 //=============================================================================================================
63 
64 using namespace RtHpiPlugin;
65 using namespace MNEX;
66 using namespace XMEASLIB;
67 
68 
69 //*************************************************************************************************************
70 //=============================================================================================================
71 // DEFINE MEMBER METHODS
72 //=============================================================================================================
73 
74 RtHpi::RtHpi()
75 : m_bIsRunning(false)
76 , m_bProcessData(false)
77 , m_pRTMSAInput(NULL)
78 , m_pRTMSAOutput(NULL)
79 , m_pRtHpiBuffer(CircularMatrixBuffer<double>::SPtr())
80 {
81 }
82 
83 
84 //*************************************************************************************************************
85 
86 RtHpi::~RtHpi()
87 {
88  if(this->isRunning())
89  stop();
90 }
91 
92 
93 //*************************************************************************************************************
94 
95 QSharedPointer<IPlugin> RtHpi::clone() const
96 {
97  QSharedPointer<RtHpi> pRtHpiClone(new RtHpi);
98  return pRtHpiClone;
99 }
100 
101 
102 //*************************************************************************************************************
103 //=============================================================================================================
104 // Creating required display instances and set configurations
105 //=============================================================================================================
106 
107 void RtHpi::init()
108 {
109  // Input
110  m_pRTMSAInput = PluginInputData<NewRealTimeMultiSampleArray>::create(this, "Rt HPI In", "RT HPI input data");
111  connect(m_pRTMSAInput.data(), &PluginInputConnector::notify, this, &RtHpi::update, Qt::DirectConnection);
112  m_inputConnectors.append(m_pRTMSAInput);
113 
114  // Output
115  m_pRTMSAOutput = PluginOutputData<NewRealTimeMultiSampleArray>::create(this, "Rt HPI Out", "RT HPI output data");
116  m_outputConnectors.append(m_pRTMSAOutput);
117 
118  m_pRTMSAOutput->data()->setMultiArraySize(100);
119  m_pRTMSAOutput->data()->setVisibility(true);
120 
121  //init channels when fiff info is available
122  connect(this, &RtHpi::fiffInfoAvailable, this, &RtHpi::initConnector);
123 
124  //Delete Buffer - will be initailzed with first incoming data
125  if(!m_pRtHpiBuffer.isNull())
126  m_pRtHpiBuffer = CircularMatrixBuffer<double>::SPtr();
127 }
128 
129 
130 //*************************************************************************************************************
131 
132 void RtHpi::unload()
133 {
134 
135 }
136 
137 
138 //*************************************************************************************************************
139 
140 void RtHpi::initConnector()
141 {
142  qDebug() << "void RtHpi::initConnector()";
143  if(m_pFiffInfo)
144  m_pRTMSAOutput->data()->initFromFiffInfo(m_pFiffInfo);
145 }
146 
147 
148 //*************************************************************************************************************
149 
150 bool RtHpi::start()
151 {
152  //Check if the thread is already or still running. This can happen if the start button is pressed immediately after the stop button was pressed. In this case the stopping process is not finished yet but the start process is initiated.
153  if(this->isRunning())
154  QThread::wait();
155 
156  m_bIsRunning = true;
157 
158  // Start threads
159  QThread::start();
160 
161  return true;
162 }
163 
164 
165 //*************************************************************************************************************
166 
167 bool RtHpi::stop()
168 {
169  //Wait until this thread is stopped
170  m_bIsRunning = false;
171 
172  if(m_bProcessData)
173  {
174  //In case the semaphore blocks the thread -> Release the QSemaphore and let it exit from the pop function (acquire statement)
175  m_pRtHpiBuffer->releaseFromPop();
176  m_pRtHpiBuffer->releaseFromPush();
177 
178  m_pRtHpiBuffer->clear();
179 
180  m_pRTMSAOutput->data()->clear();
181  }
182 
183  return true;
184 }
185 
186 
187 //*************************************************************************************************************
188 
189 IPlugin::PluginType RtHpi::getType() const
190 {
191  return _IAlgorithm;
192 }
193 
194 
195 //*************************************************************************************************************
196 
197 QString RtHpi::getName() const
198 {
199  return "RtHpi Toolbox";
200 }
201 
202 
203 //*************************************************************************************************************
204 
205 QWidget* RtHpi::setupWidget()
206 {
207  RtHpiSetupWidget* setupWidget = new RtHpiSetupWidget(this);//widget is later distroyed by CentralWidget - so it has to be created everytime new
208  return setupWidget;
209 }
210 
211 
212 //*************************************************************************************************************
213 
214 void RtHpi::update(XMEASLIB::NewMeasurement::SPtr pMeasurement)
215 {
216  QSharedPointer<NewRealTimeMultiSampleArray> pRTMSA = pMeasurement.dynamicCast<NewRealTimeMultiSampleArray>();
217 
218  if(pRTMSA)
219  {
220  //Check if buffer initialized
221  if(!m_pRtHpiBuffer)
222  m_pRtHpiBuffer = CircularMatrixBuffer<double>::SPtr(new CircularMatrixBuffer<double>(64, pRTMSA->getNumChannels(), pRTMSA->getMultiSampleArray()[0].cols()));
223 
224  //Fiff information
225  if(!m_pFiffInfo)
226  {
227  m_pFiffInfo = pRTMSA->info();
228  emit fiffInfoAvailable();
229  }
230 
231  if(m_bProcessData)
232  {
233  MatrixXd t_mat;
234 
235  for(qint32 i = 0; i < pRTMSA->getMultiArraySize(); ++i)
236  {
237  t_mat = pRTMSA->getMultiSampleArray()[i];
238  m_pRtHpiBuffer->push(&t_mat);
239  }
240  }
241  }
242 }
243 
244 
245 
246 //*************************************************************************************************************
247 
249 {
250  struct sens sensors;
251  struct coilParam coil;
252  int numCoils = 4;
253 
254  coil.pos = Eigen::MatrixXd::Zero(numCoils,3);
255  coil.mom = Eigen::MatrixXd::Zero(numCoils,3);
256 
257  int numCh = m_pFiffInfo->nchan;
258  //
259  // Read Fiff Info
260  //
261  while(!m_pFiffInfo) msleep(10);// Wait for fiff Info
262 
263  m_bProcessData = true;
264 
265  QVector<int> refind(0); // Reference Indices
266 
267  for(int j=0;j<numCh;j++) {
268  if (m_pFiffInfo->ch_names[j].indexOf("TRG013") >= 0) {
269  refind.append(j);
270  }
271  if (m_pFiffInfo->ch_names[j].indexOf("TRG014") >= 0) {
272  refind.append(j);
273  }
274  if (m_pFiffInfo->ch_names[j].indexOf("TRG015") >= 0) {
275  refind.append(j);
276  }
277  if (m_pFiffInfo->ch_names[j].indexOf("TRG016") >= 0) {
278  refind.append(j);
279  }
280  }
281 
282  QVector<int> innerind(0); // Inner layer channels
283 
284  for (int i=0;i<numCh;i++) {
285  int j = 0;
286  if(m_pFiffInfo->chs[i].coil_type == 7002) {
287  // Check if the sensor is bad, if not get the position
288  // and orientation
289  if(!(m_pFiffInfo->bads.contains(m_pFiffInfo->ch_names[i]))) {
290  innerind.append(j);
291  sensors.coilpos(j,0) = m_pFiffInfo->chs[i].loc(0);
292  sensors.coilpos(j,1) = m_pFiffInfo->chs[i].loc(1);
293  sensors.coilpos(j,2) = m_pFiffInfo->chs[i].loc(2);
294  sensors.coilori(j,0) = m_pFiffInfo->chs[i].loc(10);
295  sensors.coilori(j,1) = m_pFiffInfo->chs[i].loc(11);
296  sensors.coilori(j,2) = m_pFiffInfo->chs[i].loc(12);
297  j++;
298  }
299  }
300  }
301 
302  int samF = m_pFiffInfo->sfreq;
303 
304  int numLoc = 3, numBlock, samLoc, phase; // numLoc : Number of times to localize in a second
305  samLoc = samF/numLoc; // minimum samples required to localize numLoc times in a second
306 
307  QVector<MatrixXd> buffer;
308  Eigen::MatrixXd amp(innerind.size(),numCoils);
309 
310  while (m_bIsRunning) {
311  if(m_bProcessData) {
312  /* Dispatch the inputs */
313  MatrixXd t_mat = m_pRtHpiBuffer->pop();
314 
315  buffer.append(t_mat);
316 
317  if(buffer.size()*t_mat.cols() >= samLoc) {
318 
319  Eigen::MatrixXd alldata(t_mat.rows(),buffer.size()*t_mat.cols());
320  // Concatenate data into a matrix
321  for(int i=0;i<buffer.size();i++) alldata << buffer[i];
322 
323  // Get the data from inner layer channels
324  Eigen::MatrixXd innerdata(innerind.size(),samLoc);
325  Eigen::MatrixXd refdata(numCoils,samLoc);
326 
327  numBlock = alldata.cols()/samLoc;
328 
329  // Loop for localizing coils
330  for(int i = 0;i<numBlock;i++) {
331  for(int j = 0;j < innerind.size();j++)
332  innerdata.row(j) << alldata.block(innerind[j],i*samLoc,1,samLoc);
333  for(int j = 0;j < refind.size();j++)
334  refdata.row(j) << alldata.block(refind[j],i*samLoc,1,samLoc);
335 
336  // Once we have the required data, compute the amplitude and phase of the channels
337  for(int j = 0;j < numCoils;j++) {
338  for(int k = 0;k < innerind.size();k++) {
339  amp(k,j) = innerdata.row(k).cwiseProduct(refdata.row(j)).sum();
340  phase = amp(k,j)/abs(amp(k,j));
341  amp(k,j) = amp(k,j)*phase;
342  }
343  }
344 
345  coil = dipfit(coil, sensors, amp);
346  }
347 
348  buffer.clear();
349  }
350 
351  for(qint32 i = 0; i < t_mat.cols(); ++i)
352  m_pRTMSAOutput->data()->setValue(t_mat.col(i));
353  }
354  }
355 }
356 
357 /*********************************************************************************
358  * dipfit function is adapted from Fieldtrip Software. It has been
359  * heavily edited for use with MNE-X Software
360  *********************************************************************************/
361 
362 coilParam RtHpi::dipfit(struct coilParam coil, struct sens sensors, Eigen::MatrixXd data)
363 {
364  // Initialize variables
365  int display = 0;
366  int maxiter = 100;
367  dipError temp;
368 
369  coil = fminsearch(coil.pos, maxiter, 2 * maxiter * coil.pos.cols(), display, data, sensors);
370 
371  temp = dipfitError(coil.pos, data, sensors);
372 
373  coil.mom = temp.moment;
374 
375  return coil;
376 }
377 
378 /*********************************************************************************
379  * fminsearch Multidimensional unconstrained nonlinear minimization (Nelder-Mead).
380  * X = fminsearch(X0, maxiter, maxfun, display, data, sensors) starts at X0 and
381  * attempts to find a local minimizer
382  *********************************************************************************/
383 
384 coilParam RtHpi::fminsearch(Eigen::MatrixXd pos,int maxiter, int maxfun, int display, Eigen::MatrixXd data, struct sens sensors)
385 {
386  double tolx, tolf, rho, chi, psi, sigma, func_evals, usual_delta, zero_term_delta, temp1, temp2;
387  std::string header, how;
388  int n, itercount, prnt;
389  Eigen::MatrixXd onesn, two2np1, one2n, v, y, v1, tempX1, tempX2, xbar, xr, x, xe, xc, xcc, xin;
390  std::vector <double> fv, fv1;
391  std::vector <int> idx;
392  coilParam coil;
393  dipError tempdip, fxr, fxe, fxc, fxcc;
394 
395  tolx = tolf = 1e-4;
396 
397  switch(display)
398  {
399  case 0:
400  prnt = 0;
401  break;
402  default:
403  prnt = 1;
404  }
405 
406  header = " Iteration Func-count min f(x) Procedure";
407 
408  n = pos.cols();
409 
410  // Initialize parameters
411  rho = 1; chi = 2; psi = 0.5; sigma = 0.5;
412  onesn = Eigen::MatrixXd::Ones(1,n);
413  two2np1 = one2n = Eigen::MatrixXd::Zero(1,n);
414 
415  for(int i = 0;i < n;i++) {
416  two2np1(i) = 1 + i;
417  one2n(i) = i;
418  }
419 
420  v = v1 = Eigen::MatrixXd::Zero(n, n+1);
421  fv.resize(n+1);idx.resize(n+1);fv1.resize(n+1);
422 
423  for(int i = 0;i < n; i++) v(i,0) = pos(i);
424 
425  tempdip = dipfitError(pos, data, sensors);
426  fv[0] = tempdip.error;
427 
428 
429  func_evals = 1;itercount = 0;how = "";
430 
431  // Continue setting up the initial simplex.
432  // Following improvement suggested by L.Pfeffer at Stanford
433  usual_delta = 0.05; // 5 percent deltas for non-zero terms
434  zero_term_delta = 0.00025; // Even smaller delta for zero elements of x
435  xin = pos.transpose();
436 
437  for(int j = 0;j < n;j++) {
438  y = xin;
439 
440  if(y(j) != 0) y(j) = (1 + usual_delta) * y(j);
441  else y(j) = zero_term_delta;
442 
443  v.col(j+1).array() = y;
444  pos = y.transpose();
445  tempdip = dipfitError(pos, data, sensors);
446  fv[j+1] = tempdip.error;
447  }
448 
449  // Sort elements of fv
450  base_arr = fv;
451  for (int i = 0; i < n+1; i++) idx[i] = i;
452 
453  sort (idx.begin(), idx.end(), RtHpi::compar);
454 
455  for (int i = 0;i < n+1;i++) {
456  v1.col(i) = v.col(idx[i]);
457  fv1[i] = fv[idx[i]];
458  }
459 
460  v = v1;fv = fv1;
461 
462  how = "initial simplex";
463  itercount = itercount + 1;
464  func_evals = n + 1;
465 
466  tempX1 = Eigen::MatrixXd::Zero(1,n);
467 
468  while ((func_evals < maxfun) && (itercount < maxiter)) {
469  for (int i = 0;i < n;i++) tempX1(i) = abs(fv[0] - fv[i+1]);
470  temp1 = tempX1.maxCoeff();
471 
472  tempX2 = Eigen::MatrixXd::Zero(n,n);
473 
474  for(int i = 0;i < n;i++) tempX2.col(i) = v.col(i+1) - v.col(0);
475 
476  tempX2 = tempX2.array().abs();
477 
478  temp2 = tempX2.maxCoeff();
479 
480  if((temp1 <= tolf) && (temp2 <= tolx)) break;
481 
482  xbar = v.block(0,0,n,n).rowwise().sum();
483  xbar /= n;
484 
485  xr = (1+rho) * xbar - rho * v.block(0,n,v.rows(),1);
486 
487  x = xr.transpose();
488  std::cout << "Iteration Count: " << itercount << ":" << x << std::endl;
489 
490  fxr = dipfitError(x, data, sensors);
491 
492  func_evals = func_evals+1;
493 
494  if (fxr.error < fv[0]) {
495  // Calculate the expansion point
496  xe = (1 + rho * chi) * xbar - rho * chi * v.col(v.cols()-1);
497  x = xe.transpose();
498  fxe = dipfitError(x, data, sensors);
499  func_evals = func_evals+1;
500 
501  if(fxe.error < fxr.error) {
502  v.col(v.cols()-1) = xe;
503  fv[n] = fxe.error;
504  how = "expand";
505  }
506  else {
507  v.col(v.cols()-1) = xr;
508  fv[n] = fxr.error;
509  how = "reflect";
510  }
511  }
512  else {
513  if(fxr.error < fv[n-1]) {
514  v.col(v.cols()-1) = xr;
515  fv[n] = fxr.error;
516  how = "reflect";
517  }
518  else { // fxr.error >= fv[:,n-1]
519  // Perform contraction
520  if(fxr.error < fv[n]) {
521  // Perform an outside contraction
522  xc = (1 + psi * rho) * xbar - psi * rho * v.col(v.cols()-1);
523  x = xc.transpose();
524  fxc = dipfitError(x, data, sensors);
525  func_evals = func_evals + 1;
526 
527  if(fxc.error <= fxr.error) {
528  v.col(v.cols()-1) = xc;
529  fv[n] = fxc.error;
530  how = "contract outside";
531  }
532  else {
533  // perform a shrink
534  how = "shrink";
535  }
536  }
537  else {
538  xcc = (1 - psi) * xbar + psi * v.col(v.cols()-1);
539  x = xcc.transpose();
540  fxcc = dipfitError(x, data, sensors);
541  func_evals = func_evals+1;
542  if(fxcc.error < fv[n]) {
543  v.col(v.cols()-1) = xcc;
544  fv[n] = fxcc.error;
545  how = "contract inside";
546  }
547  else {
548  // perform a shrink
549  how = "shrink";
550  }
551  }
552  if(how.compare("shrink") == 0) {
553  for(int j = 1;j < n+1;j++) {
554  v.col(j).array() = v.col(0).array() + sigma * (v.col(j).array() - v.col(0).array());
555  x = v.col(j).array().transpose();
556  tempdip = dipfitError(x,data, sensors);
557  fv[j] = tempdip.error;
558  }
559  }
560  }
561  }
562  // Sort elements of fv
563  base_arr = fv;
564  for (int i = 0; i < n+1; i++) idx[i] = i;
565  sort (idx.begin (), idx.end (), RtHpi::compar);
566  for (int i = 0;i < n+1;i++) {
567  v1.col(i) = v.col(idx[i]);
568  fv1[i] = fv[idx[i]];
569  }
570  v = v1;fv = fv1;
571  itercount = itercount + 1;
572  }
573 
574  x = v.col(0).transpose();
575  std::cout << x << std::endl;
576  coil.pos = x;
577  return coil;
578 }
579 
580 
581 /*********************************************************************************
582  * dipfitError computes the error between measured and model data
583  * and can be used for non-linear fitting of dipole position.
584  * The function has been compared with matlab dipfit_error and it gives
585  * same output
586  *********************************************************************************/
587 
588 dipError RtHpi::dipfitError(Eigen::MatrixXd pos, Eigen::MatrixXd data, struct sens sensors)
589 {
590  // Variable Declaration
591  struct dipError e;
592  Eigen::MatrixXd lf, dif;
593 
594  // Compute lead field for a magnetic dipole in infinite vacuum
595  lf = ft_compute_leadfield(pos, sensors);
596 
597  //e.moment = pinv(lf)*data;
598  e.moment = pinv(lf) * data;
599 
600  dif = data - lf * e.moment;
601 
602  e.error = dif.array().square().sum()/data.array().square().sum();
603 
604  return e;
605 }
606 
607 /*********************************************************************************
608  * ft_compute_leadfield computes a forward solution for a dipole in a a volume
609  * conductor model. The forward solution is expressed as the leadfield
610  * matrix (Nchan*3), where each column corresponds with the potential or field
611  * distributions on all sensors for one of the x,y,z-orientations of the dipole.
612  * The function has been compared with matlab ft_compute_leadfield and it gives
613  * same output
614  *********************************************************************************/
615 
616 Eigen::MatrixXd RtHpi::ft_compute_leadfield(Eigen::MatrixXd pos, struct sens sensors)
617 {
618 
619  Eigen::MatrixXd pnt, ori, lf;
620 
621  pnt = sensors.coilpos; // position of each coil
622  ori = sensors.coilori; // orientation of each coil
623 
624  lf = magnetic_dipole(pos, pnt, ori);
625 
626  lf = sensors.tra * lf;
627 
628  return lf;
629 }
630 
631 /*********************************************************************************
632  * magnetic_dipole leadfield for a magnetic dipole in an infinite medium
633  * The function has been compared with matlab magnetic_dipole and it gives same output
634  *********************************************************************************/
635 
636 Eigen::MatrixXd RtHpi::magnetic_dipole(Eigen::MatrixXd pos, Eigen::MatrixXd pnt, Eigen::MatrixXd ori) {
637 
638  double u0 = 1e-7;
639  int nchan;
640  Eigen::MatrixXd r, r2, r5, x, y, z, mx, my, mz, Tx, Ty, Tz, lf, temp;
641 
642  nchan = pnt.rows();
643 
644  // Shift the magnetometers so that the dipole is in the origin
645  pnt.array().col(0) -= pos(0);pnt.array().col(1) -= pos(1);pnt.array().col(2) -= pos(2);
646 
647  r = pnt.array().square().rowwise().sum().sqrt();
648 
649  r2 = r5 = x = y = z = mx = my = mz = Tx = Ty = Tz = lf = Eigen::MatrixXd::Zero(nchan,3);
650 
651  for(int i = 0;i < nchan;i++) {
652  r2.row(i).array().fill(pow(r(i),2));
653  r5.row(i).array().fill(pow(r(i),5));
654  }
655 
656  for(int i = 0;i < nchan;i++) {
657  x.row(i).array().fill(pnt(i,0));
658  y.row(i).array().fill(pnt(i,1));
659  z.row(i).array().fill(pnt(i,2));
660  }
661  mx.col(0).array().fill(1);my.col(1).array().fill(1);mz.col(2).array().fill(1);
662 
663  Tx = 3 * x.cwiseProduct(pnt) - mx.cwiseProduct(r2);
664  Ty = 3 * y.cwiseProduct(pnt) - my.cwiseProduct(r2);
665  Tz = 3 * z.cwiseProduct(pnt) - mz.cwiseProduct(r2);
666 
667  for(int i = 0;i < nchan;i++) {
668  lf(i,0) = Tx.row(i).dot(ori.row(i));
669  lf(i,1) = Ty.row(i).dot(ori.row(i));
670  lf(i,2) = Tz.row(i).dot(ori.row(i));
671  }
672 
673  for(int i = 0;i < nchan;i++) {
674  for(int j = 0;j < 3;j++) {
675  lf(i,j) = u0 * lf(i,j)/(4 * M_PI * r5(i,j));
676  }
677  }
678  return lf;
679 }
680 
681 bool RtHpi::compar (int a, int b){
682  return ((base_arr)[a] < (base_arr)[b]);
683 }
684 
685 Eigen::Matrix4d computeTransformation(Eigen::MatrixXd NH, Eigen::MatrixXd BT)
686 {
687  Eigen::MatrixXd xdiff, ydiff, zdiff, C, Q;
688  Eigen::Matrix4d trans = Eigen::Matrix4d::Identity(4,4), Rot = Eigen::Matrix4d::Zero(4,4), Trans = Eigen::Matrix4d::Identity(4,4);
689  double meanx,meany,meanz,normf;
690 
691  for(int i = 0;i < 50;i++) {
692  zdiff = NH.col(2) - BT.col(2);
693  ydiff = NH.col(1) - BT.col(1);
694  xdiff = NH.col(0) - BT.col(0);
695 
696  meanx=xdiff.mean();
697  meany=ydiff.mean();
698  meanz=zdiff.mean();
699 
700  for (int j=0;j<NH.rows();j++) {
701  BT(j,0) = BT(j,0) + meanx;
702  BT(j,1) = BT(j,1) + meany;
703  BT(j,2) = BT(j,2) + meanz;
704  }
705 
706  C = BT.transpose() * NH;
707 
708  Eigen::JacobiSVD< Eigen::MatrixXd > svd(C ,Eigen::ComputeThinU | Eigen::ComputeThinV);
709 
710  Q = svd.matrixU() * svd.matrixV().transpose();
711 
712  BT = BT * Q;
713 
714  normf = (NH.transpose()-BT.transpose()).norm();
715 
716  for(int j=0;j<3;j++) {
717  for(int k=0;k<3;k++) Rot(j,k) = Q(k,j);
718  }
719 
720  Rot(3,3) = 1;
721  Trans(0,3) = meanx;Trans(1,3) = meany;Trans(2,3) = meanz;
722  trans = Rot * Trans * trans;
723  }
724  return trans;
725 }
726 
727 Eigen::MatrixXd RtHpi::pinv(Eigen::MatrixXd a)
728 {
729  double epsilon = std::numeric_limits<double>::epsilon();
730  Eigen::JacobiSVD< Eigen::MatrixXd > svd(a ,Eigen::ComputeThinU | Eigen::ComputeThinV);
731  double tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs()(0);
732  return svd.matrixV() * (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(),0).matrix().asDiagonal() * svd.matrixU().adjoint();
733 }
734 std::vector <double>RtHpi::base_arr;
OutputConnectorList m_outputConnectors
Definition: IPlugin.h:222
The RtHpiSetupWidget class provides the RtHpi configuration window.
QSharedPointer< CircularMatrixBuffer > SPtr
The circular matrix buffer.
The RtHpi class provides a RtHpi algorithm structure.
Definition: rthpi.h:117
Contains the declaration of the RTHPI class.
virtual void run()
Definition: rthpi.cpp:248
void push(const Matrix< _Tp, Dynamic, Dynamic > *pMatrix)
Contains the declaration of the RtHpiSetupWidget class.
Definition: arrow.h:75
Matrix< _Tp, Dynamic, Dynamic > pop()
QSharedPointer< NewMeasurement > SPtr
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...