36 #include "rtsssalgo_test.h"
40 RtSssAlgoTest::RtSssAlgoTest()
45 RtSssAlgoTest::~RtSssAlgoTest()
50 QList<MatrixXd> RtSssAlgoTest::buildLinearEqn()
52 QList<MatrixXd> Eqn, EqnRR;
53 QList<MatrixXd> LinEqn;
54 int LInRR, LOutRR, LIn, LOut;
60 if (MACHINE_TYPE == 1)
63 std::cout <<
"loading coil information (VectorView).....";
64 getCoilInfoVectorView4Sim();
65 std::cout <<
"loading coil information (VectorView)..... ** finished **" << endl;
67 else if (MACHINE_TYPE == 2)
70 std::cout <<
"loading coil information (BabyMEG).....";
71 getCoilInfoBabyMEG4Sim();
72 std::cout <<
"loading coil information (BabyMEG)..... ** finisehd **" << endl;
81 EqnRR = getSSSEqn(LInRR, LOutRR);
85 Eqn = getSSSEqn(LIn, LOut);
95 EqnARR.resize(NumCoil, EqnInRR.cols()+EqnOutRR.cols());
96 EqnA.resize(NumCoil, EqnIn.cols()+EqnOut.cols());
97 EqnB.resize(NumCoil,1);
99 CoilScale.setOnes(NumCoil);
100 for(
int i=0; i<NumCoil; i++)
101 if (CoilGrad(i) == 0) CoilScale(i) = MagScale;
103 EqnARR << EqnInRR, EqnOutRR;
104 EqnA << EqnIn, EqnOut;
106 EqnARR = CoilScale.asDiagonal() * EqnARR;
107 EqnA = CoilScale.asDiagonal() * EqnA;
108 EqnB = CoilScale.asDiagonal() * MEGData;
112 LinEqn.append(EqnIn);
113 LinEqn.append(EqnOut);
114 LinEqn.append(EqnARR);
126 std::cout <<
"building SSS linear equation .....finished !" << endl;
133 void RtSssAlgoTest::getCoilInfoVectorView4Sim()
138 cout <<
"Num of Coil = " << NumCoil << endl;
142 Origin << 0.0, 0.0, 0.04;
144 std::cout <<
"loading MEGIn ....";
145 ifstream inMEGIn(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGIn.txt",ios::in);
146 MEGIn.resize(NumCoil);
147 for (
int k=0; k<NumCoil; k++) inMEGIn >> MEGIn(k);
150 std::cout <<
" finished !" << endl;
152 std::cout <<
"loading MEGOut ....";
153 ifstream inMEGOut(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGOut.txt",ios::in);
154 MEGOut.resize(NumCoil);
155 for(
int k=0; k<NumCoil; k++) inMEGOut >> MEGOut(k);
158 std::cout <<
" finished !" << endl;
160 std::cout <<
"loading MEGNoise ....";
161 ifstream inMEGNoise(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGNoise.txt",ios::in);
162 MEGNoise.resize(NumCoil);
163 for(
int k=0; k<NumCoil; k++) inMEGNoise >> MEGNoise(k);
166 std::cout <<
" finished !" << endl;
168 std::cout <<
"loading MEGData ....";
169 ifstream inMEGData(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/VectorView/MEGData.txt",ios::in);
170 MEGData.resize(NumCoil);
171 for(
int k=0; k<NumCoil; k++) inMEGData >> MEGData(k);
174 std::cout <<
" finished !" << endl;
177 for (
int i=0; i<NumCoil; i++)
180 CoilName.append(
"VV_PLANAR_T1");
182 CoilName.append(
"VV_PLANAR_T1");
184 CoilName.append(
"VV_MAG_T1");
189 std::cout <<
"loading coil transformation matrix ....";
190 ifstream inCT(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_coil_def/VectorView/CoilT.txt",ios::in);
192 for(
int k=0; k<NumCoil; k++)
194 for(
int i=0; i<4; i++)
195 for(
int j=0; j<4; j++)
203 std::cout <<
" finished !" << endl;
206 CoilTk.append(
"VV_PLANAR_T1");
207 CoilTk.append(
"VV_MAG_T1");
216 MatrixXd Cintpts8(3,8), Cintpts9(3,9);
217 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;
218 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;
219 CoilRk.append(Cintpts8);
220 CoilRk.append(Cintpts9);
223 MatrixXd Wintpts8(1,8), Wintpts9(1,9);
224 Wintpts8 << 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790, 14.9790, -14.9790;
225 Wintpts9 << 0.1975, 0.0772, 0.0772, 0.0772, 0.0772, 0.1235, 0.1235, 0.1235, 0.1235;
226 CoilWk.append(Wintpts8);
227 CoilWk.append(Wintpts9);
230 CoilGrad.resize(NumCoil);
231 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;
237 void RtSssAlgoTest::getCoilInfoBabyMEG4Sim()
240 int NumCoilOut = 105;
242 NumCoil = NumCoilIn + NumCoilOut;
243 cout <<
"Num of Coil = " << NumCoil << endl;
247 Origin << 0.0, 0.0, 0.04;
249 std::cout <<
"loading MEGIn ....";
250 ifstream inMEGIn(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGIn2.txt",ios::in);
251 MEGIn.resize(NumCoil);
252 for (
int k=0; k<NumCoil; k++) inMEGIn >> MEGIn(k);
255 std::cout <<
" finished !" << endl;
257 std::cout <<
"loading MEGOut ....";
258 ifstream inMEGOut(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGOut2.txt",ios::in);
259 MEGOut.resize(NumCoil);
260 for(
int k=0; k<NumCoil; k++) inMEGOut >> MEGOut(k);
263 std::cout <<
" finished !" << endl;
265 std::cout <<
"loading MEGNoise ....";
266 ifstream inMEGNoise(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGNoise.txt",ios::in);
267 MEGNoise.resize(NumCoil);
268 for(
int k=0; k<NumCoil; k++) inMEGNoise >> MEGNoise(k);
271 std::cout <<
" finished !" << endl;
273 std::cout <<
"loading MEGData ....";
274 ifstream inMEGData(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_test_sim/BabyMEG/MEGData2.txt",ios::in);
275 MEGData.resize(NumCoil);
276 for(
int k=0; k<NumCoil; k++) inMEGData >> MEGData(k);
279 std::cout <<
" finished !" << endl;
282 for (
int i=0; i<NumCoilIn; i++)
283 CoilName.append(
"BSQ2_In_Layer");
284 for (
int i=0; i<NumCoilOut; i++)
285 CoilName.append(
"BSQ2_Out_Layer");
290 std::cout <<
"loading coil transformation matrix ....";
291 ifstream inCT(
"/autofs/cluster/fusion/slew/MEG_realtime/mne-matlab-master/matlab/DATA_coil_def/BabyMEG/CoilT.txt",ios::in);
293 for(
int k=0; k<NumCoil; k++)
295 for(
int i=0; i<4; i++)
296 for(
int j=0; j<4; j++)
304 std::cout <<
" finished !" << endl;
308 CoilTk.append(
"BSQ2_In_Layer");
309 CoilTk.append(
"BSQ2_Out_Layer");
318 MatrixXd CintptsInlayer(3,7), CintptsOutlayer(3,7);
319 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;
320 CintptsOutlayer = CintptsInlayer;
321 CoilRk.append(CintptsInlayer);
322 CoilRk.append(CintptsOutlayer);
325 MatrixXd WintptsInlayer(1,7), WintptsOutlayer(1,7);
326 WintptsInlayer << 0.25, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125;
327 WintptsOutlayer = WintptsInlayer;
328 CoilWk.append(WintptsInlayer);
329 CoilWk.append(WintptsOutlayer);
333 CoilGrad.setZero(NumCoil);
360 QList<MatrixXd> RtSssAlgoTest::getSSSEqn(
int LIn,
int LOut)
365 VectorXd coil_distance, coil_clocation(4,1), coil_vector(4,1);
366 MatrixXd coil_location;
367 MatrixXd EqnIn, EqnOut;
371 NumBIn = (LIn*LIn) + 2*LIn;
372 NumBOut = (LOut*LOut) + 2*LOut;
375 coil_distance.resize(NumCoil);
377 for(
int i = 0; i<NumCoil; i++)
379 coil_clocation = CoilT[i].block(0,3,3,1);
380 coil_distance(i) = (coil_clocation-Origin).norm();
383 RScale = exp(coil_distance.array().log().mean());
386 EqnIn.setZero(NumCoil,NumBIn);
387 EqnOut.setZero(NumCoil,NumBOut);
389 for(
int i = 0; i<NumCoil; i++)
391 if (CoilName[i] == CoilTk[0])
393 else if (CoilName[i] == CoilTk[1])
397 cout<<
"There is no matching coil" << endl;
402 coil_vector = CoilT[i].block(0,2,3,1);
405 int NumCoilPts = CoilNk(coil_index);
406 MatrixXd tmpmat; tmpmat.setOnes(4,NumCoilPts);
408 tmpmat.topRows(3) = CoilRk[coil_index];
409 coil_location = CoilT[i] * tmpmat;
410 tmpmat = coil_location.topRows(3); coil_location.resize(3,NumCoilPts); coil_location = tmpmat;
411 coil_location = coil_location - Origin.replicate(1,NumCoilPts);
412 coil_location = coil_location / RScale;
415 getSSSBasis(coil_location.row(0).transpose(), coil_location.row(1).transpose(), coil_location.row(2).transpose(), LIn, LOut);
419 MatrixXd b_in, b_out;
420 b_in = coil_vector(0)*BInX + coil_vector(1)*BInY + coil_vector(2)*BInZ;
421 b_out = coil_vector(0)*BOutX + coil_vector(1)*BOutY + coil_vector(2)*BOutZ;
425 EqnIn.block(i,0,1,NumBIn) = CoilWk[coil_index] * b_in;
426 EqnOut.block(i,0,1,NumBOut) = CoilWk[coil_index] * b_out;
455 void RtSssAlgoTest::getSSSBasis(VectorXd X, VectorXd Y, VectorXd Z,
int LIn,
int LOut)
457 int NumSample, NumBIn, NumBOut;
460 QList<MatrixXd> P, dP_dx;
461 MatrixXd cur_p, cur_p1, cur_p2;
463 QList<MatrixXcd> YP, dY_dTHETA;
466 MatrixXd BInR, BInPHI, BInTHETA;
469 VectorXcd cur_R, cur_PHI, cur_THETA;
470 MatrixXd BOutR, BOutPHI, BOutTHETA;
474 for(
int i=0; i<X.rows(); i++)
475 if ( (abs(X(i)) < 1e-30) && (abs(Y(i)) < 1e-30) )
477 std::cout <<
"Zero THETA detected!**** ";
478 std::cout <<
"X, Y: " << X(i) <<
", " << Y(i) << endl;
483 LMax = qMax(LIn,LOut);
484 NumSample = X.size();
485 NumBIn = (LIn*LIn) + 2*LIn;
486 NumBOut = (LOut*LOut) + 2*LOut;
489 getCartesianToSpherCoordinate(X,Y,Z);
494 for(
int l=1; l<=LMax; l++)
496 cur_p = legendre(l, THETA.array().cos());
497 P.append(cur_p.transpose());
506 for(
int l=0; l<LMax; l++)
508 dP_dx.append(MatrixXd::Zero(NumSample,(l+1)+1));
509 for(
int m=0; m <= l+1; m++)
511 cur_p1 = P[l].col(m);
515 cur_p2 = -factorial((l+1)-1)/factorial((l+1)+1) * P[l].col(1);
519 cur_p2 = P[l].col(m-1);
527 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);
534 for(
int l=0; l<LMax; l++)
536 YP.append(MatrixXcd::Zero(NumSample,(l+1)+1));
537 dY_dTHETA.append(MatrixXcd::Zero(NumSample,(l+1)+1));
539 for(
int m=0; m <=l+1; m++)
544 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();
547 YP[l].col(m) = cur_phi.array() * P[l].col(m).array().cast<cplxd >();
548 dY_dTHETA[l].col(m) = -cur_phi.array() * THETA.array().sin().cast<cplxd >() * dP_dx[l].col(m).array().cast<cplxd >();
555 BInR.setZero(NumSample,NumBIn);
556 BInPHI.setZero(NumSample,NumBIn);
557 BInTHETA.setZero(NumSample,NumBIn);
559 for(
int l=0; l<LIn; l++)
561 scale_r= R.array().pow((l+1)+2);
563 BInR.col(basis_index) = -((l+1)+1) * YP[l].col(0).array().real() / scale_r.array();
564 BInPHI.col(basis_index) = VectorXd::Zero(NumSample);
565 BInTHETA.col(basis_index) = dY_dTHETA[l].col(0).array().real() / scale_r.array();
566 basis_index = basis_index + 1;
570 for(
int m=0; m<=l; m++)
572 cur_R = -sqrt(2) * ((l+1)+1) * YP[l].col(m+1).array() / scale_r.array().cast<cplxd >();
573 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 >();
574 cur_THETA = sqrt(2) * dY_dTHETA[l].col(m+1).array() / scale_r.array().cast<cplxd >();
577 BInR.col(basis_index) = cur_R.real();
578 BInPHI.col(basis_index) = cur_PHI.real();
579 BInTHETA.col(basis_index) = cur_THETA.real();
580 basis_index = basis_index + 1;
582 BInR.col(basis_index) = cur_R.imag();
583 BInPHI.col(basis_index) = cur_PHI.imag();
584 BInTHETA.col(basis_index) = cur_THETA.imag();
585 basis_index = basis_index + 1;
593 BOutR.setZero(NumSample,NumBOut);
594 BOutPHI.setZero(NumSample,NumBOut);
595 BOutTHETA.setZero(NumSample,NumBOut);
597 for(
int l=0; l<LOut; l++)
599 scale_r= R.array().pow((l+1)-1);
601 BOutR.col(basis_index) = (l+1) * YP[l].col(0).array().real() * scale_r.array();
602 BOutPHI.col(basis_index) = VectorXd::Zero(NumSample);
603 BOutTHETA.col(basis_index) = dY_dTHETA[l].col(0).array().real() * scale_r.array();
604 basis_index = basis_index + 1;
606 for(
int m=0; m<=l; m++)
608 cur_R = sqrt(2) * (l+1) * YP[l].col(m+1).array() * scale_r.array().cast<cplxd >();
609 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 >();
610 cur_THETA = sqrt(2) * dY_dTHETA[l].col(m+1).array() * scale_r.array().cast<cplxd >();
611 BOutR.col(basis_index) = cur_R.real();
612 BOutPHI.col(basis_index) = cur_PHI.real();
613 BOutTHETA.col(basis_index) = cur_THETA.real();
614 basis_index = basis_index + 1;
615 BOutR.col(basis_index) = cur_R.imag();
616 BOutPHI.col(basis_index) = cur_PHI.imag();
617 BOutTHETA.col(basis_index) = cur_THETA.imag();
618 basis_index = basis_index + 1;
627 getSphereToCartesianVector();
628 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();
629 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();
630 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();
631 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();
632 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();
633 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();
665 void RtSssAlgoTest::getCartesianToSpherCoordinate(VectorXd X, VectorXd Y, VectorXd Z)
671 hypotxy = hypot(X,Y);
672 R = hypot(hypotxy,Z);
674 THETA = atan2vec(hypotxy,Z);
706 void RtSssAlgoTest::getSphereToCartesianVector()
708 int NumExp = R.size();
710 R_X = THETA.array().sin() * PHI.array().cos();
711 R_Y = THETA.array().sin() * PHI.array().sin();
712 R_Z = THETA.array().cos();
714 PHI_X = -PHI.array().sin();
715 PHI_Y = PHI.array().cos();
716 PHI_Z = VectorXd::Zero(NumExp,1);
719 THETA_X = THETA.array().cos() * PHI.array().cos();
720 THETA_Y = THETA.array().cos() * PHI.array().sin();
721 THETA_Z = -THETA.array().sin();
756 QList<MatrixXd> RtSssAlgoTest::getSSSRR(MatrixXd EqnIn, MatrixXd EqnOut, MatrixXd EqnARR, MatrixXd EqnA, MatrixXd EqnB)
758 int NumBIn, NumBOut, NumCoil, NumExp;
759 MatrixXd EqnRRInv, EqnInv, SSSIn, SSSOut, Weight;
761 double RR_K1, RR_K2, RR_K3;
762 double eqn_scale0, eqn_scale;
763 MatrixXd sol_X, sol_X_old, eqn_Y, eqn_D, temp_M, temp_N, sol_in, sol_out;
765 VectorXd eqn_err, weight_index;
766 QList<MatrixXd> RRsss;
769 double ErrTolRel = 1e-3;
772 double WeightThres = 1 - 1e-6;
775 NumBIn = EqnIn.cols();
776 NumBOut = EqnOut.cols();
777 NumCoil = EqnB.rows();
778 NumExp = EqnB.cols();
823 EqnRRInv = (EqnARR.transpose() * EqnARR).inverse();
824 EqnInv = (EqnA.transpose() * EqnA).inverse();
826 SSSIn.setZero(NumCoil,NumExp);
827 SSSOut.setZero(NumCoil,NumExp);
828 Weight.setZero(NumCoil,NumExp);
829 ErrRel.setZero(NumExp);
832 RR_K1 = qSqrt(1-qSqrt(3)/2) * RR_K2;
834 for(
int i=0; i<NumExp; i++)
837 sol_X = EqnRRInv * (EqnARR.transpose() * EqnB.col(i));
841 eqn_err = EqnARR * sol_X - EqnB.col(i);
842 eqn_scale0 = stdev(eqn_err);
843 eqn_err = eqn_err.cwiseAbs() / eqn_scale0;
846 sol_X_old.setConstant(sol_X.rows(), sol_X.cols(), 1e30);
849 while (((sol_X.array()-sol_X_old.array()).matrix().norm() / sol_X.norm()) > ErrTolRel)
856 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);
860 weight_index = eigen_LT_index(Weight.col(i), WeightThres);
865 eqn_Y.resize(weight_index.size(), EqnARR.cols());
866 eqn_D.resize(weight_index.size(), 1);
867 for(
int k=0; k<weight_index.size(); k++)
869 eqn_Y.row(k) = EqnARR.row(weight_index(k));
872 eqn_D(k) = Weight(weight_index(k),i) - 1;
876 temp_M = EqnARR.transpose() * (Weight.col(i).array() * EqnB.col(i).array()).matrix();
884 temp_N = EqnRRInv * eqn_Y.transpose();
893 diagMat = (1 / eqn_D.array()).matrix().asDiagonal();
896 sol_X = EqnRRInv * temp_M - temp_N * (diagMat + eqn_Y * temp_N).inverse() * (temp_N.transpose() * temp_M);
897 eqn_err = (EqnARR * sol_X - EqnB.col(i)).cwiseAbs();
898 eqn_scale = qMin(eqn_scale0, RR_K3 * qSqrt((Weight.col(i).array() * eqn_err.array() * eqn_err.array()).mean()));
899 eqn_err = eqn_err / eqn_scale;
909 eqn_Y.resize(weight_index.size(), NumBIn+NumBOut);
910 for(
int k=0; k<weight_index.size(); k++) eqn_Y.row(k) = EqnA.row(weight_index(k));
911 temp_M = EqnA.transpose() * (Weight.col(i).array() * EqnB.col(i).array()).matrix();
912 temp_N = EqnInv * eqn_Y.transpose();
917 diagMat = (1 / eqn_D.array()).matrix().asDiagonal();
918 sol_X = EqnInv * temp_M - temp_N * (diagMat + eqn_Y * temp_N).inverse() * (temp_N.transpose() * temp_M);
920 ErrRel(i) = (EqnA * sol_X - EqnB.col(i)).norm() / EqnB.col(i).norm();
921 sol_in = sol_X.block(0,0,NumBIn,1);
922 sol_out = sol_X.block(NumBIn,0,NumBOut,1);
925 SSSIn.col(i) = EqnIn * sol_in;
926 SSSOut.col(i) = EqnOut * sol_out;
931 RRsss.append(SSSOut);
932 RRsss.append(Weight);
933 RRsss.append(ErrRel);
961 QList<MatrixXd> RtSssAlgoTest::getSSSOLS(MatrixXd EqnIn, MatrixXd EqnOut, MatrixXd EqnA, MatrixXd EqnB)
963 int NumBIn, NumBOut, NumCoil, NumExp;
964 MatrixXd EqnRRInv, EqnInv, SSSIn, SSSOut;
965 MatrixXd sol_X, sol_in, sol_out;
967 QList<MatrixXd> OLSsss;
970 NumBIn = EqnIn.cols();
971 NumBOut = EqnOut.cols();
972 NumCoil = EqnB.rows();
973 NumExp = EqnB.cols();
1020 EqnRRInv = (EqnA.transpose() * EqnA).inverse();
1021 SSSIn.setZero(NumCoil,NumExp);
1022 SSSOut.setZero(NumCoil,NumExp);
1023 ErrRel.setZero(NumExp);
1025 for (
int i=0; i<NumExp; i++)
1028 sol_X = EqnRRInv * (EqnA.transpose() * EqnB.col(i));
1029 ErrRel(i) = (EqnA * sol_X - EqnB.col(i)).norm() / EqnB.col(i).norm();
1030 sol_in = sol_X.block(0,0,NumBIn,1);
1031 sol_out = sol_X.block(NumBIn,0, NumBOut,1);
1034 SSSIn.col(i) = EqnIn * sol_in;
1035 SSSOut.col(i) = EqnOut * sol_out;
1037 OLSsss.append(SSSIn);
1038 OLSsss.append(SSSOut);
1039 OLSsss.append(ErrRel);
1044 QList<MatrixXd> RtSssAlgoTest::getLinEqn()
1046 QList<MatrixXd> LinEqn;
1048 LinEqn.append(EqnIn);
1049 LinEqn.append(EqnOut);
1050 LinEqn.append(EqnARR);
1051 LinEqn.append(EqnA);
1052 LinEqn.append(EqnB);
Contains the declaration of the RtSssAlgo class.