C++ Im1D_REAL8::data方法代码示例

void Craig_etal_L1
     (
             Im2D_REAL8  A,
             Im1D_REAL8  B,
             REAL        TOLER,
             Im1D_REAL8  SOL,
             Im1D_REAL8  RESIDU
     )
{
    INT n = SOL.tx();
    INT m = B.tx();

    BENCH_ASSERT
    (
           (A.tx() == n+2)
        && (A.ty() == m+2)
        && (B.tx() == m)
        && (SOL.tx() == n)
        && (RESIDU.tx() == m)
    );
    Craig_Barrodale_Roberts_l1
    (
        m,n,
        A.data_lin(),
        B.data(),
        TOLER,
        SOL.data(),
        RESIDU.data()
    );
}

Im1D_REAL8 cManipOrdInc::ReordonneSol(Im1D_REAL8 aIm)
{
   Im1D_REAL8 aRes (aIm.tx());
   for (int aK=0 ; aK<aIm.tx() ; aK++)
      aRes.data()[aK] = aIm.data()[mAlloc2Solve[aK]];

   return aRes;
}

Pt3dr cGenSysSurResol::Pt3dSolInter(bool * aOk)
{
     Im1D_REAL8  aSol = GSSR_Solve(aOk);
     if (aOk && (! *aOk))
        return Pt3dr(1e33,-1e44,0);
     ELISE_ASSERT(aSol.tx()==3,"cGenSysSurResol::SolInter");

     return Pt3dr(aSol.data()[0],aSol.data()[1],aSol.data()[2]);
}

void cGenSysSurResol::GSSR_SolveEqFitDroite(REAL & aAx,REAL &aB,bool * aOk)
{
     Im1D_REAL8  aSol = GSSR_Solve(aOk);
     if (aOk && (! *aOk))
        return;
     ELISE_ASSERT(aSol.tx()==2,"cGenSysSurResol::GSSR_SolveEqFitDroite");
     aAx = aSol.data()[0];
     aB = aSol.data()[1];
}

void cGenSysSurResol::GSSR_SolveEqFitPlan(REAL & aAx,REAL &aBy,REAL & aC,bool * aOk)
{
     Im1D_REAL8  aSol = GSSR_Solve(aOk);
     if (aOk && (! *aOk))
        return;
     ELISE_ASSERT(aSol.tx()==3,"cGenSysSurResol::GSSR_SolveEqFiPlan");
     aAx = aSol.data()[0];
     aBy = aSol.data()[1];
     aC  = aSol.data()[2];
}

REAL SystLinSurResolu::Residu(Im1D_REAL8 anIm,INT iEq) const
{
   AssertIndexEqValide(iEq);
   AssertIndexGoodNbVar(anIm.tx());

   return Residu(anIm.data(),iEq);
}

void cIncEnsembleCamera::OptLineaireOnDirL2
     (
           std::list<cIncSetLiaison *> *        aListSL,
           const std::vector<cFonctrPond> &     aFoncAux
           // bool                                 CalcMatr
     )
{
    mListSl = aListSL;
    mFoncsAux = aFoncAux;

    SetOpt();
/*
    if (CalcMatr)
    {
        ResetEquation();
        ScoreCurGen(true,true);
    }
*/

    Im1D_REAL8  anIm =  mSys->GSSR_Solve((bool *)0);


    SetImState0();
    SetCurDir(anIm.data());

    cIEC_OptimCurDir anOpti(*this);
    REAL aLambda = anOpti.Optim(0,1);

    SetLambda(aLambda);
}

void cIncEnsembleCamera::OptimPowel
     (
          std::list<cIncSetLiaison *> *  aListSL,
          const std::vector<cFonctrPond> &     aFoncAux,
          REAL tol,
          INT ItMax
     )
{
    mListSl = aListSL;
    mFoncsAux = aFoncAux;
    SetOpt();
    Im1D_REAL8 I = CurVals();

    cIEC_OptimPowel aOpt(*this);
    aOpt.powel(I.data(),tol,ItMax);
    SetPtCur(I.data());
}

void SystLinSurResolu::PushEquation 
     ( 
            Im1D_REAL8    aFormLin,
            REAL          aValue,
            REAL          aPds
     )
{
   AssertIndexGoodNbVar(aFormLin.tx());
   PushEquation (aFormLin.data(),aValue,aPds);
}

void MakeFoncRepart(Im1D_REAL8 aH,int * aVMax=0)
{
    double aNbP;
    ELISE_COPY(aH.all_pts(),aH.in(),sigma(aNbP));
    REAL8 * aDH = aH.data();
    for (int aK=1 ; aK<aNbH ; aK++)
    {
        if (aDH[aK] && aVMax) *aVMax = aK;
        aDH[aK] += aDH[aK-1];
    }
    ELISE_COPY(aH.all_pts(),aH.in() * (255.0/aNbP),aH.out());
}

