From a2dc1a2c5323e53bc29016ed0d205a5475672e56 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Mon, 1 Jul 2024 17:28:54 +0200 Subject: [PATCH 1/3] In CheckBoard --- MMVII/include/MMVII_ExtractLines.h | 17 +- MMVII/include/MMVII_Geom2D.h | 4 + .../CodedTarget/cCheckBoardTargetExtract.cpp | 96 +++++++-- MMVII/src/Geom2D/GeomsBase2D.cpp | 10 +- .../src/ImagesInfoExtract/cScoreTetaLine.cpp | 185 ++++++++++++++++-- MMVII/src/Perso/cMMVII_ReduceVideo.cpp | 27 ++- 6 files changed, 293 insertions(+), 46 deletions(-) diff --git a/MMVII/include/MMVII_ExtractLines.h b/MMVII/include/MMVII_ExtractLines.h index 309e16591f..2e598d8ca0 100755 --- a/MMVII/include/MMVII_ExtractLines.h +++ b/MMVII/include/MMVII_ExtractLines.h @@ -272,9 +272,15 @@ class cScoreTetaLine : public tFunc1DReal // herit from tFunc1DReal for optimiza /// Assure that teta1->teta2 is trigo and teta1 * DIm() const; ///< Accessor + /// Current lenght can vary + void SetLengthCur(tREAL8 aL); + + tREAL8 Prolongate(tREAL8 aLMax, std::pair & aTeta,bool ReestimateTeta); + + tREAL8 Score2Teta(const std::pair & aTeta, tREAL8 aAbscMin) const; private : /// fix center of reusing the data (to avoid cost for cTabulatedDiffInterpolator) @@ -288,14 +294,21 @@ class cScoreTetaLine : public tFunc1DReal // herit from tFunc1DReal for optimiza /// Value as tFunc1DReal cPt1dr Value(const cPt1dr& aPt) const override; + tREAL8 ScoreOfTeta(const tREAL8 & aTeta, tREAL8 aAbscMin,tREAL8 aSign) const; + cDataIm2D * mDIm; - tREAL8 mLength; + + tREAL8 mLengthInit; + tREAL8 mStepAInit; + tREAL8 mLengthCur; int mNb; tREAL8 mStepAbsc; tREAL8 mStepTeta; cTabulatedDiffInterpolator mTabInt; cPt2dr mC; tREAL8 mCurSign; /// used to specify an orientaion of segment + tREAL8 mStepTetaInit; + tREAL8 mStepTetaLim; }; diff --git a/MMVII/include/MMVII_Geom2D.h b/MMVII/include/MMVII_Geom2D.h index dfa16d7142..bc9f5b5165 100755 --- a/MMVII/include/MMVII_Geom2D.h +++ b/MMVII/include/MMVII_Geom2D.h @@ -58,6 +58,10 @@ template inline T Teta(const cPtxd & aP1) ///< From x,y to To return std::atan2(aP1.y(),aP1.x()); } +/// return the "line" angle : i.e angle between 2 non oriented direction, it's always in [0,PI/2] +template T LineAngles(const cPtxd & aDir1,const cPtxd & aDir2); + + template inline cPtxd ToPolar(const cPtxd & aP1,T aDefTeta) ///< With Def value 4 teta { return IsNotNull(aP1) ? ToPolar(aP1) : cPtxd(0,aDefTeta); diff --git a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp index 5713b052bf..63b7cc7007 100755 --- a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp @@ -152,6 +152,7 @@ class cCdRadiom : public cCdSym tREAL8 mCostCorrel; // 1-Correlation of model tREAL8 mRatioBW; // ratio min/max of BW + tREAL8 mScoreTeta; // ratio min/max of BW tREAL8 mBlack; tREAL8 mWhite; }; @@ -241,7 +242,16 @@ std::pair cCdRadiomCompiled::TheorRadiom(const cPt2dr &aPt) con return std::pair(aPos,aGrayTh); } -cCdRadiom::cCdRadiom(const cCdSym & aCdSym,const cDataIm2D & aDIm,tREAL8 aTeta0,tREAL8 aTeta1,tREAL8 aThickness) : + + +cCdRadiom::cCdRadiom +( + const cCdSym & aCdSym, + const cDataIm2D & aDIm, + tREAL8 aTeta0, + tREAL8 aTeta1, + tREAL8 aThickness +) : cCdSym (aCdSym), mTetas {aTeta0,aTeta1}, mCostCorrel (2.001), // over maximal theoreticall value @@ -249,7 +259,6 @@ cCdRadiom::cCdRadiom(const cCdSym & aCdSym,const cDataIm2D & aDIm,tREAL8 { static int aCpt=0 ; aCpt++; - 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); @@ -364,6 +373,12 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli tREAL8 mThickness; ///< used for fine estimation of radiom bool mOptimSegByRadiom; ///< Do we optimize the segment on average radiom + tREAL8 mLengtSInit;// = 05.0; + tREAL8 mLengtProlong;// = 20.0; + tREAL8 mStepSeg;// = 0.5; + tREAL8 mMaxCostCorrIm;// = 0.1; + int mNumDebug; + // =========== Internal param ============ tIm mImIn; ///< Input global image cPt2di mSzIm; ///< Size of image @@ -387,12 +402,17 @@ cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector 1 + 2 cos(T)sin(T) U.V = 1 + sin(2T) U.V, ValMin -> 1 -U.V + * + */ cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL,const cCdSym & aCdSym,tREAL8 aThickness) { + bool IsMarqed = IsPtTest(aCdSym.mC); + static int aCptGlob=0 ; aCptGlob++; + static int aCptMarq=0 ; if (IsMarqed) aCptMarq++; + DebugCB = (aCptMarq == mNumDebug) && IsMarqed; + + auto aPairTeta = aSTL.Tetas_CheckBoard(aCdSym.mC,0.1,1e-3); + tREAL8 aLength = aSTL.Prolongate(mLengtProlong,aPairTeta,true); + + auto [aTeta0,aTeta1] = aPairTeta; + cPt2dr aV0 = FromPolar(1.0,aTeta0); + cPt2dr aV1 = FromPolar(1.0,aTeta1); - auto [aTeta0,aTeta1] = aSTL.Tetas_CheckBoard(aCdSym.mC,0.5,1e-3); + cAffin2D aMapEll2Ori(aCdSym.mC,aV0*aLength,aV1*aLength); + cAffin2D aMapOri2Ell = aMapEll2Ori.MapInverse(); + + /* x,y = mC + (V0 V1) y + */ cCdRadiom aCdRadiom(aCdSym,*mDImIn,aTeta0,aTeta1,aThickness); if (mOptimSegByRadiom) { - aCdRadiom.OptimSegIm(*(aSTL.DIm()),aSTL.Length()); + aCdRadiom.OptimSegIm(*(aSTL.DIm()),aSTL.LengthCur()); } - if (IsPtTest(aCdSym.mC)) + if (IsMarqed) { - - static int aCpt=0 ; aCpt++; - StdOut() << " CPT=" << aCpt << " Corrrr=" << aCdRadiom.mCostCorrel + StdOut() << " CPT=" << aCptMarq << " " << aCptGlob << " Corrrr=" << aCdRadiom.mCostCorrel << " Ratio=" << aCdRadiom.mRatioBW - << " V0="<< aCdRadiom.mBlack << " V1=" << aCdRadiom.mWhite << "\n"; + << " V0="<< aCdRadiom.mBlack << " V1=" << aCdRadiom.mWhite + << " ScTeta=" << aSTL.Score2Teta(aPairTeta,2.0) + << " LLL=" << aLength + << "\n"; // StdOut() << " CPT=" << aCpt << " TETASRef= " << aTeta0 << " " << aTeta1 << "\n"; int aZoom = 9; @@ -486,6 +530,20 @@ cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL, } } + if (1) + { + cCdRadiomCompiled aCRC(aCdRadiom,2.0); + for (const auto & aPix : aIm.ImR().DIm()) + { + cPt2dr aPtR = ToR(aPix+aDec); + aPtR = aMapOri2Ell.Value(aPtR); + if (Norm2(aPtR) < 1) + { + aIm.SetRGBPix(aPix,cRGBImage::Cyan); + } + } + } + int aKT=0; for (const auto & aTeta : aCdRadiom.mTetas) { @@ -516,7 +574,7 @@ cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL, } aIm.SetRGBPoint(aCLoc,cRGBImage::Red); - aIm.ToFile("TestCenter_" + ToStr(aCpt) + ".tif"); + aIm.ToFile("TestCenter_" + ToStr(aCptMarq) + ".tif"); // getchar(); } @@ -534,10 +592,6 @@ 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 = 05.0; - tREAL8 mStepSeg = 0.5; - tREAL8 mMaxCostCorrIm = 0.1; - // computed threshold tINT8 mDistRectInt = 20; // to see later how we compute it @@ -695,7 +749,7 @@ void cAppliCheckBoardTargetExtract::DoOneImage() if (aCdRad.mCostCorrel <= mMaxCostCorrIm) { aVCdtRad.push_back(aCdRad); - mDImLabel->SetV(ToI(aCdRad.mC),eFilterRadiom); + mDImLabel->SetV(ToI(aCdRad.mC),eFilterRadiom); } } } diff --git a/MMVII/src/Geom2D/GeomsBase2D.cpp b/MMVII/src/Geom2D/GeomsBase2D.cpp index 291b26b39a..b61082e9bb 100755 --- a/MMVII/src/Geom2D/GeomsBase2D.cpp +++ b/MMVII/src/Geom2D/GeomsBase2D.cpp @@ -1273,6 +1273,13 @@ template cTriangle ImageOfTri(const cTriangle< return cTriangle(aMap.Value(aTri.Pt(0)),aMap.Value(aTri.Pt(1)),aMap.Value(aTri.Pt(2))); } +template T LineAngles(const cPtxd & aDir1,const cPtxd & aDir2) +{ + T aRes = std::abs(Teta(aDir1* conj(aDir2))); // Angle Dir1/Dir2 + return std::min(aRes,(T)M_PI - aRes); +} + + /* ========================== */ /* INSTANTIATION */ @@ -1285,7 +1292,8 @@ template std::pair >,std::vector>> Rando template class cSegment2DCompiled;\ template class cAffin2D;\ template class cHomogr2D;\ -template cHomogr2D RandomMapId(TYPE); +template cHomogr2D RandomMapId(TYPE);\ +template TYPE LineAngles(const cPtxd & aDir1,const cPtxd & aDir2); INSTANTIATE_GEOM_REAL(tREAL4) INSTANTIATE_GEOM_REAL(tREAL8) diff --git a/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp b/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp index bfdace3f18..cb7625f664 100755 --- a/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp +++ b/MMVII/src/ImagesInfoExtract/cScoreTetaLine.cpp @@ -15,12 +15,19 @@ namespace MMVII cScoreTetaLine::cScoreTetaLine(cDataIm2D & aDIm,const cDiffInterpolator1D & anInt,tREAL8 aLength,tREAL8 aStep) : mDIm (&aDIm), - mLength (aLength), - mNb (round_up(mLength/aStep)), - mStepAbsc (mLength/mNb), + mLengthInit (aLength), + mStepAInit (aStep), mStepTeta (-1), mTabInt (anInt,100) { + SetLengthCur(mLengthInit); +} + +void cScoreTetaLine::SetLengthCur(tREAL8 aL) +{ + mLengthCur = aL; + mNb = round_up(mLengthCur/mStepAInit); + mStepAbsc = mLengthCur / mNb; } void cScoreTetaLine::SetCenter(const cPt2dr & aC) { mC = aC; } @@ -41,9 +48,8 @@ void cScoreTetaLine::SetCenter(const cPt2dr & aC) { mC = aC; } * */ -cPt1dr cScoreTetaLine::Value(const cPt1dr& aPtAngle) const +tREAL8 cScoreTetaLine::ScoreOfTeta(const tREAL8 & aTeta, tREAL8 aAbscMin,tREAL8 aSign) const { - tREAL8 aTeta = aPtAngle.x(); cPt2dr aTgt = FromPolar(1.0,aTeta); cPt2dr aNormal = aTgt * cPt2dr(0,1.0); cWeightAv aWA; // Average of difference @@ -53,24 +59,41 @@ cPt1dr cScoreTetaLine::Value(const cPt1dr& aPtAngle) const if (aKL) { tREAL8 aAbsc = mStepAbsc*aKL; - cPt2dr aPt = mC + aTgt * aAbsc; - auto [aVal,aGrad] = mDIm->GetValueAndGradInterpol(mTabInt,aPt); + if ( std::abs(aAbsc)>aAbscMin ) + { + cPt2dr aPt = mC + aTgt * aAbsc; + auto [aVal,aGrad] = mDIm->GetValueAndGradInterpol(mTabInt,aPt); - tREAL8 aN2 = Norm2(aGrad); - if (aN2 > 0) - { + 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); - } + aGrad = aGrad / (aN2 * aSign * (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 + return aWA.Average(1e10); // 1e10 => def value in case no point +} + +cPt1dr cScoreTetaLine::Value(const cPt1dr& aPtAngle) const +{ + return cPt1dr(ScoreOfTeta(aPtAngle.x(),1.0,mCurSign)); +} + +tREAL8 cScoreTetaLine::Score2Teta(const std::pair & aTeta, tREAL8 aAbscMin) const +{ + tREAL8 aScTeta0 = ScoreOfTeta(aTeta.first ,aAbscMin,-1); + tREAL8 aScTeta1 = ScoreOfTeta(aTeta.second,aAbscMin,1); + + return std::max(aScTeta0,aScTeta1); } + /** Get initial guess of direction, simply parse all possible value at given step */ tREAL8 cScoreTetaLine::GetTetasInit(tREAL8 aStepPix,int aCurSign) @@ -78,7 +101,7 @@ 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 + mStepTeta = aStepPix / mLengthCur; // Step in radian to have given step in pixel int aNbTeta = round_up(M_PI / mStepTeta); mStepTeta = M_PI / aNbTeta; @@ -87,23 +110,30 @@ tREAL8 cScoreTetaLine::GetTetasInit(tREAL8 aStepPix,int aCurSign) tREAL8 aTeta = aKTeta * mStepTeta; tREAL8 aVal = Value(cPt1dr(aTeta)).x(); aWMin.Add(aTeta,aVal); // update + // } return aWMin.IndexExtre() ; } + + /** Refine the existing value of mimimal using "cOptimByStep", suppose we are close enouh */ tREAL8 cScoreTetaLine::Refine(tREAL8 aTeta0,tREAL8 aStepPix,int aSign) { mCurSign = aSign; cOptimByStep<1> aOpt(*this,true,10.0); // 10.0 => DistMin - cPt1dr aRes = aOpt.Optim(cPt1dr(aTeta0),mStepTeta,aStepPix/mLength).second; + auto [aVal,aRes] = aOpt.Optim(cPt1dr(aTeta0),mStepTeta,aStepPix/mLengthCur); return aRes.x(); } std::pair cScoreTetaLine::Tetas_CheckBoard(const cPt2dr & aC,tREAL8 aStepInit,tREAL8 aStepLim) { + mStepTetaInit = aStepInit; + mStepTetaLim = aStepLim; + SetLengthCur(mLengthInit); + tREAL8 aVTeta[2]; SetCenter(aC); @@ -118,6 +148,125 @@ std::pair cScoreTetaLine::Tetas_CheckBoard(const cPt2dr & aC,tREA return std::pair(aVTeta[0],aVTeta[1]); } +extern bool DebugCB; + +tREAL8 cScoreTetaLine::Prolongate(tREAL8 aLMax,std::pair & aPairTeta,bool ReestimateTeta) +{ + std::vector aVTeta{aPairTeta.first,aPairTeta.second}; + std::vector aVTgt; + std::vector aVNormBlack; + std::vector aVCurOrd; + + std::vector aVBl; + std::vector aVWh; + + // [1] Compute, for the existing "small" segment , the value of bl & white + + for (size_t aKTeta=0 ; aKTeta<2 ; aKTeta++) + { + for (const auto & aSign : {-1,1}) + { + cPt2dr aTgt = FromPolar(1.0,aVTeta.at(aKTeta))*double(aSign); + aVTgt.push_back(aTgt); + cPt2dr aNormBlack = aTgt * cPt2dr(0, (aKTeta==0)?1:-1 ); + aVNormBlack.push_back(aNormBlack); + aVCurOrd.push_back(0.0); + + int aNb = round_up(mLengthInit) ; // (2/pixel) + lenght/2 + for (int aKA=0 ; aKA<= aNb ; aKA++) + { + tREAL8 aAbsc= mLengthInit * (0.5 + aKA/(2.0*aNb)) ; + cPt2dr aPt = mC + aTgt *aAbsc; + + tREAL8 aBl = mDIm->GetVBL(aPt + aNormBlack*(1.0)); + tREAL8 aWh = mDIm->GetVBL(aPt - aNormBlack*(1.0)); + aVBl.push_back(aBl); + aVWh.push_back(aWh); + } + } + } + tREAL8 aBl = NonConstMediane(aVBl); + tREAL8 aWh = NonConstMediane(aVWh); + + tREAL8 aWBound = 0.25; // 0.5 => bound=avg, 0 = initial bound + tREAL8 aMaxBl = (1-aWBound) * aBl + aWBound * aWh ; + tREAL8 aMinWh = aWBound* aBl + (1-aWBound) * aWh ; + + bool aGoOn = true; + tREAL8 aMargin = 2.0; + tREAL8 aAbsc=mLengthCur ; + + tREAL8 aAngleFaisc = LineAngles(FromPolar(1.0,aVTeta.at(0)), FromPolar(1.0,aVTeta.at(1))); + aAngleFaisc = std::min(0.1,aAngleFaisc/4.0); + + + cWeightAv aAvgD; + while (aGoOn) + { + for (int aK=0 ; aK<4 ; aK++) + { + cPt2dr aTgt = aVTgt[aK]; + cPt2dr aNormBlack = aVNormBlack[aK]; + cPt2dr aPt = mC + aTgt *aAbsc + aNormBlack * aVCurOrd[aK]; + + cWhichMax aWMax; + for (tREAL8 aKOrd =-1 ; aKOrd<=1 ; aKOrd++) + { + tREAL8 aDeltaOrd = aKOrd * 0.2; + cPt2dr aNewPt = aPt + aNormBlack * aDeltaOrd; + auto [aVal,aGrad] = mDIm->GetValueAndGradInterpol(mTabInt,aNewPt); + aWMax.Add(aDeltaOrd,Norm2(aGrad)); + } + aVCurOrd[aK] += aWMax.IndexExtre(); + aPt = mC + aTgt *aAbsc + aNormBlack * aVCurOrd[aK]; + + // tREAL8 aOrdRad = std::max(1.0, aAngleFaisc * aAbsc); + tREAL8 aOrdRad = 1.0 + aAngleFaisc * aAbsc; + + tREAL8 aBl = mDIm->GetVBL(aPt + aNormBlack*(aOrdRad)); + tREAL8 aWh = mDIm->GetVBL(aPt - aNormBlack*(aOrdRad)); + if ((aBl>aMaxBl) || (aWhGetValueAndGradInterpol(mTabInt,aPt); + tREAL8 aN2 = Norm2(aGrad); + if (aN2==0) + aGoOn = false; + else + { + aGrad = aGrad/aN2; + if (Norm2(aGrad+aNormBlack) > 0.5) + { + if (DebugCB) + { + StdOut() << "GRRRADD " << aGrad << " " << aNormBlack << " N2=" << Norm2(aGrad+aNormBlack) << "\n"; + } + aGoOn = false; + } + + aAvgD.Add(1.0,Norm2(aGrad+aNormBlack)); + } + } + + if (aAbsc > aLMax+aMargin) aGoOn = false; + if (aGoOn) + aAbsc += 0.5; + } + aAbsc = std::max(mLengthCur,aAbsc-aMargin); + + + if (0) StdOut() << " -------------- ABSC= " << aAbsc << " AVGD=" << aAvgD.Average() << "\n"; + + return aAbsc; +} + + void cScoreTetaLine::NormalizeTeta(t2Teta & aVTeta) { // the angles are define % pi, t2ddo have a unique representation @@ -130,7 +279,7 @@ void cScoreTetaLine::NormalizeTeta(t2Teta & aVTeta) } -const tREAL8 & cScoreTetaLine::Length() const {return mLength;} +const tREAL8 & cScoreTetaLine::LengthCur() const {return mLengthCur;} cDataIm2D * cScoreTetaLine::DIm() const {return mDIm;} diff --git a/MMVII/src/Perso/cMMVII_ReduceVideo.cpp b/MMVII/src/Perso/cMMVII_ReduceVideo.cpp index 87fe7bea92..4debd94eeb 100755 --- a/MMVII/src/Perso/cMMVII_ReduceVideo.cpp +++ b/MMVII/src/Perso/cMMVII_ReduceVideo.cpp @@ -128,15 +128,20 @@ class cAppli_ReduceVideo : public cMMVII_Appli int Exe() override; cCollecSpecArg2007 & ArgObl(cCollecSpecArg2007 & anArgObl) override; cCollecSpecArg2007 & ArgOpt(cCollecSpecArg2007 & anArgOpt) override; + + std::vector Samples() const override; + private : std::string mPat; ///< Pattern of input file std::string mDir; ///< Pattern of input file bool mExec; ///< Execute cat and remove file (else just create file) + bool mIsAudio; ///< For converting mp4->mp3, nor reduction cPt2di mSzReduc; std::string mPrefixRed; std::string mPrefixUnRed; std::string mPrefixRedDone; + std::string mPostFRes; ///< Postfix for result int mLevMin; int mLevMax; @@ -159,10 +164,16 @@ cCollecSpecArg2007 & cAppli_ReduceVideo::ArgOpt(cCollecSpecArg2007 & anArgOpt) return anArgOpt << AOpt2007(mExec,"Exec","Execute reduction, else only print",{eTA2007::HDV}) + << AOpt2007(mIsAudio,"IsAudio","If true, done reduce, but convert mp4,wmv... ->mP3",{eTA2007::HDV}) << AOpt2007(mSzReduc,"SzRed","Reduction size, x=-1 if no reduction",{eTA2007::HDV}) ; } +std::vector cAppli_ReduceVideo::Samples() const +{ + return {}; +// mm3d MapCmd S=1 ffmpeg -i "P=(.*)\.mp4" -vf scale=1280:720 "c=Reduced-\$1.mp4" +} cAppli_ReduceVideo::cAppli_ReduceVideo ( @@ -172,6 +183,7 @@ cAppli_ReduceVideo::cAppli_ReduceVideo cMMVII_Appli (aVArgs,aSpec), mDir (DirCur()), mExec (false), + mIsAudio (false), mSzReduc (1280,720), // 854x480 640x360 mPrefixRed ("Reduced-"), mPrefixUnRed (mPrefixRed+ "UNRED-"), @@ -183,6 +195,13 @@ cAppli_ReduceVideo::cAppli_ReduceVideo int cAppli_ReduceVideo::Exe() { + if (mIsAudio) + { + SetIfNotInit(mPrefixRed,std::string("Audio-")); + } + mPostFRes = mIsAudio ? ".mp3" : ".mp4"; + + mAllFiles = RecGetFilesFromDir(mDir,AllocRegex(mPat),mLevMin,mLevMax); std::set aDirsWithWhite;// store folder with ' ' @@ -216,10 +235,10 @@ int cAppli_ReduceVideo::Exe() RenameFiles(aCurDir+aNameInit,aCurDir+aNameCor); aNameInit = aNameCor; } - std::string aNameTmp = mPrefixRed + LastPrefix(aNameInit) + "-TmpRed" + ".mp4"; + std::string aNameTmp = mPrefixRed + LastPrefix(aNameInit) + "-TmpRed" + mPostFRes; cParamCallSys aCom("ffmpeg","-i",aCurDir+aNameInit); - if (mSzReduc.x()>0) // x=-1 => convention for conserving size + if ((mSzReduc.x()>0) && (!mIsAudio)) // x=-1 => convention for conserving size { aCom.AddArgs("-vf","scale=" + ToStr(mSzReduc.x()) + ":" + ToStr(mSzReduc.y())); } @@ -229,13 +248,13 @@ int cAppli_ReduceVideo::Exe() { if ( GlobSysCall(aCom,true) == EXIT_SUCCESS) { - std::string aNameReduc = mPrefixRed + LastPrefix(aNameInit) + ".mp4"; + std::string aNameReduc = mPrefixRed + LastPrefix(aNameInit) + mPostFRes; std::string aNameDone = mPrefixRed + LastPrefix(aNameInit) + "-DONE."+ LastPostfix(aNameInit) ; // if reduction didnt work, maintain names if (SizeFile(aFullN0) Date: Thu, 4 Jul 2024 19:03:49 +0200 Subject: [PATCH 2/3] In reengenering CodedTargetCheckBoardExtract --- MMVII/include/MMVII_ExtractLines.h | 48 ++ MMVII/include/MMVII_util_tpl.h | 3 + .../CodedTarget/cCheckBoardTargetExtract.cpp | 691 ++++++++---------- .../src/ImagesInfoExtract/cScoreTetaLine.cpp | 72 ++ 4 files changed, 439 insertions(+), 375 deletions(-) diff --git a/MMVII/include/MMVII_ExtractLines.h b/MMVII/include/MMVII_ExtractLines.h index 2e598d8ca0..31cc386a88 100755 --- a/MMVII/include/MMVII_ExtractLines.h +++ b/MMVII/include/MMVII_ExtractLines.h @@ -311,6 +311,54 @@ class cScoreTetaLine : public tFunc1DReal // herit from tFunc1DReal for optimiza tREAL8 mStepTetaLim; }; +/** Base class for optimizing of a segment it inherits of tFunc2DReal so that it can be used by " cOptimByStep<2>" + * in OptimizeSeg . The function "CostOfSeg" must be defined. +*/ + +class cOptimPosSeg : public tFunc2DReal +{ + public : + /// constructor takes initial segment + cOptimPosSeg(const tSeg2dr &); + + + tSeg2dr OptimizeSeg(tREAL8 aStepInit,tREAL8 aStepLim,bool IsMin,tREAL8 aMaxDInfInit) const; + private : + cSegment2DCompiled ModifiedSeg(const cPt2dr & aPModif) const; + /// Inteface CostOfSeg to make it a function of tFunc2DReal + cPt1dr Value(const cPt2dr &) const override; + /// Method to define for indicating the cost to optimize + virtual tREAL8 CostOfSeg(const cSegment2DCompiled &) const = 0; + + cSegment2DCompiled mSegInit; ///< initial value of segment + cPt2dr mP0Loc; ///< First point in local coordinates of segment (y=0) + cPt2dr mP1Loc; ///< Second point in local coordinates of segment (y=0) +}; + +/** Class for refine the position of a segment where the objective is that, on a given image "I" , the points have a given value "V" : + * + * ----- + * \ + * / | I(m) - V| + * Cost(Seg) = ----- + * p in Seg + * */ + +template 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; ///< Sampling step on the seg + int mNbOnSeg; ///< Number of points (computed from mStepOnSeg) + const cDataIm2D & mDataIm; ///< Image on which it is computed + tREAL8 mTargetValue; ///< +}; + + + }; #endif // _MMVII_EXTRACT_LINES_H_ diff --git a/MMVII/include/MMVII_util_tpl.h b/MMVII/include/MMVII_util_tpl.h index 724540cb26..56aa125283 100755 --- a/MMVII/include/MMVII_util_tpl.h +++ b/MMVII/include/MMVII_util_tpl.h @@ -324,6 +324,9 @@ template void ResizeUp(std::vector & aV1,size_t aSz,const Typ aV1.resize(aSz,aVal); } +template void ResizeDown(std::vector & aV1,size_t aSz) { aV1.resize(std::min(aV1.size(),aSz)); } + + template void SetAndResize(std::vector & aVec,size_t aSz,const Type &aVal,const Type & aDef) { ResizeUp(aVec,aSz,aDef); diff --git a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp index 63b7cc7007..b1bcd285a9 100755 --- a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp @@ -15,96 +15,6 @@ 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 mSegInit; - cPt2dr mP0Loc; - cPt2dr mP1Loc; -}; - - -cOptimPosSeg::cOptimPosSeg(const tSeg2dr & aSeg0) : - mSegInit (aSeg0), - mP0Loc (mSegInit.ToCoordLoc(mSegInit.P1())), - mP1Loc (mSegInit.ToCoordLoc(mSegInit.P2())) -{ -} - -cSegment2DCompiled cOptimPosSeg::ModifiedSeg(const cPt2dr & aPModif) const -{ - cPt2dr aP0Glob = mSegInit.FromCoordLoc(mP0Loc+cPt2dr(0,aPModif.x())); - cPt2dr aP1Glob = mSegInit.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; -}; - -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; - - 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 ; @@ -114,6 +24,16 @@ static constexpr tU_INT1 eTopoMaxLoc = 4 ; static constexpr tU_INT1 eFilterSym = 5 ; static constexpr tU_INT1 eFilterRadiom = 6 ; + +class cCdSadle; +class cCdSym; +class cCdRadiom; +class cTmpCdRadiomPos ; + +/* **************************************************** */ +/* cCdSadle */ +/* **************************************************** */ + /// candidate that are pre-selected on the sadle-point criteria class cCdSadle { @@ -126,6 +46,10 @@ class cCdSadle tREAL8 mSadCrit; }; +/* **************************************************** */ +/* cCdSym */ +/* **************************************************** */ + /// candidate that are pre-selected after symetry criterion class cCdSym : public cCdSadle { @@ -141,14 +65,21 @@ class cCdSym : public cCdSadle 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); + /// Cstr use the 2 direction + Thickness of transition between black & white + cCdRadiom(const cCdSym &,const cDataIm2D & aDIm,tREAL8 aTeta1,tREAL8 aTeta2,tREAL8 aLength,tREAL8 aThickness); - /// Once blac/white & teta are computed, refine seg using + /// Once black/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); + /// compute an ellipse that contain + void ComputePtsOfEllipse(std::vector & aRes) const; + + /// Make a visualisation of geometry + void ShowDetail(int aCptMarq,const cScoreTetaLine & aSTL,const std::string &) const; + tREAL8 mTetas[2]; + tREAL8 mLength; ///< length of the lines + tREAL8 mThickness; ///< thickness of the transition B/W tREAL8 mCostCorrel; // 1-Correlation of model tREAL8 mRatioBW; // ratio min/max of BW @@ -157,6 +88,8 @@ class cCdRadiom : public cCdSym tREAL8 mWhite; }; + + enum class eTPosCB { eUndef, @@ -171,10 +104,10 @@ inline bool IsInside(eTPosCB aState) {return (aState==eTPosCB::eInsideBlack) || inline bool IsOk(eTPosCB aState) {return aState!=eTPosCB::eUndef;} /// Used temporary for compilation of radiom -class cCdRadiomCompiled : public cCdRadiom +class cTmpCdRadiomPos : public cCdRadiom { public : - cCdRadiomCompiled(const cCdRadiom &,tREAL8 aThickness); + cTmpCdRadiomPos(const cCdRadiom &,tREAL8 aThickness); /// Theoreticall radiom of modelize checkboard + bool if was computed std::pair TheorRadiom(const cPt2dr &) const; @@ -184,65 +117,11 @@ class cCdRadiomCompiled : public cCdRadiom cSegment2DCompiled mSeg1 ; }; - -cCdRadiomCompiled::cCdRadiomCompiled(const cCdRadiom & aCDR,tREAL8 aThickness) : - cCdRadiom (aCDR), - mThickness (aThickness), - mSeg0 (mC,mC+FromPolar(1.0,mTetas[0])), - mSeg1 (mC,mC+FromPolar(1.0,mTetas[1])) -{ -} - -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::cCdRadiom ( @@ -250,10 +129,13 @@ cCdRadiom::cCdRadiom const cDataIm2D & aDIm, tREAL8 aTeta0, tREAL8 aTeta1, + tREAL8 aLength, tREAL8 aThickness ) : cCdSym (aCdSym), mTetas {aTeta0,aTeta1}, + mLength (aLength), + mThickness (aThickness), mCostCorrel (2.001), // over maximal theoreticall value mRatioBW (0) { @@ -261,7 +143,6 @@ cCdRadiom::cCdRadiom 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; @@ -270,11 +151,12 @@ cCdRadiom::cCdRadiom cMatIner2Var aCorGrayAll; cMatIner2Var aCorGrayInside; - cCdRadiomCompiled aCRC(*this,aThickness); + cTmpCdRadiomPos aCRC(*this,aThickness); - for (const auto & aDelta : aDisk) + std::vector aVPixEllipse; + ComputePtsOfEllipse(aVPixEllipse); + for (const auto & aPImI : aVPixEllipse) { - cPt2di aPImI = aDelta + ToI(mC); tREAL8 aValIm = aDIm.GetV(aPImI); cPt2dr aPImR = ToR(aPImI); @@ -311,7 +193,7 @@ void cCdRadiom::OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength) { cPt2dr aTgt = FromPolar(aLength,mTetas[aKTeta]); tSeg2dr aSegInit(mC-aTgt,mC+aTgt); - cOptimSeg_ValueIm aOSVI(aSegInit,0.5,aDIm,(mBlack+mWhite)/2.0); + 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); @@ -325,6 +207,140 @@ void cCdRadiom::OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength) cScoreTetaLine::NormalizeTeta(mTetas); } +void cCdRadiom::ComputePtsOfEllipse(std::vector & aRes) const +{ + aRes.clear(); + cPt2dr aV0 = FromPolar(mLength,mTetas[0]); + cPt2dr aV1 = FromPolar(mLength,mTetas[1]); + + // ---- x,y -> mC + x V0 + y V1 ------ + cAffin2D aMapEll2Ori(mC,aV0,aV1); + cAffin2D aMapOri2Ell = aMapEll2Ori.MapInverse(); + + cTplBoxOfPts aBox; + int aNbTeta = 100; + for (int aKTeta=0 ; aKTeta aPairTeta(mTetas[0],mTetas[1]); + + StdOut() << " CPT=" << aCptMarq + << " Corrrr=" << mCostCorrel + << " Ratio=" << mRatioBW + << " V0="<< mBlack << " V1=" << mWhite + << " ScTeta=" << aSTL.Score2Teta(aPairTeta,2.0) + << " LLL=" << mLength + << " ThickN=" << mThickness + << "\n"; + + int aZoom = 9; + cPt2di aSz(50,50); + cPt2di aDec = ToI(mC) - aSz/2; + cPt2dr aCLoc = mC-ToR(aDec); + + cRGBImage aIm = cRGBImage:: FromFile(aNameIm,cBox2di(aDec,aDec+aSz),aZoom); + aIm.ResetGray(); + + if (1) // generate the theoretical image + { + cTmpCdRadiomPos aCRC(*this,mThickness); + std::vector aVEllipse; + ComputePtsOfEllipse(aVEllipse); + + for (const auto & aPix : aVEllipse) + { + auto [aState,aGray] = aCRC.TheorRadiom(ToR(aPix)); + if (aState != eTPosCB::eUndef) + { + aIm.SetGrayPix(aPix-aDec,mBlack+ aGray*(mWhite-mBlack)); + } + } + } + + aIm.SetRGBPoint(aCLoc,cRGBImage::Red); + aIm.ToFile("TestCenter_" + ToStr(aCptMarq) + ".tif"); +} + +/* ***************************************************** */ +/* */ +/* cTmpCdRadiomPos */ +/* */ +/* ***************************************************** */ + + +cTmpCdRadiomPos::cTmpCdRadiomPos(const cCdRadiom & aCDR,tREAL8 aThickness) : + cCdRadiom (aCDR), + mThickness (aThickness), + mSeg0 (mC,mC+FromPolar(1.0,mTetas[0])), + mSeg1 (mC,mC+FromPolar(1.0,mTetas[1])) +{ +} + +std::pair cTmpCdRadiomPos::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); +} + + + + /* *********************************************************** */ /* */ /* cAppliCheckBoardTargetExtract */ @@ -353,7 +369,17 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli bool IsPtTest(const cPt2dr & aPt) const; ///< Is it a point marqed a test + /// Method called do each image void DoOneImage() ; + /// Read Image from files, init Label Image, init eventually masq 4 debug, compute blurred version of input image + void ReadImagesAndBlurr(); + /// compute points that are "topologicall" : memorized in label image as label "eTopoTmpCC" + void ComputeTopoSadles(); + /// start from topo point, compute the "saddle" criterion, and use it for filtering on relative max + void SaddleCritFiler() ; + /// start from saddle point, optimize position on symetry criterion then filter on sym thresholds + void SymetryFiler() ; + void MakeImageSaddlePoints(const tDIm &,const cDataIm2D & aDMasq) const; cPhotogrammetricProject mPhProj; @@ -373,26 +399,44 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli tREAL8 mThickness; ///< used for fine estimation of radiom bool mOptimSegByRadiom; ///< Do we optimize the segment on average radiom - tREAL8 mLengtSInit;// = 05.0; - tREAL8 mLengtProlong;// = 20.0; - tREAL8 mStepSeg;// = 0.5; - tREAL8 mMaxCostCorrIm;// = 0.1; - int mNumDebug; + tREAL8 mLengtSInit; ///< = 05.0; + tREAL8 mLengtProlong; ///< = 20.0; + tREAL8 mStepSeg; ///< = 0.5; + tREAL8 mMaxCostCorrIm; ///< = 0.1; + int mNbBlur1; ///< = 4, Number of initial blurring + // + // ---------------- Thresholds for Saddle point criteria -------------------- + tREAL8 mDistMaxLocSad ; ///< =10.0, for supressing sadle-point, not max loc in a neighboorhoud + int mDistRectInt; ///< = 20, insideness of points for seed detection + size_t mMaxNbSP_ML0 ; ///< = 30000 Max number of best point saddle points, before MaxLoc + size_t mMaxNbSP_ML1 ; ///< = 2000 Max number of best point saddle points, after MaxLoc + cPt2di mPtLimCalcSadle; ///< =(2,1) limit point for calc sadle neighbour , included, + + // ---------------- Thresholds for Symetry criteria -------------------- + tREAL8 mThresholdSym ; ///< = 0.5, 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 + int mNumDebug; // =========== Internal param ============ tIm mImIn; ///< Input global image cPt2di mSzIm; ///< Size of image tDIm * mDImIn; ///< Data input image + tIm mImBlur; ///< Blurred image, used in pre-detetction + tDIm * mDImBlur; ///< 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 + std::vector mVCdtSad; ///< Candidate that are selected as local max of saddle criteria + std::vector mNbSads; ///< For info, number of sadle points at different step + std::vector mVCdtSym; ///< Candidate that are selected on the symetry criteria }; /* *************************************************** */ /* */ -/* cAppliCheckBoardTargetExtract */ +/* cAppliCheckBoardTargetExtract */ /* */ /* *************************************************** */ @@ -406,9 +450,20 @@ cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector aMapEll2Ori(aCdSym.mC,aV0*aLength,aV1*aLength); - cAffin2D aMapOri2Ell = aMapEll2Ori.MapInverse(); - /* x,y = mC + (V0 V1) y - */ - - cCdRadiom aCdRadiom(aCdSym,*mDImIn,aTeta0,aTeta1,aThickness); + cCdRadiom aCdRadiom(aCdSym,*mDImIn,aTeta0,aTeta1,aLength,aThickness); if (mOptimSegByRadiom) { @@ -500,48 +547,18 @@ cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL, if (IsMarqed) { - StdOut() << " CPT=" << aCptMarq << " " << aCptGlob << " Corrrr=" << aCdRadiom.mCostCorrel - << " Ratio=" << aCdRadiom.mRatioBW - << " V0="<< aCdRadiom.mBlack << " V1=" << aCdRadiom.mWhite - << " ScTeta=" << aSTL.Score2Teta(aPairTeta,2.0) - << " LLL=" << aLength - << "\n"; - // StdOut() << " CPT=" << aCpt << " TETASRef= " << aTeta0 << " " << aTeta1 << "\n"; - - int aZoom = 9; - cPt2di aSz(50,50); - cPt2di aDec = ToI(aCdSym.mC) - aSz/2; - cPt2dr aCLoc = aCdRadiom.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); - } - } - } + aCdRadiom.ShowDetail(aCptMarq,aSTL,mNameIm); +/* if (1) { - cCdRadiomCompiled aCRC(aCdRadiom,2.0); - for (const auto & aPix : aIm.ImR().DIm()) - { - cPt2dr aPtR = ToR(aPix+aDec); - aPtR = aMapOri2Ell.Value(aPtR); - if (Norm2(aPtR) < 1) - { - aIm.SetRGBPix(aPix,cRGBImage::Cyan); - } - } + std::vector aVPix; + aCdRadiom.ComputePtsOfEllipse(aVPix); + for (const auto & aPix : aVPix) + { + // StdOut() << "D=" << aDec << " " << aPix << "\n"; + aIm.SetRGBPix(aPix -aDec,cRGBImage::Cyan); + } } int aKT=0; @@ -575,36 +592,24 @@ cCdRadiom cAppliCheckBoardTargetExtract::TestBinarization(cScoreTetaLine & aSTL, aIm.SetRGBPoint(aCLoc,cRGBImage::Red); aIm.ToFile("TestCenter_" + ToStr(aCptMarq) + ".tif"); + */ // getchar(); } return aCdRadiom; } -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 - - 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 - - // computed threshold - tINT8 mDistRectInt = 20; // to see later how we compute it - +void cAppliCheckBoardTargetExtract::ReadImagesAndBlurr() +{ /* [0] Initialise : read image and mask */ - cAutoTimerSegm aTSInit(mTimeSegm,"Init"); + cAutoTimerSegm aTSInit(mTimeSegm,"0-Init"); // [0.0] read image mImIn = tIm::FromFile(mNameIm); mDImIn = &mImIn.DIm() ; mSzIm = mDImIn->Sz(); - cRect2 aRectInt = mDImIn->Dilate(-mDistRectInt); // [0.1] initialize labeling image //mImLabel(mSzIm,nullptr,eModeInitImage::eMIA_Null); @@ -613,46 +618,40 @@ void cAppliCheckBoardTargetExtract::DoOneImage() mDImLabel->InitCste(eNone); + // [0.2] Generate potential mask for test points mHasMasqTest = mPhProj.ImageHasMask(mNameIm); if (mHasMasqTest) mMasqTest = mPhProj.MaskOfImage(mNameIm,*mDImIn); - /* [1] Compute a blurred image => less noise, less low level saddle */ - cAutoTimerSegm aTSBlur(mTimeSegm,"Blurr"); - - tIm aImBlur = mImIn.Dup(); // create image blurred with less noise - tDIm& aDImBlur = aImBlur.DIm(); + cAutoTimerSegm aTSBlur(mTimeSegm,"1-Blurr"); - SquareAvgFilter(aDImBlur,mNbBlur1,1,1); + mImBlur = mImIn.Dup(); // create image blurred with less noise + mDImBlur = &(mImBlur.DIm()); + SquareAvgFilter(*mDImBlur,mNbBlur1,1,1); // 1,1 => Nbx,Nby +} - - /* [2] Compute "topological" saddle point */ - - cAutoTimerSegm aTSTopoSad(mTimeSegm,"TopoSad"); +void cAppliCheckBoardTargetExtract::ComputeTopoSadles() +{ + cAutoTimerSegm aTSTopoSad(mTimeSegm,"2.0-TopoSad"); + cRect2 aRectInt = mDImIn->Dilate(-mDistRectInt); // Rectangle excluding point too close to border // 2.1 point with criteria on conexity of point > in neighoor - int aNbSaddle=0; - int aNbTot=0; - for (const auto & aPix : aRectInt) { - if (FlagSup8Neigh(aDImBlur,aPix).NbConComp() >=4) + if (FlagSup8Neigh(*mDImBlur,aPix).NbConComp() >=4) { 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"); + cAutoTimerSegm aTSMaxCC(mTimeSegm,"2.1-MaxCCSad"); int aNbCCSad=0; std::vector aVCC; const std::vector & aV8 = Alloc8Neighbourhood(); @@ -666,76 +665,98 @@ void cAppliCheckBoardTargetExtract::DoOneImage() cWhichMax aBestPInCC; for (const auto & aPixCC : aVCC) { - aBestPInCC.Add(aPixCC,CriterionTopoSadle(aDImBlur,aPixCC)); + aBestPInCC.Add(aPixCC,CriterionTopoSadle(*mDImBlur,aPixCC)); } cPt2di aPCC = aBestPInCC.IndexExtre(); mDImLabel->SetV(aPCC,eTopoMaxOfCC); } } +} - /* [3] Compute point that are max local */ - - - std::vector aVCdtSad; - cAutoTimerSegm aTSCritSad(mTimeSegm,"CritSad"); +/** The saddle criteria is defined by fitting a quadratic function on the image. Having computed the eigen value of quadratic function : + * + * - this criteria is 0 if they have same sign + * - else it is the smallest eigen value + * + * This fitting is done a smpothed version of the image : + * - it seem more "natural" for fitting a smooth model + * - it limits the effect of delocalization + * - it (should be) is not a problem as long as the kernel is smaller than the smallest checkbord we want to detect + * + * As it is used on a purely relative criteria, we dont have to bother how it change the value. + * + */ +void cAppliCheckBoardTargetExtract::SaddleCritFiler() +{ + cAutoTimerSegm aTSCritSad(mTimeSegm,"3.0-CritSad"); - cCalcSaddle aCalcSBlur(aRayCalcSadle+0.001,1.0); + cCalcSaddle aCalcSBlur(Norm2(mPtLimCalcSadle)+0.001,1.0); // structure for computing saddle criteria // [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); - aVCdtSad.push_back(cCdSadle(ToR(aPix),aCritS)); -// cPt3dr(aPix.x(),aPix.y(),aCritS)); + tREAL8 aCritS = aCalcSBlur.CalcSaddleCrit(*mDImBlur,aPix); + mVCdtSad.push_back(cCdSadle(ToR(aPix),aCritS)); } } - int aNbSad0 = aVCdtSad.size(); + mNbSads.push_back(mVCdtSad.size()); // memo size for info - // [3.2] select KBest + MaxLocal - cAutoTimerSegm aTSMaxLoc(mTimeSegm,"MaxLoc"); + // [3.2] sort by decreasing criteria of saddles => "-" aCdt.mSadCrit + limit size + cAutoTimerSegm aTSMaxLoc(mTimeSegm,"3.1-MaxLoc"); - // 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(); + SortOnCriteria(mVCdtSad,[](const auto & aCdt){return - aCdt.mSadCrit;}); + ResizeDown(mVCdtSad,mMaxNbSP_ML0); + mNbSads.push_back(mVCdtSad.size()); // memo size for info - // limit the number of point , a bit rough but first experiment show that sadle criterion is almost perfect on good images - aVCdtSad.resize(std::min(aVCdtSad.size(),size_t(mMaxNbMLS))); - for (const auto & aCdt : aVCdtSad) - mDImLabel->SetV(ToI(aCdt.mC),eTopoMaxLoc); + // [3.3] select MaxLocal + mVCdtSad = FilterMaxLoc((cPt2dr*)nullptr,mVCdtSad,[](const auto & aCdt) {return aCdt.mC;}, mDistMaxLocSad); + mNbSads.push_back(mVCdtSad.size()); // memo size for info - int aNbSad2 = aVCdtSad.size(); - StdOut() << "END MAXLOC \n"; + // [3.3] select KBest + MaxLocal + // limit the number of point , a bit rough but first experiment show that sadle criterion is almost perfect on good images + // mVCdtSad.resize(std::min(mVCdtSad.size(),size_t(mMaxNbMLS))); + ResizeDown(mVCdtSad,mMaxNbSP_ML1); + mNbSads.push_back(mVCdtSad.size()); // memo size for info + for (const auto & aCdt : mVCdtSad) + mDImLabel->SetV(ToI(aCdt.mC),eTopoMaxLoc); +} - /* [4] Calc Symetry criterion */ +void cAppliCheckBoardTargetExtract::SymetryFiler() +{ + cAutoTimerSegm aTSSym(mTimeSegm,"4-SYM"); + cFilterDCT * aFSym = cFilterDCT::AllocSym(mImIn,0.0,mDistCalcSym0,1.0); + cOptimByStep<2> aOptimSym(*aFSym,true,mDistDivSym); - std::vector aVCdtSym; - cAutoTimerSegm aTSSym(mTimeSegm,"SYM"); + for (auto & aCdtSad : mVCdtSad) { - cCalcSaddle aCalcSInit(1.5,1.0); - cFilterDCT * aFSym = cFilterDCT::AllocSym(mImIn,0.0,mDistCalcSym0,1.0); - cOptimByStep aOptimSym(*aFSym,true,mDistDivSym); - - 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) - { - aVCdtSym.push_back(cCdSym(aCdtSad,aValSym)); - mDImLabel->SetV(ToI(aNewP),eFilterSym); - } - } + auto [aValSym,aNewP] = aOptimSym.Optim(aCdtSad.mC,1.0,0.01); // Pos Init, Step Init, Step Lim + aCdtSad.mC = aNewP; - delete aFSym; + if (aValSym < mThresholdSym) + { + mVCdtSym.push_back(cCdSym(aCdtSad,aValSym)); + mDImLabel->SetV(ToI(aNewP),eFilterSym); + } } - int aNbSym = aVCdtSym.size(); + + delete aFSym; +} + +void cAppliCheckBoardTargetExtract::DoOneImage() +{ + /* [0] Initialise : read image , mask + Blurr */ + ReadImagesAndBlurr(); + /* [2] Compute "topological" saddle point */ + ComputeTopoSadles(); + /* [3] Compute point that are max local of saddle point criteria */ + SaddleCritFiler(); + /* [4] Calc Symetry criterion */ + SymetryFiler(); /* [5] Compute lines, radiom model & correlation */ std::vector aVCdtRad; @@ -743,7 +764,7 @@ void cAppliCheckBoardTargetExtract::DoOneImage() { cCubicInterpolator aCubI(-0.5); cScoreTetaLine aSTL(*mDImIn,aCubI,mLengtSInit,mStepSeg); - for (const auto & aCdtSym : aVCdtSym) + for (const auto & aCdtSym : mVCdtSym) { cCdRadiom aCdRad = TestBinarization(aSTL,aCdtSym,mThickness); if (aCdRad.mCostCorrel <= mMaxCostCorrIm) @@ -755,92 +776,12 @@ void cAppliCheckBoardTargetExtract::DoOneImage() } - -#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) - { - cPt2dr aP0 = Proj(aP3); - cPt2dr aP1 = aP0; - aCalcSBlur.RefineSadlePointFromIm(aImBlur,aP1,true); - cPt2dr aP2 = aP1; - aCalcSInit.RefineSadlePointFromIm(aImBlur,aP2,true); - - - tREAL8 aSymInt = aOptimSym.Optim(aP0,1.0,0.05).first; - tREAL8 aCorBin = TestBinarization(aSTL,Proj(aP3)); - - bool Ok = IsPtTest(Proj(aP3)); - - //StdOut() << "SYMM=" << aFSym->ComputeVal(aP0) << "\n"; - - if (Ok) - { - 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); - } - } - - 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"; - - 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"); - - 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 << " " < cOptimPosSeg::ModifiedSeg(const cPt2dr & aPModif) const +{ + cPt2dr aP0Glob = mSegInit.FromCoordLoc(mP0Loc+cPt2dr(0,aPModif.x())); + cPt2dr aP1Glob = mSegInit.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); +} + +/* ************************************************************************ */ +/* */ +/* cOptimSeg_ValueIm */ +/* */ +/* ************************************************************************ */ + +template 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) +{ +} + +template 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; +} + +template class cOptimSeg_ValueIm; + + + /* ************************************************************************ */ /* */ /* cScoreTetaLine */ From bb7383cd332beb1c834339bf89e3d9bb88662629 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Tue, 9 Jul 2024 18:28:47 +0200 Subject: [PATCH 3/3] New CheckB Target : ellipse+ corner --- MMVII/include/MMVII_Geom2D.h | 6 +- .../CodedTarget/cCheckBoardTargetExtract.cpp | 399 +++++++++++++++--- MMVII/src/Geom2D/cEllipse2D.cpp | 35 +- 3 files changed, 367 insertions(+), 73 deletions(-) diff --git a/MMVII/include/MMVII_Geom2D.h b/MMVII/include/MMVII_Geom2D.h index bc9f5b5165..392e3bae17 100755 --- a/MMVII/include/MMVII_Geom2D.h +++ b/MMVII/include/MMVII_Geom2D.h @@ -612,6 +612,7 @@ class cEllipse cPt2dr PtAndGradOfTeta(tREAL8 aTeta,cPt2dr &,tREAL8 aMulRho=1.0) const; /// return also the gradien of belong function cPt2dr ToCoordLoc(const cPt2dr &) const; /// in a sys when ellipse is unity circle + cPt2dr VectToCoordLoc(const cPt2dr &) const; /// for vector (dont use center) cPt2dr FromCoordLoc(const cPt2dr &) const; /// in a sys when ellipse is unity circle cPt2dr VectFromCoordLoc(const cPt2dr &) const; /// for vector (dont use center)in a sys when ellipse is unity circle cPt2dr ToRhoTeta(const cPt2dr &) const; /// Invert function of PtOfTeta @@ -622,6 +623,8 @@ class cEllipse cPt2dr Tgt(const cPt2dr &) const; cPt2dr NormalInt(const cPt2dr &) const; + cPt2dr InterSemiLine(tREAL8 aTeta) const; /// compute the intesection of 1/2 line of direction teta with the ellipse + private : void OneBenchEllispe(); cDenseVect mV; @@ -648,12 +651,13 @@ class cEllipse_Estimate cLeasSqtAA & Sys(); // indicate a rough center, for better numerical accuracy - cEllipse_Estimate(const cPt2dr & aC0); + cEllipse_Estimate(const cPt2dr & aC0,bool isCenterFree=true); void AddPt(cPt2dr aP) ; cEllipse Compute() ; ~cEllipse_Estimate(); private : + bool mIsCenterFree; cLeasSqtAA *mSys; cPt2dr mC0; diff --git a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp index b1bcd285a9..ce1a542e42 100755 --- a/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp +++ b/MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp @@ -60,7 +60,7 @@ class cCdSym : public cCdSadle }; -/// candid obtain after radiometric modelization, +/// candidate obtain after radiometric modelization, class cCdRadiom : public cCdSym { @@ -68,15 +68,31 @@ class cCdRadiom : public cCdSym /// Cstr use the 2 direction + Thickness of transition between black & white cCdRadiom(const cCdSym &,const cDataIm2D & aDIm,tREAL8 aTeta1,tREAL8 aTeta2,tREAL8 aLength,tREAL8 aThickness); + /// Theoretical threshold + tREAL8 Threshold() const ; /// Once black/white & teta are computed, refine seg using void OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength); - /// compute an ellipse that contain + /// compute an ellipse that "safely" contains the points of checkboard + void ComputePtsOfEllipse(std::vector & aRes,tREAL8 aLength) const; + /// call previous with length void ComputePtsOfEllipse(std::vector & aRes) const; + /// compute an ellipse that contain + bool FrontBlackCC(std::vector & aRes,cDataIm2D & aMarq) const; + /// Select point of front that are on ellipse + void SelEllAndRefineFront(std::vector & aRes,const std::vector &) const; + + + /// Is the possibly the point on arc of black ellipse + bool PtIsOnEll(cPt2dr &) const; + /// Is the point on the line for one of angles + bool PtIsOnLine(const cPt2dr &,tREAL8 aTeta) const; + /// Make a visualisation of geometry - void ShowDetail(int aCptMarq,const cScoreTetaLine & aSTL,const std::string &) const; + void ShowDetail(int aCptMarq,const cScoreTetaLine & aSTL,const std::string &,cDataIm2D & aMarq) const; + const cDataIm2D * mDIm; tREAL8 mTetas[2]; tREAL8 mLength; ///< length of the lines tREAL8 mThickness; ///< thickness of the transition B/W @@ -88,6 +104,23 @@ class cCdRadiom : public cCdSym tREAL8 mWhite; }; +class cCdEllipse : public cCdRadiom +{ + public : + cCdEllipse(const cCdRadiom &,cDataIm2D & aMarq); + bool IsOk() const; + const cEllipse & Ell() const; + const cPt2dr & V1() const; + const cPt2dr & V2() const; + + private : + void AssertOk() const; + + bool mIsOk; + cEllipse mEll; + cPt2dr mV1; + cPt2dr mV2; +}; enum class eTPosCB @@ -133,6 +166,7 @@ cCdRadiom::cCdRadiom tREAL8 aThickness ) : cCdSym (aCdSym), + mDIm (&aDIm), mTetas {aTeta0,aTeta1}, mLength (aLength), mThickness (aThickness), @@ -141,11 +175,6 @@ cCdRadiom::cCdRadiom { static int aCpt=0 ; aCpt++; - cSegment2DCompiled aSeg0 (mC,mC+FromPolar(1.0,mTetas[0])); - cSegment2DCompiled aSeg1 (mC,mC+FromPolar(1.0,mTetas[1])); - cStdStatRes aW0; - cStdStatRes aW1; - int aNbIn0=0,aNbIn1=0; cMatIner2Var aCorGrayAll; @@ -186,6 +215,8 @@ cCdRadiom::cCdRadiom mWhite = a+b; } +tREAL8 cCdRadiom::Threshold() const {return (mBlack+mWhite)/2.0;} + void cCdRadiom::OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength) { std::vector> aVSegOpt; @@ -193,7 +224,7 @@ void cCdRadiom::OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength) { cPt2dr aTgt = FromPolar(aLength,mTetas[aKTeta]); tSeg2dr aSegInit(mC-aTgt,mC+aTgt); - cOptimSeg_ValueIm aOSVI(aSegInit,0.5,aDIm,(mBlack+mWhite)/2.0); + cOptimSeg_ValueIm aOSVI(aSegInit,0.5,aDIm,Threshold()); tSeg2dr aSegOpt = aOSVI.OptimizeSeg(0.5,0.01,true,2.0); aVSegOpt.push_back(aSegOpt); @@ -208,10 +239,16 @@ void cCdRadiom::OptimSegIm(const cDataIm2D & aDIm,tREAL8 aLength) } void cCdRadiom::ComputePtsOfEllipse(std::vector & aRes) const +{ + ComputePtsOfEllipse(aRes,mLength); +} + + +void cCdRadiom::ComputePtsOfEllipse(std::vector & aRes,tREAL8 aLength) const { aRes.clear(); - cPt2dr aV0 = FromPolar(mLength,mTetas[0]); - cPt2dr aV1 = FromPolar(mLength,mTetas[1]); + cPt2dr aV0 = FromPolar(aLength,mTetas[0]); + cPt2dr aV1 = FromPolar(aLength,mTetas[1]); // ---- x,y -> mC + x V0 + y V1 ------ cAffin2D aMapEll2Ori(mC,aV0,aV1); @@ -233,7 +270,145 @@ void cCdRadiom::ComputePtsOfEllipse(std::vector & aRes) const } } -void cCdRadiom::ShowDetail(int aCptMarq,const cScoreTetaLine & aSTL,const std::string & aNameIm) const +bool cCdRadiom::PtIsOnLine(const cPt2dr & aPAbs,tREAL8 aTeta) const +{ + cSegment2DCompiled aSeg(mC,mC+FromPolar(1.0,aTeta)); + + cPt2dr aPLoc = aSeg.ToCoordLoc(aPAbs); + + if (std::abs(aPLoc.y()) <= 1.0 + std::abs(aPLoc.x()) /30.0) + return true; + + + + return false; +} + +bool cCdRadiom::PtIsOnEll(cPt2dr & aPtAbs) const +{ + if (Norm2(aPtAbs - mC)<3.0) + return false; + + for (const auto & aTeta : mTetas ) + if (PtIsOnLine(aPtAbs,aTeta)) + return false; + + cGetPts_ImInterp_FromValue aGIFV(*mDIm,Threshold(),0.1,aPtAbs, VUnit(aPtAbs - mC)); + cPt2dr aNewP = aPtAbs; + if (aGIFV.Ok()) + { + aNewP = aGIFV.PRes(); + if (Norm2(aPtAbs-aNewP)>2.0) + return false; + + cPt2dr aPGr = Proj(mDIm->GetGradAndVBL(aNewP)); + tREAL8 aSc = std::abs(Cos(aPGr,aNewP-mC)); + if (aSc<0.5) + return false; + + aPtAbs = aNewP; + } + else + return false; + + // cGetPts_ImInterp_FromValue aGIFV(*mDIm,aV,0.1,aPt+ToR(aDec)-aNorm, aNorm); + /* + cPt2dr aPGr = Proj(mDIm->GetGradAndVBL(aPtAbs)); + if (IsNull(aPGr)) return false; + + cPt2dr aDir = (aPtAbs-mC) / aPGr; + tREAL8 aTeta = + if (Norm2(aPGr) + + */ + + + + return true; +} + +void cCdRadiom::SelEllAndRefineFront(std::vector & aRes,const std::vector & aFrontI) const +{ + aRes.clear(); + for (const auto & aPix : aFrontI) + { + cPt2dr aRPix = ToR(aPix); + if (PtIsOnEll(aRPix)) + aRes.push_back(aRPix); + } +} + +bool cCdRadiom::FrontBlackCC(std::vector & aVFront,cDataIm2D & aDMarq) const +{ + std::vector aRes; + aVFront.clear(); + + std::vector aVPtsEll; + ComputePtsOfEllipse(aVPtsEll,5.0); + + tREAL8 aThrs = Threshold(); + for (const auto & aPix : aVPtsEll) + { + if (mDIm->GetV(aPix) & aV4 = Alloc4Neighbourhood(); + + cRect2 aImOk(aDMarq.Dilate(-10)); + bool isOk = true; + + while (aIndBot != aRes.size()) + { + for (const auto & aDelta : aV4) + { + cPt2di aPix = aRes.at(aIndBot) + aDelta; + if ((aDMarq.GetV(aPix)==0) && (mDIm->GetV(aPix) & aV8 = Alloc8Neighbourhood(); + // compute frontier points + for (const auto & aPix : aRes) + { + bool has8NeighWhite = false; + for (const auto & aDelta : aV8) + { + if (aDMarq.GetV(aPix+aDelta)==0) + has8NeighWhite = true; + } + + if (has8NeighWhite) + aVFront.push_back(aPix); + } + // StdOut() << "FFFF=" << aVFront << "\n"; + + for (const auto & aPix : aRes) + aDMarq.SetV(aPix,0); + + return isOk; +} + + + // if (has8NeighWhite && PtIsOnEll(aRPix)) + +void cCdRadiom::ShowDetail(int aCptMarq,const cScoreTetaLine & aSTL,const std::string & aNameIm,cDataIm2D & aMarq) const { std::pair aPairTeta(mTetas[0],mTetas[1]); @@ -247,14 +422,20 @@ void cCdRadiom::ShowDetail(int aCptMarq,const cScoreTetaLine & aSTL,const std:: << "\n"; int aZoom = 9; - cPt2di aSz(50,50); + cPt2di aSz(100,100); cPt2di aDec = ToI(mC) - aSz/2; cPt2dr aCLoc = mC-ToR(aDec); cRGBImage aIm = cRGBImage:: FromFile(aNameIm,cBox2di(aDec,aDec+aSz),aZoom); aIm.ResetGray(); - if (1) // generate the theoretical image + cPt3di aCol_BlRight = cRGBImage::Red; // color for teta wih blakc on right (on visualization) + cPt3di aCol_BlLeft = cRGBImage::Blue; // color for teta wih blakc on left (on visualization) + cPt3di aCol_StrL = cRGBImage::Yellow; // color for theoreticall straight line + cPt3di aCoulFront = cRGBImage::Cyan; + cPt3di aCoulEllFront = cRGBImage::Orange; + + if (0) // generate the theoretical image + the area (ellipse) of gray modelization { cTmpCdRadiomPos aCRC(*this,mThickness); std::vector aVEllipse; @@ -262,18 +443,144 @@ void cCdRadiom::ShowDetail(int aCptMarq,const cScoreTetaLine & aSTL,const std:: for (const auto & aPix : aVEllipse) { - auto [aState,aGray] = aCRC.TheorRadiom(ToR(aPix)); + auto [aState,aWeightWhite] = aCRC.TheorRadiom(ToR(aPix)); if (aState != eTPosCB::eUndef) { - aIm.SetGrayPix(aPix-aDec,mBlack+ aGray*(mWhite-mBlack)); + tREAL8 aGr = mBlack+ aWeightWhite*(mWhite-mBlack); + cPt3di aCoul (128,aGr,aGr); + // cPt3di aCoul (aGr,aGr,aGr); + aIm.SetRGBPix(aPix-aDec,aCoul); } } } + if (1) // generate the connected + { + std::vector aIFront; + FrontBlackCC(aIFront,aMarq); + if (0) + { + for (const auto & aPix : aIFront) + { + aIm.SetRGBPix(aPix-aDec,aCoulFront); + } + } + + std::vector aEllFr; + SelEllAndRefineFront(aEllFr,aIFront); + for (const auto & aPt : aEllFr) + { + aIm.SetRGBPoint(aPt-ToR(aDec),aCoulEllFront); + } + } + + if (1) + { + cCdEllipse aCDE(*this,aMarq); + int aNb= 500 ; + for (int aK=0 ; aK 1.0) + { + // The orientation of normal change : (1) when we cross intersection + // (2) when we change the segment + tREAL8 aSign = ((aAbsc>0) ? 1.0 : -1.0) * ((aKT==0) ? 1 : -1) ; + cPt2dr aNorm = FromPolar(1.0,aTeta + M_PI/2.0) * aSign; + tREAL8 aV = Threshold(); + + // note that we retract from segment position by "-aNorm" , because the "cGetPts_ImInterp_FromValue" + // is looking only in one direction ! + cGetPts_ImInterp_FromValue aGIFV(*mDIm,aV,0.1,aPt+ToR(aDec)-aNorm, aNorm); + if (aGIFV.Ok()) + { + aIm.SetRGBPoint(aGIFV.PRes()-ToR(aDec),cRGBImage::Green); + } + // StdOut() << "OKKK " << aGIFV.Ok() << " K=" << aK << "\n"; + } + /* + */ + // aIm.SetRGBPoint(aPt,cRGBImage::Green); + } + aKT++; + } + } + + aIm.SetRGBPoint(aCLoc,cRGBImage::Red); aIm.ToFile("TestCenter_" + ToStr(aCptMarq) + ".tif"); } +/* ***************************************************** */ +/* */ +/* cCdEllipse */ +/* */ +/* ***************************************************** */ + +cCdEllipse::cCdEllipse(const cCdRadiom & aCdR,cDataIm2D & aMarq) : + cCdRadiom (aCdR), + mIsOk (false), + mEll (cPt2dr(0,0),0,1,1) +{ + std::vector aIFront; + FrontBlackCC(aIFront,aMarq); + + std::vector aEllFr; + SelEllAndRefineFront(aEllFr,aIFront); + + if (aEllFr.size() < 6) + { + return; + } + + mIsOk = true; + + cEllipse_Estimate anEE(mC,false); + for (const auto & aPixFr : aEllFr) + { + anEE.AddPt(aPixFr); + } + + mEll = anEE.Compute(); + + mV1 = mEll.InterSemiLine(mTetas[0]); + mV2 = mEll.InterSemiLine(mTetas[1]); +} + +bool cCdEllipse::IsOk() const {return mIsOk;} +void cCdEllipse::AssertOk() const +{ + MMVII_INTERNAL_ASSERT_tiny(mIsOk,"No ellipse Ok in cCdEllipse"); +} + +const cEllipse & cCdEllipse::Ell() const {AssertOk(); return mEll;} +const cPt2dr & cCdEllipse::V1() const {AssertOk(); return mV1;} +const cPt2dr & cCdEllipse::V2() const {AssertOk(); return mV2;} + /* ***************************************************** */ /* */ /* cTmpCdRadiomPos */ @@ -428,6 +735,9 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli cIm2D mMasqTest; ///< Possible image of mas 4 debug, print info ... cIm2D mImLabel; ///< Image storing labels of centers cDataIm2D * mDImLabel; ///< Data Image of label + cIm2D mImTmp; ///< Temporary image for connected components + cDataIm2D * mDImTmp; ///< Data Image of "mImTmp" + std::vector mVCdtSad; ///< Candidate that are selected as local max of saddle criteria std::vector mNbSads; ///< For info, number of sadle points at different step std::vector mVCdtSym; ///< Candidate that are selected on the symetry criteria @@ -467,7 +777,9 @@ cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector aVPix; - aCdRadiom.ComputePtsOfEllipse(aVPix); - for (const auto & aPix : aVPix) - { - // StdOut() << "D=" << aDec << " " << aPix << "\n"; - aIm.SetRGBPix(aPix -aDec,cRGBImage::Cyan); - } - } - - int aKT=0; - 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); - 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(aCptMarq) + ".tif"); - */ -// getchar(); + aCdRadiom.ShowDetail(aCptMarq,aSTL,mNameIm,*mDImTmp); } return aCdRadiom; @@ -612,11 +877,13 @@ void cAppliCheckBoardTargetExtract::ReadImagesAndBlurr() mSzIm = mDImIn->Sz(); // [0.1] initialize labeling image - //mImLabel(mSzIm,nullptr,eModeInitImage::eMIA_Null); mDImLabel = &(mImLabel.DIm()); mDImLabel->Resize(mSzIm); mDImLabel->InitCste(eNone); + mDImTmp = &(mImTmp.DIm() ); + mDImTmp->Resize(mSzIm); + mDImTmp->InitCste(0); // [0.2] Generate potential mask for test points mHasMasqTest = mPhProj.ImageHasMask(mNameIm); diff --git a/MMVII/src/Geom2D/cEllipse2D.cpp b/MMVII/src/Geom2D/cEllipse2D.cpp index 9ceba07b35..fdb9260751 100755 --- a/MMVII/src/Geom2D/cEllipse2D.cpp +++ b/MMVII/src/Geom2D/cEllipse2D.cpp @@ -113,11 +113,24 @@ cPt2dr cEllipse::PtOfTeta(tREAL8 aTeta,tREAL8 aMulRho) const return mCenter+ mVGa *(cos(aTeta)*mLGa*aMulRho) + mVSa *(sin(aTeta)*mLSa*aMulRho); } +cPt2dr cEllipse::InterSemiLine(tREAL8 aTeta) const +{ + cPt2dr aPt = VectToCoordLoc(FromPolar(1.0,aTeta)); + aPt = VUnit(aPt); + return FromCoordLoc(aPt); +} cPt2dr cEllipse::ToCoordLoc(const cPt2dr & aP0) const { - cPt2dr aP = (aP0-mCenter)/mVGa; + // cPt2dr aP = (aP0-mCenter)/mVGa; + // return cPt2dr(aP.x()/mLGa,aP.y()/mLSa); + return VectToCoordLoc(aP0-mCenter); +} + +cPt2dr cEllipse::VectToCoordLoc(const cPt2dr & aP0) const +{ + cPt2dr aP = aP0/mVGa; return cPt2dr(aP.x()/mLGa,aP.y()/mLSa); } @@ -437,9 +450,10 @@ void cEllipse::BenchEllispe() /* cEllipseEstimate */ /* */ /* *********************************************************** */ -cEllipse_Estimate::cEllipse_Estimate(const cPt2dr & aC0) : - mSys (new cLeasSqtAA (5)), - mC0 (aC0) +cEllipse_Estimate::cEllipse_Estimate(const cPt2dr & aC0,bool isCenterFree) : + mIsCenterFree (isCenterFree), + mSys (new cLeasSqtAA (5)), + mC0 (aC0) { } @@ -458,8 +472,16 @@ void cEllipse_Estimate::AddPt(cPt2dr aP) aDV(0) = Square(aP.x()); aDV(1) = 2 * aP.x() * aP.y(); aDV(2) = Square(aP.y()); - aDV(3) = aP.x(); - aDV(4) = aP.y(); + if (mIsCenterFree) + { + aDV(3) = aP.x(); + aDV(4) = aP.y(); + } + else + { + aDV(3) = 0.0 ; + aDV(4) = 0.0 ; + } mSys->PublicAddObservation(1.0,aDV,1.0); @@ -470,6 +492,7 @@ cEllipse cEllipse_Estimate::Compute() { auto aSol = mSys->Solve(); + return cEllipse(aSol,mC0); /// return aRes; }