67 , m_iNumLeadFieldCombinations(0)
68 , m_ppPairIdxCombinations(NULL)
71 , m_iSamplesStcWindow(-1)
84 , m_iNumLeadFieldCombinations(0)
85 , m_ppPairIdxCombinations(NULL)
88 , m_iSamplesStcWindow(-1)
92 init(p_pFwd, p_bSparsed, p_iN, p_dThr);
111 std::cout <<
"OpenMP enabled" << std::endl;
114 std::cout <<
"OpenMP disabled (to enable it: VS2010->Project Properties->C/C++->Language, then modify OpenMP Support)" << std::endl;
117 std::cout <<
"Available Threats: " <<
m_iMaxNumThreads << std::endl << std::endl;
120 std::cout <<
"##### Initialization RAP MUSIC started ######\n\n";
139 if ( p_pFwd.
sol->data.cols() % 3 != 0 )
141 std::cout <<
"Gain matrix is not associated with a 3D grid!\n";
158 std::cout <<
"Calculate gain matrix combinations. \n";
166 std::cout <<
"Gain matrix combinations calculated. \n\n";
174 std::cout <<
"Number of sources to find: " <<
m_iN <<
"\n\n";
180 std::cout <<
"##### Initialization RAP MUSIC completed ######\n\n\n";
182 Q_UNUSED(p_bSparsed);
210 Q_UNUSED(pick_normal);
216 std::cout <<
"Number of FiffEvoked channels (" << p_fiffEvoked.
data.rows() <<
") doesn't match the number of channels (" <<
m_iNumChannels <<
") of the forward solution." << std::endl;
217 return p_sourceEstimate;
232 p_sourceEstimate.
tmin = p_fiffEvoked.
times[0];
233 p_sourceEstimate.
tstep = p_fiffEvoked.
times[1] - p_fiffEvoked.
times[0];
237 QList< DipolePair<double> > t_RapDipoles;
240 for(qint32 i = 0; i < t_RapDipoles.size(); ++i)
242 double dip1 = sqrt( pow(t_RapDipoles[i].m_Dipole1.phi_x(),2) +
243 pow(t_RapDipoles[i].m_Dipole1.phi_y(),2) +
244 pow(t_RapDipoles[i].m_Dipole1.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
246 double dip2 = sqrt( pow(t_RapDipoles[i].m_Dipole2.phi_x(),2) +
247 pow(t_RapDipoles[i].m_Dipole2.phi_y(),2) +
248 pow(t_RapDipoles[i].m_Dipole2.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
250 RowVectorXd dip1Time = RowVectorXd::Constant(p_fiffEvoked.
data.cols(), dip1);
251 RowVectorXd dip2Time = RowVectorXd::Constant(p_fiffEvoked.
data.cols(), dip2);
253 p_sourceEstimate.
data.block(t_RapDipoles[i].m_iIdx1, 0, 1, p_fiffEvoked.
data.cols()) = dip1Time;
254 p_sourceEstimate.
data.block(t_RapDipoles[i].m_iIdx2, 0, 1, p_fiffEvoked.
data.cols()) = dip2Time;
262 qint32 t_iNumSensors = p_fiffEvoked.
data.rows();
263 qint32 t_iNumSteps = p_fiffEvoked.
data.cols();
266 qint32 t_iSamplesDiscard = t_iSamplesOverlap/2;
268 MatrixXd data = MatrixXd::Zero(t_iNumSensors, m_iSamplesStcWindow);
270 qint32 curSample = 0;
271 qint32 curResultSample = 0;
272 qint32 stcWindowSize = m_iSamplesStcWindow - 2*t_iSamplesDiscard;
276 QList< DipolePair<double> > t_RapDipoles;
279 if(curSample + m_iSamplesStcWindow >= t_iNumSteps)
285 data = p_fiffEvoked.
data.block(0, curSample, t_iNumSensors, m_iSamplesStcWindow);
288 curSample += (m_iSamplesStcWindow - t_iSamplesOverlap);
290 curSample -= t_iSamplesDiscard;
297 stcWindowSize = p_sourceEstimate.
data.cols() - curResultSample;
299 for(qint32 i = 0; i < t_RapDipoles.size(); ++i)
301 double dip1 = sqrt( pow(t_RapDipoles[i].m_Dipole1.phi_x(),2) +
302 pow(t_RapDipoles[i].m_Dipole1.phi_y(),2) +
303 pow(t_RapDipoles[i].m_Dipole1.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
305 double dip2 = sqrt( pow(t_RapDipoles[i].m_Dipole2.phi_x(),2) +
306 pow(t_RapDipoles[i].m_Dipole2.phi_y(),2) +
307 pow(t_RapDipoles[i].m_Dipole2.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
309 RowVectorXd dip1Time = RowVectorXd::Constant(stcWindowSize, dip1);
310 RowVectorXd dip2Time = RowVectorXd::Constant(stcWindowSize, dip2);
312 p_sourceEstimate.
data.block(t_RapDipoles[i].m_iIdx1, curResultSample, 1, stcWindowSize) = dip1Time;
313 p_sourceEstimate.
data.block(t_RapDipoles[i].m_iIdx2, curResultSample, 1, stcWindowSize) = dip2Time;
317 curResultSample += stcWindowSize;
325 return p_sourceEstimate;
337 std::cout <<
"Number of FiffEvoked channels (" << data.rows() <<
") doesn't match the number of channels (" <<
m_iNumChannels <<
") of the forward solution." << std::endl;
338 return p_sourceEstimate;
352 p_sourceEstimate.
times = RowVectorXf::Zero(data.cols());
353 p_sourceEstimate.
times[0] = tmin;
354 for(qint32 i = 1; i < p_sourceEstimate.
times.size(); ++i)
355 p_sourceEstimate.
times[i] = p_sourceEstimate.
times[i-1] + tstep;
356 p_sourceEstimate.
tmin = tmin;
357 p_sourceEstimate.
tstep = tstep;
362 for(qint32 i = 0; i < t_RapDipoles.size(); ++i)
364 double dip1 = sqrt( pow(t_RapDipoles[i].m_Dipole1.phi_x(),2) +
365 pow(t_RapDipoles[i].m_Dipole1.phi_y(),2) +
366 pow(t_RapDipoles[i].m_Dipole1.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
368 double dip2 = sqrt( pow(t_RapDipoles[i].m_Dipole2.phi_x(),2) +
369 pow(t_RapDipoles[i].m_Dipole2.phi_y(),2) +
370 pow(t_RapDipoles[i].m_Dipole2.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
372 RowVectorXd dip1Time = RowVectorXd::Constant(data.cols(), dip1);
373 RowVectorXd dip2Time = RowVectorXd::Constant(data.cols(), dip2);
375 p_sourceEstimate.
data.block(t_RapDipoles[i].m_iIdx1, 0, 1, data.cols()) = dip1Time;
376 p_sourceEstimate.
data.block(t_RapDipoles[i].m_iIdx2, 0, 1, data.cols()) = dip2Time;
379 return p_sourceEstimate;
392 std::cout <<
"RAP MUSIC wasn't initialized!";
393 return p_SourceEstimate;
399 std::cout <<
"Lead Field channels do not fit to number of measurement channels!";
400 return p_SourceEstimate;
416 int t_r =
calcPhi_s(p_matMeasurement, t_pMatPhi_s);
418 int t_iMaxSearch =
m_iN < t_r ?
m_iN : t_r;
422 std::cout <<
"Warning: Rank " << t_r <<
" of the measurement data is smaller than the " <<
m_iN;
423 std::cout <<
" sources to find." << std::endl;
424 std::cout <<
" Searching now for " << t_iMaxSearch <<
" correlated sources.";
425 std::cout << std::endl << std::endl;
431 t_matOrthProj.setIdentity();
435 t_matA_k_1.setZero();
451 p_RapDipoles.clear();
453 std::cout <<
"##### Calculation of RAP MUSIC started ######\n\n";
455 MatrixXT t_matProj_Phi_s(t_matOrthProj.rows(), t_pMatPhi_s->cols());
459 for(
int r = 0; r < t_iMaxSearch ; ++r)
461 t_matProj_Phi_s = t_matOrthProj*(*t_pMatPhi_s);
468 Eigen::JacobiSVD< MatrixXT > t_svdProj_Phi_S(t_matProj_Phi_s, Eigen::ComputeThinU);
470 useFullRank(t_svdProj_Phi_S.matrixU(), t_svdProj_Phi_S.singularValues().asDiagonal(), t_matU_B);
478 clock_t start_subcorr, end_subcorr;
479 start_subcorr = clock();
483 #pragma omp parallel num_threads(m_iMaxNumThreads)
493 MatrixX6T t_matProj_G(t_matProj_LeadField.rows(),6);
522 end_subcorr = clock();
524 float t_fSubcorrElapsedTime = ( (float)(end_subcorr-start_subcorr) / (float)CLOCKS_PER_SEC ) * 1000.0f;
525 std::cout <<
"Time Elapsed: " << t_fSubcorrElapsedTime <<
" ms" << std::endl;
530 VectorXT::Index t_iMaxIdx;
532 t_val_roh_k = t_vecRoh.maxCoeff(&t_iMaxIdx);
539 std::cout <<
"Iteration: " << r+1 <<
" of " << t_iMaxSearch
540 <<
"; Correlation: " << t_val_roh_k<<
"; Position (Idx+1): " << t_iIdx1+1 <<
" - " << t_iIdx2+1 <<
"\n\n";
546 MatrixX6T t_matProj_G_k_1(t_matOrthProj.rows(), t_matG_k_1.cols());
547 t_matProj_G_k_1 = t_matOrthProj * t_matG_k_1;
562 std::cout <<
"Searching stopped, last correlation " << t_val_roh_k;
563 std::cout <<
" is smaller then the given threshold " <<
m_dThreshold << std::endl << std::endl;
577 std::cout <<
"##### Calculation of RAP MUSIC completed ######"<< std::endl << std::endl << std::endl;
581 float t_fElapsedTime = ( (float)(end-start) / (float)CLOCKS_PER_SEC ) * 1000.0f;
582 std::cout <<
"Total Time Elapsed: " << t_fElapsedTime <<
" ms" << std::endl << std::endl;
587 return p_SourceEstimate;
597 if (p_matMeasurement.cols() > p_matMeasurement.rows())
600 t_matF =
MatrixXT(p_matMeasurement);
602 Eigen::JacobiSVD<MatrixXT> t_svdF(t_matF, Eigen::ComputeThinU);
604 int t_r =
getRank(t_svdF.singularValues().asDiagonal());
608 if (p_pMatPhi_s != NULL)
615 memcpy(p_pMatPhi_s->data(), t_svdF.matrixU().data(),
sizeof(double) *
m_iNumChannels * t_iCols);
628 Matrix6XT t_matU_A_T(6, p_matProj_G.rows());
630 Eigen::JacobiSVD<MatrixXT> t_svdProj_G(p_matProj_G, Eigen::ComputeThinU);
632 t_matSigma_A = t_svdProj_G.singularValues().asDiagonal();
633 t_matU_A_T = t_svdProj_G.matrixU().transpose();
642 MatrixXT t_matCor(t_matU_A_T_full.rows(), p_matU_B.cols());
645 t_matCor = t_matU_A_T_full*p_matU_B;
649 if (t_matCor.cols() > t_matCor.rows())
651 MatrixXT t_matCor_H(t_matCor.cols(), t_matCor.rows());
652 t_matCor_H = t_matCor.adjoint();
655 Eigen::JacobiSVD<MatrixXT> t_svdCor_H(t_matCor_H);
657 t_vecSigma_C = t_svdCor_H.singularValues();
661 Eigen::JacobiSVD<MatrixXT> t_svdCor(t_matCor);
663 t_vecSigma_C = t_svdCor.singularValues();
667 double t_dRetSigma_C;
668 t_dRetSigma_C = t_vecSigma_C(0);
673 return t_dRetSigma_C;
687 Eigen::JacobiSVD<MatrixXT> svdOfProj_G(p_matProj_G, Eigen::ComputeThinU | Eigen::ComputeThinV);
689 sigma_A = svdOfProj_G.singularValues().asDiagonal();
690 U_A_T = svdOfProj_G.matrixU().transpose();
691 V_A = svdOfProj_G.matrixV();
700 t_matCor = U_A_T*p_matU_B;
708 if (t_matCor.cols() > t_matCor.rows())
711 Cor_H = t_matCor.adjoint();
713 Eigen::JacobiSVD<MatrixXT> svdOfCor_H(Cor_H, Eigen::ComputeThinV);
715 U_C = svdOfCor_H.matrixV();
716 sigma_C = svdOfCor_H.singularValues();
720 Eigen::JacobiSVD<MatrixXT> svdOfCor(t_matCor, Eigen::ComputeThinU);
722 U_C = svdOfCor.matrixU();
723 sigma_C = svdOfCor.singularValues();
727 sigma_a_inv = sigma_A.inverse();
730 X = (V_A*sigma_a_inv)*U_C;
735 double norm_X = 1/(X_max.norm());
738 p_vec_phi_k_1 = X_max*norm_X;
745 ret_sigma_C = sigma_C(0);
762 VectorXT t_vec_a_theta_k_1(p_matG_k_1.rows(),1);
764 t_vec_a_theta_k_1 = p_matG_k_1*p_matPhi_k_1;
766 p_matA_k_1.block(0,p_iIdxk_1,p_matA_k_1.rows(),1) = t_vec_a_theta_k_1;
776 MatrixXT t_matA_k_1_tmp(p_matA_k_1.cols(), p_matA_k_1.cols());
777 t_matA_k_1_tmp = p_matA_k_1.adjoint()*p_matA_k_1;
780 int t_size = t_matA_k_1_tmp.cols();
782 while (!t_matA_k_1_tmp.block(0,0,t_size,t_size).fullPivLu().isInvertible())
787 MatrixXT t_matA_k_1_tmp_inv(t_matA_k_1_tmp.rows(), t_matA_k_1_tmp.cols());
788 t_matA_k_1_tmp_inv.setZero();
790 t_matA_k_1_tmp_inv.block(0,0,t_size,t_size) = t_matA_k_1_tmp.block(0,0,t_size,t_size).inverse();
793 t_matA_k_1_tmp = MatrixXT::Zero(p_matA_k_1.rows(), p_matA_k_1.cols());
795 t_matA_k_1_tmp = p_matA_k_1*t_matA_k_1_tmp_inv;
798 MatrixXT t_matA_k_1_tmp2(p_matA_k_1.rows(), p_matA_k_1.rows());
800 t_matA_k_1_tmp2 = t_matA_k_1_tmp*p_matA_k_1.adjoint();
806 p_matOrthProj = I-t_matA_k_1_tmp2;
816 const int p_iNumCombinations,
817 Pair** p_ppPairIdxCombinations)
const
824 #pragma omp parallel num_threads(m_iMaxNumThreads) private(idx1, idx2)
830 for (
int i = 0; i < p_iNumCombinations; ++i)
834 Pair* t_pairCombination =
new Pair();
835 t_pairCombination->
x1 = idx1;
836 t_pairCombination->
x2 = idx2;
838 p_ppPairIdxCombinations[i] = t_pairCombination;
848 int ii = p_iPoints*(p_iPoints+1)/2-1-p_iCurIdx;
849 int K = (int)floor((sqrt((
double)(8*ii+1))-1)/2);
851 p_iIdx1 = p_iPoints-1-K;
853 p_iIdx2 = (p_iCurIdx-p_iPoints*(p_iPoints+1)/2 + (K+1)*(K+2)/2)+p_iIdx1;
861 int p_iIdx1,
int p_iIdx2)
863 p_matGainMarix_Pair.block(0,0,p_matGainMarix.rows(),3) = p_matGainMarix.block(0, p_iIdx1*3, p_matGainMarix.rows(), 3);
865 p_matGainMarix_Pair.block(0,3,p_matGainMarix.rows(),3) = p_matGainMarix.block(0, p_iIdx2*3, p_matGainMarix.rows(), 3);
878 t_pRapDipolePair.
m_iIdx1 = p_iDipoleIdx1;
879 t_pRapDipolePair.
m_iIdx2 = p_iDipoleIdx2;
889 t_pRapDipolePair.
m_Dipole1.phi_x() = p_vec_phi_k_1[0];
890 t_pRapDipolePair.
m_Dipole1.phi_y() = p_vec_phi_k_1[1];
891 t_pRapDipolePair.
m_Dipole1.phi_z() = p_vec_phi_k_1[2];
893 t_pRapDipolePair.
m_Dipole2.phi_x() = p_vec_phi_k_1[3];
894 t_pRapDipolePair.
m_Dipole2.phi_y() = p_vec_phi_k_1[4];
895 t_pRapDipolePair.
m_Dipole2.phi_z() = p_vec_phi_k_1[5];
899 p_RapDipoles.append(t_pRapDipolePair);
Eigen::Matrix< double, 6, 1 > Vector6T
virtual const MNESourceSpace & getSourceSpace() const
virtual const char * getName() const
static void calcA_k_1(const MatrixX6T &p_matG_k_1, const Vector6T &p_matPhi_k_1, const int p_iIdxk_1, MatrixXT &p_matA_k_1)
static int useFullRank(const MatrixXT &p_Mat, const MatrixXT &p_matSigma_src, MatrixXT &p_matFull_Rank, int type=NOT_TRANSPOSED)
Source Space descritpion.
Eigen::Matrix< double, Eigen::Dynamic, 6 > MatrixX6T
Eigen::Matrix< double, 6, Eigen::Dynamic > Matrix6XT
int m_iNumLeadFieldCombinations
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixXT
void setStcAttr(int p_iSampStcWin, float p_fStcOverlap)
static void getGainMatrixPair(const MatrixXT &p_matGainMarix, MatrixX6T &p_matGainMarix_Pair, int p_iIdx1, int p_iIdx2)
void calcPairCombinations(const int p_iNumPoints, const int p_iNumCombinations, Pair **p_ppPairIdxCombinations) const
bool init(MNEForwardSolution &p_pFwd, bool p_bSparsed=false, int p_iN=2, double p_dThr=0.5)
Eigen::Matrix< double, 6, 6 > Matrix6T
virtual MNESourceEstimate calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal=false)
MNEForwardSolution m_ForwardSolution
void calcOrthProj(const MatrixXT &p_matA_k_1, MatrixXT &p_matOrthProj) const
FiffNamedMatrix::SDPtr sol
RapMusic algorithm class declaration.
static MatrixXT makeSquareMat(const MatrixXT &p_matF)
static int getRank(const MatrixXT &p_matSigma)
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorXT
static void getPointPair(const int p_iPoints, const int p_iCurIdx, int &p_iIdx1, int &p_iIdx2)
static void insertSource(int p_iDipoleIdx1, int p_iDipoleIdx2, const Vector6T &p_vec_phi_k_1, double p_valCor, QList< DipolePair< double > > &p_RapDipoles)
MNEMath class declaration.
int calcPhi_s(const MatrixXT &p_matMeasurement, MatrixXT *&p_pMatPhi_s) const
Pair ** m_ppPairIdxCombinations
static double subcorr(MatrixX6T &p_matProj_G, const MatrixXT &p_pMatU_B)
static int nchoose2(int n)