void cMEPCoCentrik::OneItereRotPur(ElMatrix<REAL>  & aMat,double & anErrStd)
{
    L2SysSurResol aSysLin3(3);
    aSysLin3.GSSR_Reset(false);

    std::vector<double> aVRes;
    double aSomP=0;
    double aSomErr=0;
    for (ElPackHomologue::const_iterator itP=mPack.begin() ; itP!=mPack.end() ; itP++)
    {
         Pt3dr aQ1 = vunit(PZ1(itP->P1()));
         Pt3dr aQ2 =  aMat * vunit(PZ1(itP->P2()));
         double aVQ2[3],aVQ1[3];
         aQ2.to_tab(aVQ2);
         aQ1.to_tab(aVQ1);

         double anEcart = euclid(aQ1-aQ2);
         aVRes.push_back(anEcart);
         double aPds =  itP->Pds() / (1 + ElSquare(anEcart / (2*anErrStd)));

         aSomP += aPds;
         aSomErr += aPds * square_euclid(aQ1-aQ2);;

         ElMatrix<REAL>  aMQ2 =  MatProVect(aQ2);
         for (int aY=0 ; aY< 3 ; aY++)
         {
             double aCoeff[3];
             for (int aX=0 ; aX< 3 ; aX++)
                 aCoeff[aX] = aMQ2(aX,aY);

             aSysLin3.GSSR_AddNewEquation(aPds,aCoeff,aVQ2[aY]-aVQ1[aY],0);
         }
    }
    double anErrRobut = MoyKPPVal(aVRes,NbForEcart(aVRes.size()));
    double anErrQuad = sqrt(aSomErr/aSomP);
    if (0)
    {
       std::cout << "ERR QUAD " << anErrQuad * mFoc << " Robust " << anErrRobut * mFoc << "\n";
    }
    anErrStd =  anErrQuad;
    Im1D_REAL8   aSol = aSysLin3.GSSR_Solve (0);
    double * aData = aSol.data();

    ElMatrix<double> aMPV =  MatProVect(Pt3dr(aData[0],aData[1],aData[2]));

     // std::cout << " aData " << aData[0] <<  " " << aData[1] << " " << aData[2] << "\n";
   

    // aMat  = NearestRotation(aMat * (ElMatrix<double>(3,true) +aMPV));
    aMat  = NearestRotation((ElMatrix<double>(3,true) +aMPV) * aMat);
}

void bench_triviale_opt_sous_contrainte()
{
   // Miminise x2+y2, sous la contrainte x+y=2
     L2SysSurResol aSys(2);
     double C[2] = {1,1};
     aSys.GSSR_AddContrainte(C,3);

     double Fx[2] = {1,0};
     aSys.GSSR_AddNewEquation(1.0,Fx,0);
     double Fy[2] = {0,1};
     aSys.GSSR_AddNewEquation(1.0,Fy,0);

     Im1D_REAL8 aSol = aSys.GSSR_Solve(0);
     BENCH_ASSERT(std::abs(aSol.data()[0] -1.5)<epsilon);
     BENCH_ASSERT(std::abs(aSol.data()[1] -1.5)<epsilon);

}

void cASAMG::ComputeIncidGradProf()
{
   Im2DGen * aImProf = mStdN->ImProf();
   L2SysSurResol  aSys (3);
   Pt2di aP0;

   Im2D_Bits<1> aMasqTmp = ImMarqueurCC(mSz);
   TIm2DBits<1> aTMasqTmp(aMasqTmp);


   for (aP0.x=0 ; aP0.x<mSz.x ; aP0.x++)
   {
       for (aP0.y=0 ; aP0.y<mSz.y ; aP0.y++)
       {
            double Angle = 1.5;
            if (mTMasqN.get(aP0))
            {
                cCalcPlanImage  aCalcPl(CCDist(),aSys,aImProf);
                int aNb = OneZC(aP0,CCV4(),aTMasqTmp,1,0,mTMasqN,1,aCalcPl);
                ResetMarqueur(aTMasqTmp,aCalcPl.mVPts);
                if (aNb >= SeuimNbPtsCCDist())
                {

                    Im1D_REAL8 aSol = aSys.Solve((bool *)0);
                    double * aDS = aSol.data();
                    Pt2dr aGrad (aDS[0],aDS[1]);

                    Angle = euclid(aGrad);
                }
/*
*/
            }
            mTIncid.oset(aP0,ElMin(255,ElMax(0,round_ni(Angle*DynAng()))));
       }
   }
}

