fiff_info.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "fiff_info.h"
#include "fiff_stream.h"
#include "fiff_file.h"
 
#include <utils/ioutils.h>
#include <iostream>
#include <QDebug>
 
#include <Eigen/LU>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace FIFFLIB;
using namespace UTILSLIB;
using namespace Eigen;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
FiffInfo::FiffInfo()
: FiffInfoBase() //nchan(-1)
, sfreq(-1.0)
, linefreq(-1.0)
, highpass(-1.0)
, lowpass(-1.0)
, gantry_angle(-1)
, acq_pars("")
, acq_stim("")
{
    meas_date[0] = -1;
}
 
//=============================================================================================================
 
FiffInfo::FiffInfo(const FiffInfo& p_FiffInfo)
: FiffInfoBase(p_FiffInfo)
, file_id(p_FiffInfo.file_id)
, sfreq(p_FiffInfo.sfreq)
, linefreq(p_FiffInfo.linefreq)
, highpass(p_FiffInfo.highpass)
, lowpass(p_FiffInfo.lowpass)
, proj_id(p_FiffInfo.proj_id)
, proj_name(p_FiffInfo.proj_name)
, xplotter_layout(p_FiffInfo.xplotter_layout)
, experimenter(p_FiffInfo.experimenter)
, description(p_FiffInfo.description)
, utc_offset(p_FiffInfo.utc_offset)
, gantry_angle(p_FiffInfo.gantry_angle)
, dev_ctf_t(p_FiffInfo.dev_ctf_t)
, dig(p_FiffInfo.dig)
, dig_trans(p_FiffInfo.dig_trans)
, projs(p_FiffInfo.projs)
, comps(p_FiffInfo.comps)
, acq_pars(p_FiffInfo.acq_pars)
, acq_stim(p_FiffInfo.acq_stim)
, hpi_coil_freqs(p_FiffInfo.hpi_coil_freqs)
{
    meas_date[0] = p_FiffInfo.meas_date[0];
    meas_date[1] = p_FiffInfo.meas_date[1];
}
 
//=============================================================================================================
 
FiffInfo::~FiffInfo()
{
}
 
//=============================================================================================================
 
void FiffInfo::clear()
{
    FiffInfoBase::clear();
    file_id = FiffId();
    meas_date[0] = -1;
    sfreq = -1.0;
    linefreq = -1.0;
    highpass = -1.0;
    lowpass = -1.0;
    proj_id = -1;
    proj_name = "";
    xplotter_layout = "";
    experimenter = "";
    description = "";
    utc_offset = ""; 
    gantry_angle = -1;
    dev_ctf_t.clear();
    dig.clear();
    dig_trans.clear();
    projs.clear();
    comps.clear();
    acq_pars = "";
    acq_stim = "";
    hpi_coil_freqs.clear();
}
 
//=============================================================================================================
 
qint32 FiffInfo::get_current_comp()
{
    qint32 comp = 0;
    qint32 first_comp = -1;
 
    qint32 k = 0;
    for (k = 0; k < this->nchan; ++k) {
        if (this->chs[k].kind == FIFFV_MEG_CH) {
            comp = this->chs[k].chpos.coil_type >> 16;
            if (first_comp < 0)
                first_comp = comp;
            else if (comp != first_comp)
                qWarning("Compensation is not set equally on all MEG channels");
        }
    }
    return comp;
}
 
//=============================================================================================================
 
