compute_fwd.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "compute_fwd.h"
#include <mne/mne_forward_solution.h>
 
#include <fiff/fiff_stream.h>
#include <fiff/fiff_info.h>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QCoreApplication>
#include <QFile>
#include <QTextStream>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace FWDLIB;
using namespace MNELIB;
using namespace FIFFLIB;
using namespace Eigen;
 
//=============================================================================================================
// CONSTANTS
//=============================================================================================================
 
constexpr int FAIL = -1;
constexpr int OK = 0;
 
constexpr int X = 0;
constexpr int Y = 1;
constexpr int Z = 2;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
ComputeFwd::ComputeFwd(std::shared_ptr<ComputeFwdSettings> pSettings)
: m_pSettings(pSettings)
, m_meg_forward(new FiffNamedMatrix)
, m_meg_forward_grad(new FiffNamedMatrix)
, m_eeg_forward(new FiffNamedMatrix)
, m_eeg_forward_grad(new FiffNamedMatrix)
{
    initFwd();
}
 
//=============================================================================================================
 
ComputeFwd::~ComputeFwd()
{
}
 
//=============================================================================================================
 
//=============================================================================================================
 
void ComputeFwd::initFwd()
{
    m_bInitialized = false;
    m_spaces.clear();
    m_iNSource = 0;
 
    m_mri_head_t = FiffCoordTrans();
    m_meg_head_t = FiffCoordTrans();
 
    m_listMegChs = QList<FiffChInfo>();
    m_listEegChs = QList<FiffChInfo>();
    m_listCompChs = QList<FiffChInfo>();
 
    int iNMeg = 0;
    int iNEeg = 0;
    int iNComp = 0;
 
    m_templates.reset();
    m_megcoils.reset();
    m_compcoils.reset();
    m_eegels.reset();
    m_eegModels.reset();
    m_iNChan = 0;
 
    int k;
    m_mri_id = FiffId();
    m_meas_id.clear();
 
    QFile filteredFile;
    QTextStream* filteredStream = nullptr;
 
    m_eegModel.reset();
    m_bemModel.reset();
 
    // Report the setup
    qInfo("Source space                 : %s", m_pSettings->srcname.toUtf8().constData());
    if (!(m_pSettings->transname.isEmpty()) || !(m_pSettings->mriname.isEmpty())) {
        qInfo("MRI -> head transform source : %s", !(m_pSettings->mriname.isEmpty()) ? m_pSettings->mriname.toUtf8().constData() : m_pSettings->transname.toUtf8().constData());
    } else {
        qInfo("MRI and head coordinates are assumed to be identical.");
    }
    qInfo("Measurement data             : %s", m_pSettings->measname.toUtf8().constData());
    if (!m_pSettings->bemname.isEmpty()) {
        qInfo("BEM model                    : %s", m_pSettings->bemname.toUtf8().constData());
    } else {
        qInfo("Sphere model                 : origin at (% 7.2f % 7.2f % 7.2f) mm",
              1000.0f * m_pSettings->r0[X], 1000.0f * m_pSettings->r0[Y], 1000.0f * m_pSettings->r0[Z]);
        if (m_pSettings->include_eeg) {
            if (m_pSettings->eeg_model_file.isEmpty()) {
                qWarning("EEG model file not specified; using default.");
            }
            m_eegModels.reset(FwdEegSphereModelSet::fwd_load_eeg_sphere_models(m_pSettings->eeg_model_file, m_eegModels.release()));
            m_eegModels->fwd_list_eeg_sphere_models();
 
            if (m_pSettings->eeg_model_name.isEmpty()) {
                m_pSettings->eeg_model_name = QString("Default");
            }
            m_eegModel.reset(m_eegModels->fwd_select_eeg_sphere_model(m_pSettings->eeg_model_name));
            if (!m_eegModel) {
                return;
            }
 
            if (!m_eegModel->fwd_setup_eeg_sphere_model(m_pSettings->eeg_sphere_rad, m_pSettings->use_equiv_eeg, 3)) {
                return;
            }
 
            qInfo("Using EEG sphere model \"%s\" with scalp radius %7.1f mm",
                  m_pSettings->eeg_model_name.toUtf8().constData(), 1000 * m_pSettings->eeg_sphere_rad);
            qInfo("%s the electrode locations to scalp", m_pSettings->scale_eeg_pos ? "Scale" : "Do not scale");
 
            m_eegModel->scale_pos = m_pSettings->scale_eeg_pos;
            m_eegModel->r0 = m_pSettings->r0;
        }
    }
    qInfo("%s field computations", m_pSettings->accurate ? "Accurate" : "Standard");
    qInfo("Do computations in %s coordinates.", FiffCoordTrans::frame_name(m_pSettings->coord_frame).toUtf8().constData());
    qInfo("%s source orientations", m_pSettings->fixed_ori ? "Fixed" : "Free");
    if (m_pSettings->compute_grad) {
        qInfo("Compute derivatives with respect to source location coordinates");
    }
    qInfo("Destination for the solution : %s", m_pSettings->solname.toUtf8().constData());
    if (m_pSettings->do_all) {
        qInfo("Calculate solution for all source locations.");
    }
    if (m_pSettings->nlabel > 0) {
        qInfo("Source space will be restricted to sources in %d labels", m_pSettings->nlabel);
    }
 
    // Read the source locations
    qInfo("Reading %s...", m_pSettings->srcname.toUtf8().constData());
    if (MNESourceSpace::read_source_spaces(m_pSettings->srcname, m_spaces) != OK) {
        return;
    }
    for (k = 0, m_iNSource = 0; k < static_cast<int>(m_spaces.size()); k++) {
        if (m_pSettings->do_all) {
            m_spaces[k]->enable_all_sources();
        }
        m_iNSource += m_spaces[k]->nuse;
    }
    if (m_iNSource == 0) {
        qCritical("No sources are active in these source spaces. --all option should be used.");
        return;
    }
    qInfo("Read %d source spaces a total of %d active source locations", static_cast<int>(m_spaces.size()), m_iNSource);
    if (MNESourceSpace::restrict_sources_to_labels(m_spaces, m_pSettings->labels, m_pSettings->nlabel) == FAIL) {
        return;
    }
 
    // Read the MRI -> head coordinate transformation
    if (!m_pSettings->mriname.isEmpty()) {
        m_mri_head_t = FiffCoordTrans::readMriTransform(m_pSettings->mriname);
        if (m_mri_head_t.isEmpty()) {
            return;
        }
        {
            QFile mriFile(m_pSettings->mriname);
            FiffStream::SPtr mriStream(new FiffStream(&mriFile));
            if (mriStream->open()) {
                m_mri_id.version = mriStream->id().version;
                m_mri_id.machid[0] = mriStream->id().machid[0];
                m_mri_id.machid[1] = mriStream->id().machid[1];
                m_mri_id.time = mriStream->id().time;
                mriStream->close();
            } else {
                mriStream->close();
                m_mri_id = FiffId();
            }
        }
        if (m_mri_id.isEmpty()) {
            qCritical("Couldn't read MRI file id (How come?)");
            return;
        }
    } else if (!m_pSettings->transname.isEmpty()) {
        FiffCoordTrans t = FiffCoordTrans::readFShead2mriTransform(m_pSettings->transname.toUtf8().data());
        if (t.isEmpty()) {
            return;
        }
        m_mri_head_t = t.inverted();
    } else {
        m_mri_head_t = FiffCoordTrans::identity(FIFFV_COORD_MRI, FIFFV_COORD_HEAD);
    }
    m_mri_head_t.print();
 
    // Read the channel information and the MEG device -> head coordinate transformation
 
    if (!m_pSettings->pFiffInfo) {
        QFile measname(m_pSettings->measname);
        FIFFLIB::FiffDirNode::SPtr DirNode;
        FiffStream::SPtr pStream(new FiffStream(&measname));
        FIFFLIB::FiffInfo fiffInfo;
        if (!pStream->open()) {
            qCritical() << "Could not open Stream.";
            return;
        }
 
        if (!pStream->read_meas_info(pStream->dirtree(), fiffInfo, DirNode)) {
            qCritical() << "Could not find the channel information.";
            return;
        }
        pStream->close();
        m_pInfoBase = QSharedPointer<FIFFLIB::FiffInfo>(new FiffInfo(fiffInfo));
    } else {
        m_pInfoBase = m_pSettings->pFiffInfo;
    }
    if (!m_pInfoBase) {
        qCritical("ComputeFwd::initFwd(): no FiffInfo");
        return;
    }
    m_pInfoBase->mne_read_meg_comp_eeg_ch_info(m_listMegChs,
                                               iNMeg,
                                               m_listCompChs,
                                               iNComp,
                                               m_listEegChs,
                                               iNEeg,
                                               m_meg_head_t,
                                               m_meas_id);
    if (!m_pSettings->meg_head_t.isEmpty()) {
        m_meg_head_t = m_pSettings->meg_head_t;
    }
    if (m_meg_head_t.isEmpty()) {
        qCritical("MEG -> head coordinate transformation not found.");
        return;
    }
 
    m_iNChan = iNMeg + iNEeg;
 
    if (iNMeg > 0) {
        qInfo("Read %3d MEG channels from %s", iNMeg, m_pSettings->measname.toUtf8().constData());
    }
    if (iNComp > 0) {
        qInfo("Read %3d MEG compensation channels from %s", iNComp, m_pSettings->measname.toUtf8().constData());
    }
    if (iNEeg > 0) {
        qInfo("Read %3d EEG channels from %s", iNEeg, m_pSettings->measname.toUtf8().constData());
    }
    if (!m_pSettings->include_meg) {
        qInfo("MEG not requested. MEG channels omitted.");
        m_listMegChs.clear();
        m_listCompChs.clear();
        iNMeg = 0;
        iNComp = 0;
    } else
        m_meg_head_t.print();
    if (!m_pSettings->include_eeg) {
        qInfo("EEG not requested. EEG channels omitted.");
        m_listEegChs.clear();
        iNEeg = 0;
    } else {
        if (!FiffChInfo::checkEegLocations(m_listEegChs, iNEeg)) {
            return;
        }
    }
 
    // Create coil descriptions with transformation to head or MRI frame
 
    if (m_pSettings->include_meg) {
        m_qPath = QString(QCoreApplication::applicationDirPath() + "/../resources/general/coilDefinitions/coil_def.dat");
        if (!QCoreApplication::startingUp()) {
            m_qPath = QCoreApplication::applicationDirPath() + QString("/../resources/general/coilDefinitions/coil_def.dat");
        } else if (!QFile::exists(m_qPath)) {
            m_qPath = "../resources/general/coilDefinitions/coil_def.dat";
        }
 
        m_templates = FwdCoilSet::read_coil_defs(m_qPath);
        if (!m_templates) {
            return;
        }
 
        // Compensation data
 
        m_compData = MNECTFCompDataSet::read(m_pSettings->measname);
        if (!m_compData) {
            return;
        }
        if (m_compData->ncomp > 0) {
            qInfo("%d compensation data sets in %s", m_compData->ncomp, m_pSettings->measname.toUtf8().constData());
        } else {
            m_listCompChs.clear();
            iNComp = 0;
 
            m_compData.reset();
        }
    }
    if (m_pSettings->coord_frame == FIFFV_COORD_MRI) {
        FiffCoordTrans head_mri_t = m_mri_head_t.inverted();
        FiffCoordTrans meg_mri_t = FiffCoordTrans::combine(FIFFV_COORD_DEVICE, FIFFV_COORD_MRI, m_meg_head_t, head_mri_t);
        if (meg_mri_t.isEmpty()) {
            return;
        }
        m_megcoils = m_templates->create_meg_coils(m_listMegChs,
                                                   iNMeg,
                                                   m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
                                                   meg_mri_t);
        if (!m_megcoils) {
            return;
        }
        if (iNComp > 0) {
            m_compcoils = m_templates->create_meg_coils(m_listCompChs,
                                                        iNComp,
                                                        FWD_COIL_ACCURACY_NORMAL,
                                                        meg_mri_t);
            if (!m_compcoils) {
                return;
            }
        }
        m_eegels = FwdCoilSet::create_eeg_els(m_listEegChs,
                                              iNEeg,
                                              head_mri_t);
        if (!m_eegels) {
            return;
        }
 
        qInfo("MRI coordinate coil definitions created.");
    } else {
        m_megcoils = m_templates->create_meg_coils(m_listMegChs,
                                                   iNMeg,
                                                   m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
                                                   m_meg_head_t);
        if (!m_megcoils) {
            return;
        }
 
        if (iNComp > 0) {
            m_compcoils = m_templates->create_meg_coils(m_listCompChs,
                                                        iNComp,
                                                        FWD_COIL_ACCURACY_NORMAL, m_meg_head_t);
            if (!m_compcoils) {
                return;
            }
        }
        m_eegels = FwdCoilSet::create_eeg_els(m_listEegChs,
                                              iNEeg);
        if (!m_eegels) {
            return;
        }
        qInfo("Head coordinate coil definitions created.");
    }
 
    // Transform the source spaces into the appropriate coordinates
    {
        if (MNESourceSpace::transform_source_spaces_to(m_pSettings->coord_frame, m_mri_head_t, m_spaces) != OK) {
            return;
        }
    }
    qInfo("Source spaces are now in %s coordinates.", FiffCoordTrans::frame_name(m_pSettings->coord_frame).toUtf8().constData());
 
    // Prepare the BEM model if necessary
 
    if (!m_pSettings->bemname.isEmpty()) {
        QString bemsolname = FwdBemModel::fwd_bem_make_bem_sol_name(m_pSettings->bemname);
        m_pSettings->bemname = bemsolname;
 
        qInfo("Setting up the BEM model using %s...", m_pSettings->bemname.toUtf8().constData());
        qInfo("Loading surfaces...");
        m_bemModel = FwdBemModel::fwd_bem_load_three_layer_surfaces(m_pSettings->bemname);
 
        if (m_bemModel) {
            qInfo("Three-layer model surfaces loaded.");
        } else {
            m_bemModel = FwdBemModel::fwd_bem_load_homog_surface(m_pSettings->bemname);
            if (!m_bemModel) {
                return;
            }
            qInfo("Homogeneous model surface loaded.");
        }
        if (iNEeg > 0 && m_bemModel->nsurf == 1) {
            qCritical("Cannot use a homogeneous model in EEG calculations.");
            return;
        }
        qInfo("Loading the solution matrix...");
        if (m_bemModel->fwd_bem_load_recompute_solution(m_pSettings->bemname.toUtf8().data(), FWD_BEM_UNKNOWN, false) == FAIL) {
            return;
        }
        if (m_pSettings->coord_frame == FIFFV_COORD_HEAD) {
            qInfo("Employing the head->MRI coordinate transform with the BEM model.");
            if (m_bemModel->fwd_bem_set_head_mri_t(m_mri_head_t) == FAIL) {
                return;
            }
        }
        qInfo("BEM model %s is now set up", m_bemModel->sol_name.toUtf8().constData());
    } else {
        qInfo("Using the sphere model.");
        m_bemModel = std::make_unique<FwdBemModel>(); // no surfaces selects the sphere branch
    }
 
    // Try to circumvent numerical problems by excluding points too close or outside the inner skull surface
 
    if (m_pSettings->filter_spaces) {
        if (!m_pSettings->mindistoutname.isEmpty()) {
            filteredFile.setFileName(m_pSettings->mindistoutname);
            if (!filteredFile.open(QIODevice::WriteOnly | QIODevice::Text)) {
                qCritical() << m_pSettings->mindistoutname;
                return;
            }
            filteredStream = new QTextStream(&filteredFile);
            qInfo("Omitted source space points will be output to : %s", m_pSettings->mindistoutname.toUtf8().constData());
        }
        MNESourceSpace::filter_source_spaces(m_pSettings->mindist,
                                             m_pSettings->bemname,
                                             m_mri_head_t,
                                             m_spaces,
                                             filteredStream, m_pSettings->use_threads);
        delete filteredStream;
        filteredStream = nullptr;
    }
    m_bInitialized = true;
}
 
