mne_msh_display_surface.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "mne_msh_display_surface.h"
#include "mne_surface.h"
#include "mne_morph_map.h"
#include "mne_msh_picked.h"
#include "mne_msh_color_scale_def.h"
 
#include <fiff/fiff_stream.h>
#include <fiff/fiff_coord_trans.h>
#include <fiff/fiff_digitizer_data.h>
#include <fiff/fiff_dig_point.h>
 
#include <math/sphere.h>
 
#include <Eigen/Geometry>
 
// MSVC builds already define _USE_MATH_DEFINES globally (see src/CMakeLists.txt),
// so define it here only for the toolchains that do not.
#ifndef _USE_MATH_DEFINES
#define _USE_MATH_DEFINES
#endif
#include <math.h>
#include <QDebug>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace Eigen;
using namespace FIFFLIB;
using namespace MNELIB;
 
//============================= local macros =================================
 
constexpr int FAIL = -1;
constexpr int OK = 0;
 
// Axis indices for coordinate access (kept for documentation).
[[maybe_unused]] constexpr int X = 0;
[[maybe_unused]] constexpr int Y = 1;
[[maybe_unused]] constexpr int Z = 2;
 
constexpr int SHOW_CURVATURE_NONE = 0;
constexpr int SHOW_CURVATURE_OVERLAY = 1;
constexpr int SHOW_OVERLAY_HEAT = 1;
 
constexpr float POS_CURV_COLOR = 0.25f;
constexpr float NEG_CURV_COLOR = 0.375f;
constexpr float EVEN_CURV_COLOR = 0.375f;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
MNEMshDisplaySurface::MNEMshDisplaySurface() = default;
 
//=============================================================================================================
 
MNEMshDisplaySurface::~MNEMshDisplaySurface()
{
    if (user_data_free)
        user_data_free(user_data);
}
 
//=============================================================================================================
// Alignment functions (moved from MNESurfaceOrVolume)
//=============================================================================================================
 
int MNEMshDisplaySurface::align_fiducials(FiffDigitizerData& head_dig,
                                          const FiffDigitizerData& mri_dig,
                                          int niter,
                                          int scale_head,
                                          float omit_dist,
                                          Eigen::Vector3f& scales)
 
{
    using FidMatrix = Eigen::Matrix<float, 3, 3, Eigen::RowMajor>;
    FidMatrix head_fid, mri_fid;
    bool head_fid_found[3] = {false, false, false};
    bool mri_fid_found[3] = {false, false, false};
    int j, k;
    FiffDigPoint p;
    float nasion_weight = 5.0;
 
    for (j = 0; j < 2; j++) {
        const FiffDigitizerData& d = (j == 0) ? head_dig : mri_dig;
        FidMatrix& fid = (j == 0) ? head_fid : mri_fid;
        bool* found = (j == 0) ? head_fid_found : mri_fid_found;
 
        for (k = 0; k < d.npoint; k++) {
            p = d.points[k];
            if (p.kind == FIFFV_POINT_CARDINAL) {
                if (p.ident == FIFFV_POINT_LPA) {
                    fid.row(0) = Eigen::Map<const Eigen::RowVector3f>(d.points[k].r);
                    found[0] = true;
                } else if (p.ident == FIFFV_POINT_NASION) {
                    fid.row(1) = Eigen::Map<const Eigen::RowVector3f>(d.points[k].r);
                    found[1] = true;
                } else if (p.ident == FIFFV_POINT_RPA) {
                    fid.row(2) = Eigen::Map<const Eigen::RowVector3f>(d.points[k].r);
                    found[2] = true;
                }
            }
        }
    }
 
    for (k = 0; k < 3; k++) {
        if (!head_fid_found[k]) {
            qCritical("Some of the MEG fiducials were missing");
            return FAIL;
        }
 
        if (!mri_fid_found[k]) {
            qCritical("Some of the MRI fiducials were missing");
            return FAIL;
        }
    }
 
    if (scale_head) {
        get_head_scale(head_dig, mri_fid, scales);
        qInfo("xscale = %.3f yscale = %.3f zscale = %.3f\n", scales[0], scales[1], scales[2]);
 
        for (j = 0; j < 3; j++)
            for (k = 0; k < 3; k++)
                mri_fid(j, k) = mri_fid(j, k) * scales[k];
 
        scale(scales);
    }
 
    // Initial alignment
    head_dig.head_mri_t_adj = std::make_unique<FiffCoordTrans>(FIFFLIB::FiffCoordTrans::fromCardinalPoints(FIFFV_COORD_HEAD, FIFFV_COORD_MRI,
                                                                                                           mri_fid.row(0).data(), mri_fid.row(1).data(), mri_fid.row(2).data()));
 
    // Populate mri_fids from cardinal digitizer points transformed into MRI coords
    head_dig.pickCardinalFiducials();
 
    // Overwrite the fiducial locations with the ones from the MRI digitizer data
    for (k = 0; k < head_dig.nfids(); k++)
        Eigen::Map<Eigen::Vector3f>(head_dig.mri_fids[k].r) = mri_fid.row(k).transpose();
    head_dig.head_mri_t_adj->print();
    qInfo("After simple alignment : \n");
 
    if (omit_dist > 0)
        discard_outlier_digitizer_points(head_dig, omit_dist);
 
    // Optional iterative refinement
    if (niter > 0) {
        for (k = 0; k < niter; k++) {
            if (iterate_alignment_once(head_dig, nasion_weight, Eigen::Vector3f(mri_fid.row(1).transpose()), k == niter - 1 && niter > 1) == FAIL)
                return FAIL;
        }
 
        qInfo("%d / %d iterations done. RMS dist = %7.1f mm\n", k, niter,
              1000.0 * rms_digitizer_distance(head_dig));
        qInfo("After refinement :\n");
        head_dig.head_mri_t_adj->print();
    }
 
    return OK;
}
 
