mne_ctf_comp_data_set.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "mne_ctf_comp_data_set.h"
#include "mne_ctf_comp_data.h"
 
#include "mne_types.h"
 
#include <fiff/fiff_coord_trans.h>
 
#include <fiff/fiff_types.h>
 
#include <Eigen/Core>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QFile>
#include <QDebug>
 
constexpr int FAIL = -1;
constexpr int OK = 0;
 
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace Eigen;
using namespace FIFFLIB;
using namespace MNELIB;
 
//============================= mne_read_forward_solution.c =============================
 
int mne_read_meg_comp_eeg_ch_info_32(const QString& name,
                                     QList<FIFFLIB::FiffChInfo>& megp, /* MEG channels */
                                     int* nmegp,
                                     QList<FIFFLIB::FiffChInfo>& meg_compp,
                                     int* nmeg_compp,
                                     QList<FIFFLIB::FiffChInfo>& eegp, /* EEG channels */
                                     int* neegp,
                                     FiffCoordTrans* meg_head_t,
                                     fiffId* idp) /* The measurement ID */
/*
      * Read the channel information and split it into three arrays,
      * one for MEG, one for MEG compensation channels, and one for EEG
      */
{
    QFile file(name);
    FiffStream::SPtr stream(new FiffStream(&file));
 
    QList<FIFFLIB::FiffChInfo> chs;
    int nchan = 0;
    QList<FIFFLIB::FiffChInfo> meg;
    int nmeg = 0;
    QList<FIFFLIB::FiffChInfo> meg_comp;
    int nmeg_comp = 0;
    QList<FIFFLIB::FiffChInfo> eeg;
    int neeg = 0;
    std::unique_ptr<FiffId> id;
    QList<FiffDirNode::SPtr> nodes;
    FiffDirNode::SPtr info;
    FiffTag::UPtr t_pTag;
    FIFFLIB::FiffChInfo this_ch;
    FiffCoordTrans t;
    fiff_int_t kind, pos;
    int j, k, to_find;
 
    if (!stream->open()) {
        stream->close();
        return FIFF_FAIL;
    }
 
    nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_PARENT_MEAS_FILE);
 
    if (nodes.size() == 0) {
        nodes = stream->dirtree()->dir_tree_find(FIFFB_MEAS_INFO);
        if (nodes.size() == 0) {
            qCritical("Could not find the channel information.");
            stream->close();
            return FIFF_FAIL;
        }
    }
    info = nodes[0];
    to_find = 0;
    for (k = 0; k < info->nent(); k++) {
        kind = info->dir[k]->kind;
        pos = info->dir[k]->pos;
        switch (kind) {
            case FIFF_NCHAN:
                if (!stream->read_tag(t_pTag, pos)) {
                    stream->close();
                    return FIFF_FAIL;
                }
                nchan = *t_pTag->toInt();
 
                for (j = 0; j < nchan; j++) {
                    chs.append(FiffChInfo());
                    chs[j].scanNo = -1;
                }
                to_find = nchan;
                break;
 
            case FIFF_PARENT_BLOCK_ID:
                if (!stream->read_tag(t_pTag, pos)) {
                    stream->close();
                    return FIFF_FAIL;
                }
                id = std::make_unique<FiffId>(*(FiffId*)t_pTag->data());
                break;
 
            case FIFF_COORD_TRANS:
                if (!stream->read_tag(t_pTag, pos)) {
                    stream->close();
                    return FIFF_FAIL;
                }
                //            t = t_pTag->toCoordTrans();
                t = FiffCoordTrans::readFromTag(t_pTag);
                if (t.from != FIFFV_COORD_DEVICE || t.to != FIFFV_COORD_HEAD)
                    t = FiffCoordTrans();
                break;
 
            case FIFF_CH_INFO: /* Information about one channel */
                if (!stream->read_tag(t_pTag, pos)) {
                    stream->close();
                    return FIFF_FAIL;
                }
 
                this_ch = t_pTag->toChInfo();
                if (this_ch.scanNo <= 0 || this_ch.scanNo > nchan) {
                    qCritical("FIFF_CH_INFO : scan # out of range %d (%d)!", this_ch.scanNo, nchan);
                    stream->close();
                    return FIFF_FAIL;
                } else
                    chs[this_ch.scanNo - 1] = this_ch;
                to_find--;
                break;
        }
    }
    if (to_find != 0) {
        qCritical("Some of the channel information was missing.");
        stream->close();
        return FIFF_FAIL;
    }
    if (t.isEmpty() && meg_head_t != nullptr) {
        /*
     * Try again in a more general fashion
     */
        t = FiffCoordTrans::readMeasTransform(name);
        if (t.isEmpty()) {
            qCritical("MEG -> head coordinate transformation not found.");
            stream->close();
            return FIFF_FAIL;
        }
    }
    /*
   * Sort out the channels
   */
    for (k = 0; k < nchan; k++) {
        if (chs[k].kind == FIFFV_MEG_CH) {
            meg.append(chs[k]);
            nmeg++;
        } else if (chs[k].kind == FIFFV_REF_MEG_CH) {
            meg_comp.append(chs[k]);
            nmeg_comp++;
        } else if (chs[k].kind == FIFFV_EEG_CH) {
            eeg.append(chs[k]);
            neeg++;
        }
    }
    //    fiff_close(in);
    stream->close();
 
    megp = meg;
    if (nmegp) {
        *nmegp = nmeg;
    }
 
    meg_compp = meg_comp;
    if (nmeg_compp) {
        *nmeg_compp = nmeg_comp;
    }
 
    eegp = eeg;
    if (neegp) {
        *neegp = neeg;
    }
 
    if (idp == nullptr) {
        /* id is auto-deleted by unique_ptr */
    } else
        *idp = id.release();
    if (meg_head_t == nullptr) {
    } else
        *meg_head_t = t;
 
    return FIFF_OK;
}
 
