Merge branch 'mpd'

MMVII/include/MMVII_Geom2D.h

@@ -630,6 +630,10 @@ class cEllipse
 
        cPt2dr InterSemiLine(tREAL8 aTeta) const;    /// compute the intesection of 1/2 line of direction teta with the ellipse
 
+       /// get points on ellipse that are +- less regularly sampled at a given step
+       void GetTetasRegularSample(std::vector<tREAL8> & aVTetas,const tREAL8 & aDist);
+
+
     private :
        inline void AssertOk() const;
 

MMVII/src/CodedTarget/FilterCodedTarget.h

@@ -335,6 +335,52 @@ class cCalcSaddle
 std::pair<cIm2D<tREAL4>,cIm2D<tREAL4>> FastComputeSaddleCriterion(cIm2D<tREAL4>  aIm,double aRay);
 
 
+/**  Class for "fine" optimiization of symetry criterion on a image , the input
+ * are : the image "IM", the interpolator "i", a set of point "Pk", an initial center "C",
+ * it trie to adapt "C" so that the equation  "Eq1" is satisfied
+ *
+ *     IM(Pk) = IM(P'k) =  IM(C + C-Pk)  "Eq1"
+ *
+ *    Note that in Pk, we do not store the symetic points P'K, they are computed "on the fly"
+ *
+ *  It offers 2 "strategy" a "fast" one using gradient and another using heuristik
+ *  (which is maintained because it exists and may be safer, but not really sure ...)
+ *
+ *  Inherit from "cDataMapping" to be usable in heuristik optimization by "cOptimByStep".
+ *
+ *  Generally (at least in B/W target) the points will be located on frontier/gradient for 2 reasons :
+ *
+ *      - efficiency in computation (less point)
+ *      - accuraccy, in theoretically homogeneous region, the un-semitry in radiometry is only due to noise
+ */
+
+template <class Type> class cOptimSymetryOnImage : public cDataMapping<tREAL8,2,1>
+{
+        public :
+           typedef  cDataIm2D<Type> tDIm;
+           ///  return symetry  coefficient
+           cPt1dr Value(const cPt2dr & ) const override;
+
+           /// do one iteration of refined with gradient & least square
+           tREAL8 OneIterLeastSqGrad();
+           /// do N iteration of "OneIterLeastSqGrad" until decreasing is small enough or number max of iter reached
+           int IterLeastSqGrad(tREAL8 aGainMin,int aNbMax);
+
+           const cPt2dr & C0() const {return mC0;} ///< Accessor
+
+           /// Constructor memorize parameters and initializd mPtsOpt empty
+           cOptimSymetryOnImage(const cPt2dr & aC0,const tDIm & ,const cDiffInterpolator1D &);
+           /// Add a pts in mPtsOpt
+           void AddPts(const cPt2dr & aPts);
+        protected :
+
+           cPt2dr                       mC0;     ///<   Centers , can be updated
+           const tDIm &                 mDIm;    ///< Image on which optimization os made
+           const cDiffInterpolator1D &  mInterp; ///< Interpolator, used for sub-pixel refinement
+           std::vector<cPt2dr>          mPtsOpt; ///< set of points used for computation
+};
+
+
 };
 #endif // _FILTER_CODED_TARGET_H_
 

MMVII/src/CodedTarget/cCdts_CheckBoardTarget.cpp

@@ -1,8 +1,99 @@
 #include "cCheckBoardTargetExtract.h"