//=============================================================================================================
 
void ComputeFwd::populateMetadata(MNEForwardSolution& fwd)
{
    fwd.coord_frame = m_pSettings->coord_frame;
    fwd.source_ori = m_pSettings->fixed_ori ? FIFFV_MNE_FIXED_ORI : FIFFV_MNE_FREE_ORI;
    fwd.surf_ori = false;
    fwd.mri_filename = m_pSettings->mriname;
 
    int nmeg = m_megcoils ? m_megcoils->ncoil() : 0;
    int neeg = m_eegels ? m_eegels->ncoil() : 0;
    fwd.nchan = nmeg + neeg;
 
    fwd.mri_head_t = m_mri_head_t;
    fwd.mri_id = m_mri_id;
 
    // Source spaces: store in MRI frame (FIFF convention), keep m_spaces in computation frame
    {
        if (MNESourceSpace::transform_source_spaces_to(FIFFV_COORD_MRI, m_mri_head_t, m_spaces) != OK) {
            return;
        }
        fwd.src.clear();
        fwd.nsource = 0;
        for (int i = 0; i < static_cast<int>(m_spaces.size()); ++i) {
            fwd.src.append(*m_spaces[i]);
            fwd.nsource += m_spaces[i]->nuse;
        }
        if (MNESourceSpace::transform_source_spaces_to(m_pSettings->coord_frame, m_mri_head_t, m_spaces) != OK) {
            return;
        }
    }
 
    // Source locations and orientations in the computation frame, as MNEForwardSolution::read sets them.
    const int nOri = m_pSettings->fixed_ori ? 1 : 3;
    fwd.source_rr = MatrixX3f(fwd.nsource, 3);
    fwd.source_nn = MatrixX3f(nOri * fwd.nsource, 3);
    for (int i = 0, p = 0; i < static_cast<int>(m_spaces.size()); ++i) {
        const MNESourceSpace& s = *m_spaces[i];
        for (int k = 0; k < s.nuse; ++k, ++p) {
            fwd.source_rr.row(p) = s.rr.row(s.vertno(k));
            if (nOri == 1) {
                fwd.source_nn.row(p) = s.nn.row(s.vertno(k));
            } else {
                fwd.source_nn.middleRows(3 * p, 3) = Matrix3f::Identity();
            }
        }
    }
 
    // Measurement provenance
    fwd.info.filename = m_pSettings->measname;
    fwd.info.meas_id = m_meas_id;
    fwd.info.dev_head_t = m_meg_head_t;
    fwd.info.nchan = fwd.nchan;
 
    fwd.info.chs.clear();
    fwd.info.ch_names.clear();
    for (int i = 0; i < nmeg; ++i) {
        fwd.info.chs.append(m_listMegChs[i]);
        fwd.info.ch_names.append(m_listMegChs[i].ch_name);
    }
    for (int i = 0; i < neeg; ++i) {
        fwd.info.chs.append(m_listEegChs[i]);
        fwd.info.ch_names.append(m_listEegChs[i].ch_name);
    }
 
    // Bad channels
    if (!m_pSettings->measname.isEmpty()) {
        QFile fileBad(m_pSettings->measname);
        FiffStream::SPtr t_pStreamBads(new FiffStream(&fileBad));
        if (t_pStreamBads->open()) {
            fwd.info.bads = t_pStreamBads->read_bad_channels(t_pStreamBads->dirtree());
            t_pStreamBads->close();
        }
    }
}
 
