MNE-CPP  beta 1.0
rtsssalgo.cpp
Go to the documentation of this file.
1 //=============================================================================================================
39 #include "rtsssalgo.h"
40 #include <QFuture>
41 #include <QtConcurrent/QtConcurrentMap>
42 //#include "FormFiles/rtssssetupwidget.h"
43 
44 RtSssAlgo::RtSssAlgo()
45 {
46 
47 }
48 
49 RtSssAlgo::~RtSssAlgo()
50 {
51 
52 }
53 
54 //QList<MatrixXd> RtSssAlgo::buildLinearEqn()
55 MatrixXd RtSssAlgo::buildLinearEqn()
56 {
57  QList<MatrixXd> Eqn, EqnRR;
58  QList<MatrixXd> LinEqn;
59  qint32 LIn, LOut;
60 // MatrixXd EqnInRR, EqnOutRR, EqnIn, EqnOut;
61 
62 // int MagScale;
63 // int MACHINE_TYPE = VECTORVIEW;
64 
65 // if (MACHINE_TYPE == VECTORVIEW)
66 // {
67 // MagScale = 100;
68 // std::cout << "loading coil information (VectorView).....";
69 // getCoilInfoVectorView();
70 // std::cout << "loading coil information (VectorView)..... ** finished **" << endl;
71 // }
72 // else if (MACHINE_TYPE == BABYMEG)
73 // {
74 // MagScale = 1;
75 // std::cout << "loading coil information (BabyMEG).....";
76 // getCoilInfoBabyMEG4Sim();
77 // std::cout << "loading coil information (BabyMEG)..... ** finisehd **" << endl;
78 // }
79 
80 
81 // Compute SSS equation
82  LIn = LInRR;
83  LOut = LOutRR;
84  EqnRR = getSSSEqn(LIn, LOut);
85  EqnInRR = EqnRR[0];
86  EqnOutRR = EqnRR[1];
87 
88  LIn = LInOLS;
89  LOut = LOutOLS;
90  Eqn = getSSSEqn(LIn, LOut);
91  EqnIn = Eqn[0];
92  EqnOut = Eqn[1];
93 
94 //
95 // Vector2i LexpRR, LexpOLS;
96 // LexpRR(0) =LInRR;
97 // LexpRR(1) = LOutRR;
98 
99 // LexpOLS(0) =LInOLS;
100 // LexpOLS(1) = LOutOLS;
101 
102 // QList<Vector2i> Lexp;
103 // Lexp.append(LexpRR);
104 // Lexp.append(LexpOLS);
105 
106 // QFuture<QList> res = QtConcurrent::mapped(Lexp,getSSSEqn);
107 // res.waitForFinished();
108 
109 
110 // build linear equation
111 // MatrixXd EqnARR(NumCoil, EqnInRR.cols()+EqnOutRR.cols());
112 // MatrixXd EqnA(NumCoil, EqnIn.cols()+EqnOut.cols());
113 // MatrixXd EqnB(NumCoil,1);
114  EqnARR.resize(NumCoil, EqnInRR.cols()+EqnOutRR.cols());
115  EqnA.resize(NumCoil, EqnIn.cols()+EqnOut.cols());
116  EqnB.resize(NumCoil,1);
117 
118  // Find out if coils are all gradiometers, all magnetometers, or both.
119  // When both gradiometers and magnetometers are used,
120  // MagScale facor of 100 must be appiled to magnetomters.
121  float MagScale;
122  if ((0 < CoilGrad.sum()) && (CoilGrad.sum() < NumCoil)) MagScale = 100;
123  else MagScale = 1;
124 
125  VectorXd CoilScale;
126  CoilScale.setOnes(NumCoil);
127  for(int i=0; i<NumCoil; i++)
128  {
129  if (CoilGrad(i) == 0) CoilScale(i) = MagScale;
130 // std::cout << "i=" << i << "CoilGrad: " << CoilGrad(i) << ", CoilScale: " << CoilScale(i) << std::endl;
131  }
132 
133  EqnARR << EqnInRR, EqnOutRR;
134  EqnA << EqnIn, EqnOut;
135 
136  EqnARR = CoilScale.asDiagonal() * EqnARR;
137  EqnA = CoilScale.asDiagonal() * EqnA;
138 // std::cout << "pass 1" << std::endl;
139 // std::cout << "MEGData: " << MEGData.rows() << " x " << MEGData.cols() << std::endl;
140 // EqnB = CoilScale.asDiagonal() * MEGData;
141 // std::cout << "EqnB: " << EqnB.rows() << " x " << EqnB.cols() << std::endl;
142 // std::cout << "pass 2" << std::endl;
143 
144 // LinEqn.append(EqnInRR);
145 // LinEqn.append(EqnOutRR);
146  LinEqn.append(EqnIn);
147  LinEqn.append(EqnOut);
148  LinEqn.append(EqnARR);
149  LinEqn.append(EqnA);
150 // LinEqn.append(EqnB);
151  LinEqn.append(CoilScale.asDiagonal());
152 
153 
154 // std::cout << "EqnInRR *********************************" << endl << EqnInRR << endl << endl;
155 // std::cout << "EqnOutRR ********************************" << endl << EqnOutRR << endl << endl;
156 // std::cout << "EqnIn ***********************************" << endl << EqnIn << endl << endl;
157 // std::cout << "EqnOut **********************************" << endl << EqnOut << endl << endl;
158 // std::cout << "EqnARR **********************************" << endl << EqnARR << endl << endl;
159 // std::cout << "EqnA ************************************" << endl << EqnA << endl << endl;
160 // std::cout << "EqnB ************************************" << endl << EqnB.transpose() << endl;
161 
162 // std::cout << "building SSS linear equation .....finished !" << endl;
163 
164 // return LinEqn;
165  return CoilScale.asDiagonal();
166 }
167 
168 void RtSssAlgo::setSSSParameter(QList<int> expansionOrder)
169 {
170 // LInRR = 5;
171 // LOutRR = 4;
172 // LInOLS = 8;
173 // LOutOLS = 4;
174 
175  LInRR = expansionOrder[0];
176  LOutRR = expansionOrder[1];
177  LInOLS = expansionOrder[2];
178  LOutOLS = expansionOrder[3];
179 }
180 
181 void RtSssAlgo::setMEGInfo(FiffInfo::SPtr fiffInfo)
182 {
183 
184  // Set origin of head(?) coordinate
185  Origin.resize(3);
186  Origin << 0.0, 0.0, 0.04;
187 
188 // // Find the number of MEG channels
189 // qint32 nmegchan = 0;
190 // for (qint32 i=0; i<fiffInfo->nchan; ++i)
191 // if(fiffInfo->chs[i].kind == FIFFV_MEG_CH)
192 // nmegchan ++;
193 
194  // Find the number of MEG channels
195  NumMEGChan = 0,
196  NumBadCoil = 0;
197  BadChan.setZero(fiffInfo->nchan);
198  for (qint32 i=0; i<fiffInfo->nchan; ++i)
199  if(fiffInfo->chs[i].kind == FIFFV_MEG_CH)
200  {
201  for(qint32 j = 0; j < fiffInfo->bads.size(); j++)
202  if(fiffInfo->chs[i].ch_name == fiffInfo->bads[j])
203  {
204  BadChan(i) = 1;
205  qDebug() << "bad ch: " << i << "(" << fiffInfo->chs[i].ch_name << ")";
206  NumBadCoil ++;
207 // nmegbadchan ++;
208  }
209  NumMEGChan ++;
210  }
211 
212  NumCoil = NumMEGChan - NumBadCoil;
213 
214  qDebug() << "number of meg channels: " << NumMEGChan;
215  qDebug() << "number of bad meg channels : " << NumBadCoil;
216  qDebug() << "number of meg channels used for rtSSS: " << NumCoil;
217 
218 // CoilTk.append("VV_PLANAR_T1");
219 // CoilTk.append("VV_MAG_T1");
220 
221  // Number of coil integration points for a coil
222  CoilNk.resize(NumCoil);
223 
224  // Coordinate of coil integrating points: This is assigned to CoilRk.
225  // 7002 (babyMEG circular (10mm diameter) inlayer magnetometer): 7 pts
226  // 7003 (babyMEG square (20x20mm) outlayer magnetometer): 7 pts
227  // 3012 (VectorView planar grdiometer, Type 1): 8 pts
228  // 3024 (VectorView magnetometer,Type 3): 9 pts
229  // initialization order X1.....Xn Y1.....Yn Z1.....Zn
230  MatrixXd C7002intpts7(3,7), C7003intpts4(3,4), C7003intpts16(3,16), C3012intpts8(3,8), C3024intpts9(3,9);
231  C7002intpts7 << 0.0000000, 0.0020412, -0.0020412, -0.0040825, -0.0020412, 0.0020412, 0.0040825, 0.0000000, 0.0035356, 0.0035356, 0.0000000, -0.0035356, -0.0035356, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000;
232  C7003intpts4 << 0.0050000, -0.0050000, -0.0050000, 0.0050000, 0.0050000, 0.0050000, 0.0050000, 0.0050000, 0.0000000, 0.0000000, 0.0000000, 0.0000000;
233  C7003intpts16 << 0.006, -0.006, -0.006, 0.006, 0.006, -0.006, -0.006, 0.006, 0.002, -0.002, -0.002, 0.002, 0.002, -0.002, -0.002, 0.002, 0.006, 0.006, -0.006, -0.006, 0.002, 0.002, -0.002, -0.002, 0.006, 0.006, -0.006, -0.006, 0.002, 0.002, -0.002, -0.002, 0.000, 0.000, 0.000, 0.000,0.000, 0.000, 0.000, 0.000,0.000, 0.000, 0.000, 0.000,0.000, 0.000, 0.000, 0.000;
234 
235  C3012intpts8 << 0.0059, -0.0059, 0.0059, -0.0059, 0.0108, -0.0108, 0.0108, -0.0108, 0.0067, 0.0067, -0.0067, -0.0067, 0.0067, 0.0067, -0.0067, -0.0067, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003;
236  C3024intpts9 << 0.0000, 0.0100, -0.0100, 0.0100, -0.0100, 0.0000, 0.0000, 0.0100, -0.0100, 0.0000, 0.0100, 0.0100, -0.0100, -0.0100, 0.0100, -0.0100, 0.0000, 0.0000, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003;
237 
238  // Coil point weight: This is assinged to CoilWk
239  MatrixXd C7002Wintpts7(1,7), C7003Wintpts4(1,4), C7003Wintpts16(1,16), C3012Wintpts8(1,8), C3024Wintpts9(1,9);
240  C7002Wintpts7 << 0.25, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125;
241  C7003Wintpts4 << 0.25, 0.25, 0.25, 0.25;
242  C7003Wintpts16 << 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625, 0.0625;
243 
244  C3012Wintpts8 << 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790;
245  C3024Wintpts9 << 0.1975, 0.0772, 0.0772, 0.0772, 0.0772, 0.1235, 0.1235, 0.1235, 0.1235;
246 
247  qint32 cid = 0;
248  CoilGrad.setZero(NumCoil);
249 
250  for (qint32 i=0; i<fiffInfo->nchan; ++i)
251  {
252  if(fiffInfo->chs[i].kind == FIFFV_MEG_CH && BadChan(i) == 0)
253  {
254 
255  CoilT.append(fiffInfo->chs[i].coil_trans);
256  CoilName.append(fiffInfo->chs[i].ch_name);
257 
258 // qDebug() << "coil type= " << fiffInfo->chs[i].coil_type;
259 
260  switch (fiffInfo->chs[i].coil_type)
261  {
262  case 3012:
263  CoilGrad(cid) = 1;
264  CoilNk(cid) = 8;
265  CoilRk.append(C3012intpts8);
266  CoilWk.append(C3012Wintpts8);
267  break;
268 
269  case 3024:
270  CoilGrad(cid) = 0;
271  CoilNk(cid) = 9;
272  CoilRk.append(C3024intpts9);
273  CoilWk.append(C3024Wintpts9);
274  break;
275 
276  case 7002:
277  CoilGrad(cid) = 0;
278  CoilNk(cid) = 7;
279  CoilRk.append(C7002intpts7);
280  CoilWk.append(C7002Wintpts7);
281  break;
282 
283  case 7003:
284  CoilGrad(cid) = 0;
285  CoilNk(cid) = 4;
286  CoilRk.append(C7003intpts4);
287  CoilWk.append(C7003Wintpts4);
288 // CoilRk.append(C7003intpts16);
289 // CoilWk.append(C7003Wintpts16);
290  break;
291 
292  default:
293  qDebug() << " This coil type is NOT supported";
294  break;
295  }
296 // qDebug() << fiffInfo->chs[i].coil_type << ", CoilNk(" << cid << "): " << CoilNk(cid);
297  cid++;
298  }
299  }
300 
301 // std::cout << "loading MEGData ....";
302 // ifstream inMEGData("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGData.txt",ios::in);
303 // MEGData.resize(NumCoil,1);
304 // VectorXd tmpvec(NumCoil);
305 // for(int k=0; k<NumCoil; k++) inMEGData >> tmpvec(k);
306 // MEGData.col(0) = tmpvec;
307 // inMEGData.close();
308 //NumCoil = 275;
309 }
310 
311 
312 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
313 //% build linear equation for SSS
314 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
315 //% Origin: origin of device coordinate system [X; Y; Z];
316 //% LIn: expansion order for internal basis functions
317 //% LOut: expansion order for external basis functions
318 //% CoilTk{i,1}: coil type - name
319 //% CoilNk(i,1): coil type - number of samples
320 //% CoilWk{i,1}: coil type - weight of samples
321 //% CoilRk{i,1}: coil type - location of samples
322 //% i: i-th type
323 //% CoilName{i,1}: name per coil
324 //% CoilT{i,1}: transformation matrix per coil
325 //% i: i-th coil
326 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
327 //% EqnIn(i,1): internal basis functions
328 //% EqnOut(i,1): external basis functions
329 //% i: i-th coil
330 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
331 // function [EqnIn,EqnOut] = get_SSS_Eqn(Origin,LIn,LOut,CoilTk,CoilNk,CoilWk,CoilRk,CoilName,CoilT)
332 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
333 QList<MatrixXd> RtSssAlgo::getSSSEqn(qint32 LIn, qint32 LOut)
334 //QList<MatrixXd> RtSssAlgo::getSSSEqn(VectorXi Lexp)
335 {
336  int NumBIn, NumBOut;
337 // int coil_index=1;
338  double RScale;
339  VectorXd coil_distance, coil_clocation(4,1), coil_vector(4,1);
340  MatrixXd coil_location;
341  MatrixXd EqnIn, EqnOut;
342  QList<MatrixXd> Eqn;
343 
344 // qint32 LIn,LOut;
345 // LIn = Lexp(0);
346 // LOut = Lexp(1);
347 
348  // % initialization
349  NumBIn = (LIn*LIn) + 2*LIn;
350  NumBOut = (LOut*LOut) + 2*LOut;
351 
352 //% calculate scaling factor for distance
353  coil_distance.resize(NumCoil);
354 
355  for(int i = 0; i<NumCoil; i++)
356  {
357  coil_clocation = CoilT[i].block(0,3,3,1);
358  coil_distance(i) = (coil_clocation-Origin).norm();
359  }
360 
361  RScale = exp(coil_distance.array().log().mean());
362 
363 //% build linear equation for internal/external basis functions
364  EqnIn.setZero(NumCoil,NumBIn);
365  EqnOut.setZero(NumCoil,NumBOut);
366 
367  for(int i = 0; i<NumCoil; i++)
368  {
369 // if (CoilName[i] == CoilTk[0])
370 // coil_index = 0;
371 // else if (CoilName[i] == CoilTk[1])
372 // coil_index = 1;
373 // else
374 // {
375 // cout<< "There is no matching coil" << endl;
376 // exit(1);
377 // }
378 
379 // % calculate coil orientation
380  coil_vector = CoilT[i].block(0,2,3,1);
381 
382 // % calculate coil locations (multiple points)
383 // int NumCoilPts = CoilNk(coil_index);
384  qint32 NumCoilPts = CoilNk(i);
385  MatrixXd tmpmat; tmpmat.setOnes(4,NumCoilPts);
386 //qDebug() << "pass 2";
387 // tmpmat.topRows(3) = CoilRk[coil_index];
388 // std::cout << "CoilRk: " << CoilRk[i].rows() << " x " << CoilRk[i].cols() << std::endl;
389  tmpmat.topRows(3) = CoilRk[i];
390 //qDebug() << "pass 3";
391  coil_location = CoilT[i] * tmpmat;
392  tmpmat = coil_location.topRows(3); coil_location.resize(3,NumCoilPts); coil_location = tmpmat;
393  coil_location = coil_location - Origin.replicate(1,NumCoilPts);
394  coil_location = coil_location / RScale;
395 
396 // % build linear equation
397  getSSSBasis(coil_location.row(0).transpose(), coil_location.row(1).transpose(), coil_location.row(2).transpose(), LIn, LOut);
398 // std::cout << "NumCoilPts: " << NumCoilPts << ", NumBIn: " << NumBIn << ", NumBOut: " << NumBOut << endl;
399 // std::cout << "BInX: " << BInX.rows() << " x " << BInX.cols() << endl;
400 
401  MatrixXd b_in, b_out;
402  b_in = coil_vector(0)*BInX + coil_vector(1)*BInY + coil_vector(2)*BInZ;
403  b_out = coil_vector(0)*BOutX + coil_vector(1)*BOutY + coil_vector(2)*BOutZ;
404 // std::cout << "b_in: Coil= " << i+1 << endl << b_in << endl;
405 // std::cout << "b_out: Coil= " << i+1 << endl << b_out << endl;
406 
407 // EqnIn.block(i,0,1,NumBIn) = CoilWk[coil_index] * b_in;
408 // EqnOut.block(i,0,1,NumBOut) = CoilWk[coil_index] * b_out;
409  EqnIn.block(i,0,1,NumBIn) = CoilWk[i] * b_in;
410  EqnOut.block(i,0,1,NumBOut) = CoilWk[i] * b_out;
411  }
412 
413  Eqn.append(EqnIn);
414  Eqn.append(EqnOut);
415 // std::cout << "EqnIn: Coil 1 " << endl << Eqn[0].row(0) << endl;
416 // std::cout << "EqnIn: Coil 306 " << endl << Eqn[0].row(305) << endl;
417 // std::cout << "EqnOut: Coil 1 " << endl << Eqn[1].row(0) << endl;
418 // std::cout << "EqnOut: Coil 306 " << endl << Eqn[1].row(305) << endl;
419 
420  return Eqn;
421 }
422 
423 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
424 //% generate SSS basis vectors at given locations
425 //% -- real basis function
426 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
427 //% X/Y/Z(i,1): Cartesian coordinate
428 //% i: i-th point
429 //% LIn: expansion order for internal basis functions
430 //% LOut: expansion order for external basis functions
431 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
432 //% BInX/Y/Z(i,j): internal basis vectors
433 //% BOutX/Y/Z(i,j): external basis vectors
434 //% i: i-th point
435 //% j: j-th basis function
436 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
437 // function [BInX,BInY,BInZ,BOutX,BOutY,BOutZ] = get_SSS_Basis(X,Y,Z,LIn,LOut)
438 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
439 void RtSssAlgo::getSSSBasis(VectorXd X, VectorXd Y, VectorXd Z, qint32 LIn, qint32 LOut)
440 {
441  int NumSample, NumBIn, NumBOut;
442  int LMax;
443 
444  QList<MatrixXd> P, dP_dx;
445  MatrixXd cur_p, cur_p1, cur_p2;
446 
447  QList<MatrixXcd> YP, dY_dTHETA;
448  VectorXcd cur_phi;
449 
450  MatrixXd BInR, BInPHI, BInTHETA;
451  int basis_index;
452  VectorXd scale_r;
453  VectorXcd cur_R, cur_PHI, cur_THETA;
454  MatrixXd BOutR, BOutPHI, BOutTHETA;
455 
456 // % check input parameters
457 // % -- R ~= 0 & THETA ~= 0
458  for(int i=0; i<X.rows(); i++)
459  if ( (abs(X(i)) < 1e-30) && (abs(Y(i)) < 1e-30) )
460  {
461  std::cout << "Zero THETA detected!**** ";
462  std::cout << "X, Y: " << X(i) << ", " << Y(i) << endl;
463  exit(1);
464  }
465 
466 // % initialization
467  LMax = qMax(LIn,LOut);
468  NumSample = X.size();
469  NumBIn = (LIn*LIn) + 2*LIn;
470  NumBOut = (LOut*LOut) + 2*LOut;
471 
472 // % Coordinate transdorm: (X,Y,Z) --> (R, PHI, THETA)
473  getCartesianToSpherCoordinate(X,Y,Z);
474 // std::cout << "X, Y, Z: " << endl << X.transpose() << endl << Y.transpose() << endl << Z.transpose() << endl;
475 // std::cout << "R, PHI, THETA: " << endl << R.transpose() << endl << PHI.transpose() << endl << THETA.transpose() << endl;
476 
477 // % calculate P
478  for(int l=1; l<=LMax; l++)
479  {
480  cur_p = legendre(l, THETA.array().cos());
481  P.append(cur_p.transpose());
482  }
483 // std::cout << "Legendre Polynomial 1-----" << endl << P[0].transpose() << endl;
484 // std::cout << "Legendre Polynomial 2-----" << endl << P[1].transpose() << endl;
485 // std::cout << "Legendre Polynomial 3-----" << endl << P[2].transpose() << endl;
486 // std::cout << "Legendre Polynomial 4-----" << endl << P[3].transpose() << endl;
487 // std::cout << "Legendre Polynomial 5-----" << endl << P[4].transpose() << endl;
488 
489 // % calculate dP_dx
490  for(int l=0; l<LMax; l++)
491  {
492  dP_dx.append(MatrixXd::Zero(NumSample,(l+1)+1));
493  for(int m=0; m <= l+1; m++)
494  {
495  cur_p1 = P[l].col(m);
496 // std::cout << "cur_p1. m= " << m << endl << cur_p1.transpose() << endl;
497  if (m == 0)
498  {
499  cur_p2 = -factorial((l+1)-1)/factorial((l+1)+1) * P[l].col(1);
500  }
501  else
502  {
503  cur_p2 = P[l].col(m-1);
504  }
505 // std::cout << "cur_p2. m=." << m << endl << cur_p2.transpose() << endl;
506 
507  // std::cout << "THETA: " << THETA.rows() << " x " << THETA.cols() << endl;
508  // std::cout << "cur_p1: " << cur_p1.rows() << " x " << cur_p1.cols() << endl;
509  // std::cout << "P[0]: " << P[0].col(1) << endl;
510 // std::cout << "dp_dx[" << l << "]: " << endl << dP_dx[l] << endl;
511  dP_dx[l].col(m) = (m*THETA.array().cos() * cur_p1.array() + (((l+1)+m)*((l+1)-m+1)*THETA.array().sin()) * cur_p2.array()) / THETA.array().sin().pow(2);
512  }
513 
514 // std::cout << "dP_dx[" << l << "]" <<endl << dP_dx[l] << endl;
515  }
516 
517 // % calculate Y & dY_dTHETA
518  for(int l=0; l<LMax; l++)
519  {
520  YP.append(MatrixXcd::Zero(NumSample,(l+1)+1));
521  dY_dTHETA.append(MatrixXcd::Zero(NumSample,(l+1)+1));
522 
523  for(int m=0; m <=l+1; m++)
524  {
525 // % cur_phi = sqrt((2*l+1)*factorial(l-m)/factorial(l+m)/2/pi) * exp(sqrt(-1)*m*PHI);
526 // % Y{l,1}(:,m+1) = cur_phi .* P{l,1}(:,m+1);
527 // % dY_dTHETA{l,1}(:,m+1) = -cur_phi .* sin(THETA) .* dP_dx{l,1}(:,m+1);
528  cur_phi = sqrt((2*(l+1)+1)*factorial((l+1)-m)/factorial((l+1)+m)/2/M_PI) * (sqrt(cplxd(-1,0)) * cplxd(m,0) * PHI.cast<cplxd >()).array().exp();
529 // std::cout << "cur_phi. m=." << m << endl << cur_phi.transpose() << endl;
530 
531  YP[l].col(m) = cur_phi.array() * P[l].col(m).array().cast<cplxd >();
532  dY_dTHETA[l].col(m) = -cur_phi.array() * THETA.array().sin().cast<cplxd >() * dP_dx[l].col(m).array().cast<cplxd >();
533  }
534 // std::cout << "YP[" << l << "]" <<endl << YP[l] << endl;
535 // std::cout << "dY_dTHETA[" << l << "]" <<endl << dY_dTHETA[l] << endl;
536  }
537 
538 // % calculate BIn (spherical coordinate)
539  BInR.setZero(NumSample,NumBIn);
540  BInPHI.setZero(NumSample,NumBIn);
541  BInTHETA.setZero(NumSample,NumBIn);
542  basis_index = 0;
543  for(int l=0; l<LIn; l++)
544  {
545  scale_r= R.array().pow((l+1)+2);
546 
547  BInR.col(basis_index) = -((l+1)+1) * YP[l].col(0).array().real() / scale_r.array();
548  BInPHI.col(basis_index) = VectorXd::Zero(NumSample);
549  BInTHETA.col(basis_index) = dY_dTHETA[l].col(0).array().real() / scale_r.array();
550  basis_index = basis_index + 1;
551 
552 // std::cout << "BInR" << endl << BInR << endl;
553 
554  for(int m=0; m<=l; m++)
555  {
556  cur_R = -sqrt(2) * ((l+1)+1) * YP[l].col(m+1).array() / scale_r.array().cast<cplxd >();
557  cur_PHI = sqrt(cplxd(-2,0)) * cplxd(m+1,0) * YP[l].col(m+1).array() / THETA.array().sin().cast<cplxd >() / scale_r.array().cast<cplxd >();
558  cur_THETA = sqrt(2) * dY_dTHETA[l].col(m+1).array() / scale_r.array().cast<cplxd >();
559 // std::cout << "cur_R, m=" << m << endl << cur_R.transpose() << endl;
560 
561  BInR.col(basis_index) = cur_R.real();
562  BInPHI.col(basis_index) = cur_PHI.real();
563  BInTHETA.col(basis_index) = cur_THETA.real();
564  basis_index = basis_index + 1;
565 
566  BInR.col(basis_index) = cur_R.imag();
567  BInPHI.col(basis_index) = cur_PHI.imag();
568  BInTHETA.col(basis_index) = cur_THETA.imag();
569  basis_index = basis_index + 1;
570  }
571  }
572 // std::cout << "BInR" << endl << BInR << endl;
573 // std::cout << "BInPHI" << endl << BInPHI << endl;
574 // std::cout << "BInTHETA" << endl << BInTHETA << endl;
575 
576 // % calculate BOut (spherical coordinate)
577  BOutR.setZero(NumSample,NumBOut);
578  BOutPHI.setZero(NumSample,NumBOut);
579  BOutTHETA.setZero(NumSample,NumBOut);
580  basis_index = 0;
581  for(int l=0; l<LOut; l++)
582  {
583  scale_r= R.array().pow((l+1)-1);
584 
585  BOutR.col(basis_index) = (l+1) * YP[l].col(0).array().real() * scale_r.array();
586  BOutPHI.col(basis_index) = VectorXd::Zero(NumSample);
587  BOutTHETA.col(basis_index) = dY_dTHETA[l].col(0).array().real() * scale_r.array();
588  basis_index = basis_index + 1;
589 
590  for(int m=0; m<=l; m++)
591  {
592  cur_R = sqrt(2) * (l+1) * YP[l].col(m+1).array() * scale_r.array().cast<cplxd >();
593  cur_PHI = sqrt(cplxd(-2,0)) * cplxd(m+1,0) * YP[l].col(m+1).array() / THETA.array().sin().cast<cplxd >() * scale_r.array().cast<cplxd >();
594  cur_THETA = sqrt(2) * dY_dTHETA[l].col(m+1).array() * scale_r.array().cast<cplxd >();
595  BOutR.col(basis_index) = cur_R.real();
596  BOutPHI.col(basis_index) = cur_PHI.real();
597  BOutTHETA.col(basis_index) = cur_THETA.real();
598  basis_index = basis_index + 1;
599  BOutR.col(basis_index) = cur_R.imag();
600  BOutPHI.col(basis_index) = cur_PHI.imag();
601  BOutTHETA.col(basis_index) = cur_THETA.imag();
602  basis_index = basis_index + 1;
603  }
604  }
605 // std::cout << "BOutR" << endl << BOutR << endl;
606 // std::cout << "BOutPHI" << endl << BOutPHI << endl;
607 // std::cout << "BOutTHETA" << endl << BOutTHETA << endl;
608 
609 // % convert from spherical coordinate to Cartesian coordinate
610 // [R_X,R_Y,R_Z,PHI_X,PHI_Y,PHI_Z,THETA_X,THETA_Y,THETA_Z] = get_Sph_2_Cart_Vector(R,PHI,THETA);
611  getSphereToCartesianVector();
612  BInX = BInR.array() * (R_X.replicate(1,NumBIn)).array() + BInPHI.array() * (PHI_X.replicate(1,NumBIn)).array() + BInTHETA.array() * (THETA_X.replicate(1,NumBIn)).array();
613  BInY = BInR.array() * (R_Y.replicate(1,NumBIn)).array() + BInPHI.array() * (PHI_Y.replicate(1,NumBIn)).array() + BInTHETA.array() * (THETA_Y.replicate(1,NumBIn)).array();
614  BInZ = BInR.array() * (R_Z.replicate(1,NumBIn)).array() + BInPHI.array() * (PHI_Z.replicate(1,NumBIn)).array() + BInTHETA.array() * (THETA_Z.replicate(1,NumBIn)).array();
615  BOutX = BOutR.array() * (R_X.replicate(1,NumBOut)).array() + BOutPHI.array() * (PHI_X.replicate(1,NumBOut)).array() + BOutTHETA.array() * (THETA_X.replicate(1,NumBOut)).array();
616  BOutY = BOutR.array() * (R_Y.replicate(1,NumBOut)).array() + BOutPHI.array() * (PHI_Y.replicate(1,NumBOut)).array() + BOutTHETA.array() * (THETA_Y.replicate(1,NumBOut)).array();
617  BOutZ = BOutR.array() * (R_Z.replicate(1,NumBOut)).array() + BOutPHI.array() * (PHI_Z.replicate(1,NumBOut)).array() + BOutTHETA.array() * (THETA_Z.replicate(1,NumBOut)).array();
618 
619 // std::cout << "BInX" << endl << BInX << endl;
620 // std::cout << "BInY" << endl << BInY << endl;
621 // std::cout << "BInZ" << endl << BInZ << endl;
622 // std::cout << "BOutX" << endl << BOutX << endl;
623 // std::cout << "BOutY" << endl << BOutY << endl;
624 // std::cout << "BOutZ" << endl << BOutZ << endl;
625 }
626 
627 
628 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
629 //%
630 //% RtSssAlgo::getCartesianToSpherCeoordinate(double X, double Y, double Z
631 //%
632 //% convert point coordinate from Cartesian to spherical
633 //%
634 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
635 //% X/Y/Z(i,1): Cartesian coordinate
636 //% i: i-th point
637 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
638 //% R/PHI/THETA(i,1): spherical coordinate
639 //% i: i-th point
640 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
641 // function [R,PHI,THETA] = get_Cart_2_Sph_Coordinate(X,Y,Z)
642 // hypotxy = hypot(X,Y);
643 // R = hypot(hypotxy,Z);
644 // PHI = atan2(Y,X);
645 // THETA = atan2(hypotxy,Z);
646 // return;
647 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
648 void RtSssAlgo::getCartesianToSpherCoordinate(VectorXd X, VectorXd Y, VectorXd Z)
649 {
650  VectorXd hypotxy;
651 // std::cout << "X: " << X.rows() << " x " << X.cols() << endl;
652 // std::cout << "Y: " << Y.rows() << " x " << Y.cols() << endl;
653 
654  hypotxy = hypot(X,Y);
655  R = hypot(hypotxy,Z);
656  PHI = atan2vec(Y,X);
657  THETA = atan2vec(hypotxy,Z);
658 }
659 
660 
661 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
662 //%
663 //% RtSssAlgo::getSphereToCartesianVector(double R, double PHI, double THETA)
664 //%
665 //% convert unit vector from spherical to Cartesian
666 //%
667 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
668 //% R/PHI/THETA(i,1): spherical coordinate
669 //% i: i-th point
670 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
671 //% R_X/R_Y/R_Z(i,1): unit vector R in Cartesian coordinate
672 //% PHI_X/PHI_Y/PHI_Z(i,1): unit vector PHI in Cartesian coordinate
673 //% THETA_X/THETA_Y/THETA_Z(i,1): unit vector THETA in Cartesian coordinate
674 //% i: i-th point
675 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
676 // function [R_X,R_Y,R_Z,PHI_X,PHI_Y,PHI_Z,THETA_X,THETA_Y,THETA_Z] = get_Sph_2_Cart_Vector(R,PHI,THETA)
677 // NumExp = size(R,1);
678 // R_X = sin(THETA) .* cos(PHI);
679 // R_Y = sin(THETA) .* sin(PHI);
680 // R_Z = cos(THETA);
681 // PHI_X = -sin(PHI);
682 // PHI_Y = cos(PHI);
683 // PHI_Z = zeros(NumExp,1);
684 // THETA_X = cos(THETA) .* cos(PHI);
685 // THETA_Y = cos(THETA) .* sin(PHI);
686 // THETA_Z = -sin(THETA);
687 // return;
688 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
689 void RtSssAlgo::getSphereToCartesianVector()
690 {
691  int NumExp = R.size();
692 
693  R_X = THETA.array().sin() * PHI.array().cos();
694  R_Y = THETA.array().sin() * PHI.array().sin();
695  R_Z = THETA.array().cos();
696 
697  PHI_X = -PHI.array().sin();
698  PHI_Y = PHI.array().cos();
699  PHI_Z = VectorXd::Zero(NumExp,1);
700 
701 
702  THETA_X = THETA.array().cos() * PHI.array().cos();
703  THETA_Y = THETA.array().cos() * PHI.array().sin();
704  THETA_Z = -THETA.array().sin();
705 }
706 
707 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
708 //% SSS by robust regression (subspace iteration & low-rank update)
709 //% -- RR_K1 is determined by bi-square function (w = 0.75 at RR_K1)
710 //% -- RR_K2 is determined by bi-square function (95% statistical efficiency)
711 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
712 //% EqnIn(i,j): internal basis functions
713 //% EqnOut(i,j): external basis functions
714 //% i: i-th coil
715 //% j: j-th basis function
716 //% EqnARR(i,j): SSS equation (LHS) after scaling -- subspace
717 //% i: i-th coil
718 //% j: j-th basis function
719 //% EqnA(i,j): SSS equation (LHS) after scaling -- full
720 //% i: i-th coil
721 //% j: j-th basis function
722 //% EqnB(i,j): SSS equation (RHS) after scaling
723 //% i: i-th coil
724 //% j: j-th sample in time domain
725 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
726 //% SSSIn(i,j): internal MEG signal recovered by SSS
727 //% SSSOut(i,j): external MEG signal recovered by SSS
728 //% Weight(i,j): optimal weight for MEG coil
729 //% ErrRel(1,j): relative residual of over-determined linear equation
730 //% i: i-th coil
731 //% j: j-th sample in time domain
732 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
733 // function [SSSIn,SSSOut,Weight,ErrRel] = get_SSS_RR(EqnIn,EqnOut,EqnARR,EqnA,EqnB)
734 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
735 
736 //QList<MatrixXd> RtSssAlgo::getSSSRR(MatrixXd EqnIn, MatrixXd EqnOut, MatrixXd EqnARR, MatrixXd EqnA, MatrixXd EqnB)
737 //QList<MatrixXd> RtSssAlgo::getSSSRR(MatrixXd EqnB)
738 MatrixXd RtSssAlgo::getSSSRR(MatrixXd EqnB)
739 {
740  int NumBIn, NumBOut, NumCoil, NumExp;
741  MatrixXd EqnRRInv, EqnInv, SSSIn, SSSOut, Weight; //, ErrRel;
742  VectorXd ErrRel;
743  double RR_K1, RR_K2, RR_K3;
744  double eqn_scale0, eqn_scale;
745  MatrixXd sol_X, sol_X_old, eqn_Y, eqn_D, temp_M, temp_N, sol_in, sol_out;
746  MatrixXd diagMat;
747  VectorXd eqn_err, weight_index;
748  QList<MatrixXd> RRsss;
749 
750 // % error tolerance for robust regression
751  double ErrTolRel = 1e-3;
752 
753 // % weight threshold for robust regression
754  double WeightThres = 1 - 1e-6;
755 
756 // % initialization
757  NumBIn = EqnIn.cols();
758  NumBOut = EqnOut.cols();
759  NumCoil = EqnB.rows();
760  NumExp = EqnB.cols();
761 
762  EqnRRInv = (EqnARR.transpose() * EqnARR).inverse();
763  EqnInv = (EqnA.transpose() * EqnA).inverse();
764 
765  SSSIn.setZero(NumCoil,NumExp);
766  SSSOut.setZero(NumCoil,NumExp);
767  Weight.setZero(NumCoil,NumExp);
768  ErrRel.setZero(NumExp);
769  RR_K3 = 3;
770  RR_K2 = 4.685;
771  RR_K1 = qSqrt(1-qSqrt(3)/2) * RR_K2;
772 
773  for(int i=0; i<NumExp; i++)
774  {
775 // % solve OLS solution
776  sol_X = EqnRRInv * (EqnARR.transpose() * EqnB.col(i));
777 // std::cout << "sol_X ************************" << endl << sol_X.transpose() << endl;
778 
779 // % scale linear equation
780  eqn_err = EqnARR * sol_X - EqnB.col(i);
781  eqn_scale0 = stdev(eqn_err);
782  eqn_err = eqn_err.cwiseAbs() / eqn_scale0;
783 
784 // % solve iteratively re-weighted least squares (Bi-Square) -- subspace
785  sol_X_old.setConstant(sol_X.rows(), sol_X.cols(), 1e30);
786 // while (((sol_X-sol_X_old).norm() / sol_X.norm()) > ErrTolRel)
787  int cnt=0;
788  while (((sol_X.array()-sol_X_old.array()).matrix().norm() / sol_X.norm()) > ErrTolRel)
789 // for (int h=0; h<2; h++)
790  {
791  cnt++;
792 // cout << "while cnt: " << cnt << endl;
793  sol_X_old = sol_X;
794 // Weight(:,i) = (eqn_err <= RR_K1) + (eqn_err > RR_K1 & eqn_err <= RR_K2) .* (1-(eqn_err-RR_K1).^2/(RR_K2-RR_K1)^2).^2;
795 // Weight.col(i) = eigen_LTE(eqn_err, RR_K1) + eigen_AND(eigen_GT(eqn_err, RR_K1),eigen_LTE(eqn_err, RR_K2));
796  Weight.col(i) = eigen_LTE(eqn_err,RR_K1).array() + eigen_AND(eigen_GT(eqn_err,RR_K1),eigen_LTE(eqn_err,RR_K2)).array() * (1 - ((eqn_err.array()-RR_K1).pow(2)) / pow(RR_K2-RR_K1,2) ).pow(2);
797 // std::cout << "Weight **************************" << endl << Weight.transpose() << endl;
798 
799 // % weight_index = find(Weight(:,i) < WeightThres);
800  weight_index = eigen_LT_index(Weight.col(i), WeightThres);
801 // weight_index = eigen_LT_index_test(Weight.col(i), h);
802 // std:cout << "weight_index ********************* " << endl << weight_index << endl;
803 
804 // % eqn_Y = EqnARR(weight_index,:); eqn_D = Weight(weight_index,i) - 1;
805  eqn_Y.resize(weight_index.size(), EqnARR.cols());
806  eqn_D.resize(weight_index.size(),1);
807  for(int k=0; k<weight_index.size(); k++)
808  {
809  eqn_Y.row(k) = EqnARR.row(weight_index(k));
810 // std::cout << "eqn_Y: " << eqn_Y.rows() << " x " << eqn_Y.cols() << endl;
811 // std::cout << "eqn_D: " << eqn_D.rows() << " x " << eqn_D.cols() << endl;
812  eqn_D(k) = Weight(weight_index(k),i) - 1;
813  }
814 // std::cout << "eqn_D: " << endl << eqn_D.transpose() << endl << endl;
815  temp_M = EqnARR.transpose() * (Weight.col(i).array() * EqnB.col(i).array()).matrix();
816 // std::cout << "temp_M: " << endl << temp_M.transpose() << endl << endl;
817 
818 // std::cout << "EqnRRInv: " << EqnRRInv.rows() << " x " << EqnRRInv.cols() << endl;
819 // std::cout << "eqn_Y: " << eqn_Y.rows() << " x " << eqn_Y.cols() << endl;
820 // std::cout << "EqnRRInv: " << endl << EqnRRInv.transpose() << endl << endl;
821 // std::cout << "eqn_Y: " << endl << eqn_Y.transpose() << endl << endl;
822 
823  temp_N = EqnRRInv * eqn_Y.transpose();
824 // std::cout << "temp_N: " << endl << temp_N.transpose() << endl << endl;
825 
826 
827 // diagMat = eqn_D.asDiagonal();
828 // std::cout << "diagMat ************************" << endl << diagMat.transpose() << endl;
829 // diagMat = 1 / diagMat.array();
830 // std::cout << "diagMat ************************" << endl << diagMat.transpose() << endl;
831 
832  diagMat = (1 / eqn_D.array()).matrix().asDiagonal();
833 // cout << "diagMat ************************" << i+1 << endl << diagMat.rows() << " x" << diagMat.cols() << endl;
834 
835  sol_X = EqnRRInv * temp_M - temp_N * (diagMat + eqn_Y * temp_N).inverse() * (temp_N.transpose() * temp_M);
836  eqn_err = (EqnARR * sol_X - EqnB.col(i)).cwiseAbs();
837  eqn_scale = qMin(eqn_scale0, RR_K3 * qSqrt((Weight.col(i).array() * eqn_err.array() * eqn_err.array()).mean()));
838  eqn_err = eqn_err / eqn_scale;
839 // std::cout << "temp_M: " << endl << temp_M.transpose() << endl << endl;
840 
841 // std::cout << "(sol_X.array()-sol_X_old.array()).matrix().norm(): " << (sol_X.array()-sol_X_old.array()).matrix().norm() << endl << endl;
842 // std::cout << "sol_X.norm(): " << sol_X.norm() << endl << endl;
843 // std::cout << "sol_X ************************" << endl << sol_X.transpose() << endl << endl;
844 
845  }
846 // % solve weighted SSS - full
847 // % eqn_Y = EqnA(weight_index,:); temp_M = EqnA' * (Weight(:,i).*EqnB(:,i)); temp_N = EqnInv * eqn_Y';
848  eqn_Y.resize(weight_index.size(), NumBIn+NumBOut);
849  for(int k=0; k<weight_index.size(); k++) eqn_Y.row(k) = EqnA.row(weight_index(k));
850  temp_M = EqnA.transpose() * (Weight.col(i).array() * EqnB.col(i).array()).matrix();
851  temp_N = EqnInv * eqn_Y.transpose();
852 
853 
854 // % sol_X = EqnInv * temp_M - temp_N * ((diag(1./eqn_D) + eqn_Y * temp_N); // \ (temp_N'*temp_M));
855 // diagMat = eqn_D.asDiagonal();
856 // diagMat = 1 / diagMat.array();
857  diagMat = (1 / eqn_D.array()).matrix().asDiagonal();
858  sol_X = EqnInv * temp_M - temp_N * (diagMat + eqn_Y * temp_N).inverse() * (temp_N.transpose() * temp_M);
859 
860  ErrRel(i) = (EqnA * sol_X - EqnB.col(i)).norm() / EqnB.col(i).norm();
861 
862  sol_in = sol_X.block(0,0,NumBIn,1);
863  sol_out = sol_X.block(NumBIn,0,NumBOut,1); // sol_X(NumBIn+1:NumBIn+NumBOut,1);
864 
865 // % recover internal/external MEG siganl
866  SSSIn.col(i) = EqnIn * sol_in;
867  SSSOut.col(i) = EqnOut * sol_out;
868 
869  }
870  RRsss.append(SSSIn);
871  RRsss.append(SSSOut);
872  RRsss.append(Weight);
873  RRsss.append(ErrRel);
874 
875 // return RRsss;
876  return SSSIn;
877 }
878 
879 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
880 //% SSS by OLS
881 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
882 //% EqnIn(i,j): internal basis functions
883 //% EqnOut(i,j): external basis functions
884 //% i: i-th coil
885 //% j: j-th basis function
886 //% EqnA(i,j): SSS equation (LHS) after scaling
887 //% i: i-th coil
888 //% j: j-th basis function
889 //% EqnB(i,j): SSS equation (RHS) after scaling
890 //% i: i-th coil
891 //% j: j-th sample in time domain
892 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
893 //% SSSIn(i,j): internal MEG signal recovered by SSS
894 //% SSSOut(i,j): external MEG signal recovered by SSS
895 //% ErrRel(1,j): relative residual of over-determined linear equation
896 //% i: i-th coil
897 //% j: j-th sample in time domain
898 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
899 //function [SSSIn,SSSOut,ErrRel] = get_SSS_OLS(EqnIn,EqnOut,EqnA,EqnB)
900 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
901 //QList<MatrixXd> RtSssAlgo::getSSSOLS(MatrixXd EqnIn, MatrixXd EqnOut, MatrixXd EqnA, MatrixXd EqnB)
902 //QList<MatrixXd> RtSssAlgo::getSSSOLS(MatrixXd EqnB)
903 MatrixXd RtSssAlgo::getSSSOLS(MatrixXd EqnB)
904 {
905  int NumBIn, NumBOut, NumCoil, NumExp;
906  MatrixXd EqnRRInv, EqnInv, SSSIn, SSSOut;
907  MatrixXd sol_X, sol_in, sol_out;
908  VectorXd ErrRel;
909  QList<MatrixXd> OLSsss;
910 
911 // % initialization
912  NumBIn = EqnIn.cols();
913  NumBOut = EqnOut.cols();
914  NumCoil = EqnB.rows();
915  NumExp = EqnB.cols();
916 
917  EqnRRInv = (EqnA.transpose() * EqnA).inverse();
918  SSSIn.setZero(NumCoil,NumExp);
919  SSSOut.setZero(NumCoil,NumExp);
920  ErrRel.setZero(NumExp);
921 
922  for (int i=0; i<NumExp; i++)
923  {
924 // % solve OLS solution
925  sol_X = EqnRRInv * (EqnA.transpose() * EqnB.col(i));
926 
927  ErrRel(i) = (EqnA * sol_X - EqnB.col(i)).norm() / EqnB.col(i).norm();
928 
929  sol_in = sol_X.block(0,0,NumBIn,1);
930 
931  sol_out = sol_X.block(NumBIn,0, NumBOut,1); // probably NumBIn+0 is correct. Please debug !!!!
932 
933 // % recover internal/external MEG siganl
934  SSSIn.col(i) = EqnIn * sol_in;
935  SSSOut.col(i) = EqnOut * sol_out;
936  }
937 
938  OLSsss.append(SSSIn);
939  OLSsss.append(SSSOut);
940  OLSsss.append(ErrRel);
941 
942 // return OLSsss;
943  return SSSIn;
944 }
945 
946 // Return number of meg channels
947 qint32 RtSssAlgo::getNumMEGChanUsed()
948 {
949  return NumCoil;
950 }
951 
952 qint32 RtSssAlgo::getNumMEGChan()
953 {
954  return NumMEGChan;
955 }
956 
957 VectorXi RtSssAlgo::getBadChan()
958 {
959  return BadChan;
960 }
961 
962 qint32 RtSssAlgo::getNumMEGBadChan()
963 {
964  return NumBadCoil;
965 }
966 
967 // Set MEG signal
968 //void RtSssAlgo::setMEGsignal(MatrixXd megfrombuffer)
969 //{
970 // MEGData = megfrombuffer;
971 // EqnB = megfrombuffer;
973 //}
974 
975 QList<MatrixXd> RtSssAlgo::getLinEqn()
976 {
977  QList<MatrixXd> LinEqn;
978 
979  LinEqn.append(EqnIn);
980  LinEqn.append(EqnOut);
981  LinEqn.append(EqnARR);
982  LinEqn.append(EqnA);
983  LinEqn.append(EqnB);
984 
985  return LinEqn;
986 }
987 
988 
989 // Legendre Polyn0mial
990 MatrixXd legendre(int l, VectorXd x)
991 {
992  MatrixXd outP;
993 // std::cout << "x: " << x.rows() << " x " << x.cols() << endl;
994 
995  outP.setZero(l+1, x.size());
996 // std::cout << " outP: " << outP.rows() << " x " << outP.cols() << endl;
997 
998  for(int m=0; m<=l; m++)
999  {
1000  for(int k=0; k<x.size(); k++)
1001  outP(m,k) = plgndr(l,m,x(k));
1002  }
1003 
1004  return outP;
1005 }
1006 
1007 
1008 float plgndr(int l, int m, float x)
1009 {
1010 // void nrerror(char error_text[]);
1011  float fact,pll,pmm,pmmp1,somx2;
1012  int i,ll;
1013 
1014  if (m < 0 || m > l || fabs(x) > 1.0)
1015  std::cout << "Bad arguments in routine plgndr";
1016 // nrerror("Bad arguments in routine plgndr");
1017  pmm=1.0;
1018  if (m > 0)
1019  {
1020  somx2=sqrt((1.0-x)*(1.0+x));
1021  fact=1.0;
1022  for (i=1;i<=m;i++)
1023  {
1024  pmm *= -fact*somx2;
1025  fact += 2.0;
1026  }
1027  }
1028  if (l == m)
1029  return pmm;
1030  else {
1031  pmmp1=x*(2*m+1)*pmm;
1032  if (l == (m+1))
1033  return pmmp1;
1034  else
1035  {
1036  for (ll=m+2;ll<=l;ll++)
1037  {
1038  pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
1039  pmm=pmmp1;
1040  pmmp1=pll;
1041  }
1042  return pll;
1043  }
1044  }
1045 }
1046 
1047 //---------------------------------------------------------------------
1048 /* MATLAB: hypot.m
1049  Robust computation of the square root of the sum of squares.
1050  C = hypot(A,B) returns SQRT(ABS(A).^2+ABS(B).^2) carefully computed to
1051  avoid underflow and overflow.*/
1052 VectorXd hypot(VectorXd X, VectorXd Y)
1053 {
1054  VectorXd rst;
1055  rst = (X.array().abs().pow(2) + Y.array().abs().pow(2)).sqrt();
1056  return rst;
1057 }
1058 
1059 //---------------------------------------------------------------------
1060 // atan2() for vector
1061 // since no
1062 VectorXd atan2vec(VectorXd a, VectorXd b)
1063 {
1064  int n = a.size();
1065  VectorXd at2(n);
1066 
1067  for(int i=0; i<n; i++)
1068  {
1069  at2(i) = qAtan2(a(i),b(i));
1070  }
1071 
1072  return at2;
1073 }
1074 
1075 VectorXd find(MatrixXd M, int I)
1076 {
1077  VectorXd outIDX;
1078  int k=0;
1079  for(int i=0; i<M.cols(); i++)
1080  for(int j=0; j<M.rows(); j++)
1081  if (M(i,j) == I)
1082  {
1083  outIDX(k) = i*M.rows() + j;
1084  k++;
1085  }
1086  return outIDX;
1087 }
1088 
1089 
1090 
1091 double factorial(int n)
1092 {
1093  double fac=1;
1094 
1095  for(int i=2; i<=n; i++) fac *= i;
1096 
1097  return fac;
1098 }
1099 
1100 
1101 //----------------------------------------------------------------------
1102 // standard deviation of vector
1103 double stdev(VectorXd V)
1104 {
1105  int size = V.size();
1106  double ave, sum=0;
1107 
1108  ave = V.mean();
1109  for(int i=0; i<size; i++)
1110  {
1111  sum += pow((V(i)-ave),2);
1112  }
1113  return sqrt(sum/size);
1114 }
1115 
1116 //
1117 // Check matrix elements to see if they are less than or equal to the given tolerance
1118 // Returns a vector whose elements consists of 1 if true, or 0 if false
1119 VectorXd eigen_LTE(VectorXd V, double tol)
1120 {
1121  int row = V.rows();
1122  VectorXd outTrue = VectorXd::Zero(row);
1123 
1124  for(int i=0; i<row; i++)
1125  if(V(i) <= tol) outTrue(i) = 1;
1126 // std::cout << "LTE **********************" << outTrue.transpose() << endl;
1127  return outTrue;
1128 }
1129 
1130 //
1131 // Check matrix elements to see if they are greater than the given tolerance
1132 // Returns a vector whose elements consists of 1 if true, or 0 if false
1133 VectorXd eigen_GT(VectorXd V, double tol)
1134 {
1135  int row = V.rows();
1136  VectorXd outTrue = VectorXd::Zero(row);
1137 
1138  for(int i=0; i<row; i++)
1139  if(V(i) > tol) outTrue(i) = 1;
1140  return outTrue;
1141 }
1142 
1143 //
1144 // Check if the elements of two matrices are both true
1145 // Returns a vector whose elements consists of 1 if true, or 0 if false
1146 VectorXd eigen_AND(VectorXd V1, VectorXd V2)
1147 {
1148  int row = V1.rows();
1149  VectorXd outTrue = VectorXd::Zero(row);
1150 
1151  for(int i=0; i<row; i++)
1152  if(V1(i) && V2(i)) outTrue(i) = 1;
1153 // std::cout << "AND **********************" << outTrue.transpose() << endl;
1154  return outTrue;
1155 }
1156 
1157 // Check matrix elements to see if they are less than the given tolerance
1158 // Returns a vector whose elements consists index of true
1159 VectorXd eigen_LT_index(VectorXd V, double tol)
1160 {
1161  int cnt = 0;
1162  int size = V.size();
1163  VectorXd outIdx; outIdx.setZero(size);
1164 
1165  for(int i=0; i<size; i++)
1166  if(V(i) < tol)
1167  {
1168  outIdx(cnt) = i;
1169  cnt++;
1170  }
1171  return outIdx.head(cnt);
1172 }
1173 
1174 //VectorXd eigen_LT_index_test(VectorXd V, int cnt)
1175 VectorXd eigen_LT_index_test(int cnt)
1176 {
1177  VectorXd outIdx = VectorXd::Zero(cnt+1);
1178 
1179  for (int i=0; i<cnt+1; i++)
1180  outIdx(i) = i+1;
1181 
1182 // std::cout << "outIdx: " << outIdx.rows() << " x " << outIdx.cols() << ", cnt: " << cnt << endl;
1183  return outIdx;
1184 }
1185 
1186 //******************************************************************************************************************
1187 //******************************************************************************************************************
1188 //******************************************************************************************************************
1189 
1190 // Read FIFF coil configuration and assign them to the Coil class variables
1191 // such as CoilName,CoilT,CoilTk, CoilNk, CoilRk, CoilWk, CoilGrad.
1192 void RtSssAlgo::getCoilInfoBabyMEG4Sim()
1193 {
1194  int NumCoilIn = 270;
1195  int NumCoilOut = 105;
1196 
1197  NumCoil = NumCoilIn + NumCoilOut;
1198  cout << "Num of Coil = " << NumCoil << endl;
1199 
1200 // simulation data from test1_sim.mat
1201  Origin.resize(3);
1202  Origin << 0.0, 0.0, 0.04;
1203 
1204  std::cout << "loading MEGIn ....";
1205  ifstream inMEGIn("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGIn2.txt",ios::in);
1206  MEGIn.resize(NumCoil);
1207  for (int k=0; k<NumCoil; k++) inMEGIn >> MEGIn(k);
1208  inMEGIn.close();
1209 // std::cout << MEGIn << endl;
1210  std::cout << " finished !" << endl;
1211 
1212  std::cout << "loading MEGOut ....";
1213  ifstream inMEGOut("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGOut2.txt",ios::in);
1214  MEGOut.resize(NumCoil);
1215  for(int k=0; k<NumCoil; k++) inMEGOut >> MEGOut(k);
1216  inMEGOut.close();
1217 // std::cout << MEGOut << endl;
1218  std::cout << " finished !" << endl;
1219 
1220  std::cout << "loading MEGNoise ....";
1221  ifstream inMEGNoise("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGNoise.txt",ios::in);
1222  MEGNoise.resize(NumCoil);
1223  for(int k=0; k<NumCoil; k++) inMEGNoise >> MEGNoise(k);
1224  inMEGNoise.close();
1225 // std::cout << MEGNoise << endl;
1226  std::cout << " finished !" << endl;
1227 
1228  std::cout << "loading MEGData ....";
1229  ifstream inMEGData("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGData2.txt",ios::in);
1230  MEGData.resize(NumCoil,1);
1231  VectorXd tmpvec(NumCoil);
1232  for(int k=0; k<NumCoil; k++) inMEGData >> tmpvec(k);
1233  MEGData.col(0) = tmpvec;
1234  inMEGData.close();
1235 // std::cout << MEGData << endl;
1236  std::cout << " finished !" << endl;
1237 
1238  // name per coil
1239  for (int i=0; i<NumCoilIn; i++)
1240  CoilName.append("BSQ2_In_Layer");
1241  for (int i=0; i<NumCoilOut; i++)
1242  CoilName.append("BSQ2_Out_Layer");
1243 
1244 // std:cout << CoilName[0] << ", " << CoilName[1] << ", " << CoilName[2] << endl;
1245 
1246  // transformation matrix per coil
1247  std::cout << "loading coil transformation matrix ....";
1248  ifstream inCT("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_coil_def/BabyMEG/CoilT.txt",ios::in);
1249  MatrixXd tmat(4,4);
1250  for(int k=0; k<NumCoil; k++)
1251  {
1252  for(int i=0; i<4; i++)
1253  for(int j=0; j<4; j++)
1254  {
1255  inCT >> tmat(i,j);
1256  }
1257  CoilT.append(tmat);
1258  }
1259  inCT.close();
1260 // std::cout << CoilT[305] << endl;
1261  std::cout << " finished !" << endl;
1262 
1263 
1264  // coil type - name
1265  CoilTk.append("BSQ2_In_Layer");
1266  CoilTk.append("BSQ2_Out_Layer");
1267 // std::cout << CoilTk[0] << ", " << CoilTk[1] << endl;
1268 
1269 
1270  // coil type - Number of integration points for each coil type
1271  CoilNk.resize(2);
1272  CoilNk << 7, 7;
1273 
1274  // coil type - Coordinate of integrating points
1275  MatrixXd CintptsInlayer(3,7), CintptsOutlayer(3,7);
1276  CintptsInlayer << 0.000000, 0.000000, 0.000000, 0.0040825, 0.000000, 0.000000, -0.0040825, 0.000000, 0.000000, 0.0020412, 0.0035356, 0.000000, 0.0020412, -0.0035356, 0.000000, -0.0020412, 0.0035356, 0.000000, -0.0020412, -0.0035356, 0.000000;
1277  CintptsOutlayer << 0.000000, 0.000000, 0.000000, 0.0081650, 0.000000, 0.000000, -0.0081650, 0.000000, 0.000000, 0.0040824, 0.0070712, 0.000000, 0.0040824, -0.0070712, 0.000000, -0.0040824, 0.0070712, 0.000000, -0.0040824, -0.0070712, 0.000000;
1278  CoilRk.append(CintptsInlayer);
1279  CoilRk.append(CintptsOutlayer);
1280 
1281  // Coil type - Coil point weight
1282  MatrixXd CWintptsInlayer(1,7), CWintptsOutlayer(1,7);
1283  CWintptsInlayer << 0.25, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125;
1284  CWintptsOutlayer = CWintptsInlayer;
1285  CoilWk.append(CWintptsInlayer);
1286  CoilWk.append(CWintptsOutlayer);
1287 
1288  // BabyMEG consists of magnetometers only
1289 // CoilGrad.resize(NumCoil);
1290  CoilGrad.setZero(NumCoil);
1291 // std::cout << CoilGrad << endl;
1292 }
1293 
1294 void RtSssAlgo::getCoilInfoVectorView()
1295 {
1296  Origin.resize(3);
1297  Origin << 0.0, 0.0, 0.04;
1298 
1299  // Load Fiff info
1300  QFile t_fileRaw("/autofs/cluster/fusion/slew/GitHub/mne-cpp/bin/MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
1301  FiffRawData raw(t_fileRaw);
1302 
1303  cout << "number of chs: " << raw.info.nchan << endl;
1304 
1305  // coil type - name
1306  CoilTk.append(QString::number(3012));
1307  CoilTk.append(QString::number(3024));
1308  std::cout << CoilTk[0].toStdString() << " " << CoilTk[1].toStdString() << endl;
1309 
1310  CoilGrad.resize(306);
1311  int nmegcoil=0;
1312  for (int i=0; i<raw.info.nchan; i++ )
1313  {
1314  // kind = 1(MEG), 2(EEG), 3(STI...)
1315  // coil_type = 3012(planar gradiometer), 3024(magnetometer)
1316  if ((raw.info.chs[i].kind == 1 ) && ((raw.info.chs[i].coil_type == 3012) || (raw.info.chs[i].coil_type == 3024)))
1317  {
1318  CoilName.append(QString::number(raw.info.chs[i].coil_type));
1319 // std::cout << "coil no[" << nmegcoil <<"] "<< CoilName[nmegcoil].toStdString() << " ";
1320 
1321  CoilT.append(raw.info.chs[i].coil_trans);
1322 
1323  if (raw.info.chs[i].coil_type == 3012) CoilGrad(nmegcoil) = 1;
1324  else CoilGrad(nmegcoil) = 0;
1325 // std::cout << CoilGrad(nmegcoil) << " ,";
1326 
1327  nmegcoil++;
1328  }
1329  }
1330  cout << endl << "meg coil cnt: " << nmegcoil << endl;
1331  NumCoil = nmegcoil;
1332 
1333  // Number of integration points for each coil type: 8 pts(3012), 9 pts(3024)
1334  CoilNk.resize(2);
1335  CoilNk << 8, 9;
1336 
1337  // Coordinate of coil integrating points: 8 pts(3012), 9 pts(3024)
1338  MatrixXd Cintpts8(3,8), Cintpts9(3,9);
1339  Cintpts8 << 0.0059, -0.0059, 0.0059, -0.0059, 0.0108, -0.0108, 0.0108, -0.0108, 0.0067, 0.0067, -0.0067, -0.0067, 0.0067, 0.0067, -0.0067, -0.0067, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003;
1340  Cintpts9 << 0.0000, 0.0100, -0.0100, 0.0100, -0.0100, 0.0000, 0.0000, 0.0100, -0.0100, 0.0000, 0.0100, 0.0100, -0.0100, -0.0100, 0.0100, -0.0100, 0.0000, 0.0000, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003;
1341  CoilRk.append(Cintpts8);
1342  CoilRk.append(Cintpts9);
1343 
1344  // Coil point weight: 8 pts(3012), 9 pts(3024)
1345  MatrixXd CWintpts8(1,8), CWintpts9(1,9);
1346  CWintpts8 << 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790;
1347  CWintpts9 << 0.1975, 0.0772, 0.0772, 0.0772, 0.0772, 0.1235, 0.1235, 0.1235, 0.1235;
1348  CoilWk.append(CWintpts8);
1349  CoilWk.append(CWintpts9);
1350 
1351  std::cout << "loading MEGData ....";
1352  ifstream inMEGData("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGData.txt",ios::in);
1353  MEGData.resize(NumCoil,1);
1354  VectorXd tmpvec(NumCoil);
1355  for(int k=0; k<NumCoil; k++) inMEGData >> tmpvec(k);
1356  MEGData.col(0) = tmpvec;
1357  inMEGData.close();
1358 // std::cout << MEGData << endl;
1359  std::cout << " finished !" << endl;
1360 }
1361 
1362 // Read FIFF coil configuration and assign them to the Coil class variables
1363 // such as CoilName,CoilT,CoilTk, CoilNk, CoilRk, CoilWk, CoilGrad.
1364 void RtSssAlgo::getCoilInfoVectorView4Sim()
1365 {
1366  QString str;
1367 
1368  NumCoil = 306;
1369  cout << "Num of Coil = " << NumCoil << endl;
1370 
1371 // simulation data from test1_sim.mat
1372  Origin.resize(3);
1373  Origin << 0.0, 0.0, 0.04;
1374 
1375  std::cout << "loading MEGIn ....";
1376  ifstream inMEGIn("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGIn.txt",ios::in);
1377  MEGIn.resize(NumCoil);
1378  for (int k=0; k<NumCoil; k++) inMEGIn >> MEGIn(k);
1379  inMEGIn.close();
1380 // std::cout << MEGIn << endl;
1381  std::cout << " finished !" << endl;
1382 
1383  std::cout << "loading MEGOut ....";
1384  ifstream inMEGOut("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGOut.txt",ios::in);
1385  MEGOut.resize(NumCoil);
1386  for(int k=0; k<NumCoil; k++) inMEGOut >> MEGOut(k);
1387  inMEGOut.close();
1388 // std::cout << MEGOut << endl;
1389  std::cout << " finished !" << endl;
1390 
1391  std::cout << "loading MEGNoise ....";
1392  ifstream inMEGNoise("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGNoise.txt",ios::in);
1393  MEGNoise.resize(NumCoil);
1394  for(int k=0; k<NumCoil; k++) inMEGNoise >> MEGNoise(k);
1395  inMEGNoise.close();
1396 // std::cout << MEGNoise << endl;
1397  std::cout << " finished !" << endl;
1398 
1399  std::cout << "loading MEGData ....";
1400  ifstream inMEGData("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGData.txt",ios::in);
1401  MEGData.resize(NumCoil,1);
1402  VectorXd tmpvec(NumCoil);
1403  for(int k=0; k<NumCoil; k++) inMEGData >> tmpvec(k);
1404  MEGData.col(0) = tmpvec;
1405  inMEGData.close();
1406 // std::cout << MEGData << endl;
1407  std::cout << " finished !" << endl;
1408 
1409  // name per coil
1410  for (int i=0; i<NumCoil; i++)
1411  {
1412  if (i%3 == 0)
1413  CoilName.append("VV_PLANAR_T1");
1414  else if (i%3 ==1)
1415  CoilName.append("VV_PLANAR_T1");
1416  else
1417  CoilName.append("VV_MAG_T1");
1418  }
1419 // std:cout << CoilName[0] << ", " << CoilName[1] << ", " << CoilName[2] << endl;
1420 
1421  // transformation matrix per coil
1422  std::cout << "loading coil transformation matrix ....";
1423  ifstream inCT("/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_coil_def/VectorView/CoilT.txt",ios::in);
1424  MatrixXd tmat(4,4);
1425  for(int k=0; k<NumCoil; k++)
1426  {
1427  for(int i=0; i<4; i++)
1428  for(int j=0; j<4; j++)
1429  {
1430  inCT >> tmat(i,j);
1431  }
1432  CoilT.append(tmat);
1433  }
1434  inCT.close();
1435 // std::cout << CoilT[305] << endl;
1436  std::cout << " finished !" << endl;
1437 
1438  // coil type - name
1439  CoilTk.append("VV_PLANAR_T1");
1440  CoilTk.append("VV_MAG_T1");
1441 // std::cout << CoilTk[0] << ", " << CoilTk[1] << endl;
1442 
1443 
1444  // coil type - Number of integration points for each coil type
1445  CoilNk.resize(2);
1446  CoilNk << 8, 9;
1447 
1448  // coil type - Coordinate of integrating points
1449  MatrixXd Cintpts8(3,8), Cintpts9(3,9);
1450  Cintpts8 << 0.0059, -0.0059, 0.0059, -0.0059, 0.0108, -0.0108, 0.0108, -0.0108, 0.0067, 0.0067, -0.0067, -0.0067, 0.0067, 0.0067, -0.0067, -0.0067, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003;
1451  Cintpts9 << 0.0000, 0.0100, -0.0100, 0.0100, -0.0100, 0.0000, 0.0000, 0.0100, -0.0100, 0.0000, 0.0100, 0.0100, -0.0100, -0.0100, 0.0100, -0.0100, 0.0000, 0.0000, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003, 0.0003;
1452  CoilRk.append(Cintpts8);
1453  CoilRk.append(Cintpts9);
1454 
1455  // Coil type - Coil point weight
1456  MatrixXd CWintpts8(1,8), CWintpts9(1,9);
1457  CWintpts8 << 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790;
1458  CWintpts9 << 0.1975, 0.0772, 0.0772, 0.0772, 0.0772, 0.1235, 0.1235, 0.1235, 0.1235;
1459  CoilWk.append(CWintpts8);
1460  CoilWk.append(CWintpts9);
1461 
1462  // Gradiometer ?
1463  CoilGrad.resize(NumCoil);
1464  CoilGrad << 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0;
1465 // std::cout << CoilGrad << endl;
1466 }
1467 
Contains the declaration of the RtSssAlgo class.
QSharedPointer< FiffInfo > SPtr
Definition: fiff_info.h:99
FIFF raw measurement data.
Definition: fiff_raw_data.h:94