inv_dics.cpp

Go to the documentation of this file.

//=============================================================================================================
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "inv_dics.h"
#include "inv_beamformer_compute.h"
 
#include <mne/mne_forward_solution.h>
#include <fiff/fiff_cov.h>
#include <fiff/fiff_info.h>
 
#include <QDebug>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Dense>
 
//=============================================================================================================
// USED NAMESPACES
//=============================================================================================================
 
using namespace Eigen;
using namespace INVLIB;
using namespace MNELIB;
using namespace FIFFLIB;
 
//=============================================================================================================
// DEFINE MEMBER METHODS
//=============================================================================================================
 
InvBeamformer InvDICS::makeDICS([[maybe_unused]] const FiffInfo& info,
                                const MNEForwardSolution& forward,
                                const std::vector<MatrixXd>& csdMatrices,
                                const VectorXd& frequencies,
                                double reg,
                                bool realFilter,
                                const FiffCov& noiseCov,
                                BeamformerPickOri pickOri,
                                BeamformerWeightNorm weightNorm,
                                bool reduceRank,
                                BeamformerInversion invMethod)
{
    InvBeamformer result;
    result.kind = "DICS";
 
    const int nFreqs = static_cast<int>(csdMatrices.size());
    if (nFreqs == 0) {
        qWarning("InvDICS::makeDICS - No CSD matrices provided!");
        return result;
    }
    if (frequencies.size() != nFreqs) {
        qWarning("InvDICS::makeDICS - Frequency vector size mismatch with CSD count!");
        return result;
    }
 
    // -----------------------------------------------------------------------
    // Extract leadfield
    // -----------------------------------------------------------------------
    if (!forward.sol || forward.sol->data.size() == 0) {
        qWarning("InvDICS::makeDICS - Forward solution has no gain matrix!");
        return result;
    }
 
    MatrixXd G = forward.sol->data;
    const int nChannels = static_cast<int>(G.rows());
    const int nOrient = (forward.source_ori == FIFFV_MNE_FREE_ORI) ? 3 : 1;
    const int nSources = static_cast<int>(G.cols()) / nOrient;
 
    qInfo("InvDICS::makeDICS - Leadfield: %d channels x %d sources (n_orient=%d), %d frequencies",
          nChannels, nSources, nOrient, nFreqs);
 
    // -----------------------------------------------------------------------
    // Whitening matrix
    // -----------------------------------------------------------------------
    MatrixXd whitener;
    if (noiseCov.data.size() > 0 && noiseCov.eig.size() > 0 && noiseCov.eigvec.size() > 0) {
        VectorXd invSqrtEig(noiseCov.eig.size());
        for (int i = 0; i < noiseCov.eig.size(); ++i) {
            invSqrtEig(i) = (noiseCov.eig(i) > 1e-30)
                ? 1.0 / std::sqrt(noiseCov.eig(i))
                : 0.0;
        }
        // Rows of FiffCov::eigvec are the eigenvectors.
        whitener = invSqrtEig.asDiagonal() * noiseCov.eigvec;
    } else {
        whitener = MatrixXd::Identity(nChannels, nChannels);
    }
 
    MatrixXd projMat = MatrixXd::Identity(nChannels, nChannels);
 
    // Whiten leadfield (shared across all frequencies)
    MatrixXd Gw = whitener * G;
 
    // Source normals
    MatrixX3d nn = forward.source_nn.cast<double>();
    if (nOrient == 3 && nn.rows() == 3 * nSources) {
        // Free orientation stores three rows per source; mne-python uses nn[2::3].
        MatrixX3d perSource(nSources, 3);
        for (int s = 0; s < nSources; ++s)
            perSource.row(s) = nn.row(3 * s + 2);
        nn = perSource;
    }
 
    // -----------------------------------------------------------------------
    // Compute filter for each frequency
    // -----------------------------------------------------------------------
    for (int fi = 0; fi < nFreqs; ++fi) {
        MatrixXd Cm = csdMatrices[fi];
 
        if (Cm.rows() != nChannels || Cm.cols() != nChannels) {
            qWarning("InvDICS::makeDICS - CSD[%d] dimension mismatch!", fi);
            return InvBeamformer();
        }
 
        // Optional: take real part of CSD
        if (realFilter) {
            // CSD is provided as real-valued after user extracts real part,
            // or we ensure it here
            Cm = Cm.real();
        }
 
        // Whiten CSD
        MatrixXd CmW = whitener * Cm * whitener.transpose();
        CmW = (CmW + CmW.transpose()) * 0.5; // Ensure symmetry
 
        // Compute filter
        MatrixXd W;
        MatrixX3d mpOri;
 
        bool ok = InvBeamformerCompute::computeBeamformer(
            Gw, CmW, reg, nOrient,
            weightNorm, pickOri, reduceRank, invMethod,
            nn, W, mpOri);
 
        if (!ok) {
            qWarning("InvDICS::makeDICS - Filter computation failed at frequency %d (%.1f Hz)!",
                     fi, frequencies(fi));
            return InvBeamformer();
        }
 
        result.weights.push_back(W);
 
        // Store max-power orientation from first frequency
        if (fi == 0 && pickOri == BeamformerPickOri::MaxPower) {
            result.maxPowerOri = mpOri;
        }
    }
 
    // -----------------------------------------------------------------------
    // Populate metadata
    // -----------------------------------------------------------------------
    result.whitener = whitener;
    result.proj = projMat;
    result.chNames = forward.sol->row_names;
    result.isFreOri = (nOrient == 3 && pickOri != BeamformerPickOri::Normal && pickOri != BeamformerPickOri::MaxPower);
    result.nSourcesTotal = nSources;
    result.srcType = "surface";
    result.weightNorm = weightNorm;
    result.pickOri = pickOri;
    result.inversion = invMethod;
    result.reg = reg;
    result.rank = static_cast<int>(nChannels);
    result.sourceNn = forward.source_nn;
    result.frequencies = frequencies;
 
    VectorXi verts(0);
    if (forward.src.size() >= 2) {
        verts.resize(forward.src[0].vertno.size() + forward.src[1].vertno.size());
        verts << forward.src[0].vertno, forward.src[1].vertno;
        // Keep where the left hemisphere ends. The concatenation above is the
        // only place that knows it, and without it the result cannot be
        // written as the -lh/-rh pair mne-python expects.
        result.nVerticesLh = static_cast<int>(forward.src[0].vertno.size());
    } else if (forward.src.size() == 1) {
        verts = forward.src[0].vertno;
    }
    result.vertices = verts;
 
    qInfo("InvDICS::makeDICS - Done. %d frequency filters computed.", nFreqs);
 
    return result;
}
 
