61 #include <QtConcurrent>
72 using namespace FSLIB;
90 , source_rr(MatrixX3f::Zero(0,3))
91 , source_nn(MatrixX3f::Zero(0,3))
109 , source_rr(MatrixX3f::Zero(0,3))
110 , source_nn(MatrixX3f::Zero(0,3))
112 if(!
read(p_IODevice, *
this, force_fixed, surf_ori, include, exclude, bExcludeBads))
114 printf(
"\tForward solution not found.\n");
123 : info(p_MNEForwardSolution.info)
124 , source_ori(p_MNEForwardSolution.source_ori)
125 , surf_ori(p_MNEForwardSolution.surf_ori)
126 , coord_frame(p_MNEForwardSolution.coord_frame)
127 , nsource(p_MNEForwardSolution.nsource)
128 , nchan(p_MNEForwardSolution.nchan)
129 , sol(p_MNEForwardSolution.sol)
130 , sol_grad(p_MNEForwardSolution.sol_grad)
131 , mri_head_t(p_MNEForwardSolution.mri_head_t)
132 , src(p_MNEForwardSolution.src)
133 , source_rr(p_MNEForwardSolution.source_rr)
134 , source_nn(p_MNEForwardSolution.source_nn)
173 printf(
"Cluster forward solution using %s.\n", p_sMethod.toUtf8().constData());
182 printf(
"Error: Fixed orientation not implemented jet!\n");
195 MatrixXd t_G_Whitened(0,0);
196 bool t_bUseWhitened =
false;
204 MatrixXd p_outWhitener;
205 qint32 p_outNumNonZero;
207 this->
prepare_forward(p_pInfo, p_pNoise_cov,
false, p_outFwdInfo, t_G_Whitened, p_outNoiseCov, p_outWhitener, p_outNumNonZero);
208 printf(
"\tWhitening the forward solution.\n");
210 t_G_Whitened = p_outWhitener*t_G_Whitened;
211 t_bUseWhitened =
true;
223 for(qint32 h = 0; h < this->
src.
size(); ++h )
231 for(qint32 j = 0; j < h; ++j)
232 offset += this->
src[j].nuse;
235 printf(
"Cluster Left Hemisphere\n");
237 printf(
"Cluster Right Hemisphere\n");
239 Colortable t_CurrentColorTable = p_AnnotationSet[h].getColortable();
240 VectorXi label_ids = t_CurrentColorTable.
getLabelIds();
243 VectorXi vertno_labeled = VectorXi::Zero(this->
src[h].vertno.rows());
246 for(qint32 i = 0; i < vertno_labeled.rows(); ++i)
247 vertno_labeled[i] = p_AnnotationSet[h].getLabelIds()[this->
src[h].vertno[i]];
250 QList<RegionData> m_qListRegionDataIn;
255 for (qint32 i = 0; i < label_ids.rows(); ++i)
257 if (label_ids[i] != 0)
259 QString curr_name = t_CurrentColorTable.
struct_names[i];
260 printf(
"\tCluster %d / %li %s...", i+1, label_ids.rows(), curr_name.toUtf8().constData());
265 VectorXi idcs = VectorXi::Zero(vertno_labeled.rows());
269 for(qint32 j = 0; j < vertno_labeled.rows(); ++j)
271 if(vertno_labeled[j] == label_ids[i])
277 idcs.conservativeResize(c);
280 MatrixXd t_G(this->
sol->data.rows(), idcs.rows()*3);
281 MatrixXd t_G_Whitened_Roi(t_G_Whitened.rows(), idcs.rows()*3);
283 for(qint32 j = 0; j < idcs.rows(); ++j)
285 t_G.block(0, j*3, t_G.rows(), 3) = this->
sol->data.block(0, (idcs[j]+offset)*3, t_G.rows(), 3);
287 t_G_Whitened_Roi.block(0, j*3, t_G_Whitened_Roi.rows(), 3) = t_G_Whitened.block(0, (idcs[j]+offset)*3, t_G_Whitened_Roi.rows(), 3);
290 qint32 nSens = t_G.rows();
291 qint32 nSources = t_G.cols()/3;
299 t_sensG.
nClusters = ceil((
double)nSources/(
double)p_iClusterSize);
305 printf(
"%d Cluster(s)... ", t_sensG.
nClusters);
308 t_sensG.
matRoiG = MatrixXd(t_G.cols()/3, 3*nSens);
310 t_sensG.
matRoiGWhitened = MatrixXd(t_G_Whitened_Roi.cols()/3, 3*nSens);
312 for(qint32 j = 0; j < nSens; ++j)
314 for(qint32 k = 0; k < t_sensG.
matRoiG.rows(); ++k)
315 t_sensG.
matRoiG.block(k,j*3,1,3) = t_G.block(j,k*3,1,3);
318 t_sensG.
matRoiGWhitened.block(k,j*3,1,3) = t_G_Whitened_Roi.block(j,k*3,1,3);
325 m_qListRegionDataIn.append(t_sensG);
331 printf(
"failed! Label contains no sources.\n");
340 printf(
"Clustering... ");
341 QFuture< RegionDataOut > res;
342 res = QtConcurrent::mapped(m_qListRegionDataIn, &RegionData::cluster);
343 res.waitForFinished();
348 MatrixXd t_G_partial;
352 QList<RegionData>::const_iterator itIn;
353 itIn = m_qListRegionDataIn.begin();
354 QFuture<RegionDataOut>::const_iterator itOut;
355 for (itOut = res.constBegin(); itOut != res.constEnd(); ++itOut)
357 nClusters = itOut->ctrs.rows();
358 nSens = itOut->ctrs.cols()/3;
359 t_G_partial = MatrixXd::Zero(nSens, nClusters*3);
367 for(qint32 j = 0; j < nSens; ++j)
368 for(qint32 k = 0; k < nClusters; ++k)
369 t_G_partial.block(j, k*3, 1, 3) = itOut->ctrs.block(k,j*3,1,3);
374 for(qint32 j = 0; j < nClusters; ++j)
376 VectorXi clusterIdcs = VectorXi::Zero(itOut->roiIdx.rows());
377 VectorXd clusterDistance = VectorXd::Zero(itOut->roiIdx.rows());
378 MatrixX3f clusterSource_rr = MatrixX3f::Zero(itOut->roiIdx.rows(), 3);
379 qint32 nClusterIdcs = 0;
380 for(qint32 k = 0; k < itOut->roiIdx.rows(); ++k)
382 if(itOut->roiIdx[k] == j)
384 clusterIdcs[nClusterIdcs] = itIn->idcs[k];
386 qint32 offset = h == 0 ? 0 : this->
src[0].nuse;
387 clusterSource_rr.row(nClusterIdcs) = this->
source_rr.row(offset + itIn->idcs[k]);
388 clusterDistance[nClusterIdcs] = itOut->D(k,j);
392 clusterIdcs.conservativeResize(nClusterIdcs);
393 clusterSource_rr.conservativeResize(nClusterIdcs,3);
394 clusterDistance.conservativeResize(nClusterIdcs);
396 VectorXi clusterVertnos = VectorXi::Zero(clusterIdcs.size());
397 for(qint32 k = 0; k < clusterVertnos.size(); ++k)
398 clusterVertnos(k) = this->
src[h].vertno[clusterIdcs(k)];
401 p_fwdOut.
src[h].cluster_info.clusterVertnos.append(clusterVertnos);
402 p_fwdOut.
src[h].cluster_info.clusterSource_rr.append(clusterSource_rr);
403 p_fwdOut.
src[h].cluster_info.clusterDistances.append(clusterDistance);
404 p_fwdOut.
src[h].cluster_info.clusterLabelIds.append(label_ids[itOut->iLabelIdxOut]);
405 p_fwdOut.
src[h].cluster_info.clusterLabelNames.append(t_CurrentColorTable.
getNames()[itOut->iLabelIdxOut]);
412 if(t_G_partial.rows() > 0 && t_G_partial.cols() > 0)
414 t_G_new.conservativeResize(t_G_partial.rows(), t_G_new.cols() + t_G_partial.cols());
415 t_G_new.block(0, t_G_new.cols() - t_G_partial.cols(), t_G_new.rows(), t_G_partial.cols()) = t_G_partial;
418 for(qint32 k = 0; k < nClusters; ++k)
422 double sqec = sqrt((itIn->matRoiGOrig.block(0, j*3, itIn->matRoiGOrig.rows(), 3) - t_G_partial.block(0, k*3, t_G_partial.rows(), 3)).array().pow(2).sum());
423 double sqec_min = sqec;
426 for(qint32 j = 1; j < itIn->idcs.rows(); ++j)
428 sqec = sqrt((itIn->matRoiGOrig.block(0, j*3, itIn->matRoiGOrig.rows(), 3) - t_G_partial.block(0, k*3, t_G_partial.rows(), 3)).array().pow(2).sum());
441 qint32 sel_idx = itIn->idcs[j_min];
443 p_fwdOut.
src[h].cluster_info.centroidVertno.append(this->
src[h].vertno[sel_idx]);
444 p_fwdOut.
src[h].cluster_info.centroidSource_rr.append(this->
src[h].rr.row(sel_idx));
450 p_fwdOut.
src[h].vertno[count] = p_fwdOut.
src[h].cluster_info.clusterLabelIds[count];
469 p_fwdOut.
src[h].vertno.conservativeResize(count);
478 qint32 totalNumOfClust = 0;
479 for (qint32 h = 0; h < 2; ++h)
480 totalNumOfClust += p_fwdOut.
src[h].cluster_info.clusterVertnos.
size();
483 p_D = MatrixXd::Zero(this->
sol->data.cols(), totalNumOfClust);
485 p_D = MatrixXd::Zero(this->
sol->data.cols(), totalNumOfClust*3);
492 qint32 currentCluster = 0;
493 for (qint32 h = 0; h < 2; ++h)
495 int hemiOffset = h == 0 ? 0 : t_vertnos[0].size();
496 for(qint32 i = 0; i < p_fwdOut.
src[h].cluster_info.clusterVertnos.
size(); ++i)
505 idx_sel.array() += hemiOffset;
511 double selectWeight = 1.0/idx_sel.size();
514 for(qint32 j = 0; j < idx_sel.size(); ++j)
515 p_D.col(currentCluster)[idx_sel(j)] = selectWeight;
519 qint32 clustOffset = currentCluster*3;
520 for(qint32 j = 0; j < idx_sel.size(); ++j)
522 qint32 idx_sel_Offset = idx_sel(j)*3;
524 p_D(idx_sel_Offset,clustOffset) = selectWeight;
526 p_D(idx_sel_Offset+1, clustOffset+1) = selectWeight;
528 p_D(idx_sel_Offset+2, clustOffset+2) = selectWeight;
692 p_fwdOut.
sol->data = t_G_new;
693 p_fwdOut.
sol->ncol = t_G_new.cols();
708 qint32 np = isFixed ? p_fwdOut.
sol->data.cols() : p_fwdOut.
sol->data.cols()/3;
710 if(p_iNumDipoles > np)
713 VectorXi sel(p_iNumDipoles);
715 float t_fStep = (float)np/(
float)p_iNumDipoles;
717 for(qint32 i = 0; i < p_iNumDipoles; ++i)
719 float t_fCurrent = ((float)i)*t_fStep;
720 sel[i] = (quint32)floor(t_fCurrent);
725 p_D = MatrixXd::Zero(p_fwdOut.
sol->data.cols(), p_iNumDipoles);
726 for(qint32 i = 0; i < p_iNumDipoles; ++i)
731 p_D = MatrixXd::Zero(p_fwdOut.
sol->data.cols(), p_iNumDipoles*3);
732 for(qint32 i = 0; i < p_iNumDipoles; ++i)
733 for(qint32 j = 0; j < 3; ++j)
734 p_D((sel[i]*3)+j, (i*3)+j) = 1;
785 p_fwdOut.
sol->data = this->
sol->data * p_D;
787 MatrixX3f rr(p_iNumDipoles,3);
789 MatrixX3f nn(p_iNumDipoles,3);
792 for(qint32 i = 0; i < p_iNumDipoles; ++i)
801 p_fwdOut.
sol->ncol = p_fwdOut.
sol->data.cols();
803 p_fwdOut.
nsource = p_iNumDipoles;
813 printf(
"\tCreating the depth weighting matrix...\n");
824 d = (G.array().square()).rowwise().sum();
829 qint32 n_pos = G.cols() / 3;
830 d = VectorXd::Zero(n_pos);
832 for (qint32 k = 0; k < n_pos; ++k)
834 Gk = G.block(0,3*k, G.rows(), 3);
835 JacobiSVD<MatrixXd> svd(Gk.transpose()*Gk);
836 d[k] = svd.singularValues().maxCoeff();
841 if(patch_areas.cols() > 0)
844 printf(
"\tToDo!!!!! >>> Patch areas taken into account in the depth weighting\n");
848 VectorXd w = d.cwiseInverse();
851 MNEMath::sort<double>(ws,
false);
852 double weight_limit = pow(limit, 2);
853 if (!limit_depth_chs)
859 limit = ws[ind] * weight_limit;
864 limit = ws[ws.size()-1];
867 if (ws[ws.size()-1] > weight_limit * ws[0])
869 double th = weight_limit * ws[0];
870 for(qint32 i = 0; i < ws.size(); ++i)
883 printf(
"\tlimit = %d/%li = %f", n_limit + 1, d.size(), sqrt(limit / ws[0]));
884 double scale = 1.0 / limit;
885 printf(
"\tscale = %g exp = %g", scale, exp);
887 VectorXd t_w = w.array() / limit;
888 for(qint32 i = 0; i < t_w.size(); ++i)
889 t_w[i] = t_w[i] > 1 ? 1 : t_w[i];
890 wpp = t_w.array().pow(exp);
894 depth_prior.
data = wpp;
897 depth_prior.
data.resize(wpp.rows()*3, 1);
900 for(qint32 i = 0; i < wpp.rows(); ++i)
904 depth_prior.
data(idx, 0) = v;
905 depth_prior.
data(idx+1, 0) = v;
906 depth_prior.
data(idx+2, 0) = v;
911 depth_prior.
diag =
true;
912 depth_prior.
dim = depth_prior.
data.rows();
913 depth_prior.
nfree = 1;
924 qint32 n_sources = this->
sol->data.cols();
926 if (0 <= loose && loose <= 1)
930 printf(
"\tForward operator is not oriented in surface coordinates. loose parameter should be None not %f.", loose);
932 printf(
"\tSetting loose to %f.\n", loose);
937 printf(
"\tIgnoring loose parameter with forward operator with fixed orientation.\n");
943 if(loose < 0 || loose > 1)
945 qWarning(
"Warning: Loose value should be in interval [0,1] not %f.\n", loose);
946 loose = loose > 1 ? 1 : 0;
947 printf(
"Setting loose to %f.\n", loose);
952 orient_prior.
data = VectorXd::Ones(n_sources);
953 if(!is_fixed_ori && (0 <= loose && loose <= 1))
955 printf(
"\tApplying loose dipole orientations. Loose value of %f.\n", loose);
956 for(qint32 i = 0; i < n_sources; i+=3)
957 orient_prior.
data.block(i,0,2,1).array() *= loose;
960 orient_prior.
diag =
true;
961 orient_prior.
dim = orient_prior.
data.size();
962 orient_prior.
nfree = 1;
974 if(include.size() == 0 && exclude.size() == 0)
980 quint32 nuse = sel.size();
984 printf(
"Nothing remains after picking. Returning original forward solution.\n");
987 printf(
"\t%d out of %d channels remain after picking\n", nuse, fwd.
nchan);
990 MatrixXd newData(nuse, fwd.
sol->data.cols());
991 for(quint32 i = 0; i < nuse; ++i)
992 newData.row(i) = fwd.
sol->data.row(sel[i]);
994 fwd.
sol->data = newData;
995 fwd.
sol->nrow = nuse;
997 QStringList ch_names;
998 for(qint32 i = 0; i < sel.cols(); ++i)
999 ch_names << fwd.
sol->row_names[sel(i)];
1001 fwd.
sol->row_names = ch_names;
1003 QList<FiffChInfo> chs;
1004 for(qint32 i = 0; i < sel.cols(); ++i)
1005 chs.append(fwd.
info.
chs[sel(i)]);
1010 for(qint32 i = 0; i < fwd.
info.
bads.size(); ++i)
1011 if(ch_names.contains(fwd.
info.
bads[i]))
1017 newData.resize(nuse, fwd.
sol_grad->data.cols());
1018 for(quint32 i = 0; i < nuse; ++i)
1019 newData.row(i) = fwd.
sol_grad->data.row(sel[i]);
1022 QStringList row_names;
1023 for(qint32 i = 0; i < sel.cols(); ++i)
1024 row_names << fwd.
sol_grad->row_names[sel(i)];
1025 fwd.
sol_grad->row_names = row_names;
1036 VectorXi selVertices;
1039 for(qint32 i = 0; i < p_qListLabels.size(); ++i)
1041 VectorXi currentSelection;
1044 selVertices.conservativeResize(iSize+currentSelection.size());
1045 selVertices.block(iSize,0,currentSelection.size(),1) = currentSelection;
1046 iSize = selVertices.size();
1053 MatrixX3f rr(selVertices.size(),3);
1054 MatrixX3f nn(selVertices.size(),3);
1056 for(qint32 i = 0; i < selVertices.size(); ++i)
1058 rr.block(i, 0, 1, 3) = selectedFwd.
source_rr.row(selVertices[i]);
1059 nn.block(i, 0, 1, 3) = selectedFwd.
source_nn.row(selVertices[i]);
1066 MatrixXd G(selectedFwd.
sol->data.rows(),selSolIdcs.size());
1068 qint32 rows = G.rows();
1070 for(qint32 i = 0; i < selSolIdcs.size(); ++i)
1071 G.block(0, i, rows, 1) = selectedFwd.
sol->data.col(selSolIdcs[i]);
1073 selectedFwd.
sol->data = G;
1074 selectedFwd.
sol->nrow = selectedFwd.
sol->data.rows();
1075 selectedFwd.
sol->ncol = selectedFwd.
sol->data.cols();
1076 selectedFwd.
nsource = selectedFwd.
sol->ncol / 3;
1088 RowVectorXi sel =
info.
pick_types(meg, eeg,
false, include, exclude);
1090 QStringList include_ch_names;
1091 for(qint32 i = 0; i < sel.cols(); ++i)
1094 return this->pick_channels(include_ch_names);
1102 QStringList fwd_ch_names, ch_names;
1103 for(qint32 i = 0; i < this->
info.
chs.size(); ++i)
1104 fwd_ch_names << this->
info.
chs[i].ch_name;
1107 for(qint32 i = 0; i < p_info.
chs.size(); ++i)
1108 if( !p_info.
bads.contains(p_info.
chs[i].ch_name)
1109 && !p_noise_cov.
bads.contains(p_info.
chs[i].ch_name)
1110 && fwd_ch_names.contains(p_info.
chs[i].ch_name))
1111 ch_names << p_info.
chs[i].ch_name;
1113 qint32 n_chan = ch_names.size();
1114 printf(
"Computing inverse operator with %d channels.\n", n_chan);
1122 p_outNumNonZero = 0;
1123 VectorXi t_vecNonZero = VectorXi::Zero(n_chan);
1124 for(qint32 i = 0; i < p_outNoiseCov.
eig.rows(); ++i)
1126 if(p_outNoiseCov.
eig[i] > 0)
1128 t_vecNonZero[p_outNumNonZero] = i;
1132 if(p_outNumNonZero > 0)
1133 t_vecNonZero.conservativeResize(p_outNumNonZero);
1135 if(p_outNumNonZero > 0)
1139 qWarning(
"Warning in MNEForwardSolution::prepare_forward: if (p_pca) havent been debugged.");
1140 p_outWhitener = MatrixXd::Zero(n_chan, p_outNumNonZero);
1142 for(qint32 i = 0; i < p_outNumNonZero; ++i)
1143 p_outWhitener.col(t_vecNonZero[i]) = p_outNoiseCov.
eigvec.col(t_vecNonZero[i]).array() / sqrt(p_outNoiseCov.
eig(t_vecNonZero[i]));
1144 printf(
"\tReducing data rank to %d.\n", p_outNumNonZero);
1148 printf(
"Creating non pca whitener.\n");
1149 p_outWhitener = MatrixXd::Zero(n_chan, n_chan);
1150 for(qint32 i = 0; i < p_outNumNonZero; ++i)
1151 p_outWhitener(t_vecNonZero[i],t_vecNonZero[i]) = 1.0 / sqrt(p_outNoiseCov.
eig(t_vecNonZero[i]));
1153 p_outWhitener *= p_outNoiseCov.
eigvec;
1157 VectorXi fwd_idx = VectorXi::Zero(ch_names.size());
1158 VectorXi info_idx = VectorXi::Zero(ch_names.size());
1160 qint32 count_fwd_idx = 0;
1161 qint32 count_info_idx = 0;
1162 for(qint32 i = 0; i < ch_names.size(); ++i)
1164 idx = fwd_ch_names.indexOf(ch_names[i]);
1167 fwd_idx[count_fwd_idx] = idx;
1170 idx = p_info.
ch_names.indexOf(ch_names[i]);
1173 info_idx[count_info_idx] = idx;
1177 fwd_idx.conservativeResize(count_fwd_idx);
1178 info_idx.conservativeResize(count_info_idx);
1180 gain.resize(count_fwd_idx, this->
sol->data.cols());
1181 for(qint32 i = 0; i < count_fwd_idx; ++i)
1182 gain.row(i) = this->
sol->data.row(fwd_idx[i]);
1184 p_outFwdInfo = p_info.
pick_info(info_idx);
1186 printf(
"\tTotal rank is %d\n", p_outNumNonZero);
1196 QList<FiffDirEntry> t_Dir;
1198 printf(
"Reading forward solution from %s...\n", t_pStream->streamName().toUtf8().constData());
1199 if(!t_pStream->open(t_Tree, t_Dir))
1204 QList<FiffDirTree> fwds = t_Tree.
dir_tree_find(FIFFB_MNE_FORWARD_SOLUTION);
1206 if (fwds.size() == 0)
1208 t_pStream->device()->close();
1209 std::cout <<
"No forward solutions in " << t_pStream->streamName().toUtf8().constData();
1215 QList<FiffDirTree> parent_mri = t_Tree.
dir_tree_find(FIFFB_MNE_PARENT_MRI_FILE);
1216 if (parent_mri.size() == 0)
1218 t_pStream->device()->close();
1219 std::cout <<
"No parent MRI information in " << t_pStream->streamName().toUtf8().constData();
1226 t_pStream->device()->close();
1227 std::cout <<
"Could not read the source spaces\n";
1232 for(qint32 k = 0; k < t_SourceSpace.
size(); ++k)
1241 bads = t_pStream->read_bad_channels(t_Tree);
1244 printf(
"\t%d bad channels ( ",bads.size());
1245 for(qint32 i = 0; i < bads.size(); ++i)
1246 printf(
"\"%s\" ", bads[i].toLatin1().constData());
1257 for(qint32 k = 0; k < fwds.size(); ++k)
1259 if(!fwds[k].find_tag(t_pStream.data(), FIFF_MNE_INCLUDED_METHODS, t_pTag))
1261 t_pStream->device()->close();
1262 std::cout <<
"Methods not listed for one of the forward solutions\n";
1265 if (*t_pTag->toInt() == FIFFV_MNE_MEG)
1267 printf(
"MEG solution found\n");
1270 else if(*t_pTag->toInt() == FIFFV_MNE_EEG)
1272 printf(
"EEG solution found\n");
1273 eegnode = fwds.at(k);
1279 if (read_one(t_pStream.data(), megnode, megfwd))
1281 if (megfwd.
source_ori == FIFFV_MNE_FIXED_ORI)
1282 ori = QString(
"fixed");
1284 ori = QString(
"free");
1285 printf(
"\tRead MEG forward solution (%d sources, %d channels, %s orientations)\n", megfwd.
nsource,megfwd.
nchan,ori.toUtf8().constData());
1288 if (read_one(t_pStream.data(), eegnode, eegfwd))
1290 if (eegfwd.
source_ori == FIFFV_MNE_FIXED_ORI)
1291 ori = QString(
"fixed");
1293 ori = QString(
"free");
1294 printf(
"\tRead EEG forward solution (%d sources, %d channels, %s orientations)\n", eegfwd.
nsource,eegfwd.
nchan,ori.toUtf8().constData());
1304 if (megfwd.
sol->data.cols() != eegfwd.
sol->data.cols() ||
1309 t_pStream->device()->close();
1310 std::cout <<
"The MEG and EEG forward solutions do not match\n";
1315 fwd.
sol->data = MatrixXd(megfwd.
sol->nrow + eegfwd.
sol->nrow, megfwd.
sol->ncol);
1317 fwd.
sol->data.block(0,0,megfwd.
sol->nrow,megfwd.
sol->ncol) = megfwd.
sol->data;
1318 fwd.
sol->data.block(megfwd.
sol->nrow,0,eegfwd.
sol->nrow,eegfwd.
sol->ncol) = eegfwd.
sol->data;
1319 fwd.
sol->nrow = megfwd.
sol->nrow + eegfwd.
sol->nrow;
1320 fwd.
sol->row_names.append(eegfwd.
sol->row_names);
1333 printf(
"\tMEG and EEG forward solutions combined\n");
1343 if(!parent_mri[0].find_tag(t_pStream.data(), FIFF_COORD_TRANS, t_pTag))
1345 t_pStream->device()->close();
1346 std::cout <<
"MRI/head coordinate transformation not found\n";
1358 t_pStream->device()->close();
1359 std::cout <<
"MRI/head coordinate transformation not found\n";
1368 t_pStream->read_meas_info_base(t_Tree, fwd.
info);
1371 t_pStream->device()->close();
1379 std::cout <<
"Only forward solutions computed in MRI or head coordinates are acceptable";
1386 for(qint32 k = 0; k < t_SourceSpace.
size(); ++k)
1387 nuse += t_SourceSpace[k].nuse;
1390 throw(
"Source spaces do not match the forward solution.\n");
1392 printf(
"\tSource spaces transformed to the forward solution coordinate frame\n");
1393 fwd.
src = t_SourceSpace;
1402 for(qint32 k = 0; k < t_SourceSpace.
size();++k)
1404 for(qint32 q = 0; q < t_SourceSpace[k].nuse; ++q)
1406 fwd.
source_rr.block(q,0,1,3) = t_SourceSpace[k].rr.block(t_SourceSpace[k].vertno(q),0,1,3);
1407 fwd.
source_nn.block(q,0,1,3) = t_SourceSpace[k].nn.block(t_SourceSpace[k].vertno(q),0,1,3);
1409 nuse += t_SourceSpace[k].nuse;
1416 printf(
"\tChanging to fixed-orientation forward solution...");
1418 MatrixXd tmp = fwd.
source_nn.transpose().cast<
double>();
1420 fwd.
sol->data *= (*fix_rot);
1426 SparseMatrix<double> t_matKron;
1427 SparseMatrix<double> t_eye(3,3);
1428 for (qint32 i = 0; i < 3; ++i)
1429 t_eye.insert(i,i) = 1.0f;
1430 t_matKron = kroneckerProduct(*fix_rot,t_eye);
1443 printf(
"\tConverting to surface-based source orientations...");
1445 bool use_ave_nn =
false;
1446 if(t_SourceSpace[0].patch_inds.
size() > 0)
1449 printf(
"\tAverage patch normals will be employed in the rotation to the local surface coordinates...\n");
1457 qWarning(
"Warning source_ori: Rotating the source coordinate system haven't been verified --> Singular Vectors U are different from MATLAB!");
1459 for(qint32 k = 0; k < t_SourceSpace.
size();++k)
1462 for (qint32 q = 0; q < t_SourceSpace[k].nuse; ++q)
1463 fwd.
source_rr.block(q+nuse,0,1,3) = t_SourceSpace[k].rr.block(t_SourceSpace[k].vertno(q),0,1,3);
1465 for (qint32 p = 0; p < t_SourceSpace[k].nuse; ++p)
1473 VectorXi t_vIdx = t_SourceSpace[k].pinfo[t_SourceSpace[k].patch_inds[p]];
1474 Matrix3Xf t_nn(3, t_vIdx.size());
1475 for(qint32 i = 0; i < t_vIdx.size(); ++i)
1476 t_nn.col(i) = t_SourceSpace[k].nn.block(t_vIdx[i],0,1,3).transpose();
1477 nn = t_nn.rowwise().sum();
1478 nn.array() /= nn.norm();
1481 nn = t_SourceSpace[k].nn.block(t_SourceSpace[k].vertno(p),0,1,3).transpose();
1483 Matrix3f tmp = Matrix3f::Identity(nn.rows(), nn.rows()) - nn*nn.transpose();
1485 JacobiSVD<MatrixXf> t_svd(tmp, Eigen::ComputeThinU);
1487 VectorXf t_s = t_svd.singularValues();
1488 MatrixXf U = t_svd.matrixU();
1489 MNEMath::sort<float>(t_s, U);
1494 if ((nn.transpose() * U.block(0,2,3,1))(0,0) < 0)
1496 fwd.
source_nn.block(pp, 0, 3, 3) = U.transpose();
1499 nuse += t_SourceSpace[k].nuse;
1501 MatrixXd tmp = fwd.
source_nn.transpose().cast<
double>();
1504 fwd.
sol->data *= *surf_rot;
1508 SparseMatrix<double> t_matKron;
1509 SparseMatrix<double> t_eye(3,3);
1510 for (qint32 i = 0; i < 3; ++i)
1511 t_eye.insert(i,i) = 1.0f;
1512 t_matKron = kroneckerProduct(*surf_rot,t_eye);
1520 printf(
"\tCartesian source orientations...");
1523 for(qint32 k = 0; k < t_SourceSpace.
size(); ++k)
1525 for (qint32 q = 0; q < t_SourceSpace[k].nuse; ++q)
1526 fwd.
source_rr.block(q+nuse,0,1,3) = t_SourceSpace[k].rr.block(t_SourceSpace[k].vertno(q),0,1,3);
1528 nuse += t_SourceSpace[k].nuse;
1531 MatrixXf t_ones = MatrixXf::Ones(fwd.
nsource,1);
1532 Matrix3f t_eye = Matrix3f::Identity();
1533 fwd.
source_nn = kroneckerProduct(t_ones,t_eye);
1541 QStringList exclude_bads = exclude;
1542 if (bads.size() > 0)
1544 for(qint32 k = 0; k < bads.size(); ++k)
1545 if(!exclude_bads.contains(bads[k],Qt::CaseInsensitive))
1546 exclude_bads << bads[k];
1674 t_pStream->device()->close();
1695 p_pStream->device()->close();
1696 std::cout <<
"Source orientation tag not found.";
1704 p_pStream->device()->close();
1705 std::cout <<
"Coordinate frame tag not found.";
1713 p_pStream->device()->close();
1714 std::cout <<
"Number of sources not found.";
1718 one.
nsource = *t_pTag->toInt();
1720 if(!p_Node.
find_tag(p_pStream, FIFF_NCHAN, t_pTag))
1722 p_pStream->device()->close();
1723 printf(
"Number of channels not found.");
1727 one.
nchan = *t_pTag->toInt();
1730 one.
sol->transpose_named_matrix();
1733 p_pStream->device()->close();
1734 printf(
"Forward solution data not found .");
1740 one.
sol_grad->transpose_named_matrix();
1745 if (one.
sol->data.rows() != one.
nchan ||
1748 p_pStream->device()->close();
1749 printf(
"Forward solution matrix has wrong dimensions.\n");
1758 p_pStream->device()->close();
1759 printf(
"Forward solution gradient matrix has wrong dimensions.\n");
1772 if(info.
chs.size() != G.rows())
1774 printf(
"Error G.rows() and length of info.chs do not match: %li != %i", G.rows(), info.
chs.size());
1778 RowVectorXi sel = info.
pick_types(QString(
"grad"));
1781 for(qint32 i = 0; i < sel.size(); ++i)
1782 G.row(i) = G.row(sel[i]);
1783 G.conservativeResize(sel.size(), G.cols());
1784 printf(
"\t%li planar channels", sel.size());
1791 for(qint32 i = 0; i < sel.size(); ++i)
1792 G.row(i) = G.row(sel[i]);
1793 G.conservativeResize(sel.size(), G.cols());
1794 printf(
"\t%li magnetometer or axial gradiometer channels", sel.size());
1801 for(qint32 i = 0; i < sel.size(); ++i)
1802 G.row(i) = G.row(sel[i]);
1803 G.conservativeResize(sel.size(), G.cols());
1804 printf(
"\t%li EEG channels\n", sel.size());
1807 printf(
"Could not find MEG or EEG channels\n");
1819 qWarning(
"Warning: Only surface-oriented, free-orientation forward solutions can be converted to fixed orientaton.\n");
1823 for(qint32 i = 2; i < this->
sol->data.cols(); i += 3)
1824 this->
sol->data.col(count) = this->
sol->data.col(i);
1825 this->
sol->data.conservativeResize(this->
sol->data.rows(), count);
1826 this->
sol->ncol = this->
sol->ncol / 3;
1828 printf(
"\tConverted the forward solution into the fixed-orientation mode.\n");
static qint32 find_source_space_hemi(MNEHemisphere &p_Hemisphere)
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
MNEForwardSolution pick_types(bool meg, bool eeg, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
#define FIFF_MNE_SOURCE_ORIENTATION
FIFF measurement file information.
Source Space descritpion.
bool find_tag(FiffStream *p_pStream, fiff_int_t findkind, QSharedPointer< FiffTag > &p_pTag) const
bool transform_source_space_to(fiff_int_t dest, FiffCoordTrans &trans)
MNEForwardSolution cluster_forward_solution(const AnnotationSet &p_AnnotationSet, qint32 p_iClusterSize, MatrixXd &p_D=defaultD, const FiffCov &p_pNoise_cov=defaultCov, const FiffInfo &p_pInfo=defaultInfo, QString p_sMethod="cityblock") const
QList< FiffDirTree > dir_tree_find(fiff_int_t p_kind) const
QSharedDataPointer< FiffNamedMatrix > SDPtr
VectorXi getLabelIds() const
FiffInfo pick_info(const RowVectorXi &sel=defaultVectorXi) const
FiffCoordTrans mri_head_t
QSharedPointer< FiffTag > SPtr
Colortable class declaration.
static void restrict_gain_matrix(MatrixXd &G, const FiffInfo &info)
VectorXi tripletSelection(const VectorXi &p_vecIdxSelection) const
static bool readFromStream(FiffStream::SPtr &p_pStream, bool add_geom, FiffDirTree &p_Tree, MNESourceSpace &p_SourceSpace)
static SparseMatrix< double > * make_block_diag(const MatrixXd &A, qint32 n)
QSharedPointer< FiffStream > SPtr
bool read_named_matrix(const FiffDirTree &p_Node, fiff_int_t matkind, FiffNamedMatrix &mat)
static FiffCov compute_depth_prior(const MatrixXd &Gain, const FiffInfo &gain_info, bool is_fixed_ori, double exp=0.8, double limit=10.0, const MatrixXd &patch_areas=defaultConstMatrixXd, bool limit_depth_chs=false)
FiffNamedMatrix::SDPtr sol_grad
MNEForwardSolution class declaration, which provides the forward solution including the source space ...
static VectorXi intersect(const VectorXi &v1, const VectorXi &v2, VectorXi &idx_sel)
Directory tree structure.
static VectorXi sort(Matrix< T, Dynamic, 1 > &v, bool desc=true)
#define FIFFV_MNE_ORIENT_PRIOR_COV
FiffNamedMatrix::SDPtr sol
#define FIFF_MNE_COORD_FRAME
MNEForwardSolution pick_regions(const QList< Label > &p_qListLabels) const
MNEForwardSolution pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
RowVectorXi pick_types(const QString meg, bool eeg=false, bool stim=false, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
Vertices label based lookup table.
FiffCov prepare_noise_cov(const FiffInfo &p_info, const QStringList &p_chNames) const
static bool read(QIODevice &p_IODevice, MNEForwardSolution &fwd, bool force_fixed=false, bool surf_ori=false, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList, bool bExcludeBads=true)
KMeans class declaration.
QList< VectorXi > get_vertno() const
MNESourceSpace pick_regions(const QList< Label > &p_qListLabels) const
bool isFixedOrient() const
void prepare_forward(const FiffInfo &p_info, const FiffCov &p_noise_cov, bool p_pca, FiffInfo &p_outFwdInfo, MatrixXd &gain, FiffCov &p_outNoiseCov, MatrixXd &p_outWhitener, qint32 &p_outNumNonZero) const
FiffCov compute_orient_prior(float loose=0.2)
MNEMath class declaration.
MNEForwardSolution reduce_forward_solution(qint32 p_iNumDipoles, MatrixXd &p_D) const
static RowVectorXi pick_channels(const QStringList &ch_names, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList)
#define FIFFV_MNE_DEPTH_PRIOR_COV
QStringList getNames() const
QList< VectorXi > label_src_vertno_sel(const Label &p_label, VectorXi &src_sel) const