MneICP

Namespace: MNELIB · Library: MNE Library

Module

This page documents a header-level module — a collection of free functions that share an algorithmic topic. There is no enclosing C++ class; the functions live directly in the library namespace.

Python equivalent

mne.coreg.Coregistration in MNE-Python.

#include <mne/mne_icp.h>

MNELIB::MNEIcp wraps the classical point-to-surface ICP used by mne-cpp's coregistration tools to refine an initial head <-> MRI transformation. It projects each digitiser point onto the nearest triangle of an MNEBemSurface (via MNEProjectToSurface), solves a rigid-body least squares with SVD and iterates until convergence.

Functions

performIcp(mneSurfacePoints, matPointCloud, transFromTo, fRMSE, bScale, iMaxIter, fTol, vecWeights)

bool performIcp(const QSharedPointer<[MNEProjectToSurface](/docs/api/mne/mne-project-to-surface)> mneSurfacePoints, const Eigen::MatrixXf & matPointCloud, [FiffCoordTrans](/docs/api/fiff/fiff-coord-trans) & transFromTo, float & fRMSE, bool bScale, int iMaxIter, float fTol, const Eigen::VectorXf & vecWeights);

The ICP algorithm to register a point cloud with a surface.

Parameters:

mneSurfacePoints : const QSharedPointer<MNEProjectToSurface> The MNEProjectToSurface object that contains the surface triangles etc. (To).
matPointCloud : const Eigen::MatrixXf & The point cloud to be registered (From).
transFromTo : FiffCoordTrans & The forward transformation matrix. It can contain an initial transformation (e.g. from fiducial alignment).
fRMSE : float & The resulting Root-Mean-Square-Error in m.
bScale : bool Whether to apply scaling or not. Should be false for matching data sets, defaults to false.
iMaxIter : int The maximum number of iterations for the ICP algorithm, defaults to 20.
fTol : float The convergence tolerance, defaults to 0.001.
vecWeights : const Eigen::VectorXf & The weights to apply, defaults to zeros.

Returns:

bool — Whether the registration was successful.

fitMatchedPoints(matSrcPoint, matDstPoint, matTrans, fScale, bScale, vecWeights)

bool fitMatchedPoints(const Eigen::MatrixXf & matSrcPoint, const Eigen::MatrixXf & matDstPoint, Eigen::Matrix4f & matTrans, float fScale, bool bScale, const Eigen::VectorXf & vecWeights);

Corresponding point set registration using quaternions.

Parameters:

matSrcPoint : const Eigen::MatrixXf & The source point set.
matDstPoint : const Eigen::MatrixXf & The destination point set.
matTrans : Eigen::Matrix4f & The forward transformation matrix.
fScale : float The scaling parameter, defaults to 1.0.
bScale : bool Whether to apply scaling or not. Should be false for matching data sets, defaults to false.
vecWeights : const Eigen::VectorXf & The weights to apply, defaults to zeros.

Returns:

bool — Whether the matching was successful.

discard3DPointOutliers(mneSurfacePoints, matPointCloud, transFromTo, vecTake, matTakePoint, fMaxDist)

bool discard3DPointOutliers(const QSharedPointer<[MNEProjectToSurface](/docs/api/mne/mne-project-to-surface)> mneSurfacePoints, const Eigen::MatrixXf & matPointCloud, const [FiffCoordTrans](/docs/api/fiff/fiff-coord-trans) & transFromTo, Eigen::VectorXi & vecTake, Eigen::MatrixXf & matTakePoint, float fMaxDist);

Discard outliers compared to a given 3D surface.

Parameters:

mneSurfacePoints : const QSharedPointer<MNEProjectToSurface> The MNEProjectToSurface object that contains the surface triangles etc. (To).
matPointCloud : const Eigen::MatrixXf & The point cloud to be registered (From).
transFromTo : const FiffCoordTrans & The forward transformation matrix.
vecTake : Eigen::VectorXi & The index of taken digitizers.
matTakePoint : Eigen::MatrixXf & The digitizer points to take.
fMaxDist : float The maximum distance to the surface in mm, defaults to 0 mm.

Returns:

bool — Whether the discarding was successful.

Example

Source: src/examples/ex_coreg/main.cpp

#include <iostream>

#include <disp3D/view/brainview.h>
#include <disp3D/model/braintreemodel.h>
#include <disp3D/model/items/digitizersettreeitem.h>
#include <disp3D/model/items/bemtreeitem.h>

#include "fiff/fiff_dig_point_set.h"
#include "fiff/fiff_dig_point.h"
#include "fiff/fiff_coord_trans.h"

