mne_bem.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "mne_bem.h"
 
#include <math/warp.h>
#include <fs/fs_label.h>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QFile>
 
#include <stdexcept>
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace UTILSLIB;
using namespace FSLIB;
using namespace MNELIB;
using namespace FIFFLIB;
using namespace Eigen;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
MNEBem::MNEBem()
{
}
 
//=============================================================================================================
 
MNEBem::MNEBem(QIODevice& p_IODevice) //const MNEBem &p_MNEBem
//: m_qListBemSurface()
{
    FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
 
    if (!MNEBem::readFromStream(t_pStream, true, *this)) {
        t_pStream->close();
        throw std::runtime_error("Could not read the bem surfaces\n");
        //ToDo error(me,'Could not read the bem surfaces (%s)',mne_omit_first_line(lasterr));
        //        return false;
    }
 
    //    bool testStream =t_pStream->device()->isOpen();
}
 
//=============================================================================================================
 
MNEBem::~MNEBem()
{
}
 
//=============================================================================================================
 
void MNEBem::clear()
{
    m_qListBemSurface.clear();
}
 
//=============================================================================================================
 
bool MNEBem::readFromStream(FiffStream::SPtr& p_pStream, bool add_geom, MNEBem& p_Bem)
{
    //
    //   Open the file, create directory
    //
    bool open_here = false;
    QFile t_file; //ToDo TCPSocket;
 
    if (!p_pStream->device()->isOpen()) {
        QString t_sFileName = p_pStream->streamName();
 
        t_file.setFileName(t_sFileName);
        p_pStream = FiffStream::SPtr(new FiffStream(&t_file));
        if (!p_pStream->open()) {
            return false;
        }
        open_here = true;
        //        if(t_pDir)
        //            delete t_pDir;
    }
 
    //
    //   Find all BEM surfaces
    //
 
    QList<FiffDirNode::SPtr> bem = p_pStream->dirtree()->dir_tree_find(FIFFB_BEM);
    if (bem.isEmpty()) {
        qCritical() << "No BEM block found!";
        if (open_here) {
            p_pStream->close();
        }
        return false;
    }
 
    QList<FiffDirNode::SPtr> bemsurf = p_pStream->dirtree()->dir_tree_find(FIFFB_BEM_SURF);
    if (bemsurf.isEmpty()) {
        qCritical() << "No BEM surfaces found!";
        if (open_here) {
            p_pStream->close();
        }
        return false;
    }
 
    QList<MNEBemSurface> surfaces;
    for (int k = 0; k < bemsurf.size(); ++k) {
        MNEBemSurface p_BemSurface;
        qInfo("\tReading a BEM surface...");
        if (!MNEBem::readBemSurface(p_pStream, bemsurf[k], p_BemSurface)) {
            if (open_here) {
                p_pStream->close();
            }
            return false;
        }
        p_BemSurface.addTriangleData();
        if (add_geom || p_BemSurface.nn.rows() == 0) {
            p_BemSurface.addVertexNormals();
        }
        qInfo("\t[done]\n");
 
        surfaces.append(p_BemSurface);
    }
    p_Bem.m_qListBemSurface.append(surfaces);
 
    qInfo("\t%lld bem surfaces read\n", static_cast<long long>(bemsurf.size()));
 
    if (open_here) {
        p_pStream->close();
    }
    return true;
}
 
//=============================================================================================================
 
