inv_ecd_set.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "inv_ecd_set.h"
#include <fiff/fiff_byte_swap.h>
#include <fiff/fiff_types.h>
 
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QDebug>
#include <QFile>
#include <QTextStream>
#include <QString>
#include <QRegularExpression>
#include <QDebug>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace INVLIB;
using namespace FIFFLIB;
 
//=============================================================================================================
// DEFINES
//=============================================================================================================
 
constexpr int X = 0;
constexpr int Y = 1;
constexpr int Z = 2;
 
namespace INVLIB
{

typedef struct
{
    int dipole;               /* Which dipole in a multi-dipole set */
    float begin, end;         /* Fitting time range */
    float r0[3];              /* Sphere model origin */
    float rd[3];              /* Dipole location */
    float Q[3];               /* InvDipole amplitude */
    float goodness;           /* Goodness-of-fit */
    int errors_computed;      /* Have we computed the errors */
    float noise_level;        /* Noise level used for error computations */
    float single_errors[5];   /* Single parameter error limits */
    float error_matrix[5][5]; /* This fully describes the conf. ellipsoid */
    float conf_vol;           /* The xyz confidence volume */
    float khi2;               /* The khi^2 value */
    float prob;               /* Probability to exceed khi^2 by chance */
    float noise_est;          /* Total noise estimate */
} bdipEcdRec;
using bdipEcd = bdipEcdRec*;
 
} // Namespace
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
InvEcdSet::InvEcdSet()
{
}
 
//=============================================================================================================
 
InvEcdSet::InvEcdSet(const InvEcdSet& p_ECDSet)
: dataname(p_ECDSet.dataname)
, m_qListDips(p_ECDSet.m_qListDips)
{
}
 
//=============================================================================================================
 
InvEcdSet::~InvEcdSet()
{
}
 
//=============================================================================================================
 
void InvEcdSet::addEcd(const InvEcd& p_ecd)
{
    m_qListDips.append(p_ecd);
}
 
//=============================================================================================================
 
InvEcdSet InvEcdSet::read_dipoles_dip(const QString& fileName)
{
    InvEcdSet set;
 
    QFile inputFile(fileName);
    if (inputFile.open(QIODevice::ReadOnly | QIODevice::Text)) {
        QTextStream in(&inputFile);
        while (!in.atEnd()) {
            QString line = in.readLine();
            QStringList list = line.split(QRegularExpression("\\s+"));
 
            if (list[0].contains("#") || list.size() != 11) {
                continue;
            } else {
                InvEcd one;
                one.valid = true;
                one.time = list[1].toFloat() / 1000.0f;
                one.rd[X] = list[3].toFloat() / 1000.0f;
                one.rd[Y] = list[4].toFloat() / 1000.0f;
                one.rd[Z] = list[5].toFloat() / 1000.0f;
                one.Q[X] = list[7].toFloat() / 1e9f;
                one.Q[Y] = list[8].toFloat() / 1e9f;
                one.Q[Z] = list[9].toFloat() / 1e9f;
                one.good = list[10].toFloat() / 100.0f;
                set << one;
            }
        }
        inputFile.close();
 
        qInfo("Read %d dipoles in dip format from %s\n", set.size(), fileName.toUtf8().data());
    } else {
        qCritical("Not able to read from: %s\n", fileName.toUtf8().data());
    }
 
    return set;
}
 
//=============================================================================================================
 
