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