bool MNEBem::readBemSurface(FiffStream::SPtr& p_pStream, const FiffDirNode::SPtr& p_Tree, MNEBemSurface& p_BemSurface)
{
    p_BemSurface.clear();
 
    FiffTag::UPtr t_pTag;
 
    //=====================================================================
    if (!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_ID, t_pTag)) {
        p_BemSurface.id = FIFFV_BEM_SURF_ID_UNKNOWN;
    } else {
        p_BemSurface.id = *t_pTag->toInt();
    }
 
    //    qDebug() << "Read BemSurface ID; type:" << t_pTag->getType() << "value:" << *t_pTag->toInt();
 
    //=====================================================================
    if (!p_Tree->find_tag(p_pStream, FIFF_BEM_SIGMA, t_pTag)) {
        p_BemSurface.sigma = 1.0;
    } else {
        p_BemSurface.sigma = *t_pTag->toFloat();
    }
 
    //    qDebug() <<
 
    //=====================================================================
    if (!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NNODE, t_pTag)) {
        qWarning() << "np not found!";
        return false;
    } else {
        p_BemSurface.np = *t_pTag->toInt();
    }
 
    //    qDebug() <<
 
    //=====================================================================
    if (!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NTRI, t_pTag)) {
        qWarning() << "ntri not found!";
        return false;
    } else {
        p_BemSurface.ntri = *t_pTag->toInt();
    }
 
    //    qDebug() <<
 
    //=====================================================================
    if (!p_Tree->find_tag(p_pStream, FIFF_MNE_COORD_FRAME, t_pTag)) {
        qWarning() << "FIFF_MNE_COORD_FRAME not found, trying FIFF_BEM_COORD_FRAME.";
        if (!p_Tree->find_tag(p_pStream, FIFF_BEM_COORD_FRAME, t_pTag)) {
            qWarning() << "Coordinate frame information not found.";
            return false;
        } else {
            p_BemSurface.coord_frame = *t_pTag->toInt();
        }
    } else {
        p_BemSurface.coord_frame = *t_pTag->toInt();
    }
 
    //    qDebug() <<
 
    //=====================================================================
    //
    //   Vertices, normals, and triangles
    //
    //=====================================================================
    if (!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NODES, t_pTag)) {
        qWarning() << "Vertex data not found.";
        return false;
    }
 
    p_BemSurface.rr = t_pTag->toFloatMatrix().transpose();
    qint32 rows_rr = p_BemSurface.rr.rows();
 
    if (rows_rr != p_BemSurface.np) {
        qWarning() << "Vertex information is incorrect.";
        return false;
    }
 
    //    qDebug() << "Surf Nodes; type:" << t_pTag->getType();
 
    //=====================================================================
    // Stored normals are optional (mne-python, MNE-C); readFromStream computes missing ones.
    if (p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_NORMALS, t_pTag) || p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_NORMALS, t_pTag)) {
        p_BemSurface.nn = t_pTag->toFloatMatrix().transpose();
        if (p_BemSurface.nn.rows() != p_BemSurface.np) {
            qWarning() << "Vertex normal information is incorrect.";
            return false;
        }
    } else {
        p_BemSurface.nn.resize(0, 3);
    }
 
    //    qDebug() << "Bem Vertex Normals; type:" << t_pTag->getType();
 
    //=====================================================================
    if (p_BemSurface.ntri > 0) {
        if (!p_Tree->find_tag(p_pStream, FIFF_BEM_SURF_TRIANGLES, t_pTag)) {
            if (!p_Tree->find_tag(p_pStream, FIFF_MNE_SOURCE_SPACE_TRIANGLES, t_pTag)) {
                qWarning() << "Triangulation not found.";
                return false;
            } else {
                p_BemSurface.itris = t_pTag->toIntMatrix().transpose();
                p_BemSurface.itris -= MatrixXi::Constant(p_BemSurface.itris.rows(), 3, 1); //0 based indizes
            }
        } else {
            p_BemSurface.itris = t_pTag->toIntMatrix().transpose();
            p_BemSurface.itris -= MatrixXi::Constant(p_BemSurface.itris.rows(), 3, 1); //0 based indizes
        }
 
        if (p_BemSurface.itris.rows() != p_BemSurface.ntri) {
            qWarning() << "Triangulation information is incorrect.";
            return false;
        }
    } else {
        p_BemSurface.itris.resize(0, 3);
    }
 
    return true;
}
 
//=============================================================================================================
 
bool MNEBem::write(QIODevice& p_IODevice)
{
    FiffStream::SPtr t_pStream = FiffStream::start_file(p_IODevice);
    if (!t_pStream) {
        return false;
    }
    qInfo("Write BEM surface in %s...\n", t_pStream->streamName().toUtf8().constData());
    this->writeToStream(t_pStream.data());
    t_pStream->end_file();
    return true;
}
 
