fs_surface.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "fs_surface.h"
#include <fiff/fiff_byte_swap.h>
 
#include <iostream>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QFile>
#include <QDataStream>
#include <QTextStream>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace FSLIB;
using namespace Eigen;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
FsSurface::FsSurface()
: m_sFilePath("")
, m_sFileName("")
, m_iHemi(-1)
, m_sSurf("")
, m_vecOffset(Vector3f::Zero(3))
{
}
 
//=============================================================================================================
 
FsSurface::FsSurface(const QString& p_sFile)
: m_sFilePath("")
, m_sFileName("")
, m_iHemi(-1)
, m_sSurf("")
, m_vecOffset(Vector3f::Zero(3))
{
    FsSurface::read(p_sFile, *this);
}
 
//=============================================================================================================
 
FsSurface::FsSurface(const QString& subject_id, qint32 hemi, const QString& surf, const QString& subjects_dir)
: m_sFilePath("")
, m_sFileName("")
, m_iHemi(-1)
, m_sSurf("")
, m_vecOffset(Vector3f::Zero(3))
{
    FsSurface::read(subject_id, hemi, surf, subjects_dir, *this);
}
 
//=============================================================================================================
 
FsSurface::FsSurface(const QString& path, qint32 hemi, const QString& surf)
: m_sFilePath("")
, m_sFileName("")
, m_iHemi(-1)
, m_sSurf("")
, m_vecOffset(Vector3f::Zero(3))
{
    FsSurface::read(path, hemi, surf, *this);
}
 
//=============================================================================================================
 
FsSurface::~FsSurface()
{
}
 
//=============================================================================================================
 
void FsSurface::clear()
{
    m_sFilePath.clear();
    m_sFileName.clear();
    m_iHemi = -1;
    m_sSurf.clear();
    m_matRR.resize(0, 3);
    m_matTris.resize(0, 3);
    m_matNN.resize(0, 3);
    m_vecCurv.resize(0);
}
 
//=============================================================================================================
 
MatrixX3f FsSurface::compute_normals(const MatrixX3f& rr, const MatrixX3i& tris)
{
    qInfo("\tcomputing normals\n");
    // first, compute triangle normals
    MatrixX3f r1(tris.rows(), 3);
    MatrixX3f r2(tris.rows(), 3);
    MatrixX3f r3(tris.rows(), 3);
 
    for (qint32 i = 0; i < tris.rows(); ++i) {
        r1.row(i) = rr.row(tris(i, 0));
        r2.row(i) = rr.row(tris(i, 1));
        r3.row(i) = rr.row(tris(i, 2));
    }
 
    MatrixX3f x = r2 - r1;
    MatrixX3f y = r3 - r1;
    MatrixX3f tri_nn(x.rows(), y.cols());
    tri_nn.col(0) = x.col(1).cwiseProduct(y.col(2)) - x.col(2).cwiseProduct(y.col(1));
    tri_nn.col(1) = x.col(2).cwiseProduct(y.col(0)) - x.col(0).cwiseProduct(y.col(2));
    tri_nn.col(2) = x.col(0).cwiseProduct(y.col(1)) - x.col(1).cwiseProduct(y.col(0));
 
    //   Triangle normals and areas
    MatrixX3f tmp = tri_nn.cwiseProduct(tri_nn);
    VectorXf normSize = tmp.rowwise().sum();
    normSize = normSize.cwiseSqrt();
 
    for (qint32 i = 0; i < normSize.size(); ++i)
        if (normSize(i) != 0)
            tri_nn.row(i) /= normSize(i);
 
    MatrixX3f nn = MatrixX3f::Zero(rr.rows(), 3);
 
    for (qint32 p = 0; p < tris.rows(); ++p) {
        Vector3i verts = tris.row(p);
        for (qint32 j = 0; j < verts.size(); ++j)
            nn.row(verts(j)) += tri_nn.row(p);
    }
 
    tmp = nn.cwiseProduct(nn);
    normSize = tmp.rowwise().sum();
    normSize = normSize.cwiseSqrt();
 
    for (qint32 i = 0; i < normSize.size(); ++i)
        if (normSize(i) != 0)
            nn.row(i) /= normSize(i);
 
    return nn;
}
 
//=============================================================================================================
 
bool FsSurface::read(const QString& subject_id, qint32 hemi, const QString& surf, const QString& subjects_dir, FsSurface& p_Surface, bool p_bLoadCurvature)
{
    if (hemi != 0 && hemi != 1)
        return false;
 
    QString p_sFile = QString("%1/%2/surf/%3.%4").arg(subjects_dir).arg(subject_id).arg(hemi == 0 ? "lh" : "rh").arg(surf);
 
    return read(p_sFile, p_Surface, p_bLoadCurvature);
}
 
//=============================================================================================================
 
bool FsSurface::read(const QString& path, qint32 hemi, const QString& surf, FsSurface& p_Surface, bool p_bLoadCurvature)
{
    if (hemi != 0 && hemi != 1)
        return false;
 
    QString p_sFile = QString("%1/%2.%3").arg(path).arg(hemi == 0 ? "lh" : "rh").arg(surf);
 
    return read(p_sFile, p_Surface, p_bLoadCurvature);
}
 
