53 using namespace FSLIB;
62 MNESourceSpace::MNESourceSpace()
70 : m_qListHemispheres(p_MNESourceSpace.m_qListHemispheres)
88 m_qListHemispheres.clear();
96 QList<VectorXi> p_vertices;
97 for(qint32 i = 0; i < m_qListHemispheres.size(); ++i)
98 p_vertices.push_back(m_qListHemispheres[i].vertno);
110 QList<VectorXi> vertno;
111 vertno << this->m_qListHemispheres[0].vertno << this->m_qListHemispheres[1].vertno;
113 if (p_label.
hemi == 0)
116 vertno[0] = vertno_sel;
117 vertno[1] = VectorXi();
119 else if (p_label.
hemi == 1)
122 src_sel.array() += p_label.
vertices.size();
123 vertno[0] = VectorXi();
124 vertno[1] = vertno_sel;
153 qWarning(
"Unknown hemisphere type\n");
154 vertno[0] = VectorXi::Zero(0);
155 vertno[1] = VectorXi::Zero(0);
166 Q_UNUSED(p_qListLabels);
170 for(qint32 h = 0; h < 2; ++h)
172 VectorXi selVertices;
176 for(qint32 i = 0; i < p_qListLabels.size(); ++i)
178 if(p_qListLabels[i].hemi == h)
180 VectorXi currentSelection;
182 MNEMath::intersect(m_qListHemispheres[h].vertno, p_qListLabels[i].vertices, currentSelection);
184 selVertices.conservativeResize(iSize+currentSelection.size());
185 selVertices.block(iSize,0,currentSelection.size(),1) = currentSelection;
186 iSize = selVertices.size();
192 VectorXi newVertno(selVertices.size());
194 selectedSrc.m_qListHemispheres[h].inuse = VectorXi::Zero(selectedSrc.m_qListHemispheres[h].np);
196 for(qint32 i = 0; i < selVertices.size(); ++i)
198 selectedSrc.m_qListHemispheres[h].inuse[selVertices[i]] = 1;
199 newVertno[i] = this->m_qListHemispheres[h].vertno[selVertices[i]];
202 selectedSrc.m_qListHemispheres[h].nuse = selVertices.size();
203 selectedSrc.m_qListHemispheres[h].vertno = newVertno;
208 VectorXi idx_select = VectorXi::Zero(this->m_qListHemispheres[h].use_tris.rows());
209 for(qint32 i = 0; i < 3; ++i)
211 VectorXi tri_dim = this->m_qListHemispheres[h].use_tris.col(i);
215 for(qint32 j = 0; j < idx_dim.size(); ++j)
216 idx_select[idx_dim[j]] = 1;
220 for(qint32 i = 0; i < idx_select.size(); ++i)
221 if(idx_select[i] == 1)
224 selectedSrc.m_qListHemispheres[h].nuse_tri = countSel;
226 MatrixX3i use_tris_new(countSel,3);
227 MatrixX3d use_tri_cent_new(countSel,3);
228 MatrixX3d use_tri_nn_new(countSel,3);
229 VectorXd use_tri_area_new(countSel);
232 for(qint32 i = 0; i < idx_select.size(); ++i)
234 if(idx_select[i] == 1)
236 use_tris_new.row(countSel) = this->m_qListHemispheres[h].use_tris.row(i);
237 use_tri_cent_new.row(countSel) = this->m_qListHemispheres[h].use_tri_cent.row(i);
238 use_tri_nn_new.row(countSel) = this->m_qListHemispheres[h].use_tri_nn.row(i);
239 use_tri_area_new[countSel] = this->m_qListHemispheres[h].use_tri_area[i];
244 selectedSrc.m_qListHemispheres[h].use_tris = use_tris_new;
245 selectedSrc.m_qListHemispheres[h].use_tri_cent = use_tri_cent_new;
246 selectedSrc.m_qListHemispheres[h].use_tri_nn = use_tri_nn_new;
247 selectedSrc.m_qListHemispheres[h].use_tri_area = use_tri_area_new;
266 bool open_here =
false;
267 if (!p_pStream->device()->isOpen())
269 QList<FiffDirEntry> t_Dir;
270 QString t_sFileName = p_pStream->streamName();
272 QFile t_file(t_sFileName);
274 if(!p_pStream->open(p_Tree, t_Dir))
283 QList<FiffDirTree> spaces = p_Tree.
dir_tree_find(FIFFB_MNE_SOURCE_SPACE);
284 if (spaces.size() == 0)
287 p_pStream->device()->close();
288 std::cout <<
"No source spaces found";
292 for(
int k = 0; k < spaces.size(); ++k)
295 printf(
"\tReading a source space...");
296 MNESourceSpace::read_source_space(p_pStream.data(), spaces[k], p_Hemisphere);
297 printf(
"\t[done]\n" );
299 complete_source_space_info(p_Hemisphere);
301 p_SourceSpace.m_qListHemispheres.append(p_Hemisphere);
306 printf(
"\t%d source spaces read\n", spaces.size());
309 p_pStream->device()->close();
319 double xave = p_Hemisphere.
rr.col(0).sum();
323 hemi = FIFFV_MNE_SURF_LEFT_HEMI;
325 hemi = FIFFV_MNE_SURF_RIGHT_HEMI;
335 for(
int k = 0; k < this->m_qListHemispheres.size(); ++k)
337 if(!this->m_qListHemispheres[k].transform_hemisphere_to(dest,trans))
339 printf(
"Could not transform source space.\n");
351 p_Hemisphere.
clear();
357 p_Hemisphere.
id = FIFFV_MNE_SURF_UNKNOWN;
359 p_Hemisphere.
id = *t_pTag->toInt();
366 p_pStream->device()->close();
367 std::cout <<
"error: Number of vertices not found.";
371 p_Hemisphere.
np = *t_pTag->toInt();
378 p_Hemisphere.
ntri = 0;
380 p_Hemisphere.
ntri = *t_pTag->toInt();
384 p_Hemisphere.
ntri = *t_pTag->toInt();
392 p_pStream->device()->close();
393 std::cout <<
"Coordinate frame information not found.";
405 p_pStream->device()->close();
406 std::cout <<
"Vertex data not found.";
410 p_Hemisphere.
rr = t_pTag->toFloatMatrix().transpose();
411 qint32 rows_rr = p_Hemisphere.
rr.rows();
414 if (rows_rr != p_Hemisphere.
np)
416 p_pStream->device()->close();
417 std::cout <<
"Vertex information is incorrect.";
426 p_pStream->device()->close();
427 std::cout <<
"Vertex normals not found.";
431 p_Hemisphere.
nn = t_pTag->toFloatMatrix().transpose();
432 qint32 rows_nn = p_Hemisphere.
nn.rows();
434 if (rows_nn != p_Hemisphere.
np)
436 p_pStream->device()->close();
437 std::cout <<
"Vertex normal information is incorrect.";
444 if (p_Hemisphere.
ntri > 0)
450 p_pStream->device()->close();
451 std::cout <<
"Triangulation not found.";
456 p_Hemisphere.
tris = t_pTag->toIntMatrix().transpose();
457 p_Hemisphere.
tris -= MatrixXi::Constant(p_Hemisphere.
tris.rows(),3,1);
462 p_Hemisphere.
tris = t_pTag->toIntMatrix().transpose();
463 p_Hemisphere.
tris -= MatrixXi::Constant(p_Hemisphere.
tris.rows(),3,1);
465 if (p_Hemisphere.
tris.rows() != p_Hemisphere.
ntri)
467 p_pStream->device()->close();
468 std::cout <<
"Triangulation information is incorrect.";
474 MatrixXi p_defaultMatrix(0, 0);
475 p_Hemisphere.
tris = p_defaultMatrix;
485 p_Hemisphere.
nuse = 0;
486 p_Hemisphere.
inuse = VectorXi::Zero(p_Hemisphere.
nuse);
487 VectorXi p_defaultVector;
488 p_Hemisphere.
vertno = p_defaultVector;
492 p_Hemisphere.
nuse = *t_pTag->toInt();
495 p_pStream->device()->close();
496 std::cout <<
"Source selection information missing.";
499 p_Hemisphere.
inuse = VectorXi(Map<VectorXi>(t_pTag->toInt(), t_pTag->size()/4, 1));
501 p_Hemisphere.
vertno = VectorXi::Zero(p_Hemisphere.
nuse);
502 if (p_Hemisphere.
inuse.rows() != p_Hemisphere.
np)
504 p_pStream->device()->close();
505 std::cout <<
"Incorrect number of entries in source space selection.";
509 for (
int p = 0; p < p_Hemisphere.
np; ++p)
511 if(p_Hemisphere.
inuse(p) == 1)
513 p_Hemisphere.
vertno(pp) = p;
527 MatrixX3i p_defaultMatrix;
529 p_Hemisphere.
use_tris = p_defaultMatrix;
533 p_Hemisphere.
nuse_tri = *t_pTag1->toInt();
534 p_Hemisphere.
use_tris = t_pTag2->toIntMatrix().transpose();
535 p_Hemisphere.
use_tris -= MatrixXi::Constant(p_Hemisphere.
use_tris.rows(),3,1);
544 VectorXi p_defaultVector;
545 p_Hemisphere.
nearest = p_defaultVector;
546 VectorXd p_defaultFloatVector;
552 p_Hemisphere.
nearest = VectorXi(Map<VectorXi>(t_pTag1->toInt(), t_pTag1->size()/4, 1));
553 p_Hemisphere.
nearest_dist = VectorXd((Map<VectorXf>(t_pTag2->toFloat(), t_pTag1->size()/4, 1)).cast<
double>());
558 printf(
"\tPatch information added...");
567 p_Hemisphere.
dist = SparseMatrix<double>();
572 p_Hemisphere.
dist = t_pTag1->toSparseFloatMatrix();
573 p_Hemisphere.
dist_limit = *t_pTag2->toFloat();
577 SparseMatrix<double> distT = p_Hemisphere.
dist.transpose();
578 p_Hemisphere.
dist += distT;
589 if (p_Hemisphere.
nearest.rows() == 0)
591 p_Hemisphere.
pinfo.clear();
596 printf(
"\tComputing patch statistics...");
598 std::vector< std::pair<int,int> > t_vIndn;
600 for(qint32 i = 0; i < p_Hemisphere.
nearest.rows(); ++i)
602 std::pair<int,int> t_pair(i, p_Hemisphere.
nearest(i));
603 t_vIndn.push_back(t_pair);
605 std::sort(t_vIndn.begin(),t_vIndn.end(), MNEMath::compareIdxValuePairSmallerThan<int> );
607 VectorXi nearest_sorted(t_vIndn.size());
610 std::vector<qint32> t_vfirsti;
611 t_vfirsti.push_back(current);
612 std::vector<qint32> t_vlasti;
614 for(quint32 i = 0; i < t_vIndn.size(); ++i)
616 nearest_sorted[i] = t_vIndn[i].second;
617 if (t_vIndn[current].second != t_vIndn[i].second)
620 t_vlasti.push_back(i-1);
621 t_vfirsti.push_back(current);
624 t_vlasti.push_back(t_vIndn.size()-1);
626 for(quint32 k = 0; k < t_vfirsti.size(); ++k)
628 std::vector<int> t_vIndex;
630 for(
int l = t_vfirsti[k]; l <= t_vlasti[k]; ++l)
631 t_vIndex.push_back(t_vIndn[l].first);
633 std::sort(t_vIndex.begin(),t_vIndex.end());
635 int* t_pV = &t_vIndex[0];
636 Eigen::Map<Eigen::VectorXi> t_vPInfo(t_pV, t_vIndex.size());
638 p_Hemisphere.
pinfo.append(t_vPInfo);
642 std::vector<qint32> patch_verts;
643 patch_verts.reserve(t_vlasti.size());
644 for(quint32 i = 0; i < t_vlasti.size(); ++i)
645 patch_verts.push_back(nearest_sorted[t_vlasti[i]]);
648 std::vector<qint32>::iterator it;
649 for(qint32 i = 0; i < p_Hemisphere.
vertno.size(); ++i)
651 it = std::find(patch_verts.begin(), patch_verts.end(), p_Hemisphere.
vertno[i]);
652 p_Hemisphere.
patch_inds[i] = it-patch_verts.begin();
661 bool MNESourceSpace::complete_source_space_info(
MNEHemisphere& p_Hemisphere)
666 printf(
"\tCompleting triangulation info...");
667 p_Hemisphere.
tri_cent = MatrixX3d::Zero(p_Hemisphere.
ntri,3);
668 p_Hemisphere.
tri_nn = MatrixX3d::Zero(p_Hemisphere.
ntri,3);
669 p_Hemisphere.
tri_area = VectorXd::Zero(p_Hemisphere.
ntri);
675 for (qint32 i = 0; i < p_Hemisphere.
ntri; ++i)
677 for ( qint32 j = 0; j < 3; ++j)
679 k = p_Hemisphere.
tris(i, j);
681 r(j,0) = p_Hemisphere.
rr(k, 0);
682 r(j,1) = p_Hemisphere.
rr(k, 1);
683 r(j,2) = p_Hemisphere.
rr(k, 2);
685 p_Hemisphere.
tri_cent(i, 0) += p_Hemisphere.
rr(k, 0);
686 p_Hemisphere.
tri_cent(i, 1) += p_Hemisphere.
rr(k, 1);
687 p_Hemisphere.
tri_cent(i, 2) += p_Hemisphere.
rr(k, 2);
689 p_Hemisphere.
tri_cent.row(i) /= 3.0f;
692 a = r.row(1) - r.row(0 );
693 b = r.row(2) - r.row(0);
694 p_Hemisphere.
tri_nn(i,0) = a(1)*b(2)-a(2)*b(1);
695 p_Hemisphere.
tri_nn(i,1) = a(2)*b(0)-a(0)*b(2);
696 p_Hemisphere.
tri_nn(i,2) = a(0)*b(1)-a(1)*b(0);
699 size = p_Hemisphere.
tri_nn.row(i)*p_Hemisphere.
tri_nn.row(i).transpose();
700 size = std::pow(size, 0.5f );
702 p_Hemisphere.
tri_area(i) = size/2.0f;
719 printf(
"\tCompleting selection triangulation info...");
727 for (qint32 i = 0; i < p_Hemisphere.
nuse_tri; ++i)
729 for ( qint32 j = 0; j < 3; ++j)
733 r(j,0) = p_Hemisphere.
rr(k, 0);
734 r(j,1) = p_Hemisphere.
rr(k, 1);
735 r(j,2) = p_Hemisphere.
rr(k, 2);
744 a = r.row(1) - r.row(0 );
745 b = r.row(2) - r.row(0);
746 p_Hemisphere.
use_tri_nn(i,0) = a(1)*b(2)-a(2)*b(1);
747 p_Hemisphere.
use_tri_nn(i,1) = a(2)*b(0)-a(0)*b(2);
748 p_Hemisphere.
use_tri_nn(i,2) = a(0)*b(1)-a(1)*b(0);
752 size = std::pow(size, 0.5f );
766 printf(
"\tCompleting triangle and vertex neighboring info...");
767 add_geometry_info(p_Hemisphere);
776 bool MNESourceSpace::add_geometry_info(
MNEHemisphere& p_Hemisphere)
785 for (p = 0; p < p_Hemisphere.
tris.rows(); p++)
786 for (k = 0; k < 3; k++)
789 std::cout<<
"MNESourceSpace::add_geometry_info() - p_Hemisphere.neighbor_tri[100].size(): "<<p_Hemisphere.
neighbor_tri[100].size()<<std::endl;
794 for (k = 0; k < p_Hemisphere.
np; k++) {
795 for (p = 0; p < p_Hemisphere.
neighbor_tri[k].size(); p++) {
797 for (c = 0; c < 3; c++) {
822 std::cout<<
"MNESourceSpace::add_geometry_info() - p_Hemisphere.neighbor_vert[100].size(): "<<p_Hemisphere.
neighbor_tri[100].size()<<std::endl;
831 for(qint32 h = 0; h < m_qListHemispheres.size(); ++h)
833 printf(
"\tWrite a source space... ");
835 m_qListHemispheres[h].writeToStream(p_pStream);
836 p_pStream->
end_block(FIFFB_MNE_SOURCE_SPACE);
839 printf(
"\t%d source spaces written\n", m_qListHemispheres.size());
847 if(m_qListHemispheres.size() > idx)
848 return m_qListHemispheres[idx];
851 qWarning(
"Warning: Index out of bound! Returning last element.");
852 return m_qListHemispheres[m_qListHemispheres.size()-1];
861 if(m_qListHemispheres.size() > idx)
862 return m_qListHemispheres[idx];
865 qWarning(
"Warning: Index out of bound! Returning last element.");
866 return m_qListHemispheres[m_qListHemispheres.size()-1];
875 if(idt.compare(
"lh") == 0)
876 return m_qListHemispheres[0];
877 else if(idt.compare(
"rh") == 0)
878 return m_qListHemispheres[1];
881 qWarning(
"Warning: Identifier is not 'lh' or 'rh'! Returning 'lh'.");
882 return m_qListHemispheres[0];
891 if(idt.compare(
"lh") == 0)
892 return m_qListHemispheres[0];
893 else if(idt.compare(
"rh") == 0)
894 return m_qListHemispheres[1];
897 qWarning(
"Warning: Identifier is not 'lh' or 'rh'! Returning 'lh'.");
898 return m_qListHemispheres[0];
static qint32 find_source_space_hemi(MNEHemisphere &p_Hemisphere)
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
#define FIFF_MNE_SOURCE_SPACE_NORMALS
#define FIFF_MNE_SOURCE_SPACE_NUSE
void writeToStream(FiffStream *p_pStream)
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)
QList< FiffDirTree > dir_tree_find(fiff_int_t p_kind) const
#define FIFF_MNE_SOURCE_SPACE_NTRI
#define FIFF_MNE_SOURCE_SPACE_ID
#define FIFF_MNE_SOURCE_SPACE_TRIANGLES
QSharedPointer< FiffTag > SPtr
#define FIFF_MNE_SOURCE_SPACE_POINTS
static bool readFromStream(FiffStream::SPtr &p_pStream, bool add_geom, FiffDirTree &p_Tree, MNESourceSpace &p_SourceSpace)
QMap< int, QVector< int > > neighbor_vert
QSharedPointer< FiffStream > SPtr
SparseMatrix< double > dist
static VectorXi intersect(const VectorXi &v1, const VectorXi &v2, VectorXi &idx_sel)
void start_block(fiff_int_t kind)
#define FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
MNEHemisphere & operator[](qint32 idx)
Hemisphere provides geometry information.
Directory tree structure.
static VectorXi sort(Matrix< T, Dynamic, 1 > &v, bool desc=true)
#define FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
static bool patch_info(MNEHemisphere &p_Hemisphere)
#define FIFF_MNE_COORD_FRAME
QList< VectorXi > get_vertno() const
MNESourceSpace pick_regions(const QList< Label > &p_qListLabels) const
#define FIFF_MNE_SOURCE_SPACE_NUSE_TRI
#define FIFF_MNE_SOURCE_SPACE_DIST
Coordinate transformation description.
MNESourceSpace class declaration.
#define FIFF_BEM_SURF_TRIANGLES
MNEMath class declaration.
#define FIFF_BEM_SURF_NTRI
QMap< int, QVector< int > > neighbor_tri
#define FIFF_MNE_SOURCE_SPACE_NEAREST
void end_block(fiff_int_t kind)
#define FIFF_MNE_SOURCE_SPACE_SELECTION
#define FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
QList< VectorXi > label_src_vertno_sel(const Label &p_label, VectorXi &src_sel) const