C++ MNESourceEstimate::write方法代码示例

//.........这里部分代码省略.........
        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");
//    AnnotationSet t_annotationSet("/home/chdinh/sl_data/subjects/mind006/label/lh.aparc.a2009s.annot", "/home/chdinh/sl_data/subjects/mind006/label/rh.aparc.a2009s.annot");
    AnnotationSet t_annotationSet("E:/Data/sl_data/subjects/mind006/label/lh.aparc.a2009s.annot", "E:/Data/sl_data/subjects/mind006/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(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;

    //
    // 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);
    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;

    //Condition Numbers
//    MatrixXd mags(102, t_Fwd.sol->data.cols());
//    qint32 count = 0;
//    for(qint32 i = 2; i < 306; i += 3)
//    {
//        mags.row(count) = t_Fwd.sol->data.row(i);
//        ++count;
//    }
//    MatrixXd magsClustered(102, t_clusteredFwd.sol->data.cols());
//    count = 0;
//    for(qint32 i = 2; i < 306; i += 3)
//    {
//        magsClustered.row(count) = t_clusteredFwd.sol->data.row(i);
//        ++count;

//.........这里部分代码省略.........
    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();

    qint32 t_iSamplesOverlap = (qint32)floor(((float)samplesStcWindow)*stcOverlap);
    qint32 t_iSamplesDiscard = t_iSamplesOverlap/2;

    MatrixXd measData = MatrixXd::Zero(t_iNumSensors, samplesStcWindow);

    qint32 curSample = 0;
    qint32 curResultSample = 0;
    qint32 stcWindowSize = samplesStcWindow - 2*t_iSamplesDiscard;

    while(!last)
    {
        //Data
        if(curSample + samplesStcWindow >= t_iNumSteps) //last
        {
            last = true;
            measData = data.block(0, data.cols()-samplesStcWindow, t_iNumSensors, samplesStcWindow);
        }
        else
            measData = data.block(0, curSample, t_iNumSensors, samplesStcWindow);


        curSample += (samplesStcWindow - t_iSamplesOverlap);
        if(first)
            curSample -= t_iSamplesDiscard; //shift on start t_iSamplesDiscard backwards

        //Calculate
        MNESourceEstimate stcData = t_pwlRapMusic.calculateInverse(measData, 0.0f, tstep);

        //Assign Result
        if(last)
            stcWindowSize = stcData.data.cols() - curResultSample;

        sourceEstimate.data.block(0,curResultSample,sourceEstimate.data.rows(),stcWindowSize) =
                                            stcData.data.block(0,0,stcData.data.rows(),stcWindowSize);

        curResultSample += stcWindowSize;

        if(first)
            first = false;
    }

    if(sourceEstimate.isEmpty())
        return 1;

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

    InverseView view(t_pwlRapMusic.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;
        }
        int h = args.at(h_pos + 1).toInt(&ok);
        if (!ok)
        {
            qWarning() << "Could not parse height argument:" << args;
            return 1;
        }
        view.resize(w, h);
    }
    else
    {
        view.resize(800, 600);
    }
    view.show();

    //Push Estimate
    view.pushSourceEstimate(sourceEstimate);

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

    return a.exec();
}

//.........这里部分代码省略.........
    //
    // Load data
    //
    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;
    //
    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;
    }

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

    //Only take first finger movement of each tapping session (each tapping session consists of 4 seperate finger movement)