bool InvEcdSet::save_dipoles_bdip(const QString& fileName)
/*
   * Save dipoles in the bdip format employed by xfit
   */
{
    QFile out(fileName);
    bdipEcdRec one_out;
    InvEcd one;
    int k, p;
    int nsave;
 
    if (fileName.isEmpty() || this->size() == 0)
        return true;
 
    if (!out.open(QIODevice::WriteOnly)) {
        qInfo("%s", fileName.toUtf8().constData());
        return false;
    }
 
    for (k = 0, nsave = 0; k < this->size(); k++) {
        one = m_qListDips[k];
        if (one.valid) {
            one_out.dipole = FIFFLIB::swap_int(1);
            one_out.begin = FIFFLIB::swap_float(one.time);
            for (p = 0; p < 3; p++) {
                one_out.r0[p] = FIFFLIB::swap_float(0.0);
                one_out.rd[p] = FIFFLIB::swap_float(one.rd[p]);
                one_out.Q[p] = FIFFLIB::swap_float(one.Q[p]);
            }
            one_out.goodness = FIFFLIB::swap_float(one.good);
            one_out.errors_computed = FIFFLIB::swap_int(0);
            one_out.khi2 = FIFFLIB::swap_float(one.khi2);
            if (out.write(reinterpret_cast<const char*>(&one_out), sizeof(bdipEcdRec)) != sizeof(bdipEcdRec)) {
                qCritical("Failed to write a dipole");
                out.close();
                QFile::remove(fileName);
                return false;
            }
            nsave++;
        }
    }
    if (!out.flush()) {
        qInfo("%s", fileName.toUtf8().constData());
        return false;
    }
    out.close();
    qInfo("Save %d dipoles in bdip format to %s\n", nsave, fileName.toUtf8().data());
    return true;
}
 
//=============================================================================================================
 
bool InvEcdSet::save_dipoles_dip(const QString& fileName) const
{
    QFile out(fileName);
    int k, nsave;
    InvEcd one;
 
    if (fileName.isEmpty() || this->size() == 0)
        return true;
    if (!out.open(QIODevice::WriteOnly | QIODevice::Text)) {
        qInfo("%s", fileName.toUtf8().constData());
        return false;
    }
    //
    // The column layout is read back by other tools, so keep printf style
    // formatting rather than rebuilding the widths with stream manipulators.
    //
    QTextStream stream(&out);
    stream << "# CoordinateSystem \"Head\"\n";
    stream << QString::asprintf("# %7s %7s %8s %8s %8s %8s %8s %8s %8s %6s\n",
                                "begin", "end", "X (mm)", "Y (mm)", "Z (mm)", "Q(nAm)", "Qx(nAm)", "Qy(nAm)", "Qz(nAm)", "g/%");
    for (k = 0, nsave = 0; k < this->size(); k++) {
        one = this->m_qListDips[k];
        if (one.valid) {
            stream << QString::asprintf("  %7.1f %7.1f %8.2f %8.2f %8.2f %8.3f %8.3f %8.3f %8.3f %6.1f\n",
                                        1000 * one.time, 1000 * one.time,
                                        1000 * one.rd[X], 1000 * one.rd[Y], 1000 * one.rd[Z],
                                        1e9 * one.Q.norm(), 1e9 * one.Q[X], 1e9 * one.Q[Y], 1e9 * one.Q[Z], 100.0 * one.good);
            nsave++;
        }
    }
    stream << QString::asprintf("## Name \"%s dipoles\" Style \"Dipoles\"\n", "ALL");
    stream.flush();
    if (!out.flush()) {
        qInfo("%s", fileName.toUtf8().constData());
        return false;
    }
    out.close();
    qInfo("Save %d dipoles in dip format to %s\n", nsave, fileName.toUtf8().data());
    return true;
}
 
//=============================================================================================================
 
const InvEcd& InvEcdSet::operator[](int idx) const
{
    if (idx >= m_qListDips.length()) {
        qWarning("Warning: Required InvEcd doesn't exist! Returning InvEcd '0'.");
        idx = 0;
    }
    return m_qListDips[idx];
}
 
//=============================================================================================================
 
InvEcd& InvEcdSet::operator[](int idx)
{
    if (idx >= m_qListDips.length()) {
        qWarning("Warning: Required InvEcd doesn't exist! Returning InvEcd '0'.");
        idx = 0;
    }
    return m_qListDips[idx];
}
 
//=============================================================================================================
 
InvEcdSet& InvEcdSet::operator<<(const InvEcd& p_ecd)
{
    this->m_qListDips.append(p_ecd);
    return *this;
}