C++ MNESourceEstimate类代码示例

bool MNESourceEstimate::read(QIODevice &p_IODevice, MNESourceEstimate& p_stc)
{
    QSharedPointer<QDataStream> t_pStream(new QDataStream(&p_IODevice));

    t_pStream->setFloatingPointPrecision(QDataStream::SinglePrecision);
    t_pStream->setByteOrder(QDataStream::BigEndian);
    t_pStream->setVersion(QDataStream::Qt_5_0);

    if(!t_pStream->device()->open(QIODevice::ReadOnly))
        return false;

    QFile* t_pFile = qobject_cast<QFile*>(&p_IODevice);
    if(t_pFile)
        printf("Reading source estimate from %s...", t_pFile->fileName().toUtf8().constData());
    else
        printf("Reading source estimate...");

    // read start time in ms
    *t_pStream >> p_stc.tmin;
    p_stc.tmin /= 1000;
    // read sampling rate in ms
    *t_pStream >> p_stc.tstep;
    p_stc.tstep /= 1000;
    // read number of vertices
    quint32 t_nVertices;
    *t_pStream >> t_nVertices;
    p_stc.vertices = VectorXi(t_nVertices);
    // read the vertex indices
    for(quint32 i = 0; i < t_nVertices; ++i)
        *t_pStream >> p_stc.vertices[i];
    // read the number of timepts
    quint32 t_nTimePts;
    *t_pStream >> t_nTimePts;
    //
    // read the data
    //
    p_stc.data = MatrixXd(t_nVertices, t_nTimePts);
    for(qint32 i = 0; i < p_stc.data.array().size(); ++i)
    {
        float value;
        *t_pStream >> value;
        p_stc.data.array()(i) = value;
    }

    //Update time vector
    p_stc.update_times();

    // close the file
    t_pStream->device()->close();

    printf("[done]\n");

    return true;
}

//.........这里部分代码省略.........
    wrtFWD_28_left.close();
    wrtFWD_29_left.close();
    wrtFWD_45_left.close();
    wrtFWD_28_right.close();
    wrtFWD_29_right.close();
    wrtFWD_45_right.close();
    wrtFWD.close();

    //
    // make an inverse operators
    //
    FiffInfo info = evoked.info;

    MNEInverseOperator inverse_operator(info, t_clusteredFwd, noise_cov, 0.2f, 0.8f);