//=============================================================================================================
 
// Simple head size fit
void MNEMshDisplaySurface::get_head_scale(FIFFLIB::FiffDigitizerData& dig,
                                          const Eigen::Matrix<float, 3, 3, Eigen::RowMajor>& mri_fid,
                                          Eigen::Vector3f& scales)
{
    int k, ndig, nhead;
    float simplex_size = 2e-2f;
    Eigen::VectorXf r0(3);
    float Rdig, Rscalp;
 
    scales[0] = scales[1] = scales[2] = 1.0;
 
    Eigen::MatrixXf dig_rr(dig.npoint, 3);
    Eigen::MatrixXf head_rr(np, 3);
 
    // Pick only the points with positive z
    for (k = 0, ndig = 0; k < dig.npoint; k++) {
        if (dig.points[k].r[2] > 0) {
            dig_rr.row(ndig++) = Eigen::Map<const Eigen::RowVector3f>(dig.points[k].r);
        }
    }
 
    if (!UTILSLIB::Sphere::fit_sphere_to_points(dig_rr.topRows(ndig), simplex_size, r0, Rdig)) {
        return;
    }
 
    qInfo("Polhemus : (%.1f %.1f %.1f) mm R = %.1f mm\n", 1000 * r0[0], 1000 * r0[1], 1000 * r0[2], 1000 * Rdig);
 
    // Pick only the points above the fiducial plane
    Eigen::Vector3f LR = mri_fid.row(2).transpose() - mri_fid.row(0).transpose();
    Eigen::Vector3f LN = mri_fid.row(1).transpose() - mri_fid.row(0).transpose();
    Eigen::Vector3f norm = LR.cross(LN);
    norm.normalize();
 
    for (k = 0, nhead = 0; k < np; k++) {
        Eigen::Vector3f diff_vec = rr.row(k).transpose() - mri_fid.row(0).transpose();
        if (diff_vec.dot(norm) > 0) {
            head_rr.row(nhead++) = rr.row(k);
        }
    }
 
    if (!UTILSLIB::Sphere::fit_sphere_to_points(head_rr.topRows(nhead), simplex_size, r0, Rscalp)) {
        return;
    }
 
    qInfo("Scalp : (%.1f %.1f %.1f) mm R = %.1f mm\n", 1000 * r0[0], 1000 * r0[1], 1000 * r0[2], 1000 * Rscalp);
 
    scales[0] = scales[1] = scales[2] = Rdig / Rscalp;
}
 