#define MNE_CTFV_COMP_UNKNOWN -1
#define MNE_CTFV_COMP_NONE 0
#define MNE_CTFV_COMP_G1BR 0x47314252
#define MNE_CTFV_COMP_G2BR 0x47324252
#define MNE_CTFV_COMP_G3BR 0x47334252
#define MNE_CTFV_COMP_G2OI 0x47324f49
#define MNE_CTFV_COMP_G3OI 0x47334f49
 
static struct
{
    int grad_comp;
    int ctf_comp;
} compMap[] = {{MNE_CTFV_NOGRAD, MNE_CTFV_COMP_NONE},
               {MNE_CTFV_GRAD1, MNE_CTFV_COMP_G1BR},
               {MNE_CTFV_GRAD2, MNE_CTFV_COMP_G2BR},
               {MNE_CTFV_GRAD3, MNE_CTFV_COMP_G3BR},
               {MNE_4DV_COMP1, MNE_4DV_COMP1}, /* One-to-one mapping for 4D data */
               {MNE_CTFV_COMP_UNKNOWN, MNE_CTFV_COMP_UNKNOWN}};
 
int mne_unmap_ctf_comp_kind(int ctf_comp)
 
{
    int k;
 
    for (k = 0; compMap[k].grad_comp >= 0; k++)
        if (ctf_comp == compMap[k].ctf_comp)
            return compMap[k].grad_comp;
    return ctf_comp;
}
 
std::unique_ptr<FiffSparseMatrix> mne_convert_to_sparse(const Eigen::MatrixXf& dense, /* The dense matrix to be converted */
                                                        int stor_type,                /* Either FIFFTS_MC_CCS or FIFFTS_MC_RCS */
                                                        float small)                  /* How small elements should be ignored? */