//=============================================================================================================
 
void MNEBem::writeToStream(FiffStream* p_pStream)
{
    p_pStream->start_block(FIFFB_BEM);
    for (qint32 h = 0; h < m_qListBemSurface.size(); ++h) {
        qInfo("\tWrite a bem surface... ");
        p_pStream->start_block(FIFFB_BEM_SURF);
        m_qListBemSurface[h].writeToStream(p_pStream);
        p_pStream->end_block(FIFFB_BEM_SURF);
        qInfo("[done]\n");
    }
    qInfo("\t%lld bem surfaces written\n", static_cast<long long>(m_qListBemSurface.size()));
    p_pStream->end_block(FIFFB_BEM);
}
 
//=============================================================================================================
 
const MNEBemSurface& MNEBem::operator[](qint32 idx) const
{
    //
    // Falling back to surface 0 is no help when there are no surfaces at all:
    // indexing an empty QList asserts in a debug build and is undefined
    // otherwise. Hand back a default surface instead, which callers can spot
    // through its empty vertex list.
    //
    if (m_qListBemSurface.isEmpty()) {
        qWarning("Warning: No BEM surfaces available! Returning an empty surface.");
 
        static const MNEBemSurface defaultSurface;
 
        return defaultSurface;
    }
 
    if (idx < 0 || idx >= m_qListBemSurface.length()) {
        qWarning("Warning: Required surface doesn't exist! Returning surface '0'.");
        idx = 0;
    }
    return m_qListBemSurface[idx];
}
 
//=============================================================================================================
 
MNEBemSurface& MNEBem::operator[](qint32 idx)
{
    //
    // See the const overload: an empty list has no surface 0 to fall back to,
    // so keep one default surface around rather than indexing out of bounds.
    //
    if (m_qListBemSurface.isEmpty()) {
        qWarning("Warning: No BEM surfaces available! Returning an empty surface.");
 
        static MNEBemSurface defaultSurface;
        defaultSurface = MNEBemSurface();
 
        return defaultSurface;
    }
 
    if (idx < 0 || idx >= m_qListBemSurface.length()) {
        qWarning("Warning: Required surface doesn't exist! Returning surface '0'.");
        idx = 0;
    }
    return m_qListBemSurface[idx];
}
 
//=============================================================================================================
 
MNEBem& MNEBem::operator<<(const MNEBemSurface& surf)
{
    this->m_qListBemSurface.append(surf);
    return *this;
}
 
//=============================================================================================================
 
MNEBem& MNEBem::operator<<(const MNEBemSurface* surf)
{
    this->m_qListBemSurface.append(*surf);
    return *this;
}
 
//=============================================================================================================
 
void MNEBem::warp(const MatrixXf& sLm, const MatrixXf& dLm)
{
    Warp help;
    QList<MatrixXf> vertList;
    for (int i = 0; i < this->m_qListBemSurface.size(); i++) {
        vertList.append(this->m_qListBemSurface[i].rr);
    }
 
    help.calculate(sLm, dLm, vertList);
 
    for (int i = 0; i < this->m_qListBemSurface.size(); i++) {
        this->m_qListBemSurface[i].rr = vertList.at(i);
    }
    return;
}
 
//=============================================================================================================
 
void MNEBem::transform(const FiffCoordTrans& trans)
{
    MatrixX3f vert;
    for (int i = 0; i < this->m_qListBemSurface.size(); i++) {
        vert = this->m_qListBemSurface[i].rr;
        vert = trans.apply_trans(vert);
        this->m_qListBemSurface[i].rr = vert;
    }
    return;
}
 
//=============================================================================================================
 
void MNEBem::invtransform(const FiffCoordTrans& trans)
{
    MatrixX3f vert;
    for (int i = 0; i < this->m_qListBemSurface.size(); i++) {
        vert = this->m_qListBemSurface[i].rr;
        vert = trans.apply_inverse_trans(vert);
        this->m_qListBemSurface[i].rr = vert;
    }
    return;
}