Skip to content

Commit

Permalink
Prop Cov , transfert matrix written
Browse files Browse the repository at this point in the history
  • Loading branch information
deseilligny committed Jul 9, 2022
1 parent 3a67db2 commit 95a640d
Show file tree
Hide file tree
Showing 5 changed files with 96 additions and 27 deletions.
2 changes: 1 addition & 1 deletion MMVII/include/MMVII_Geom2D.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ template <class Type> class cSim2D


static cSim2D FromExample(const tPt&aP0In,const tPt&aP1In,const tPt&aP0Out,const tPt&aP1Out);
static cSim2D FromExample(const std::vector<tPt>& aVIn,const std::vector<tPt>& aVOut);
static cSim2D FromExample(const std::vector<tPt>& aVIn,const std::vector<tPt>& 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);

Expand Down
22 changes: 16 additions & 6 deletions MMVII/src/Geom2D/GeomsBase2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,13 +122,13 @@ template <class TypeMap> class cLeastSquareEstimate
typedef std::vector<tPt> 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 <class TypeMap>
TypeMap cLeastSquareEstimate<TypeMap>::Estimate(const tVPts & aVIn,const tVPts & aVOut)
TypeMap cLeastSquareEstimate<TypeMap>::Estimate(const tVPts & aVIn,const tVPts & aVOut,tTypeElem * aRes2)
{
cLeasSqtAA<tTypeElem> aSys(TypeMap::NbDOF());
cDenseVect<tTypeElem> aVX(TypeMap::NbDOF());
Expand All @@ -144,18 +144,28 @@ template <class TypeMap>
aSys.AddObservation(1,aVY,aVOut[aK].y());
}
cDenseVect<tTypeElem> aSol = aSys.Solve();
TypeMap aMap = TypeMap::FromParam(aSol);

if (aRes2)
{
*aRes2 = 0.0;
for (int aK=0; aK<int(aVIn.size()) ; aK++)
{
*aRes2 += SqN2(aVOut[aK]-aMap.Value(aVIn[aK]));
}
///StdOut() << "NOOrrrmSol= " << aSomR2 << "\n";
}
/*
cDenseVect<tTypeElem> aTest = aSys.tAA() * aSol ;
StdOut() << "NOOrrrmSol= " << aTest.L2Dist(aSys.tARhs()) << "\n";
*/

return TypeMap::FromParam(aSol);
return aMap;
}

template <class Type> cSim2D<Type> cSim2D<Type>::FromExample(const std::vector<tPt>& aVIn,const std::vector<tPt>& aVOut)
template <class Type> cSim2D<Type> cSim2D<Type>::FromExample(const std::vector<tPt>& aVIn,const std::vector<tPt>& aVOut,Type * aRes2)
{
return cLeastSquareEstimate<cSim2D<Type>>::Estimate(aVIn,aVOut);
return cLeastSquareEstimate<cSim2D<Type>>::Estimate(aVIn,aVOut,aRes2);
}


Expand Down Expand Up @@ -304,7 +314,7 @@ template cSim2D<TYPE> cSim2D<TYPE>::RandomSimInv(const TYPE & AmplTr,const TYPE
template cSim2D<TYPE> cSim2D<TYPE>::FromExample(const tPt & aP0In,const tPt & aP1In,const tPt & aP0Out,const tPt & aP1Out ) ;\
template cSimilitud3D<TYPE> cSim2D<TYPE>::Ext3D() const;\
template cSim2D<TYPE> cSim2D<TYPE>::FromParam(const cDenseVect<TYPE> & aVec) ;\
template cSim2D<TYPE> cSim2D<TYPE>::FromExample(const std::vector<tPt>& aVIn,const std::vector<tPt>& aVOut);\
template cSim2D<TYPE> cSim2D<TYPE>::FromExample(const std::vector<tPt>&,const std::vector<tPt>&,TYPE*);\
template cDenseMatrix<TYPE> MatOfMul (const cPtxd<TYPE,2> & aP);


Expand Down
4 changes: 2 additions & 2 deletions MMVII/src/TutoBenchTrianguRSNL/TrianguRSNL.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
11 changes: 11 additions & 0 deletions MMVII/src/TutoBenchTrianguRSNL/cMainNetwork.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,8 @@ template <class Type> cMainNetwork <Type>::cMainNetwork
mMatrP (cMemManager::AllocMat<tPNetPtr>(mX_SzM,mY_SzM)), // alloc a grid of pointers
// Amplitude of scale muste
mSimInd2G (cSim2D<Type>::RandomSimInv(5.0,3.0,0.1))
/*
*/
{

// generate in VPix a regular grid, put them in random order for testing more config in matrix
Expand All @@ -58,6 +60,15 @@ template <class Type> cMainNetwork <Type>::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<int(mVPts.size()) ; aK++)
{
mVPts[aK].mNumPt = aK;
}
}
// Put adress of points in a grid so that they are accessible by indexes
for (auto & aPNet : mVPts)
PNetPtrOfGrid(aPNet.mInd) = & aPNet;
Expand Down
84 changes: 66 additions & 18 deletions MMVII/src/TutoBenchTrianguRSNL/cNetPropCov.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include "TrianguRSNL.h"
#include "include/MMVII_Tpl_Images.h"


namespace MMVII
Expand Down Expand Up @@ -49,15 +50,10 @@ template <class Type> cElemNetwork<Type>::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<Type>::RandomSimInv(5.0,3.0,0.3))
#endif

mSimM2This (cSim2D<Type>::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() ;
Expand All @@ -84,25 +80,77 @@ template <class Type> void cElemNetwork<Type>::CalcCov(int aNbIter)

template <class Type> void cElemNetwork<Type>::PropagCov()
{
bool Bug = (mDebugN==3);

std::vector<tPt> aVLoc;
std::vector<tPt> 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<Type> aSim = cSim2D<Type>::FromExample(aVLoc,aVMain);
Type aSqResidual;
cSim2D<Type> aSimM2L = cSim2D<Type>::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<int> aVIndTransf(this->mNum,-1);
cDenseMatrix<Type> aMatrixTranf(aNbVar,eModeInitImage::eMIA_Null); ///< Square
cDenseVect<Type> 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<Type> aVecLoc(aNbVar,eModeInitImage::eMIA_Null); ///< Square
cDenseVect<Type> 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<Type> aVLoc2 = (aMatrixTranf * aVecGlob) + aVecTranf;
cDenseVect<Type> aVDif = aVLoc2 - aVecLoc;
StdOut() << "DIF " << aVDif.L2Norm() << "\n";
*/

}

StdOut() << mDebugN << " SSS " << aSim.Tr() << " " << aSim.Sc() << "\n";
getchar();
}

/* *************************************** */
Expand Down

0 comments on commit 95a640d

Please sign in to comment.