/*
 * Convert a dense matrix to sparse using Eigen's sparseView.
 */
{
    Q_UNUSED(stor_type);
 
    if (small < 0) { /* Automatic scaling */
        float maxval = dense.cwiseAbs().maxCoeff();
        if (maxval > 0)
            small = maxval * std::fabs(small);
        else
            small = std::fabs(small);
    }
 
    Eigen::SparseMatrix<float> eigenSparse = dense.sparseView(small, 1.0f);
    eigenSparse.makeCompressed();
 
    if (eigenSparse.nonZeros() <= 0) {
        qWarning("No nonzero elements found.");
        return nullptr;
    }
 
    return std::make_unique<FiffSparseMatrix>(std::move(eigenSparse), FIFFTS_MC_RCS);
}
 
int mne_sparse_vec_mult2_32(FiffSparseMatrix* mat, /* The sparse matrix */
                            float* vector,         /* Vector to be multiplied */
                            float* res)            /* Result of the multiplication */
/*
      * Multiply a vector by a sparse matrix using Eigen.
      */
{
    Eigen::Map<const Eigen::VectorXf> vecIn(vector, mat->cols());
    Eigen::Map<Eigen::VectorXf> vecOut(res, mat->rows());
    vecOut = mat->eigen() * vecIn;
    return 0;
}
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
MNECTFCompDataSet::MNECTFCompDataSet()
: ncomp(0)
, nch(0)
, undo(nullptr)
, current(nullptr)
{
}
 
//=============================================================================================================
 
MNECTFCompDataSet::MNECTFCompDataSet(const MNECTFCompDataSet& set)
: ncomp(0)
, nch(set.nch)
, undo(nullptr)
, current(nullptr)
{
    if (set.ncomp > 0) {
        for (int k = 0; k < set.ncomp; k++)
            if (set.comps[k])
                this->comps.push_back(std::make_unique<MNECTFCompData>(*set.comps[k]));
        this->ncomp = static_cast<int>(this->comps.size());
    }
 
    this->chs = set.chs;
 
    if (set.undo)
        this->undo = std::make_unique<MNECTFCompData>(*set.undo);
 
    if (set.current)
        this->current = std::make_unique<MNECTFCompData>(*set.current);
}
 
//=============================================================================================================
 
MNECTFCompDataSet::~MNECTFCompDataSet()
{
}
 
//=============================================================================================================
 
std::unique_ptr<MNECTFCompDataSet> MNECTFCompDataSet::read(const QString& name)
/*
     * Read all CTF compensation data from a given file
     */
{
    QFile file(name);
    FiffStream::SPtr stream(new FiffStream(&file));
 
    std::unique_ptr<MNECTFCompDataSet> set;
    QList<FiffDirNode::SPtr> nodes;
    QList<FiffDirNode::SPtr> comps;
    int ncomp;
    int kind, k;
    FiffTag::UPtr t_pTag;
    QList<FiffChInfo> chs;
    int nch = 0;
    int calibrated;
    /*
        * Read the channel information
        */
    {
        QList<FiffChInfo> comp_chs, temp;
        int ncompch = 0;
 
        if (mne_read_meg_comp_eeg_ch_info_32(name,
                                             chs,
                                             &nch,
                                             comp_chs,
                                             &ncompch,
                                             temp,
                                             nullptr,
                                             nullptr,
                                             nullptr) == FAIL)
            return nullptr;
        if (ncompch > 0) {
            for (k = 0; k < ncompch; k++)
                chs.append(comp_chs[k]);
            nch = nch + ncompch;
        }
    }
    /*
        * Read the rest of the stuff
        */
    if (!stream->open()) {
        stream->close();
        return nullptr;
    }
    set = std::make_unique<MNECTFCompDataSet>();
    /*
        * Locate the compensation data sets
        */
    nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_CTF_COMP);
    if (nodes.size() == 0) {
        stream->close();
        return set;
    }
    comps = nodes[0]->dir_tree_find(FIFFB_MNE_CTF_COMP_DATA);
    if (comps.size() == 0) {
        stream->close();
        return set;
    }
    ncomp = comps.size();
    /*
        * Set the channel info
        */
    set->chs = chs;
    set->nch = nch;
    /*
        * Read each data set
        */
    for (k = 0; k < ncomp; k++) {
        auto mat = MNENamedMatrix::read(stream, comps[k], FIFF_MNE_CTF_COMP_DATA);
        if (!mat) {
            stream->close();
            return nullptr;
        }
        comps[k]->find_tag(stream, FIFF_MNE_CTF_COMP_KIND, t_pTag);
        if (t_pTag) {
            kind = *t_pTag->toInt();
        } else {
            stream->close();
            return nullptr;
        }
        comps[k]->find_tag(stream, FIFF_MNE_CTF_COMP_CALIBRATED, t_pTag);
        if (t_pTag) {
            calibrated = *t_pTag->toInt();
        } else
            calibrated = 0;
        /*
            * Add these data to the set
            */
        auto one = std::make_unique<MNECTFCompData>();
        one->data = std::move(mat);
        one->kind = kind;
        one->mne_kind = mne_unmap_ctf_comp_kind(one->kind);
        one->calibrated = calibrated;
 
        if (one->calibrate(set->chs, set->nch, true) == FAIL) {
            qWarning("Warning: Compensation data for '%s' omitted\n", explain_comp(one->kind).toUtf8().constData());
        } else {
            set->comps.push_back(std::move(one));
            set->ncomp++;
        }
    }