//=============================================================================================================
 
int MNEMshDisplaySurface::discard_outlier_digitizer_points(FIFFLIB::FiffDigitizerData& d,
                                                           float maxdist) const
/*
      * Discard outlier digitizer points
      */
{
    int discarded = 0;
    int k;
 
    d.dist_valid = false;
    calculate_digitizer_distances(d, true, true);
    for (k = 0; k < d.npoint; k++) {
        d.discard[k] = 0;
        /*
        * Discard unless cardinal landmark or HPI coil
        */
        if (std::fabs(d.dist(k)) > maxdist &&
            d.points[k].kind != FIFFV_POINT_CARDINAL &&
            d.points[k].kind != FIFFV_POINT_HPI) {
            discarded++;
            d.discard[k] = 1;
        }
    }
    qInfo("%d points discarded (maxdist = %6.1f mm).\n", discarded, 1000 * maxdist);
 
    return discarded;
}
 
//=============================================================================================================
 
void MNEMshDisplaySurface::calculate_digitizer_distances(FIFFLIB::FiffDigitizerData& dig,
                                                         bool do_all, bool do_approx) const
/*
 * Calculate the distances from the scalp surface
 */
{
    int k, nactive;
    FiffDigPoint point;
    Q_ASSERT(dig.head_mri_t);
    const FiffCoordTrans& t = (dig.head_mri_t_adj && !dig.head_mri_t_adj->isEmpty()) ? *dig.head_mri_t_adj : *dig.head_mri_t;
    int nstep = 4;
 
    if (dig.dist_valid)
        return;
 
    PointsT digPoints(dig.npoint, 3);
 
    dig.dist.conservativeResize(dig.npoint);
    if (dig.closest.size() == 0) {
        /*
        * Ensure that all closest values are initialized correctly
        */
        dig.closest = Eigen::VectorXi::Constant(dig.npoint, -1);
    }
 
    dig.closest_point.setZero(dig.npoint, 3);
    Eigen::VectorXi closest(dig.npoint);
    Eigen::VectorXf dists(dig.npoint);
 
    for (k = 0, nactive = 0; k < dig.npoint; k++) {
        if ((dig.active[k] && !dig.discard[k]) || do_all) {
            point = dig.points.at(k);
            digPoints.row(nactive) = Eigen::Map<const Eigen::RowVector3f>(point.r);
            FiffCoordTrans::apply_trans(digPoints.row(nactive).data(), t, FIFFV_MOVE);
            if (do_approx) {
                closest[nactive] = dig.closest(k);
                if (closest[nactive] < 0)
                    do_approx = false;
            } else
                closest[nactive] = -1;
            nactive++;
        }
    }
 
    find_closest_on_surface_approx(digPoints, nactive, closest, dists, nstep);
    /*
     * Project the points on the triangles
     */
    if (!do_approx)
        qInfo("Inside or outside for %d points...", nactive);
    for (k = 0, nactive = 0; k < dig.npoint; k++) {
        if ((dig.active[k] && !dig.discard[k]) || do_all) {
            dig.dist(k) = dists[nactive];
            dig.closest(k) = closest[nactive];
            {
                Eigen::Vector3f pt = Eigen::Map<const Eigen::Vector3f>(digPoints.row(nactive).data());
                Eigen::Vector3f proj = project_to_triangle(dig.closest(k), pt);
                dig.closest_point.row(k) = proj.transpose();
            }
            // The distance above is with respect to the closest triangle only,
            // so its sign says which side of that one triangle the point is on
            // rather than whether the point is inside the surface. The solid
            // angle criterion below answers the latter, and is the check MNE-C
            // describes here.
            //
            // It is off on purpose, and not because it is wrong or unfinished:
            //
            //  - The criterion is correct. sum_solids adds 2*atan2 per
            //    triangle, so a closed surface gives 4*pi for a point inside
            //    and 0 for one outside, which makes the 0.9 threshold exact.
            //    Verified on a closed cube and on a 9 cm sphere, where it
            //    still separates cleanly 2 mm either side of the surface.
            //    mne-python computes the same quantity for its inside/outside
            //    test, only normalised to 2*pi because its per triangle term
            //    is -arctan2 rather than 2*atan2.
            //
            //  - Nothing reads the sign. Both consumers of dig.dist discard
            //    it: discard_outlier_digitizer_points compares
            //    std::fabs(dist) against maxdist, and rms_digitizer_distance
            //    squares it.
            //
            // Enabling it therefore costs a full pass over every triangle for
            // each digitizer point, thousands of atan2 calls per point, and
            // changes no observable behaviour. Flip the constant if a caller
            // ever needs a signed distance.
            constexpr bool bUseSolidAngleSign = false;
            if (bUseSolidAngleSign && !do_approx) {
                Eigen::Vector3f pt = Eigen::Map<const Eigen::Vector3f>(digPoints.row(nactive).data());
                if (sum_solids(pt) / (4 * M_PI) > 0.9)
                    dig.dist(k) = -std::fabs(dig.dist(k));
                else
                    dig.dist(k) = std::fabs(dig.dist(k));
            }
            nactive++;
        }
    }
 
    if (!do_approx)
        qInfo("[done]\n");
 
    dig.dist_valid = true;
 
    return;
}
 