//    QFile t_fileInv("D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_02_07_Lorenz_Esch_007/Processed/inverse operators/Lorenz-140128-Duke128-eeg-inv.fif");
//    MNEInverseOperator inverse_operator(t_fileInv);

    //
    // save clustered inverse
    //
    if(!t_sFileNameClusteredInv.isEmpty())
    {
        QFile t_fileClusteredInverse(t_sFileNameClusteredInv);
        inverse_operator.write(t_fileClusteredInverse);
    }

    //
    // Compute inverse solution
    //
    // Calculate stc for averaged data
    MinimumNorm minimumNorm(inverse_operator, lambda2, method);

    MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);

    if(sourceEstimate.isEmpty())
        return 1;

    // Calculate stc for all trials and write to file - dirty
    if(doRawTrialLocalization)
    {
        MNESourceEstimate sourceEstimate_trial;

        for(int i = 0; i<data.size(); i++)
        {
            QString file_trial_stc_orig = t_sFileNameStc;

            QString file_trial_stc_tmp = t_sFileNameStc;
            QString stc_trial_sub_folder = file_trial_stc_tmp.remove(0,t_sFileNameStc.indexOf("/So")+1);

            stc_trial_sub_folder.remove(".stc");
            stc_trial_sub_folder.append("-raw-trials/");

            QString dirCheck = data_location;
            dirCheck.append("/Processed/stc/");
            dirCheck.append(stc_trial_sub_folder);

            if(QDir(dirCheck).exists() == false)
                QDir().mkdir(dirCheck);

            stc_trial_sub_folder.prepend("stc/");

            file_trial_stc_orig.replace("stc/", stc_trial_sub_folder);

            QString temp;
            temp.append("_");

/**
* 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[])
{
    QApplication a(argc, argv);

    //########################################################################################
    //
    // Source Estimate START
    //
    //########################################################################################

    QFile t_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
    QFile t_fileCov("./MNE-sample-data/MEG/sample/sample_audvis-cov.fif");
    QFile t_fileEvoked("./MNE-sample-data/MEG/sample/sample_audvis-ave.fif");

    QString t_sFileClusteredInverse("");//QFile t_fileClusteredInverse("./clusteredInverse-inv.fif");
    AnnotationSet t_annotationSet ("sample", 2, "aparc.a2009s", "./MNE-sample-data/subjects");

    double snr = 1.0;
    double lambda2 = 1.0 / pow(snr, 2);
    QString method("dSPM"); //"MNE" | "dSPM" | "sLORETA"

    // Load data
    fiff_int_t setno = 1;
    QPair<QVariant, QVariant> baseline(QVariant(), 0);
    FiffEvoked evoked(t_fileEvoked, setno, baseline);
    if(evoked.isEmpty())
        return 1;

    std::cout << "Evoked description: " << evoked.comment.toLatin1().constData() << std::endl;

    MNEForwardSolution t_Fwd(t_fileFwd);
    if(t_Fwd.isEmpty())
        return 1;

    FiffCov noise_cov(t_fileCov);

    // regularize noise covariance
    noise_cov = noise_cov.regularize(evoked.info, 0.05, 0.05, 0.1, true);

    //
    // Cluster forward solution;
    //
    MNEForwardSolution t_clusteredFwd = t_Fwd;//.cluster_forward_solution(t_annotationSet, 40);

    //
    // make an inverse operators
    //
    FiffInfo info = evoked.info;

    MNEInverseOperator inverse_operator(info, t_clusteredFwd, noise_cov, 0.2f, 0.8f);

    if(!t_sFileClusteredInverse.isEmpty())
    {
        QFile t_fileClusteredInverse(t_sFileClusteredInverse);
        inverse_operator.write(t_fileClusteredInverse);
    }

    //
    // Compute inverse solution
    //
    MinimumNorm minimumNorm(inverse_operator, lambda2, method);
    MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);

    if(sourceEstimate.isEmpty())
        return 1;

//    // View activation time-series
//    std::cout << "\nsourceEstimate:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
//    std::cout << "time\n" << sourceEstimate.times.block(0,0,1,10) << std::endl;
//    std::cout << "timeMin\n" << sourceEstimate.times[0] << std::endl;
//    std::cout << "timeMax\n" << sourceEstimate.times[sourceEstimate.times.size()-1] << std::endl;
//    std::cout << "time step\n" << sourceEstimate.tstep << std::endl;

    //########################################################################################
    //
    //Source Estimate END
    //
    //########################################################################################

    //########################################################################################
    //
    // Create the test view START
    //
    //########################################################################################

    std::cout<<"Creating BrainView"<<std::endl;

    SurfaceSet tSurfSet ("sample", 2, "inflated", "./MNE-sample-data/subjects");
    AnnotationSet tAnnotSet ("sample", 2, "aparc.a2009s", "./MNE-sample-data/subjects");
    //Surface tSurfRight ("sample", 1, "inflated", "./MNE-sample-data/subjects");
    //Annotation tAnnotRight ("sample", 1, "aparc.a2009s", "./MNE-sample-data/subjects");
    //Surface tSurfLeft ("sample", 0, "inflated", "./MNE-sample-data/subjects");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    //
    // Cluster forward solution;
    //
    MatrixXd D;
    MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution(t_annotationSet, 20, D, noise_cov, evoked.info);

    //
    // make an inverse operators
    //
    FiffInfo info = evoked.info;

    MNEInverseOperator inverse_operator(info, t_clusteredFwd, noise_cov, 0.2f, 0.8f);

    //
    // save clustered inverse
    //
    if(!t_sFileNameClusteredInv.isEmpty())
    {
        QFile t_fileClusteredInverse(t_sFileNameClusteredInv);
        inverse_operator.write(t_fileClusteredInverse);
    }

    //
    // Compute inverse solution
    //
    MinimumNorm minimumNorm(inverse_operator, lambda2, method);

#ifdef BENCHMARK
    //
    //   Set up the inverse according to the parameters
    //
    minimumNorm.doInverseSetup(vecSel.size(),false);

    MNESourceEstimate sourceEstimate;
    QList<qint64> qVecElapsedTime;
    for(qint32 i = 0; i < 100; ++i)
    {
        //Benchmark time
        QElapsedTimer timer;
        timer.start();
        sourceEstimate = minimumNorm.calculateInverse(evoked.data, evoked.times(0), evoked.times(1)-evoked.times(0));
        qVecElapsedTime.append(timer.elapsed());
    }

    double meanTime = 0.0;
    qint32 offset = 19;
    qint32 c = 0;
    for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
    {
        meanTime += qVecElapsedTime[i];
        ++c;
    }

    meanTime /= (double)c;

    double varTime = 0;
    for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
        varTime += pow(qVecElapsedTime[i] - meanTime,2);

    varTime /= (double)c - 1.0f;
    varTime = sqrt(varTime);

    qDebug() << "MNE calculation took" << meanTime << "+-" << varTime << "ms in average";

#else
    MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);

/**
* 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[])
{
    QGuiApplication a(argc, argv);

    //########################################################################################
    // Source Estimate

    QFile t_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
    QFile t_fileEvoked("./MNE-sample-data/MEG/sample/sample_audvis-ave.fif");
    AnnotationSet t_annotationSet("sample", 2, "aparc.a2009s", "./MNE-sample-data/subjects");
    SurfaceSet t_surfSet("sample", 2, "white", "./MNE-sample-data/subjects");

    QString t_sFileNameStc("");//"RapMusic.stc");

    qint32 numDipolePairs = 7;

    bool doMovie = false;//true;

    // Parse command line parameters
    for(qint32 i = 0; i < argc; ++i)
    {
        if(strcmp(argv[i], "-stc") == 0 || strcmp(argv[i], "--stc") == 0)
        {
            if(i + 1 < argc)
                t_sFileNameStc = QString::fromUtf8(argv[i+1]);
        }else if(strcmp(argv[i], "-num") == 0 || strcmp(argv[i], "--num") == 0)
        {
            if(i + 1 < argc)
                numDipolePairs = atof(argv[i+1]);
        }
    }

    qDebug() << "Start calculation with stc:" << t_sFileNameStc;

    // Load data
    fiff_int_t setno = 0;
    QPair<QVariant, QVariant> baseline(QVariant(), 0);
    FiffEvoked evoked(t_fileEvoked, setno, baseline);
    if(evoked.isEmpty())
        return 1;

    std::cout << "evoked first " << evoked.first << "; last " << evoked.last << std::endl;

    MNEForwardSolution t_Fwd(t_fileFwd);
    if(t_Fwd.isEmpty())
        return 1;

    QStringList ch_sel_names = t_Fwd.info.ch_names;
    FiffEvoked pickedEvoked = evoked.pick_channels(ch_sel_names);

    //
    // Cluster forward solution;
    //
    MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution(t_annotationSet, 20);//40);

//    std::cout << "Size " << t_clusteredFwd.sol->data.rows() << " x " << t_clusteredFwd.sol->data.cols() << std::endl;
//    std::cout << "Clustered Fwd:\n" << t_clusteredFwd.sol->data.row(0) << std::endl;


    RapMusic t_rapMusic(t_clusteredFwd, false, numDipolePairs);

    if(doMovie)
        t_rapMusic.setStcAttr(200,0.5);


    MNESourceEstimate sourceEstimate = t_rapMusic.calculateInverse(pickedEvoked);
    if(sourceEstimate.isEmpty())
        return 1;


    QList<Label> t_qListLabels;
    QList<RowVector4i> t_qListRGBAs;

    //ToDo overload toLabels using instead of t_surfSet rr of MNESourceSpace
    t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);

    InverseView view(t_rapMusic.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24, true, false, false);//true);


    if (view.stereoType() != QGLView::RedCyanAnaglyph)
        view.camera()->setEyeSeparation(0.3f);
    QStringList args = QCoreApplication::arguments();
    int w_pos = args.indexOf("-width");
    int h_pos = args.indexOf("-height");
    if (w_pos >= 0 && h_pos >= 0)
    {
        bool ok = true;
        int w = args.at(w_pos + 1).toInt(&ok);
        if (!ok)
        {
            qWarning() << "Could not parse width argument:" << args;
            return 1;
//.........这里部分代码省略.........

/**
* 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[])
{
    QGuiApplication a(argc, argv);

    //########################################################################################
    // Source Estimate

    QFile t_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
    QFile t_fileCov("./MNE-sample-data/MEG/sample/sample_audvis-cov.fif");
    QFile t_fileEvoked("./MNE-sample-data/MEG/sample/sample_audvis-ave.fif");

    QFile t_fileClusteredInverse("./clusteredInverse-inv.fif");

    double snr = 3.0;
    double lambda2 = 1.0 / pow(snr, 2);
    QString method("sLORETA"); //"MNE" | "dSPM" | "sLORETA"

    // Load data
    fiff_int_t setno = 0;
    QPair<QVariant, QVariant> baseline(QVariant(), 0);
    FiffEvoked evoked(t_fileEvoked, setno, baseline);
    if(evoked.isEmpty())
        return 1;

    std::cout << "Evoked description: " << evoked.comment.toLatin1().constData() << std::endl;

    MNEForwardSolution t_Fwd(t_fileFwd);
    if(t_Fwd.isEmpty())
        return 1;

    AnnotationSet t_annotationSet("./MNE-sample-data/subjects/sample/label/lh.aparc.a2009s.annot", "./MNE-sample-data/subjects/sample/label/rh.aparc.a2009s.annot");

    FiffCov noise_cov(t_fileCov);

    // regularize noise covariance
    noise_cov = noise_cov.regularize(evoked.info, 0.05, 0.05, 0.1, true);

    //
    // Cluster forward solution;
    //
    MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution_ccr(t_annotationSet, 40);

    //
    // make an inverse operators
    //
    FiffInfo info = evoked.info;

    MNEInverseOperator inverse_operator(info, t_clusteredFwd, noise_cov, 0.2f, 0.8f);

    inverse_operator.write(t_fileClusteredInverse);

    //
    // Compute inverse solution
    //
    MinimumNorm minimumNorm(inverse_operator, lambda2, method);
    MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);

    if(sourceEstimate.isEmpty())
        return 1;

    // View activation time-series
    std::cout << "\nsourceEstimate:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
    std::cout << "time\n" << sourceEstimate.times.block(0,0,1,10) << std::endl;
    std::cout << "timeMin\n" << sourceEstimate.times[0] << std::endl;
    std::cout << "timeMax\n" << sourceEstimate.times[sourceEstimate.times.size()-1] << std::endl;
    std::cout << "time step\n" << sourceEstimate.tstep << std::endl;

    //Source Estimate end
    //########################################################################################

    AnnotationSet t_annotSet("./MNE-sample-data/subjects/sample/label/lh.aparc.a2009s.annot","./MNE-sample-data/subjects/sample/label/rh.aparc.a2009s.annot");
    SurfaceSet t_surfSet("./MNE-sample-data/subjects/sample/surf/lh.white", "./MNE-sample-data/subjects/sample/surf/rh.white");

    QList<Label> t_qListLabels;
    QList<RowVector4i> t_qListRGBAs;

    //ToDo overload toLabels using instead of t_surfSet rr of MNESourceSpace
    t_annotSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);

    InverseView view(minimumNorm.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24, true);

    if (view.stereoType() != QGLView::RedCyanAnaglyph)
        view.camera()->setEyeSeparation(0.3f);
    QStringList args = QCoreApplication::arguments();
    int w_pos = args.indexOf("-width");
    int h_pos = args.indexOf("-height");
    if (w_pos >= 0 && h_pos >= 0)
    {
        bool ok = true;
        int w = args.at(w_pos + 1).toInt(&ok);
        if (!ok)
        {
//.........这里部分代码省略.........

//.........这里部分代码省略.........

    FiffCov noise_cov(t_fileCov);

    // regularize noise covariance
    noise_cov = noise_cov.regularize(evoked.info, 0.05, 0.05, 0.1, true);

    //
    // Cluster forward solution;
    //
    MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution_ccr(t_annotationSet, 20);//40);

    //
    // make an inverse operators
    //
    FiffInfo info = evoked.info;

    MNEInverseOperator inverse_operator(info, t_clusteredFwd, noise_cov, 0.2f, 0.8f);

    //
    // save clustered inverse
    //
    if(!t_sFileNameClusteredInv.isEmpty())
    {
        QFile t_fileClusteredInverse(t_sFileNameClusteredInv);
        inverse_operator.write(t_fileClusteredInverse);
    }

    //
    // Compute inverse solution
    //
    MinimumNorm minimumNorm(inverse_operator, lambda2, method);

#ifdef BENCHMARK
    MNESourceEstimate sourceEstimate;
    QList<qint64> qVecElapsedTime;
    for(qint32 i = 0; i < 20; ++i)
    {
        //Benchmark time
        QElapsedTimer timer;
        timer.start();
        sourceEstimate = minimumNorm.calculateInverse(evoked);
        qVecElapsedTime.append(timer.elapsed());
    }

    double meanTime = 0.0;
    for(qint32 i = 0; i < qVecElapsedTime.size(); ++i)
        meanTime += qVecElapsedTime[i];

    meanTime /= qVecElapsedTime.size();
    qDebug() << "MNE calculation took" << meanTime << "ms in average";

#else
    MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);
#endif

    if(sourceEstimate.isEmpty())
        return 1;

    // View activation time-series
    std::cout << "\nsourceEstimate:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
    std::cout << "time\n" << sourceEstimate.times.block(0,0,1,10) << std::endl;
    std::cout << "timeMin\n" << sourceEstimate.times[0] << std::endl;
    std::cout << "timeMax\n" << sourceEstimate.times[sourceEstimate.times.size()-1] << std::endl;
    std::cout << "time step\n" << sourceEstimate.tstep << std::endl;

    //Condition Numbers

//.........这里部分代码省略.........

//    for(qint32 i = 0; i < vecSel.size(); ++i)
//    {
//        vecSel(i) = i;
//    }

//    //Option 3 - Manual selection
//    VectorXi vecSel(20);

//    vecSel << 76, 74, 13, 61, 97, 94, 75, 71, 60, 56, 26, 57, 56, 0, 52, 72, 33, 86, 96, 67;


    std::cout << "Select following epochs to average:\n" << vecSel << std::endl;

    FiffEvoked evoked = data.average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);

    QStringList ch_sel_names = t_Fwd.info.ch_names;
    FiffEvoked pickedEvoked = evoked.pick_channels(ch_sel_names);

    //########################################################################################
    // RAP MUSIC Source Estimate

    //
    // Cluster forward solution;
    //
    MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution(t_annotationSet, 20);//40);

    //
    // Compute inverse solution
    //
    PwlRapMusic t_pwlRapMusic(t_clusteredFwd, false, numDipolePairs);

#ifdef BENCHMARK
    MNESourceEstimate sourceEstimate;
    QList<qint64> qVecElapsedTime;
    for(qint32 i = 0; i < 100; ++i)
    {
        //Benchmark time
        QElapsedTimer timer;
        timer.start();
        sourceEstimate = t_pwlRapMusic.calculateInverse(pickedEvoked);
        qVecElapsedTime.append(timer.elapsed());
    }

    double meanTime = 0.0;
    qint32 offset = 19;
    qint32 c = 0;
    for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
    {
        meanTime += qVecElapsedTime[i];
        ++c;
    }

    meanTime /= (double)c;

    double varTime = 0;
    for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
        varTime += pow(qVecElapsedTime[i] - meanTime,2);

    varTime /= (double)c - 1.0f;
    varTime = sqrt(varTime);

    qDebug() << "RAP-MUSIC calculation took" << meanTime << "+-" << varTime << "ms in average";

#else
    int iWinSize = 200;

//.........这里部分代码省略.........
        }
    }

    if(data.size() > 0)
    {
        printf("Read %d epochs, %d samples each.\n",data.size(),(qint32)data[0]->epoch.cols());

//        //DEBUG
//        std::cout << data[0]->epoch.block(0,0,10,10) << std::endl;
//        qDebug() << data[0]->epoch.rows() << " x " << data[0]->epoch.cols();

//        std::cout << times.block(0,0,1,10) << std::endl;
//        qDebug() << times.rows() << " x " << times.cols();
    }




    //
    // Init RAP-MUSIC
    //
    RapMusic t_rapMusic(t_clusteredFwd, false, numDipolePairs);




    //
    // calculate the average
    //
    for(qint32 numAverages = 1; numAverages <= 20; numAverages += 1)
    {
        for(qint32 it = 0; it <= 30; ++it)
        {
            //
            // calculate the average
            //
            VectorXi vecSel(numAverages);
            srand (time(NULL)); // initialize random seed

            for(qint32 i = 0; i < vecSel.size(); ++i)
            {
                qint32 val = rand() % data.size();
                vecSel(i) = val;
            }

//            std::cout << "Select following epochs to average:\n" << vecSel << std::endl;

//            QString sSelFile = QString("aveInfo_%1_%2.txt").arg(numAverages).arg(it);
//            std::ofstream selFile(sSelFile.toLatin1().constData());
//            if (selFile.is_open())
//            {
//              selFile << vecSel << '\n';
//            }


            FiffEvoked evoked = data.average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);



            QStringList ch_sel_names = t_Fwd.info.ch_names;
            FiffEvoked pickedEvoked = evoked.pick_channels(ch_sel_names);


            //########################################################################################
            // RAP MUSIC Source Estimate

//            if(doMovie)
//                t_pwlRapMusic.setStcAttr(200,0.5);

            MNESourceEstimate sourceEstimate = t_rapMusic.calculateInverse(pickedEvoked);

//            std::cout << "Source Estimate:\n" << sourceEstimate.data << std::endl;

//            std::cout << "Source Estimate vertices:\n" << sourceEstimate.vertices << std::endl;



            if(!sourceEstimate.isEmpty())
            {
                QString t_sFileNameStc = sTargetDir+QString("%1_%2_ave_it_%3.stc").arg(sTargetPrefix).arg(numAverages).arg(it);

                qDebug() << "Write to:" << t_sFileNameStc;


                QDir dir(sTargetDir);
                if (!dir.exists()) {
                    dir.mkpath(".");
                }

                if(!t_sFileNameStc.isEmpty())
                {
                    QFile t_fileClusteredStc(t_sFileNameStc);
                    sourceEstimate.write(t_fileClusteredStc);
                }
            }
        }
    }

    return 0;//app.exec();
}

/**
* 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[])
{
    QCoreApplication app(argc, argv);

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

    QCommandLineOption evokedFileOption("ave", "Path to evoked <file>.", "file", "./MNE-sample-data/MEG/sample/sample_audvis-ave.fif");
    QCommandLineOption invFileOption("inv", "Path to inverse operator <file>.", "file", "./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif");
    QCommandLineOption snrOption("snr", "The <snr> value used for computation.", "snr", "1.0");
    QCommandLineOption methodOption("method", "Inverse estimation <method>, i.e., 'MNE', 'dSPM' or 'sLORETA'.", "method", "dSPM");
    QCommandLineOption stcFileOption("stcOut", "Path to stc <file>, which is to be written.", "file", "");

    parser.addOption(evokedFileOption);
    parser.addOption(invFileOption);
    parser.addOption(snrOption);
    parser.addOption(methodOption);
    parser.addOption(stcFileOption);
    parser.process(app);

    QFile t_fileEvoked(parser.value(evokedFileOption));
    QFile t_fileInv(parser.value(invFileOption));

    float snr = parser.value(snrOption).toFloat();
    QString method(parser.value(methodOption));
    QString t_sFileNameStc(parser.value(stcFileOption));

    double lambda2 = 1.0 / pow(snr, 2);
    qDebug() << "Start calculation with: SNR" << snr << "; Lambda" << lambda2 << "; Method" << method << "; stc:" << t_sFileNameStc;

    //
    //   Read the data first
    //
    fiff_int_t setno = 0;
    QPair<QVariant, QVariant> baseline(QVariant(), 0);
    FiffEvoked evoked(t_fileEvoked, setno, baseline);
    if(evoked.isEmpty())
        return 1;

    //
    //   Then the inverse operator
    //
    MNEInverseOperator inverse_operator(t_fileInv);

    //
    // Compute inverse solution
    //
    MinimumNorm minimumNorm(inverse_operator, lambda2, method);
    MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);

    //
    //Results
    //
    std::cout << "\npart ( block( 0, 0, 10, 10) ) of the inverse solution:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
    printf("tmin = %f s\n", sourceEstimate.tmin);
    printf("tstep = %f s\n", sourceEstimate.tstep);


    if(!t_sFileNameStc.isEmpty())
    {
        QFile t_fileStc(t_sFileNameStc);
        sourceEstimate.write(t_fileStc);

        //test if everything was written correctly
        MNESourceEstimate readSourceEstimate(t_fileStc);

        std::cout << "\npart ( block( 0, 0, 10, 10) ) of the inverse solution:\n" << readSourceEstimate.data.block(0,0,10,10) << std::endl;
        printf("tmin = %f s\n", readSourceEstimate.tmin);
        printf("tstep = %f s\n", readSourceEstimate.tstep);
    }

    return 0;//app.exec();
}

/**
* 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[])
{
    QCoreApplication a(argc, argv);

//  fname_data  - Name of the data file
//  setno       - Data set number
//  fname_inv   - Inverse operator file name
//  nave        - Number of averages (scales the noise covariance)
//             If negative, the number of averages in the data will be
//             used
//  lambda2     - The regularization factor
//  dSPM        - do dSPM?
//  sLORETA     - do sLORETA?

//    QFile t_fileEvoked("./MNE-sample-data/MEG/sample/sample_audvis-ave.fif");
//    QFile t_fileInv("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif");

    QFile t_fileEvoked("E:/Data/sl_data/MEG/mind006/mind006_051209_auditory01_raw-ave.fif");
    QFile t_fileInv("E:/Data/sl_data/MEG/mind006/mind006_051209_auditory01_raw-oct-6p-meg-inv.fif");

    float snr = 1.0f;
    QString method("dSPM"); //"MNE" | "dSPM" | "sLORETA"
    QString t_sFileNameStc("");

    // Parse command line parameters
    for(qint32 i = 0; i < argc; ++i)
    {
        if(strcmp(argv[i], "-snr") == 0 || strcmp(argv[i], "--snr") == 0)
        {
            if(i + 1 < argc)
                snr = atof(argv[i+1]);
        }
        else if(strcmp(argv[i], "-method") == 0 || strcmp(argv[i], "--method") == 0)
        {
            if(i + 1 < argc)
                method = QString::fromUtf8(argv[i+1]);
        }
        else if(strcmp(argv[i], "-stc") == 0 || strcmp(argv[i], "--stc") == 0)
        {
            if(i + 1 < argc)
                t_sFileNameStc = QString::fromUtf8(argv[i+1]);
        }
    }

    double lambda2 = 1.0 / pow(snr, 2);
    qDebug() << "Start calculation with: SNR" << snr << "; Lambda" << lambda2 << "; Method" << method << "; stc:" << t_sFileNameStc;

    //
    //   Read the data first
    //
    fiff_int_t setno = 0;
    QPair<QVariant, QVariant> baseline(QVariant(), 0);
    FiffEvoked evoked(t_fileEvoked, setno, baseline);
    if(evoked.isEmpty())
        return 1;

    //
    //   Then the inverse operator
    //
    MNEInverseOperator inverse_operator(t_fileInv);

    //
    // Compute inverse solution
    //
    MinimumNorm minimumNorm(inverse_operator, lambda2, method);
    MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);

    //
    //Results
    //
    std::cout << "\npart ( block( 0, 0, 10, 10) ) of the inverse solution:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
    printf("tmin = %f s\n", sourceEstimate.tmin);
    printf("tstep = %f s\n", sourceEstimate.tstep);


    if(!t_sFileNameStc.isEmpty())
    {
        QFile t_fileStc(t_sFileNameStc);
        sourceEstimate.write(t_fileStc);

        //test if everything was written correctly
        MNESourceEstimate readSourceEstimate(t_fileStc);

        std::cout << "\npart ( block( 0, 0, 10, 10) ) of the inverse solution:\n" << readSourceEstimate.data.block(0,0,10,10) << std::endl;
        printf("tmin = %f s\n", readSourceEstimate.tmin);
        printf("tstep = %f s\n", readSourceEstimate.tstep);
    }

    return 0;//a.exec();
}

//.........这里部分代码省略.........
    //   Read a data segment
    //   times output argument is optional
    //
    bool readSuccessful = false;
    MatrixXd data;
    MatrixXd times;
    if (in_samples)
        readSuccessful = raw.read_raw_segment(data, times, (qint32)from, (qint32)to, picks);
    else
        readSuccessful = raw.read_raw_segment_times(data, times, from, to, picks);

    if (!readSuccessful)
    {
        printf("Could not read raw segment.\n");
        return -1;
    }

    printf("Read %d samples.\n",(qint32)data.cols());

    //########################################################################################
    // RAP MUSIC Source Estimate

    //
    // Cluster forward solution;
    //
    MNEForwardSolution t_clusteredFwd = t_SelectFwd.cluster_forward_solution(t_annotationSet, 20);//t_Fwd.cluster_forward_solution_ccr(t_annotationSet, 20);//40);


    //
    // Compute inverse solution
    //
    PwlRapMusic t_pwlRapMusic(t_clusteredFwd, false, numDipolePairs);

    MNESourceEstimate sourceEstimate;

    float tstep = 1.0f/raw.info.sfreq;

    float tmin;
    if(in_samples)
        tmin = from * tstep;
    else
        tmin = from;

    //
    // Rap MUSIC Source estimate
    //
    sourceEstimate.data = MatrixXd::Zero(t_clusteredFwd.nsource, data.cols());

    //Results
    sourceEstimate.vertices = VectorXi(t_clusteredFwd.src[0].vertno.size() + t_clusteredFwd.src[1].vertno.size());
    sourceEstimate.vertices << t_clusteredFwd.src[0].vertno, t_clusteredFwd.src[1].vertno;

    sourceEstimate.times = RowVectorXf::Zero(data.cols());
    sourceEstimate.times[0] = tmin;
    for(qint32 i = 1; i < sourceEstimate.times.size(); ++i)
        sourceEstimate.times[i] = sourceEstimate.times[i-1] + tstep;
    sourceEstimate.tmin = tmin;
    sourceEstimate.tstep = tstep;


    bool first = true;
    bool last = false;

    qint32 t_iNumSensors = data.rows();
    qint32 t_iNumSteps = data.cols();