bool FiffInfo::make_compensator(fiff_int_t from, fiff_int_t to, FiffCtfComp& ctf_comp, bool exclude_comp_chs) const
{
    MatrixXd C1, C2, comp_tmp;
 
    //    if(ctf_comp.data)
    //        delete ctf_comp.data;
    ctf_comp.data->clear();
 
    if (from == to) {
        ctf_comp.data->data = MatrixXd::Identity(this->nchan, this->nchan);
        return false;
    }
 
    if (from == 0)
        C1 = MatrixXd::Zero(this->nchan, this->nchan);
    else {
        if (!this->make_compensator(from, C1)) {
            qWarning("Cannot create compensator C1\n");
            qWarning("Desired compensation matrix (kind = %d) not found\n", from);
            return false;
        }
    }
 
    if (to == 0)
        C2 = MatrixXd::Zero(this->nchan, this->nchan);
    else {
        if (!this->make_compensator(to, C2)) {
            qWarning("Cannot create compensator C2\n");
            qWarning("Desired compensation matrix (kind = %d) not found\n", to);
            return false;
        }
    }
    //
    //   s_from = (I - C1)*s_orig  =>  s_orig = (I - C1)^-1 * s_from
    //   s_to   = (I - C2)*s_orig
    //   (I + C1) is only the inverse when C1*C1 = 0, which fails as soon as
    //   reference channels are themselves compensated (e.g. CTF grade 1).
    //
    const MatrixXd eye = MatrixXd::Identity(this->nchan, this->nchan);
    comp_tmp = (eye - C2) * (eye - C1).partialPivLu().inverse();
 
    qint32 k;
    if (exclude_comp_chs) {
        VectorXi pick = VectorXi::Zero(this->nchan);
        qint32 npick = 0;
        for (k = 0; k < this->nchan; ++k) {
            if (this->chs[k].kind != FIFFV_REF_MEG_CH) {
                pick(npick) = k;
                ++npick;
            }
        }
        if (npick == 0) {
            qWarning("Nothing remains after excluding the compensation channels\n");
            return false;
        }
 
        ctf_comp.data->data.resize(npick, this->nchan);
        for (k = 0; k < npick; ++k)
            ctf_comp.data->data.row(k) = comp_tmp.block(pick(k), 0, 1, this->nchan);
    } else {
        ctf_comp.data->data = comp_tmp;
    }
    // FiffRawData applies the compensator only when kind is set.
    ctf_comp.kind = to;
 
    return true;
}
 
//=============================================================================================================
 
bool FiffInfo::make_compensator(fiff_int_t kind, MatrixXd& this_comp) const //private method
{
    FiffNamedMatrix::SDPtr this_data;
    MatrixXd presel, postsel;
    qint32 k, col, c, ch = 0, row, row_ch = 0, channelAvailable;
 
    for (k = 0; k < this->comps.size(); ++k) {
        if (this->comps[k].kind == kind) {
            this_data = this->comps[k].data;
 
            //
            //   Create the preselector
            //
            presel = MatrixXd::Zero(this_data->ncol, this->nchan);
 
            for (col = 0; col < this_data->ncol; ++col) {
                channelAvailable = 0;
                for (c = 0; c < this->ch_names.size(); ++c) {
                    if (QString::compare(this_data->col_names.at(col), this->ch_names.at(c)) == 0) {
                        ++channelAvailable;
                        ch = c;
                    }
                }
                if (channelAvailable == 0) {
                    qWarning("Channel %s is not available in data\n", this_data->col_names.at(col).toUtf8().constData());
                    return false;
                } else if (channelAvailable > 1) {
                    qWarning("Ambiguous channel %s", this_data->col_names.at(col).toUtf8().constData());
                    return false;
                }
                presel(col, ch) = 1.0;
            }
            //
            //   Create the postselector
            //
            postsel = MatrixXd::Zero(this->nchan, this_data->nrow);
 
            for (c = 0; c < this->nchan; ++c) {
                channelAvailable = 0;
                for (row = 0; row < this_data->row_names.size(); ++row) {
                    if (QString::compare(this->ch_names.at(c), this_data->row_names.at(row)) == 0) {
                        ++channelAvailable;
                        row_ch = row;
                    }
                }
                if (channelAvailable > 1) {
                    qWarning("Ambiguous channel %s", this->ch_names.at(c).toUtf8().constData());
                    return false;
                } else if (channelAvailable == 1) {
                    postsel(c, row_ch) = 1.0;
                }
            }
            this_comp = postsel * this_data->data * presel;
            return true;
        }
    }
    this_comp = defaultMatrixXd;
    return false;
}
 
//=============================================================================================================
 