#ifdef DEBUG
    qInfo("%d CTF compensation data sets read from %s\n", set->ncomp, name);
#endif
    stream->close();
    return set;
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::make_comp(const QList<FiffChInfo>& chList,
                                 int nChan,
                                 QList<FiffChInfo> compchs,
                                 int nCompChan) /* How many of these */
/*
     * Make compensation data to apply to a set of channels to yield (or uncompensated) compensated data
     */
{
    Eigen::VectorXi compGrades;
    int need_comp;
    int first_comp;
    MNECTFCompData* this_comp;
    Eigen::VectorXi comp_sel;
    QStringList names;
    QString name;
    int j, k, p;
 
    std::unique_ptr<FiffSparseMatrix> presel;
    std::unique_ptr<FiffSparseMatrix> postsel;
    std::unique_ptr<MNENamedMatrix> data;
 
    QStringList emptyList;
 
    if (compchs.isEmpty()) {
        compchs = chList;
        nCompChan = nChan;
    }
    qInfo("Setting up compensation data...\n");
    if (nChan == 0)
        return OK;
    current.reset();
    compGrades.resize(nChan);
    for (k = 0, need_comp = 0, first_comp = MNE_CTFV_COMP_NONE; k < nChan; k++) {
        if (chList[k].kind == FIFFV_MEG_CH) {
            compGrades[k] = chList[k].chpos.coil_type >> 16;
            if (compGrades[k] != MNE_CTFV_COMP_NONE) {
                if (first_comp == MNE_CTFV_COMP_NONE)
                    first_comp = compGrades[k];
                else {
                    if (compGrades[k] != first_comp) {
                        qCritical("We do not support nonuniform compensation yet.");
                        return FAIL;
                    }
                }
                need_comp++;
            }
        } else
            compGrades[k] = MNE_CTFV_COMP_NONE;
    }
    if (need_comp == 0) {
        qInfo("\tNo compensation set. Nothing more to do.\n");
        return OK;
    }
    qInfo("\t%d out of %d channels have the compensation set.\n", need_comp, nChan);
    /*
        * Find the desired compensation data matrix
        */
    for (k = 0, this_comp = nullptr; k < this->ncomp; k++) {
        if (this->comps[k]->mne_kind == first_comp) {
            this_comp = this->comps[k].get();
            break;
        }
    }
    if (!this_comp) {
        qCritical("Did not find the desired compensation data : %s",
                  explain_comp(map_comp_kind(first_comp)).toUtf8().constData());
        return FAIL;
    }
    qInfo("\tDesired compensation data (%s) found.\n", explain_comp(map_comp_kind(first_comp)).toUtf8().constData());
    /*
        * Find the compensation channels
        */
    comp_sel.resize(this_comp->data->ncol);
    for (k = 0; k < this_comp->data->ncol; k++) {
        comp_sel[k] = -1;
        name = this_comp->data->collist[k];
        for (p = 0; p < nCompChan; p++)
            if (QString::compare(name, compchs[p].ch_name) == 0) {
                comp_sel[k] = p;
                break;
            }
        if (comp_sel[k] < 0) {
            qCritical("Compensation channel %s not found", name.toUtf8().constData());
            return FAIL;
        }
    }
    qInfo("\tAll compensation channels found.\n");
    /*
        * Create the preselector
        */
    {
        Eigen::MatrixXf sel = Eigen::MatrixXf::Zero(this_comp->data->ncol, nCompChan);
        for (j = 0; j < this_comp->data->ncol; j++)
            sel(j, comp_sel[j]) = 1.0f;
        presel = mne_convert_to_sparse(sel, FIFFTS_MC_RCS, 1e-30f);
        if (!presel)
            return FAIL;
        qInfo("\tPreselector created.\n");
    }
    /*
     * Pick the desired channels
     */
    for (k = 0; k < nChan; k++) {
        if (compGrades[k] != MNE_CTFV_COMP_NONE)
            names.append(chList[k].ch_name);
    }
 
    {
        auto d = this_comp->data->pick(names, need_comp, emptyList, 0);
        if (!d)
            return FAIL;
        data = std::move(d);
    }
    qInfo("\tCompensation data matrix created.\n");
    /*
        * Create the postselector
        */
    {
        Eigen::MatrixXf sel = Eigen::MatrixXf::Zero(nChan, data->nrow);
        for (j = 0, p = 0; j < nChan; j++) {
            if (compGrades[j] != MNE_CTFV_COMP_NONE)
                sel(j, p++) = 1.0f;
        }
        postsel = mne_convert_to_sparse(sel, FIFFTS_MC_RCS, 1e-30f);
        if (!postsel)
            return FAIL;
        qInfo("\tPostselector created.\n");
    }
    current = std::make_unique<MNECTFCompData>();
    current->kind = this_comp->kind;
    current->mne_kind = this_comp->mne_kind;
    current->data = std::move(data);
    current->presel = std::move(presel);
    current->postsel = std::move(postsel);
 
    qInfo("\tCompensation set up.\n");
    return OK;
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::set_comp(QList<FIFFLIB::FiffChInfo>& chs,
                                int nch,
                                int comp)
/*
     * Set the compensation bits to the desired value
     */
{
    int k;
    int nset;
    for (k = 0, nset = 0; k < nch; k++) {
        if (chs[k].kind == FIFFV_MEG_CH) {
            chs[k].chpos.coil_type = (chs[k].chpos.coil_type & 0xFFFF) | (comp << 16);
            nset++;
        }
    }
    qInfo("A new compensation value (%s) was assigned to %d MEG channels.\n",
          explain_comp(map_comp_kind(comp)).toUtf8().constData(), nset);
    return nset;
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::apply(bool do_it, Eigen::Ref<Eigen::VectorXf> data)
{
    return apply(do_it, data, data);
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::apply(bool do_it, Eigen::Ref<Eigen::VectorXf> data, Eigen::Ref<const Eigen::VectorXf> compdata)
/*
     * Apply compensation or revert to uncompensated data
     */
{
    MNECTFCompData* this_comp;
    int ndata = static_cast<int>(data.size());
    int ncompdata = static_cast<int>(compdata.size());
 
    if (!current)
        return OK;
    this_comp = current.get();
    /*
       * Dimension checks
       */
    if (this_comp->presel) {
        if (this_comp->presel->cols() != ncompdata) {
            qCritical("Compensation data dimension mismatch. Expected %d, got %d channels.",
                      this_comp->presel->cols(), ncompdata);
            return FAIL;
        }
    } else if (this_comp->data->ncol != ncompdata) {
        qCritical("Compensation data dimension mismatch. Expected %d, got %d channels.",
                  this_comp->data->ncol, ncompdata);
        return FAIL;
    }
    if (this_comp->postsel) {
        if (this_comp->postsel->rows() != ndata) {
            qCritical("Data dimension mismatch. Expected %d, got %d channels.",
                      this_comp->postsel->rows(), ndata);
            return FAIL;
        }
    } else if (this_comp->data->nrow != ndata) {
        qCritical("Data dimension mismatch. Expected %d, got %d channels.",
                  this_comp->data->nrow, ndata);
        return FAIL;
    }
    /*
        * Preselection is optional
        */
    const float* presel;
    if (this_comp->presel) {
        if (this_comp->presel_data.size() == 0)
            this_comp->presel_data.resize(this_comp->presel->rows());
        if (mne_sparse_vec_mult2_32(this_comp->presel.get(), const_cast<float*>(compdata.data()), this_comp->presel_data.data()) != OK)
            return FAIL;
        presel = this_comp->presel_data.data();
    } else
        presel = compdata.data();
    /*
        * This always happens
        */
    if (this_comp->comp_data.size() == 0)
        this_comp->comp_data.resize(this_comp->data->nrow);
    {
        Eigen::Map<const Eigen::VectorXf> preselVec(presel, this_comp->data->ncol);
        Eigen::Map<Eigen::VectorXf> compVec(this_comp->comp_data.data(), this_comp->data->nrow);
        compVec = this_comp->data->data * preselVec;
    }
    /*
        * Optional postselection
        */
    const float* comp;
    if (!this_comp->postsel)
        comp = this_comp->comp_data.data();
    else {
        if (this_comp->postsel_data.size() == 0)
            this_comp->postsel_data.resize(this_comp->postsel->rows());
        if (mne_sparse_vec_mult2_32(this_comp->postsel.get(), this_comp->comp_data.data(), this_comp->postsel_data.data()) != OK)
            return FAIL;
        comp = this_comp->postsel_data.data();
    }
    /*
        * Compensate or revert compensation?
        */
    Eigen::Map<const Eigen::VectorXf> compVec(comp, ndata);
    if (do_it)
        data -= compVec;
    else
        data += compVec;
    return OK;
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::apply_transpose(bool do_it, Eigen::MatrixXf& data)
/*
     * Apply compensation or revert to uncompensated data
     */
{
    MNECTFCompData* this_comp;
    int ndata = static_cast<int>(data.rows());
    int ncompdata = ndata;
 
    if (!current)
        return OK;
    this_comp = current.get();
    /*
        * Dimension checks
        */
    if (this_comp->presel) {
        if (this_comp->presel->cols() != ncompdata) {
            qCritical("Compensation data dimension mismatch. Expected %d, got %d channels.",
                      this_comp->presel->cols(), ncompdata);
            return FAIL;
        }
    } else if (this_comp->data->ncol != ncompdata) {
        qCritical("Compensation data dimension mismatch. Expected %d, got %d channels.",
                  this_comp->data->ncol, ncompdata);
        return FAIL;
    }
    if (this_comp->postsel) {
        if (this_comp->postsel->rows() != ndata) {
            qCritical("Data dimension mismatch. Expected %d, got %d channels.",
                      this_comp->postsel->rows(), ndata);
            return FAIL;
        }
    } else if (this_comp->data->nrow != ndata) {
        qCritical("Data dimension mismatch. Expected %d, got %d channels.",
                  this_comp->data->nrow, ndata);
        return FAIL;
    }
    /*
        * Preselection is optional — sparse matrix * data
        */
    Eigen::MatrixXf preselMat;
    if (this_comp->presel) {
        preselMat = this_comp->presel->eigen() * data;
    } else {
        preselMat = data;
    }
    /*
        * Compensation: comp = data_matrix * preselMat
        * data_matrix is float** (contiguous row-major via mne_cmatrix)
        */
    Eigen::MatrixXf comp = this_comp->data->data * preselMat;
    /*
        * Optional postselection — sparse matrix * comp
        */
    if (this_comp->postsel) {
        comp = this_comp->postsel->eigen() * comp;
    }
    /*
        * Compensate or revert compensation?
        */
    if (do_it)
        data -= comp;
    else
        data += comp;
    return OK;
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::get_comp(const QList<FIFFLIB::FiffChInfo>& chs, int nch)
{
    int res = MNE_CTFV_NOGRAD;
    int first_comp, comp;
    int k;
 
    for (k = 0, first_comp = -1; k < nch; k++) {
        if (chs[k].kind == FIFFV_MEG_CH) {
            comp = chs[k].chpos.coil_type >> 16;
            if (first_comp < 0)
                first_comp = comp;
            else if (first_comp != comp) {
                qCritical("Non uniform compensation not supported.");
                return FAIL;
            }
        }
    }
    if (first_comp >= 0)
        res = first_comp;
    return res;
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::map_comp_kind(int grad)
/*
     * Simple mapping
     */
{
    int k;
 
    for (k = 0; compMap[k].grad_comp >= 0; k++)
        if (grad == compMap[k].grad_comp)
            return compMap[k].ctf_comp;
    return grad;
}
 
//=============================================================================================================
 
QString MNECTFCompDataSet::explain_comp(int kind)
{
    static const struct
    {
        int kind;
        const char* expl;
    } explain[] = {{MNE_CTFV_COMP_NONE, "uncompensated"},
                   {MNE_CTFV_COMP_G1BR, "first order gradiometer"},
                   {MNE_CTFV_COMP_G2BR, "second order gradiometer"},
                   {MNE_CTFV_COMP_G3BR, "third order gradiometer"},
                   {MNE_4DV_COMP1, "4D comp 1"},
                   {MNE_CTFV_COMP_UNKNOWN, "unknown"}};
    int k;
 
    for (k = 0; explain[k].kind != MNE_CTFV_COMP_UNKNOWN; k++)
        if (explain[k].kind == kind)
            return QString::fromLatin1(explain[k].expl);
    return QString::fromLatin1(explain[k].expl);
}
 
//=============================================================================================================
 
int MNECTFCompDataSet::set_compensation(int compensate_to,
                                        QList<FiffChInfo>& chList,
                                        int nchan,
                                        QList<FiffChInfo> comp_chs,
                                        int ncomp_chan) /* How many */
/*
     * Make data which has the third-order gradient compensation applied
     */
{
    int k;
    int have_comp_chs;
    int comp_was = MNE_CTFV_COMP_UNKNOWN;
 
    if (comp_chs.isEmpty()) {
        comp_chs = chList;
        ncomp_chan = nchan;
    }
    undo.reset();
    current.reset();
    for (k = 0, have_comp_chs = 0; k < ncomp_chan; k++)
        if (comp_chs[k].kind == FIFFV_REF_MEG_CH)
            have_comp_chs++;
    if (have_comp_chs == 0 && compensate_to != MNE_CTFV_NOGRAD) {
        qWarning("No compensation channels in these data.");
        return FAIL;
    }
    /*
        * Update the 'current' field to reflect the compensation possibly present in the data now
        */
    if (make_comp(chList, nchan, comp_chs, ncomp_chan) == FAIL)
        return FAIL;
    /*
        * Are we there already?
        */
    if (current && current->mne_kind == compensate_to) {
        qInfo("No further compensation necessary (comp = %s)\n", explain_comp(current->kind).toUtf8().constData());
        current.reset();
        return OK;
    }
    undo = std::move(current);
    if (compensate_to == MNE_CTFV_NOGRAD) {
        qInfo("No compensation was requested.\n");
        set_comp(chList, nchan, compensate_to);
        return OK;
    }
    if (set_comp(chList, nchan, compensate_to) > 0) {
        if (undo)
            comp_was = undo->mne_kind;
        else
            comp_was = MNE_CTFV_NOGRAD;
        if (make_comp(chList, nchan, comp_chs, ncomp_chan) == FAIL) {
            if (comp_was != MNE_CTFV_COMP_UNKNOWN)
                set_comp(chList, nchan, comp_was);
            return FAIL;
        }
        qInfo("Compensation set up as requested (%s -> %s).\n",
              explain_comp(map_comp_kind(comp_was)).toUtf8().constData(),
              explain_comp(current->kind).toUtf8().constData());
    }
    return OK;
}