69 :
RapMusic(p_pFwd, p_bSparsed, p_iN, p_dThr)
72 init(p_pFwd, p_bSparsed, p_iN, p_dThr);
78 PwlRapMusic::~PwlRapMusic()
88 return "Powell RAP MUSIC";
117 std::cout <<
"RAP MUSIC wasn't initialized!";
118 return p_SourceEstimate;
124 std::cout <<
"Lead Field channels do not fit to number of measurement channels!";
125 return p_SourceEstimate;
141 int t_r =
calcPhi_s(p_matMeasurement, t_pMatPhi_s);
143 int t_iMaxSearch =
m_iN < t_r ?
m_iN : t_r;
147 std::cout <<
"Warning: Rank " << t_r <<
" of the measurement data is smaller than the " <<
m_iN;
148 std::cout <<
" sources to find." << std::endl;
149 std::cout <<
" Searching now for " << t_iMaxSearch <<
" correlated sources.";
150 std::cout << std::endl << std::endl;
156 t_matOrthProj.setIdentity();
160 t_matA_k_1.setZero();
176 p_RapDipoles.clear();
178 std::cout <<
"##### Calculation of PWL RAP MUSIC started ######\n\n";
180 MatrixXT t_matProj_Phi_s(t_matOrthProj.rows(), t_pMatPhi_s->cols());
184 for(
int r = 0; r < t_iMaxSearch ; ++r)
186 t_matProj_Phi_s = t_matOrthProj*(*t_pMatPhi_s);
193 Eigen::JacobiSVD< MatrixXT > t_svdProj_Phi_S(t_matProj_Phi_s, Eigen::ComputeThinU);
195 useFullRank(t_svdProj_Phi_S.matrixU(), t_svdProj_Phi_S.singularValues().asDiagonal(), t_matU_B);
203 clock_t start_subcorr, end_subcorr;
204 start_subcorr = clock();
209 int t_iCurrentRow = 2;
214 int t_iMaxIdx_old = -1;
224 while(t_iMaxFound == 0)
229 #pragma omp parallel num_threads(m_iMaxNumThreads)
235 for(
int i = 0; i < t_iNumVecElements; i++)
237 int k = t_pVecIdxElements(i);
240 MatrixX6T t_matProj_G(t_matProj_LeadField.rows(),6);
269 VectorXT::Index t_iMaxIdx;
271 t_val_roh_k = t_vecRoh.maxCoeff(&t_iMaxIdx);
273 if((
int)t_iMaxIdx == t_iMaxIdx_old)
280 t_iMaxIdx_old = t_iMaxIdx;
288 if(t_iIdx1 == t_iCurrentRow)
289 t_iCurrentRow = t_iIdx2;
291 t_iCurrentRow = t_iIdx1;
297 end_subcorr = clock();
299 float t_fSubcorrElapsedTime = ( (float)(end_subcorr-start_subcorr) / (float)CLOCKS_PER_SEC ) * 1000.0f;
300 std::cout <<
"Time Elapsed: " << t_fSubcorrElapsedTime <<
" ms" << std::endl;
304 std::cout <<
"Iteration: " << r+1 <<
" of " << t_iMaxSearch
305 <<
"; Correlation: " << t_val_roh_k<<
"; Position (Idx+1): " << t_iIdx1+1 <<
" - " << t_iIdx2+1 <<
"\n\n";
311 MatrixX6T t_matProj_G_k_1(t_matOrthProj.rows(), t_matG_k_1.cols());
312 t_matProj_G_k_1 = t_matOrthProj * t_matG_k_1;
327 std::cout <<
"Searching stopped, last correlation " << t_val_roh_k;
328 std::cout <<
" is smaller then the given threshold " <<
m_dThreshold << std::endl << std::endl;
342 std::cout <<
"##### Calculation of PWL RAP MUSIC completed ######"<< std::endl << std::endl << std::endl;
346 float t_fElapsedTime = ( (float)(end-start) / (float)CLOCKS_PER_SEC ) * 1000.0f;
347 std::cout <<
"Total Time Elapsed: " << t_fElapsedTime <<
" ms" << std::endl << std::endl;
352 return p_SourceEstimate;
358 int PwlRapMusic::PowellOffset(
int p_iRow,
int p_iNumPoints)
361 return p_iRow*p_iNumPoints - (( (p_iRow-1)*p_iRow) / 2);
368 void PwlRapMusic::PowellIdxVec(
int p_iRow,
int p_iNumPoints, Eigen::VectorXi& p_pVecElements)
378 for(
int i = 0; i <= p_iRow; ++i)
379 p_pVecElements(i) = PwlRapMusic::PowellOffset(i+1,p_iNumPoints)-(p_iNumPoints-p_iRow);
383 int off = PwlRapMusic::PowellOffset(p_iRow,p_iNumPoints);
384 int length = p_iNumPoints - p_iRow;
386 for(
int i = p_iRow; i < p_iRow+length; ++i)
388 p_pVecElements(i) = off+k;
Eigen::Matrix< double, 6, 1 > Vector6T
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)
Eigen::Matrix< double, Eigen::Dynamic, 6 > MatrixX6T
int m_iNumLeadFieldCombinations
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixXT
static void getGainMatrixPair(const MatrixXT &p_matGainMarix, MatrixX6T &p_matGainMarix_Pair, int p_iIdx1, int p_iIdx2)
bool init(MNEForwardSolution &p_pFwd, bool p_bSparsed=false, int p_iN=2, double p_dThr=0.5)
PwlRapMusic algorithm class declaration.
virtual MNESourceEstimate calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal=false)
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
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorXT
The RapMusic class provides the RAP MUSIC Algorithm CPU implementation. ToDo: Paper references...
static void insertSource(int p_iDipoleIdx1, int p_iDipoleIdx2, const Vector6T &p_vec_phi_k_1, double p_valCor, QList< DipolePair< double > > &p_RapDipoles)
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)
virtual const char * getName() const