#include "mne/mne_bem.h"
#include "mne/mne_bem_surface.h"
#include "mne/mne_project_to_surface.h"

#include <utils/generics/mne_logger.h>
#include "mne/mne_icp.h"

//=============================================================================================================
// QT INCLUDES
//=============================================================================================================

#include <QApplication>
#include <QMainWindow>
#include <QCommandLineParser>
#include <QDebug>
#include <QFile>

//=============================================================================================================
// Eigen
//=============================================================================================================

#include <Eigen/Core>

//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================

using namespace Eigen;
using namespace UTILSLIB;
using namespace FIFFLIB;
using namespace MNELIB;

//=============================================================================================================
// MAIN
//=============================================================================================================

//=============================================================================================================
/**
 * The function main marks the entry point of the program.
 * By default, main has the storage class extern.
 *
 * @param[in] argc (argument count) is an integer that indicates how many arguments were entered on the command line when the program was started.
 * @param[in] argv (argument vector) is an array of pointers to arrays of character objects. The array objects are null-terminated strings, representing the arguments that were entered on the command line when the program was started.
 * @return the value that was set to exit() (which is 0 if exit() is called via quit()).
 */
int main(int argc, char *argv[])
{
    #ifdef STATICBUILD
    // Q_INIT_RESOURCE(mne_disp3d);
    #endif

    qInstallMessageHandler(MNELogger::customLogWriter);
    QApplication a(argc, argv);

    // Command Line Parser
    QCommandLineParser parser;
    parser.setApplicationDescription("Example Coregistration");
    parser.addHelpOption();

    QCommandLineOption digOption("dig", "The destination point set", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis-ave.fif");
    QCommandLineOption bemOption("bem", "The bem file", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/subjects/sample/bem/sample-head.fif");
    QCommandLineOption transOption("trans", "The MRI-Head transformation file", "file", QCoreApplication::applicationDirPath() + "/../resources/data/MNE-sample-data/MEG/sample/sample_audvis_raw-trans.fif");
    QCommandLineOption scaleOption("scale", "Weather to scale during the registration or not", "bool", "false");
    QCommandLineOption tolOption("tol", "The convergence limit for the icp algorithm.", "float", "0.001");
    QCommandLineOption distOption("dist", "The maximum distance between digitizer and head shape in mm.", "float", "0.02");
    QCommandLineOption iterOption("iter", "The maximum number of icp iterations.", "int", "20");

    parser.addOption(digOption);
    parser.addOption(bemOption);
    parser.addOption(scaleOption);
    parser.addOption(transOption);
    parser.addOption(tolOption);
    parser.addOption(distOption);
    parser.addOption(iterOption);

    parser.process(a);

    // get cli parameters
    QFile t_fileDig(parser.value(digOption));
    QFile t_fileBem(parser.value(bemOption));
    QFile t_fileTrans(parser.value(transOption));

    bool bScale = false;
    if(parser.value(scaleOption) == "false" || parser.value(scaleOption) == "0") {
        bScale = false;
    } else if(parser.value(scaleOption) == "true" || parser.value(scaleOption) == "1") {
        bScale = true;
    }

    float fTol = parser.value(tolOption).toFloat();
    float fMaxDist = parser.value(distOption).toFloat();
    int iMaxIter = parser.value(iterOption).toInt();

    // read Trans
    FiffCoordTrans transHeadMriRef(t_fileTrans);

    // read Bem
    MNEBem bemHead(t_fileBem);
    MNEBemSurface::SPtr bemSurface = MNEBemSurface::SPtr::create(bemHead[0]);
    MNEProjectToSurface::SPtr mneSurfacePoints = MNEProjectToSurface::SPtr::create(*bemSurface);

    // read digitizer data
    QList<int> lPickFiducials({FIFFV_POINT_CARDINAL});
    QList<int> lPickHSP({FIFFV_POINT_CARDINAL,FIFFV_POINT_HPI,FIFFV_POINT_EXTRA,FIFFV_POINT_EEG});
    FiffDigPointSet digSetSrc = FiffDigPointSet(t_fileDig).pickTypes(lPickFiducials);   // Fiducials Head-Space
    FiffDigPointSet digSetDst = FiffDigPointSet(t_fileDig).pickTypes(lPickFiducials);
    digSetDst.applyTransform(transHeadMriRef, false);
    FiffDigPointSet digSetHsp = FiffDigPointSet(t_fileDig).pickTypes(lPickHSP);         // Head shape points Head-Space

    // Initial Fiducial Alignment
    // Declare variables
    Matrix3f matSrc(digSetSrc.size(),3);
    Matrix3f matDst(digSetDst.size(),3);
    Matrix4f matTrans;
    Vector3f vecWeights(digSetSrc.size()); // LPA, Nasion, RPA
    float fScale;

    // get coordinates
    for(int i = 0; i< digSetSrc.size(); ++i) {
        matSrc(i,0) = digSetSrc[i].r[0]; matSrc(i,1) = digSetSrc[i].r[1]; matSrc(i,2) = digSetSrc[i].r[2];
        matDst(i,0) = digSetDst[i].r[0]; matDst(i,1) = digSetDst[i].r[1]; matDst(i,2) = digSetDst[i].r[2];

        // set standart weights
        if(digSetSrc[i].ident == FIFFV_POINT_NASION) {
            vecWeights(i) = 10.0;
        } else {
            vecWeights(i) = 1.0;
        }
    }

    // align fiducials
    if(!MNELIB::fitMatchedPoints(matSrc,matDst,matTrans,fScale,bScale,vecWeights)) {
        qWarning() << "Point cloud registration not succesfull.";
    }

    fiff_int_t iFrom = digSetSrc[0].coord_frame;
    fiff_int_t iTo = bemSurface->coord_frame;
    FiffCoordTrans transHeadMri(iFrom, iTo, matTrans);

    // Prepare Icp:
    VectorXf vecWeightsICP(digSetHsp.size()); // Weigths vector
    MatrixXf matHsp(digSetHsp.size(),3);

    for(int i = 0; i < digSetHsp.size(); ++i) {
        matHsp(i,0) = digSetHsp[i].r[0]; matHsp(i,1) = digSetHsp[i].r[1]; matHsp(i,2) = digSetHsp[i].r[2];
        // set standart weights
        if((digSetHsp[i].kind == FIFFV_POINT_CARDINAL) && (digSetHsp[i].ident == FIFFV_POINT_NASION)) {
            vecWeightsICP(i) = 10.0;
        } else {
            vecWeightsICP(i) = 1.0;
        }
    }

    MatrixXf matHspClean;
    VectorXi vecTake;
    float fRMSE = 0.0;

    // discard outliers
    if(!MNELIB::discard3DPointOutliers(mneSurfacePoints, matHsp, transHeadMri, vecTake, matHspClean, fMaxDist)) {
        qWarning() << "Discard outliers was not succesfull.";
    }
    VectorXf vecWeightsICPClean(vecTake.size());

    for(int i = 0; i < vecTake.size(); ++i) {
        vecWeightsICPClean(i) = vecWeightsICP(vecTake(i));
    }

    // icp
    if(!MNELIB::performIcp(mneSurfacePoints, matHspClean, transHeadMri, fRMSE, bScale, iMaxIter, fTol, vecWeightsICPClean)) {
        qWarning() << "ICP was not succesfull.";
    }
    qInfo() << "transHeadMri:";
    transHeadMri.print();
    qInfo() << "transHeadMriRef:";
    transHeadMriRef.print();

    BrainView *pBrainView = new BrainView();
    BrainTreeModel *pModel = new BrainTreeModel();
    pBrainView->setModel(pModel);
    pModel->addDigitizerData(digSetSrc.getList());
    pModel->addDigitizerData(digSetHsp.getList());

    // Add BEM surfaces
    for (int i = 0; i < bemHead.size(); ++i) {
        pModel->addBemSurface("Sample", "Head", bemHead[i]);
    }

    pBrainView->show();

    return a.exec();
}

Authors of this file

Christoph Dinh <christoph.dinh@mne-cpp.org>

Functions​

performIcp(mneSurfacePoints, matPointCloud, transFromTo, fRMSE, bScale, iMaxIter, fTol, vecWeights)​

fitMatchedPoints(matSrcPoint, matDstPoint, matTrans, fScale, bScale, vecWeights)​

discard3DPointOutliers(mneSurfacePoints, matPointCloud, transFromTo, vecTake, matTakePoint, fMaxDist)​

Example​

Authors of this file​

Functions

performIcp(mneSurfacePoints, matPointCloud, transFromTo, fRMSE, bScale, iMaxIter, fTol, vecWeights)

fitMatchedPoints(matSrcPoint, matDstPoint, matTrans, fScale, bScale, vecWeights)

discard3DPointOutliers(mneSurfacePoints, matPointCloud, transFromTo, vecTake, matTakePoint, fMaxDist)

Example

Authors of this file