void cIncEnsembleCamera::OptimJacobi
     (
          std::list<cIncSetLiaison *> *  aListSL,
          const std::vector<cFonctrPond> &     aFoncAux
     )
{
    mListSl = aListSL;
    mFoncsAux = aFoncAux;
    SetOpt();
    ELISE_ASSERT(mL2Opt,"Need L2 Sys for OptimJacobi");

    mSysL2->GetMatr(mMatrL2,mMatrtB);
    jacobi_diag(mMatrL2,mMatrValP,mMatrVecP);

    mtBVecP = mMatrtB * mMatrVecP;

cElRanGen aR;
Im1D_REAL8 P0 = CurVals();
Im1D_REAL8 mImBest  = CurVals();


REAL aScInit = ScoreCur(false);
REAL aScMin = aScInit;
REAL aScAmMin = aScInit;
for (INT aTest =0 ; aTest < 100 ; aTest ++)
{
    SetPtCur(P0.data());
    Im1D_REAL8 aDir( NbVal(),0.0);
    

    for (INT aK = 0 ; aK < NbVal() ; aK++)
    {
        ELISE_ASSERT(mMatrValP(aK,aK) != 0,"Jcobi ");
        REAL Val = (mtBVecP(aK,0) /  ElAbs(mMatrValP(aK,aK)))  ;

	REAL aRan = aR.cNRrandom3()  ;

        if (aRan < 0.25)
        {
             Val = 0;
        }
        else if (aRan < 0.5)
	;
	/*
        {
             Val = Val; 
        }
	*/
        else
            Val *=  3 * aR.cNRrandom3() -1;
     

        for (INT aY =0 ; aY< NbVal() ; aY++)
            aDir.data()[aY] += mMatrVecP(aK,aY) * Val;

    }

    ELISE_COPY(aDir.all_pts(),P0.in() + aDir.in(),aDir.out());
    SetPtCur(aDir.data());
    REAL aSc = ScoreCur(false);

if (aSc < aScMin)
{
    ElSetMin(aScMin,aSc);

    for (INT aK= 0 ; aK< 10 ; aK++)
        OneItereDevL1( aK < 3);

    REAL aSc = ScoreCur(false);
    if (aSc < aScAmMin)
    {
         aScAmMin = aSc;
         ELISE_COPY(mImBest.all_pts(),CurVals().in(),mImBest.out());
    }
   
}
    cout  <<  aScMin << " " << aScInit  << " " << aSc <<  " " << aScAmMin << "\n";
}
    SetPtCur(mImBest.data());
}

int Luc_main_XAlign(int argc,char ** argv)
{
    //MMD_InitArgcArgv(argc,argv,3);

    std::string aFilePtsIn;
    //Reading the arguments
    ElInitArgMain
    (
        argc,argv,
        LArgMain()  << EAMC(aFilePtsIn,"Input file"),
        LArgMain()
    );

    std::string aFilePtsOut="GCP_xAligned.xml";

    std::ifstream file(aFilePtsIn.c_str(), ios::in);
    int nbIm;
    file >> nbIm;
    std::vector<Pt3dr> aVPts(nbIm);
    std::vector<Pt3dr> aVInc(nbIm);
    std::vector<std::string> aVName(nbIm,"");
    for(int i=0 ; i<nbIm ; i++)
    {
        string name;
        file >> aVName[i] >> aVPts[i].x >> aVPts[i].y >> aVPts[i].z >> aVInc[i].x >> aVInc[i].y >> aVInc[i].z;
    }

    file.close();
    //Least Square

    // Create L2SysSurResol to solve least square equation with 3 unknown
    L2SysSurResol aSys(2);

    //For Each SIFT point
    double sumX=0, sumY=0;
    for(int i=0;i<int(aVPts.size());i++){
        double aPds[2]={aVPts[i].x,1};
        double poids=1;
        aSys.AddEquation(poids,aPds,aVPts[i].y);
        sumX=sumX+aVPts[i].x;
        sumY=sumY+aVPts[i].y;
    }

    Pt2dr aRotCenter; aRotCenter.x=sumX/aVPts.size();aRotCenter.y=sumY/aVPts.size();

    bool Ok;
    Im1D_REAL8 aSol = aSys.GSSR_Solve(&Ok);

    double aAngle;
    if (Ok)
    {
        double* aData = aSol.data();
        aAngle=atan(aData[0]);
        cout<<"Angle = "<<aAngle<<endl<<"Rot Center = "<<aRotCenter<<endl;

    for(int i=0;i<int(aVPts.size());i++){
        Pt2dr aPt; aPt.x=aVPts[i].x; aPt.y=aVPts[i].y;
        aPt=Rot2D(aAngle, aPt, aRotCenter);aVPts[i].x=aPt.x;aVPts[i].y=aPt.y;
    }
    }

//End Least Square

    cDicoAppuisFlottant  aDico;
    for (int aKP=0 ; aKP<int(aVPts.size()) ; aKP++)
    {
        cOneAppuisDAF aOAD;
        aOAD.Pt() = aVPts[aKP];
        aOAD.NamePt() = aVName[aKP];
        aOAD.Incertitude() = aVInc[aKP];

        aDico.OneAppuisDAF().push_back(aOAD);
    }


    MakeFileXML(aDico,aFilePtsOut);

    return 0;
}