From f3ece253935b7b3035dac965a57ff43e8f7a8a97 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Thu, 30 May 2024 19:03:37 +0200 Subject: [PATCH 1/5] In topol-sadle detect 4 checkboard target --- MMVII/include/MMVII_Ptxd.h | 5 + .../CodedTarget/cCheckBoardTargetExtract.cpp | 207 ++++++++++++++++-- MMVII/src/Geom2D/GeomsBase2D.cpp | 3 + 3 files changed, 195 insertions(+), 20 deletions(-) diff --git a/MMVII/include/MMVII_Ptxd.h b/MMVII/include/MMVII_Ptxd.h index e12b32eefc..20dc4bbce6 100755 --- a/MMVII/include/MMVII_Ptxd.h +++ b/MMVII/include/MMVII_Ptxd.h @@ -204,6 +204,11 @@ typedef cPtxd cPt3df ; // If Dim=3 1-> 6 2-> 3->26 ( 3^3 -1) template const std::vector> & AllocNeighbourhood(int aNbVois); +// classical 8-Neighbourhood +const std::vector & Alloc8Neighbourhood(); +// classical 4-Neighbourhood +const std::vector & Alloc4Neighbourhood(); + // Create a tab where K entrie represent vectors having NormInf equal to K // !! => Entry go from 0 to aDistMax included // !! => the size can be larger (but obviously not smaller) than dist required, as function remumber previous calls .... diff --git a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp index bdd954feb4..92b2737912 100755 --- a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp @@ -3,6 +3,7 @@ #include "MMVII_Tpl_Images.h" #include "MMVII_Interpolators.h" #include "MMVII_Linear2DFiltering.h" +#include // Test git branch @@ -36,8 +37,7 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli void DoOneImage() ; int IsTopoSaddlePoint(const cPt2di & aPt) const; - - + void MakeImageSaddlePoints(const tDIm &) const; // =========== Mandatory args ============ std::string mNameIm; ///< Name of background image @@ -51,6 +51,7 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli // =========== Internal param ============ tIm mImIn; ///< Input global image + cPt2di mSzIm; tDIm * mDImIn; }; @@ -86,24 +87,149 @@ cCollecSpecArg2007 & cAppliCheckBoardTargetExtract::ArgOpt(cCollecSpecArg2007 & << AOpt2007(mShowQuickSaddleF,"ShowQSF","Vector to show quick saddle filters") ; } -#if (0) -#endif -// cPt2di FreemanV8[8] = {{ + + +/** Compute, as a flagt of bit, the set of Fremaan-8 neighboor that are > over a pixel, + * for value==, the comparison is done on Y then X + */ + + +template tU_INT1 FlagSup8Neigh(const cDataIm2D & aDIm,const cPt2di & aPt) +{ + Type aV0 = aDIm.GetV(aPt); + tU_INT1 aResult = 0; + // for freeman 1,4 Y is positive, for freeman 0, X is positive + for (int aK=0 ; aK<4 ; aK++) + { + if (aDIm.GetV(aPt+FreemanV8[aK]) >= aV0) + aResult |= 1<< aK; + } + for (int aK=4 ; aK<8 ; aK++) + { + if (aDIm.GetV(aPt+FreemanV8[aK]) > aV0) + aResult |= 1<< aK; + } + return aResult; +} + +class cCCOfOrdNeigh +{ + public : + bool mIsSup; // is it a component with value > pixel + int mBit0; // First neighboorr + int mBit1; // last neighboor +}; + +/* +void ComputeCCOfeNeigh(std::vector> & aVCC) +{ + for (size_t aFlag=0 ; aFlag<256 ; aFlag++) + { + for (size_t aBit=0 ; aBit<8 ; aBit++) + { + int aBitPrec = (7+aBit) % 8; + bool ThisIs0 = (aFlag& (1< & aVBool,size_t aFlag) +{ +} +*/ + +void ComputeNbbCCOfFlag(tU_INT1 * aTabFlag) +{ + for (size_t aFlag=0 ; aFlag<256 ; aFlag++) + { + tU_INT1 aNbCC = 0; + for (size_t aBit=0 ; aBit<8 ; aBit++) + { + int aNextBit = (1+aBit) % 8; + bool ThisIs0 = (aFlag& (1< tREAL8 SadleDecision(const cDataIm2D & aDIm,const cPt2di & aPt) +{ + +} +*/ + + +void ConnectedComponent + ( + std::vector & aVPts, + cDataIm2D & aDIm, + const std::vector & aNeighbourhood, + const cPt2di& aSeed, + int aMarqInit=1, + int aNewMarq=0 + ) +{ + aVPts.clear(); + if (aDIm.GetV(aSeed) != aMarqInit) + return; + + aDIm.SetV(aSeed,aNewMarq); + aVPts.push_back(aSeed); + size_t aIndBottom = 0; + + while (aIndBottom!=aVPts.size()) + { + cPt2di aP0 = aVPts[aIndBottom]; + for (const auto & aDelta : aNeighbourhood) + { + cPt2di aNeigh = aP0 + aDelta; + if (aDIm.GetV(aNeigh)==aMarqInit) + { + aDIm.SetV(aNeigh,aNewMarq); + aVPts.push_back(aNeigh); + } + } + aIndBottom++; + } +} + + int cAppliCheckBoardTargetExtract::IsTopoSaddlePoint(const cPt2di & aPt) const { tREAL4 aV0 = mDImIn->GetV(aPt); - bool Sup0[8]; + std::vector Sup0(8); for (int aK=0 ; aK<4 ; aK++) { - Sup0[aK] = mDImIn->GetV(aPt+FreemanV8[aK]) > aV0; + Sup0[aK] = mDImIn->GetV(aPt+FreemanV8[aK]) >= aV0; } for (int aK=4 ; aK<8 ; aK++) { - Sup0[aK] = mDImIn->GetV(aPt+FreemanV8[aK]) >= aV0; + Sup0[aK] = mDImIn->GetV(aPt+FreemanV8[aK]) > aV0; } int aNbDist = 0; @@ -117,35 +243,76 @@ int cAppliCheckBoardTargetExtract::IsTopoSaddlePoint(const cPt2di & aPt) const return aNbDist; } +void cAppliCheckBoardTargetExtract::MakeImageSaddlePoints(const tDIm & aDIm) const +{ + cRGBImage aRGB = RGBImFromGray(aDIm); + // cRect2 aRectInt = aDIm->Dilate(-1); + + for (const auto & aPix : cRect2(aDIm.Dilate(-1))) + { + if (NbbCCOfFlag(FlagSup8Neigh(aDIm,aPix))>=4) + { + aRGB.SetRGBPix(aPix,cRGBImage::Red); + } + } + aRGB.ToFile("Saddles.tif"); +} void cAppliCheckBoardTargetExtract::DoOneImage() { mImIn = tIm::FromFile(mNameIm); mDImIn = &mImIn.DIm() ; + mSzIm = mDImIn->Sz(); + StdOut() << "BITSET , 4 " << sizeof(std::bitset<4>) << " 8:" << sizeof(std::bitset<8>) << " 9:" << sizeof(std::bitset<8>) << "\n"; + StdOut() << "END READ IMAGE \n"; - new cIm2D(cPt2di(2,2)); - new cIm2D(cPt2di(2,2)); + SquareAvgFilter(*mDImIn,4,1,1); - // ExpFilterOfStdDev(*mDImIn,4,2.0); + StdOut() << "END FILTER \n"; - SquareAvgFilter(*mDImIn,4,1,1); - mDImIn->ToFile("tooooottttoto.tif"); + // mDImIn->ToFile("tooooottttoto.tif"); cRect2 aRectInt = mDImIn->Dilate(-1); - int aNbExtre=0; int aNbSaddle=0; - int aNbStd=0; int aNbTot=0; + + // MakeImageSaddlePoints(*mDImIn); + + cIm2D aImMasq(mSzIm,nullptr,eModeInitImage::eMIA_Null); + cDataIm2D & aDMasq = aImMasq.DIm(); for (const auto & aPix : aRectInt) { + if (NbbCCOfFlag(FlagSup8Neigh(*mDImIn,aPix)) >=4) + { + aImMasq.DIm().SetV(aPix,1); + aNbSaddle++; + } aNbTot++; - int aNb = IsTopoSaddlePoint(aPix); - if (aNb==0) aNbExtre++; - if (aNb==2) aNbStd++; - if (aNb==4) aNbSaddle++; } - StdOut() << "Saddle = " << aNbSaddle << " %=" << (100.0*aNbSaddle)/double(aNbTot) << "\n"; + int aNbCCSad=0; + std::vector aVCC; + const std::vector & aV8 = Alloc8Neighbourhood(); + for (const auto& aPix : aDMasq) + { + if (aDMasq.GetV(aPix)==1) + { + aNbCCSad++; + ConnectedComponent(aVCC,aDMasq,aV8,aPix,1,0); + for (const auto & aPixCC : aVCC) + { + FakeUseIt(aPixCC); + } + } + } + + + StdOut() << "NBS " << (100.0*aNbSaddle)/aNbTot << " " << (100.0*aNbCCSad)/aNbTot << "\n"; + /* + for (const auto & aPix : aRectInt) + { + } + */ } diff --git a/MMVII/src/Geom2D/GeomsBase2D.cpp b/MMVII/src/Geom2D/GeomsBase2D.cpp index 7713258c47..dc88b927b7 100755 --- a/MMVII/src/Geom2D/GeomsBase2D.cpp +++ b/MMVII/src/Geom2D/GeomsBase2D.cpp @@ -8,6 +8,9 @@ namespace MMVII const cPt2di FreemanV8[8] {{1,0},{1,1},{0,1},{-1,1},{-1,0},{-1,-1},{0,-1},{1,-1}}; const cPt2di FreemanV10[10] {{1,0},{1,1},{0,1},{-1,1},{-1,0},{-1,-1},{0,-1},{1,-1},{1,0},{1,1}}; +const std::vector & Alloc4Neighbourhood() { return AllocNeighbourhood<2>(1); } +const std::vector & Alloc8Neighbourhood() { return AllocNeighbourhood<2>(2); } + /* ========================== */ /* cBox2di */ /* ========================== */ From 2ba15153fa80f5e3563401b1508dc3d7bea9aaaf Mon Sep 17 00:00:00 2001 From: deseilligny Date: Tue, 4 Jun 2024 19:54:38 +0200 Subject: [PATCH 2/5] Checkboard : center +- OK (saddle/Sym) --- MMVII/include/MMVII_2Include_Tiling.h | 78 +++- MMVII/include/MMVII_HeuristikOpt.h | 55 +++ MMVII/include/MMVII_ImageMorphoMath.h | 109 +++++ MMVII/include/MMVII_Mappings.h | 2 + MMVII/include/MMVII_Matrix.h | 1 + MMVII/include/MMVII_SysSurR.h | 6 +- MMVII/include/MMVII_nums.h | 7 +- MMVII/src/Appli/cMMVII_Ap_CPU.cpp | 2 +- MMVII/src/CodedTarget/FilterCodedTarget.h | 21 +- .../CodedTarget/cCheckBoardTargetExtract.cpp | 417 ++++++++++-------- MMVII/src/CodedTarget/cFilterTarget.cpp | 14 + MMVII/src/CodedTarget/cSaddleFiter.cpp | 121 ++++- MMVII/src/Geoms/TilesIndexGeom.cpp | 35 +- MMVII/src/Mappings/MappingGen.cpp | 1 + MMVII/src/Matrix/BaseMatrixes.cpp | 4 + MMVII/src/MorphoMat/ConnectedComponents.cpp | 64 +++ MMVII/src/TieP/cToTiePMultiple.cpp | 2 +- MMVII/src/UtiMaths/cOptimByStep.cpp | 135 ++++++ 18 files changed, 860 insertions(+), 214 deletions(-) create mode 100755 MMVII/include/MMVII_HeuristikOpt.h create mode 100755 MMVII/include/MMVII_ImageMorphoMath.h create mode 100644 MMVII/src/MorphoMat/ConnectedComponents.cpp create mode 100644 MMVII/src/UtiMaths/cOptimByStep.cpp diff --git a/MMVII/include/MMVII_2Include_Tiling.h b/MMVII/include/MMVII_2Include_Tiling.h index cc9c537af6..4e4d6a3a29 100755 --- a/MMVII/include/MMVII_2Include_Tiling.h +++ b/MMVII/include/MMVII_2Include_Tiling.h @@ -116,17 +116,17 @@ template bool EqualPt(const cPtxd & aP1,const cPtxd */ -template class cTiling : public cTilingIndex +template class cTiling : public cTilingIndex { public : - typedef cTilingIndex tTI; + typedef cTilingIndex tTI; typedef typename tTI::tRBox tRBox; typedef typename tTI::tRPt tRPt; typedef typename tTI::tIPt tIPt; typedef typename Type::tPrimGeom tPrimGeom; typedef typename Type::tArgPG tArgPG; - static constexpr int Dim = Type::Dim; + static constexpr int Dim = Type::TheDim; typedef std::list tCont1Tile; typedef std::vector tVectTiles; @@ -248,14 +248,21 @@ FakeUseIt(aDistWMargin); tArgPG mArgPG; }; + +/** We consider 2d point that are valuated with z() + * Extract point that do not contain any point > in a neighbourhood of size aDist, + * Non Const vect because they are sorted at the begining of process + */ +std::vector FilterMaxLoc(std::vector & aVpt,tREAL8 aDist); + /** Sometime we only need to put points in a spatial index, this * class make the interface allowing to use simple point in spatial indexing */ -template class cPointSpInd +template class cPointSpInd { public : - static constexpr int Dim = TheDim; + static constexpr int TheDim = Dim; typedef cPtxd tPrimGeom; typedef int tArgPG; /// unused here @@ -293,6 +300,67 @@ template class cGeneratePointDiff }; +template + void IndexeseFilterMaxLoc + ( + TypePrim *, // Type specifier + std::vector & aVInd, + const std::vector & aVPt, + const TypeCalcP2 & aCalc, + tREAL8 aDist + ) +{ + // compute box of points + cTplBoxOfPts aBoxOfPt; + for (const auto & aPt : aVPt) + aBoxOfPt.Add(aCalc(aPt)); + + cTplBox aBox = aBoxOfPt.CurBox(); + + // estimate number of case + int aNbCase = std::min + ( + round_up(aVPt.size()/30.0) , + round_up( 0.1*(aBox.NbElem()/pow(aDist,TypePrim::TheDim)) ) + ); + + + // creat tiling + clear indexe + cTiling> aTiling(aBox,true,aNbCase,-1); + aVInd.clear(); + + // parse and add + for (size_t aKPt=0 ; aKPt + std::vector + FilterMaxLoc + ( + TypePrim *, // Type specifier + const std::vector & aVPt, + const TypeCalcP2 & aCalc, + tREAL8 aDist + ) +{ + std::vector aVInd; + IndexeseFilterMaxLoc((TypePrim*) nullptr,aVInd,aVPt,aCalc,aDist); + + std::vector aRes; + for (const auto & aInd : aVInd) + aRes.push_back(aVPt.at(aInd)); + return aRes; +} + + }; #endif // _MMVII_Tiling_H_ diff --git a/MMVII/include/MMVII_HeuristikOpt.h b/MMVII/include/MMVII_HeuristikOpt.h new file mode 100755 index 0000000000..4c6bb76dd8 --- /dev/null +++ b/MMVII/include/MMVII_HeuristikOpt.h @@ -0,0 +1,55 @@ +#ifndef _MMVII_Heuristik_Opt_H_ +#define _MMVII_Heuristik_Opt_H_ + +#include "MMVII_Ptxd.h" + + +namespace MMVII +{ + +/** \file "MMVII_HeuristikOpt.h" + \brief Declaration functionnalities for optimization that are not based on local linearization, + typicall used on image matching +*/ + +/** Class for optimizing by local descent function that are not differenciable. Method : + * + * - multi-step approach + * - iteratively research the optimum in a neigourhood (for example can be V4 or V8 if dim=2) + * - take cautious not to visit twice the same point (this is the only "non trivial part") + */ + +typedef std::vector> tVLInt; + +template class cOptimByStep +{ + public : + typedef cPtxd tPtR; + typedef cPtxd tPtI; + typedef cDataMapping tMap; + + cOptimByStep(const tMap & aMap,bool IsMin,tREAL8 aMaxDInfInit,int aDist1Max=Dim); + std::pair Optim(const tPtR & ,tREAL8 aStepInit,tREAL8 aStepLim,tREAL8 aMul=0.5); + + private : + std::pair CurValue() const; + bool DoOneStep(tREAL8 aStep); ///< return false if point go too far from init + // void OneStep(const tPtR + tREAL8 Value(const tPtR& aPt) const; + + const tMap & mFunc; ///< contain the func to optimize + tREAL8 mSign; ///< +1 for min, -1 for max + cWhichMin mWMin; ///< store best result and arg-best + tREAL8 mMaxDInfInit; ///< maximal dist to initial value we accept + int mDist1Max; ///< max dist-1 of neighboorhooud + std::vector> mINeigh; ///< list of integer neighboorhoud + const tVLInt & mIndNextN; ///< given last neigh, will give fast access to next neigh + tPtR mPt0; ///< memorize init point to avoid we go too far + +}; + + + +}; // namespace MMVII + +#endif // _MMVII_Heuristik_Opt_H_ diff --git a/MMVII/include/MMVII_ImageMorphoMath.h b/MMVII/include/MMVII_ImageMorphoMath.h new file mode 100755 index 0000000000..8c110d6cea --- /dev/null +++ b/MMVII/include/MMVII_ImageMorphoMath.h @@ -0,0 +1,109 @@ +#ifndef _MMVII_MorphoMat_H_ +#define _MMVII_MorphoMat_H_ + +#include "MMVII_Image2D.h" + + +namespace MMVII +{ + +/** \file "MMVII_ImageMorphoMath.h" + \brief Declaration functionnalities for morpho-mat & graph like processing +*/ + +/* *********************************************** */ +/* */ +/* Graph & Flag of bits */ +/* */ +/* *********************************************** */ + +/** Class for representing the connected component of 8-neigboorhoud, using Freeman code */ + +class cCCOfOrdNeigh +{ + public : + bool mIs1; // is it a component with 1 + int mFirstBit; // First neighboor (include) + int mLastBit; // last neighboor (exclude) +}; + + +/// class for storing a subset of neighboors as flag of bit on Freman neighbourhood + + +class cFlag8Neigh +{ + public : + + cFlag8Neigh() : mFlag(0) {} + cFlag8Neigh(tU_INT1 aFlag) : mFlag(aFlag) {} + + /// If bit to 1 , safe because work with bit<0 or bit>=8 (a bit slower ...) + inline bool SafeIsIn(size_t aBit) const { return (mFlag& (1< & ConComp() const + { + static std::vector > aVCC; + static bool First=true; + if (First) + { + First = false; + ComputeCCOfeNeigh(aVCC); + } + return aVCC.at(mFlag); + } + + /// compute the number of connected component of a given neighbourhood flag + inline int NbConComp() const { return ConComp().size(); } + + private : + tU_INT1 mFlag; + + static void ComputeCCOfeNeigh(std::vector> & aVCC); +}; + +/** Compute, as a flag of bit, the set of Fremaan-8 neighboor that are > over a pixel, + * to have a strict order, for value==, the comparison is done on Y then X */ +template cFlag8Neigh FlagSup8Neigh(const cDataIm2D & aDIm,const cPt2di & aPt); + + +/** Given a point, that is a "topologicall saddle" (i.e multiple connected component of points >) , + * carecterize the "importance" of saddle */ +template tREAL8 CriterionTopoSadle(const cDataIm2D & aDIm,const cPt2di & aPixC); + + +/** Put in vector "VPts", the connected component of a seed , of point having the same colour in Image, + * update with "NewMarq" */ + +void ConnectedComponent + ( + std::vector & aVPts, // Vector of results + cDataIm2D & aDIm, // Image marquing point of CC to compute + const std::vector & aNeighbourhood, // generally 4 or 8 neighbourhood + const std::vector& aSeed, // seeds of the connected component + int aMarqInit=1, // will extract CC of point having this colour + int aNewMarq=0, // will marq points reached with this colour + bool Restore=false // if true restore initial marq at end + ); + +/** Idem, current and simplified case where the seed is a single point */ +void ConnectedComponent + ( + std::vector & , cDataIm2D & , const std::vector & aNeighbourhood, + const cPt2di& aSeed, // seed of the connected component + int aMarqInit=1, int aNewMarq=0, bool Restore=false + ); + +}; + + + + + +#endif // _MMVII_MorphoMat_H_ diff --git a/MMVII/include/MMVII_Mappings.h b/MMVII/include/MMVII_Mappings.h index c51c011eff..afdc45116e 100755 --- a/MMVII/include/MMVII_Mappings.h +++ b/MMVII/include/MMVII_Mappings.h @@ -373,6 +373,8 @@ template class cDataMapping : publ #endif // MAP_STATIC_BUF }; +typedef cDataMapping tFunc2DReal; + /** Specialization for DimIn=DimOut , introduce because we want to force invertible mapping to have DimIn==DimOut */ diff --git a/MMVII/include/MMVII_Matrix.h b/MMVII/include/MMVII_Matrix.h index 9638d311e4..5f01b285a4 100755 --- a/MMVII/include/MMVII_Matrix.h +++ b/MMVII/include/MMVII_Matrix.h @@ -154,6 +154,7 @@ template class cDenseVect MMVII_INTERNAL_ASSERT_medium(aV.IsInside(Sz()) ,"Sparse Vector out dense vect"); } void WeightedAddIn(Type aWeight,const tSpV & aColLine); + void WeightedAddIn(Type aWeight,const tDV & aColLine); /* ========= Othognalization & projection stuff =========== */ diff --git a/MMVII/include/MMVII_SysSurR.h b/MMVII/include/MMVII_SysSurR.h index 31e33a5c69..aa5deefe25 100755 --- a/MMVII/include/MMVII_SysSurR.h +++ b/MMVII/include/MMVII_SysSurR.h @@ -584,12 +584,12 @@ template class cLeasSqtAA : public cLeasSq // ================ Accessor used in Schur elim ======== : - const cDenseMatrix & tAA () const; ///< Accessor + const cDenseMatrix & tAA () const; ///< Accessor , warn not symetrized const cDenseVect & tARhs () const; ///< Accessor - cDenseMatrix & tAA () ; ///< Accessor + cDenseMatrix & tAA () ; ///< Accessor , warn not symetrized cDenseVect & tARhs () ; ///< Accessor - /// access to tAA via virtual interface + /// access to tAA via virtual interface, duplicate then do the symetrization cDenseMatrix V_tAA() const override; /// access to tARhs via virtual interface cDenseVect V_tARhs() const override; diff --git a/MMVII/include/MMVII_nums.h b/MMVII/include/MMVII_nums.h index 9da123fca8..e43a81bff6 100755 --- a/MMVII/include/MMVII_nums.h +++ b/MMVII/include/MMVII_nums.h @@ -657,14 +657,17 @@ template class cWhichExtrem } bool IsInit() const {return mIsInit;} - void Add(const TypeIndex & anIndex,const TypeVal & aNewVal) + // return value indicate if modif was done + bool Add(const TypeIndex & anIndex,const TypeVal & aNewVal) { if ( (IsMin?(aNewVal=mValExtre)) || (!mIsInit)) { mValExtre = aNewVal; mIndexExtre = anIndex; + mIsInit = true; + return true; } - mIsInit = true; + return false; } const TypeIndex & IndexExtre() const {AssertIsInit();return mIndexExtre;} const TypeVal & ValExtre () const {AssertIsInit();return mValExtre;} diff --git a/MMVII/src/Appli/cMMVII_Ap_CPU.cpp b/MMVII/src/Appli/cMMVII_Ap_CPU.cpp index 57a036fedc..3f7e4dcd08 100755 --- a/MMVII/src/Appli/cMMVII_Ap_CPU.cpp +++ b/MMVII/src/Appli/cMMVII_Ap_CPU.cpp @@ -127,7 +127,7 @@ cAutoTimerSegm::~cAutoTimerSegm() /* */ /***********************************/ -static const std::string DefTime(" OTHERS"); +static const std::string DefTime("OTHERS"); cTimerSegm::cTimerSegm(cMMVII_Ap_CPU * anAppli) : mLastIndex (DefTime), diff --git a/MMVII/src/CodedTarget/FilterCodedTarget.h b/MMVII/src/CodedTarget/FilterCodedTarget.h index ac125b3e21..52ac329a8a 100755 --- a/MMVII/src/CodedTarget/FilterCodedTarget.h +++ b/MMVII/src/CodedTarget/FilterCodedTarget.h @@ -4,6 +4,7 @@ #include "MMVII_Linear2DFiltering.h" #include "MMVII_Geom2D.h" #include "MMVII_SysSurR.h" +#include "MMVII_Mappings.h" /** \file FilterCodedTarget.h @@ -119,7 +120,7 @@ class cParamAllFilterDCT */ -template class cFilterDCT : public cMemCheck +template class cFilterDCT : public tFunc2DReal { public : typedef cIm2D tIm; @@ -129,6 +130,8 @@ template class cFilterDCT : public cMemCheck // =============== Allocator ============= static cFilterDCT * AllocSym(tIm anIm,const cParamAllFilterDCT &); + static cFilterDCT * AllocSym(tIm anIm,double aR0,double aR1,double aEpsilon); + // static cFilterDCT * AllocBin(tIm anIm,double aR0,double aR1); static cFilterDCT * AllocRad(const tImGr & aImGr,const cParamAllFilterDCT &); static cFilterDCT * AllocBin(tIm anIm,const cParamAllFilterDCT &); @@ -144,12 +147,14 @@ template class cFilterDCT : public cMemCheck virtual void UpdateSelected(cDCT & aDC) const ; double ComputeVal(const cPt2dr & aP); + cPt1dr Value(const cPt2dr&) const; tIm ComputeIm(); double ComputeValMaxCrown(const cPt2dr & aP,const double& aThreshold); tIm ComputeImMaxCrown(const double& aThreshold); eDCTFilters ModeF() const; + protected : cFilterDCT(bool IsCumul,eDCTFilters aMode,tIm anIm,bool IsSym,double aR0,double aR1,double aThickN=1.5); cFilterDCT(tIm anIm,const cParamAllFilterDCT & aGlob,const cParam1FilterDCT *); @@ -263,6 +268,8 @@ class cCompiledNeighBasisFunc * Typically a faster implemantation will store for once the var/covar matrix and its inverse * (depend only of the neighbourhood) */ cDenseVect SlowCalc(const std::vector & ); + + const cDenseVect & FastCalc(const std::vector & ); private : inline void AssertValsIsOk(const std::vector & aVV) { @@ -275,6 +282,10 @@ class cCompiledNeighBasisFunc std::vector> mVFuncs; ///< for each neighboor contain value on the basis size_t mNbFunc; ///< number of functions, commodity cLeasSqtAA mSys; ///< least square system + bool mFisrtFast; ///< First time we do Fast => compute mMatInvCov + cDenseMatrix mMatInvCov; ///< Invers of cov as it does not vary + cDenseVect mRHS; + cDenseVect mSolFast; }; /** Class specific to saddle point computation using a least square fitting of neighbourood by quadratic @@ -287,16 +298,20 @@ class cCalcSaddle cCalcSaddle(double aRay,double aStep); /// not use 4 now, compute criteria/eigen values - tREAL8 CalcSaddleCrit(const std::vector &,bool Show); + tREAL8 CalcSaddleCrit(const std::vector &,bool Show=false); + + /// Dont use interpol, require step = 1 + tREAL8 CalcSaddleCrit(cDataIm2D& aIm,cPt2di); /// compute refined displacement a saddle point using values of neihboor cPt2dr RefineSadlPtFromVals(const std::vector & aVVals,bool Show); /// optimize position by iteration on RefineSadlPtFromVals + bool RefineSadlePointFromIm(cIm2D aIm,cPt2dr & aPt,bool ResetIfDivg=false); void RefineSadlePointFromIm(cIm2D aIm,cDCT & aDCT); double Ray() const { return mRay;} //CM: avoid mRay unused double Step() const { return mStep;} //CM: avoid mStep unused - + const std::vector & VINeigh() const ; ///< Accessor private : double mRay; ///< Radius of neighbourhood double mStep; ///< Step for discretization diff --git a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp index 92b2737912..19b30f622f 100755 --- a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp @@ -3,7 +3,10 @@ #include "MMVII_Tpl_Images.h" #include "MMVII_Interpolators.h" #include "MMVII_Linear2DFiltering.h" -#include +#include "MMVII_ImageMorphoMath.h" +#include "MMVII_2Include_Tiling.h" +#include "MMVII_Sensor.h" +#include "MMVII_HeuristikOpt.h" // Test git branch @@ -11,6 +14,14 @@ namespace MMVII { + +static constexpr tU_INT1 eNone = 0 ; +static constexpr tU_INT1 eTopo0 = 1 ; +static constexpr tU_INT1 eTopoTmpCC = 2 ; +static constexpr tU_INT1 eTopoMaxOfCC = 3 ; +static constexpr tU_INT1 eTopoMaxLoc = 4 ; +static constexpr tU_INT1 eFilterSym = 5 ; + /* *********************************************************** */ /* */ /* cAppliCheckBoardTargetExtract */ @@ -35,9 +46,17 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli cCollecSpecArg2007 & ArgObl(cCollecSpecArg2007 & anArgObl) override ; cCollecSpecArg2007 & ArgOpt(cCollecSpecArg2007 & anArgOpt) override ; + bool IsPtTest(const cPt2dr & aPt) const; ///< Is it a point marqed a test + + void DoOneImage() ; - int IsTopoSaddlePoint(const cPt2di & aPt) const; - void MakeImageSaddlePoints(const tDIm &) const; + void MakeImageSaddlePoints(const tDIm &,const cDataIm2D & aDMasq) const; + + cWhichMin OptimFilter(cFilterDCT *,const cPt2dr &) const; + + cPhotogrammetricProject mPhProj; + cTimerSegm mTimeSegm; + // =========== Mandatory args ============ std::string mNameIm; ///< Name of background image @@ -50,9 +69,13 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli std::vector mShowQuickSaddleF; // =========== Internal param ============ - tIm mImIn; ///< Input global image - cPt2di mSzIm; - tDIm * mDImIn; + tIm mImIn; ///< Input global image + cPt2di mSzIm; ///< Size of image + tDIm * mDImIn; ///< Data input image + bool mHasMasqTest; ///< Do we have a test image 4 debuf (with masq) + cIm2D mMasqTest; ///< Possible image of mas 4 debug, print info ... + cIm2D mImLabel; ///< Image storing labels of centers + cDataIm2D * mDImLabel; ///< Data Image of label }; @@ -64,8 +87,14 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector & aVArgs,const cSpecMMVII_Appli & aSpec) : cMMVII_Appli (aVArgs,aSpec), + mPhProj (*this), + mTimeSegm (this), mImIn (cPt2di(1,1)), - mDImIn (nullptr) + mDImIn (nullptr), + mHasMasqTest (false), + mMasqTest (cPt2di(1,1)), + mImLabel (cPt2di(1,1)), + mDImLabel (nullptr) { } @@ -84,241 +113,279 @@ cCollecSpecArg2007 & cAppliCheckBoardTargetExtract::ArgOpt(cCollecSpecArg2007 & { return anArgOpt - << AOpt2007(mShowQuickSaddleF,"ShowQSF","Vector to show quick saddle filters") + << mPhProj.DPMask().ArgDirInOpt("TestMask","Mask for selecting point used in detailed mesg/output") + << AOpt2007(mShowQuickSaddleF,"ShowQSF","Vector to show quick saddle filters") ; } - -/** Compute, as a flagt of bit, the set of Fremaan-8 neighboor that are > over a pixel, - * for value==, the comparison is done on Y then X - */ - - -template tU_INT1 FlagSup8Neigh(const cDataIm2D & aDIm,const cPt2di & aPt) +void cAppliCheckBoardTargetExtract::MakeImageSaddlePoints(const tDIm & aDIm,const cDataIm2D & aDMasq) const { - Type aV0 = aDIm.GetV(aPt); - tU_INT1 aResult = 0; - // for freeman 1,4 Y is positive, for freeman 0, X is positive - for (int aK=0 ; aK<4 ; aK++) - { - if (aDIm.GetV(aPt+FreemanV8[aK]) >= aV0) - aResult |= 1<< aK; - } - for (int aK=4 ; aK<8 ; aK++) - { - if (aDIm.GetV(aPt+FreemanV8[aK]) > aV0) - aResult |= 1<< aK; - } - return aResult; -} - -class cCCOfOrdNeigh -{ - public : - bool mIsSup; // is it a component with value > pixel - int mBit0; // First neighboorr - int mBit1; // last neighboor -}; + cRGBImage aRGB = RGBImFromGray(aDIm); -/* -void ComputeCCOfeNeigh(std::vector> & aVCC) -{ - for (size_t aFlag=0 ; aFlag<256 ; aFlag++) - { - for (size_t aBit=0 ; aBit<8 ; aBit++) - { - int aBitPrec = (7+aBit) % 8; - bool ThisIs0 = (aFlag& (1<= (int) eTopoMaxLoc) + { + aRGB.SetRGBPix(aPix,(aDMasq.GetV(aPix)==eFilterSym) ? cRGBImage::Red : cRGBImage::Green ); + } + } + aRGB.ToFile("Saddles.tif"); } -void FlagToVect(std::vector & aVBool,size_t aFlag) +bool cAppliCheckBoardTargetExtract::IsPtTest(const cPt2dr & aPt) const { + return mHasMasqTest && (mMasqTest.DIm().GetV(ToI(aPt)) != 0); } -*/ -void ComputeNbbCCOfFlag(tU_INT1 * aTabFlag) +cWhichMin cAppliCheckBoardTargetExtract::OptimFilter(cFilterDCT * aFilter,const cPt2dr & aP0) const { - for (size_t aFlag=0 ; aFlag<256 ; aFlag++) + cWhichMin aRes; + for (tREAL8 aDx=-1.0 ; aDx<1.0 ; aDx += 0.1) { - tU_INT1 aNbCC = 0; - for (size_t aBit=0 ; aBit<8 ; aBit++) - { - int aNextBit = (1+aBit) % 8; - bool ThisIs0 = (aFlag& (1<ComputeValMaxCrown(aPt,1e10); + tREAL8 aVal = aFilter->ComputeVal(aPt); + aRes.Add(aPt,aVal); + // UpdateMin(aMinSym,aFSym->ComputeValMaxCrown(aP2+cPt2dr(aDx,aDy),10)); } - aTabFlag[aFlag] = aNbCC; } + return aRes; } -int NbbCCOfFlag(tU_INT1 aFlag) -{ - static tU_INT1 TabFlag[256]; - static bool First=true; - if (First) - { - First = false; - ComputeNbbCCOfFlag(TabFlag); - } - return TabFlag[aFlag]; -} -/* -template tREAL8 SadleDecision(const cDataIm2D & aDIm,const cPt2di & aPt) -{ +void cAppliCheckBoardTargetExtract::DoOneImage() +{ + int mNbBlur1 = 4; // Number of iteration of initial blurring + tREAL8 mDistMaxLocSad = 10.0; // for supressing sadle-point, not max loc in a neighboorhoud + int mMaxNbMLS = 2000; // Max number of point in best saddle points + tREAL8 aRayCalcSadle = sqrt(4+1); // limit point 2,1 -} -*/ - - -void ConnectedComponent - ( - std::vector & aVPts, - cDataIm2D & aDIm, - const std::vector & aNeighbourhood, - const cPt2di& aSeed, - int aMarqInit=1, - int aNewMarq=0 - ) -{ - aVPts.clear(); - if (aDIm.GetV(aSeed) != aMarqInit) - return; + tREAL8 mThresholdSym = 0.50; + tREAL8 mDistCalcSym0 = 8.0; + tREAL8 mDistDivSym = 2.0; + - aDIm.SetV(aSeed,aNewMarq); - aVPts.push_back(aSeed); - size_t aIndBottom = 0; + // computed threshold + tINT8 mDistRectInt = 20; // to see later how we compute it - while (aIndBottom!=aVPts.size()) - { - cPt2di aP0 = aVPts[aIndBottom]; - for (const auto & aDelta : aNeighbourhood) - { - cPt2di aNeigh = aP0 + aDelta; - if (aDIm.GetV(aNeigh)==aMarqInit) - { - aDIm.SetV(aNeigh,aNewMarq); - aVPts.push_back(aNeigh); - } - } - aIndBottom++; - } -} + /* [0] Initialise : read image and mask */ + cAutoTimerSegm aTSInit(mTimeSegm,"Init"); + // [0.0] read image + mImIn = tIm::FromFile(mNameIm); + mDImIn = &mImIn.DIm() ; + mSzIm = mDImIn->Sz(); + cRect2 aRectInt = mDImIn->Dilate(-mDistRectInt); -int cAppliCheckBoardTargetExtract::IsTopoSaddlePoint(const cPt2di & aPt) const -{ - tREAL4 aV0 = mDImIn->GetV(aPt); - std::vector Sup0(8); + // [0.1] initialize labeling image + //mImLabel(mSzIm,nullptr,eModeInitImage::eMIA_Null); + mDImLabel = &(mImLabel.DIm()); + mDImLabel->Resize(mSzIm); + mDImLabel->InitCste(eNone); - for (int aK=0 ; aK<4 ; aK++) - { - Sup0[aK] = mDImIn->GetV(aPt+FreemanV8[aK]) >= aV0; - } - for (int aK=4 ; aK<8 ; aK++) - { - Sup0[aK] = mDImIn->GetV(aPt+FreemanV8[aK]) > aV0; - } - int aNbDist = 0; - for (int aK=0 ; aK<8 ; aK++) - { - if (Sup0[aK]!=Sup0[(aK+1)%8]) - aNbDist++; - } + mHasMasqTest = mPhProj.ImageHasMask(mNameIm); + if (mHasMasqTest) + mMasqTest = mPhProj.MaskOfImage(mNameIm,*mDImIn); - return aNbDist; -} + /* [1] Compute a blurred image => less noise, less low level saddle */ -void cAppliCheckBoardTargetExtract::MakeImageSaddlePoints(const tDIm & aDIm) const -{ - cRGBImage aRGB = RGBImFromGray(aDIm); - // cRect2 aRectInt = aDIm->Dilate(-1); + cAutoTimerSegm aTSBlur(mTimeSegm,"Blurr"); + + tIm aImBlur = mImIn.Dup(); // create image blurred with less noise + tDIm& aDImBlur = aImBlur.DIm(); + + SquareAvgFilter(aDImBlur,mNbBlur1,1,1); - for (const auto & aPix : cRect2(aDIm.Dilate(-1))) - { - if (NbbCCOfFlag(FlagSup8Neigh(aDIm,aPix))>=4) - { - aRGB.SetRGBPix(aPix,cRGBImage::Red); - } - } - aRGB.ToFile("Saddles.tif"); -} -void cAppliCheckBoardTargetExtract::DoOneImage() -{ - mImIn = tIm::FromFile(mNameIm); - mDImIn = &mImIn.DIm() ; - mSzIm = mDImIn->Sz(); - StdOut() << "BITSET , 4 " << sizeof(std::bitset<4>) << " 8:" << sizeof(std::bitset<8>) << " 9:" << sizeof(std::bitset<8>) << "\n"; - StdOut() << "END READ IMAGE \n"; - SquareAvgFilter(*mDImIn,4,1,1); + /* [2] Compute "topological" saddle point */ - StdOut() << "END FILTER \n"; + cAutoTimerSegm aTSTopoSad(mTimeSegm,"TopoSad"); - // mDImIn->ToFile("tooooottttoto.tif"); + // 2.1 point with criteria on conexity of point > in neighoor - cRect2 aRectInt = mDImIn->Dilate(-1); int aNbSaddle=0; int aNbTot=0; - // MakeImageSaddlePoints(*mDImIn); - - cIm2D aImMasq(mSzIm,nullptr,eModeInitImage::eMIA_Null); - cDataIm2D & aDMasq = aImMasq.DIm(); for (const auto & aPix : aRectInt) { - if (NbbCCOfFlag(FlagSup8Neigh(*mDImIn,aPix)) >=4) + if (FlagSup8Neigh(aDImBlur,aPix).NbConComp() >=4) { - aImMasq.DIm().SetV(aPix,1); + mDImLabel->SetV(aPix,eTopo0); aNbSaddle++; } aNbTot++; } + + // 2.2 as often there 2 "touching" point with this criteria + // select 1 point in conected component + + cAutoTimerSegm aTSMaxCC(mTimeSegm,"MaxCCSad"); int aNbCCSad=0; std::vector aVCC; const std::vector & aV8 = Alloc8Neighbourhood(); - for (const auto& aPix : aDMasq) + + for (const auto& aPix : *mDImLabel) { - if (aDMasq.GetV(aPix)==1) + if (mDImLabel->GetV(aPix)==eTopo0) { aNbCCSad++; - ConnectedComponent(aVCC,aDMasq,aV8,aPix,1,0); + ConnectedComponent(aVCC,*mDImLabel,aV8,aPix,eTopo0,eTopoTmpCC); + cWhichMax aBestPInCC; for (const auto & aPixCC : aVCC) { - FakeUseIt(aPixCC); + aBestPInCC.Add(aPixCC,CriterionTopoSadle(aDImBlur,aPixCC)); } + + cPt2di aPCC = aBestPInCC.IndexExtre(); + mDImLabel->SetV(aPCC,eTopoMaxOfCC); } } + /* [3] Compute point that are max local */ - StdOut() << "NBS " << (100.0*aNbSaddle)/aNbTot << " " << (100.0*aNbCCSad)/aNbTot << "\n"; - /* - for (const auto & aPix : aRectInt) + cAutoTimerSegm aTSCritSad(mTimeSegm,"CritSad"); + + std::vector aVSad0; + cCalcSaddle aCalcSBlur(aRayCalcSadle+0.001,1.0); + + for (const auto& aPix : *mDImLabel) { + if (mDImLabel->GetV(aPix)==eTopoMaxOfCC) + { + tREAL8 aCritS = aCalcSBlur.CalcSaddleCrit(aDImBlur,aPix); + aVSad0.push_back(cPt3dr(aPix.x(),aPix.y(),aCritS)); + } } - */ + + // [3.2] select KBest + MaxLocal + cAutoTimerSegm aTSMaxLoc(mTimeSegm,"MaxLoc"); + + SortOnCriteria(aVSad0,[](const auto & aPt){return - aPt.z();} ); + std::vector aVMaxLoc = FilterMaxLoc((cPt2dr*)nullptr,aVSad0,[](const auto & aP) {return Proj(aP);}, mDistMaxLocSad); + + // limit the number of point , a bit rough but first experiment show that sadle criterion is almost perfect on good images + aVMaxLoc.resize(std::min(aVMaxLoc.size(),size_t(mMaxNbMLS))); + + for (const auto & aP3 : aVMaxLoc) + mDImLabel->SetV(ToI(Proj(aP3)),eTopoMaxLoc); + + StdOut() << "END MAXLOC \n"; + + + + /* [4] Calc Symetry criterion */ + + cAutoTimerSegm aTSSym(mTimeSegm,"SYM"); + if (1) + { + cCalcSaddle aCalcSInit(1.5,1.0); + std::vector aNewP3; + cFilterDCT * aFSym = cFilterDCT::AllocSym(mImIn,0.0,mDistCalcSym0,1.0); + cOptimByStep aOptimSym(*aFSym,true,mDistDivSym); + for (const auto & aP3 : aVMaxLoc) + { + cPt2dr aP0 = Proj(aP3); + + auto [aValSym,aNewP] = aOptimSym.Optim(aP0,1.0,0.01); + + if (aValSym< mThresholdSym) + { + aNewP3.push_back(cPt3dr(aNewP.x(),aNewP.y(),aValSym)); + mDImLabel->SetV(ToI(aNewP),eFilterSym); + } + } + + delete aFSym; + aVMaxLoc = aNewP3; + } + + if (0) + { + cFilterDCT * aFSym = cFilterDCT::AllocSym(mImIn,0.0,mDistCalcSym0,1.0); + cOptimByStep aOptimSym(*aFSym,true,mDistDivSym); +// tPtR Optim(const tPtR & ,tREAL8 aStepInit,tREAL8 aStepLim,tREAL8 aMul=0.5); + + + cStdStatRes aSSad1; + cStdStatRes aSSad0; + cStdStatRes aSSymInt_1; + cStdStatRes aSSymInt_0; + + cCalcSaddle aCalcSInit(1.5,1.0); + for (const auto & aP3 : aVMaxLoc) + { + cPt2dr aP0 = Proj(aP3); + cPt2dr aP1 = aP0; + aCalcSBlur.RefineSadlePointFromIm(aImBlur,aP1,true); + cPt2dr aP2 = aP1; + aCalcSInit.RefineSadlePointFromIm(aImBlur,aP2,true); + + // tREAL8 aSymInt = aFSym->ComputeValMaxCrown(aP2,10); + // tREAL8 aSymInt = OptimFilter(aFSym,aP2).ValExtre(); + + tREAL8 aSymInt = aOptimSym.Optim(aP0,1.0,0.05).first; + + bool Ok = IsPtTest(Proj(aP3)); + + //StdOut() << "SYMM=" << aFSym->ComputeVal(aP0) << "\n"; + + if (Ok) + { + aSSad1.Add(aP3.z()); + aSSymInt_1.Add(aSymInt); + } + else + { + aSSad0.Add(aP3.z()); + aSSymInt_0.Add(aSymInt); + } + } + + if (mHasMasqTest) + { + StdOut() << " ================ STAT LOC CRITERIA ==================== \n"; + StdOut() << " * Saddle , #Ok Min=" << aSSad1.Min() + << " #NotOk 90% " << aSSad0.ErrAtProp(0.9) + << " 99% " << aSSad0.ErrAtProp(0.99) + << " 99.9% " << aSSad0.ErrAtProp(0.999) + << "\n"; + + StdOut() << " * SYM , #Ok=" << aSSymInt_1.Max() << " %75=" << aSSymInt_1.ErrAtProp(0.75) + << " NotOk 50% " << aSSymInt_0.ErrAtProp(0.5) + << " 10% " << aSSymInt_0.ErrAtProp(0.10) + << "\n"; + } + delete aFSym; + } + + cAutoTimerSegm aTSMakeIm(mTimeSegm,"OTHERS"); + + MakeImageSaddlePoints(*mDImIn,*mDImLabel); + StdOut() << "NBS " << (100.0*aNbSaddle)/aNbTot << " " << (100.0*aNbCCSad)/aNbTot + << " " << (100.0*aVMaxLoc.size())/aNbTot << " NB=" << aVMaxLoc.size() << "\n"; + aDImBlur.ToFile("Blurred.tif"); + } +/* + * Dist= sqrt(5) + * T=6.3751 + * 2 sqrt->6.42483 + */ int cAppliCheckBoardTargetExtract::Exe() { + mPhProj.FinishInit(); + if (RunMultiSet(0,0)) { return ResultMultiSet(); diff --git a/MMVII/src/CodedTarget/cFilterTarget.cpp b/MMVII/src/CodedTarget/cFilterTarget.cpp index 0b1fc7e535..5a8ed46f89 100755 --- a/MMVII/src/CodedTarget/cFilterTarget.cpp +++ b/MMVII/src/CodedTarget/cFilterTarget.cpp @@ -237,6 +237,13 @@ template double cFilterDCT::ComputeVal(const cPt2dr & aP) return Compute(); } + +template cPt1dr cFilterDCT::Value(const cPt2dr& aP) const +{ + return cPt1dr(const_cast&>(*this).ComputeVal(aP)); +} + + //static bool BUGF = false; template double cFilterDCT::ComputeValMaxCrown(const cPt2dr & aP,const double & aThreshold) @@ -687,6 +694,13 @@ template return new cSymFilterCT(anIm,aParam); } +template + cFilterDCT * cFilterDCT::AllocSym(tIm anIm,double aR0,double aR1,double aEpsilon) +{ + return new cSymFilterCT(anIm,aR0,aR1,aEpsilon); +} + + template cFilterDCT * cFilterDCT::AllocBin(tIm anIm,const cParamAllFilterDCT & aParam) { diff --git a/MMVII/src/CodedTarget/cSaddleFiter.cpp b/MMVII/src/CodedTarget/cSaddleFiter.cpp index f03025b9dc..03de5af33c 100755 --- a/MMVII/src/CodedTarget/cSaddleFiter.cpp +++ b/MMVII/src/CodedTarget/cSaddleFiter.cpp @@ -55,10 +55,13 @@ std::pair> cBasisFuncQuad::WeightAndVals(const cPt2dr cCompiledNeighBasisFunc::cCompiledNeighBasisFunc(const cBasisFunc & aBaseF ,const tVNeigh & aVNeigh) : - mVNeigh (aVNeigh), - mNbNeigh (mVNeigh.size()), - mSys (1) - + mVNeigh (aVNeigh), + mNbNeigh (mVNeigh.size()), + mSys (1), + mFisrtFast (true), + mMatInvCov (cPt2di(1,1)), + mRHS (1), + mSolFast (1) { MMVII_INTERNAL_ASSERT_tiny(!mVNeigh.empty(),"Emty neigh cCompiledNeighBasisFunc"); for (const auto & aNeigh : mVNeigh) @@ -85,9 +88,37 @@ cDenseVect cCompiledNeighBasisFunc::SlowCalc(const std::vector mSys.PublicAddObservation(mVWeight[aK],mVFuncs[aK],aVV[aK]); } - return mSys.Solve(); + cDenseVect aRes = mSys.Solve(); + + return aRes; +} + +const cDenseVect & cCompiledNeighBasisFunc::FastCalc(const std::vector & aVV) +{ + if (mFisrtFast) + { + mFisrtFast = false; + cLeasSqtAA aSys(mNbFunc); + for (size_t aK=0 ; aK cCompiledNeighBasisFunc::SlowCalc(const std::vector cCalcSaddle::cCalcSaddle(double aRay,double aStep) : mRay (aRay), mStep (aStep), - mVINeigh (SortedVectOfRadius(-1,aRay/aStep)), + mVINeigh (SortedVectOfRadius(-1,mRay/aStep)), mNbNeigh (mVINeigh.size()), - mRNeigh (ToR(mVINeigh,aStep)), + mRNeigh (ToR(mVINeigh,mStep)), mCalcQuad (cBasisFuncQuad(),mRNeigh), mVVals (mNbNeigh) { } +const std::vector & cCalcSaddle::VINeigh() const {return mVINeigh;} + tREAL8 cCalcSaddle::CalcSaddleCrit(const std::vector & aVVals,bool Show) { - cDenseVect aVect = mCalcQuad.SlowCalc(aVVals); + const cDenseVect & aVect = mCalcQuad.FastCalc(aVVals); + // cDenseVect aVect = mCalcQuad.SlowCalc(aVVals); tREAL8 aXX = aVect(0); tREAL8 aXY = aVect(1); tREAL8 aYY = aVect(2); + if (Show) StdOut() << " dxx=" << aXX << " dxy=" << aXY << " dyy=" << aYY << std::endl; // (A -X B) @@ -122,10 +157,22 @@ tREAL8 cCalcSaddle::CalcSaddleCrit(const std::vector & aVVals,bool Show tREAL8 aLapl = aXX+aYY; tREAL8 aDiscr= std::sqrt(Square(aXX-aYY) + 4*Square(aXY)); - tREAL8 aL1 = (aLapl+aDiscr)/2.0; - tREAL8 aL2 = (aLapl-aDiscr)/2.0; + tREAL8 aL1 = (aLapl+aDiscr)/2.0; // Highest Eigen Val + tREAL8 aL2 = (aLapl-aDiscr)/2.0; // Smallest Eigen Val + // aL1 aL2 = (aLapl^2 - aDiscr^2)/4 = (aXX^2 + aYY^2 + 2 aXX aYY -(aXX^2+aYY^2-aXX aYY+4 aXY^2)) + // = aXX aYY - aXY ^2 = det if (Show) StdOut() << " L1=" << aL1 << " L2=" << aL2 << std::endl; + + tREAL8 aRes = std::min(std::abs(aL1),std::abs(aL2)); + + if ((aL1>0)!=(aL2>0)) + return aRes; // we are saddle , min abs show how safe we are + else + return -aRes; // we are not saddle -min abs show how far we are + + + /* if (aL1>0) { aL2 = -aL2; @@ -137,12 +184,28 @@ tREAL8 cCalcSaddle::CalcSaddleCrit(const std::vector & aVVals,bool Show } return aL2; + */ +} + + +tREAL8 cCalcSaddle::CalcSaddleCrit(cDataIm2D& aDIm,cPt2di aPix) +{ + MMVII_INTERNAL_ASSERT_tiny(mStep==1.0,"Require step 1 for integer CalcSaddleCrit"); + + std::vector aVVals; + for (const auto & aDelta : mVINeigh) + { + aVVals.push_back(aDIm.GetV(aPix+aDelta)); + } + return CalcSaddleCrit(aVVals); } + + cPt2dr cCalcSaddle::RefineSadlPtFromVals(const std::vector & aVVals,bool Show) { - cDenseVect aVect = mCalcQuad.SlowCalc(aVVals); + const cDenseVect& aVect = mCalcQuad.FastCalc(aVVals); tREAL8 dXX = aVect(0); tREAL8 dXY = aVect(1); @@ -163,32 +226,44 @@ cPt2dr cCalcSaddle::RefineSadlPtFromVals(const std::vector & aVVals,bo return -cPt2dr(dYY*dX -dXY*dY,-dXY*dX +dXX*dY) / aDet; } -void cCalcSaddle::RefineSadlePointFromIm(cIm2D aIm,cDCT & aDCT) +bool cCalcSaddle::RefineSadlePointFromIm(cIm2D aIm,cPt2dr & aPt,bool ResetIfDivg) { double aThrDiv = 3.0; - cPt2dr aP0 = aDCT.mPt; + cPt2dr aP0 = aPt; cDataIm2D & aDIm = aIm.DIm(); for (int aK=0 ; aK<4 ; aK++) { for (size_t aKNeigh=0; aKNeigh aThrDiv) || (aDIm.Interiority(ToI(aDCT.mPt))<20)) + aPt += aDPt; + if ((Norm2(aPt-aP0) > aThrDiv) || (aDIm.Interiority(ToI(aPt))<20)) { - if (aDCT.mGT) - { - StdOut() << "DIVG " << aDCT.mGT->mC << " DPT " << aDPt << std::endl; - } - aDCT.mState = eResDCT::Divg; - return; + if (ResetIfDivg) + aPt = aP0; + return false; } if (Norm2(aDPt) < 1e-2) - return; + return true; + } + return true; +} + +void cCalcSaddle::RefineSadlePointFromIm(cIm2D aIm,cDCT & aDCT) +{ + bool Ok = RefineSadlePointFromIm(aIm,aDCT.mPt); + + if (!Ok) + { + if (aDCT.mGT) + { + StdOut() << "DIVG " << aDCT.mGT->mC << std::endl; + } + aDCT.mState = eResDCT::Divg; } } diff --git a/MMVII/src/Geoms/TilesIndexGeom.cpp b/MMVII/src/Geoms/TilesIndexGeom.cpp index a7343c156b..a513f0bb8d 100755 --- a/MMVII/src/Geoms/TilesIndexGeom.cpp +++ b/MMVII/src/Geoms/TilesIndexGeom.cpp @@ -40,6 +40,39 @@ template template class cGeneratePointDiff<2>; +/* **************************************** */ +/* */ +/* cTilingIndex */ +/* */ +/* **************************************** */ + +class cPoint2DValuated +{ + public : + static constexpr int Dim = 2; + typedef cPt2dr tPrimGeom; + typedef int tArgPG; /// unused here + + tPrimGeom GetPrimGeom(int Arg=-1) const {return Proj(mPt);} + + cPoint2DValuated(const cPt3dr & aPt) : + mPt (aPt) + { + } + + private : + cPt3dr mPt; +}; + + +std::vector FilterMaxLoc(std::vector & aVPt,tREAL8 aDist) +{ + SortOnCriteria( aVPt,[](const auto & aPt){return - aPt.z();} ); + + return FilterMaxLoc((cPt2dr*)nullptr,aVPt,[](const auto & aP) {return Proj(aP);}, aDist); +} + + /* **************************************** */ @@ -203,7 +236,7 @@ void OneBenchSpatialIndex() cBox2dr aBox(aP0,aP1); // Box of the tiling cBox2dr aBoxMargin(aP0-aSzMargin,aP1+aSzMargin); // box slightly bigger - cTiling> aSI(aBox,true,1000,-1); // The tiling we want to check + cTiling> aSI(aBox,true,1000,-1); // The tiling we want to check , -1 is the fake arge // Test the function GetObjAtPos std::list aLPt; diff --git a/MMVII/src/Mappings/MappingGen.cpp b/MMVII/src/Mappings/MappingGen.cpp index d1a82bb2f2..2ac6eb38df 100755 --- a/MMVII/src/Mappings/MappingGen.cpp +++ b/MMVII/src/Mappings/MappingGen.cpp @@ -509,6 +509,7 @@ INSTANCE_TWO_DIM_MAPPING(DIM,2);\ INSTANCE_TWO_DIM_MAPPING(DIM,3); INSTANCE_TWO_DIM_MAPPING(1,1); +INSTANCE_TWO_DIM_MAPPING(2,1); INSTANCE_ONE_DIM_MAPPING(1) INSTANCE_ONE_DIM_MAPPING(2) diff --git a/MMVII/src/Matrix/BaseMatrixes.cpp b/MMVII/src/Matrix/BaseMatrixes.cpp index c91fb56649..29bab51242 100755 --- a/MMVII/src/Matrix/BaseMatrixes.cpp +++ b/MMVII/src/Matrix/BaseMatrixes.cpp @@ -284,6 +284,10 @@ template void cDenseVect::WeightedAddIn(Type aW,const tSpV & aD[aP.mInd] += aW * aP.mVal; } +template void cDenseVect::WeightedAddIn(Type aW,const tDV & aVect) +{ + MMVII::WeightedAddIn(mIm.DIm(),aW,aVect.mIm.DIm()); +} template std::ostream & operator << (std::ostream & OS,const cDenseVect &aV) { diff --git a/MMVII/src/MorphoMat/ConnectedComponents.cpp b/MMVII/src/MorphoMat/ConnectedComponents.cpp new file mode 100644 index 0000000000..e3558dbe65 --- /dev/null +++ b/MMVII/src/MorphoMat/ConnectedComponents.cpp @@ -0,0 +1,64 @@ +#include "MMVII_ImageMorphoMath.h" +namespace MMVII +{ + +void ConnectedComponent + ( + std::vector & aVPts, + cDataIm2D & aDIm, + const std::vector & aNeighbourhood, + const std::vector & aVecSeed, + int aMarqInit, + int aNewMarq, + bool Restore + ) +{ + aVPts.clear(); + for (const auto & aPSeed : aVecSeed) + { + if (aDIm.GetV(aPSeed) == aMarqInit) + { + aDIm.SetV(aPSeed,aNewMarq); + aVPts.push_back(aPSeed); + } + } + size_t aIndBottom = 0; + + while (aIndBottom!=aVPts.size()) + { + cPt2di aP0 = aVPts[aIndBottom]; + for (const auto & aDelta : aNeighbourhood) + { + cPt2di aNeigh = aP0 + aDelta; + if (aDIm.GetV(aNeigh)==aMarqInit) + { + aDIm.SetV(aNeigh,aNewMarq); + aVPts.push_back(aNeigh); + } + } + aIndBottom++; + } + + if (Restore) + { + for (const auto & aPix : aVPts) + aDIm.SetV(aPix,aMarqInit); + } +} + +void ConnectedComponent + ( + std::vector & aVPts, + cDataIm2D & aDIm, + const std::vector & aNeighbourhood, + const cPt2di& aPtSeed, + int aMarqInit, + int aNewMarq, + bool Restore + ) +{ + std::vector aVSeed {aPtSeed}; + ConnectedComponent(aVPts,aDIm,aNeighbourhood,aVSeed,aMarqInit,aNewMarq,Restore); +} + +}; diff --git a/MMVII/src/TieP/cToTiePMultiple.cpp b/MMVII/src/TieP/cToTiePMultiple.cpp index ce08a030b3..f93a1a4212 100644 --- a/MMVII/src/TieP/cToTiePMultiple.cpp +++ b/MMVII/src/TieP/cToTiePMultiple.cpp @@ -555,7 +555,7 @@ void cOneImMEff2MP::ComputeIndexPts(cInterfImportHom & anImport,const cMemoryEf { public : // static constexpr int Dim = 2; => dont work with local class, hapilly enum works - enum {Dim=2}; + enum {TheDim=2}; typedef cPt2dr tPrimGeom; // geometric primitives indexed are points // type of arg that we will used in call back "GetPrimGeom", we need to refer to images typedef cOneImMEff2MP * tArgPG; diff --git a/MMVII/src/UtiMaths/cOptimByStep.cpp b/MMVII/src/UtiMaths/cOptimByStep.cpp new file mode 100644 index 0000000000..0fd32d86d6 --- /dev/null +++ b/MMVII/src/UtiMaths/cOptimByStep.cpp @@ -0,0 +1,135 @@ +#include "MMVII_HeuristikOpt.h" +#include "MMVII_Mappings.h" + +namespace MMVII +{ + + +/** Given a dim of point, and a Dist1Max corresponding to a neighboor (for ex : 1 for V4, 2 for V8 with Dim=2 ...), + * let P1 in a neighboord, we want to have the indexe of point P2 such that P1+P2 is not in a neighboorhoud, + * + * It is used in "combinatory" optim, once we have the best point P1 in the neighboor, we want to explorate the points P2, + * neighboor of P1 that , have not been already explored at previus step (i.e were not in the neighboorhoud). + * + * For ex, with Dim=2, DistMax= 2 , neigboorhood is something like (not sure numbering) + * + * 3 2 1 + * 4 x 0 + * 5 6 7 + * + * Res[0] = (1,0,7) + * Res[1] = (0,1,2,3,7) + * Also conventionnaly we add at the end a special entry with all neighboors + * + * Res[8] = (0,1,2 ... 7) + * + */ + + +template const tVLInt & VLNeighOptim(int aDist1Max) +{ + static std::map TheMap; // Buffer of result + + tVLInt & aVInt = TheMap[aDist1Max]; + if (aVInt.empty()) // if was not alerady computed + { + // tVLInt & aVInt = TheMap[aDist1Max]; + const auto & aVNeigh = AllocNeighbourhood(aDist1Max); + aVInt.resize(aVNeigh.size()+1); + + for (size_t aK1=0 ; aK1 & aP1 = aVNeigh.at(aK1); + for (size_t aK2=0 ; aK2 aV12 = aP1 + aVNeigh.at(aK2); + // P in neigboors such : NormInf(P) <=1 && Norm1(P) <= aDist1Max + if ( (NormInf(aV12) >1) || (Norm1(aV12)>aDist1Max) ) + aVInt.at(aK1).push_back(aK2); + } + aVInt.at(aVNeigh.size()).push_back(aK1); // at the end we push all the neighbours + } + } + + return aVInt; +} + +/* ******************************************* */ +/* */ +/* cOptimByStep */ +/* */ +/* ******************************************* */ + + + +template tREAL8 cOptimByStep::Value(const tPtR & aPt) const +{ + return mSign * mFunc.Value(aPt).x(); +} + +template cOptimByStep::cOptimByStep(const tMap & aMap,bool IsMin,tREAL8 aMaxDInfInit,int aDist1Max) : + mFunc (aMap), + mSign (IsMin ? 1.0 : -1.0), + mWMin (), + mMaxDInfInit (aMaxDInfInit), + mDist1Max (aDist1Max), + mINeigh (AllocNeighbourhood(aDist1Max)), + mIndNextN (VLNeighOptim(aDist1Max)) +{ + VLNeighOptim(aDist1Max); +} + +template + std::pair> cOptimByStep::Optim(const tPtR & aP0 ,tREAL8 aStepInit,tREAL8 aStepLim,tREAL8 aMul) +{ + mPt0 = aP0; + mWMin = cWhichMin(aP0,Value(aP0)); + + for (tREAL8 aStep = aStepInit ; aStep>= aStepLim ; aStep *= aMul) + { + if (! DoOneStep(aStep) ) + return CurValue(); + } + + return CurValue(); +} + +template std::pair> cOptimByStep::CurValue() const +{ + return std::pair(mWMin.ValExtre()*mSign,mWMin.IndexExtre()); +} +template bool cOptimByStep::DoOneStep(tREAL8 aStep) +{ + int aLastN = mINeigh.size(); // initially set to all neighboors + for (;;) // whil we can ameliorate + { + tPtR aPOpt = mWMin.IndexExtre(); // current value of optimal point + int aNewLast = -1; // by default no new index + for (const auto & aKN : mIndNextN.at(aLastN)) // parse all indexes not already explored + { + tPtR aNewP = aPOpt + ToR(mINeigh.at(aKN)) * aStep; // new point to test + bool aNewBest = mWMin.Add(aNewP,Value(aNewP)); // update optimum + if (aNewBest) // if some update was made + { + aNewLast = aKN; // the new best last indexe + // precaution usefull with image processing + if (NormInf(mPt0-mWMin.IndexExtre()) > mMaxDInfInit) + return false; + } + } + // if no update, finish w/o problem + if (aNewLast<0) + return true; + aLastN = aNewLast; + } +} + + +template class cOptimByStep<1>; +template class cOptimByStep<2>; +template class cOptimByStep<3>; + + + +}; + From 73d39ce901cca8fac7f89ab142bfca7e8299c4b8 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Tue, 4 Jun 2024 20:12:48 +0200 Subject: [PATCH 3/5] Forgoten file --- MMVII/src/MorphoMat/FlagBitsOperation.cpp | 109 ++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 MMVII/src/MorphoMat/FlagBitsOperation.cpp diff --git a/MMVII/src/MorphoMat/FlagBitsOperation.cpp b/MMVII/src/MorphoMat/FlagBitsOperation.cpp new file mode 100644 index 0000000000..7e01b74e0d --- /dev/null +++ b/MMVII/src/MorphoMat/FlagBitsOperation.cpp @@ -0,0 +1,109 @@ +#include "MMVII_ImageMorphoMath.h" + +namespace MMVII +{ + + +/* ******************************************************** */ +/* */ +/* cFlag8Neigh */ +/* */ +/* ******************************************************** */ + +void cFlag8Neigh::ComputeCCOfeNeigh(std::vector> & aVCC) +{ + aVCC = std::vector> (256); + + for (size_t aFlag=0 ; aFlag<256 ; aFlag++) // parse all possible flag combinaisons + { + cFlag8Neigh aNeigh(aFlag); + for (size_t aBit=0 ; aBit<8 ; aBit++) // parse all bit + { + if (aNeigh.SafeIsIn(aBit) != aNeigh.SafeIsIn(aBit-1)) // select bit of transition + { + int aBitNext = aBit+1; + while (aNeigh.SafeIsIn(aBit) == aNeigh.SafeIsIn(aBitNext)) // go to next transition + aBitNext++; + + // create a new connected component + cCCOfOrdNeigh aCC; + aCC.mIs1 = aNeigh.SafeIsIn(aBit); + aCC.mFirstBit = aBit; + aCC.mLastBit = aBitNext; + // add it + aVCC.at(aFlag).push_back(aCC); + } + } + } +} + +/* ******************************************************** */ +/* */ +/* :: Global func, based on value ordering & neigboors */ +/* */ +/* ******************************************************** */ + +template cFlag8Neigh FlagSup8Neigh(const cDataIm2D & aDIm,const cPt2di & aPt) +{ + Type aV0 = aDIm.GetV(aPt); + cFlag8Neigh aResult ; + // for freeman 1,4 Y is positive, for freeman 0, X is positive + for (int aK=0 ; aK<4 ; aK++) + { + if (aDIm.GetV(aPt+FreemanV8[aK]) >= aV0) + aResult.AddNeigh(aK); + } + for (int aK=4 ; aK<8 ; aK++) + { + if (aDIm.GetV(aPt+FreemanV8[aK]) > aV0) + aResult.AddNeigh(aK); + } + return aResult; +} + +/* For this "topo" criteria, we compute : + * + * - aMaxOfMin = max for all connected component of the min of the connected component + * typically we should have 2 CC, they should be both low + * - idem aMinOfMax + * + * Value(Point) - aMaxOfMin : reflect it is a sadlle point for low value + * aMinOfMax - Value(Point) : idem + * + * we select the worst criteria + */ + + + +template tREAL8 CriterionTopoSadle(const cDataIm2D & aDIm,const cPt2di & aPixC) +{ + cFlag8Neigh aNeigh = FlagSup8Neigh(aDIm,aPixC); + const std::vector & aVCC = aNeigh.ConComp() ; + + tREAL8 aMaxOfMin(-1e10); + tREAL8 aMinOfMax(1e10); + for (const auto & aCC : aVCC) + { + cBoundVals aBounds; + for (int aBit = aCC.mFirstBit ; aBit & aDIm,const cPt2di & aPt); +template tREAL8 CriterionTopoSadle(const cDataIm2D & aDIm,const cPt2di & aPixC); + + + +}; From f0cc142bae44835f741b877c3ce4f86613cd4518 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Fri, 7 Jun 2024 18:46:47 +0200 Subject: [PATCH 4/5] In checkboard , radiom modelization --- MMVII/include/MMVII_ExtractLines.h | 37 ++ MMVII/include/MMVII_Image2D.h | 2 + MMVII/include/MMVII_Mappings.h | 1 + MMVII/include/MMVII_Matrix.h | 4 + MMVII/include/MMVII_nums.h | 3 + .../CodedTarget/cCheckBoardTargetExtract.cpp | 427 ++++++++++++++++-- MMVII/src/ImagesBase/cRGBImage.cpp | 14 + .../src/ImagesInfoExtract/cScoreTetaLine.cpp | 128 ++++++ MMVII/src/UtiMaths/uti_stat.cpp | 33 ++ 9 files changed, 605 insertions(+), 44 deletions(-) create mode 100755 MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp diff --git a/MMVII/include/MMVII_ExtractLines.h b/MMVII/include/MMVII_ExtractLines.h index cfc57981b9..95680a4599 100755 --- a/MMVII/include/MMVII_ExtractLines.h +++ b/MMVII/include/MMVII_ExtractLines.h @@ -251,6 +251,43 @@ template class cExtractLines std::vector mPtsCont; ///< List of point in mImMasqCont }; + +/** Class for extracting a line arround a point, for now its very specialized in the extraction of + * direction in checkboard, may evolve to things more general ? + */ + +class cScoreTetaLine : public tFunc1DReal // herit from tFunc1DReal for optimization in cOptimByStep +{ + public : + /// constructor : aL length of the segment, aStep step of discretization in a segment + cScoreTetaLine(cDataIm2D &,const cDiffInterpolator1D & ,tREAL8 aL,tREAL8 aStep); + + /// extract the 2 angle of line in checkboar, aStepInit & aStepLim => used in cOptimByStep + std::pair Tetas_CheckBoard(const cPt2dr& aC,tREAL8 aStepInit,tREAL8 aStepLim); + private : + + /// fix center of reusing the data (to avoid cost for cTabulatedDiffInterpolator) + void SetCenter(const cPt2dr & aC); + + /// Extract the initial values of tetas by testing all at a given step (in pixel) + tREAL8 GetTetasInit(tREAL8 aStepPix,int aCurSign); + /// Refine the value with a limit step in pixel + tREAL8 Refine(tREAL8 aTeta0,tREAL8 aStepPix,int aSign); + + /// Value as tFunc1DReal + cPt1dr Value(const cPt1dr& aPt) const override; + + cDataIm2D * mDIm; + tREAL8 mLength; + int mNb; + tREAL8 mStepAbsc; + tREAL8 mStepTeta; + cTabulatedDiffInterpolator mTabInt; + cPt2dr mC; + tREAL8 mCurSign; /// used to specify an orientaion of segment +}; + + }; #endif // _MMVII_EXTRACT_LINES_H_ // diff --git a/MMVII/include/MMVII_Image2D.h b/MMVII/include/MMVII_Image2D.h index 04fb5177c1..df639b64d9 100755 --- a/MMVII/include/MMVII_Image2D.h +++ b/MMVII/include/MMVII_Image2D.h @@ -499,6 +499,8 @@ class cRGBImage void Write(const std::string &,const cPt2di & aP0,double aDyn=1,const cRect2& =cRect2::TheEmptyBox) const; // 1 to 1 /* */ + // transformate the RGB internal image in gray + void ResetGray(); /// set values iff param are OK, RGB image are made for visu, not for intensive computation void SetRGBPix(const cPt2di & aPix,int aR,int aG,int aB); diff --git a/MMVII/include/MMVII_Mappings.h b/MMVII/include/MMVII_Mappings.h index afdc45116e..288fbd1d44 100755 --- a/MMVII/include/MMVII_Mappings.h +++ b/MMVII/include/MMVII_Mappings.h @@ -373,6 +373,7 @@ template class cDataMapping : publ #endif // MAP_STATIC_BUF }; +typedef cDataMapping tFunc1DReal; typedef cDataMapping tFunc2DReal; /** Specialization for DimIn=DimOut , introduce because we want to force invertible mapping diff --git a/MMVII/include/MMVII_Matrix.h b/MMVII/include/MMVII_Matrix.h index 5f01b285a4..a14f1ab7e1 100755 --- a/MMVII/include/MMVII_Matrix.h +++ b/MMVII/include/MMVII_Matrix.h @@ -728,6 +728,9 @@ template class cMatIner2Var Type CorrelNotC(const Type &aEpsilon=1e-10) const; // Non centered correl Type StdDev1() const; Type StdDev2() const; + + /// [A B] as least-square solution of V1 = A V2 + B + std::pair FitLineDirect() const; private : Type mS0; ///< Som of W Type mS1; ///< Som of W * V1 @@ -744,6 +747,7 @@ template class cWeightAv cWeightAv(); void Add(const TypeWeight & aWeight,const TypeVal & aVal); TypeVal Average() const; + TypeVal Average(const TypeVal & aDef) const; const TypeVal & SVW() const; /// Accessor to sum weighted vals const TypeWeight & SW() const; /// Accessor to sum weighted vals private : diff --git a/MMVII/include/MMVII_nums.h b/MMVII/include/MMVII_nums.h index e43a81bff6..6d3168dbc2 100755 --- a/MMVII/include/MMVII_nums.h +++ b/MMVII/include/MMVII_nums.h @@ -725,6 +725,9 @@ template class cWhichMinMax const cWhichMin & Min() const {return mMin;} const cWhichMax & Max() const {return mMax;} + const TypeIndex & IndMin() const {return mMin.IndexExtre();} + const TypeIndex & IndMax() const {return mMax.IndexExtre();} + private : cWhichMin mMin; cWhichMax mMax; diff --git a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp index 19b30f622f..f34711fb4b 100755 --- a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp @@ -7,20 +7,272 @@ #include "MMVII_2Include_Tiling.h" #include "MMVII_Sensor.h" #include "MMVII_HeuristikOpt.h" +#include "MMVII_ExtractLines.h" -// Test git branch - namespace MMVII { +class cOptimPosSeg : public tFunc2DReal +{ + public : + cOptimPosSeg(const tSeg2dr &); + cPt1dr Value(const cPt2dr &) const override; + virtual tREAL8 CostOfSeg(const cSegment2DCompiled &) const = 0; + + cSegment2DCompiled ModifiedSeg(const cPt2dr & aPModif) const; + + tSeg2dr OptimizeSeg(tREAL8 aStepInit,tREAL8 aStepLim,bool IsMin,tREAL8 aMaxDInfInit) const; + + private : + cSegment2DCompiled mSeg0; + cPt2dr mP0Loc; + cPt2dr mP1Loc; +}; + + +cOptimPosSeg::cOptimPosSeg(const tSeg2dr & aSeg0) : + mSeg0 (aSeg0), + mP0Loc (mSeg0.ToCoordLoc(mSeg0.P1())), + mP1Loc (mSeg0.ToCoordLoc(mSeg0.P2())) +{ +} + +cSegment2DCompiled cOptimPosSeg::ModifiedSeg(const cPt2dr & aPModif) const +{ + cPt2dr aP0Glob = mSeg0.FromCoordLoc(mP0Loc+cPt2dr(0,aPModif.x())); + cPt2dr aP1Glob = mSeg0.FromCoordLoc(mP1Loc+cPt2dr(0,aPModif.y())); + + return cSegment2DCompiled(aP0Glob,aP1Glob); +} + +cPt1dr cOptimPosSeg::Value(const cPt2dr & aPt) const +{ + return cPt1dr(CostOfSeg(ModifiedSeg(aPt))); +} + +tSeg2dr cOptimPosSeg::OptimizeSeg(tREAL8 aStepInit,tREAL8 aStepLim,bool IsMin,tREAL8 aMaxDInfInit) const +{ + cOptimByStep<2> aOpt(*this,true,IsMin,aMaxDInfInit); + auto [aScore,aSol] = aOpt.Optim(cPt2dr(0,0),aStepInit,aStepLim); + + return ModifiedSeg(aSol); +} + +class cOptimSeg_ValueIm : public cOptimPosSeg +{ + public : + cOptimSeg_ValueIm(const tSeg2dr &,tREAL8 aStepOnSeg,const cDataIm2D & aDIm,tREAL8 aTargetValue); + + tREAL8 CostOfSeg(const cSegment2DCompiled &) const override; + private : + tREAL8 mStepOnSeg; + int mNbOnSeg; + const cDataIm2D & mDataIm; + tREAL8 mTargetValue; +}; + +tREAL8 cOptimSeg_ValueIm::CostOfSeg(const cSegment2DCompiled & aSeg) const +{ + tREAL8 aSum=0; + + for (int aK=0 ; aK<= mNbOnSeg ; aK++) + { + cPt2dr aPt = Centroid(aK/double(mNbOnSeg),aSeg.P1(),aSeg.P2()); + aSum += std::abs(mDataIm.GetVBL(aPt) - mTargetValue); + } + + return aSum / mNbOnSeg; +} + + + + static constexpr tU_INT1 eNone = 0 ; static constexpr tU_INT1 eTopo0 = 1 ; static constexpr tU_INT1 eTopoTmpCC = 2 ; static constexpr tU_INT1 eTopoMaxOfCC = 3 ; static constexpr tU_INT1 eTopoMaxLoc = 4 ; static constexpr tU_INT1 eFilterSym = 5 ; +static constexpr tU_INT1 eFilterRadiom = 6 ; + +/// candidate that are pre-selected on the sadle-point criteria +class cCdSadle +{ + public : + cCdSadle (const cPt2dr & aC,tREAL8 aCrit) : mC (aC) , mSadCrit (aCrit) {} + cCdSadle (){} + cPt2dr mC; + + /// Criterion of sadle point, obtain by fiting a quadric on gray level + tREAL8 mSadCrit; +}; + +/// candidate that are pre-selected after symetry criterion +class cCdSym : public cCdSadle +{ + public : + cCdSym(const cCdSadle & aCdSad,tREAL8 aCrit) : cCdSadle (aCdSad), mSymCrit (aCrit) { } + /// Criterion of symetry, obtain after optimisation of center + tREAL8 mSymCrit; + +}; + +/// candid obtain after radiometric modelization, + +class cCdRadiom : public cCdSym +{ + public : + cCdRadiom(const cCdSym &,const cDataIm2D & aDIm,tREAL8 aTeta1,tREAL8 aTeta2,tREAL8 aThickness); + + cMatIner2Var StatGray(const cDataIm2D & aDIm,tREAL8 aThickness,bool IncludeSegs); + + tREAL8 mTeta0; + tREAL8 mTeta1; + + tREAL8 mCostCorrel; // 1-Correlation of model + tREAL8 mBlack; + tREAL8 mWhite; +}; + +enum class eTPosCB +{ + eUndef, + eInsideBlack, + eInsideWhite, + eBorderLeft, + eBorderRight + +}; + +inline bool IsInside(eTPosCB aState) {return (aState==eTPosCB::eInsideBlack) || (aState==eTPosCB::eInsideWhite) ;} +inline bool IsOk(eTPosCB aState) {return aState!=eTPosCB::eUndef;} + +/// Used temporary for compilation of radiom +class cCdRadiomCompiled : public cCdRadiom +{ + public : + cCdRadiomCompiled(const cCdRadiom &,tREAL8 aThickness); + + /// Theoreticall radiom of modelize checkboard + bool if was computed + std::pair TheorRadiom(const cPt2dr &) const; + + tREAL8 mThickness; + cSegment2DCompiled mSeg0 ; + cSegment2DCompiled mSeg1 ; +}; + +cCdRadiomCompiled::cCdRadiomCompiled(const cCdRadiom & aCDR,tREAL8 aThickness) : + cCdRadiom (aCDR), + mThickness (aThickness), + mSeg0 (mC,mC+FromPolar(1.0,mTeta0)), + mSeg1 (mC,mC+FromPolar(1.0,mTeta1)) +{ +} + +std::pair cCdRadiomCompiled::TheorRadiom(const cPt2dr &aPt) const +{ + eTPosCB aPos = eTPosCB::eUndef; + tREAL8 aGrayTh = -1; + + // we compute locacl coordinates because the sign of y indicate if we are left/right of the oriented segment + // and sign of x indicate if we are before/after the centre + cPt2dr aLoc0 = mSeg0.ToCoordLoc(aPt); + tREAL8 aY0 = aLoc0.y(); + + cPt2dr aLoc1 = mSeg1.ToCoordLoc(aPt); + tREAL8 aY1 = aLoc1.y(); + + // compute if we are far enough of S0/S1 because the computation of gray will change + // black/white if far enough, else interpolation + bool FarS0 = std::abs(aY0)> mThickness; + bool FarS1 = std::abs(aY1)> mThickness; + + if ( FarS0 && FarS1) + { + if ((aY0>0)!=(aY1>0)) + { + aPos = eTPosCB::eInsideBlack; + aGrayTh = 0.0; + } + else + { + aPos = eTPosCB::eInsideWhite; + aGrayTh = 1.0; + } + } + else if ((!FarS0) && FarS1) + { + // (! FarS0) => teta1 + // Red = teta1 , black on left on image, right on left in coord oriented + aPos = eTPosCB::eBorderRight; + int aSignX = (aLoc0.x() >0) ? -1 : 1; + aGrayTh = (mThickness+aSignX*aY0) / (2.0*mThickness); + } + else if (FarS0 && (!FarS1)) + { + aPos = eTPosCB::eBorderLeft; + int aSignX = (aLoc1.x() <0) ? -1 : 1; + aGrayTh = (mThickness+aSignX*aY1) / (2.0 * mThickness); + } + + return std::pair(aPos,aGrayTh); +} + +cCdRadiom::cCdRadiom(const cCdSym & aCdSym,const cDataIm2D & aDIm,tREAL8 aTeta0,tREAL8 aTeta1,tREAL8 aThickness) : + cCdSym (aCdSym), + mTeta0 (aTeta0), + mTeta1 (aTeta1), + mCostCorrel (2.001) // over maximal theoreticall value +{ + static int aCpt=0 ; aCpt++; + + + cSegment2DCompiled aSeg0 (mC,mC+FromPolar(1.0,mTeta0)); + cSegment2DCompiled aSeg1 (mC,mC+FromPolar(1.0,mTeta1)); + static std::vector aDisk = VectOfRadius(0.0,5); + cStdStatRes aW0; + cStdStatRes aW1; + + int aNbIn0=0,aNbIn1=0; + + cMatIner2Var aCorGrayAll; + cMatIner2Var aCorGrayInside; + + cCdRadiomCompiled aCRC(*this,aThickness); + + for (const auto & aDelta : aDisk) + { + cPt2di aPImI = aDelta + ToI(mC); + tREAL8 aValIm = aDIm.GetV(aPImI); + cPt2dr aPImR = ToR(aPImI); + + auto [aState,aGrayTh] = aCRC.TheorRadiom(aPImR); + + if (IsInside(aState)) + { + aCorGrayInside.Add(aGrayTh,aValIm); + aNbIn0 += (aState == eTPosCB::eInsideBlack); + aNbIn1 += (aState == eTPosCB::eInsideWhite); + } + if (IsOk(aState)) + { + aCorGrayAll.Add(aGrayTh,aValIm); + } + } + + tREAL8 aRatioNb = std::min(aNbIn0,aNbIn1) / (tREAL8) std::max(aNbIn0,aNbIn1); + if (aRatioNb <0.2) + { + return ; + } + + mCostCorrel = 1-aCorGrayAll.Correl(); + auto [a,b] = aCorGrayInside.FitLineDirect(); + mBlack = b ; + mWhite = a+b; +} /* *********************************************************** */ /* */ @@ -28,6 +280,7 @@ static constexpr tU_INT1 eFilterSym = 5 ; /* */ /* *********************************************************** */ +class cScoreTetaLine; class cAppliCheckBoardTargetExtract : public cMMVII_Appli { @@ -52,11 +305,11 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli void DoOneImage() ; void MakeImageSaddlePoints(const tDIm &,const cDataIm2D & aDMasq) const; - cWhichMin OptimFilter(cFilterDCT *,const cPt2dr &) const; - cPhotogrammetricProject mPhProj; cTimerSegm mTimeSegm; + cCdRadiom TestBinarization(cScoreTetaLine&,const cCdSym &,tREAL8 aThickness); + // =========== Mandatory args ============ std::string mNameIm; ///< Name of background image @@ -66,7 +319,7 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli // -- - std::vector mShowQuickSaddleF; + tREAL8 mThickness; ///< used for fine estimation of radiom // =========== Internal param ============ tIm mImIn; ///< Input global image @@ -89,6 +342,7 @@ cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector cAppliCheckBoardTargetExtract::OptimFilter(cFilterDCT * aFilter,const cPt2dr & aP0) const + + +cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL,const cCdSym & aCdSym,tREAL8 aThickness) { - cWhichMin aRes; - for (tREAL8 aDx=-1.0 ; aDx<1.0 ; aDx += 0.1) - { - for (tREAL8 aDy=-1.0 ; aDy<1.0 ; aDy += 0.1) - { - cPt2dr aPt = aP0+cPt2dr(aDx,aDy); - // tREAL8 aVal = aFilter->ComputeValMaxCrown(aPt,1e10); - tREAL8 aVal = aFilter->ComputeVal(aPt); - aRes.Add(aPt,aVal); - // UpdateMin(aMinSym,aFSym->ComputeValMaxCrown(aP2+cPt2dr(aDx,aDy),10)); - } - } - return aRes; + + auto [aTeta0,aTeta1] = aSTL.Tetas_CheckBoard(aCdSym.mC,0.5,1e-3); + + cCdRadiom aCdRadiom(aCdSym,*mDImIn,aTeta0,aTeta1,aThickness); + + if (IsPtTest(aCdSym.mC)) + { + static int aCpt=0 ; aCpt++; + StdOut() << " CPT=" << aCpt << " Corrrr=" << aCdRadiom.mCostCorrel << " V0="<< aCdRadiom.mBlack << " V1=" << aCdRadiom.mWhite << "\n"; + // StdOut() << " CPT=" << aCpt << " TETASRef= " << aTeta0 << " " << aTeta1 << "\n"; + + int aZoom = 9; + cPt2di aSz(50,50); + cPt2di aDec = ToI(aCdSym.mC) - aSz/2; + cPt2dr aCLoc = aCdSym.mC-ToR(aDec); + + cRGBImage aIm = cRGBImage:: FromFile(mNameIm,cBox2di(aDec,aDec+aSz),aZoom); + aIm.ResetGray(); + + if (0) + { + cCdRadiomCompiled aCRC(aCdRadiom,2.0); + for (const auto & aPix : aIm.ImR().DIm()) + { + cPt2dr aPtR = ToR(aPix+aDec); + auto [aState,aGray] = aCRC.TheorRadiom(aPtR); + if (aState != eTPosCB::eUndef) + { + aIm.SetGrayPix(aPix,aGray*255); + } + } + } + + int aKT=0; + for (const auto & aTeta : {aTeta0,aTeta1}) + { + for (int aK= -aZoom * 20 ; aK<=aZoom*20 ; aK++) + { + tREAL8 aAbsc= aK/ (2.0 * aZoom); + cPt2dr aPt = aCLoc + FromPolar(aAbsc,aTeta); + aIm.SetRGBPoint(aPt,cRGBImage::Green); + if (aK==6*aZoom) + aIm.DrawCircle((aKT==0) ? cRGBImage::Blue : cRGBImage::Red,aPt,0.5); + // aIm.SetRGBPoint(aPt, (aKT==0) ? cRGBImage::Blue : cRGBImage::Red); + } + aKT++; + } + + aIm.SetRGBPoint(aCLoc,cRGBImage::Red); + aIm.ToFile("TestCenter_" + ToStr(aCpt) + ".tif"); + } + + return aCdRadiom; } void cAppliCheckBoardTargetExtract::DoOneImage() @@ -162,10 +458,13 @@ void cAppliCheckBoardTargetExtract::DoOneImage() int mMaxNbMLS = 2000; // Max number of point in best saddle points tREAL8 aRayCalcSadle = sqrt(4+1); // limit point 2,1 - tREAL8 mThresholdSym = 0.50; - tREAL8 mDistCalcSym0 = 8.0; - tREAL8 mDistDivSym = 2.0; + tREAL8 mThresholdSym = 0.50; // threshlod for symetry criteria + tREAL8 mDistCalcSym0 = 8.0; // distance for evaluating symetry criteria + tREAL8 mDistDivSym = 2.0; // maximal distance to initial value in symetry opt + tREAL8 mLengtSInit = 5.0; + tREAL8 mStepSeg = 0.5; + tREAL8 mMaxCostCorrIm = 0.1; // computed threshold tINT8 mDistRectInt = 20; // to see later how we compute it @@ -251,74 +550,103 @@ void cAppliCheckBoardTargetExtract::DoOneImage() /* [3] Compute point that are max local */ + + std::vector aVCdtSad; cAutoTimerSegm aTSCritSad(mTimeSegm,"CritSad"); - std::vector aVSad0; cCalcSaddle aCalcSBlur(aRayCalcSadle+0.001,1.0); + // [3.1] compute for each point the saddle criteria for (const auto& aPix : *mDImLabel) { if (mDImLabel->GetV(aPix)==eTopoMaxOfCC) { tREAL8 aCritS = aCalcSBlur.CalcSaddleCrit(aDImBlur,aPix); - aVSad0.push_back(cPt3dr(aPix.x(),aPix.y(),aCritS)); + aVCdtSad.push_back(cCdSadle(ToR(aPix),aCritS)); +// cPt3dr(aPix.x(),aPix.y(),aCritS)); } } + int aNbSad0 = aVCdtSad.size(); // [3.2] select KBest + MaxLocal cAutoTimerSegm aTSMaxLoc(mTimeSegm,"MaxLoc"); - SortOnCriteria(aVSad0,[](const auto & aPt){return - aPt.z();} ); - std::vector aVMaxLoc = FilterMaxLoc((cPt2dr*)nullptr,aVSad0,[](const auto & aP) {return Proj(aP);}, mDistMaxLocSad); + // sort by decreasing criteria of saddles => "-" aCdt.mSadCrit + SortOnCriteria(aVCdtSad,[](const auto & aCdt){return - aCdt.mSadCrit;}); + aVCdtSad = FilterMaxLoc((cPt2dr*)nullptr,aVCdtSad,[](const auto & aCdt) {return aCdt.mC;}, mDistMaxLocSad); + int aNbSad1 = aVCdtSad.size(); // limit the number of point , a bit rough but first experiment show that sadle criterion is almost perfect on good images - aVMaxLoc.resize(std::min(aVMaxLoc.size(),size_t(mMaxNbMLS))); + aVCdtSad.resize(std::min(aVCdtSad.size(),size_t(mMaxNbMLS))); - for (const auto & aP3 : aVMaxLoc) - mDImLabel->SetV(ToI(Proj(aP3)),eTopoMaxLoc); + for (const auto & aCdt : aVCdtSad) + mDImLabel->SetV(ToI(aCdt.mC),eTopoMaxLoc); + int aNbSad2 = aVCdtSad.size(); StdOut() << "END MAXLOC \n"; - /* [4] Calc Symetry criterion */ + std::vector aVCdtSym; cAutoTimerSegm aTSSym(mTimeSegm,"SYM"); - if (1) { cCalcSaddle aCalcSInit(1.5,1.0); - std::vector aNewP3; cFilterDCT * aFSym = cFilterDCT::AllocSym(mImIn,0.0,mDistCalcSym0,1.0); cOptimByStep aOptimSym(*aFSym,true,mDistDivSym); - for (const auto & aP3 : aVMaxLoc) - { - cPt2dr aP0 = Proj(aP3); - auto [aValSym,aNewP] = aOptimSym.Optim(aP0,1.0,0.01); + for (auto & aCdtSad : aVCdtSad) + { + auto [aValSym,aNewP] = aOptimSym.Optim(aCdtSad.mC,1.0,0.01); // Pos Init, Step Init, Step Lim + aCdtSad.mC = aNewP; - if (aValSym< mThresholdSym) + if (aValSym < mThresholdSym) { - aNewP3.push_back(cPt3dr(aNewP.x(),aNewP.y(),aValSym)); + aVCdtSym.push_back(cCdSym(aCdtSad,aValSym)); mDImLabel->SetV(ToI(aNewP),eFilterSym); } } delete aFSym; - aVMaxLoc = aNewP3; } + int aNbSym = aVCdtSym.size(); - if (0) + /* [5] Compute lines, radiom model & correlation */ + std::vector aVCdtRad; + cAutoTimerSegm aTSRadiom(mTimeSegm,"Radiom"); + { + cCubicInterpolator aCubI(-0.5); + cScoreTetaLine aSTL(*mDImIn,aCubI,mLengtSInit,mStepSeg); + for (const auto & aCdtSym : aVCdtSym) + { + cCdRadiom aCdRad = TestBinarization(aSTL,aCdtSym,mThickness); + if (aCdRad.mCostCorrel <= mMaxCostCorrIm) + aVCdtRad.push_back(aCdRad); + } + } + + + +#if (0) + + + if (1) { cFilterDCT * aFSym = cFilterDCT::AllocSym(mImIn,0.0,mDistCalcSym0,1.0); cOptimByStep aOptimSym(*aFSym,true,mDistDivSym); // tPtR Optim(const tPtR & ,tREAL8 aStepInit,tREAL8 aStepLim,tREAL8 aMul=0.5); + cCubicInterpolator aCubI(-0.5); + cScoreTetaLine aSTL(*mDImIn,aCubI,5.0,0.5); cStdStatRes aSSad1; cStdStatRes aSSad0; cStdStatRes aSSymInt_1; cStdStatRes aSSymInt_0; + cStdStatRes aSSCor_0; + cStdStatRes aSSCor_1; + cCalcSaddle aCalcSInit(1.5,1.0); for (const auto & aP3 : aVMaxLoc) { @@ -328,10 +656,9 @@ void cAppliCheckBoardTargetExtract::DoOneImage() cPt2dr aP2 = aP1; aCalcSInit.RefineSadlePointFromIm(aImBlur,aP2,true); - // tREAL8 aSymInt = aFSym->ComputeValMaxCrown(aP2,10); - // tREAL8 aSymInt = OptimFilter(aFSym,aP2).ValExtre(); tREAL8 aSymInt = aOptimSym.Optim(aP0,1.0,0.05).first; + tREAL8 aCorBin = TestBinarization(aSTL,Proj(aP3)); bool Ok = IsPtTest(Proj(aP3)); @@ -341,11 +668,13 @@ void cAppliCheckBoardTargetExtract::DoOneImage() { aSSad1.Add(aP3.z()); aSSymInt_1.Add(aSymInt); + aSSCor_1.Add(aCorBin); } else { aSSad0.Add(aP3.z()); aSSymInt_0.Add(aSymInt); + aSSCor_0.Add(aCorBin); } } @@ -362,17 +691,27 @@ void cAppliCheckBoardTargetExtract::DoOneImage() << " NotOk 50% " << aSSymInt_0.ErrAtProp(0.5) << " 10% " << aSSymInt_0.ErrAtProp(0.10) << "\n"; + + StdOut() << " * CORR , #Max=" << aSSCor_1.Max() << " %75=" << aSSCor_1.ErrAtProp(0.75) + << " NotOk 50% " << aSSCor_0.ErrAtProp(0.5) + << " 10% " << aSSCor_0.ErrAtProp(0.10) + << "\n"; } delete aFSym; } cAutoTimerSegm aTSMakeIm(mTimeSegm,"OTHERS"); - MakeImageSaddlePoints(*mDImIn,*mDImLabel); StdOut() << "NBS " << (100.0*aNbSaddle)/aNbTot << " " << (100.0*aNbCCSad)/aNbTot << " " << (100.0*aVMaxLoc.size())/aNbTot << " NB=" << aVMaxLoc.size() << "\n"; aDImBlur.ToFile("Blurred.tif"); +#endif + cAutoTimerSegm aTSMakeIm(mTimeSegm,"OTHERS"); + MakeImageSaddlePoints(*mDImIn,*mDImLabel); + + StdOut() << "NB Cd, SAD: " << aNbSad0 << " " << aNbSad1 << " " < & aDIm,const cDiffInterpolator1D & anInt,tREAL8 aLength,tREAL8 aStep) : + mDIm (&aDIm), + mLength (aLength), + mNb (round_up(mLength/aStep)), + mStepAbsc (mLength/mNb), + mStepTeta (-1), + mTabInt (anInt,100) +{ +} + +void cScoreTetaLine::SetCenter(const cPt2dr & aC) { mC = aC; } + +/** For an "oriented" segemet S : + * - C the center, Teta the direction, L the lengnt + * We return how the segement is one of the direction of the checkboard pattern. Let + * - T the tangent and N the normal + * - A be the absicsse of a point -L < A < L + * - P = C + A T a point of + * - let G be the gradient of image in P, G' = G/|G| the direction of G + * The score depends of two "signs" s1 and s2 (s1,s2 = "+-1") : + * - s1 depends if we have black ofr white a the left of S (which is oriented) + * - s2 depends if A>0 or A<0 (because the left colour change when cross C) + * So the score is finally : + * + * | G' - s1s2 N| + * + */ + +cPt1dr cScoreTetaLine::Value(const cPt1dr& aPtAngle) const +{ + tREAL8 aTeta = aPtAngle.x(); + cPt2dr aTgt = FromPolar(1.0,aTeta); + cPt2dr aNormal = aTgt * cPt2dr(0,1.0); + cWeightAv aWA; // Average of difference + + for (int aKL=-mNb ; aKL<=mNb ; aKL++) // samples all the point on the line + { + if (aKL) + { + tREAL8 aAbsc = mStepAbsc*aKL; + cPt2dr aPt = mC + aTgt * aAbsc; + auto [aVal,aGrad] = mDIm->GetValueAndGradInterpol(mTabInt,aPt); + + tREAL8 aN2 = Norm2(aGrad); + if (aN2 > 0) + { + // Orientation depend of Sign & position in the segment + aGrad = aGrad / (aN2 * mCurSign * (aKL>0 ? 1 : -1)); + tREAL8 aDist = Norm2(aNormal-aGrad); + aWA.Add(aN2,aDist); + } + } + } + + return cPt1dr(aWA.Average(1e10)); // 1e10 => def value in case no point +} + + +/** Get initial guess of direction, simply parse all possible value at given step */ + +tREAL8 cScoreTetaLine::GetTetasInit(tREAL8 aStepPix,int aCurSign) +{ + mCurSign = aCurSign; + cWhichMin aWMin; // minimal score + + mStepTeta = aStepPix / mLength; // Step in radian to have given step in pixel + int aNbTeta = round_up(M_PI / mStepTeta); + mStepTeta = M_PI / aNbTeta; + + for (int aKTeta=0 ; aKTeta aOpt(*this,true,10.0); // 10.0 => DistMin + cPt1dr aRes = aOpt.Optim(cPt1dr(aTeta0),mStepTeta,aStepPix/mLength).second; + + return aRes.x(); +} + +std::pair cScoreTetaLine::Tetas_CheckBoard(const cPt2dr & aC,tREAL8 aStepInit,tREAL8 aStepLim) +{ + std::vector aVTeta; + SetCenter(aC); + + for (const auto & aSign : {-1,1}) + { + tREAL8 aTetaInit = GetTetasInit(aStepInit,aSign); + tREAL8 aTetaRefine = Refine(aTetaInit,aStepLim,aSign); + + // at this step the angles are define % pi, theoretically in [0,Pi], btw after due to optim it can be slightly outside + aTetaRefine = mod_real(aTetaRefine,M_PI); + aVTeta.push_back(aTetaRefine); + } + + // the angle are defined %pi, to have a single representation we need that teta1/teta2 + // correpond to a direct repair + if (aVTeta.at(0)> aVTeta.at(1)) + aVTeta.at(1) += M_PI; + + return std::pair(aVTeta.at(0),aVTeta.at(1)); +} + +}; diff --git a/MMVII/src/UtiMaths/uti_stat.cpp b/MMVII/src/UtiMaths/uti_stat.cpp index 112b5aafb9..cbfc73545d 100755 --- a/MMVII/src/UtiMaths/uti_stat.cpp +++ b/MMVII/src/UtiMaths/uti_stat.cpp @@ -367,6 +367,32 @@ template cMatIner2Var StatFromImageDist(const cDataIm2D Inv = (-Wx Wx2) + * ----------- + * (Wx2 W1 -WxXx) + * + * + * ( W1 -Wx ) (Wxy) ( W1 Wxy -Wx Wy) + * (-Wx Wx2) (Wy ) = (-Wx Wxy+ Wx2 Wy) + * + */ + +template std::pair cMatIner2Var::FitLineDirect() const +{ + Type aDelta = mS11 * mS0 - Square(mS1); + + Type A = (mS12*mS0 - mS1*mS2) / aDelta; + Type B = (-mS1*mS12 + mS11 * mS2) / aDelta; + + return std::pair(A,B); +} + + /* *********************************************** */ /* */ /* cWeightAv */ @@ -412,6 +438,13 @@ template TypeVal cWeightAv: return mSVW / mSW; } +template TypeVal cWeightAv::Average(const TypeVal & aDef) const +{ + if (! ValidInvertibleFloatValue(mSW)) + return aDef; + return Average(); +} + template const TypeVal & cWeightAv::SVW () const {return mSVW;} template const TypeWeight & cWeightAv::SW () const {return mSW;} From 0197e01ddba9cebd4fd13ff69b83e3ccc2b23b74 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Mon, 10 Jun 2024 19:26:17 +0200 Subject: [PATCH 5/5] New option in codetarget recognition to avoid error in multiple target --- MMVII/Doc/Methods/CodedTarget-Theory.tex | 11 ++ MMVII/include/MMVII_ExtractLines.h | 11 ++ MMVII/include/MMVII_Geom2D.h | 7 + MMVII/src/Appli/cMMVII_Appli_MakeReport.cpp | 37 ++--- MMVII/src/Bench/BenchGeom.cpp | 20 +++ .../CodedTarget/cCheckBoardTargetExtract.cpp | 135 ++++++++++++++---- MMVII/src/CodedTarget/cCircTargetExtract.cpp | 7 +- MMVII/src/Geom2D/GeomsBase2D.cpp | 34 +++++ .../src/ImagesInfoExtract/cScoreTetaLine.cpp | 29 ++-- 9 files changed, 232 insertions(+), 59 deletions(-) diff --git a/MMVII/Doc/Methods/CodedTarget-Theory.tex b/MMVII/Doc/Methods/CodedTarget-Theory.tex index e2833bb47d..6cabd87520 100644 --- a/MMVII/Doc/Methods/CodedTarget-Theory.tex +++ b/MMVII/Doc/Methods/CodedTarget-Theory.tex @@ -1,6 +1,17 @@ \chapter{Coded target} \label{Chap:CodedTarget} +%----------------------------------------------------------------------- +%----------------------------------------------------------------------- +%----------------------------------------------------------------------- + +\section{Todo} + +\begin{itemize} + \item {\tt 2023-12-14-CurTestClino/Calib-1/Processing/} target partially occluded are recognized as others, + for example {\tt 043\_0078.JPG}, target {\tt 66 (?)} is interpreted as {\tt 34} ; +\end{itemize} + %----------------------------------------------------------------------- %----------------------------------------------------------------------- diff --git a/MMVII/include/MMVII_ExtractLines.h b/MMVII/include/MMVII_ExtractLines.h index 95680a4599..309e16591f 100755 --- a/MMVII/include/MMVII_ExtractLines.h +++ b/MMVII/include/MMVII_ExtractLines.h @@ -1,6 +1,9 @@ #ifndef _MMVII_EXTRACT_LINES_H_ #define _MMVII_EXTRACT_LINES_H_ + #include "MMVII_Image2D.h" +#include "MMVII_Interpolators.h" +#include "MMVII_Mappings.h" namespace MMVII @@ -264,6 +267,14 @@ class cScoreTetaLine : public tFunc1DReal // herit from tFunc1DReal for optimiza /// extract the 2 angle of line in checkboar, aStepInit & aStepLim => used in cOptimByStep std::pair Tetas_CheckBoard(const cPt2dr& aC,tREAL8 aStepInit,tREAL8 aStepLim); + + typedef tREAL8 t2Teta[2]; + /// Assure that teta1->teta2 is trigo and teta1 * DIm() const; ///< Accessor + private : /// fix center of reusing the data (to avoid cost for cTabulatedDiffInterpolator) diff --git a/MMVII/include/MMVII_Geom2D.h b/MMVII/include/MMVII_Geom2D.h index efa6831dc1..dfa16d7142 100755 --- a/MMVII/include/MMVII_Geom2D.h +++ b/MMVII/include/MMVII_Geom2D.h @@ -12,6 +12,7 @@ namespace MMVII typedef cSegment tSeg2dr; + /** \file MMVII_Geom2D.h \brief contain classes for geometric manipulation, specific to 2D space : 2D line, 2D plane, rotation, ... @@ -95,10 +96,16 @@ template class cSegment2DCompiled : public cSegmentCompiled Type Dist(const tPt& aPt) const; ///< Faster than upper class const tPt & Normal() const {return mNorm;} + + tPt InterSeg(const cSegment2DCompiled &,tREAL8 aMinAngle=1e-5,bool *IsOk=nullptr); private : tPt mNorm; }; +/** this class a represent a "closed" segment , it has same data than cSegment2DCompiled, + * but as a set/geometric primitive, it is limited by extremities + */ + class cClosedSeg2D { public : diff --git a/MMVII/src/Appli/cMMVII_Appli_MakeReport.cpp b/MMVII/src/Appli/cMMVII_Appli_MakeReport.cpp index 301330b700..a772fa0d45 100755 --- a/MMVII/src/Appli/cMMVII_Appli_MakeReport.cpp +++ b/MMVII/src/Appli/cMMVII_Appli_MakeReport.cpp @@ -139,26 +139,29 @@ void cMMVII_Appli::DoMergeReport() { if (BoolFind(mReport2Merge,anIt.first)) { - cMMVII_Ofs aFileGlob(anIt.second, eFileModeOut::AppendText); - const std::string & anId = anIt.first; - int aNbLines = 0; - if (mRMSWasUsed) - { - for (const auto & aNameIm : VectMainSet(0)) - { - std::string aNameIn = DirSubPReport(anId) + FileOfPath(aNameIm,false) + "." + mMapIdPostReport[anId]; - cMMVII_Ifs aIn(aNameIn, eFileModeIn::Text); - - std::string aLine; - while (std::getline(aIn.Ifs(), aLine)) + // Put aFileGlob in {} to create destruction before OnCloseReport that may generat error + { + cMMVII_Ofs aFileGlob(anIt.second, eFileModeOut::AppendText); + const std::string & anId = anIt.first; + + if (mRMSWasUsed) + { + for (const auto & aNameIm : VectMainSet(0)) { - aFileGlob.Ofs() << aLine<< "\n"; - aNbLines++; - } + std::string aNameIn = DirSubPReport(anId) + FileOfPath(aNameIm,false) + "." + mMapIdPostReport[anId]; + cMMVII_Ifs aIn(aNameIn, eFileModeIn::Text); + + std::string aLine; + while (std::getline(aIn.Ifs(), aLine)) + { + aFileGlob.Ofs() << aLine<< "\n"; + aNbLines++; + } + } } - } - RemoveRecurs(DirSubPReport(anId),false,false); + RemoveRecurs(DirSubPReport(anId),false,false); + } OnCloseReport(aNbLines,anIt.first,anIt.second); } } diff --git a/MMVII/src/Bench/BenchGeom.cpp b/MMVII/src/Bench/BenchGeom.cpp index 51d4a529fa..3e74e006ca 100755 --- a/MMVII/src/Bench/BenchGeom.cpp +++ b/MMVII/src/Bench/BenchGeom.cpp @@ -585,11 +585,31 @@ void BenchMap2D() void BenchPlane3D(); void BenchHomogr2D(); +void BenchSeg2D() +{ + for (int aK=0 ; aK<100 ; aK++) + { + cPt2dr aP1 = cPt2dr::PRandC(); + cPt2dr aT1 = cPt2dr::PRandUnit(); + cSegment2DCompiled aSeg1(aP1,aP1+aT1); + + cPt2dr aP2 = cPt2dr::PRandC(); + cPt2dr aT2 = cPt2dr::PRandUnitNonAligned(aT1); + cSegment2DCompiled aSeg2(aP2,aP2+aT2); + + cPt2dr aI = aSeg1.InterSeg(aSeg2); + + MMVII_INTERNAL_ASSERT_bench(aSeg1.DistLine(aI)<1e-8,"InterSeg"); + MMVII_INTERNAL_ASSERT_bench(aSeg2.DistLine(aI)<1e-8,"InterSeg"); + } +} void BenchGeom(cParamExeBench & aParam) { if (! aParam.NewBench("Geom")) return; + BenchSeg2D(); + BenchSampleQuat(); BenchHomogr2D(); diff --git a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp index f34711fb4b..5713b052bf 100755 --- a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp @@ -8,6 +8,7 @@ #include "MMVII_Sensor.h" #include "MMVII_HeuristikOpt.h" #include "MMVII_ExtractLines.h" +#include "MMVII_TplImage_PtsFromValue.h" namespace MMVII @@ -26,23 +27,23 @@ class cOptimPosSeg : public tFunc2DReal tSeg2dr OptimizeSeg(tREAL8 aStepInit,tREAL8 aStepLim,bool IsMin,tREAL8 aMaxDInfInit) const; private : - cSegment2DCompiled mSeg0; + cSegment2DCompiled mSegInit; cPt2dr mP0Loc; cPt2dr mP1Loc; }; cOptimPosSeg::cOptimPosSeg(const tSeg2dr & aSeg0) : - mSeg0 (aSeg0), - mP0Loc (mSeg0.ToCoordLoc(mSeg0.P1())), - mP1Loc (mSeg0.ToCoordLoc(mSeg0.P2())) + mSegInit (aSeg0), + mP0Loc (mSegInit.ToCoordLoc(mSegInit.P1())), + mP1Loc (mSegInit.ToCoordLoc(mSegInit.P2())) { } cSegment2DCompiled cOptimPosSeg::ModifiedSeg(const cPt2dr & aPModif) const { - cPt2dr aP0Glob = mSeg0.FromCoordLoc(mP0Loc+cPt2dr(0,aPModif.x())); - cPt2dr aP1Glob = mSeg0.FromCoordLoc(mP1Loc+cPt2dr(0,aPModif.y())); + cPt2dr aP0Glob = mSegInit.FromCoordLoc(mP0Loc+cPt2dr(0,aPModif.x())); + cPt2dr aP1Glob = mSegInit.FromCoordLoc(mP1Loc+cPt2dr(0,aPModif.y())); return cSegment2DCompiled(aP0Glob,aP1Glob); } @@ -63,16 +64,31 @@ tSeg2dr cOptimPosSeg::OptimizeSeg(tREAL8 aStepInit,tREAL8 aStepLim,bool IsMin,tR class cOptimSeg_ValueIm : public cOptimPosSeg { public : - cOptimSeg_ValueIm(const tSeg2dr &,tREAL8 aStepOnSeg,const cDataIm2D & aDIm,tREAL8 aTargetValue); + cOptimSeg_ValueIm(const tSeg2dr &,tREAL8 aStepOnSeg,const cDataIm2D & aDIm,tREAL8 aTargetValue); tREAL8 CostOfSeg(const cSegment2DCompiled &) const override; private : tREAL8 mStepOnSeg; int mNbOnSeg; - const cDataIm2D & mDataIm; + const cDataIm2D & mDataIm; tREAL8 mTargetValue; }; +cOptimSeg_ValueIm::cOptimSeg_ValueIm +( + const tSeg2dr & aSegInit, + tREAL8 aStepOnSeg, + const cDataIm2D & aDIm, + tREAL8 aTargetValue +) : + cOptimPosSeg (aSegInit), + mStepOnSeg (aStepOnSeg), + mNbOnSeg (round_up( Norm2(aSegInit.V12()) / mStepOnSeg )) , + mDataIm (aDIm), + mTargetValue (aTargetValue) +{ +} + tREAL8 cOptimSeg_ValueIm::CostOfSeg(const cSegment2DCompiled & aSeg) const { tREAL8 aSum=0; @@ -89,6 +105,7 @@ tREAL8 cOptimSeg_ValueIm::CostOfSeg(const cSegment2DCompiled & aSeg) con + static constexpr tU_INT1 eNone = 0 ; static constexpr tU_INT1 eTopo0 = 1 ; static constexpr tU_INT1 eTopoTmpCC = 2 ; @@ -125,13 +142,16 @@ class cCdRadiom : public cCdSym { public : cCdRadiom(const cCdSym &,const cDataIm2D & aDIm,tREAL8 aTeta1,tREAL8 aTeta2,tREAL8 aThickness); - cMatIner2Var StatGray(const cDataIm2D & aDIm,tREAL8 aThickness,bool IncludeSegs); - tREAL8 mTeta0; - tREAL8 mTeta1; + /// Once blac/white & teta are computed, refine seg using + void OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength); + +// cOptimSeg_ValueIm(const tSeg2dr &,tREAL8 aStepOnSeg,const cDataIm2D & aDIm,tREAL8 aTargetValue); + tREAL8 mTetas[2]; tREAL8 mCostCorrel; // 1-Correlation of model + tREAL8 mRatioBW; // ratio min/max of BW tREAL8 mBlack; tREAL8 mWhite; }; @@ -163,11 +183,12 @@ class cCdRadiomCompiled : public cCdRadiom cSegment2DCompiled mSeg1 ; }; + cCdRadiomCompiled::cCdRadiomCompiled(const cCdRadiom & aCDR,tREAL8 aThickness) : cCdRadiom (aCDR), mThickness (aThickness), - mSeg0 (mC,mC+FromPolar(1.0,mTeta0)), - mSeg1 (mC,mC+FromPolar(1.0,mTeta1)) + mSeg0 (mC,mC+FromPolar(1.0,mTetas[0])), + mSeg1 (mC,mC+FromPolar(1.0,mTetas[1])) { } @@ -222,15 +243,15 @@ std::pair cCdRadiomCompiled::TheorRadiom(const cPt2dr &aPt) con cCdRadiom::cCdRadiom(const cCdSym & aCdSym,const cDataIm2D & aDIm,tREAL8 aTeta0,tREAL8 aTeta1,tREAL8 aThickness) : cCdSym (aCdSym), - mTeta0 (aTeta0), - mTeta1 (aTeta1), - mCostCorrel (2.001) // over maximal theoreticall value + mTetas {aTeta0,aTeta1}, + mCostCorrel (2.001), // over maximal theoreticall value + mRatioBW (0) { static int aCpt=0 ; aCpt++; - cSegment2DCompiled aSeg0 (mC,mC+FromPolar(1.0,mTeta0)); - cSegment2DCompiled aSeg1 (mC,mC+FromPolar(1.0,mTeta1)); + cSegment2DCompiled aSeg0 (mC,mC+FromPolar(1.0,mTetas[0])); + cSegment2DCompiled aSeg1 (mC,mC+FromPolar(1.0,mTetas[1])); static std::vector aDisk = VectOfRadius(0.0,5); cStdStatRes aW0; cStdStatRes aW1; @@ -262,8 +283,8 @@ cCdRadiom::cCdRadiom(const cCdSym & aCdSym,const cDataIm2D & aDIm,tREAL8 } } - tREAL8 aRatioNb = std::min(aNbIn0,aNbIn1) / (tREAL8) std::max(aNbIn0,aNbIn1); - if (aRatioNb <0.2) + mRatioBW = std::min(aNbIn0,aNbIn1) / (tREAL8) std::max(aNbIn0,aNbIn1); + if (mRatioBW <0.05) { return ; } @@ -274,6 +295,27 @@ cCdRadiom::cCdRadiom(const cCdSym & aCdSym,const cDataIm2D & aDIm,tREAL8 mWhite = a+b; } +void cCdRadiom::OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength) +{ + std::vector> aVSegOpt; + for (int aKTeta=0 ; aKTeta<2 ; aKTeta++) + { + cPt2dr aTgt = FromPolar(aLength,mTetas[aKTeta]); + tSeg2dr aSegInit(mC-aTgt,mC+aTgt); + cOptimSeg_ValueIm aOSVI(aSegInit,0.5,aDIm,(mBlack+mWhite)/2.0); + tSeg2dr aSegOpt = aOSVI.OptimizeSeg(0.5,0.01,true,2.0); + + aVSegOpt.push_back(aSegOpt); + mTetas[aKTeta] = Teta(aSegOpt.V12()); + // mTetas[aKTeta] = aSegOpt.I// + } + + cPt2dr aC = aVSegOpt.at(0).InterSeg(aVSegOpt.at(1)); + + mC = aC; + cScoreTetaLine::NormalizeTeta(mTetas); +} + /* *********************************************************** */ /* */ /* cAppliCheckBoardTargetExtract */ @@ -320,6 +362,7 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli // -- tREAL8 mThickness; ///< used for fine estimation of radiom + bool mOptimSegByRadiom; ///< Do we optimize the segment on average radiom // =========== Internal param ============ tIm mImIn; ///< Input global image @@ -343,6 +386,7 @@ cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector= (int) eTopoMaxLoc) { - aRGB.SetRGBPix(aPix,(aDMasq.GetV(aPix)==eFilterSym) ? cRGBImage::Red : cRGBImage::Green ); + cPt3di aCoul = cRGBImage::Yellow; + if (aDMasq.GetV(aPix)== eFilterSym) aCoul = cRGBImage::Green; + if (aDMasq.GetV(aPix)== eFilterRadiom) aCoul = cRGBImage::Red; + aRGB.SetRGBPix(aPix,aCoul); } } aRGB.ToFile("Saddles.tif"); @@ -401,16 +449,25 @@ cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL, cCdRadiom aCdRadiom(aCdSym,*mDImIn,aTeta0,aTeta1,aThickness); + if (mOptimSegByRadiom) + { + aCdRadiom.OptimSegIm(*(aSTL.DIm()),aSTL.Length()); + } + if (IsPtTest(aCdSym.mC)) { + + static int aCpt=0 ; aCpt++; - StdOut() << " CPT=" << aCpt << " Corrrr=" << aCdRadiom.mCostCorrel << " V0="<< aCdRadiom.mBlack << " V1=" << aCdRadiom.mWhite << "\n"; + StdOut() << " CPT=" << aCpt << " Corrrr=" << aCdRadiom.mCostCorrel + << " Ratio=" << aCdRadiom.mRatioBW + << " V0="<< aCdRadiom.mBlack << " V1=" << aCdRadiom.mWhite << "\n"; // StdOut() << " CPT=" << aCpt << " TETASRef= " << aTeta0 << " " << aTeta1 << "\n"; int aZoom = 9; cPt2di aSz(50,50); cPt2di aDec = ToI(aCdSym.mC) - aSz/2; - cPt2dr aCLoc = aCdSym.mC-ToR(aDec); + cPt2dr aCLoc = aCdRadiom.mC-ToR(aDec); cRGBImage aIm = cRGBImage:: FromFile(mNameIm,cBox2di(aDec,aDec+aSz),aZoom); aIm.ResetGray(); @@ -430,22 +487,37 @@ cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL, } int aKT=0; - for (const auto & aTeta : {aTeta0,aTeta1}) + for (const auto & aTeta : aCdRadiom.mTetas) { for (int aK= -aZoom * 20 ; aK<=aZoom*20 ; aK++) { - tREAL8 aAbsc= aK/ (2.0 * aZoom); - cPt2dr aPt = aCLoc + FromPolar(aAbsc,aTeta); - aIm.SetRGBPoint(aPt,cRGBImage::Green); - if (aK==6*aZoom) - aIm.DrawCircle((aKT==0) ? cRGBImage::Blue : cRGBImage::Red,aPt,0.5); + tREAL8 aAbsc= aK/ (2.0 * aZoom); + cPt2dr aPt = aCLoc + FromPolar(aAbsc,aTeta); + if (aK==6*aZoom) + aIm.DrawCircle((aKT==0) ? cRGBImage::Blue : cRGBImage::Red,aPt,0.5); // aIm.SetRGBPoint(aPt, (aKT==0) ? cRGBImage::Blue : cRGBImage::Red); + tREAL8 aSign = ((aAbsc>0) ? 1.0 : -1.0) * ((aKT==0) ? 1 : -1) ; + cPt2dr aNorm = FromPolar(1.0,aTeta + M_PI/2.0) * aSign; + aIm.SetRGBPoint(aPt,cRGBImage::Yellow); + if (std::abs(aAbsc) > 1.0) + { + tREAL8 aV = (aCdRadiom.mBlack+aCdRadiom.mWhite) /2.0; + + cGetPts_ImInterp_FromValue aGIFV(*mDImIn,aV,0.1,aPt+ToR(aDec)-aNorm, aNorm); + if (aGIFV.Ok()) + { + aIm.SetRGBPoint(aGIFV.PRes()-ToR(aDec),cRGBImage::Blue); + } + // StdOut() << "OKKK " << aGIFV.Ok() << " K=" << aK << "\n"; + } + // aIm.SetRGBPoint(aPt,cRGBImage::Green); } aKT++; } aIm.SetRGBPoint(aCLoc,cRGBImage::Red); aIm.ToFile("TestCenter_" + ToStr(aCpt) + ".tif"); +// getchar(); } return aCdRadiom; @@ -462,7 +534,7 @@ void cAppliCheckBoardTargetExtract::DoOneImage() tREAL8 mDistCalcSym0 = 8.0; // distance for evaluating symetry criteria tREAL8 mDistDivSym = 2.0; // maximal distance to initial value in symetry opt - tREAL8 mLengtSInit = 5.0; + tREAL8 mLengtSInit = 05.0; tREAL8 mStepSeg = 0.5; tREAL8 mMaxCostCorrIm = 0.1; @@ -621,7 +693,10 @@ void cAppliCheckBoardTargetExtract::DoOneImage() { cCdRadiom aCdRad = TestBinarization(aSTL,aCdtSym,mThickness); if (aCdRad.mCostCorrel <= mMaxCostCorrIm) + { aVCdtRad.push_back(aCdRad); + mDImLabel->SetV(ToI(aCdRad.mC),eFilterRadiom); + } } } diff --git a/MMVII/src/CodedTarget/cCircTargetExtract.cpp b/MMVII/src/CodedTarget/cCircTargetExtract.cpp index c87ea7334e..a96ef36aa9 100755 --- a/MMVII/src/CodedTarget/cCircTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCircTargetExtract.cpp @@ -671,6 +671,7 @@ class cAppliExtractCircTarget : public cMMVII_Appli, std::string mReportMutipleDetec; // Name for report of multiple detection in on target double mRatioDMML; + int mNbMaxMulTarget; cThresholdCircTarget mThresh; std::vector mGTMissed; @@ -697,7 +698,8 @@ cAppliExtractCircTarget::cAppliExtractCircTarget mPatHihlight ("XXXXX"), mUseSimul (false), mDoReportSimul (false), - mRatioDMML (1.5) + mRatioDMML (1.5), + mNbMaxMulTarget (2) { } @@ -724,6 +726,7 @@ cCollecSpecArg2007 & cAppliExtractCircTarget::ArgOpt(cCollecSpecArg2007 & anArgO << AOpt2007(mZoomVisuSeed,"ZoomVisuSeed","Make a visualisation of seed point",{eTA2007::HDV}) << AOpt2007(mZoomVisuElFinal,"ZoomVisuEllipse","Make a visualisation extracted ellispe & target",{eTA2007::HDV}) << AOpt2007(mPatHihlight,"PatHL","Pattern for highliting targets in visu",{eTA2007::HDV}) + << AOpt2007(mNbMaxMulTarget,"NbMMT","Nb max of multiple target acceptable",{eTA2007::HDV}) << mPhProj.DPPointsMeasures().ArgDirOutOptWithDef("Std") ); } @@ -1060,7 +1063,7 @@ void cAppliExtractCircTarget::OnCloseReport(int aNbLine,const std::string & anId { if (anIdent==mReportMutipleDetec) { - if (aNbLine>4) // Each target is detected 2 times + if (aNbLine>(mNbMaxMulTarget*2)) // Each target is detected 2 times { MMVII_UsersErrror(eTyUEr::eMultipleTargetInOneImage,"Nb Multiple Target = " + ToStr(aNbLine/2)); } diff --git a/MMVII/src/Geom2D/GeomsBase2D.cpp b/MMVII/src/Geom2D/GeomsBase2D.cpp index dc88b927b7..291b26b39a 100755 --- a/MMVII/src/Geom2D/GeomsBase2D.cpp +++ b/MMVII/src/Geom2D/GeomsBase2D.cpp @@ -97,6 +97,40 @@ template Type cSegment2DCompiled::DistClosedSeg(const tPt& aP return std::abs(aPL.y()); } +template + cPtxd cSegment2DCompiled::InterSeg + ( + const cSegment2DCompiled & aSeg2, + tREAL8 aMinAngle, + bool *IsOk + ) +{ + // (aM1 + L1 T1) . N2 = M2 . N2 + // L1 (T1.N2) = N2. (aM2-aM1) + tPt aM1 = this->PMil(); + tPt aM2 = aSeg2.PMil(); + + Type aScalT1N2 = Scal(this->mTgt,aSeg2.mNorm); + + if (std::abs(aScalT1N2) <= aMinAngle) + { + if (IsOk) + { + *IsOk = false; + return cPtxd::Dummy(); + } + MMVII_INTERNAL_ERROR("Segment2DCompiled::InterSeg : almost parallel"); + } + + if (IsOk) + *IsOk = true; + + Type aLambda = Scal(aM2-aM1,aSeg2.mNorm) / aScalT1N2 ; + + return aM1 + this->mTgt * aLambda; +} + + /* ========================== */ /* cClosedSeg2D */ /* ========================== */ diff --git a/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp b/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp index bc138b08a4..bfdace3f18 100755 --- a/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp +++ b/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp @@ -104,25 +104,34 @@ tREAL8 cScoreTetaLine::Refine(tREAL8 aTeta0,tREAL8 aStepPix,int aSign) std::pair cScoreTetaLine::Tetas_CheckBoard(const cPt2dr & aC,tREAL8 aStepInit,tREAL8 aStepLim) { - std::vector aVTeta; + tREAL8 aVTeta[2]; SetCenter(aC); for (const auto & aSign : {-1,1}) { tREAL8 aTetaInit = GetTetasInit(aStepInit,aSign); tREAL8 aTetaRefine = Refine(aTetaInit,aStepLim,aSign); - - // at this step the angles are define % pi, theoretically in [0,Pi], btw after due to optim it can be slightly outside - aTetaRefine = mod_real(aTetaRefine,M_PI); - aVTeta.push_back(aTetaRefine); + aVTeta[(1+aSign)/2] = aTetaRefine; } + NormalizeTeta(aVTeta); + + return std::pair(aVTeta[0],aVTeta[1]); +} - // the angle are defined %pi, to have a single representation we need that teta1/teta2 - // correpond to a direct repair - if (aVTeta.at(0)> aVTeta.at(1)) - aVTeta.at(1) += M_PI; +void cScoreTetaLine::NormalizeTeta(t2Teta & aVTeta) +{ + // the angles are define % pi, t2ddo have a unique representation + for (int aK=0 ; aK<2 ; aK++) + aVTeta[aK] = mod_real(aVTeta[aK],M_PI); + + // we need that teta1/teta2 correpond to a direct repair + if (aVTeta[0]> aVTeta[1]) + aVTeta[1] += M_PI; - return std::pair(aVTeta.at(0),aVTeta.at(1)); } +const tREAL8 & cScoreTetaLine::Length() const {return mLength;} +cDataIm2D * cScoreTetaLine::DIm() const {return mDIm;} + + };