//    VectorXi vecSel(1);
//    vecSel << 60;

    if(count<averageIterator)
        averageIterator = count;

    VectorXi vecSel(count/averageIterator);

    int i;
    for(i=0; i<count/averageIterator; i++)
        vecSel[i] = i*averageIterator;

    // Add number of averages read from the eventfile to the file name
    t_sAvrFileNameRel.append(t_sAvrFileNameRel.number(tmin));
    t_sAvrFileNameRel.append("_");
    t_sAvrFileNameRel.append(t_sAvrFileNameRel.number(tmax));
    t_sAvrFileNameRel.append("_");
    t_sAvrFileNameRel.append(t_sAvrFileNameRel.number(i));
    t_sAvrFileNameRel.append(".fif");
    QString t_sAvrFileName = t_sAvrFileNameRel.prepend(data_location);

    t_sStcFileNameRel.append(t_sStcFileNameRel.number(tmin));
    t_sStcFileNameRel.append("_");
    t_sStcFileNameRel.append(t_sStcFileNameRel.number(tmax));
    t_sStcFileNameRel.append("_");
    t_sStcFileNameRel.append(t_sStcFileNameRel.number(i));
    t_sStcFileNameRel.append("_");
    t_sStcFileNameRel.append(method.append(".stc"));
    QString t_sFileNameStc = t_sStcFileNameRel.prepend(data_location);

    std::cout << "Select following epochs to average:\n" << vecSel << std::endl;
    raw.info.filename = QString(""); //this need to be done in order to write a new file. otherwise some of the info contents are not getting copied correctly from raw.info (i.e. digitizer data)
    FiffEvoked evoked = data.average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);

    //Write averaged data to file
    QFile m_fileOut(t_sAvrFileName);
    RowVectorXd m_cals;

    FiffStream::SPtr m_pOutfid = Fiff::start_writing_raw(m_fileOut, raw.info, m_cals);
    fiff_int_t first = 0;
    m_pOutfid->write_int(FIFF_FIRST_SAMPLE, &first);
    //m_pOutfid->finish_writing_raw();

    Eigen::MatrixXd eData = evoked.data;
    Eigen::MatrixXd matValue = MatrixXd::Zero(eData.rows()+10, eData.cols());

    matValue.block(0,0,128,matValue.cols()) = evoked.data;

    m_pOutfid->write_raw_buffer(matValue, m_cals);

    m_pOutfid->finish_writing_raw();

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

    double snr = 0.1f;//1.0f;//3.0f;//0.1f;//3.0f;

    QString t_sFileNameClusteredInv("");

    // 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)

//.........这里部分代码省略.........
    // 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;
        }
        int h = args.at(h_pos + 1).toInt(&ok);
        if (!ok)
        {
            qWarning() << "Could not parse height argument:" << args;
            return 1;
        }
        view.resize(w, h);
    }
    else
    {
        view.resize(800, 600);
    }
    view.show();

    //Push Estimate
    view.pushSourceEstimate(sourceEstimate);

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

    return a.exec();
}

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

    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();
}

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

    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(doMovie) {
        t_pwlRapMusic.setStcAttr(iWinSize, 0.6f);
    }

    MNESourceEstimate sourceEstimate = t_pwlRapMusic.calculateInverse(pickedEvoked);

    if(doMovie) {
        //Select only the activations once
        MatrixXd dataPicked(sourceEstimate.data.rows(), int(std::floor(sourceEstimate.data.cols()/iWinSize)));

        for(int i = 0; i < dataPicked.cols(); ++i) {
            dataPicked.col(i) = sourceEstimate.data.col(i*iWinSize);
        }

        sourceEstimate.data = dataPicked;
    }

    if(sourceEstimate.isEmpty()) {
        return 1;
    }

#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;

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

//    //only one time point - P100
//    qint32 sample = 0;
//    for(qint32 i = 0; i < sourceEstimate.times.size(); ++i)
//    {
//        if(sourceEstimate.times(i) >= 0)
//        {
//            sample = i;
//            break;
//        }
//    }
//    sample += (qint32)ceil(0.106/sourceEstimate.tstep); //100ms
//    sourceEstimate = sourceEstimate.reduce(sample, 1);

    AbstractView::SPtr p3DAbstractView = AbstractView::SPtr(new AbstractView());
    Data3DTreeModel::SPtr p3DDataModel = p3DAbstractView->getTreeModel();

    p3DDataModel->addSurfaceSet(parser.value(subjectOption), evoked.comment, t_surfSet, t_annotationSet);

    //Add rt source loc data and init some visualization values
    if(MneEstimateTreeItem* pRTDataItem = p3DDataModel->addSourceData(parser.value(subjectOption),
                                                                      evoked.comment,
                                                                      sourceEstimate,
                                                                      t_clusteredFwd,
                                                                      t_surfSet,
                                                                      t_annotationSet)) {
        pRTDataItem->setLoopState(true);
        pRTDataItem->setTimeInterval(17);
        pRTDataItem->setNumberAverages(1);
        pRTDataItem->setStreamingState(true);
        pRTDataItem->setThresholds(QVector3D(0.01f,0.5f,1.0f));
        pRTDataItem->setVisualizationType("Annotation based");
        pRTDataItem->setColormapType("Hot");
    }

    p3DAbstractView->show();

    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);

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

    return a.exec();//1;//a.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();
}

//.........这里部分代码省略.........
//    QPair<QVariant, QVariant> baseline(QVariant(), 0);
//    FiffEvoked evoked(t_fileEvoked, setno, baseline);
//    if(evoked.isEmpty())
//        return 1;

    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_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;