From 95a640d546b813230c208f40d6fcb3d1a2812ae0 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Sat, 9 Jul 2022 17:18:26 +0200 Subject: [PATCH] Prop Cov , transfert matrix written --- MMVII/include/MMVII_Geom2D.h | 2 +- MMVII/src/Geom2D/GeomsBase2D.cpp | 22 +++-- MMVII/src/TutoBenchTrianguRSNL/TrianguRSNL.h | 4 +- .../src/TutoBenchTrianguRSNL/cMainNetwork.cpp | 11 +++ .../src/TutoBenchTrianguRSNL/cNetPropCov.cpp | 84 +++++++++++++++---- 5 files changed, 96 insertions(+), 27 deletions(-) diff --git a/MMVII/include/MMVII_Geom2D.h b/MMVII/include/MMVII_Geom2D.h index 5e9faa104e..12d7ebeaf3 100755 --- a/MMVII/include/MMVII_Geom2D.h +++ b/MMVII/include/MMVII_Geom2D.h @@ -142,7 +142,7 @@ template class cSim2D static cSim2D FromExample(const tPt&aP0In,const tPt&aP1In,const tPt&aP0Out,const tPt&aP1Out); - static cSim2D FromExample(const std::vector& aVIn,const std::vector& aVOut); + static cSim2D FromExample(const std::vector& aVIn,const std::vector& aVOut,Type * aRes2=nullptr); /// Compute a random similitude, assuring that Amplitude of scale has a minimal value static cSim2D RandomSimInv(const Type&AmplTr,const Type&AmplSc,const Type&AmplMinSc); diff --git a/MMVII/src/Geom2D/GeomsBase2D.cpp b/MMVII/src/Geom2D/GeomsBase2D.cpp index f4b384d95f..34573b0189 100755 --- a/MMVII/src/Geom2D/GeomsBase2D.cpp +++ b/MMVII/src/Geom2D/GeomsBase2D.cpp @@ -122,13 +122,13 @@ template class cLeastSquareEstimate typedef std::vector tVPts; /// Estimate the map M such that M(aVIn[aK]) = aVOut[aK] - static inline TypeMap Estimate(const tVPts& aVIn,const tVPts & aVOut); + static inline TypeMap Estimate(const tVPts& aVIn,const tVPts & aVOut,tTypeElem * aRes2); private : }; template - TypeMap cLeastSquareEstimate::Estimate(const tVPts & aVIn,const tVPts & aVOut) + TypeMap cLeastSquareEstimate::Estimate(const tVPts & aVIn,const tVPts & aVOut,tTypeElem * aRes2) { cLeasSqtAA aSys(TypeMap::NbDOF()); cDenseVect aVX(TypeMap::NbDOF()); @@ -144,18 +144,28 @@ template aSys.AddObservation(1,aVY,aVOut[aK].y()); } cDenseVect aSol = aSys.Solve(); + TypeMap aMap = TypeMap::FromParam(aSol); + if (aRes2) + { + *aRes2 = 0.0; + for (int aK=0; aK aTest = aSys.tAA() * aSol ; StdOut() << "NOOrrrmSol= " << aTest.L2Dist(aSys.tARhs()) << "\n"; */ - return TypeMap::FromParam(aSol); + return aMap; } -template cSim2D cSim2D::FromExample(const std::vector& aVIn,const std::vector& aVOut) +template cSim2D cSim2D::FromExample(const std::vector& aVIn,const std::vector& aVOut,Type * aRes2) { - return cLeastSquareEstimate>::Estimate(aVIn,aVOut); + return cLeastSquareEstimate>::Estimate(aVIn,aVOut,aRes2); } @@ -304,7 +314,7 @@ template cSim2D cSim2D::RandomSimInv(const TYPE & AmplTr,const TYPE template cSim2D cSim2D::FromExample(const tPt & aP0In,const tPt & aP1In,const tPt & aP0Out,const tPt & aP1Out ) ;\ template cSimilitud3D cSim2D::Ext3D() const;\ template cSim2D cSim2D::FromParam(const cDenseVect & aVec) ;\ -template cSim2D cSim2D::FromExample(const std::vector& aVIn,const std::vector& aVOut);\ +template cSim2D cSim2D::FromExample(const std::vector&,const std::vector&,TYPE*);\ template cDenseMatrix MatOfMul (const cPtxd & aP); diff --git a/MMVII/src/TutoBenchTrianguRSNL/TrianguRSNL.h b/MMVII/src/TutoBenchTrianguRSNL/TrianguRSNL.h index 750a370287..6a57a81a06 100755 --- a/MMVII/src/TutoBenchTrianguRSNL/TrianguRSNL.h +++ b/MMVII/src/TutoBenchTrianguRSNL/TrianguRSNL.h @@ -35,8 +35,8 @@ namespace NS_Bench_RSNL { #define DEBUG_RSNL true #if (DEBUG_RSNL) -constexpr double AMPL_Grid2Real= 0.0; // Amplitude of random differerence between real position and regular grid -constexpr double AMPL_Real2Init = 0.0; // Amplitude of random and syst differerence betwen real an init position +constexpr double AMPL_Grid2Real= 0.00; // Amplitude of random differerence between real position and regular grid +constexpr double AMPL_Real2Init = 0.00; // Amplitude of random and syst differerence betwen real an init position #else constexpr double AMPL_Grid2Real= 0.1; // Amplitude of random differerence between real position and regular grid constexpr double AMPL_Real2Init = 0.1; // Amplitude of random and syst differerence betwen real an init position diff --git a/MMVII/src/TutoBenchTrianguRSNL/cMainNetwork.cpp b/MMVII/src/TutoBenchTrianguRSNL/cMainNetwork.cpp index 940b95658d..3dd26780cd 100755 --- a/MMVII/src/TutoBenchTrianguRSNL/cMainNetwork.cpp +++ b/MMVII/src/TutoBenchTrianguRSNL/cMainNetwork.cpp @@ -37,6 +37,8 @@ template cMainNetwork ::cMainNetwork mMatrP (cMemManager::AllocMat(mX_SzM,mY_SzM)), // alloc a grid of pointers // Amplitude of scale muste mSimInd2G (cSim2D::RandomSimInv(5.0,3.0,0.1)) +/* +*/ { // generate in VPix a regular grid, put them in random order for testing more config in matrix @@ -58,6 +60,15 @@ template cMainNetwork ::cMainNetwork aVCoord0.push_back(aPNet.mPosInit.y()); } } + if (1) // Eventually put in random order to check implicit order of NumX,NumY is not used + { + mVPts = RandomOrder(mVPts); + // But need to reset the num of points which is used in link construction + for (int aK=0 ; aK cElemNetwork::cElemNetwork(tMainNW & aMainNW,const c tMainNW (eModeSSR::eSSR_LsqDense,cRect2(cPt2di(0,0),aBoxM.Sz()),false), mMainNW (&aMainNW), mBoxM (aBoxM), -#if (DEBUG_RSNL) - mSimM2This (tPt(0,0),tPt(1,0)) -#else - (cSim2D::RandomSimInv(5.0,3.0,0.3)) -#endif - + mSimM2This (cSim2D::RandomSimInv(5.0,3.0,0.3)) { -static int TheNumDebug=0; -mDebugN = ++TheNumDebug; + static int TheNumDebug=0; + mDebugN = ++TheNumDebug; // To have the right scale compute mSimInd2G from mSimM2This this->mSimInd2G = mSimM2This * mMainNW->SimInd2G() ; @@ -84,25 +80,77 @@ template void cElemNetwork::CalcCov(int aNbIter) template void cElemNetwork::PropagCov() { -bool Bug = (mDebugN==3); - std::vector aVLoc; std::vector aVMain; + for (const auto & aPNet : this->mVPts) { const tPNet & aHomMain = MainHom(aPNet); aVLoc.push_back(aPNet.PCur()); aVMain.push_back(aHomMain.PCur()); - if (1|| Bug) - { - StdOut() << "PTS " << aPNet.PCur() << " " << aHomMain.PCur() << "\n"; - StdOut() << " I " << aPNet.mInd << aPNet.mInd - aHomMain.mInd << "\n"; - } } - cSim2D aSim = cSim2D::FromExample(aVLoc,aVMain); + Type aSqResidual; + cSim2D aSimM2L = cSim2D::FromExample(aVMain,aVLoc,&aSqResidual); + tPt aTr = aSimM2L.Tr(); + tPt aSc = aSimM2L.Sc(); +/* + Loc = aSimM2L * Main + + X_loc (Trx) (Sx -Sy) (X_Main) + Y_loc = (Try) + (Sx -Sy) * (Y_Main) + +*/ + + int aNbVar = this->mNum; + std::vector aVIndTransf(this->mNum,-1); + cDenseMatrix aMatrixTranf(aNbVar,eModeInitImage::eMIA_Null); ///< Square + cDenseVect aVecTranf(aNbVar,eModeInitImage::eMIA_Null); ///< Square + + for (const auto & aPNet : this->mVPts) + { + const tPNet & aHomMain = MainHom(aPNet); + int aKx = aPNet.mNumX; + int aKy = aPNet.mNumY; + aVIndTransf.at(aKx) = aHomMain.mNumX; + aVIndTransf.at(aKy) = aHomMain.mNumY; + + aVecTranf(aKx) = aTr.x(); + aVecTranf(aKy) = aTr.y(); + + aMatrixTranf.SetElem(aKx,aKx,aSc.x()); + aMatrixTranf.SetElem(aKy,aKx,-aSc.y()); + aMatrixTranf.SetElem(aKx,aKy,aSc.y()); + aMatrixTranf.SetElem(aKy,aKy,aSc.x()); + } + + if (DEBUG_RSNL) + { + cDenseVect aVecLoc(aNbVar,eModeInitImage::eMIA_Null); ///< Square + cDenseVect aVecGlob(aNbVar,eModeInitImage::eMIA_Null); ///< Square + for (const auto & aPNet : this->mVPts) + { + const tPNet & aHomMain = MainHom(aPNet); + int aKx = aPNet.mNumX; + int aKy = aPNet.mNumY; + tPt aPLoc = aPNet.PCur(); + tPt aPGlob = aHomMain.PCur(); + + aVecLoc(aKx) = aPLoc.x(); + aVecLoc(aKy) = aPLoc.y(); + aVecGlob(aKx) = aPGlob.x(); + aVecGlob(aKy) = aPGlob.y(); + } + // aVecGlob + aVecTranf; +/* + + cDenseVect aVLoc2 = (aMatrixTranf * aVecGlob) + aVecTranf; + cDenseVect aVDif = aVLoc2 - aVecLoc; + + StdOut() << "DIF " << aVDif.L2Norm() << "\n"; +*/ + + } - StdOut() << mDebugN << " SSS " << aSim.Tr() << " " << aSim.Sc() << "\n"; - getchar(); } /* *************************************** */