//=============================================================================================================
 
int MNEMshDisplaySurface::iterate_alignment_once(FIFFLIB::FiffDigitizerData& dig,                  /* The digitizer data */
                                                 int nasion_weight,                                /* Weight for the nasion */
                                                 const std::optional<Eigen::Vector3f>& nasion_mri, /* Fixed correspondence point for the nasion (optional) */
                                                 int last_step) const                              /* Is this the last iteration step */
/*
 * Find the best alignment of the coordinate frames
 */
{
    Eigen::MatrixXf rr_head(dig.npoint, 3);
    Eigen::MatrixXf rr_mri(dig.npoint, 3);
    Eigen::VectorXf w = Eigen::VectorXf::Zero(dig.npoint);
    int k, nactive;
    FiffDigPoint point;
    FiffCoordTrans t;
    float max_diff = 40e-3f;
 
    if (!dig.head_mri_t_adj) {
        qCritical() << "Not adjusting the transformation";
        return FAIL;
    }
    /*
     * Calculate initial distances
     */
    calculate_digitizer_distances(dig, false, true);
 
    for (k = 0, nactive = 0; k < dig.npoint; k++) {
        if (dig.active[k] && !dig.discard[k]) {
            point = dig.points.at(k);
            rr_head.row(nactive) = Eigen::Map<const Eigen::RowVector3f>(point.r);
            rr_mri.row(nactive) = dig.closest_point.row(k);
            /*
            * Special handling for the nasion
            */
            if (point.kind == FIFFV_POINT_CARDINAL &&
                point.ident == FIFFV_POINT_NASION) {
                w[nactive] = nasion_weight;
                if (nasion_mri) {
                    rr_mri.row(nactive) = nasion_mri->transpose();
                    // rr_head is column-major, so its rows are not contiguous; transform a copy.
                    Eigen::Vector3f nasion_head = *nasion_mri;
                    Q_ASSERT(dig.head_mri_t || dig.head_mri_t_adj);
                    FiffCoordTrans::apply_inverse_trans(nasion_head.data(),
                                                        dig.head_mri_t_adj ? *dig.head_mri_t_adj : *dig.head_mri_t,
                                                        FIFFV_MOVE);
                    rr_head.row(nactive) = nasion_head.transpose();
                }
            } else
                w[nactive] = 1.0;
            nactive++;
        }
    }
    if (nactive < 3) {
        qCritical() << "Not enough points to do the alignment";
        return FAIL;
    }
    if ((t = FiffCoordTrans::procrustesAlign(FIFFV_COORD_HEAD, FIFFV_COORD_MRI,
                                             rr_head.topRows(nactive),
                                             rr_mri.topRows(nactive),
                                             w.head(nactive),
                                             max_diff))
            .isEmpty())
        return FAIL;
 
    if (dig.head_mri_t_adj)
        dig.head_mri_t_adj = std::make_unique<FiffCoordTrans>(t);
    /*
     * Calculate final distances
     */
    dig.dist_valid = false;
    calculate_digitizer_distances(dig, false, !last_step);
    return OK;
}
 