+#include "FilterCodedTarget.h"
 
 namespace MMVII
 {
 
+/* ********************************************* */
+/*                                               */
+/*           cOptimSymetryOnImage                */
+/*                                               */
+/* ********************************************* */
+
+template <class Type> 
+    cOptimSymetryOnImage<Type>::cOptimSymetryOnImage(const cPt2dr & aC0,const tDIm & aDIm,const cDiffInterpolator1D & anInt) :
+         mC0     (aC0),
+         mDIm    (aDIm),
+	 mInterp (anInt)
+{
+}
+
+template <class Type> cPt1dr cOptimSymetryOnImage<Type>::Value(const cPt2dr & aDelta ) const 
+{
+     cSymMeasure<tREAL8> aSymM; // Structure to compute symetry coeff
+     cPt2dr a2NewC = (mC0 + aDelta) * 2.0;  // twice the center actualized
+
+     for (const auto & aP1 : mPtsOpt)
+     {
+          cPt2dr aP2 = a2NewC - aP1;
+	  if (mDIm.InsideInterpolator(mInterp,aP1) && mDIm.InsideInterpolator(mInterp,aP2))
+             aSymM.Add(mDIm.GetValueInterpol(mInterp,aP1),mDIm.GetValueInterpol(mInterp,aP2));
+     }
+
+     return cPt1dr(aSymM.Sym(1e-5));
+}
+
+template <class Type> void cOptimSymetryOnImage<Type>::AddPts(const cPt2dr & aPt)
+{
+   mPtsOpt.push_back(aPt);
+}
+
+template <class Type> tREAL8 cOptimSymetryOnImage<Type>::OneIterLeastSqGrad()
+{
+     cLeasSqtAA<tREAL8>  aSys(2);
+
+     cSymMeasure<tREAL8> aSymM; //
+     for (const auto & aP1 : mPtsOpt)
+     {
+          cPt2dr aP2 = mC0 * 2.0 - aP1;
+	  if (mDIm.InsideInterpolator(mInterp,aP1) && mDIm.InsideInterpolator(mInterp,aP2))
+	  {
+              auto [aV1,aG1] = mDIm.GetValueAndGradInterpol(mInterp,aP1);
+              auto [aV2,aG2] = mDIm.GetValueAndGradInterpol(mInterp,aP2);
+
+	      aSymM.Add(aV1,aV2);
+
+	      cPt2dr a2G2= (aG2*2.0);
+	      cDenseVect<tREAL8> aDV (a2G2.ToVect());
+
+	      //  NewV2 = aV2 + 2 * aG2 . delta  = aV1
+	      aSys.PublicAddObservation(1.0,aDV,aV1-aV2);
+
+
+	      //  StdOut() << aG1 - aG2 << "\n";
+	  }
+     }
+
+     cDenseVect<tREAL8> aSol = aSys.Solve();
+     cPt2dr aDelta = cPt2dr::FromVect(aSol);
+
+     mC0 += aDelta;
+     for (auto & aP1 : mPtsOpt)
+         aP1 += aDelta;
+
+     return  aSymM.Sym(1e-5);
+
+     //  StdOut() << " Ggggg " << aDelta << " " << aSymM.Sym(1e-5) << "\n";
+
+}
+
+template <class Type> int cOptimSymetryOnImage<Type>::IterLeastSqGrad(tREAL8 aGainMin,int aNbMax)
+{
+   tREAL8 aLastScore =  OneIterLeastSqGrad();
+
+   for (int aK= 1 ; aK<aNbMax ; aK++)
+   {
+       tREAL8 aNewScore = OneIterLeastSqGrad();
+       if (aNewScore>aLastScore-aGainMin)
+          return aK+1;
+       aLastScore = aNewScore;
+   }
+   return aNbMax;
+}
+
+template class cOptimSymetryOnImage<tREAL4>;
+
+
 namespace NS_CHKBRD_TARGET_EXTR { 
 
 
@@ -714,57 +805,40 @@ std::pair<eTPosCB,tREAL8>  cTmpCdRadiomPos::TheorRadiom(const cPt2dr &aPt) const
 	return TheorRadiom(aPt,mThickness,0.0);
 }
 
+
 /* ********************************************* */
 /*                                               */
 /*               cOptimPosCdM                    */
 /*                                               */
 /* ********************************************* */
 
-class cOptimPosCdM : public cDataMapping<tREAL8,2,1>
+
+
+class cOptimPosCdM : public cOptimSymetryOnImage<tREAL4>
 {
 	public :
-           cOptimPosCdM(const cCdMerged & aCdM,const cInterpolator1D & );
+           cOptimPosCdM(const cCdMerged & aCdM,const cDiffInterpolator1D & );
 
-           cPt1dr Value(const cPt2dr & ) const override;
+           // cPt1dr Value(const cPt2dr & ) const override;
 	   typedef cSegment2DCompiled<tREAL8> tSeg;
 
 	private :
-	    void AddPts(const cPt2dr & aMaster,  const cPt2dr & aSecond,bool toAvoid2);
-
+	    void AddPts1Seg(const cPt2dr & aMaster,  const cPt2dr & aSecond,bool toAvoid2);
             const cCdMerged&        mCdM;
-	    const cInterpolator1D & mCurInt;
-	    std::vector<cPt2dr>     mPtsOpt;
 };
 
 
 
-
-cPt1dr cOptimPosCdM::Value(const cPt2dr & aDelta ) const 
+cOptimPosCdM::cOptimPosCdM(const cCdMerged & aCdM,const cDiffInterpolator1D & aInt)  :
+        cOptimSymetryOnImage<tREAL4>(aCdM.mC0,(*aCdM.mDIm0),aInt),
+	mCdM      (aCdM)
+	// mCurInt   (aInt)
 {
-     cSymMeasure<tREAL8> aSymM; //
-     cPt2dr a2NewC = (mCdM.mC0 + aDelta) * 2.0;
-
-    const cDataIm2D<tREAL4> & aDIm = *(mCdM.mDIm0);
-     for (const auto & aP1 : mPtsOpt)
-     {
-          cPt2dr aP2 = a2NewC - aP1;
-	  if (aDIm.InsideInterpolator(mCurInt,aP1) && aDIm.InsideInterpolator(mCurInt,aP2))
-             aSymM.Add(aDIm.GetValueInterpol(mCurInt,aP1),aDIm.GetValueInterpol(mCurInt,aP2));
-     }
-
-     return cPt1dr(aSymM.Sym(1e-5));
+	AddPts1Seg(aCdM.CornerlEl_WB(), aCdM.CornerlEl_BW(),true);
+	AddPts1Seg(aCdM.CornerlEl_BW(), aCdM.CornerlEl_WB(),false);
 }
 
-
-cOptimPosCdM::cOptimPosCdM(const cCdMerged & aCdM,const cInterpolator1D & aInt)  :
-	mCdM      (aCdM),
-	mCurInt   (aInt)
-{
-	AddPts(mCdM.CornerlEl_WB(), mCdM.CornerlEl_BW(),true);
-	AddPts(mCdM.CornerlEl_BW(), mCdM.CornerlEl_WB(),false);
-}
-
-void cOptimPosCdM::AddPts(const cPt2dr & aSCorn1, const cPt2dr & aSCorn2,bool toAvoid2)
+void cOptimPosCdM::AddPts1Seg(const cPt2dr & aSCorn1, const cPt2dr & aSCorn2,bool toAvoid2)
 {
      cPt2dr  aCorn1 = aSCorn1  * mCdM.mScale;
      cPt2dr  aCorn2 = aSCorn2  * mCdM.mScale;
@@ -815,13 +889,26 @@ cCdMerged::cCdMerged(const cDataIm2D<tREAL4> * aDIm0,const cCdEllipse & aCDE,tRE
 {
 }
 
-void  cCdMerged::OptimizePosition(const cInterpolator1D & anInt,tREAL8 aStepEnd)
+void cCdMerged::GradOptimizePosition(const cDiffInterpolator1D & anInt,tREAL8 aStepEnd)
+{
+    cOptimPosCdM aCdGrad(*this,anInt);
+         // StdOut() << "-----------  TEST GRAD ----------------  V0=" << aCdGrad.Value(cPt2dr(0,0)) << " \n";
+
+    aCdGrad.IterLeastSqGrad(aStepEnd,5);
+    mC0 = aCdGrad.C0();
+}
+
+void  cCdMerged::HeuristikOptimizePosition(const cDiffInterpolator1D & anInt,tREAL8 aStepEnd)
 {
-     cOptimPosCdM aCdtOpt(*this,anInt);
-     cOptimByStep anOpt(aCdtOpt,true,1.0);
-     auto [aVal,aDelta] =   anOpt.Optim(cPt2dr(0,0),0.02,aStepEnd);
+      cOptimPosCdM aCdtOpt(*this,anInt);
+      cOptimByStep anOpt(aCdtOpt,true,1.0);
+
+
+      // tREAL8  aV0 = aCdtOpt.Value(cPt2dr(0,0)).x();
+      auto [aVal,aDelta] =   anOpt.Optim(cPt2dr(0,0),0.02,aStepEnd);
 
-     mC0 = mC0 + aDelta;
+     //  StdOut() << "  HEURISTIK =" << aV0 << " => " <<  aVal << "\n";
+      mC0 = mC0 + aDelta;
 }
 
 

MMVII/src/CodedTarget/cCheckBoardTargetExtract.cpp

@@ -45,7 +45,7 @@ cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector<s
    mMaxCostCorrIm    (0.1),
    mNbMaxBlackCB     (10000),
    mPropGrayDCD      (2.0/3.0),
-   mNbBlur1          (4),
+   mNbBlur1          (1),
    mStrShow          (""),
    mScales           {1.0},
    mDistMaxLocSad    (10.0),
@@ -54,13 +54,14 @@ cAppliCheckBoardTargetExtract::cAppliCheckBoardTargetExtract(const std::vector<s
    mMaxNbSP_ML1      (2000),
    mPtLimCalcSadle   (2,1),
    mThresholdSym     (0.5),
-   mRayCalcSym0      (8.0),
+   mRayCalcSym0      (4.0),
    mDistDivSym       (2.0),
    mNumDebugMT       (-1),
    mNumDebugSaddle   (-1),
    mNbMinPtEllipse   (6),
    mTryC             (true),
-   mStepRefinePos    (1e-3),
+   mStepHeuristikRefinePos    (-1),
+   mStepGradRefinePos (1e-4),
    mZoomVisuDetec    (9),
    mDefSzVisDetec    (150),
    mSpecif           (nullptr),
@@ -109,7 +110,9 @@ cCollecSpecArg2007 & cAppliCheckBoardTargetExtract::ArgOpt(cCollecSpecArg2007 &
              <<  AOpt2007(mLInitProl,"LSIP","Length Segment Init, for prolongation",{eTA2007::HDV})
              <<  AOpt2007(mNbMinPtEllipse,"NbMinPtEl","Number minimal of point for ellipse estimation",{eTA2007::HDV})
              <<  AOpt2007(mTryC,"TryC","Try also circle when ellipse fails",{eTA2007::HDV})
-             <<  AOpt2007(mStepRefinePos,"StepRefinePos","Step Refine final position with SinC interpol & over sampling (<0 : no refine)",{eTA2007::HDV})
+             <<  AOpt2007(mStepHeuristikRefinePos,"HeuristikStepRefinePos","Step Gradient-Refine final position with SinC interpol & over sampling (<0 : no refine)",{eTA2007::HDV,eTA2007::Tuning})
+
+             <<  AOpt2007(mStepGradRefinePos,"GradStepRefinePos","Step Gradient-Refine final position with SinC interpol & over sampling (<0 : no refine)",{eTA2007::HDV})
 	     <<  AOpt2007(mScales,"Scales","Diff scales of compute (! 0.5 means bigger)",{eTA2007::HDV})
              <<  AOpt2007(mOptimSegByRadiom,"OSBR","Optimize segement by radiometry",{eTA2007::HDV})
              <<  AOpt2007(mNbMaxBlackCB,"NbMaxBlackCB","Number max of point in black part of check-board ",{eTA2007::HDV})
@@ -579,13 +582,22 @@ void cAppliCheckBoardTargetExtract::DoOneImage()
 
     for (const auto & aScale : mScales)
     {
-        DoOneImageAndScale(aScale,mImIn0.Scale(aScale));
+        cAutoTimerSegm aTSRefine(mTimeSegm,"00-Scaling");
+        auto aScIm = mImIn0.Scale(aScale);
+        DoOneImageAndScale(aScale,aScIm);
     }
 
-    if (mStepRefinePos>0)
+    if (mStepHeuristikRefinePos>0)
+    {
+        cAutoTimerSegm aTSRefine(mTimeSegm,"RefineHeuristik");
+        for (auto & aCdtM : mVCdtMerged)
+            aCdtM.HeuristikOptimizePosition(*mInterpol,mStepHeuristikRefinePos);
+    }
+    if (mStepGradRefinePos>0)
     {
+        cAutoTimerSegm aTSRefine(mTimeSegm,"RefineGradient");
         for (auto & aCdtM : mVCdtMerged)
-            aCdtM.OptimizePosition(*mInterpol,mStepRefinePos);
+            aCdtM.GradOptimizePosition(*mInterpol,mStepGradRefinePos);
     }
 
 

MMVII/src/CodedTarget/cCheckBoardTargetExtract.h

@@ -191,7 +191,9 @@ class cCdMerged : public  cCdEllipse
 
 	    tREAL8 mScale;
 	    cPt2dr mC0;  // center at initial image scale
-            void  OptimizePosition(const cInterpolator1D &,tREAL8 aStepEnd);
+            void  HeuristikOptimizePosition(const cDiffInterpolator1D &,tREAL8 aStepEnd);
+
+            void  GradOptimizePosition(const cDiffInterpolator1D &,tREAL8 aStepEnd);
 
 	    const cDataIm2D<tREAL4> * mDIm0;
 };
@@ -345,7 +347,10 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli
         // ---------------- Thresholds for Ellipse  criteria --------------------
         int                   mNbMinPtEllipse;
 	bool                  mTryC;
-	tREAL8                mStepRefinePos;
+
+	tREAL8                mStepHeuristikRefinePos;
+	tREAL8                mStepGradRefinePos;
+	bool                  mDoGradRefine;
 
         // =========== Internal param ============
 
@@ -378,7 +383,7 @@ class cAppliCheckBoardTargetExtract : public cMMVII_Appli
 	std::vector<cCdMerged> mVCdtMerged; // Candidate merged form various scales
 	tREAL8                mCurScale;    /// Memorize the current value of the scale
 	bool                  mMainScale;   /// Is it the first/main scale 
-	cInterpolator1D *     mInterpol;
+	cDiffInterpolator1D * mInterpol;
 };