v2.0.0
Loading...
Searching...
No Matches
mne_source_space.cpp
Go to the documentation of this file.
1//=============================================================================================================
36
37//=============================================================================================================
38// INCLUDES
39//=============================================================================================================
40
41#include "mne_source_space.h"
42#include "mne_nearest.h"
43#include "mne_patch_info.h"
44#include "mne_mgh_tag_group.h"
45#include "mne_surface.h"
46#include "mne_hemisphere.h"
47#include "filter_thread_arg.h"
48
49#include <fiff/fiff_coord_trans.h>
50#include <fiff/fiff_constants.h>
51#include <fiff/fiff_sparse_matrix.h>
52#include <fiff/fiff_stream.h>
53#include <fiff/fiff_tag.h>
54
55#include <fiff/fiff_byte_swap.h>
56
57#include <QFile>
58#include <QTextStream>
59#include <QtConcurrent>
60#include <QDebug>
61
62#include <cstring>
63#include <memory>
64
65#define _USE_MATH_DEFINES
66#include <math.h>
67
68using FIFFLIB::FiffCoordTrans;
69
70constexpr int X = 0;
71constexpr int Y = 1;
72constexpr int Z = 2;
73
74constexpr int FAIL = -1;
75constexpr int OK = 0;
76
77constexpr int NNEIGHBORS = 26;
78
79constexpr int CURVATURE_FILE_MAGIC_NUMBER = 16777215;
80
81constexpr int TAG_OLD_MGH_XFORM = 30;
82constexpr int TAG_OLD_COLORTABLE = 1;
83constexpr int TAG_OLD_USEREALRAS = 2;
84constexpr int TAG_USEREALRAS = 4;
85
86//=============================================================================================================
87// USED NAMESPACES
88//=============================================================================================================
89
90using namespace Eigen;
91using namespace FIFFLIB;
92using namespace MNELIB;
93
94//=============================================================================================================
95// FreeSurfer I/O helpers (file-scope, used only from mne_source_space.cpp)
96//=============================================================================================================
97
98namespace {
99
100using PointsT = MNESurfaceOrVolume::PointsT;
101using TrianglesT = MNESurfaceOrVolume::TrianglesT;
102
103//=========================================================================
104// Leaf I/O functions
105//=========================================================================
106
107int read_int3(QFile &in, int &ival)
108/*
109 * Read the strange 3-byte integer
110 */
111{
112 unsigned int s = 0;
113
114 if (in.read(reinterpret_cast<char*>(&s), 3) != 3) {
115 qCritical("read_int3 could not read data");
116 return FAIL;
117 }
118 s = (unsigned int)FIFFLIB::swap_int(s);
119 ival = ((s >> 8) & 0xffffff);
120 return OK;
121}
122
123int read_int(QFile &in, qint32 &ival)
124/*
125 * Read a 32-bit integer
126 */
127{
128 qint32 s ;
129 if (in.read(reinterpret_cast<char*>(&s), sizeof(qint32)) != static_cast<qint64>(sizeof(qint32))) {
130 qCritical("read_int could not read data");
131 return FAIL;
132 }
133 ival = FIFFLIB::swap_int(s);
134 return OK;
135}
136
137int read_int2(QFile &in, int &ival)
138/*
139 * Read int from short
140 */
141{
142 short s ;
143 if (in.read(reinterpret_cast<char*>(&s), sizeof(short)) != static_cast<qint64>(sizeof(short))) {
144 qCritical("read_int2 could not read data");
145 return FAIL;
146 }
147 ival = FIFFLIB::swap_short(s);
148 return OK;
149}
150
151int read_float(QFile &in, float &fval)
152/*
153 * Read float
154 */
155{
156 float f ;
157 if (in.read(reinterpret_cast<char*>(&f), sizeof(float)) != static_cast<qint64>(sizeof(float))) {
158 qCritical("read_float could not read data");
159 return FAIL;
160 }
161 fval = FIFFLIB::swap_float(f);
162 return OK;
163}
164
165int read_long(QFile &in, long long &lval)
166/*
167 * Read a 64-bit integer
168 */
169{
170 long long s ;
171 if (in.read(reinterpret_cast<char*>(&s), sizeof(long long)) != static_cast<qint64>(sizeof(long long))) {
172 qCritical("read_long could not read data");
173 return FAIL;
174 }
175 lval = FIFFLIB::swap_long(s);
176 return OK;
177}
178
179//=========================================================================
180// Leaf validation / copy
181//=========================================================================
182
183int check_vertex(int no, int maxno)
184{
185 if (no < 0 || no > maxno-1) {
186 qCritical("Illegal vertex number %d (max %d).",no,maxno);
187 return FAIL;
188 }
189 return OK;
190}
191
192MNEVolGeom dup_vol_geom(const MNEVolGeom& g)
193{
194 MNEVolGeom dup;
195 dup = g;
196 dup.filename = g.filename;
197 return dup;
198}
199
200//=========================================================================
201// read_vol_geom
202//=========================================================================
203
204std::unique_ptr<MNEVolGeom> read_vol_geom(QFile &fp)
205/*
206 * This the volume geometry reading code from FreeSurfer
207 */
208{
209 char param[64];
210 char eq[2];
211 char buf[256];
212 int vgRead = 0;
213 int counter = 0;
214 qint64 pos = 0;
215
216 auto vg = std::make_unique<MNEVolGeom>();
217
218 while (!fp.atEnd() && counter < 8)
219 {
220 QByteArray lineData = fp.readLine(256);
221 if (lineData.isEmpty())
222 break;
223 const char *line = lineData.constData();
224 if (strlen(line) == 0)
225 break ;
226 sscanf(line, "%s %s %*s", param, eq);
227 if (!strcmp(param, "valid")) {
228 sscanf(line, "%s %s %d", param, eq, &vg->valid);
229 vgRead = 1;
230 counter++;
231 }
232 else if (!strcmp(param, "filename")) {
233 if (sscanf(line, "%s %s %s", param, eq, buf) >= 3)
234 vg->filename = QString::fromUtf8(buf);
235 counter++;
236 }
237 else if (!strcmp(param, "volume")) {
238 sscanf(line, "%s %s %d %d %d",
239 param, eq, &vg->width, &vg->height, &vg->depth);
240 counter++;
241 }
242 else if (!strcmp(param, "voxelsize")) {
243 sscanf(line, "%s %s %f %f %f",
244 param, eq, &vg->xsize, &vg->ysize, &vg->zsize);
245 /*
246 * We like these to be in meters
247 */
248 vg->xsize = vg->xsize/1000.0;
249 vg->ysize = vg->ysize/1000.0;
250 vg->zsize = vg->zsize/1000.0;
251 counter++;
252 }
253 else if (!strcmp(param, "xras")) {
254 sscanf(line, "%s %s %f %f %f",
255 param, eq, vg->x_ras, vg->x_ras+1, vg->x_ras+2);
256 counter++;
257 }
258 else if (!strcmp(param, "yras")) {
259 sscanf(line, "%s %s %f %f %f",
260 param, eq, vg->y_ras, vg->y_ras+1, vg->y_ras+2);
261 counter++;
262 }
263 else if (!strcmp(param, "zras")) {
264 sscanf(line, "%s %s %f %f %f",
265 param, eq, vg->z_ras, vg->z_ras+1, vg->z_ras+2);
266 counter++;
267 }
268 else if (!strcmp(param, "cras")) {
269 sscanf(line, "%s %s %f %f %f",
270 param, eq, vg->c_ras, vg->c_ras+1, vg->c_ras+2);
271 vg->c_ras[0] = vg->c_ras[0]/1000.0;
272 vg->c_ras[1] = vg->c_ras[1]/1000.0;
273 vg->c_ras[2] = vg->c_ras[2]/1000.0;
274 counter++;
275 }
276 /* remember the current position */
277 pos = fp.pos();
278 };
279 if (!fp.atEnd()) { /* we read one more line */
280 if (pos > 0 ) /* if success in getting pos, then */
281 fp.seek(pos); /* restore the position */
282 /* note that this won't allow compression using pipe */
283 }
284 if (!vgRead) {
285 vg = std::make_unique<MNEVolGeom>();
286 }
287 return vg;
288}
289
290//=========================================================================
291// read_tag_data
292//=========================================================================
293
294int read_tag_data(QFile &fp, int tag, long long nbytes, unsigned char *&val, long long &nbytesp)
295/*
296 * Read the data of one tag
297 */
298{
299 size_t snbytes = nbytes;
300
301 val = nullptr;
302 if (nbytes > 0) {
303 auto dum = std::make_unique<unsigned char[]>(nbytes+1);
304 if (fp.read(reinterpret_cast<char*>(dum.get()), nbytes) != static_cast<qint64>(snbytes)) {
305 qCritical("Failed to read %d bytes of tag data",static_cast<int>(nbytes));
306 return FAIL;
307 }
308 dum[nbytes] = '\0'; /* Ensure null termination */
309 val = dum.release();
310 nbytesp = nbytes;
311 }
312 else { /* Need to handle special cases */
313 if (tag == TAG_OLD_SURF_GEOM) {
314 auto g = read_vol_geom(fp);
315 if (!g)
316 return FAIL;
317 /*
318 * Serialize MNEVolGeom as POD fields followed by a
319 * null-terminated UTF-8 filename. We must NOT memcpy the
320 * entire MNEVolGeom because it contains a QString member
321 * whose internal pointers become dangling when *g is
322 * destroyed at the end of this block.
323 */
324 struct VolGeomPOD {
325 int valid;
326 int width, height, depth;
327 float xsize, ysize, zsize;
328 float x_ras[3], y_ras[3], z_ras[3];
329 float c_ras[3];
330 };
331 QByteArray fn = g->filename.toUtf8();
332 size_t totalSize = sizeof(VolGeomPOD) + fn.size() + 1;
333 auto buf = std::make_unique<unsigned char[]>(totalSize);
334 VolGeomPOD pod;
335 pod.valid = g->valid;
336 pod.width = g->width; pod.height = g->height; pod.depth = g->depth;
337 pod.xsize = g->xsize; pod.ysize = g->ysize; pod.zsize = g->zsize;
338 std::memcpy(pod.x_ras, g->x_ras, 3 * sizeof(float));
339 std::memcpy(pod.y_ras, g->y_ras, 3 * sizeof(float));
340 std::memcpy(pod.z_ras, g->z_ras, 3 * sizeof(float));
341 std::memcpy(pod.c_ras, g->c_ras, 3 * sizeof(float));
342 std::memcpy(buf.get(), &pod, sizeof(VolGeomPOD));
343 std::memcpy(buf.get() + sizeof(VolGeomPOD), fn.constData(), fn.size() + 1);
344 val = buf.release();
345 nbytesp = static_cast<long long>(totalSize);
346 }
347 else if (tag == TAG_OLD_USEREALRAS || tag == TAG_USEREALRAS) {
348 auto vi = std::make_unique<int[]>(1);
349 if (read_int(fp, vi[0]) == FAIL)
350 return FAIL;
351 val = reinterpret_cast<unsigned char *>(vi.release());
352 nbytesp = sizeof(int);
353 }
354 else {
355 qWarning("Encountered an unknown tag with no length specification : %d\n",tag);
356 val = nullptr;
357 nbytesp = 0;
358 }
359 }
360 return OK;
361}
362
363//=========================================================================
364// add_mgh_tag_to_group
365//=========================================================================
366
367void add_mgh_tag_to_group(std::optional<MNEMghTagGroup>& g, int tag, long long len, unsigned char *data)
368{
369 if (!g)
370 g = MNEMghTagGroup();
371 auto new_tag = std::make_unique<MNEMghTag>();
372 new_tag->tag = tag;
373 new_tag->len = len;
374 new_tag->data = QByteArray(reinterpret_cast<const char*>(data), static_cast<int>(len));
375 delete[] data;
376 g->tags.push_back(std::move(new_tag));
377}
378
379//=========================================================================
380// read_next_tag
381//=========================================================================
382
383int read_next_tag(QFile &fp, int &tagp, long long &lenp, unsigned char *&datap)
384/*
385 * Read the next tag in the file
386 */
387{
388 int ilen,tag;
389 long long len;
390
391 if (read_int(fp,tag) == FAIL) {
392 tagp = 0;
393 return OK;
394 }
395 if (fp.atEnd()) {
396 tagp = 0;
397 return OK;
398 }
399 switch (tag) {
400 case TAG_OLD_MGH_XFORM: /* This is obviously a burden of the past */
401 if (read_int(fp,ilen) == FAIL)
402 return FAIL;
403 len = ilen - 1;
404 break ;
405 case TAG_OLD_SURF_GEOM:
406 case TAG_OLD_USEREALRAS:
407 case TAG_OLD_COLORTABLE:
408 len = 0 ;
409 break ;
410 default:
411 if (read_long(fp,len) == FAIL)
412 return FAIL;
413 break;
414 }
415 lenp = len;
416 tagp = tag;
417 if (read_tag_data(fp,tag,len,datap,lenp) == FAIL)
418 return FAIL;
419 return OK;
420}
421
422//=========================================================================
423// read_mgh_tags
424//=========================================================================
425
426int read_mgh_tags(QFile &fp, std::optional<MNEMghTagGroup>& tagsp)
427/*
428 * Read all the tags from the file
429 */
430{
431 long long len;
432 int tag;
433 unsigned char *tag_data;
434
435 while (1) {
436 if (read_next_tag(fp,tag,len,tag_data) == FAIL)
437 return FAIL;
438 if (tag == 0)
439 break;
440 add_mgh_tag_to_group(tagsp,tag,len,tag_data);
441 }
442 return OK;
443}
444
445//=========================================================================
446// read_curvature_file
447//=========================================================================
448
449int read_curvature_file(const QString& fname,
450 Eigen::VectorXf& curv)
451
452{
453 QFile fp(fname);
454 int magic;
455
456 float curvmin,curvmax;
457 int ncurv = 0;
458 int nface,val_pervert;
459 int val,k;
460 float fval;
461
462 if (!fp.open(QIODevice::ReadOnly)) {
463 qCritical() << fname;
464 curv.resize(0); return FAIL;
465 }
466 if (read_int3(fp,magic) != 0) {
467 qCritical() << "Bad magic in" << fname;
468 curv.resize(0); return FAIL;
469 }
470 if (magic == CURVATURE_FILE_MAGIC_NUMBER) { /* A new-style curvature file */
471 /*
472 * How many and faces
473 */
474 if (read_int(fp,ncurv) != 0) {
475 curv.resize(0); return FAIL;
476 }
477 if (read_int(fp,nface) != 0) {
478 curv.resize(0); return FAIL;
479 }
480#ifdef DEBUG
481 qInfo("nvert = %d nface = %d\n",ncurv,nface);
482#endif
483 if (read_int(fp,val_pervert) != 0) {
484 curv.resize(0); return FAIL;
485 }
486 if (val_pervert != 1) {
487 qCritical("Values per vertex not equal to one.");
488 curv.resize(0); return FAIL;
489 }
490 /*
491 * Read the curvature values
492 */
493 curv.resize(ncurv);
494 curvmin = curvmax = 0.0;
495 for (k = 0; k < ncurv; k++) {
496 if (read_float(fp,fval) != 0) {
497 curv.resize(0); return FAIL;
498 }
499 curv[k] = fval;
500 if (curv[k] > curvmax)
501 curvmax = curv[k];
502 if (curv[k] < curvmin)
503 curvmin = curv[k];
504 }
505 }
506 else { /* An old-style curvature file */
507 ncurv = magic;
508 /*
509 * How many vertices
510 */
511 if (read_int3(fp,nface) != 0) {
512 curv.resize(0); return FAIL;
513 }
514#ifdef DEBUG
515 qInfo("nvert = %d nface = %d\n",ncurv,nface);
516#endif
517 /*
518 * Read the curvature values
519 */
520 curv.resize(ncurv);
521 curvmin = curvmax = 0.0;
522 for (k = 0; k < ncurv; k++) {
523 if (read_int2(fp,val) != 0) {
524 curv.resize(0); return FAIL;
525 }
526 curv[k] = static_cast<float>(val)/100.0;
527 if (curv[k] > curvmax)
528 curvmax = curv[k];
529 if (curv[k] < curvmin)
530 curvmin = curv[k];
531
532 }
533 }
534#ifdef DEBUG
535 qInfo("Curvature range: %f...%f\n",curvmin,curvmax);
536#endif
537 return OK;
538}
539
540//=========================================================================
541// read_triangle_file
542//=========================================================================
543
544int read_triangle_file(const QString& fname,
545 PointsT& vertices,
546 TrianglesT& triangles,
547 std::optional<MNEMghTagGroup>* tagsp)
548/*
549 * Read the FS triangulated surface
550 */
551{
552 QFile fp(fname);
553 int magic;
554 char c;
555
556 qint32 nvert,ntri,nquad;
557 PointsT vert;
558 TrianglesT tri;
559 int k,p;
560 int quad[4];
561 int val;
562 int which;
563
564 if (!fp.open(QIODevice::ReadOnly)) {
565 qCritical() << fname;
566 return FAIL;
567 }
568 if (read_int3(fp,magic) != 0) {
569 qCritical() << "Bad magic in" << fname;
570 return FAIL;
571 }
572 if (magic != TRIANGLE_FILE_MAGIC_NUMBER &&
573 magic != QUAD_FILE_MAGIC_NUMBER &&
574 magic != NEW_QUAD_FILE_MAGIC_NUMBER) {
575 qCritical() << "Bad magic in" << fname;
576 return FAIL;
577 }
578 if (magic == TRIANGLE_FILE_MAGIC_NUMBER) {
579 /*
580 * Get the comment
581 */
582 qInfo("Triangle file : ");
583 for (fp.getChar(&c); c != '\n'; fp.getChar(&c)) {
584 if (fp.atEnd()) {
585 qCritical()<<"Bad triangle file.";
586 return FAIL;
587 }
588 putc(c,stderr);
589 }
590 fp.getChar(&c);
591 /*
592 * How many vertices and triangles?
593 */
594 if (read_int(fp,nvert) != 0)
595 return FAIL;
596 if (read_int(fp,ntri) != 0)
597 return FAIL;
598 qInfo(" nvert = %d ntri = %d\n",nvert,ntri);
599 vert.resize(nvert, 3);
600 tri.resize(ntri, 3);
601 /*
602 * Read the vertices
603 */
604 for (k = 0; k < nvert; k++) {
605 if (read_float(fp,vert(k,0)) != 0)
606 return FAIL;
607 if (read_float(fp,vert(k,1)) != 0)
608 return FAIL;
609 if (read_float(fp,vert(k,2)) != 0)
610 return FAIL;
611 }
612 /*
613 * Read the triangles
614 */
615 for (k = 0; k < ntri; k++) {
616 if (read_int(fp,tri(k,0)) != 0)
617 return FAIL;
618 if (check_vertex(tri(k,0),nvert) != OK)
619 return FAIL;
620 if (read_int(fp,tri(k,1)) != 0)
621 return FAIL;
622 if (check_vertex(tri(k,1),nvert) != OK)
623 return FAIL;
624 if (read_int(fp,tri(k,2)) != 0)
625 return FAIL;
626 if (check_vertex(tri(k,2),nvert) != OK)
627 return FAIL;
628 }
629 }
630 else if (magic == QUAD_FILE_MAGIC_NUMBER ||
631 magic == NEW_QUAD_FILE_MAGIC_NUMBER) {
632 if (read_int3(fp,nvert) != 0)
633 return FAIL;
634 if (read_int3(fp,nquad) != 0)
635 return FAIL;
636 qInfo("%s file : nvert = %d nquad = %d\n",
637 magic == QUAD_FILE_MAGIC_NUMBER ? "Quad" : "New quad",
638 nvert,nquad);
639 vert.resize(nvert, 3);
640 if (magic == QUAD_FILE_MAGIC_NUMBER) {
641 for (k = 0; k < nvert; k++) {
642 if (read_int2(fp,val) != 0)
643 return FAIL;
644 vert(k,0) = val/100.0;
645 if (read_int2(fp,val) != 0)
646 return FAIL;
647 vert(k,1) = val/100.0;
648 if (read_int2(fp,val) != 0)
649 return FAIL;
650 vert(k,2) = val/100.0;
651 }
652 }
653 else { /* NEW_QUAD_FILE_MAGIC_NUMBER */
654 for (k = 0; k < nvert; k++) {
655 if (read_float(fp,vert(k,0)) != 0)
656 return FAIL;
657 if (read_float(fp,vert(k,1)) != 0)
658 return FAIL;
659 if (read_float(fp,vert(k,2)) != 0)
660 return FAIL;
661 }
662 }
663 ntri = 2*nquad;
664 tri.resize(ntri, 3);
665 for (k = 0, ntri = 0; k < nquad; k++) {
666 for (p = 0; p < 4; p++) {
667 if (read_int3(fp,quad[p]) != 0)
668 return FAIL;
669 }
670
671 /*
672 * The randomization is borrowed from FreeSurfer code
673 * Strange...
674 */
675#define EVEN(n) ((((n) / 2) * 2) == n)
676#ifdef FOO
677#define WHICH_FACE_SPLIT(vno0, vno1) \
678 (1*nearbyint(sqrt(1.9*vno0) + sqrt(3.5*vno1)))
679
680 which = WHICH_FACE_SPLIT(quad[0], quad[1]) ;
681#endif
682 which = quad[0];
683 /*
684 qInfo("%f ",sqrt(1.9*quad[0]) + sqrt(3.5*quad[1]));
685 */
686
687 if (EVEN(which)) {
688 tri(ntri,0) = quad[0];
689 tri(ntri,1) = quad[1];
690 tri(ntri,2) = quad[3];
691 ntri++;
692
693 tri(ntri,0) = quad[2];
694 tri(ntri,1) = quad[3];
695 tri(ntri,2) = quad[1];
696 ntri++;
697 }
698 else {
699 tri(ntri,0) = quad[0];
700 tri(ntri,1) = quad[1];
701 tri(ntri,2) = quad[2];
702 ntri++;
703
704 tri(ntri,0) = quad[0];
705 tri(ntri,1) = quad[2];
706 tri(ntri,2) = quad[3];
707 ntri++;
708 }
709 }
710 }
711 /*
712 * Optionally read the tags
713 */
714 if (tagsp) {
715 std::optional<MNEMghTagGroup> tags;
716 if (read_mgh_tags(fp, tags) == FAIL) {
717 return FAIL;
718 }
719 *tagsp = std::move(tags);
720 }
721 /*
722 * Convert mm to m and store as Eigen matrices
723 */
724 vert /= 1000.0f;
725 vertices = std::move(vert);
726 triangles = std::move(tri);
727 return OK;
728}
729
730//=========================================================================
731// get_volume_geom_from_tag
732//=========================================================================
733
734std::optional<MNEVolGeom> get_volume_geom_from_tag(const MNEMghTagGroup *tagsp)
735{
736 if (!tagsp)
737 return std::nullopt;
738
739 struct VolGeomPOD {
740 int valid;
741 int width, height, depth;
742 float xsize, ysize, zsize;
743 float x_ras[3], y_ras[3], z_ras[3];
744 float c_ras[3];
745 };
746
747 for (const auto &t : tagsp->tags) {
748 if (t->tag == TAG_OLD_SURF_GEOM) {
749 if (t->len < static_cast<long long>(sizeof(VolGeomPOD)))
750 return std::nullopt;
751
752 const unsigned char *d = reinterpret_cast<const unsigned char *>(t->data.constData());
753 VolGeomPOD pod;
754 std::memcpy(&pod, d, sizeof(VolGeomPOD));
755
756 MNEVolGeom result;
757 result.valid = pod.valid;
758 result.width = pod.width; result.height = pod.height; result.depth = pod.depth;
759 result.xsize = pod.xsize; result.ysize = pod.ysize; result.zsize = pod.zsize;
760 std::memcpy(result.x_ras, pod.x_ras, 3 * sizeof(float));
761 std::memcpy(result.y_ras, pod.y_ras, 3 * sizeof(float));
762 std::memcpy(result.z_ras, pod.z_ras, 3 * sizeof(float));
763 std::memcpy(result.c_ras, pod.c_ras, 3 * sizeof(float));
764
765 if (t->len > static_cast<long long>(sizeof(VolGeomPOD)))
766 result.filename = QString::fromUtf8(
767 reinterpret_cast<const char *>(d + sizeof(VolGeomPOD)));
768
769 return result;
770 }
771 }
772 return std::nullopt;
773}
774
775} // anonymous namespace
776
777//=============================================================================================================
778// DEFINE MEMBER METHODS
779//=============================================================================================================
780
781MNESourceSpace::MNESourceSpace(int np)
782{
783 this->np = np;
784 if (np > 0) {
785 rr = PointsT::Zero(np, 3);
786 nn = NormalsT::Zero(np, 3);
787 inuse = VectorXi::Zero(np);
788 vertno = VectorXi::Zero(np);
789 }
790 nuse = 0;
791 ntri = 0;
792 tot_area = 0.0;
793
794 nuse_tri = 0;
795
796 // tris, use_tris are std::vector<MNETriangle> — default-constructed empty
797
798 // neighbor_tri, nneighbor_tri, curv, val,
799 // neighbor_vert, nneighbor_vert, vert_dist
800 // are Eigen/std::vector types — default-constructed empty
801
802 coord_frame = FIFFV_COORD_MRI;
803 id = FIFFV_MNE_SURF_UNKNOWN;
804 subject = "";
805 type = FIFFV_MNE_SPACE_SURFACE;
806
807 // nearest is std::vector<MNENearest> — default-constructed empty
808 // patches is std::vector<optional<MNEPatchInfo>> — default-constructed empty
809
810 dist = FIFFLIB::FiffSparseMatrix();
811 dist_limit = -1.0;
812
813 voxel_surf_RAS_t.reset();
814 vol_dims[0] = vol_dims[1] = vol_dims[2] = 0;
815
816 MRI_volume = "";
817 MRI_surf_RAS_RAS_t.reset();
818 MRI_voxel_surf_RAS_t.reset();
819 MRI_vol_dims[0] = MRI_vol_dims[1] = MRI_vol_dims[2] = 0;
820 interpolator.reset();
821
822 vol_geom.reset();
823 mgh_tags.reset();
824
825 cm[0] = cm[1] = cm[2] = 0.0;
826}
827
828//=============================================================================================================
829
833
834//=============================================================================================================
835
836MNESourceSpace::SPtr MNESourceSpace::clone() const
837{
838 // Base class clone — creates a MNESourceSpace with the same fields.
839 // Derived classes (e.g., MNEHemisphere) override this to preserve their type.
840 auto copy = std::make_shared<MNESourceSpace>(this->np);
841 copy->type = this->type;
842 copy->id = this->id;
843 copy->np = this->np;
844 copy->ntri = this->ntri;
845 copy->coord_frame = this->coord_frame;
846 copy->rr = this->rr;
847 copy->nn = this->nn;
848 copy->nuse = this->nuse;
849 copy->inuse = this->inuse;
850 copy->vertno = this->vertno;
851 copy->itris = this->itris;
852 copy->use_itris = this->use_itris;
853 copy->nuse_tri = this->nuse_tri;
854 copy->dist_limit = this->dist_limit;
855 copy->dist = this->dist;
856 copy->nearest = this->nearest;
857 copy->neighbor_tri = this->neighbor_tri;
858 copy->neighbor_vert = this->neighbor_vert;
859 return copy;
860}
861
862//=============================================================================================================
863
864void MNESourceSpace::enable_all_sources()
865{
866 int k;
867 for (k = 0; k < np; k++)
868 inuse[k] = 1;
869 nuse = np;
870 return;
871}
872
873//=============================================================================================================
874
875bool MNESourceSpace::is_left_hemi() const
876/*
877 * Left or right hemisphere?
878 */
879{
880 int k;
881 float xave;
882
883 for (k = 0, xave = 0.0; k < np; k++)
884 xave += rr(k,0);
885 if (xave < 0.0)
886 return true;
887 else
888 return false;
889}
890
891//=============================================================================================================
892
893qint32 MNESourceSpace::find_source_space_hemi() const
894{
895 double xave = rr.col(0).sum();
896 if (xave < 0)
897 return FIFFV_MNE_SURF_LEFT_HEMI;
898 else
899 return FIFFV_MNE_SURF_RIGHT_HEMI;
900}
901
902//=============================================================================================================
903
904void MNESourceSpace::update_inuse(Eigen::VectorXi new_inuse)
905/*
906 * Update the active vertices
907 */
908{
909 int k,p,nuse_count;
910
911 inuse = std::move(new_inuse);
912
913 for (k = 0, nuse_count = 0; k < np; k++)
914 if (inuse[k])
915 nuse_count++;
916
917 nuse = nuse_count;
918 if (nuse > 0) {
919 vertno.conservativeResize(nuse);
920 for (k = 0, p = 0; k < np; k++)
921 if (inuse[k])
922 vertno[p++] = k;
923 }
924 else {
925 vertno.resize(0);
926 }
927 return;
928}
929
930//=============================================================================================================
931
932int MNESourceSpace::transform_source_space(const FiffCoordTrans& t)
933/*
934 * Transform source space data into another coordinate frame
935 */
936{
937 int k;
938 if (coord_frame == t.to)
939 return OK;
940 if (coord_frame != t.from) {
941 qCritical("Coordinate transformation does not match with the source space coordinate system.");
942 return FAIL;
943 }
944 for (k = 0; k < np; k++) {
945 FiffCoordTrans::apply_trans(&rr(k,0),t,FIFFV_MOVE);
946 FiffCoordTrans::apply_trans(&nn(k,0),t,FIFFV_NO_MOVE);
947 }
948 if (!tris.empty()) {
949 for (k = 0; k < ntri; k++)
950 FiffCoordTrans::apply_trans(tris[k].nn.data(),t,FIFFV_NO_MOVE);
951 }
952 coord_frame = t.to;
953 return OK;
954}
955
956//=============================================================================================================
957
958int MNESourceSpace::add_patch_stats()
959{
960 MNENearest* nearest_data = nearest.data();
961 MNENearest* this_patch;
962 std::vector<std::optional<MNEPatchInfo>> pinfo(nuse);
963 int nave,p,q,k;
964
965 qInfo("Computing patch statistics...\n");
966 if (neighbor_tri.empty())
967 if (add_geometry_info(false) != OK)
968 return FAIL;
969
970 if (nearest.empty()) {
971 qCritical("The patch information is not available.");
972 return FAIL;
973 }
974 if (nuse == 0) {
975 patches.clear();
976 return OK;
977 }
978 /*
979 * Calculate the average normals and the patch areas
980 */
981 qInfo("\tareas, average normals, and mean deviations...");
982 std::sort(nearest.begin(), nearest.end(),
983 [](const MNENearest& a, const MNENearest& b) { return a.nearest < b.nearest; });
984 nearest_data = nearest.data(); // refresh after sort
985 nave = 1;
986 for (p = 1, q = 0; p < np; p++) {
987 if (nearest_data[p].nearest != nearest_data[p-1].nearest) {
988 if (nave == 0) {
989 qCritical("No vertices belong to the patch of vertex %d",nearest_data[p-1].nearest);
990 return FAIL;
991 }
992 if (q < nuse && vertno[q] == nearest_data[p-1].nearest) { /* Some source space points may have been omitted since
993 * the patch information was computed */
994 pinfo[q] = MNEPatchInfo();
995 pinfo[q]->vert = nearest_data[p-1].nearest;
996 this_patch = nearest_data+p-nave;
997 pinfo[q]->memb_vert.resize(nave);
998 for (k = 0; k < nave; k++) {
999 pinfo[q]->memb_vert[k] = this_patch[k].vert;
1000 this_patch[k].patch = &(*pinfo[q]);
1001 }
1002 pinfo[q]->calculate_area(this);
1003 pinfo[q]->calculate_normal_stats(this);
1004 q++;
1005 }
1006 nave = 0;
1007 }
1008 nave++;
1009 }
1010 if (nave == 0) {
1011 qCritical("No vertices belong to the patch of vertex %d",nearest_data[p-1].nearest);
1012 return FAIL;
1013 }
1014 if (q < nuse && vertno[q] == nearest_data[p-1].nearest) {
1015 pinfo[q] = MNEPatchInfo();
1016 pinfo[q]->vert = nearest_data[p-1].nearest;
1017 this_patch = nearest_data+p-nave;
1018 pinfo[q]->memb_vert.resize(nave);
1019 for (k = 0; k < nave; k++) {
1020 pinfo[q]->memb_vert[k] = this_patch[k].vert;
1021 this_patch[k].patch = &(*pinfo[q]);
1022 }
1023 pinfo[q]->calculate_area(this);
1024 pinfo[q]->calculate_normal_stats(this);
1025 q++;
1026 }
1027 qInfo(" %d/%d [done]\n",q,nuse);
1028
1029 patches = std::move(pinfo);
1030
1031 return OK;
1032}
1033
1034//=============================================================================================================
1035
1036void MNESourceSpace::rearrange_source_space()
1037{
1038 int k,p;
1039
1040 for (k = 0, nuse = 0; k < np; k++)
1041 if (inuse[k])
1042 nuse++;
1043
1044 if (nuse == 0) {
1045 vertno.resize(0);
1046 }
1047 else {
1048 vertno.conservativeResize(nuse);
1049 for (k = 0, p = 0; k < np; k++)
1050 if (inuse[k])
1051 vertno[p++] = k;
1052 }
1053 if (!nearest.empty())
1054 add_patch_stats();
1055 return;
1056}
1057
1058//=============================================================================================================
1059
1060std::unique_ptr<MNESourceSpace> MNESourceSpace::create_source_space(int np)
1061/*
1062 * Create a new source space and all associated data
1063 */
1064{
1065 auto res = std::make_unique<MNESourceSpace>();
1066 res->np = np;
1067 if (np > 0) {
1068 res->rr = PointsT::Zero(np, 3);
1069 res->nn = NormalsT::Zero(np, 3);
1070 res->inuse = VectorXi::Zero(np);
1071 res->vertno = VectorXi::Zero(np);
1072 }
1073 res->nuse = 0;
1074 res->ntri = 0;
1075 res->tot_area = 0.0;
1076
1077 res->nuse_tri = 0;
1078
1079 res->sigma = -1.0;
1080 res->coord_frame = FIFFV_COORD_MRI;
1081 res->id = FIFFV_MNE_SURF_UNKNOWN;
1082 res->subject.clear();
1083 res->type = FIFFV_MNE_SPACE_SURFACE;
1084
1085 res->dist = FIFFLIB::FiffSparseMatrix();
1086 res->dist_limit = -1.0;
1087
1088 res->voxel_surf_RAS_t.reset();
1089 res->vol_dims[0] = res->vol_dims[1] = res->vol_dims[2] = 0;
1090
1091 res->MRI_volume.clear();
1092 res->MRI_surf_RAS_RAS_t.reset();
1093 res->MRI_voxel_surf_RAS_t.reset();
1094 res->MRI_vol_dims[0] = res->MRI_vol_dims[1] = res->MRI_vol_dims[2] = 0;
1095 res->interpolator.reset();
1096
1097 res->vol_geom.reset();
1098 res->mgh_tags.reset();
1099
1100 res->cm[0] = res->cm[1] = res->cm[2] = 0.0;
1101
1102 return res;
1103}
1104
1105//=============================================================================================================
1106
1107std::unique_ptr<MNESourceSpace> MNESourceSpace::load_surface(const QString& surf_file,
1108 const QString& curv_file)
1109{
1110 return load_surface_geom(surf_file,curv_file,true,true);
1111}
1112
1113//=============================================================================================================
1114
1115std::unique_ptr<MNESourceSpace> MNESourceSpace::load_surface_geom(const QString& surf_file,
1116 const QString& curv_file,
1117 bool add_geometry,
1118 bool check_too_many_neighbors)
1119 /*
1120 * Load the surface and add the geometry information
1121 */
1122{
1123 int k;
1124 std::unique_ptr<MNESourceSpace> s;
1125 std::optional<MNEMghTagGroup> tags;
1126 Eigen::VectorXf curvs;
1127 PointsT verts;
1128 TrianglesT tris;
1129
1130 if (read_triangle_file(surf_file,
1131 verts,
1132 tris,
1133 &tags) == -1)
1134 return nullptr;
1135
1136 if (!curv_file.isEmpty()) {
1137 if (read_curvature_file(curv_file, curvs) == -1)
1138 return nullptr;
1139 if (curvs.size() != verts.rows()) {
1140 qCritical()<<"Incorrect number of vertices in the curvature file.";
1141 return nullptr;
1142 }
1143 }
1144
1145 s = std::make_unique<MNESourceSpace>(0);
1146 s->rr = std::move(verts);
1147 s->itris = std::move(tris);
1148 s->ntri = s->itris.rows();
1149 s->np = s->rr.rows();
1150 if (curvs.size() > 0) {
1151 s->curv = std::move(curvs);
1152 }
1153 s->val = Eigen::VectorXf::Zero(s->np);
1154 if (add_geometry) {
1155 if (check_too_many_neighbors) {
1156 if (s->add_geometry_info(true) != OK)
1157 return nullptr;
1158 }
1159 else {
1160 if (s->add_geometry_info2(true) != OK)
1161 return nullptr;
1162 }
1163 }
1164 else if (s->nn.rows() == 0) { /* Normals only */
1165 if (s->add_vertex_normals() != OK)
1166 return nullptr;
1167 }
1168 else
1169 s->add_triangle_data();
1170 s->nuse = s->np;
1171 s->inuse = Eigen::VectorXi::Ones(s->np);
1172 s->vertno = Eigen::VectorXi::LinSpaced(s->np, 0, s->np - 1);
1173 s->mgh_tags = std::move(tags);
1174 s->vol_geom = get_volume_geom_from_tag(s->mgh_tags ? &(*s->mgh_tags) : nullptr);
1175
1176 return s;
1177}
1178
1179//=============================================================================================================
1180
1181static std::optional<FiffCoordTrans> make_voxel_ras_trans(const Eigen::Vector3f& r0,
1182 const Eigen::Vector3f& x_ras,
1183 const Eigen::Vector3f& y_ras,
1184 const Eigen::Vector3f& z_ras,
1185 const Eigen::Vector3f& voxel_size)
1186{
1187 Eigen::Matrix3f rot;
1188 rot.row(0) = x_ras.transpose() * voxel_size[0];
1189 rot.row(1) = y_ras.transpose() * voxel_size[1];
1190 rot.row(2) = z_ras.transpose() * voxel_size[2];
1191
1192 return FiffCoordTrans(FIFFV_MNE_COORD_MRI_VOXEL, FIFFV_COORD_MRI, rot, r0);
1193}
1194
1195MNESourceSpace* MNESourceSpace::make_volume_source_space(const MNESurface& surf, float grid, float exclude, float mindist)
1196/*
1197 * Make a source space which covers the volume bounded by surf
1198 */
1199{
1200 Eigen::Vector3f minV, maxV, cm;
1201 int minn[3],maxn[3];
1202 float maxdist,dist;
1203 int k,c;
1204 std::unique_ptr<MNESourceSpace> sp;
1205 int np,nplane,nrow;
1206 int nneigh;
1207 int x,y,z;
1208 /*
1209 * Figure out the grid size
1210 */
1211 cm.setZero();
1212 minV = maxV = surf.rr.row(0).transpose();
1213
1214 for (k = 0; k < surf.np; k++) {
1215 Eigen::Vector3f node = surf.rr.row(k).transpose();
1216 cm += node;
1217 minV = minV.cwiseMin(node);
1218 maxV = maxV.cwiseMax(node);
1219 }
1220 cm /= static_cast<float>(surf.np);
1221 /*
1222 * Define the sphere which fits the surface
1223 */
1224 maxdist = 0.0;
1225 for (k = 0; k < surf.np; k++) {
1226 dist = (surf.rr.row(k).transpose() - cm).norm();
1227 if (dist > maxdist)
1228 maxdist = dist;
1229 }
1230 qInfo("FsSurface CM = (%6.1f %6.1f %6.1f) mm\n",
1231 1000*cm[X], 1000*cm[Y], 1000*cm[Z]);
1232 qInfo("FsSurface fits inside a sphere with radius %6.1f mm\n",1000*maxdist);
1233 qInfo("FsSurface extent:\n"
1234 "\tx = %6.1f ... %6.1f mm\n"
1235 "\ty = %6.1f ... %6.1f mm\n"
1236 "\tz = %6.1f ... %6.1f mm\n",
1237 1000*minV[X],1000*maxV[X],
1238 1000*minV[Y],1000*maxV[Y],
1239 1000*minV[Z],1000*maxV[Z]);
1240 for (c = 0; c < 3; c++) {
1241 if (maxV[c] > 0)
1242 maxn[c] = floor(std::fabs(maxV[c])/grid)+1;
1243 else
1244 maxn[c] = -floor(std::fabs(maxV[c])/grid)-1;
1245 if (minV[c] > 0)
1246 minn[c] = floor(std::fabs(minV[c])/grid)+1;
1247 else
1248 minn[c] = -floor(std::fabs(minV[c])/grid)-1;
1249 }
1250 qInfo("Grid extent:\n"
1251 "\tx = %6.1f ... %6.1f mm\n"
1252 "\ty = %6.1f ... %6.1f mm\n"
1253 "\tz = %6.1f ... %6.1f mm\n",
1254 1000*(minn[0]*grid),1000*(maxn[0]*grid),
1255 1000*(minn[1]*grid),1000*(maxn[1]*grid),
1256 1000*(minn[2]*grid),1000*(maxn[2]*grid));
1257 /*
1258 * Now make the initial grid
1259 */
1260 np = 1;
1261 for (c = 0; c < 3; c++)
1262 np = np*(maxn[c]-minn[c]+1);
1263 nplane = (maxn[0]-minn[0]+1)*(maxn[1]-minn[1]+1);
1264 nrow = (maxn[0]-minn[0]+1);
1265 sp = MNESourceSpace::create_source_space(np);
1266 sp->type = MNE_SOURCE_SPACE_VOLUME;
1267 sp->nneighbor_vert = Eigen::VectorXi::Constant(sp->np, NNEIGHBORS);
1268 sp->neighbor_vert.resize(sp->np);
1269 for (k = 0; k < sp->np; k++) {
1270 sp->inuse[k] = 1;
1271 sp->vertno[k] = k;
1272 sp->nn(k,0) = sp->nn(k,1) = 0.0; /* Source orientation is immaterial */
1273 sp->nn(k,2) = 1.0;
1274 sp->neighbor_vert[k] = Eigen::VectorXi::Constant(NNEIGHBORS, -1);
1275 sp->nuse++;
1276 }
1277 for (k = 0, z = minn[2]; z <= maxn[2]; z++) {
1278 for (y = minn[1]; y <= maxn[1]; y++) {
1279 for (x = minn[0]; x <= maxn[0]; x++, k++) {
1280 sp->rr(k,0) = x*grid;
1281 sp->rr(k,1) = y*grid;
1282 sp->rr(k,2) = z*grid;
1283 /*
1284 * Figure out the neighborhood:
1285 * 6-neighborhood first
1286 */
1287 Eigen::VectorXi& neigh = sp->neighbor_vert[k];
1288 if (z > minn[2])
1289 neigh[0] = k - nplane;
1290 if (x < maxn[0])
1291 neigh[1] = k + 1;
1292 if (y < maxn[1])
1293 neigh[2] = k + nrow;
1294 if (x > minn[0])
1295 neigh[3] = k - 1;
1296 if (y > minn[1])
1297 neigh[4] = k - nrow;
1298 if (z < maxn[2])
1299 neigh[5] = k + nplane;
1300 /*
1301 * Then the rest to complete the 26-neighborhood
1302 * First the plane below
1303 */
1304 if (z > minn[2]) {
1305 if (x < maxn[0]) {
1306 neigh[6] = k + 1 - nplane;
1307 if (y < maxn[1])
1308 neigh[7] = k + 1 + nrow - nplane;
1309 }
1310 if (y < maxn[1])
1311 neigh[8] = k + nrow - nplane;
1312 if (x > minn[0]) {
1313 if (y < maxn[1])
1314 neigh[9] = k - 1 + nrow - nplane;
1315 neigh[10] = k - 1 - nplane;
1316 if (y > minn[1])
1317 neigh[11] = k - 1 - nrow - nplane;
1318 }
1319 if (y > minn[1]) {
1320 neigh[12] = k - nrow - nplane;
1321 if (x < maxn[0])
1322 neigh[13] = k + 1 - nrow - nplane;
1323 }
1324 }
1325 /*
1326 * Then the same plane
1327 */
1328 if (x < maxn[0] && y < maxn[1])
1329 neigh[14] = k + 1 + nrow;
1330 if (x > minn[0]) {
1331 if (y < maxn[1])
1332 neigh[15] = k - 1 + nrow;
1333 if (y > minn[1])
1334 neigh[16] = k - 1 - nrow;
1335 }
1336 if (y > minn[1] && x < maxn[0])
1337 neigh[17] = k + 1 - nrow - nplane;
1338 /*
1339 * Finally one plane above
1340 */
1341 if (z < maxn[2]) {
1342 if (x < maxn[0]) {
1343 neigh[18] = k + 1 + nplane;
1344 if (y < maxn[1])
1345 neigh[19] = k + 1 + nrow + nplane;
1346 }
1347 if (y < maxn[1])
1348 neigh[20] = k + nrow + nplane;
1349 if (x > minn[0]) {
1350 if (y < maxn[1])
1351 neigh[21] = k - 1 + nrow + nplane;
1352 neigh[22] = k - 1 + nplane;
1353 if (y > minn[1])
1354 neigh[23] = k - 1 - nrow + nplane;
1355 }
1356 if (y > minn[1]) {
1357 neigh[24] = k - nrow + nplane;
1358 if (x < maxn[0])
1359 neigh[25] = k + 1 - nrow + nplane;
1360 }
1361 }
1362 }
1363 }
1364 }
1365 qInfo("%d sources before omitting any.\n",sp->nuse);
1366 /*
1367 * Exclude infeasible points
1368 */
1369 for (k = 0; k < sp->np; k++) {
1370 dist = (sp->rr.row(k).transpose() - cm).norm();
1371 if (dist < exclude || dist > maxdist) {
1372 sp->inuse[k] = 0;
1373 sp->nuse--;
1374 }
1375 }
1376 qInfo("%d sources after omitting infeasible sources.\n",sp->nuse);
1377 {
1378 std::vector<std::unique_ptr<MNESourceSpace>> sp_vec;
1379 sp_vec.push_back(std::move(sp));
1380 if (filter_source_spaces(surf,mindist,FiffCoordTrans(),sp_vec,nullptr) != OK) {
1381 return nullptr;
1382 }
1383 sp = std::move(sp_vec[0]);
1384 }
1385 qInfo("%d sources remaining after excluding the sources outside the surface and less than %6.1f mm inside.\n",sp->nuse,1000*mindist);
1386 /*
1387 * Omit unused vertices from the neighborhoods
1388 */
1389 qInfo("Adjusting the neighborhood info...");
1390 for (k = 0; k < sp->np; k++) {
1391 Eigen::VectorXi& neigh = sp->neighbor_vert[k];
1392 nneigh = sp->nneighbor_vert[k];
1393 if (sp->inuse[k]) {
1394 for (c = 0; c < nneigh; c++)
1395 if (!sp->inuse[neigh[c]])
1396 neigh[c] = -1;
1397 }
1398 else {
1399 for (c = 0; c < nneigh; c++)
1400 neigh[c] = -1;
1401 }
1402 }
1403 qInfo("[done]\n");
1404 /*
1405 * Set up the volume data (needed for creating the interpolation matrix)
1406 */
1407 {
1408 Eigen::Vector3f r0(minn[0]*grid, minn[1]*grid, minn[2]*grid);
1409 Eigen::Vector3f voxel_size(grid, grid, grid);
1410 Eigen::Vector3f x_ras = Eigen::Vector3f::UnitX();
1411 Eigen::Vector3f y_ras = Eigen::Vector3f::UnitY();
1412 Eigen::Vector3f z_ras = Eigen::Vector3f::UnitZ();
1413 int width = (maxn[0]-minn[0]+1);
1414 int height = (maxn[1]-minn[1]+1);
1415 int depth = (maxn[2]-minn[2]+1);
1416
1417 sp->voxel_surf_RAS_t = make_voxel_ras_trans(r0,x_ras,y_ras,z_ras,voxel_size);
1418 if (!sp->voxel_surf_RAS_t || sp->voxel_surf_RAS_t->isEmpty())
1419 return nullptr;
1420
1421 sp->vol_dims[0] = width;
1422 sp->vol_dims[1] = height;
1423 sp->vol_dims[2] = depth;
1424 Eigen::Map<Eigen::Vector3f>(sp->voxel_size) = voxel_size;
1425 }
1426
1427 return sp.release();
1428}
1429
1430//=============================================================================================================
1431
1432int MNESourceSpace::filter_source_spaces(const MNESurface& surf, float limit, const FiffCoordTrans& mri_head_t, std::vector<std::unique_ptr<MNESourceSpace>>& spaces, QTextStream *filtered) /* Provide a list of filtered points here */
1433/*
1434 * Remove all source space points closer to the surface than a given limit
1435 */
1436{
1437 MNESourceSpace* s;
1438 int k,p1,p2;
1439 Eigen::Vector3f r1;
1440 float mindist,dist;
1441 int minnode;
1442 int omit,omit_outside;
1443 double tot_angle;
1444 int nspace = static_cast<int>(spaces.size());
1445
1446 if (spaces[0]->coord_frame == FIFFV_COORD_HEAD && mri_head_t.isEmpty()) {
1447 qCritical("Source spaces are in head coordinates and no coordinate transform was provided!");
1448 return FAIL;
1449 }
1450 /*
1451 * How close are the source points to the surface?
1452 */
1453 qInfo("Source spaces are in ");
1454 if (spaces[0]->coord_frame == FIFFV_COORD_HEAD)
1455 qInfo("head coordinates.\n");
1456 else if (spaces[0]->coord_frame == FIFFV_COORD_MRI)
1457 qInfo("MRI coordinates.\n");
1458 else
1459 qWarning("unknown (%d) coordinates.\n",spaces[0]->coord_frame);
1460 qInfo("Checking that the sources are inside the bounding surface ");
1461 if (limit > 0.0)
1462 qInfo("and at least %6.1f mm away",1000*limit);
1463 qInfo(" (will take a few...)\n");
1464 omit = 0;
1465 omit_outside = 0;
1466 for (k = 0; k < nspace; k++) {
1467 s = spaces[k].get();
1468 for (p1 = 0; p1 < s->np; p1++)
1469 if (s->inuse[p1]) {
1470 r1 = s->rr.row(p1).transpose(); /* Transform the point to MRI coordinates */
1471 if (s->coord_frame == FIFFV_COORD_HEAD)
1472 FiffCoordTrans::apply_inverse_trans(r1.data(),mri_head_t,FIFFV_MOVE);
1473 /*
1474 * Check that the source is inside the inner skull surface
1475 */
1476 tot_angle = surf.sum_solids(r1)/(4*M_PI);
1477 if (std::fabs(tot_angle-1.0) > 1e-5) {
1478 omit_outside++;
1479 s->inuse[p1] = 0;
1480 s->nuse--;
1481 if (filtered)
1482 *filtered << qSetFieldWidth(10) << qSetRealNumberPrecision(3) << Qt::fixed
1483 << 1000*r1[X] << " " << 1000*r1[Y] << " " << 1000*r1[Z] << "\n" << qSetFieldWidth(0);
1484 }
1485 else if (limit > 0.0) {
1486 /*
1487 * Check the distance limit
1488 */
1489 mindist = 1.0;
1490 minnode = 0;
1491 for (p2 = 0; p2 < surf.np; p2++) {
1492 dist = (surf.rr.row(p2).transpose() - r1).norm();
1493 if (dist < mindist) {
1494 mindist = dist;
1495 minnode = p2;
1496 }
1497 }
1498 if (mindist < limit) {
1499 omit++;
1500 s->inuse[p1] = 0;
1501 s->nuse--;
1502 if (filtered)
1503 *filtered << qSetFieldWidth(10) << qSetRealNumberPrecision(3) << Qt::fixed
1504 << 1000*r1[X] << " " << 1000*r1[Y] << " " << 1000*r1[Z] << "\n" << qSetFieldWidth(0);
1505 }
1506 }
1507 }
1508 }
1509 (void)minnode; // squash compiler warning, this is unused
1510 if (omit_outside > 0)
1511 qInfo("%d source space points omitted because they are outside the inner skull surface.\n",
1512 omit_outside);
1513 if (omit > 0)
1514 qInfo("%d source space points omitted because of the %6.1f-mm distance limit.\n",
1515 omit,1000*limit);
1516 qInfo("Thank you for waiting.\n");
1517 return OK;
1518}
1519
1520//=============================================================================================================
1521
1522void MNESourceSpace::filter_source_space(FilterThreadArg *arg)
1523{
1524 FilterThreadArg* a = arg;
1525 int p1,p2;
1526 double tot_angle;
1527 int omit,omit_outside;
1528 Eigen::Vector3f r1;
1529 float mindist,dist;
1530 int minnode;
1531
1532 QSharedPointer<MNESurface> surf = a->surf.toStrongRef();
1533 if (!surf) {
1534 a->stat = FAIL;
1535 return;
1536 }
1537
1538 omit = 0;
1539 omit_outside = 0;
1540
1541 for (p1 = 0; p1 < a->s->np; p1++) {
1542 if (a->s->inuse[p1]) {
1543 r1 = a->s->rr.row(p1).transpose(); /* Transform the point to MRI coordinates */
1544 if (a->s->coord_frame == FIFFV_COORD_HEAD) {
1545 Q_ASSERT(a->mri_head_t);
1546 FiffCoordTrans::apply_inverse_trans(r1.data(),*a->mri_head_t,FIFFV_MOVE);
1547 }
1548 /*
1549 * Check that the source is inside the inner skull surface
1550 */
1551 tot_angle = surf->sum_solids(r1)/(4*M_PI);
1552 if (std::fabs(tot_angle-1.0) > 1e-5) {
1553 omit_outside++;
1554 a->s->inuse[p1] = 0;
1555 a->s->nuse--;
1556 if (a->filtered)
1557 *a->filtered << qSetFieldWidth(10) << qSetRealNumberPrecision(3) << Qt::fixed
1558 << 1000*r1[X] << " " << 1000*r1[Y] << " " << 1000*r1[Z] << "\n" << qSetFieldWidth(0);
1559 }
1560 else if (a->limit > 0.0) {
1561 /*
1562 * Check the distance limit
1563 */
1564 mindist = 1.0;
1565 minnode = 0;
1566 for (p2 = 0; p2 < surf->np; p2++) {
1567 dist = (surf->rr.row(p2).transpose() - r1).norm();
1568 if (dist < mindist) {
1569 mindist = dist;
1570 minnode = p2;
1571 }
1572 }
1573 if (mindist < a->limit) {
1574 omit++;
1575 a->s->inuse[p1] = 0;
1576 a->s->nuse--;
1577 if (a->filtered)
1578 *a->filtered << qSetFieldWidth(10) << qSetRealNumberPrecision(3) << Qt::fixed
1579 << 1000*r1[X] << " " << 1000*r1[Y] << " " << 1000*r1[Z] << "\n" << qSetFieldWidth(0);
1580 }
1581 }
1582 }
1583 }
1584 (void)minnode; // squash compiler warning, set but unused
1585 if (omit_outside > 0)
1586 qInfo("%d source space points omitted because they are outside the inner skull surface.\n",
1587 omit_outside);
1588 if (omit > 0)
1589 qInfo("%d source space points omitted because of the %6.1f-mm distance limit.\n",
1590 omit,1000*a->limit);
1591 a->stat = OK;
1592 return;
1593}
1594
1595//=============================================================================================================
1596
1597int MNESourceSpace::filter_source_spaces(float limit, const QString& bemfile, const FiffCoordTrans& mri_head_t, std::vector<std::unique_ptr<MNESourceSpace>>& spaces, QTextStream *filtered, bool use_threads)
1598/*
1599 * Remove all source space points closer to the surface than a given limit
1600 */
1601{
1602 QSharedPointer<MNESurface> surf;
1603 int k;
1604 int nproc = QThread::idealThreadCount();
1605 int nspace = static_cast<int>(spaces.size());
1606
1607 if (bemfile.isEmpty())
1608 return OK;
1609
1610 {
1611 auto rawSurf = MNESurface::read_bem_surface(bemfile,FIFFV_BEM_SURF_ID_BRAIN,false);
1612 if (!rawSurf) {
1613 qCritical("BEM model does not have the inner skull triangulation!");
1614 return FAIL;
1615 }
1616 surf.reset(rawSurf.release());
1617 }
1618 /*
1619 * How close are the source points to the surface?
1620 */
1621 qInfo("Source spaces are in ");
1622 if (spaces[0]->coord_frame == FIFFV_COORD_HEAD)
1623 qInfo("head coordinates.\n");
1624 else if (spaces[0]->coord_frame == FIFFV_COORD_MRI)
1625 qInfo("MRI coordinates.\n");
1626 else
1627 qWarning("unknown (%d) coordinates.\n",spaces[0]->coord_frame);
1628 qInfo("Checking that the sources are inside the inner skull ");
1629 if (limit > 0.0)
1630 qInfo("and at least %6.1f mm away",1000*limit);
1631 qInfo(" (will take a few...)\n");
1632 if (nproc < 2 || nspace == 1 || !use_threads) {
1633 /*
1634 * This is the conventional calculation
1635 */
1636 for (k = 0; k < nspace; k++) {
1637 auto a_ptr = std::make_unique<FilterThreadArg>();
1638 a_ptr->s = spaces[k].get();
1639 a_ptr->mri_head_t = std::make_unique<FiffCoordTrans>(mri_head_t);
1640 a_ptr->surf = surf;
1641 a_ptr->limit = limit;
1642 a_ptr->filtered = filtered;
1643 filter_source_space(a_ptr.get());
1644 spaces[k]->rearrange_source_space();
1645 }
1646 }
1647 else {
1648 /*
1649 * Calculate all (both) source spaces simultaneously
1650 */
1651 QList<FilterThreadArg*> args;
1652
1653 std::vector<std::unique_ptr<FilterThreadArg>> arg_owners;
1654 for (k = 0; k < nspace; k++) {
1655 auto a_ptr = std::make_unique<FilterThreadArg>();
1656 a_ptr->s = spaces[k].get();
1657 a_ptr->mri_head_t = std::make_unique<FiffCoordTrans>(mri_head_t);
1658 a_ptr->surf = surf;
1659 a_ptr->limit = limit;
1660 a_ptr->filtered = filtered;
1661 args.append(a_ptr.get());
1662 arg_owners.push_back(std::move(a_ptr));
1663 }
1664 /*
1665 * Ready to start the threads & Wait for them to complete
1666 */
1667 QtConcurrent::blockingMap(args, filter_source_space);
1668
1669 for (k = 0; k < nspace; k++) {
1670 spaces[k]->rearrange_source_space();
1671 }
1672 }
1673 qInfo("Thank you for waiting.\n\n");
1674
1675 return OK;
1676}
1677
1678//=============================================================================================================
1679
1680int MNESourceSpace::read_source_spaces(const QString &name, std::vector<std::unique_ptr<MNESourceSpace>>& spaces)
1681/*
1682 * Read source spaces from a FIFF file
1683 */
1684{
1685 QFile file(name);
1686 FiffStream::SPtr stream(new FiffStream(&file));
1687
1688 std::vector<std::unique_ptr<MNESourceSpace>> local_spaces;
1689 std::unique_ptr<MNESourceSpace> new_space;
1690 QList<FiffDirNode::SPtr> sources;
1691 FiffDirNode::SPtr node;
1692 FiffTag::UPtr t_pTag;
1693 int j,k,p,q;
1694 int ntri;
1695
1696 if(!stream->open()) {
1697 stream->close();
1698 return FIFF_FAIL;
1699 }
1700
1701 sources = stream->dirtree()->dir_tree_find(FIFFB_MNE_SOURCE_SPACE);
1702 if (sources.size() == 0) {
1703 qCritical("No source spaces available here");
1704 stream->close();
1705 return FIFF_FAIL;
1706 }
1707 for (j = 0; j < sources.size(); j++) {
1708 new_space = MNESourceSpace::create_source_space(0);
1709 node = sources[j];
1710 /*
1711 * Get the mandatory data first
1712 */
1713 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NPOINTS, t_pTag)) {
1714 stream->close();
1715 return FIFF_FAIL;
1716 }
1717 new_space->np = *t_pTag->toInt();
1718 if (new_space->np == 0) {
1719 qCritical("No points in this source space");
1720 stream->close();
1721 return FIFF_FAIL;
1722 }
1723 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_POINTS, t_pTag)) {
1724 stream->close();
1725 return FIFF_FAIL;
1726 }
1727 MatrixXf tmp_rr = t_pTag->toFloatMatrix().transpose();
1728 new_space->rr = tmp_rr;
1729 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NORMALS, t_pTag)) {
1730 stream->close();
1731 return FIFF_FAIL;
1732 }
1733 MatrixXf tmp_nn = t_pTag->toFloatMatrix().transpose();
1734 new_space->nn = tmp_nn;
1735 if (!node->find_tag(stream, FIFF_MNE_COORD_FRAME, t_pTag)) {
1736 new_space->coord_frame = FIFFV_COORD_MRI;
1737 }
1738 else {
1739 new_space->coord_frame = *t_pTag->toInt();
1740 }
1741 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_ID, t_pTag)) {
1742 new_space->id = *t_pTag->toInt();
1743 }
1744 if (node->find_tag(stream, FIFF_SUBJ_HIS_ID, t_pTag)) {
1745 new_space->subject = t_pTag->toString();
1746 }
1747 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_TYPE, t_pTag)) {
1748 new_space->type = *t_pTag->toInt();
1749 }
1750 ntri = 0;
1751 if (node->find_tag(stream, FIFF_BEM_SURF_NTRI, t_pTag)) {
1752 ntri = *t_pTag->toInt();
1753 }
1754 else if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NTRI, t_pTag)) {
1755 ntri = *t_pTag->toInt();
1756 }
1757 if (ntri > 0) {
1758
1759 if (!node->find_tag(stream, FIFF_BEM_SURF_TRIANGLES, t_pTag)) {
1760 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_TRIANGLES, t_pTag)) {
1761 stream->close();
1762 return FIFF_FAIL;
1763 }
1764 }
1765
1766 MatrixXi tmp_itris = t_pTag->toIntMatrix().transpose();
1767 tmp_itris.array() -= 1;
1768 new_space->itris = tmp_itris;
1769 new_space->ntri = ntri;
1770 }
1771 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NUSE, t_pTag)) {
1772 if (new_space->type == FIFFV_MNE_SPACE_VOLUME) {
1773 /*
1774 * Use all
1775 */
1776 new_space->nuse = new_space->np;
1777 new_space->inuse = Eigen::VectorXi::Ones(new_space->nuse);
1778 new_space->vertno = Eigen::VectorXi::LinSpaced(new_space->nuse, 0, new_space->nuse - 1);
1779 }
1780 else {
1781 /*
1782 * None in use
1783 * NOTE: The consequences of this change have to be evaluated carefully
1784 */
1785 new_space->nuse = 0;
1786 new_space->inuse = Eigen::VectorXi::Zero(new_space->np);
1787 new_space->vertno.resize(0);
1788 }
1789 }
1790 else {
1791 new_space->nuse = *t_pTag->toInt();
1792 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_SELECTION, t_pTag)) {
1793 stream->close();
1794 return FIFF_FAIL;
1795 }
1796
1797 {
1798 Eigen::Map<Eigen::VectorXi> inuseMap(t_pTag->toInt(), new_space->np);
1799 new_space->inuse = inuseMap;
1800 }
1801 if (new_space->nuse > 0) {
1802 new_space->vertno = Eigen::VectorXi::Zero(new_space->nuse);
1803 for (k = 0, p = 0; k < new_space->np; k++) {
1804 if (new_space->inuse[k])
1805 new_space->vertno[p++] = k;
1806 }
1807 }
1808 else {
1809 new_space->vertno.resize(0);
1810 }
1811 /*
1812 * Selection triangulation
1813 */
1814 ntri = 0;
1815 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NUSE_TRI, t_pTag)) {
1816 ntri = *t_pTag->toInt();
1817 }
1818 if (ntri > 0) {
1819
1820 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES, t_pTag)) {
1821 stream->close();
1822 return FIFF_FAIL;
1823 }
1824
1825 MatrixXi tmp_itris = t_pTag->toIntMatrix().transpose();
1826 tmp_itris.array() -= 1;
1827 new_space->use_itris = tmp_itris;
1828 new_space->nuse_tri = ntri;
1829 }
1830 /*
1831 * The patch information becomes relevant here
1832 */
1833 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEAREST, t_pTag)) {
1834 Eigen::Map<Eigen::VectorXi> nearestMap(t_pTag->toInt(), new_space->np);
1835 new_space->nearest.resize(new_space->np);
1836 for (k = 0; k < new_space->np; k++) {
1837 new_space->nearest[k].vert = k;
1838 new_space->nearest[k].nearest = nearestMap[k];
1839 new_space->nearest[k].patch = nullptr;
1840 }
1841
1842 if (!node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEAREST_DIST, t_pTag)) {
1843 stream->close();
1844 return FIFF_FAIL;
1845 }
1846 Eigen::Map<const Eigen::VectorXf> nearestDistMap(t_pTag->toFloat(), new_space->np);
1847 for (k = 0; k < new_space->np; k++) {
1848 new_space->nearest[k].dist = nearestDistMap[k];
1849 }
1850 }
1851 /*
1852 * We may have the distance matrix
1853 */
1854 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_DIST_LIMIT, t_pTag)) {
1855 new_space->dist_limit = *t_pTag->toFloat();
1856 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_DIST, t_pTag)) {
1857 // SparseMatrix<double> tmpSparse = t_pTag->toSparseFloatMatrix();
1858 auto dist_lower = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag);
1859 if (!dist_lower) {
1860 stream->close();
1861 return FIFF_FAIL;
1862 }
1863 auto dist_full = dist_lower->mne_add_upper_triangle_rcs();
1864 if (!dist_full) {
1865 stream->close();
1866 return FIFF_FAIL;
1867 }
1868 new_space->dist = std::move(*dist_full);
1869 }
1870 else
1871 new_space->dist_limit = 0.0;
1872 }
1873 }
1874 /*
1875 * For volume source spaces we might have the neighborhood information
1876 */
1877 if (new_space->type == FIFFV_MNE_SPACE_VOLUME) {
1878 int ntot,nvert,ntot_count,nneigh;
1879
1880 Eigen::VectorXi neighborsVec;
1881 Eigen::VectorXi nneighborsVec;
1882 ntot = nvert = 0;
1883 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NEIGHBORS, t_pTag)) {
1884 ntot = t_pTag->size()/sizeof(fiff_int_t);
1885 neighborsVec = Eigen::Map<Eigen::VectorXi>(t_pTag->toInt(), ntot);
1886 }
1887 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_NNEIGHBORS, t_pTag)) {
1888 nvert = t_pTag->size()/sizeof(fiff_int_t);
1889 nneighborsVec = Eigen::Map<Eigen::VectorXi>(t_pTag->toInt(), nvert);
1890 }
1891 if (neighborsVec.size() > 0 && nneighborsVec.size() > 0) {
1892 if (nvert != new_space->np) {
1893 qCritical("Inconsistent neighborhood data in file.");
1894 stream->close();
1895 return FIFF_FAIL;
1896 }
1897 for (k = 0, ntot_count = 0; k < nvert; k++)
1898 ntot_count += nneighborsVec[k];
1899 if (ntot_count != ntot) {
1900 qCritical("Inconsistent neighborhood data in file.");
1901 stream->close();
1902 return FIFF_FAIL;
1903 }
1904 new_space->nneighbor_vert = Eigen::VectorXi::Zero(nvert);
1905 new_space->neighbor_vert.resize(nvert);
1906 for (k = 0, q = 0; k < nvert; k++) {
1907 new_space->nneighbor_vert[k] = nneigh = nneighborsVec[k];
1908 new_space->neighbor_vert[k] = Eigen::VectorXi(nneigh);
1909 for (p = 0; p < nneigh; p++,q++)
1910 new_space->neighbor_vert[k][p] = neighborsVec[q];
1911 }
1912 }
1913 /*
1914 * There might be a coordinate transformation and dimensions
1915 */
1916 new_space->voxel_surf_RAS_t = FiffCoordTrans(FiffCoordTrans::readTransformFromNode(stream, node, FIFFV_MNE_COORD_MRI_VOXEL, FIFFV_MNE_COORD_SURFACE_RAS));
1917 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS, t_pTag)) {
1918 Eigen::Map<Eigen::Vector3i> volDimsMap(t_pTag->toInt());
1919 Eigen::Map<Eigen::Vector3i>(new_space->vol_dims) = volDimsMap;
1920 }
1921 {
1922 QList<FiffDirNode::SPtr> mris = node->dir_tree_find(FIFFB_MNE_PARENT_MRI_FILE);
1923
1924 if (mris.size() == 0) { /* The old way */
1925 new_space->MRI_surf_RAS_RAS_t = FiffCoordTrans(FiffCoordTrans::readTransformFromNode(stream, node, FIFFV_MNE_COORD_SURFACE_RAS, FIFFV_MNE_COORD_RAS));
1926 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_MRI_FILE, t_pTag)) {
1927 new_space->MRI_volume = t_pTag->toString();
1928 }
1929 if (node->find_tag(stream, FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, t_pTag)) {
1930 new_space->interpolator = std::move(*FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag));
1931 }
1932 }
1933 else {
1934 if (node->find_tag(stream, FIFF_MNE_FILE_NAME, t_pTag)) {
1935 new_space->MRI_volume = t_pTag->toString();
1936 }
1937 new_space->MRI_surf_RAS_RAS_t = FiffCoordTrans(FiffCoordTrans::readTransformFromNode(stream, mris[0], FIFFV_MNE_COORD_SURFACE_RAS, FIFFV_MNE_COORD_RAS));
1938 new_space->MRI_voxel_surf_RAS_t = FiffCoordTrans(FiffCoordTrans::readTransformFromNode(stream, mris[0], FIFFV_MNE_COORD_MRI_VOXEL, FIFFV_MNE_COORD_SURFACE_RAS));
1939
1940 if (mris[0]->find_tag(stream, FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, t_pTag)) {
1941 new_space->interpolator = std::move(*FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag));
1942 }
1943 if (mris[0]->find_tag(stream, FIFF_MRI_WIDTH, t_pTag)) {
1944 new_space->MRI_vol_dims[0] = *t_pTag->toInt();
1945 }
1946 if (mris[0]->find_tag(stream, FIFF_MRI_HEIGHT, t_pTag)) {
1947 new_space->MRI_vol_dims[1] = *t_pTag->toInt();
1948 }
1949 if (mris[0]->find_tag(stream, FIFF_MRI_DEPTH, t_pTag)) {
1950 new_space->MRI_vol_dims[2] = *t_pTag->toInt();
1951 }
1952 }
1953 }
1954 }
1955 new_space->add_triangle_data();
1956 local_spaces.push_back(std::move(new_space));
1957 }
1958 stream->close();
1959
1960 spaces = std::move(local_spaces);
1961
1962 return FIFF_OK;
1963}
1964
1965//=============================================================================================================
1966
1967int MNESourceSpace::transform_source_spaces_to(int coord_frame, const FiffCoordTrans& t, std::vector<std::unique_ptr<MNESourceSpace>>& spaces)
1968/*
1969 * Facilitate the transformation of the source spaces
1970 */
1971{
1972 MNESourceSpace* s;
1973 int k;
1974 int nspace = static_cast<int>(spaces.size());
1975
1976 for (k = 0; k < nspace; k++) {
1977 s = spaces[k].get();
1978 if (s->coord_frame != coord_frame) {
1979 if (!t.isEmpty()) {
1980 if (s->coord_frame == t.from && t.to == coord_frame) {
1981 if (s->transform_source_space(t) != OK)
1982 return FAIL;
1983 }
1984 else if (s->coord_frame == t.to && t.from == coord_frame) {
1985 FiffCoordTrans my_t = t.inverted();
1986 if (s->transform_source_space(my_t) != OK) {
1987 return FAIL;
1988 }
1989 }
1990 else {
1991 qCritical("Could not transform a source space because of transformation incompatibility.");
1992 return FAIL;
1993 }
1994 }
1995 else {
1996 qCritical("Could not transform a source space because of missing coordinate transformation.");
1997 return FAIL;
1998 }
1999 }
2000 }
2001 return OK;
2002}
2003
2004//=============================================================================================================
2005
2006#define LH_LABEL_TAG "-lh.label"
2007#define RH_LABEL_TAG "-rh.label"
2008
2009int MNESourceSpace::restrict_sources_to_labels(std::vector<std::unique_ptr<MNESourceSpace>>& spaces, const QStringList& labels, int nlabel)
2010/*
2011 * Pick only sources within a label
2012 */
2013{
2014 MNESourceSpace* lh = nullptr;
2015 MNESourceSpace* rh = nullptr;
2016 MNESourceSpace* sp;
2017 Eigen::VectorXi lh_inuse;
2018 Eigen::VectorXi rh_inuse;
2019 Eigen::VectorXi sel;
2020 Eigen::VectorXi *inuse = nullptr;
2021 int k,p;
2022 int nspace = static_cast<int>(spaces.size());
2023
2024 if (nlabel == 0)
2025 return OK;
2026
2027 for (k = 0; k < nspace; k++) {
2028 if (spaces[k]->is_left_hemi()) {
2029 lh = spaces[k].get();
2030 lh_inuse = Eigen::VectorXi::Zero(lh->np);
2031 }
2032 else {
2033 rh = spaces[k].get();
2034 rh_inuse = Eigen::VectorXi::Zero(rh->np);
2035 }
2036 }
2037 /*
2038 * Go through each label file
2039 */
2040 for (k = 0; k < nlabel; k++) {
2041 /*
2042 * Which hemi?
2043 */
2044 if (labels[k].contains(LH_LABEL_TAG)){ //strstr(labels[k],LH_LABEL_TAG) != NULL) {
2045 sp = lh;
2046 inuse = &lh_inuse;
2047 }
2048 else if (labels[k].contains(RH_LABEL_TAG)){ //strstr(labels[k],RH_LABEL_TAG) != NULL) {
2049 sp = rh;
2050 inuse = &rh_inuse;
2051 }
2052 else {
2053 qWarning("\tWarning: cannot assign label file %s to a hemisphere.\n",labels[k].toUtf8().constData());
2054 continue;
2055 }
2056 if (sp) {
2057 if (read_label(labels[k],sel) == FAIL)
2058 return FAIL;
2059 for (p = 0; p < sel.size(); p++) {
2060 if (sel[p] >= 0 && sel[p] < sp->np)
2061 (*inuse)[sel[p]] = sp->inuse[sel[p]];
2062 else
2063 qWarning("vertex number out of range in %s (%d vs %d)\n",
2064 labels[k].toUtf8().constData(),sel[p],sp->np);
2065 }
2066 qInfo("Processed label file %s\n",labels[k].toUtf8().constData());
2067 }
2068 }
2069 if (lh) lh->update_inuse(std::move(lh_inuse));
2070 if (rh) rh->update_inuse(std::move(rh_inuse));
2071 return OK;
2072}
2073
2074//=============================================================================================================
2075
2076int MNESourceSpace::read_label(const QString& label, Eigen::VectorXi& sel)
2077/*
2078 * Find the source points within a label
2079 */
2080{
2081 int k,p,nlabel;
2082 char c;
2083 float fdum;
2084 /*
2085 * Read the label file
2086 */
2087 QFile inFile(label);
2088 if (!inFile.open(QIODevice::ReadOnly | QIODevice::Text)) {
2089 qCritical() << label;//err_set_sys_error(label);
2090 sel.resize(0);
2091 return FAIL;
2092 }
2093 inFile.getChar(&c);
2094 if (c !='#') {
2095 qCritical("FsLabel file does not start correctly.");
2096 sel.resize(0);
2097 return FAIL;
2098 }
2099 /*
2100 * Skip the comment line
2101 */
2102 while (inFile.getChar(&c) && c != '\n')
2103 ;
2104 {
2105 QTextStream in(&inFile);
2106 in >> nlabel;
2107 if (in.status() != QTextStream::Ok) {
2108 qCritical("Could not read the number of labelled points.");
2109 sel.resize(0);
2110 return FAIL;
2111 }
2112 sel.resize(nlabel);
2113 for (k = 0; k < nlabel; k++) {
2114 in >> p >> fdum >> fdum >> fdum >> fdum;
2115 if (in.status() != QTextStream::Ok) {
2116 qCritical("Could not read label point # %d",k+1);
2117 sel.resize(0);
2118 return FAIL;
2119 }
2120 sel[k] = p;
2121 }
2122 }
2123
2124 return OK;
2125}
2126
2127//=============================================================================================================
2128
2129int MNESourceSpace::writeVolumeInfo(FiffStream::SPtr& stream, bool selected_only) const
2130{
2131 int ntot,nvert;
2132 int nneigh;
2133 int k,p;
2134
2135 if (type != FIFFV_MNE_SPACE_VOLUME)
2136 return OK;
2137 if (neighbor_vert.empty() || nneighbor_vert.size() == 0)
2138 return OK;
2139
2140 Eigen::VectorXi nneighbors;
2141 Eigen::VectorXi neighbors;
2142
2143 if (selected_only) {
2144 Eigen::VectorXi inuse_map = Eigen::VectorXi::Constant(np, -1);
2145 for (k = 0,p = 0, ntot = 0; k < np; k++) {
2146 if (inuse[k]) {
2147 ntot += nneighbor_vert[k];
2148 inuse_map[k] = p++;
2149 }
2150 }
2151 nneighbors.resize(nuse);
2152 neighbors.resize(ntot);
2153 for (k = 0, nvert = 0, ntot = 0; k < np; k++) {
2154 if (inuse[k]) {
2155 const Eigen::VectorXi& neigh = neighbor_vert[k];
2156 nneigh = nneighbor_vert[k];
2157 nneighbors[nvert++] = nneigh;
2158 for (p = 0; p < nneigh; p++)
2159 neighbors[ntot++] = neigh[p] < 0 ? -1 : inuse_map[neigh[p]];
2160 }
2161 }
2162 }
2163 else {
2164 for (k = 0, ntot = 0; k < np; k++)
2165 ntot += nneighbor_vert[k];
2166 nneighbors.resize(np);
2167 neighbors.resize(ntot);
2168 nvert = np;
2169 for (k = 0, ntot = 0; k < np; k++) {
2170 const Eigen::VectorXi& neigh = neighbor_vert[k];
2171 nneigh = nneighbor_vert[k];
2172 nneighbors[k] = nneigh;
2173 for (p = 0; p < nneigh; p++)
2174 neighbors[ntot++] = neigh[p];
2175 }
2176 }
2177
2178 stream->write_int(FIFF_MNE_SOURCE_SPACE_NNEIGHBORS,nneighbors.data(),nvert);
2179 stream->write_int(FIFF_MNE_SOURCE_SPACE_NEIGHBORS,neighbors.data(),ntot);
2180
2181 if (!selected_only) {
2182 if (voxel_surf_RAS_t && !voxel_surf_RAS_t->isEmpty()) {
2183 stream->write_coord_trans(*voxel_surf_RAS_t);
2184 stream->write_int(FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS,vol_dims,3);
2185 }
2186 if (interpolator && !MRI_volume.isEmpty()) {
2187 stream->start_block(FIFFB_MNE_PARENT_MRI_FILE);
2188 if (MRI_surf_RAS_RAS_t && !MRI_surf_RAS_RAS_t->isEmpty())
2189 stream->write_coord_trans(*MRI_surf_RAS_RAS_t);
2190 if (MRI_voxel_surf_RAS_t && !MRI_voxel_surf_RAS_t->isEmpty())
2191 stream->write_coord_trans(*MRI_voxel_surf_RAS_t);
2192 stream->write_string(FIFF_MNE_FILE_NAME,MRI_volume);
2193 if (interpolator)
2194 stream->write_float_sparse_rcs(FIFF_MNE_SOURCE_SPACE_INTERPOLATOR, interpolator->eigen());
2195 if (MRI_vol_dims[0] > 0 && MRI_vol_dims[1] > 0 && MRI_vol_dims[2] > 0) {
2196 stream->write_int(FIFF_MRI_WIDTH,&MRI_vol_dims[0]);
2197 stream->write_int(FIFF_MRI_HEIGHT,&MRI_vol_dims[1]);
2198 stream->write_int(FIFF_MRI_DEPTH,&MRI_vol_dims[2]);
2199 }
2200 stream->end_block(FIFFB_MNE_PARENT_MRI_FILE);
2201 }
2202 }
2203 else {
2204 if (interpolator && !MRI_volume.isEmpty()) {
2205 stream->write_string(FIFF_MNE_SOURCE_SPACE_MRI_FILE,MRI_volume);
2206 qCritical("Cannot write the interpolator for selection yet");
2207 return FAIL;
2208 }
2209 }
2210 return OK;
2211}
2212
2213//=============================================================================================================
2214
2215int MNESourceSpace::writeToStream(FiffStream::SPtr& stream, bool selected_only) const
2216{
2217 int p, pp;
2218
2219 if (np <= 0) {
2220 qCritical("No points in the source space being saved");
2221 return FIFF_FAIL;
2222 }
2223
2224 stream->start_block(FIFFB_MNE_SOURCE_SPACE);
2225
2226 if (type != FIFFV_MNE_SPACE_UNKNOWN)
2227 stream->write_int(FIFF_MNE_SOURCE_SPACE_TYPE, &type);
2228 if (id != FIFFV_MNE_SURF_UNKNOWN)
2229 stream->write_int(FIFF_MNE_SOURCE_SPACE_ID, &id);
2230 if (!subject.isEmpty() && subject.size() > 0) {
2231 QString subj(subject);
2232 stream->write_string(FIFF_SUBJ_HIS_ID, subj);
2233 }
2234
2235 stream->write_int(FIFF_MNE_COORD_FRAME, &coord_frame);
2236
2237 if (selected_only) {
2238 if (nuse == 0) {
2239 qCritical("No vertices in use. Cannot write active-only vertices from this source space");
2240 return FIFF_FAIL;
2241 }
2242
2243 Eigen::MatrixXf sel(nuse, 3);
2244 stream->write_int(FIFF_MNE_SOURCE_SPACE_NPOINTS, &nuse);
2245
2246 for (p = 0, pp = 0; p < np; p++) {
2247 if (inuse[p]) {
2248 sel.row(pp) = rr.row(p);
2249 pp++;
2250 }
2251 }
2252 stream->write_float_matrix(FIFF_MNE_SOURCE_SPACE_POINTS, sel);
2253
2254 for (p = 0, pp = 0; p < np; p++) {
2255 if (inuse[p]) {
2256 sel.row(pp) = nn.row(p);
2257 pp++;
2258 }
2259 }
2260 stream->write_float_matrix(FIFF_MNE_SOURCE_SPACE_NORMALS, sel);
2261 }
2262 else {
2263 stream->write_int(FIFF_MNE_SOURCE_SPACE_NPOINTS, &np);
2264 stream->write_float_matrix(FIFF_MNE_SOURCE_SPACE_POINTS, Eigen::MatrixXf(rr));
2265 stream->write_float_matrix(FIFF_MNE_SOURCE_SPACE_NORMALS, Eigen::MatrixXf(nn));
2266
2267 if (nuse > 0 && inuse.size() > 0) {
2268 stream->write_int(FIFF_MNE_SOURCE_SPACE_SELECTION, inuse.data(), np);
2269 stream->write_int(FIFF_MNE_SOURCE_SPACE_NUSE, &nuse);
2270 }
2271
2272 if (ntri > 0) {
2273 stream->write_int(FIFF_MNE_SOURCE_SPACE_NTRI, &ntri);
2274 Eigen::MatrixXi file_tris = itris.array() + 1;
2275 stream->write_int_matrix(FIFF_MNE_SOURCE_SPACE_TRIANGLES, file_tris);
2276 }
2277
2278 if (nuse_tri > 0) {
2279 stream->write_int(FIFF_MNE_SOURCE_SPACE_NUSE_TRI, &nuse_tri);
2280 Eigen::MatrixXi file_use_tris = use_itris.array() + 1;
2281 stream->write_int_matrix(FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES, file_use_tris);
2282 }
2283
2284 if (!nearest.empty()) {
2285 Eigen::VectorXi nearest_v(np);
2286 Eigen::VectorXf nearest_dist_v(np);
2287
2288 std::sort(const_cast<std::vector<MNENearest>&>(nearest).begin(),
2289 const_cast<std::vector<MNENearest>&>(nearest).end(),
2290 [](const MNENearest& a, const MNENearest& b) { return a.vert < b.vert; });
2291 for (p = 0; p < np; p++) {
2292 nearest_v[p] = nearest[p].nearest;
2293 nearest_dist_v[p] = nearest[p].dist;
2294 }
2295
2296 stream->write_int(FIFF_MNE_SOURCE_SPACE_NEAREST, nearest_v.data(), np);
2297 stream->write_float(FIFF_MNE_SOURCE_SPACE_NEAREST_DIST, nearest_dist_v.data(), np);
2298 }
2299
2300 if (!dist.is_empty()) {
2301 auto m = dist.pickLowerTriangleRcs();
2302 if (!m)
2303 return FIFF_FAIL;
2304 stream->write_float_sparse_rcs(FIFF_MNE_SOURCE_SPACE_DIST, m->eigen());
2305 stream->write_float(FIFF_MNE_SOURCE_SPACE_DIST_LIMIT, &dist_limit);
2306 }
2307 }
2308
2309 if (writeVolumeInfo(stream, selected_only) != OK)
2310 return FIFF_FAIL;
2311
2312 stream->end_block(FIFFB_MNE_SOURCE_SPACE);
2313 return FIFF_OK;
2314}
2315
2316//=============================================================================================================
2317
2322static Eigen::MatrixX3f generateIcoVertices(int grade)
2323{
2324 // Base icosahedron vertices
2325 const float t = (1.0f + std::sqrt(5.0f)) / 2.0f;
2326 std::vector<Eigen::Vector3f> verts = {
2327 {-1, t, 0}, { 1, t, 0}, {-1, -t, 0}, { 1, -t, 0},
2328 { 0, -1, t}, { 0, 1, t}, { 0, -1, -t}, { 0, 1, -t},
2329 { t, 0, -1}, { t, 0, 1}, {-t, 0, -1}, {-t, 0, 1}
2330 };
2331 // Normalize to unit sphere
2332 for (auto& v : verts)
2333 v.normalize();
2334
2335 // Base icosahedron faces
2336 std::vector<std::array<int,3>> faces = {
2337 {0,11,5}, {0,5,1}, {0,1,7}, {0,7,10}, {0,10,11},
2338 {1,5,9}, {5,11,4}, {11,10,2}, {10,7,6}, {7,1,8},
2339 {3,9,4}, {3,4,2}, {3,2,6}, {3,6,8}, {3,8,9},
2340 {4,9,5}, {2,4,11}, {6,2,10}, {8,6,7}, {9,8,1}
2341 };
2342
2343 // Subdivide
2344 for (int g = 0; g < grade; ++g) {
2345 std::map<std::pair<int,int>, int> midpointCache;
2346 std::vector<std::array<int,3>> newFaces;
2347
2348 auto getMidpoint = [&](int i1, int i2) -> int {
2349 auto key = std::make_pair(std::min(i1, i2), std::max(i1, i2));
2350 auto it = midpointCache.find(key);
2351 if (it != midpointCache.end())
2352 return it->second;
2353 Eigen::Vector3f mid = (verts[i1] + verts[i2]).normalized();
2354 int idx = static_cast<int>(verts.size());
2355 verts.push_back(mid);
2356 midpointCache[key] = idx;
2357 return idx;
2358 };
2359
2360 for (const auto& f : faces) {
2361 int a = getMidpoint(f[0], f[1]);
2362 int b = getMidpoint(f[1], f[2]);
2363 int c = getMidpoint(f[2], f[0]);
2364 newFaces.push_back({f[0], a, c});
2365 newFaces.push_back({f[1], b, a});
2366 newFaces.push_back({f[2], c, b});
2367 newFaces.push_back({a, b, c});
2368 }
2369 faces = newFaces;
2370 }
2371
2372 Eigen::MatrixX3f result(static_cast<int>(verts.size()), 3);
2373 for (int i = 0; i < static_cast<int>(verts.size()); ++i)
2374 result.row(i) = verts[i];
2375 return result;
2376}
2377
2378//=============================================================================================================
2379
2380MNEHemisphere MNESourceSpace::icoDownsample(const MNEHemisphere& hemi, int icoGrade)
2381{
2382 MNEHemisphere result(hemi);
2383
2384 if (hemi.np <= 0 || hemi.rr.rows() == 0) {
2385 qWarning("MNESourceSpace::icoDownsample - Hemisphere has no vertices.");
2386 return result;
2387 }
2388
2389 // Generate icosahedral surface at the requested grade
2390 Eigen::MatrixX3f icoVerts = generateIcoVertices(icoGrade);
2391
2392 // Project hemisphere vertices to unit sphere for matching
2393 Eigen::MatrixX3f hemiNorm(hemi.np, 3);
2394 for (int i = 0; i < hemi.np; ++i) {
2395 Eigen::Vector3f v = hemi.rr.row(i);
2396 float len = v.norm();
2397 if (len > 0.0f)
2398 hemiNorm.row(i) = (v / len).transpose();
2399 else
2400 hemiNorm.row(i) = v.transpose();
2401 }
2402
2403 // Clear inuse
2404 result.inuse = Eigen::VectorXi::Zero(result.np);
2405
2406 // For each icosahedral vertex, find the nearest hemisphere vertex
2407 for (int i = 0; i < icoVerts.rows(); ++i) {
2408 Eigen::Vector3f icoV = icoVerts.row(i);
2409 float bestDist = std::numeric_limits<float>::max();
2410 int bestIdx = -1;
2411 for (int j = 0; j < hemi.np; ++j) {
2412 float d = (hemiNorm.row(j).transpose() - icoV).squaredNorm();
2413 if (d < bestDist) {
2414 bestDist = d;
2415 bestIdx = j;
2416 }
2417 }
2418 if (bestIdx >= 0)
2419 result.inuse[bestIdx] = 1;
2420 }
2421
2422 // Recount nuse and rebuild vertno
2423 result.nuse = result.inuse.sum();
2424 result.vertno.resize(result.nuse);
2425 int k = 0;
2426 for (int i = 0; i < result.np; ++i) {
2427 if (result.inuse[i])
2428 result.vertno[k++] = i;
2429 }
2430
2431 return result;
2432}
2433
M_PI
#define M_PI
Definition coregsettingsview.cpp:35
FAIL
constexpr int FAIL
Definition compute_fwd.cpp:69
Y
constexpr int Y
Definition compute_fwd.cpp:73
Z
constexpr int Z
Definition compute_fwd.cpp:74
OK
constexpr int OK
Definition compute_fwd.cpp:70
X
constexpr int X
Definition compute_fwd.cpp:72
fiff_coord_trans.h
FiffCoordTrans class declaration.
FIFF_MRI_DEPTH
#define FIFF_MRI_DEPTH
Definition fiff_file.h:657
FIFF_BEM_SURF_TRIANGLES
#define FIFF_BEM_SURF_TRIANGLES
Definition fiff_file.h:736
FIFF_MRI_HEIGHT
#define FIFF_MRI_HEIGHT
Definition fiff_file.h:655
FIFF_BEM_SURF_NTRI
#define FIFF_BEM_SURF_NTRI
Definition fiff_file.h:734
FIFFV_BEM_SURF_ID_BRAIN
#define FIFFV_BEM_SURF_ID_BRAIN
Definition fiff_file.h:749
FIFF_SUBJ_HIS_ID
#define FIFF_SUBJ_HIS_ID
Definition fiff_file.h:573
FIFF_MRI_WIDTH
#define FIFF_MRI_WIDTH
Definition fiff_file.h:653
fiff_stream.h
FiffStream class declaration.
fiff_tag.h
FiffTag class declaration, which provides fiff tag I/O and processing methods.
fiff_sparse_matrix.h
FiffSparseMatrix class declaration.
fiff_constants.h
Fiff constants.
FIFFV_MNE_SURF_RIGHT_HEMI
#define FIFFV_MNE_SURF_RIGHT_HEMI
Definition fiff_constants.h:296
FIFF_MNE_SOURCE_SPACE_ID
#define FIFF_MNE_SOURCE_SPACE_ID
Definition fiff_constants.h:359
FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
#define FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
Definition fiff_constants.h:372
FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS
#define FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS
Definition fiff_constants.h:367
FIFF_MNE_COORD_FRAME
#define FIFF_MNE_COORD_FRAME
Definition fiff_constants.h:342
FIFF_OK
#define FIFF_OK
Definition fiff_constants.h:60
FIFF_MNE_SOURCE_SPACE_TYPE
#define FIFF_MNE_SOURCE_SPACE_TYPE
Definition fiff_constants.h:360
FIFFV_MNE_SURF_UNKNOWN
#define FIFFV_MNE_SURF_UNKNOWN
Definition fiff_constants.h:294
FIFF_MNE_SOURCE_SPACE_SELECTION
#define FIFF_MNE_SOURCE_SPACE_SELECTION
Definition fiff_constants.h:355
FIFFV_MNE_SURF_LEFT_HEMI
#define FIFFV_MNE_SURF_LEFT_HEMI
Definition fiff_constants.h:295
FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
#define FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
Definition fiff_constants.h:365
FIFF_MNE_SOURCE_SPACE_NUSE_TRI
#define FIFF_MNE_SOURCE_SPACE_NUSE_TRI
Definition fiff_constants.h:364
FIFF_MNE_SOURCE_SPACE_INTERPOLATOR
#define FIFF_MNE_SOURCE_SPACE_INTERPOLATOR
Definition fiff_constants.h:368
FIFF_MNE_SOURCE_SPACE_NORMALS
#define FIFF_MNE_SOURCE_SPACE_NORMALS
Definition fiff_constants.h:353
FIFFV_MNE_COORD_MRI_VOXEL
#define FIFFV_MNE_COORD_MRI_VOXEL
Definition fiff_constants.h:493
FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
#define FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
Definition fiff_constants.h:358
FIFF_MNE_SOURCE_SPACE_DIST
#define FIFF_MNE_SOURCE_SPACE_DIST
Definition fiff_constants.h:371
FIFF_MNE_SOURCE_SPACE_POINTS
#define FIFF_MNE_SOURCE_SPACE_POINTS
Definition fiff_constants.h:352
FIFF_FAIL
#define FIFF_FAIL
Definition fiff_constants.h:59
FIFFV_MNE_SPACE_SURFACE
#define FIFFV_MNE_SPACE_SURFACE
Definition fiff_constants.h:470
FIFF_MNE_SOURCE_SPACE_NTRI
#define FIFF_MNE_SOURCE_SPACE_NTRI
Definition fiff_constants.h:362
FIFF_MNE_SOURCE_SPACE_MRI_FILE
#define FIFF_MNE_SOURCE_SPACE_MRI_FILE
Definition fiff_constants.h:369
FIFFB_MNE_SOURCE_SPACE
#define FIFFB_MNE_SOURCE_SPACE
Definition fiff_constants.h:315
FIFF_MNE_SOURCE_SPACE_NPOINTS
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
Definition fiff_constants.h:354
FIFFV_MNE_SPACE_VOLUME
#define FIFFV_MNE_SPACE_VOLUME
Definition fiff_constants.h:471
FIFFV_NO_MOVE
#define FIFFV_NO_MOVE
Definition fiff_constants.h:231
FIFFV_COORD_HEAD
#define FIFFV_COORD_HEAD
Definition fiff_constants.h:216
FIFFV_COORD_MRI
#define FIFFV_COORD_MRI
Definition fiff_constants.h:217
FIFF_MNE_SOURCE_SPACE_TRIANGLES
#define FIFF_MNE_SOURCE_SPACE_TRIANGLES
Definition fiff_constants.h:363
FIFFV_MNE_SPACE_UNKNOWN
#define FIFFV_MNE_SPACE_UNKNOWN
Definition fiff_constants.h:469
FIFF_MNE_SOURCE_SPACE_NUSE
#define FIFF_MNE_SOURCE_SPACE_NUSE
Definition fiff_constants.h:356
FIFFV_MOVE
#define FIFFV_MOVE
Definition fiff_constants.h:230
FIFFV_MNE_COORD_RAS
#define FIFFV_MNE_COORD_RAS
Definition fiff_constants.h:494
FIFF_MNE_FILE_NAME
#define FIFF_MNE_FILE_NAME
Definition fiff_constants.h:347
FIFF_MNE_SOURCE_SPACE_NEAREST
#define FIFF_MNE_SOURCE_SPACE_NEAREST
Definition fiff_constants.h:357
FIFFB_MNE_PARENT_MRI_FILE
#define FIFFB_MNE_PARENT_MRI_FILE
Definition fiff_constants.h:317
fiff_byte_swap.h
Byte-swap utility functions for FIFF binary I/O.
MNE_SOURCE_SPACE_VOLUME
#define MNE_SOURCE_SPACE_VOLUME
Definition mne_types.h:104
mne_source_space.h
MNESourceSpace class declaration.
CURVATURE_FILE_MAGIC_NUMBER
constexpr int CURVATURE_FILE_MAGIC_NUMBER
Definition mne_source_space.cpp:79
TAG_USEREALRAS
constexpr int TAG_USEREALRAS
Definition mne_source_space.cpp:84
NNEIGHBORS
constexpr int NNEIGHBORS
Definition mne_source_space.cpp:77
TAG_OLD_MGH_XFORM
constexpr int TAG_OLD_MGH_XFORM
Definition mne_source_space.cpp:81
TAG_OLD_USEREALRAS
constexpr int TAG_OLD_USEREALRAS
Definition mne_source_space.cpp:83
LH_LABEL_TAG
#define LH_LABEL_TAG
Definition mne_source_space.cpp:2006
TAG_OLD_COLORTABLE
constexpr int TAG_OLD_COLORTABLE
Definition mne_source_space.cpp:82
EVEN
#define EVEN(n)
RH_LABEL_TAG
#define RH_LABEL_TAG
Definition mne_source_space.cpp:2007
mne_patch_info.h
MNE Patch Information (MNEPatchInfo) class declaration.
mne_nearest.h
MNENearest class declaration.
QUAD_FILE_MAGIC_NUMBER
#define QUAD_FILE_MAGIC_NUMBER
Definition mne_surface_or_volume.h:81
NEW_QUAD_FILE_MAGIC_NUMBER
#define NEW_QUAD_FILE_MAGIC_NUMBER
Definition mne_surface_or_volume.h:80
FIFFV_MNE_COORD_SURFACE_RAS
#define FIFFV_MNE_COORD_SURFACE_RAS
Definition mne_surface_or_volume.h:77
FIFF_MNE_SOURCE_SPACE_NEIGHBORS
#define FIFF_MNE_SOURCE_SPACE_NEIGHBORS
Definition mne_surface_or_volume.h:75
FIFF_MNE_SOURCE_SPACE_NNEIGHBORS
#define FIFF_MNE_SOURCE_SPACE_NNEIGHBORS
Definition mne_surface_or_volume.h:74
TAG_OLD_SURF_GEOM
#define TAG_OLD_SURF_GEOM
Definition mne_surface_or_volume.h:83
TRIANGLE_FILE_MAGIC_NUMBER
#define TRIANGLE_FILE_MAGIC_NUMBER
Definition mne_surface_or_volume.h:79
mne_mgh_tag_group.h
MNEMghTagGroup class declaration.
mne_hemisphere.h
MNEHemisphere class declaration.
filter_thread_arg.h
Filter Thread Argument (FilterThreadArg) class declaration.
mne_surface.h
MNESurface class declaration.
MNELIB
Core MNE data structures (source spaces, source estimates, hemispheres).
Definition connectivitysettings.h:63
FIFFLIB
FIFF file I/O and data structures (raw, epochs, evoked, covariance, forward).
Definition connectivitysettings.h:71
FIFFLIB::swap_float
float swap_float(float source)
Definition fiff_byte_swap.h:131
FIFFLIB::swap_int
qint32 swap_int(qint32 source)
Definition fiff_byte_swap.h:69
FIFFLIB::fiff_int_t
qint32 fiff_int_t
Definition fiff_types.h:89
FIFFLIB::swap_long
qint64 swap_long(qint64 source)
Definition fiff_byte_swap.h:97
FIFFLIB::swap_short
qint16 swap_short(qint16 source)
Definition fiff_byte_swap.h:55
FIFFLIB::FiffCoordTrans
Coordinate transformation description.
Definition fiff_coord_trans.h:81
FIFFLIB::FiffCoordTrans::readTransformFromNode
static FiffCoordTrans readTransformFromNode(FiffStream::SPtr &stream, const FiffDirNode::SPtr &node, int from, int to)
Definition fiff_coord_trans.cpp:592
FIFFLIB::FiffCoordTrans::inverted
FiffCoordTrans inverted() const
Definition fiff_coord_trans.cpp:365
FIFFLIB::FiffCoordTrans::apply_inverse_trans
Eigen::MatrixX3f apply_inverse_trans(const Eigen::MatrixX3f &rr, bool do_move=true) const
FIFFLIB::FiffCoordTrans::apply_trans
Eigen::MatrixX3f apply_trans(const Eigen::MatrixX3f &rr, bool do_move=true) const
FIFFLIB::FiffDirNode::SPtr
QSharedPointer< FiffDirNode > SPtr
Definition fiff_dir_node.h:77
FIFFLIB::FiffSparseMatrix
FIFF sparse matrix storage backed by Eigen.
Definition fiff_sparse_matrix.h:86
FIFFLIB::FiffSparseMatrix::fiff_get_float_sparse_matrix
static FiffSparseMatrix::UPtr fiff_get_float_sparse_matrix(const FIFFLIB::FiffTag::UPtr &tag)
Definition fiff_sparse_matrix.cpp:168
FIFFLIB::FiffStream
FIFF File I/O routines.
Definition fiff_stream.h:107
FIFFLIB::FiffStream::SPtr
QSharedPointer< FiffStream > SPtr
Definition fiff_stream.h:109
FIFFLIB::FiffTag::UPtr
std::unique_ptr< FiffTag > UPtr
Definition fiff_tag.h:158
MNELIB::FilterThreadArg
Thread-local arguments for parallel raw data filtering (channel range, filter kernel,...
Definition filter_thread_arg.h:81
MNELIB::FilterThreadArg::surf
QWeakPointer< MNESurface > surf
Definition filter_thread_arg.h:103
MNELIB::FilterThreadArg::mri_head_t
std::unique_ptr< FIFFLIB::FiffCoordTrans > mri_head_t
Definition filter_thread_arg.h:102
MNELIB::MNEHemisphere
Hemisphere provides geometry information.
Definition mne_hemisphere.h:83
MNELIB::MNEMghTagGroup
Collection of MNEMghTag entries from a FreeSurfer MGH/MGZ file footer.
Definition mne_mgh_tag_group.h:73
MNELIB::MNEMghTagGroup::tags
std::vector< std::unique_ptr< MNEMghTag > > tags
Definition mne_mgh_tag_group.h:112
MNELIB::MNENearest
This is used in the patch definitions.
Definition mne_nearest.h:78
MNELIB::MNENearest::patch
MNEPatchInfo * patch
Definition mne_nearest.h:100
MNELIB::MNEPatchInfo
Patch information for a single source space point including vertex members and area.
Definition mne_patch_info.h:78
MNELIB::MNESourceSpace::clone
virtual MNESourceSpace::SPtr clone() const
Definition mne_source_space.cpp:836
MNELIB::MNESourceSpace::load_surface
static std::unique_ptr< MNESourceSpace > load_surface(const QString &surf_file, const QString &curv_file)
Definition mne_source_space.cpp:1107
MNELIB::MNESourceSpace::find_source_space_hemi
qint32 find_source_space_hemi() const
Definition mne_source_space.cpp:893
MNELIB::MNESourceSpace::SPtr
std::shared_ptr< MNESourceSpace > SPtr
Definition mne_source_space.h:81
MNELIB::MNESourceSpace::restrict_sources_to_labels
static int restrict_sources_to_labels(std::vector< std::unique_ptr< MNESourceSpace > > &spaces, const QStringList &labels, int nlabel)
Definition mne_source_space.cpp:2009
MNELIB::MNESourceSpace::filter_source_spaces
static int filter_source_spaces(const MNESurface &surf, float limit, const FIFFLIB::FiffCoordTrans &mri_head_t, std::vector< std::unique_ptr< MNESourceSpace > > &spaces, QTextStream *filtered)
Definition mne_source_space.cpp:1432
MNELIB::MNESourceSpace::read_source_spaces
static int read_source_spaces(const QString &name, std::vector< std::unique_ptr< MNESourceSpace > > &spaces)
Definition mne_source_space.cpp:1680
MNELIB::MNESourceSpace::create_source_space
static std::unique_ptr< MNESourceSpace > create_source_space(int np)
Definition mne_source_space.cpp:1060
MNELIB::MNESourceSpace::load_surface_geom
static std::unique_ptr< MNESourceSpace > load_surface_geom(const QString &surf_file, const QString &curv_file, bool add_geometry, bool check_too_many_neighbors)
Definition mne_source_space.cpp:1115
MNELIB::MNESourceSpace::icoDownsample
static MNEHemisphere icoDownsample(const MNEHemisphere &hemi, int icoGrade)
Definition mne_source_space.cpp:2380
MNELIB::MNESourceSpace::filter_source_space
static void filter_source_space(FilterThreadArg *arg)
Definition mne_source_space.cpp:1522
MNELIB::MNESourceSpace::writeToStream
int writeToStream(FIFFLIB::FiffStream::SPtr &stream, bool selected_only) const
Definition mne_source_space.cpp:2215
MNELIB::MNESourceSpace::make_volume_source_space
static MNESourceSpace * make_volume_source_space(const MNESurface &surf, float grid, float exclude, float mindist)
Definition mne_source_space.cpp:1195
MNELIB::MNESourceSpace::transform_source_space
int transform_source_space(const FIFFLIB::FiffCoordTrans &t)
Definition mne_source_space.cpp:932
MNELIB::MNESourceSpace::read_label
static int read_label(const QString &label, Eigen::VectorXi &sel)
Definition mne_source_space.cpp:2076
MNELIB::MNESourceSpace::transform_source_spaces_to
static int transform_source_spaces_to(int coord_frame, const FIFFLIB::FiffCoordTrans &t, std::vector< std::unique_ptr< MNESourceSpace > > &spaces)
Definition mne_source_space.cpp:1967
MNELIB::MNESourceSpace::update_inuse
void update_inuse(Eigen::VectorXi new_inuse)
Definition mne_source_space.cpp:904
MNELIB::MNESurface
This defines a surface.
Definition mne_surface.h:79
MNELIB::MNESurface::sum_solids
double sum_solids(const Eigen::Vector3f &from) const
Definition mne_surface.cpp:92
MNELIB::MNESurface::read_bem_surface
static std::unique_ptr< MNESurface > read_bem_surface(const QString &name, int which, bool add_geometry)
Definition mne_surface.cpp:376
MNELIB::MNESurfaceOrVolume::neighbor_tri
std::vector< Eigen::VectorXi > neighbor_tri
Definition mne_surface_or_volume.h:321
MNELIB::MNESurfaceOrVolume::neighbor_vert
std::vector< Eigen::VectorXi > neighbor_vert
Definition mne_surface_or_volume.h:304
MNELIB::MNESurfaceOrVolume::MRI_surf_RAS_RAS_t
std::optional< FIFFLIB::FiffCoordTrans > MRI_surf_RAS_RAS_t
Definition mne_surface_or_volume.h:341
MNELIB::MNESurfaceOrVolume::add_geometry_info
int add_geometry_info(bool do_normals, bool check_too_many_neighbors)
Definition mne_surface_or_volume.cpp:369
MNELIB::MNESurfaceOrVolume::dist
FIFFLIB::FiffSparseMatrix dist
Definition mne_surface_or_volume.h:327
MNELIB::MNESurfaceOrVolume::nearest
std::vector< MNENearest > nearest
Definition mne_surface_or_volume.h:324
MNELIB::MNESurfaceOrVolume::interpolator
std::optional< FIFFLIB::FiffSparseMatrix > interpolator
Definition mne_surface_or_volume.h:338
MNELIB::MNESurfaceOrVolume::TrianglesT
Eigen::Matrix< int, Eigen::Dynamic, 3, Eigen::RowMajor > TrianglesT
Definition mne_surface_or_volume.h:130
MNELIB::MNESurfaceOrVolume::vol_geom
std::optional< MNEVolGeom > vol_geom
Definition mne_surface_or_volume.h:283
MNELIB::MNESurfaceOrVolume::patches
std::vector< std::optional< MNEPatchInfo > > patches
Definition mne_surface_or_volume.h:325
MNELIB::MNESurfaceOrVolume::MRI_voxel_surf_RAS_t
std::optional< FIFFLIB::FiffCoordTrans > MRI_voxel_surf_RAS_t
Definition mne_surface_or_volume.h:340
MNELIB::MNESurfaceOrVolume::tris
std::vector< MNETriangle > tris
Definition mne_surface_or_volume.h:313
MNELIB::MNESurfaceOrVolume::mgh_tags
std::optional< MNEMghTagGroup > mgh_tags
Definition mne_surface_or_volume.h:284
MNELIB::MNESurfaceOrVolume::voxel_surf_RAS_t
std::optional< FIFFLIB::FiffCoordTrans > voxel_surf_RAS_t
Definition mne_surface_or_volume.h:335
MNELIB::MNESurfaceOrVolume::PointsT
Eigen::Matrix< float, Eigen::Dynamic, 3, Eigen::RowMajor > PointsT
Definition mne_surface_or_volume.h:128
MNELIB::MNEVolGeom
MRI data volume geometry information like FreeSurfer keeps it.
Definition mne_vol_geom.h:77
Generated on for MNE-CPP by doxygen 1.16.1