FiffInfo FiffInfo::pick_info(const RowVectorXi& sel) const
{
    FiffInfo res = *this; //new FiffInfo(this);
    if (sel.size() == 0)
        return res;
 
    //ToDo when pointer List do delation
    res.chs.clear();
    res.ch_names.clear();
 
    qint32 idx;
    for (qint32 i = 0; i < sel.size(); ++i) {
        idx = sel[i];
        res.chs.append(this->chs[idx]);
        res.ch_names.append(this->ch_names[idx]);
    }
    res.nchan = static_cast<int>(sel.size());
 
    return res;
}
 
//=============================================================================================================
 
QList<FiffChInfo> FiffInfo::set_current_comp(QList<FiffChInfo>& listFiffChInfo, fiff_int_t value)
{
    QList<FiffChInfo> newList;
    qint32 k;
    fiff_int_t coil_type;
 
    for (k = 0; k < listFiffChInfo.size(); ++k)
        newList.append(listFiffChInfo[k]);
 
    qint32 lower_half = 65535; // hex2dec('FFFF');
    for (k = 0; k < listFiffChInfo.size(); ++k) {
        if (listFiffChInfo[k].kind == FIFFV_MEG_CH) {
            coil_type = listFiffChInfo[k].chpos.coil_type & lower_half;
            newList[k].chpos.coil_type = (coil_type | (value << 16));
        }
    }
    return newList;
}
 
//=============================================================================================================
 
bool FiffInfo::readMegEegChannels(const QString& name,
                                  bool do_meg,
                                  bool do_eeg,
                                  const QStringList& bads,
                                  QList<FiffChInfo>& chsp,
                                  int& nmegp,
                                  int& neegp)
{
    QFile file(name);
    FiffStream::SPtr stream(new FiffStream(&file));
 
    if (!stream->open())
        return false;
 
    FiffInfo info;
    FiffDirNode::SPtr infoNode;
    if (!stream->read_meas_info(stream->dirtree(), info, infoNode)) {
        qCritical("%s : could not read measurement info", name.toUtf8().data());
        stream->close();
        return false;
    }
    stream->close();
 
    QList<FiffChInfo> meg;
    QList<FiffChInfo> eeg;
 
    for (int k = 0; k < info.chs.size(); k++) {
        if (bads.contains(info.chs[k].ch_name))
            continue;
        if (do_meg && info.chs[k].kind == FIFFV_MEG_CH)
            meg.append(info.chs[k]);
        else if (do_eeg && info.chs[k].isValidEeg())
            eeg.append(info.chs[k]);
    }
 
    chsp.clear();
    chsp.reserve(meg.size() + eeg.size());
    chsp.append(meg);
    chsp.append(eeg);
 
    nmegp = meg.size();
    neegp = eeg.size();
    return true;
}
 
//=============================================================================================================
 
