fwd_comp_data.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "fwd_comp_data.h"
 
#include <mne/mne_ctf_comp_data_set.h>
#include <fiff/fiff_types.h>
 
namespace
{
constexpr int FAIL = -1;
constexpr int OK = 0;
}
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace Eigen;
using namespace FIFFLIB;
using namespace MNELIB;
using namespace FWDLIB;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
FwdCompData::FwdCompData()
: field(nullptr)
, vec_field(nullptr)
, field_grad(nullptr)
, client(nullptr)
{
}
 
//=============================================================================================================
 
FwdCompData::~FwdCompData() = default;
 
//=============================================================================================================
 
int FwdCompData::fwd_comp_field(const Eigen::Vector3f& rd, const Eigen::Vector3f& Q, FwdCoilSet& coils, Eigen::Ref<Eigen::VectorXf> res, void* client)
{
    FwdCompData* comp = static_cast<FwdCompData*>(client);
 
    if (!comp->field) {
        qWarning("Field computation function is missing in fwd_comp_field");
        return FAIL;
    }
    /*
       * First compute the field in the primary set of coils
       */
    if (comp->field(rd, Q, coils, res, comp->client) == FAIL)
        return FAIL;
    /*
       * Compensation needed?
       */
    if (!comp->comp_coils || comp->comp_coils->ncoil() <= 0 || !comp->set || !comp->set->current)
        return OK;
    /*
       * Workspace needed?
       */
    if (comp->work.size() == 0)
        comp->work.resize(comp->comp_coils->ncoil());
    /*
       * Compute the field in the compensation coils
       */
    if (comp->field(rd, Q, *comp->comp_coils, comp->work, comp->client) == FAIL)
        return FAIL;
    /*
       * Compute the compensated field
       */
    return comp->set->apply(true, res, comp->work);
}
 
//=============================================================================================================
 
int FwdCompData::fwd_make_ctf_comp_coils(MNECTFCompDataSet* set,
                                         FwdCoilSet* coils,
                                         FwdCoilSet* comp_coils)
{
    QList<FiffChInfo> chs;
    QList<FiffChInfo> compchs;
    int nchan = 0;
    int ncomp = 0;
    int k, res;
 
    if (!set)
        return OK;
    if (!coils || coils->ncoil() <= 0) {
        qWarning("Coil data missing in fwd_make_ctf_comp_coils");
        return FAIL;
    }
 
    for (k = 0; k < coils->ncoil(); k++) {
        chs.append(FiffChInfo());
        FwdCoil* coil = coils->coils[k].get();
        chs[k].ch_name = coil->chname;
        chs[k].chpos.coil_type = coil->type;
        chs[k].kind = (coil->coil_class == FWD_COILC_EEG) ? FIFFV_EEG_CH : FIFFV_MEG_CH;
    }
    nchan = coils->ncoil();
    if (comp_coils && comp_coils->ncoil() > 0) {
        for (k = 0; k < comp_coils->ncoil(); k++) {
            compchs.append(FiffChInfo());
            FwdCoil* coil = comp_coils->coils[k].get();
            compchs[k].ch_name = coil->chname;
            compchs[k].chpos.coil_type = coil->type;
            compchs[k].kind = (coil->coil_class == FWD_COILC_EEG) ? FIFFV_EEG_CH : FIFFV_MEG_CH;
        }
        ncomp = comp_coils->ncoil();
    }
    res = set->make_comp(chs, nchan, compchs, ncomp);
 
    return res;
}
 
//=============================================================================================================
 
std::unique_ptr<FwdCompData> FwdCompData::fwd_make_comp_data(MNECTFCompDataSet* set,
                                                             FwdCoilSet* coils,
                                                             FwdCoilSet* comp_coils,
                                                             fwdFieldFunc field,
                                                             fwdVecFieldFunc vec_field,
                                                             fwdFieldGradFunc field_grad,
                                                             void* client)
{
    auto comp = std::make_unique<FwdCompData>();
 
    if (set)
        comp->set = std::make_shared<MNECTFCompDataSet>(*set);
 
    if (comp_coils)
        comp->comp_coils = comp_coils->dup_coil_set();
    else
        qWarning("No coils to duplicate");
    comp->field = field;
    comp->vec_field = vec_field;
    comp->field_grad = field_grad;
    comp->client = client;
 
    if (fwd_make_ctf_comp_coils(comp->set.get(), coils, comp->comp_coils.get()) != OK)
        return nullptr;
    return comp;
}
 