//=============================================================================================================
 
InvSourceEstimate InvDICS::applyDICSCsd(const std::vector<MatrixXd>& csdMatrices,
                                        const VectorXd& frequencies,
                                        const InvBeamformer& filters)
{
    if (!filters.isValid() || filters.kind != "DICS") {
        qWarning("InvDICS::applyDICSCsd - Invalid or non-DICS filters!");
        return InvSourceEstimate();
    }
 
    const int nFreqs = static_cast<int>(csdMatrices.size());
    const int nFilterFreqs = filters.nFreqs();
 
    if (nFreqs != nFilterFreqs) {
        qWarning("InvDICS::applyDICSCsd - CSD count (%d) does not match filter count (%d)!",
                 nFreqs, nFilterFreqs);
        return InvSourceEstimate();
    }
 
    const int nOrient = filters.nOrient();
    const int nSources = filters.nSources();
 
 
    // Power matrix: (nSources, nFreqs)
    MatrixXd powerMat(nSources, nFreqs);
 
    for (int fi = 0; fi < nFreqs; ++fi) {
        MatrixXd Cm = csdMatrices[fi];
 
        // Whiten CSD
        if (filters.whitener.size() > 0) {
            Cm = filters.whitener * Cm * filters.whitener.transpose();
        }
 
        VectorXd power = InvBeamformerCompute::computePower(Cm, filters.weights[fi], nOrient);
        powerMat.col(fi) = power;
    }
 
    // Use frequency as "time" axis for the source estimate
    float fmin = (frequencies.size() > 0) ? static_cast<float>(frequencies(0)) : 0.0f;
    float fstep = (frequencies.size() > 1)
        ? static_cast<float>(frequencies(1) - frequencies(0))
        : 1.0f;
 
    InvSourceEstimate stc(powerMat, filters.vertices, fmin, fstep);
    stc.method = InvEstimateMethod::DICS;
    stc.sourceSpaceType = InvSourceSpaceType::Surface;
    stc.nVerticesLh = filters.nVerticesLh;
    stc.orientationType = filters.isFreOri ? InvOrientationType::Free : InvOrientationType::Fixed;
 
    return stc;
}
 