//=============================================================================================================
 
bool FsSurface::read(const QString& p_sFile, FsSurface& p_Surface, bool p_bLoadCurvature)
{
    p_Surface.clear();
 
    QFile t_File(p_sFile);
 
    if (!t_File.open(QIODevice::ReadOnly)) {
        qWarning("\tError: Couldn't open the surface file");
        return false;
    }
 
    qInfo("Reading surface...\n");
 
    //Strip file name and path
    qint32 t_NameIdx = 0;
    if (p_sFile.contains("lh."))
        t_NameIdx = p_sFile.indexOf("lh.");
    else if (p_sFile.contains("rh."))
        t_NameIdx = p_sFile.indexOf("rh.");
    else
        return false;
 
    p_Surface.m_sFilePath = p_sFile.mid(0, t_NameIdx);
    p_Surface.m_sFileName = p_sFile.mid(t_NameIdx, p_sFile.size() - t_NameIdx);
 
    QDataStream t_DataStream(&t_File);
    t_DataStream.setByteOrder(QDataStream::BigEndian);
 
    //
    //   Magic numbers to identify QUAD and TRIANGLE files
    //
    //   QUAD_FILE_MAGIC_NUMBER =  (-1 & 0x00ffffff) ;
    //   NEW_QUAD_FILE_MAGIC_NUMBER =  (-3 & 0x00ffffff) ;
    //
    qint32 NEW_QUAD_FILE_MAGIC_NUMBER = 16777213;
    qint32 TRIANGLE_FILE_MAGIC_NUMBER = 16777214;
    qint32 QUAD_FILE_MAGIC_NUMBER = 16777215;
 
    qint32 magic = FsSurface::fread3(t_DataStream);
 
    qint32 nvert = 0;
    qint32 nquad = 0;
    qint32 nface = 0;
    MatrixXf verts;
    MatrixXi faces;
 
    if (magic == QUAD_FILE_MAGIC_NUMBER || magic == NEW_QUAD_FILE_MAGIC_NUMBER) {
        nvert = FsSurface::fread3(t_DataStream);
        nquad = FsSurface::fread3(t_DataStream);
        if (magic == QUAD_FILE_MAGIC_NUMBER)
            qInfo("\t%s is a quad file (nvert = %d nquad = %d)\n", p_sFile.toUtf8().constData(), nvert, nquad);
        else
            qInfo("\t%s is a new quad file (nvert = %d nquad = %d)\n", p_sFile.toUtf8().constData(), nvert, nquad);
 
        // vertices, one column per vertex like the triangle branch (transposed below)
        verts.resize(3, nvert);
        if (magic == QUAD_FILE_MAGIC_NUMBER) {
            qint16 iVal;
            for (qint32 i = 0; i < nvert; ++i) {
                for (qint32 j = 0; j < 3; ++j) {
                    t_DataStream >> iVal;
                    verts(j, i) = static_cast<float>(iVal) / 100;
                }
            }
        } else {
            t_DataStream.readRawData(reinterpret_cast<char*>(verts.data()), nvert * 3 * sizeof(float));
            for (qint32 i = 0; i < nvert; ++i)
                for (qint32 j = 0; j < 3; ++j)
                    FIFFLIB::swap_floatp(&verts(j, i));
        }
 
        MatrixXi quads = FsSurface::fread3_many(t_DataStream, nquad * 4);
        MatrixXi quads_new(4, nquad);
        qint32 count = 0;
        for (qint32 j = 0; j < quads.cols(); ++j) {
            for (qint32 i = 0; i < quads.rows(); ++i) {
                quads_new(i, j) = quads(count, 0);
                ++count;
            }
        }
        quads = quads_new.transpose();
        //
        //  Face splitting follows
        //
        faces = MatrixXi::Zero(2 * nquad, 3);
        for (qint32 k = 0; k < nquad; ++k) {
            RowVectorXi quad = quads.row(k);
            if ((quad[0] % 2) == 0) {
                faces(nface, 0) = quad[0];
                faces(nface, 1) = quad[1];
                faces(nface, 2) = quad[3];
                ++nface;
 
                faces(nface, 0) = quad[2];
                faces(nface, 1) = quad[3];
                faces(nface, 2) = quad[1];
                ++nface;
            } else {
                faces(nface, 0) = quad(0);
                faces(nface, 1) = quad(1);
                faces(nface, 2) = quad(2);
                ++nface;
 
                faces(nface, 0) = quad(0);
                faces(nface, 1) = quad(2);
                faces(nface, 2) = quad(3);
                ++nface;
            }
        }
    } else if (magic == TRIANGLE_FILE_MAGIC_NUMBER) {
        QString s = t_File.readLine();
        t_File.readLine();
 
        t_DataStream >> nvert;
        t_DataStream >> nface;
 
        qInfo("\t%s is a triangle file (nvert = %d ntri = %d)\n", p_sFile.toUtf8().constData(), nvert, nface);
        qInfo("\t%s", s.toUtf8().constData());
 
        //vertices
        verts.resize(3, nvert);
        t_DataStream.readRawData(reinterpret_cast<char*>(verts.data()), nvert * 3 * sizeof(float));
        for (qint32 i = 0; i < 3; ++i)
            for (qint32 j = 0; j < nvert; ++j)
                FIFFLIB::swap_floatp(&verts(i, j));
 
        //faces
        faces.resize(nface, 3);
        qint32 iVal;
        for (qint32 i = 0; i < nface; ++i) {
            for (qint32 j = 0; j < 3; ++j) {
                t_DataStream >> iVal;
                faces(i, j) = iVal;
            }
        }
    } else {
        qWarning("Bad magic number (%d) in surface file %s", magic, p_sFile.toUtf8().constData());
        return false;
    }
 
    verts.transposeInPlace();
    verts.array() *= 0.001f;
 
    p_Surface.m_matRR = verts.block(0, 0, verts.rows(), 3);
    p_Surface.m_matTris = faces.block(0, 0, faces.rows(), 3);
 
    //-> not needed since qglbuilder is doing that for us
    p_Surface.m_matNN = compute_normals(p_Surface.m_matRR, p_Surface.m_matTris);
 
    // hemi info
    if (t_File.fileName().contains("lh."))
        p_Surface.m_iHemi = 0;
    else if (t_File.fileName().contains("rh."))
        p_Surface.m_iHemi = 1;
    else {
        p_Surface.m_iHemi = -1;
        return false;
    }
 
    //Loaded surface
    p_Surface.m_sSurf = t_File.fileName().mid((t_NameIdx + 3), t_File.fileName().size() - (t_NameIdx + 3));
 
    //Load curvature
    if (p_bLoadCurvature) {
        QString t_sCurvatureFile = QString("%1%2.curv").arg(p_Surface.m_sFilePath).arg(p_Surface.m_iHemi == 0 ? "lh" : "rh");
        qInfo("\t");
        p_Surface.m_vecCurv = FsSurface::read_curv(t_sCurvatureFile);
    }
 
    t_File.close();
    qInfo("\tRead a surface with %d vertices from %s\n[done]\n", nvert, p_sFile.toUtf8().constData());
 
    return true;
}
 
