41 #include <QtConcurrent/QtConcurrentMap>
44 RtSssAlgo::RtSssAlgo()
49 RtSssAlgo::~RtSssAlgo()
55 MatrixXd RtSssAlgo::buildLinearEqn()
57 QList<MatrixXd> Eqn, EqnRR;
58 QList<MatrixXd> LinEqn;
84 EqnRR = getSSSEqn(LIn, LOut);
90 Eqn = getSSSEqn(LIn, LOut);
114 EqnARR.resize(NumCoil, EqnInRR.cols()+EqnOutRR.cols());
115 EqnA.resize(NumCoil, EqnIn.cols()+EqnOut.cols());
116 EqnB.resize(NumCoil,1);
122 if ((0 < CoilGrad.sum()) && (CoilGrad.sum() < NumCoil)) MagScale = 100;
126 CoilScale.setOnes(NumCoil);
127 for(
int i=0; i<NumCoil; i++)
129 if (CoilGrad(i) == 0) CoilScale(i) = MagScale;
133 EqnARR << EqnInRR, EqnOutRR;
134 EqnA << EqnIn, EqnOut;
136 EqnARR = CoilScale.asDiagonal() * EqnARR;
137 EqnA = CoilScale.asDiagonal() * EqnA;
146 LinEqn.append(EqnIn);
147 LinEqn.append(EqnOut);
148 LinEqn.append(EqnARR);
151 LinEqn.append(CoilScale.asDiagonal());
165 return CoilScale.asDiagonal();
168 void RtSssAlgo::setSSSParameter(QList<int> expansionOrder)
175 LInRR = expansionOrder[0];
176 LOutRR = expansionOrder[1];
177 LInOLS = expansionOrder[2];
178 LOutOLS = expansionOrder[3];
186 Origin << 0.0, 0.0, 0.04;
197 BadChan.setZero(fiffInfo->nchan);
198 for (qint32 i=0; i<fiffInfo->nchan; ++i)
199 if(fiffInfo->chs[i].kind == FIFFV_MEG_CH)
201 for(qint32 j = 0; j < fiffInfo->bads.size(); j++)
202 if(fiffInfo->chs[i].ch_name == fiffInfo->bads[j])
205 qDebug() <<
"bad ch: " << i <<
"(" << fiffInfo->chs[i].ch_name <<
")";
212 NumCoil = NumMEGChan - NumBadCoil;
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;
222 CoilNk.resize(NumCoil);
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;
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;
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;
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;
248 CoilGrad.setZero(NumCoil);
250 for (qint32 i=0; i<fiffInfo->nchan; ++i)
252 if(fiffInfo->chs[i].kind == FIFFV_MEG_CH && BadChan(i) == 0)
255 CoilT.append(fiffInfo->chs[i].coil_trans);
256 CoilName.append(fiffInfo->chs[i].ch_name);
260 switch (fiffInfo->chs[i].coil_type)
265 CoilRk.append(C3012intpts8);
266 CoilWk.append(C3012Wintpts8);
272 CoilRk.append(C3024intpts9);
273 CoilWk.append(C3024Wintpts9);
279 CoilRk.append(C7002intpts7);
280 CoilWk.append(C7002Wintpts7);
286 CoilRk.append(C7003intpts4);
287 CoilWk.append(C7003Wintpts4);
293 qDebug() <<
" This coil type is NOT supported";
333 QList<MatrixXd> RtSssAlgo::getSSSEqn(qint32 LIn, qint32 LOut)
339 VectorXd coil_distance, coil_clocation(4,1), coil_vector(4,1);
340 MatrixXd coil_location;
341 MatrixXd EqnIn, EqnOut;
349 NumBIn = (LIn*LIn) + 2*LIn;
350 NumBOut = (LOut*LOut) + 2*LOut;
353 coil_distance.resize(NumCoil);
355 for(
int i = 0; i<NumCoil; i++)
357 coil_clocation = CoilT[i].block(0,3,3,1);
358 coil_distance(i) = (coil_clocation-Origin).norm();
361 RScale = exp(coil_distance.array().log().mean());
364 EqnIn.setZero(NumCoil,NumBIn);
365 EqnOut.setZero(NumCoil,NumBOut);
367 for(
int i = 0; i<NumCoil; i++)
380 coil_vector = CoilT[i].block(0,2,3,1);
384 qint32 NumCoilPts = CoilNk(i);
385 MatrixXd tmpmat; tmpmat.setOnes(4,NumCoilPts);
389 tmpmat.topRows(3) = CoilRk[i];
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;
397 getSSSBasis(coil_location.row(0).transpose(), coil_location.row(1).transpose(), coil_location.row(2).transpose(), LIn, LOut);
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;
409 EqnIn.block(i,0,1,NumBIn) = CoilWk[i] * b_in;
410 EqnOut.block(i,0,1,NumBOut) = CoilWk[i] * b_out;
439 void RtSssAlgo::getSSSBasis(VectorXd X, VectorXd Y, VectorXd Z, qint32 LIn, qint32 LOut)
441 int NumSample, NumBIn, NumBOut;
444 QList<MatrixXd> P, dP_dx;
445 MatrixXd cur_p, cur_p1, cur_p2;
447 QList<MatrixXcd> YP, dY_dTHETA;
450 MatrixXd BInR, BInPHI, BInTHETA;
453 VectorXcd cur_R, cur_PHI, cur_THETA;
454 MatrixXd BOutR, BOutPHI, BOutTHETA;
458 for(
int i=0; i<X.rows(); i++)
459 if ( (abs(X(i)) < 1e-30) && (abs(Y(i)) < 1e-30) )
461 std::cout <<
"Zero THETA detected!**** ";
462 std::cout <<
"X, Y: " << X(i) <<
", " << Y(i) << endl;
467 LMax = qMax(LIn,LOut);
468 NumSample = X.size();
469 NumBIn = (LIn*LIn) + 2*LIn;
470 NumBOut = (LOut*LOut) + 2*LOut;
473 getCartesianToSpherCoordinate(X,Y,Z);
478 for(
int l=1; l<=LMax; l++)
480 cur_p = legendre(l, THETA.array().cos());
481 P.append(cur_p.transpose());
490 for(
int l=0; l<LMax; l++)
492 dP_dx.append(MatrixXd::Zero(NumSample,(l+1)+1));
493 for(
int m=0; m <= l+1; m++)
495 cur_p1 = P[l].col(m);
499 cur_p2 = -factorial((l+1)-1)/factorial((l+1)+1) * P[l].col(1);
503 cur_p2 = P[l].col(m-1);
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);
518 for(
int l=0; l<LMax; l++)
520 YP.append(MatrixXcd::Zero(NumSample,(l+1)+1));
521 dY_dTHETA.append(MatrixXcd::Zero(NumSample,(l+1)+1));
523 for(
int m=0; m <=l+1; m++)
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();
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 >();
539 BInR.setZero(NumSample,NumBIn);
540 BInPHI.setZero(NumSample,NumBIn);
541 BInTHETA.setZero(NumSample,NumBIn);
543 for(
int l=0; l<LIn; l++)
545 scale_r= R.array().pow((l+1)+2);
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;
554 for(
int m=0; m<=l; m++)
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 >();
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;
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;
577 BOutR.setZero(NumSample,NumBOut);
578 BOutPHI.setZero(NumSample,NumBOut);
579 BOutTHETA.setZero(NumSample,NumBOut);
581 for(
int l=0; l<LOut; l++)
583 scale_r= R.array().pow((l+1)-1);
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;
590 for(
int m=0; m<=l; m++)
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;
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();
648 void RtSssAlgo::getCartesianToSpherCoordinate(VectorXd X, VectorXd Y, VectorXd Z)
654 hypotxy = hypot(X,Y);
655 R = hypot(hypotxy,Z);
657 THETA = atan2vec(hypotxy,Z);
689 void RtSssAlgo::getSphereToCartesianVector()
691 int NumExp = R.size();
693 R_X = THETA.array().sin() * PHI.array().cos();
694 R_Y = THETA.array().sin() * PHI.array().sin();
695 R_Z = THETA.array().cos();
697 PHI_X = -PHI.array().sin();
698 PHI_Y = PHI.array().cos();
699 PHI_Z = VectorXd::Zero(NumExp,1);
702 THETA_X = THETA.array().cos() * PHI.array().cos();
703 THETA_Y = THETA.array().cos() * PHI.array().sin();
704 THETA_Z = -THETA.array().sin();
738 MatrixXd RtSssAlgo::getSSSRR(MatrixXd EqnB)
740 int NumBIn, NumBOut, NumCoil, NumExp;
741 MatrixXd EqnRRInv, EqnInv, SSSIn, SSSOut, Weight;
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;
747 VectorXd eqn_err, weight_index;
748 QList<MatrixXd> RRsss;
751 double ErrTolRel = 1e-3;
754 double WeightThres = 1 - 1e-6;
757 NumBIn = EqnIn.cols();
758 NumBOut = EqnOut.cols();
759 NumCoil = EqnB.rows();
760 NumExp = EqnB.cols();
762 EqnRRInv = (EqnARR.transpose() * EqnARR).inverse();
763 EqnInv = (EqnA.transpose() * EqnA).inverse();
765 SSSIn.setZero(NumCoil,NumExp);
766 SSSOut.setZero(NumCoil,NumExp);
767 Weight.setZero(NumCoil,NumExp);
768 ErrRel.setZero(NumExp);
771 RR_K1 = qSqrt(1-qSqrt(3)/2) * RR_K2;
773 for(
int i=0; i<NumExp; i++)
776 sol_X = EqnRRInv * (EqnARR.transpose() * EqnB.col(i));
780 eqn_err = EqnARR * sol_X - EqnB.col(i);
781 eqn_scale0 = stdev(eqn_err);
782 eqn_err = eqn_err.cwiseAbs() / eqn_scale0;
785 sol_X_old.setConstant(sol_X.rows(), sol_X.cols(), 1e30);
788 while (((sol_X.array()-sol_X_old.array()).matrix().norm() / sol_X.norm()) > ErrTolRel)
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);
800 weight_index = eigen_LT_index(Weight.col(i), WeightThres);
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++)
809 eqn_Y.row(k) = EqnARR.row(weight_index(k));
812 eqn_D(k) = Weight(weight_index(k),i) - 1;
815 temp_M = EqnARR.transpose() * (Weight.col(i).array() * EqnB.col(i).array()).matrix();
823 temp_N = EqnRRInv * eqn_Y.transpose();
832 diagMat = (1 / eqn_D.array()).matrix().asDiagonal();
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;
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();
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);
860 ErrRel(i) = (EqnA * sol_X - EqnB.col(i)).norm() / EqnB.col(i).norm();
862 sol_in = sol_X.block(0,0,NumBIn,1);
863 sol_out = sol_X.block(NumBIn,0,NumBOut,1);
866 SSSIn.col(i) = EqnIn * sol_in;
867 SSSOut.col(i) = EqnOut * sol_out;
871 RRsss.append(SSSOut);
872 RRsss.append(Weight);
873 RRsss.append(ErrRel);
903 MatrixXd RtSssAlgo::getSSSOLS(MatrixXd EqnB)
905 int NumBIn, NumBOut, NumCoil, NumExp;
906 MatrixXd EqnRRInv, EqnInv, SSSIn, SSSOut;
907 MatrixXd sol_X, sol_in, sol_out;
909 QList<MatrixXd> OLSsss;
912 NumBIn = EqnIn.cols();
913 NumBOut = EqnOut.cols();
914 NumCoil = EqnB.rows();
915 NumExp = EqnB.cols();
917 EqnRRInv = (EqnA.transpose() * EqnA).inverse();
918 SSSIn.setZero(NumCoil,NumExp);
919 SSSOut.setZero(NumCoil,NumExp);
920 ErrRel.setZero(NumExp);
922 for (
int i=0; i<NumExp; i++)
925 sol_X = EqnRRInv * (EqnA.transpose() * EqnB.col(i));
927 ErrRel(i) = (EqnA * sol_X - EqnB.col(i)).norm() / EqnB.col(i).norm();
929 sol_in = sol_X.block(0,0,NumBIn,1);
931 sol_out = sol_X.block(NumBIn,0, NumBOut,1);
934 SSSIn.col(i) = EqnIn * sol_in;
935 SSSOut.col(i) = EqnOut * sol_out;
938 OLSsss.append(SSSIn);
939 OLSsss.append(SSSOut);
940 OLSsss.append(ErrRel);
947 qint32 RtSssAlgo::getNumMEGChanUsed()
952 qint32 RtSssAlgo::getNumMEGChan()
957 VectorXi RtSssAlgo::getBadChan()
962 qint32 RtSssAlgo::getNumMEGBadChan()
975 QList<MatrixXd> RtSssAlgo::getLinEqn()
977 QList<MatrixXd> LinEqn;
979 LinEqn.append(EqnIn);
980 LinEqn.append(EqnOut);
981 LinEqn.append(EqnARR);
990 MatrixXd legendre(
int l, VectorXd x)
995 outP.setZero(l+1, x.size());
998 for(
int m=0; m<=l; m++)
1000 for(
int k=0; k<x.size(); k++)
1001 outP(m,k) = plgndr(l,m,x(k));
1008 float plgndr(
int l,
int m,
float x)
1011 float fact,pll,pmm,pmmp1,somx2;
1014 if (m < 0 || m > l || fabs(x) > 1.0)
1015 std::cout <<
"Bad arguments in routine plgndr";
1020 somx2=sqrt((1.0-x)*(1.0+x));
1031 pmmp1=x*(2*m+1)*pmm;
1036 for (ll=m+2;ll<=l;ll++)
1038 pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m);
1052 VectorXd hypot(VectorXd X, VectorXd Y)
1055 rst = (X.array().abs().pow(2) + Y.array().abs().pow(2)).sqrt();
1062 VectorXd atan2vec(VectorXd a, VectorXd b)
1067 for(
int i=0; i<n; i++)
1069 at2(i) = qAtan2(a(i),b(i));
1075 VectorXd find(MatrixXd M,
int I)
1079 for(
int i=0; i<M.cols(); i++)
1080 for(
int j=0; j<M.rows(); j++)
1083 outIDX(k) = i*M.rows() + j;
1091 double factorial(
int n)
1095 for(
int i=2; i<=n; i++) fac *= i;
1103 double stdev(VectorXd V)
1105 int size = V.size();
1109 for(
int i=0; i<size; i++)
1111 sum += pow((V(i)-ave),2);
1113 return sqrt(sum/size);
1119 VectorXd eigen_LTE(VectorXd V,
double tol)
1122 VectorXd outTrue = VectorXd::Zero(row);
1124 for(
int i=0; i<row; i++)
1125 if(V(i) <= tol) outTrue(i) = 1;
1133 VectorXd eigen_GT(VectorXd V,
double tol)
1136 VectorXd outTrue = VectorXd::Zero(row);
1138 for(
int i=0; i<row; i++)
1139 if(V(i) > tol) outTrue(i) = 1;
1146 VectorXd eigen_AND(VectorXd V1, VectorXd V2)
1148 int row = V1.rows();
1149 VectorXd outTrue = VectorXd::Zero(row);
1151 for(
int i=0; i<row; i++)
1152 if(V1(i) && V2(i)) outTrue(i) = 1;
1159 VectorXd eigen_LT_index(VectorXd V,
double tol)
1162 int size = V.size();
1163 VectorXd outIdx; outIdx.setZero(size);
1165 for(
int i=0; i<size; i++)
1171 return outIdx.head(cnt);
1175 VectorXd eigen_LT_index_test(
int cnt)
1177 VectorXd outIdx = VectorXd::Zero(cnt+1);
1179 for (
int i=0; i<cnt+1; i++)
1192 void RtSssAlgo::getCoilInfoBabyMEG4Sim()
1194 int NumCoilIn = 270;
1195 int NumCoilOut = 105;
1197 NumCoil = NumCoilIn + NumCoilOut;
1198 cout <<
"Num of Coil = " << NumCoil << endl;
1202 Origin << 0.0, 0.0, 0.04;
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);
1210 std::cout <<
" finished !" << endl;
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);
1218 std::cout <<
" finished !" << endl;
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);
1226 std::cout <<
" finished !" << endl;
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;
1236 std::cout <<
" finished !" << endl;
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");
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);
1250 for(
int k=0; k<NumCoil; k++)
1252 for(
int i=0; i<4; i++)
1253 for(
int j=0; j<4; j++)
1261 std::cout <<
" finished !" << endl;
1265 CoilTk.append(
"BSQ2_In_Layer");
1266 CoilTk.append(
"BSQ2_Out_Layer");
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);
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);
1290 CoilGrad.setZero(NumCoil);
1294 void RtSssAlgo::getCoilInfoVectorView()
1297 Origin << 0.0, 0.0, 0.04;
1300 QFile t_fileRaw(
"/autofs/cluster/fusion/slew/GitHub/mne-cpp/bin/MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
1303 cout <<
"number of chs: " << raw.info.nchan << endl;
1306 CoilTk.append(QString::number(3012));
1307 CoilTk.append(QString::number(3024));
1308 std::cout << CoilTk[0].toStdString() <<
" " << CoilTk[1].toStdString() << endl;
1310 CoilGrad.resize(306);
1312 for (
int i=0; i<raw.info.nchan; i++ )
1316 if ((raw.info.chs[i].kind == 1 ) && ((raw.info.chs[i].coil_type == 3012) || (raw.info.chs[i].coil_type == 3024)))
1318 CoilName.append(QString::number(raw.info.chs[i].coil_type));
1321 CoilT.append(raw.info.chs[i].coil_trans);
1323 if (raw.info.chs[i].coil_type == 3012) CoilGrad(nmegcoil) = 1;
1324 else CoilGrad(nmegcoil) = 0;
1330 cout << endl <<
"meg coil cnt: " << nmegcoil << endl;
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);
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);
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;
1359 std::cout <<
" finished !" << endl;
1364 void RtSssAlgo::getCoilInfoVectorView4Sim()
1369 cout <<
"Num of Coil = " << NumCoil << endl;
1373 Origin << 0.0, 0.0, 0.04;
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);
1381 std::cout <<
" finished !" << endl;
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);
1389 std::cout <<
" finished !" << endl;
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);
1397 std::cout <<
" finished !" << endl;
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;
1407 std::cout <<
" finished !" << endl;
1410 for (
int i=0; i<NumCoil; i++)
1413 CoilName.append(
"VV_PLANAR_T1");
1415 CoilName.append(
"VV_PLANAR_T1");
1417 CoilName.append(
"VV_MAG_T1");
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);
1425 for(
int k=0; k<NumCoil; k++)
1427 for(
int i=0; i<4; i++)
1428 for(
int j=0; j<4; j++)
1436 std::cout <<
" finished !" << endl;
1439 CoilTk.append(
"VV_PLANAR_T1");
1440 CoilTk.append(
"VV_MAG_T1");
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);
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);
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;
Contains the declaration of the RtSssAlgo class.
QSharedPointer< FiffInfo > SPtr
FIFF raw measurement data.