//=============================================================================================================
 
float MNEMshDisplaySurface::rms_digitizer_distance(FIFFLIB::FiffDigitizerData& dig) const
{
    float rms;
    int k, nactive;
 
    calculate_digitizer_distances(dig, false, true);
 
    for (k = 0, rms = 0.0, nactive = 0; k < dig.npoint; k++)
        if (dig.active[k] && !dig.discard[k]) {
            rms = rms + dig.dist(k) * dig.dist(k);
            nactive++;
        }
    if (nactive > 1)
        rms = rms / (nactive - 1);
    return sqrt(rms);
}
 
//=============================================================================================================
 
void MNEMshDisplaySurface::scale(const Eigen::Vector3f& scales)
/*
 * Not quite complete yet
 */
{
    int j, k;
 
    for (k = 0; k < 3; k++) {
        minv[k] = scales[k] * minv[k];
        maxv[k] = scales[k] * maxv[k];
    }
    for (j = 0; j < np; j++)
        for (k = 0; k < 3; k++)
            rr(j, k) = rr(j, k) * scales[k];
    return;
}
 
//=============================================================================================================
 
void MNEMshDisplaySurface::decide_surface_extent([[maybe_unused]] const QString& tag)
{
    Eigen::Vector3f mn = rr.row(0).transpose();
    Eigen::Vector3f mx = mn;
 
    for (int k = 1; k < np; k++) {
        Eigen::Vector3f r = rr.row(k).transpose();
        mn = mn.cwiseMin(r);
        mx = mx.cwiseMax(r);
    }
 
#ifdef DEBUG
    qInfo("%s:\n", tag.toUtf8().constData());
    qInfo("\tx = %f ... %f mm\n", 1000 * mn[0], 1000 * mx[0]);
    qInfo("\ty = %f ... %f mm\n", 1000 * mn[1], 1000 * mx[1]);
    qInfo("\tz = %f ... %f mm\n", 1000 * mn[2], 1000 * mx[2]);
#endif
 
    fov = std::max(mn.cwiseAbs().maxCoeff(), mx.cwiseAbs().maxCoeff());
 
    minv = mn;
    maxv = mx;
    fov_scale = 1.1f;
}
 
//=============================================================================================================
 
void MNEMshDisplaySurface::decide_curv_display(const QString& name)
{
    if (name.startsWith("inflated") || name.startsWith("sphere") || name.startsWith("white"))
        curvature_color_mode = SHOW_CURVATURE_OVERLAY;
    else
        curvature_color_mode = SHOW_CURVATURE_NONE;
    overlay_color_mode = SHOW_OVERLAY_HEAT;
}
 
//=============================================================================================================
 
void MNEMshDisplaySurface::setup_curvature_colors()
{
    if (np == 0)
        return;
 
    const int ncolor = nvertex_colors;
    const int totalSize = ncolor * np;
 
    if (vertex_colors.size() == 0)
        vertex_colors.resize(totalSize);
 
    [[maybe_unused]] float curv_sum = 0.0f;
    if (curvature_color_mode == SHOW_CURVATURE_OVERLAY) {
        for (int k = 0; k < np; k++) {
            const int base = k * ncolor;
            curv_sum += std::fabs(curv[k]);
            const float c = (curv[k] > 0) ? POS_CURV_COLOR : NEG_CURV_COLOR;
            for (int j = 0; j < 3; j++)
                vertex_colors[base + j] = c;
            if (ncolor == 4)
                vertex_colors[base + 3] = 1.0f;
        }
    } else {
        for (int k = 0; k < np; k++) {
            const int base = k * ncolor;
            curv_sum += std::fabs(curv[k]);
            for (int j = 0; j < 3; j++)
                vertex_colors[base + j] = EVEN_CURV_COLOR;
            if (ncolor == 4)
                vertex_colors[base + 3] = 1.0f;
        }
    }
#ifdef DEBUG
    qInfo("Average curvature : %f\n", curv_sum / np);
#endif
}