//=============================================================================================================
 
int FwdCompData::fwd_comp_field_vec(const Eigen::Vector3f& rd, FwdCoilSet& coils, Eigen::Ref<Eigen::MatrixXf> res, void* client)
{
    FwdCompData* comp = static_cast<FwdCompData*>(client);
 
    if (!comp->vec_field) {
        qWarning("Field computation function is missing in fwd_comp_field_vec");
        return FAIL;
    }
    /*
       * First compute the field in the primary set of coils
       */
    if (comp->vec_field(rd, coils, res, comp->client) == FAIL)
        return FAIL;
    /*
       * Compensation needed?
       */
    if (!comp->comp_coils || comp->comp_coils->ncoil() <= 0 || !comp->set || !comp->set->current)
        return OK;
    /*
       * Need workspace?
       */
    if (comp->vec_work.size() == 0)
        comp->vec_work.resize(3, comp->comp_coils->ncoil());
    /*
       * Compute the field at the compensation sensors
       */
    if (comp->vec_field(rd, *comp->comp_coils, comp->vec_work, comp->client) == FAIL)
        return FAIL;
    /*
       * Compute the compensated field of three orthogonal dipoles
       */
    for (int k = 0; k < 3; k++) {
        Eigen::VectorXf resRow = res.row(k).transpose();
        Eigen::VectorXf workRow = comp->vec_work.row(k).transpose();
        if (comp->set->apply(true, resRow, workRow) == FAIL)
            return FAIL;
        res.row(k) = resRow.transpose();
    }
    return OK;
}
 
//=============================================================================================================
 
int FwdCompData::fwd_comp_field_grad(const Eigen::Vector3f& rd, const Eigen::Vector3f& Q, FwdCoilSet& coils, Eigen::Ref<Eigen::VectorXf> res, Eigen::Ref<Eigen::VectorXf> xgrad, Eigen::Ref<Eigen::VectorXf> ygrad, Eigen::Ref<Eigen::VectorXf> zgrad, void* client)
{
    FwdCompData* comp = static_cast<FwdCompData*>(client);
 
    if (!comp->field_grad) {
        qCritical("Field and gradient computation function is missing in fwd_comp_field_grad");
        return FAIL;
    }
    /*
     * First compute the field in the primary set of coils
     */
    if (comp->field_grad(rd, Q, coils, res, xgrad, ygrad, zgrad, comp->client) == FAIL)
        return FAIL;
    /*
     * Compensation needed?
     */
    if (!comp->comp_coils || comp->comp_coils->ncoil() <= 0 || !comp->set || !comp->set->current)
        return OK;
    /*
     * Workspace needed?
     */
    if (comp->work.size() == 0)
        comp->work.resize(comp->comp_coils->ncoil());
    if (comp->vec_work.size() == 0)
        comp->vec_work.resize(3, comp->comp_coils->ncoil());
    /*
     * Compute the field in the compensation coils
     */
    Eigen::VectorXf vw0 = comp->vec_work.row(0).transpose();
    Eigen::VectorXf vw1 = comp->vec_work.row(1).transpose();
    Eigen::VectorXf vw2 = comp->vec_work.row(2).transpose();
    if (comp->field_grad(rd, Q, *comp->comp_coils, comp->work, vw0, vw1, vw2, comp->client) == FAIL)
        return FAIL;
    comp->vec_work.row(0) = vw0.transpose();
    comp->vec_work.row(1) = vw1.transpose();
    comp->vec_work.row(2) = vw2.transpose();
    /*
     * Compute the compensated field
     */
    if (comp->set->apply(true, res, comp->work) != OK)
        return FAIL;
 
    vw0 = comp->vec_work.row(0).transpose();
    if (comp->set->apply(true, xgrad, vw0) != OK)
        return FAIL;
 
    vw1 = comp->vec_work.row(1).transpose();
    if (comp->set->apply(true, ygrad, vw1) != OK)
        return FAIL;
 
    vw2 = comp->vec_work.row(2).transpose();
    if (comp->set->apply(true, zgrad, vw2) != OK)
        return FAIL;
 
    return OK;
}