//=============================================================================================================
 
std::unique_ptr<MNEForwardSolution> ComputeFwd::calculateFwd()
{
    if (!m_bInitialized) {
        qCritical("ComputeFwd::calculateFwd - the forward computation could not be set up.");
        return nullptr;
    }
    auto fwdSolution = std::make_unique<MNEForwardSolution>();
    populateMetadata(*fwdSolution);
    int iNMeg = 0;
    int iNEeg = 0;
 
    if (m_megcoils) {
        iNMeg = m_megcoils->ncoil();
    }
    if (m_eegels) {
        iNEeg = m_eegels->ncoil();
    }
    if (!m_bemModel || m_bemModel->nsurf == 0) {
        m_pSettings->use_threads = false;
    }
 
    // check if source spaces are still in computation frame
    if (m_spaces[0]->coord_frame != m_pSettings->coord_frame) {
        if (MNESourceSpace::transform_source_spaces_to(m_pSettings->coord_frame, m_mri_head_t, m_spaces) != OK) {
            return nullptr;
        }
    }
 
    // Do the actual computation
    if (iNMeg > 0) {
        if ((m_bemModel->compute_forward_meg(m_spaces,
                                             m_megcoils.get(),
                                             m_compcoils.get(),
                                             m_compData.get(),
                                             m_pSettings->fixed_ori,
                                             m_pSettings->r0,
                                             m_pSettings->use_threads,
                                             *m_meg_forward.data(),
                                             *m_meg_forward_grad.data(),
                                             m_pSettings->compute_grad)) == FAIL) {
            return nullptr;
        }
    }
    if (iNEeg > 0) {
        if ((m_bemModel->compute_forward_eeg(m_spaces,
                                             m_eegels.get(),
                                             m_pSettings->fixed_ori,
                                             m_eegModel.get(),
                                             m_pSettings->use_threads,
                                             *m_eeg_forward.data(),
                                             *m_eeg_forward_grad.data(),
                                             m_pSettings->compute_grad)) == FAIL) {
            return nullptr;
        }
    }
 
    // Assemble combined sol
    if (iNMeg > 0 && iNEeg > 0) {
        if (m_meg_forward->data.cols() != m_eeg_forward->data.cols()) {
            qWarning() << "The MEG and EEG forward solutions do not match";
            return nullptr;
        }
        fwdSolution->sol->clear();
        fwdSolution->sol->nrow = m_meg_forward->nrow + m_eeg_forward->nrow;
        fwdSolution->sol->ncol = m_meg_forward->ncol;
        fwdSolution->sol->data = MatrixXd(fwdSolution->sol->nrow, fwdSolution->sol->ncol);
        fwdSolution->sol->data.block(0, 0, m_meg_forward->nrow, m_meg_forward->ncol) = m_meg_forward->data;
        fwdSolution->sol->data.block(m_meg_forward->nrow, 0, m_eeg_forward->nrow, m_eeg_forward->ncol) = m_eeg_forward->data;
        fwdSolution->sol->row_names = m_meg_forward->row_names;
        fwdSolution->sol->row_names.append(m_eeg_forward->row_names);
        fwdSolution->sol->col_names = m_meg_forward->col_names;
    } else if (iNMeg > 0) {
        fwdSolution->sol = m_meg_forward;
    } else {
        fwdSolution->sol = m_eeg_forward;
    }
 
    if (m_pSettings->compute_grad) {
        if (iNMeg > 0 && iNEeg > 0) {
            if (m_meg_forward_grad->data.cols() != m_eeg_forward_grad->data.cols()) {
                qWarning() << "The MEG and EEG forward solutions do not match";
                return nullptr;
            }
            fwdSolution->sol_grad->clear();
            fwdSolution->sol_grad->nrow = m_meg_forward_grad->nrow + m_eeg_forward_grad->nrow;
            fwdSolution->sol_grad->ncol = m_meg_forward_grad->ncol;
            fwdSolution->sol_grad->data = MatrixXd(fwdSolution->sol_grad->nrow, fwdSolution->sol_grad->ncol);
            fwdSolution->sol_grad->data.block(0, 0, m_meg_forward_grad->nrow, m_meg_forward_grad->ncol) = m_meg_forward_grad->data;
            fwdSolution->sol_grad->data.block(m_meg_forward_grad->nrow, 0, m_eeg_forward_grad->nrow, m_eeg_forward_grad->ncol) = m_eeg_forward_grad->data;
            fwdSolution->sol_grad->row_names = m_meg_forward_grad->row_names;
            fwdSolution->sol_grad->row_names.append(m_eeg_forward_grad->row_names);
            fwdSolution->sol_grad->col_names = m_meg_forward_grad->col_names;
        } else if (iNMeg > 0) {
            fwdSolution->sol_grad = m_meg_forward_grad;
        } else {
            fwdSolution->sol_grad = m_eeg_forward_grad;
        }
    }
 
    return fwdSolution;
}
 