//=============================================================================================================
 
VectorXf FsSurface::read_curv(const QString& p_sFileName)
{
    VectorXf curv;
 
    qInfo("Reading curvature...");
    QFile t_File(p_sFileName);
 
    if (!t_File.open(QIODevice::ReadOnly)) {
        qWarning("\tError: Couldn't open the curvature file");
        return curv;
    }
 
    QDataStream t_DataStream(&t_File);
    t_DataStream.setByteOrder(QDataStream::BigEndian);
 
    qint32 vnum = FsSurface::fread3(t_DataStream);
    qint32 NEW_VERSION_MAGIC_NUMBER = 16777215;
 
    if (vnum == NEW_VERSION_MAGIC_NUMBER) {
        qint32 fnum, vals_per_vertex;
        t_DataStream >> vnum;
 
        t_DataStream >> fnum;
        t_DataStream >> vals_per_vertex;
 
        curv.resize(vnum, 1);
        t_DataStream.readRawData(reinterpret_cast<char*>(curv.data()), vnum * sizeof(float));
        for (qint32 i = 0; i < vnum; ++i)
            FIFFLIB::swap_floatp(&curv(i));
    } else {
        qint32 fnum = FsSurface::fread3(t_DataStream);
        Q_UNUSED(fnum)
        qint16 iVal;
        curv.resize(vnum, 1);
        for (qint32 i = 0; i < vnum; ++i) {
            t_DataStream >> iVal;
            curv(i) = static_cast<float>(iVal) / 100;
        }
    }
    t_File.close();
 
    qInfo("[done]\n");
 
    return curv;
}
 
//=============================================================================================================
 
qint32 FsSurface::fread3(QDataStream& stream)
{
    char bytes[3];
    stream.readRawData(bytes, 3);
    return (static_cast<unsigned char>(bytes[0]) << 16) + (static_cast<unsigned char>(bytes[1]) << 8) + static_cast<unsigned char>(bytes[2]);
}
 
//=============================================================================================================
 
qint32 FsSurface::fread3(std::iostream& stream)
{
    char bytes[3];
    stream.read(bytes, 3);
    return (static_cast<unsigned char>(bytes[0]) << 16) + (static_cast<unsigned char>(bytes[1]) << 8) + static_cast<unsigned char>(bytes[2]);
}
 
//=============================================================================================================
 
VectorXi FsSurface::fread3_many(QDataStream& stream, qint32 count)
{
    VectorXi res(count);
    for (qint32 i = 0; i < count; ++i)
        res[i] = FsSurface::fread3(stream);
    return res;
}
 
//=============================================================================================================
 
VectorXi FsSurface::fread3_many(std::iostream& stream, qint32 count)
{
    VectorXi res(count);
    for (qint32 i = 0; i < count; ++i)
        res[i] = FsSurface::fread3(stream);
    return res;
}