//=============================================================================================================
 
InvSourceEstimate InvDICS::applyDICS(const MatrixXd& data,
                                     float tmin,
                                     float tstep,
                                     const InvBeamformer& filters,
                                     int freqIdx)
{
    if (!filters.isValid() || filters.kind != "DICS") {
        qWarning("InvDICS::applyDICS - Invalid or non-DICS filters!");
        return InvSourceEstimate();
    }
    if (freqIdx < 0 || freqIdx >= filters.nFreqs()) {
        qWarning("InvDICS::applyDICS - freqIdx %d out of range (0..%d)!",
                 freqIdx, filters.nFreqs() - 1);
        return InvSourceEstimate();
    }
 
    // Apply projection + whitening + spatial filter
    MatrixXd processed = data;
    if (filters.proj.size() > 0 && filters.proj.rows() == data.rows()) {
        processed = filters.proj * processed;
    }
    if (filters.whitener.size() > 0 && filters.whitener.rows() == processed.rows()) {
        processed = filters.whitener * processed;
    }
 
    MatrixXd sol = filters.weights[freqIdx] * processed;
 
    // Combine XYZ if needed
    const int nOrient = filters.nOrient();
    if (nOrient == 3 && filters.pickOri != BeamformerPickOri::Vector) {
        const int nSources = static_cast<int>(sol.rows()) / 3;
        const int nTimes = static_cast<int>(sol.cols());
        MatrixXd combined(nSources, nTimes);
        for (int s = 0; s < nSources; ++s) {
            combined.row(s) = sol.middleRows(s * 3, 3).colwise().norm();
        }
        sol = combined;
    }
 
    InvSourceEstimate stc(sol, filters.vertices, tmin, tstep);
    stc.method = InvEstimateMethod::DICS;
    stc.sourceSpaceType = InvSourceSpaceType::Surface;
    stc.nVerticesLh = filters.nVerticesLh;
    stc.orientationType = filters.isFreOri ? InvOrientationType::Free : InvOrientationType::Fixed;
 
    return stc;
}
 
//=============================================================================================================
 
QList<InvSourceEstimate> InvDICS::applyDICSEpochs(const QList<MatrixXd>& epochs,
                                                  float tmin,
                                                  float tstep,
                                                  const InvBeamformer& filters,
                                                  int freqIdx)
{
    QList<InvSourceEstimate> results;
 
    if (epochs.isEmpty()) {
        qWarning("InvDICS::applyDICSEpochs - No epochs provided.");
        return results;
    }
    if (!filters.isValid() || filters.kind != "DICS") {
        qWarning("InvDICS::applyDICSEpochs - Invalid or non-DICS filters!");
        return results;
    }
 
    for (int i = 0; i < epochs.size(); ++i) {
        InvSourceEstimate stc = applyDICS(epochs[i], tmin, tstep, filters, freqIdx);
        if (stc.isEmpty()) {
            qWarning("InvDICS::applyDICSEpochs - Epoch %d produced empty source estimate.", i);
        }
        results.append(stc);
    }
 
    return results;
}