void FiffInfo::writeToStream(FiffStream* p_pStream) const
{
    //
    //   We will always write floats
    //
    fiff_int_t data_type = 4;
    QList<FiffChInfo> chsToWrite;
 
    for (qint32 k = 0; k < this->nchan; ++k)
        chsToWrite << this->chs[k];
 
    fiff_int_t nchanToWrite = chsToWrite.size();
 
    //
    // write the essentials
    //
    p_pStream->start_block(FIFFB_MEAS); //4
    p_pStream->write_id(FIFF_BLOCK_ID); //5
    if (this->meas_id.version != -1) {
        p_pStream->write_id(FIFF_PARENT_BLOCK_ID, this->meas_id); //6
    }
    //
    //    Measurement info
    //
    p_pStream->start_block(FIFFB_MEAS_INFO); //7
 
    //
    //    Blocks from the original -> skip this
    //
    //        QList<fiff_int_t> blocks;
    //        blocks << FIFFB_SUBJECT << FIFFB_HPI_MEAS << FIFFB_HPI_RESULT << FIFFB_ISOTRAK << FIFFB_PROCESSING_HISTORY;
    bool have_hpi_result = false;
    bool have_isotrak = false;
    //
    //    megacq parameters
    //
    if (!this->acq_pars.isEmpty() || !this->acq_stim.isEmpty()) {
        p_pStream->start_block(FIFFB_DACQ_PARS);
        if (!this->acq_pars.isEmpty())
            p_pStream->write_string(FIFF_DACQ_PARS, this->acq_pars);
 
        if (!this->acq_stim.isEmpty())
            p_pStream->write_string(FIFF_DACQ_STIM, this->acq_stim);
 
        p_pStream->end_block(FIFFB_DACQ_PARS);
    }
    //
    //    Coordinate transformations if the HPI result block was not there
    //
    if (!have_hpi_result) {
        if (!this->dev_head_t.isEmpty())
            p_pStream->write_coord_trans(this->dev_head_t);
 
        if (!this->ctf_head_t.isEmpty())
            p_pStream->write_coord_trans(this->ctf_head_t);
    }
    //
    //    Polhemus data
    //
    if (this->dig.size() > 0 && !have_isotrak) {
        p_pStream->start_block(FIFFB_ISOTRAK);
        for (qint32 k = 0; k < this->dig.size(); ++k)
            p_pStream->write_dig_point(this->dig[k]);
 
        p_pStream->end_block(FIFFB_ISOTRAK);
    }
    //
    //    Projectors
    //
    p_pStream->write_proj(this->projs);
    //
    //    CTF compensation info
    //
    p_pStream->write_ctf_comp(this->comps);
    //
    //    Bad channels
    //
    if (this->bads.size() > 0) {
        p_pStream->start_block(FIFFB_MNE_BAD_CHANNELS);
        p_pStream->write_name_list(FIFF_MNE_CH_NAME_LIST, this->bads);
        p_pStream->end_block(FIFFB_MNE_BAD_CHANNELS);
    }
    //
    //    General
    //
    p_pStream->write_float(FIFF_SFREQ, &this->sfreq);
    p_pStream->write_float(FIFF_LINE_FREQ, &this->linefreq);
    p_pStream->write_float(FIFF_HIGHPASS, &this->highpass);
    p_pStream->write_float(FIFF_LOWPASS, &this->lowpass);
    p_pStream->write_string(FIFF_EXPERIMENTER, this->experimenter);
    p_pStream->write_string(FIFF_DESCRIPTION, this->description);
    p_pStream->write_string(FIFF_UTC_OFFSET, this->utc_offset);
    p_pStream->write_string(FIFF_PROJ_NAME, this->proj_name);
    p_pStream->write_int(FIFF_PROJ_ID, &this->proj_id);
    p_pStream->write_int(FIFF_GANTRY_ANGLE, &this->gantry_angle);
    p_pStream->write_int(FIFF_NCHAN, &nchanToWrite);
    p_pStream->write_int(FIFF_DATA_PACK, &data_type);
    if (this->meas_date[0] != -1)
        p_pStream->write_int(FIFF_MEAS_DATE, this->meas_date, 2);
    //
    //    Channel info
    //
    MatrixXd cals(1, nchanToWrite);
 
    for (qint32 k = 0; k < nchanToWrite; ++k) {
        //
        //    Scan numbers may have been messed up
        //
        chsToWrite[k].scanNo = k + 1;                        //+1 because
                                                             //        chs[k].range  = 1.0f;//Why? -> cause its already calibrated through reading
        cals(0, k) = static_cast<double>(chsToWrite[k].cal); //ToDo whats going on with cals?
        p_pStream->write_ch_info(chsToWrite[k]);
    }
    //
    //
    p_pStream->end_block(FIFFB_MEAS_INFO);
}
 
//=============================================================================================================
 
void FiffInfo::print() const
{
    std::cout << "Sample frequency: " << sfreq << "\n";
    std::cout << "LineFreq: " << linefreq << " | Highpass: " << highpass << " | Lowpass: " << lowpass << "\n";
    std::cout << "Number of digitizer points: " << dig.size() << "\n";
    for (auto& point : dig) {
        if (point.kind == FIFFV_POINT_HPI) {
            std::cout << "HPI Point " << point.ident << " - " << point.r[0] << ", " << point.r[1] << ", " << point.r[2] << "\n";
        }
    }
}