//=============================================================================================================
 
bool ComputeFwd::updateHeadPos(const FiffCoordTrans& transDevHead, MNEForwardSolution& fwd)
{
    if (!m_bInitialized)
        return false;
    int iNMeg = 0;
    if (m_megcoils) {
        iNMeg = m_megcoils->ncoil();
    }
 
    int iNComp = 0;
    if (m_compcoils) {
        iNComp = m_compcoils->ncoil();
    }
 
    // create new coilset with updated head position
    if (m_pSettings->coord_frame == FIFFV_COORD_MRI) {
        FiffCoordTrans head_mri_t = m_mri_head_t.inverted();
        FiffCoordTrans meg_mri_t = FiffCoordTrans::combine(FIFFV_COORD_DEVICE, FIFFV_COORD_MRI, transDevHead, head_mri_t);
        if (meg_mri_t.isEmpty()) {
            return false;
        }
        m_megcoils = m_templates->create_meg_coils(m_listMegChs,
                                                   iNMeg,
                                                   m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
                                                   meg_mri_t);
        if (!m_megcoils) {
            return false;
        }
        if (iNComp > 0) {
            m_compcoils = m_templates->create_meg_coils(m_listCompChs,
                                                        iNComp,
                                                        FWD_COIL_ACCURACY_NORMAL,
                                                        meg_mri_t);
            if (!m_compcoils) {
                return false;
            }
        }
    } else {
        m_megcoils = m_templates->create_meg_coils(m_listMegChs,
                                                   iNMeg,
                                                   m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
                                                   transDevHead);
        if (!m_megcoils) {
            return false;
        }
 
        if (iNComp > 0) {
            m_compcoils = m_templates->create_meg_coils(m_listCompChs,
                                                        iNComp,
                                                        FWD_COIL_ACCURACY_NORMAL, transDevHead);
            if (!m_compcoils) {
                return false;
            }
        }
    }
 
    // check if source spaces are still in computation frame
    if (m_spaces[0]->coord_frame != m_pSettings->coord_frame) {
        if (MNESourceSpace::transform_source_spaces_to(m_pSettings->coord_frame, m_mri_head_t, m_spaces) != OK) {
            return false;
        }
    }
 
    // recompute meg forward
    if ((m_bemModel->compute_forward_meg(m_spaces,
                                         m_megcoils.get(),
                                         m_compcoils.get(),
                                         m_compData.get(),
                                         m_pSettings->fixed_ori,
                                         m_pSettings->r0,
                                         m_pSettings->use_threads,
                                         *m_meg_forward.data(),
                                         *m_meg_forward_grad.data(),
                                         m_pSettings->compute_grad)) == FAIL) {
        return false;
    }
 
    // Update transformation matrix and info
    m_meg_head_t = transDevHead;
    fwd.info.dev_head_t = transDevHead;
 
    // update solution
    fwd.sol->data.block(0, 0, m_meg_forward->nrow, m_meg_forward->ncol) = m_meg_forward->data;
    if (m_pSettings->compute_grad) {
        fwd.sol_grad->data.block(0, 0, m_meg_forward_grad->nrow, m_meg_forward_grad->ncol) = m_meg_forward_grad->data;
    }
 
    return true;
}