diff --git a/MMVII/include/MMVII_DeclareAllCmd.h b/MMVII/include/MMVII_DeclareAllCmd.h index b8b98c1c31..9f06889263 100755 --- a/MMVII/include/MMVII_DeclareAllCmd.h +++ b/MMVII/include/MMVII_DeclareAllCmd.h @@ -96,6 +96,7 @@ extern cSpecMMVII_Appli TheSpecTestSensor; extern cSpecMMVII_Appli TheSpecParametrizeSensor; extern cSpecMMVII_Appli TheSpec_TutoSerial; extern cSpecMMVII_Appli TheSpec_TutoFormalDeriv; +extern cSpecMMVII_Appli TheSpecAppliExtractLine; }; #endif // _MMVII_DeclareAllCmd_H_ diff --git a/MMVII/include/MMVII_ExtractLines.h b/MMVII/include/MMVII_ExtractLines.h new file mode 100755 index 0000000000..9d7ea5c525 --- /dev/null +++ b/MMVII/include/MMVII_ExtractLines.h @@ -0,0 +1,195 @@ +#ifndef _MMVII_EXTRACT_LINES_H_ +#define _MMVII_EXTRACT_LINES_H_ +#include "MMVII_Image2D.h" + + +namespace MMVII +{ + +class cHoughTransform; +template class cImGradWithN; +template class cExtractLines; + + +class cHoughPS : public cMemCheck +{ + public : + typedef cSegment2DCompiled tSeg; + + cHoughPS(const cHoughTransform *,const cPt2dr & aPS,tREAL8 aCumul,const cPt2dr & aP1,const cPt2dr & aP2); + /// Angular distance for line anti parallel + tREAL8 DistAnglAntiPar(const cHoughPS& aPS2) const; + tREAL8 DY(const cHoughPS&) const; + tREAL8 Dist(const cHoughPS&,const tREAL8 &aFactTeta=1.0) const; + + const cPt2dr & TetaRho() const; ///< Accessor + const tREAL8 & Teta() const; ///< Accessor + const tREAL8 & Rho() const; ///< Accessor + const tSeg & Seg() const ; ///< Accessor + cHoughPS * Matched() const; ///< Accessor + const tREAL8 & Cumul() const; ///< Accessor + + cPt2dr IndTetaRho() const; ///< Teta/Rho in hough accum dynamic + + void Test(const cHoughPS & ) const; + + bool Match(const cHoughPS &,bool IsDark,tREAL8 aMaxTeta,tREAL8 aDMin,tREAL8 aDMax) const; + + static void SetMatch(std::vector& mVPS,bool IsLight,tREAL8 aMaxTeta,tREAL8 aDMin,tREAL8 aDMax); + + + private : + void InitMatch(); + void UpdateMatch(cHoughPS *,tREAL8 aDist); + + const cHoughTransform * mHT; + cPt2dr mTetaRho; + tREAL8 mCumul; + tSeg mSegE; + cHoughPS * mMatched; + tREAL8 mDistM; +}; + + +/** cHoughTransform + + for a line L (rho=R,teta=T) with , vector orthog to (cos T,sin T) : + the equation is : + + (x,y) in L <=> x cos(T) + y Sin(T) = R + + For now we consider oriented line, typically when the point on a line comes from gradient, + in this case we have : + + T in [0,2Pi] R in [-RMax,+RMax] +*/ +class cHoughTransform +{ + public : + cHoughTransform + ( + const cPt2dr & aSzIn, // Sz of the space + const cPt2dr & aMulTetaRho, // Multiplicator of Rho & Teta + const tREAL8 & aSigmTeta, // Incertitude on gradient direction + cPerspCamIntrCalib * aCalib = nullptr // possible internal calib for distorsion correction + ); + + /// Add a point with a given direction + void AccumulatePtAndDir(const cPt2dr & aPt,tREAL8 aTeta0,tREAL8 aWeight); + cIm2D Accum() const; ///< Accessor + + tREAL8 AvgS2() const; ///< Avg of square, possibly use to measure "compactness" + // tREAL8 Max() const; ///< Max of val, possibly use to measure "compactness" + + + /// Extract the local maxima retur point/line in hough space + value of max + std::vector ExtractLocalMax + ( + size_t aNbMax, // bound to number of max-loc + tREAL8 aDist, // min distance between maxima + tREAL8 aThrAvg, // threshold on Average, select if Accum > aThrAvg * Avg + tREAL8 aThrMax // threshold on Max, select if Accum > Max * Avg + ) const; + + /// max the conversion houg-point -> euclidian line + rho teta + cHoughPS * PtToLine(const cPt3dr &) const; + + const tREAL8 & RhoMax() const; ///< Accessor + /// return the angle teta, of a given index/position in hough accumulator + inline tREAL8 RInd2Teta(tREAL8 aIndTeta) const {return aIndTeta *mFactI2T;} + /// idem, return teta for "integer" index + inline tREAL8 Ind2Teta(int aK) const {return RInd2Teta(aK);} + /// for a given teta, return the index (RInd2Teta o Teta2RInd = Identity) + inline tREAL8 Teta2RInd(const tREAL8 & aTeta) const {return aTeta /mFactI2T;} + + /// return the rho of a given index/position in hough accumulator + inline tREAL8 RInd2Rho(const tREAL8 & aRInd) const { return (aRInd-1.0) / mMulRho - mRhoMax; } + /// return the index of a given rho (Rho2RInd o RInd2Rho = Identity) + inline tREAL8 Rho2RInd(const tREAL8 & aRho) const {return 1.0+ (aRho+mRhoMax) * mMulRho;} + + tREAL8 GetValueBlob(cPt2di aP,int aMaxNeigh) const; + private : + + // inline tREAL8 R2Teta(tREAL8 aIndTeta) const {return aIndTeta *mFactI2T;} + + cPt2dr mMiddle; ///< Middle point, use a origin of Rho + tREAL8 mRhoMax; ///< Max of distance to middle point + tREAL8 mMulTeta; ///< Teta multiplier, if =1 , 1 pix teta ~ 1 pix init (in worst case) + tREAL8 mMulRho; ///< Rho multiplier , if =1, 1 pix-rho ~ 1 pix init + tREAL8 mSigmTeta; ///< incertitude on teta + cPerspCamIntrCalib* mCalib; ///< Potential calibration for distorsion + int mNbTeta; ///< Number of Teta for hough-accum + tREAL8 mFactI2T ; ///< Ratio Teta-Radian / Teta-Index + int mNbRho; ///< Number of Rho for hough-accum + + cIm1D mTabSin; ///< Tabulation of sinus for a given index of teta + cDataIm1D& mDTabSin; ///< Data Image of "mTabSin" + cIm1D mTabCos; ///< Tabulation of co-sinus for a given index of teta + cDataIm1D& mDTabCos; ///< Data Image of "mTabCos" + cIm2D mAccum; ///< Accumulator of Hough + cDataIm2D& mDAccum; ///< Data Image of "mAccum" +}; + +/** Class for storing Grad + its norm + */ +template class cImGradWithN : public cImGrad +{ + public : + /// Constructor using size + cImGradWithN(const cPt2di & aSz); + /// Constructor using image & Alpha-Deriche parameters + cImGradWithN(const cDataIm2D & aImIn,Type aAlphaDeriche); + + + + /// Is it a local-maxima in the direction of the gradient + bool IsMaxLocDirGrad(const cPt2di& aPix,const std::vector &,tREAL8 aRatioXY = 1.0) const; + /// Allocat the neighbourhood use for computing local-maxima + static std::vector NeighborsForMaxLoc(tREAL8 aRay); + cIm2D NormG() {return mNormG;} ///< Accessor + cPt2dr RefinePos(const cPt2dr &) const; ///< Refine a sub-pixelar position of the contour + private : + cPt2dr OneRefinePos(const cPt2dr &) const; ///< One iteration of refinement + + cIm2D mNormG; ///< Image of norm of gradient + cDataIm2D& mDataNG; ///< Data Image of "mNormG" +}; +/// Compute the deriche + its norm +template void ComputeDericheAndNorm(cImGradWithN & aResGrad,const cDataIm2D & aImIn,double aAlpha) ; + + +/** Class for extracting line using gradient & hough transform*/ + +template class cExtractLines +{ + public : + typedef cIm2D tIm; + typedef cDataIm2D tDIm; + + cExtractLines(tIm anIm); ///< constructor , memorize image + ~cExtractLines(); + + /// initialize the gradient + void SetDericheGradAndMasq(tREAL8 aAlphaDerich,tREAL8 aRayMaxLoc,int aBorder,bool Show=false); + /// Initialize the hough transform + void SetHough(const cPt2dr & aMulTetaRho,tREAL8 aSigmTeta,cPerspCamIntrCalib *,bool AffineMax,bool Show=false); + + /// Generate an image for visualizing the contour, + cRGBImage MakeImageMaxLoc(tREAL8 aAlphaTransparency); + + cHoughTransform & Hough(); ///< Accessor + cImGradWithN & Grad(); ///< Acessor + private : + cPt2di mSz; ///< Size of the image + tIm mIm; ///< Memorize the image + cIm2D mImMasqCont; ///< Masq of point selected as contour + int mNbPtsCont; ///< Number of point detected as contour + + cImGradWithN * mGrad; ///< Structure allocated for computing gradient + cHoughTransform * mHough; ///< Structure allocatedf or computing hough + cPerspCamIntrCalib * mCalib; ///< (Optional) calibration for distorsion correction +}; + +}; +#endif // _MMVII_EXTRACT_LINES_H_ + // diff --git a/MMVII/include/MMVII_Geom2D.h b/MMVII/include/MMVII_Geom2D.h index db70974724..2eb1f23271 100755 --- a/MMVII/include/MMVII_Geom2D.h +++ b/MMVII/include/MMVII_Geom2D.h @@ -49,6 +49,12 @@ template inline cPtxd ToPolar(const cPtxd & aP1) ///< Fro AssertNonNul(aP1); return cPtxd(std::hypot(aP1.x(),aP1.y()),std::atan2(aP1.y(),aP1.x())); } +template inline T Teta(const cPtxd & aP1) ///< From x,y to To rho,teta +{ + AssertNonNul(aP1); + return std::atan2(aP1.y(),aP1.x()); +} + template inline cPtxd ToPolar(const cPtxd & aP1,T aDefTeta) ///< With Def value 4 teta { return IsNotNull(aP1) ? ToPolar(aP1) : cPtxd(0,aDefTeta); @@ -70,6 +76,9 @@ template inline cPtxd PSymXY (const cPtxd & aP) /// matrix of linear function q -> q * aP template cDenseMatrix MatOfMul (const cPtxd & aP); +/** This specialization is specific to dim 2, as the normal to a vector is + * specific to d2 + */ template class cSegment2DCompiled : public cSegmentCompiled { public : diff --git a/MMVII/include/MMVII_Linear2DFiltering.h b/MMVII/include/MMVII_Linear2DFiltering.h index 2e1778b52e..122356b0e3 100755 --- a/MMVII/include/MMVII_Linear2DFiltering.h +++ b/MMVII/include/MMVII_Linear2DFiltering.h @@ -362,8 +362,18 @@ template class cGaussianPyramid : public cMemCheck template class cImGrad { public : - cIm2D mGx; - cIm2D mGy; + cIm2D mGx; + cDataIm2D* mDGx; + cIm2D mGy; + cDataIm2D* mDGy; + + const Type & Gx(const cPt2di & aPix) const {return mDGx->GetV(aPix);} + const Type & Gy(const cPt2di & aPix) const {return mDGy->GetV(aPix);} + cPtxd Grad(const cPt2di & aPix) const {return cPtxd(Gx(aPix),Gy(aPix));} + + Type GxBL(const cPt2dr & aPix) const {return mDGx->GetVBL(aPix);} + Type GyBL(const cPt2dr & aPix) const {return mDGy->GetVBL(aPix);} + cPtxd GradBL(const cPt2dr & aPix) const {return cPtxd(GxBL(aPix),GyBL(aPix));} /// memorize cImGrad(const cIm2D & aGx,const cIm2D & aGy) ; /// Use size to allocate @@ -373,6 +383,8 @@ template class cImGrad }; template cImGrad Deriche(const cDataIm2D &aImIn,double aAlpha); +template void ComputeDeriche(cImGrad & aResGrad,const cDataIm2D & aImIn,double aAlpha); + }; diff --git a/MMVII/include/MMVII_PCSens.h b/MMVII/include/MMVII_PCSens.h index 4667fea779..10735f840f 100755 --- a/MMVII/include/MMVII_PCSens.h +++ b/MMVII/include/MMVII_PCSens.h @@ -250,6 +250,9 @@ class cPerspCamIntrCalib : public cObj2DelAtEnd, dont see what can be done ... */ tPtOut Undist(const tPtOut &) const; + /** Inverse function of Undist ... */ + tPtOut Redist(const tPtOut &) const; + // ================== Accessors & Modifiers =================== diff --git a/MMVII/include/MMVII_enums.h b/MMVII/include/MMVII_enums.h index 4a6156a506..a539829d91 100755 --- a/MMVII/include/MMVII_enums.h +++ b/MMVII/include/MMVII_enums.h @@ -31,6 +31,10 @@ enum class eTA2007 FileImage, ///< File containing an image FileCloud, ///< File containing a cloud file (ply ?) File3DRegion, ///< File containing a 3D region + FileTagged, ///< File containing a "xml" or "json" extension + FileTxt, ///< Text file, no extension specified + FileAny, ///< Any file, no more specificiation can be given + FolderAny, ///< Any folder, no more specificiation can be given MPatFile, ///< Major PaternIm => "" or "0" in sem for set1, "1" or other for set2 FFI, ///< File Filter Interval Orient, ///< Orientation diff --git a/MMVII/include/MMVII_util_tpl.h b/MMVII/include/MMVII_util_tpl.h index 1f67dc7462..7e4420466d 100755 --- a/MMVII/include/MMVII_util_tpl.h +++ b/MMVII/include/MMVII_util_tpl.h @@ -91,6 +91,34 @@ template cSelector Str2Interv(const std::string & aStr); /* */ /* ============================================= */ +template class cTransformator +{ + public : + virtual Type Transfo(const Type &) const = 0; +}; + +template class cIdTransformator : public cTransformator +{ + public : + Type Transfo(const Type &) const override; +}; + +/** Transformation by pattern of regular expression */ +class cPatternTransfo : public cTransformator +{ + public : + cPatternTransfo(const std::string & aPat,const std::string & aSubst); + /// Can be 2 vect or empty + cPatternTransfo(const std::vector & aPat); + + std::string Transfo(const std::string &) const override; + private : + std::string mPat; + std::string mSubst; +}; + + + /// Bench some sets functionnalities void BenchSet(const std::string & aDir); @@ -125,6 +153,7 @@ template class cExtSet : public cSelector bool Suppress(const Type &) ; void clear() ; int size() const ; + void Filter(const cSelector &,const cTransformator&); void Filter(const cSelector &); virtual void PutInVect(std::vector &,bool Sorted) const ; ///< Some type requires iteration diff --git a/MMVII/include/V1VII.h b/MMVII/include/V1VII.h index 90ece7673f..77a276a9f5 100755 --- a/MMVII/include/V1VII.h +++ b/MMVII/include/V1VII.h @@ -70,6 +70,9 @@ std::string V1NameMasqOfIm(const std::string & aName); // Call V1 Fast kth value extraction double KthVal(std::vector &, double aProportion); +// Idem but indicate a number and not a proportion +double IKthVal(std::vector & aV, int aK); + // Call V1 for roots of polynomials template std::vector V1RealRoots(const std::vector & aVCoef, Type aTol,int aNbMaxIter); diff --git a/MMVII/include/cMMVII_Appli.h b/MMVII/include/cMMVII_Appli.h index 629242a1ba..39c71db68a 100755 --- a/MMVII/include/cMMVII_Appli.h +++ b/MMVII/include/cMMVII_Appli.h @@ -647,6 +647,7 @@ class cMMVII_Appli : public cMMVII_Ap_NameManip, bool mRMSWasUsed; ///< Indicate if MultiCall was used std::string mIntervFilterMS[NbMaxMainSets]; ///< Filterings interval + std::vector mTransfoFFI[NbMaxMainSets]; ///< Pattern of transformation for FFI // Variable for setting num of mm version for output int mNumOutPut; ///< specified by user diff --git a/MMVII/src/Appli/cMMVII_Appli.cpp b/MMVII/src/Appli/cMMVII_Appli.cpp index 350256d0aa..0957483da0 100755 --- a/MMVII/src/Appli/cMMVII_Appli.cpp +++ b/MMVII/src/Appli/cMMVII_Appli.cpp @@ -430,13 +430,6 @@ void cMMVII_Appli::InitParam(cGenArgsSpecContext *aArgsSpecs) cSpecOneArg2007::tAllSemPL aGlob{eTA2007::Global}; // just to make shorter lines cSpecOneArg2007::tAllSemPL aGlobHDV{eTA2007::Global,eTA2007::HDV}; // just to make shorter lines - - /* Decoding AOpt2007(mIntervFilterMS[0],GOP_Int0,"File Filter Interval, Main Set" ,{eTA2007::Common,{eTA2007::FFI,"0"}}) - mIntervFilterMS[0] => string member, will store the value - GOP_Int0 => const name, Global Optionnal Interval , num 0, declared in MMVII_DeclareCste.h - {eTA2007::Common,{eTA2007::FFI,"0"}} attibute, it's common, it's intervall with attribute "0" - */ - if (HasSharedSPO(eSharedPO::eSPO_CarPO)) { mArgFac << AOpt2007(mCarPPrefOut,"CarPOut","Name for Output caracteristic points",{eTA2007::HDV}); @@ -457,13 +450,17 @@ void cMMVII_Appli::InitParam(cGenArgsSpecContext *aArgsSpecs) TestMainSet(anArgObl,HasMain0,HasMain1); TestMainSet(anArgFac,HasMain0,HasMain1); if (HasMain0) + { mArgFac << AOpt2007(mIntervFilterMS[0],GOP_Int0,"File Filter Interval, Main Set" ,{eTA2007::Shared,{eTA2007::FFI,"0"}}); + mArgFac << AOpt2007(mTransfoFFI[0],"Pat"+GOP_Int0,"Pattern Transfo File Filter Interval, Main Set" ,{eTA2007::Shared}); + } if (HasMain1) + { mArgFac << AOpt2007(mIntervFilterMS[1],GOP_Int1,"File Filter Interval, Second Set",{eTA2007::Shared,{eTA2007::FFI,"1"}}); + mArgFac << AOpt2007(mTransfoFFI[1],"Pat"+GOP_Int1,"Pattern Transfo File Filter Interval, Main Set" ,{eTA2007::Shared}); + } } mArgFac - // << AOpt2007(mIntervFilterMS[0],GOP_Int0,"File Filter Interval, Main Set" ,{eTA2007::Common,{eTA2007::FFI,"0"}}) - // << AOpt2007(mIntervFilterMS[1],GOP_Int1,"File Filter Interval, Second Set",{eTA2007::Common,{eTA2007::FFI,"1"}}) << AOpt2007(mNumOutPut,GOP_NumVO,"Num version for output format (1 or 2)",{eTA2007::Global,{eTA2007::Range,"[1,2]"}}) << AOpt2007(mSeedRand,GOP_SeedRand,"Seed for random,if <=0 init from time",aGlobHDV) << AOpt2007(mExtandPattern,"ExtPatFile","Do we extang patterns for files (or interpret them literally)",aGlobHDV) @@ -781,10 +778,11 @@ void cMMVII_Appli::InitParam(cGenArgsSpecContext *aArgsSpecs) // Filter with interval { - std::string & aNameInterval = mIntervFilterMS[aNum]; + const std::string & aNameInterval = mIntervFilterMS[aNum]; if (IsInit(&aNameInterval)) { - mVMainSets.at(aNum).Filter(Str2Interv(aNameInterval)); + cPatternTransfo aPat(mTransfoFFI[aNum]); + mVMainSets.at(aNum).Filter(Str2Interv(aNameInterval),aPat); } } // Test non empty @@ -806,13 +804,6 @@ void cMMVII_Appli::InitParam(cGenArgsSpecContext *aArgsSpecs) { MMVII_INTERNAL_ASSERT_always(false,"Multiple main set im for num:"+ToStr(aNum)); } -/* - std::string & aNameInterval = mIntervFilterMS[aNum]; - if (IsInit(&aNameInterval)) - { - mVMainSets.at(aNum).Filter(Str2Interv(aNameInterval)); - } -*/ } } // Check validity of main set initialization diff --git a/MMVII/src/Appli/cSpecMMVII_Appli.cpp b/MMVII/src/Appli/cSpecMMVII_Appli.cpp index 98b1c97eb0..cff9e4b909 100755 --- a/MMVII/src/Appli/cSpecMMVII_Appli.cpp +++ b/MMVII/src/Appli/cSpecMMVII_Appli.cpp @@ -240,6 +240,7 @@ std::vector & cSpecMMVII_Appli::InternVecAll() TheVecAll.push_back(&TheSpec_ChSysCoGCP); TheVecAll.push_back(&TheSpec_TutoSerial); TheVecAll.push_back(&TheSpec_TutoFormalDeriv); + TheVecAll.push_back(&TheSpecAppliExtractLine); std::sort(TheVecAll.begin(),TheVecAll.end(),CmpCmd); } diff --git a/MMVII/src/ConvertFormat/ImportGCP.cpp b/MMVII/src/ConvertFormat/ImportGCP.cpp index 879d53b407..405bcaabd3 100644 --- a/MMVII/src/ConvertFormat/ImportGCP.cpp +++ b/MMVII/src/ConvertFormat/ImportGCP.cpp @@ -44,6 +44,7 @@ class cAppli_ImportGCP : public cMMVII_Appli int mComment; int mNbDigName; std::vector mPatternTransfo; + double mMulCoord; }; cAppli_ImportGCP::cAppli_ImportGCP(const std::vector & aVArgs,const cSpecMMVII_Appli & aSpec) : @@ -51,7 +52,8 @@ cAppli_ImportGCP::cAppli_ImportGCP(const std::vector & aVArgs,const mPhProj (*this), mL0 (0), mLLast (-1), - mComment (-1) + mComment (-1), + mMulCoord (1.0) { } @@ -73,6 +75,7 @@ cCollecSpecArg2007 & cAppli_ImportGCP::ArgOpt(cCollecSpecArg2007 & anArgObl) << AOpt2007(mLLast,"NumLast","Num of last line to read (-1 if at end of file)",{eTA2007::HDV}) << AOpt2007(mPatternTransfo,"PatName","Pattern for transforming name (first sub-expr)",{{eTA2007::ISizeV,"[2,2]"}}) << mPhProj.ArgChSys(true) // true => default init with None + << AOpt2007(mMulCoord,"MulCoord","Coordinate multiplier, used to change unity as meter to mm") ; } @@ -127,7 +130,7 @@ int cAppli_ImportGCP::Exe() std::string aAdditionalInfo = ""; if (aHasAdditionalInfo) aAdditionalInfo = aVNames.at(aK).at(aRankA); - aSetM.AddMeasure(cMes1GCP(aChSys.Value(aVXYZ[aK]),aName,1.0,aAdditionalInfo)); + aSetM.AddMeasure(cMes1GCP(aChSys.Value(aVXYZ[aK]*mMulCoord),aName,1.0,aAdditionalInfo)); } mPhProj.SaveGCP(aSetM); diff --git a/MMVII/src/ConvertFormat/V1ImportGCP.cpp b/MMVII/src/ConvertFormat/V1ImportGCP.cpp index ada479134c..f917df7f0a 100644 --- a/MMVII/src/ConvertFormat/V1ImportGCP.cpp +++ b/MMVII/src/ConvertFormat/V1ImportGCP.cpp @@ -41,8 +41,8 @@ cAppli_ConvertV1V2_GCPIM::cAppli_ConvertV1V2_GCPIM(const std::vector & aDI1,int aNbIt,double aStdDev) template cImGrad::cImGrad(const cIm2D & aGx,const cIm2D & aGy) : - mGx (aGx), - mGy (aGy) + mGx (aGx), + mDGx (&mGx.DIm()), + mGy (aGy), + mDGy (&mGy.DIm()) { } template cImGrad:: cImGrad(const cPt2di & aSz) : diff --git a/MMVII/src/ImagesInfoExtract/cAppliExtractLine.cpp b/MMVII/src/ImagesInfoExtract/cAppliExtractLine.cpp new file mode 100755 index 0000000000..f07bd40ce2 --- /dev/null +++ b/MMVII/src/ImagesInfoExtract/cAppliExtractLine.cpp @@ -0,0 +1,399 @@ +#include "MMVII_PCSens.h" +#include "MMVII_ImageInfoExtract.h" +#include "MMVII_ExtractLines.h" + + +namespace MMVII +{ + +cHoughPS::cHoughPS(const cHoughTransform * aHT,const cPt2dr & aTR,tREAL8 aCumul,const cPt2dr & aP1,const cPt2dr & aP2) : + mHT (aHT), + mTetaRho (aTR), + mCumul (aCumul), + mSegE (aP1,aP2) +{ + InitMatch(); +} + + + + +tREAL8 cHoughPS::DistAnglAntiPar(const cHoughPS& aPS2) const +{ + return diff_circ(Teta()+M_PI,aPS2.Teta(),2*M_PI); +} + +tREAL8 cHoughPS::DY(const cHoughPS & aHPS) const +{ + return mSegE.ToCoordLoc(aHPS.mSegE.PMil()) .y() ; +} + +const cSegment2DCompiled & cHoughPS::Seg() const {return mSegE;} +const cPt2dr & cHoughPS::TetaRho() const {return mTetaRho;} +const tREAL8 & cHoughPS::Teta() const {return mTetaRho.x();} +const tREAL8 & cHoughPS::Rho() const {return mTetaRho.y();} +const tREAL8 & cHoughPS::Cumul() const {return mCumul;} +cHoughPS * cHoughPS::Matched() const {return mMatched;} + +cPt2dr cHoughPS::IndTetaRho() const +{ + return cPt2dr(mHT->Teta2RInd(Teta()),mHT->Rho2RInd(Rho())); +} + +tREAL8 cHoughPS::Dist(const cHoughPS& aPS2,const tREAL8 &aFactTeta) const +{ + return std::abs(Rho()-aPS2.Rho()) + + DistAnglAntiPar(aPS2) * (aFactTeta * mHT->RhoMax()) + ; +} + +void cHoughPS::Test(const cHoughPS & aHPS) const +{ + StdOut() << "LARG " << mSegE.ToCoordLoc(aHPS.mSegE.PMil()) .y() + << " RMDistAng " << DistAnglAntiPar(aHPS) * mHT->RhoMax() + << " DistAng " << DistAnglAntiPar(aHPS) + // << " T1=" << Teta()<< " T2=" aHPS.Teta() << " "<< diff_circ(Teta(),aHPS.Teta(),2*M_PI) + << "\n"; +} + +bool cHoughPS::Match(const cHoughPS & aPS2,bool IsLight,tREAL8 aMaxTeta,tREAL8 aDMin,tREAL8 aDMax) const +{ + tREAL8 aDY1 = DY(aPS2 ); + tREAL8 aDY2 = aPS2.DY(*this); + + // Test the coherence of left/right position + if ( ((aDY1>0) ==IsLight) || ((aDY2>0) ==IsLight)) return false; + + if ( DistAnglAntiPar(aPS2) > aMaxTeta) return false; + + tREAL8 aDYAbs1 = std::abs(aDY1); + tREAL8 aDYAbs2 = std::abs(aDY2); + + return (aDYAbs1>aDMin) && (aDYAbs1aDMin) && (aDYAbs2 & mVPS,bool IsLight,tREAL8 aMaxTeta,tREAL8 aDMin,tREAL8 aDMax) +{ + // Reset matches + for (auto & aPtr : mVPS) + aPtr->InitMatch(); + + // compute best match + for (size_t aK1=0 ; aK1Match(*mVPS[aK2],IsLight,aMaxTeta,aDMin,aDMax)) + { + tREAL8 aD12 = mVPS[aK1]->Dist(*mVPS[aK2],2.0); + mVPS[aK1]->UpdateMatch(mVPS[aK2],aD12); + mVPS[aK2]->UpdateMatch(mVPS[aK1],aD12); + } + } + } + + // Test if reciproc match + for (auto & aPtr : mVPS) + { + if (aPtr->mMatched && (aPtr->mMatched->mMatched!=aPtr)) + { + aPtr->mMatched->mMatched =nullptr; + aPtr->mMatched =nullptr; + } + } +} + +/* =============================================== */ +/* */ +/* cAppliExtractLine */ +/* */ +/* =============================================== */ + + +/** An application for computing line extraction. Created for CERN + * porject, and parametrization is more or less optimized for this purpose. + */ + +class cAppliExtractLine : public cMMVII_Appli +{ + public : + cAppliExtractLine(const std::vector & aVArgs,const cSpecMMVII_Appli & aSpec); + private : + typedef tREAL4 tIm; + + cPt2dr Redist(const cPt2dr &) const; + + int Exe() override; + cCollecSpecArg2007 & ArgObl(cCollecSpecArg2007 & anArgObl) override ; + cCollecSpecArg2007 & ArgOpt(cCollecSpecArg2007 & anArgOpt) override ; + std::vector Samples() const override; + virtual ~cAppliExtractLine(); + + void DoOneImage(const std::string & aNameIm) ; + + void MakeVisu(const std::string & aNameIm); + + cPhotogrammetricProject mPhProj; + std::string mPatImage; + bool mLineIsWhite; + bool mShowSteps; + std::vector mVParams; + cPerspCamIntrCalib * mCalib; + bool mAffineMax; + std::vector mThreshCpt; + tREAL8 mTransparencyCont; + cExtractLines* mExtrL; + int mZoomImL; + std::vector mVPS; + std::string mNameReport; + tREAL8 mRelThrsCum; +}; + + +cAppliExtractLine::cAppliExtractLine(const std::vector & aVArgs,const cSpecMMVII_Appli & aSpec) : + cMMVII_Appli (aVArgs,aSpec), + mPhProj (*this), + mShowSteps (true), + mVParams { 4e-4,2.0,5.0}, + mCalib (nullptr), + mAffineMax (true), + mThreshCpt {100,200,400,600}, + mTransparencyCont (0.5), + mExtrL (nullptr), + mZoomImL (1), + mNameReport ("LineExtract"), + mRelThrsCum (0.3) +{ +} + +cAppliExtractLine::~cAppliExtractLine() +{ + delete mExtrL; + DeleteAllAndClear(mVPS); +} + +cCollecSpecArg2007 & cAppliExtractLine::ArgObl(cCollecSpecArg2007 & anArgObl) +{ + return anArgObl + << Arg2007(mPatImage,"Name of input Image", {eTA2007::FileDirProj,{eTA2007::MPatFile,"0"}}) + << Arg2007(mLineIsWhite," True : its a light line , false dark ") + ; +} + +cCollecSpecArg2007 & cAppliExtractLine::ArgOpt(cCollecSpecArg2007 & anArgOpt) +{ + return anArgOpt + << mPhProj.DPOrient().ArgDirInOpt("","Folder for calibration to integrate distorsion") + << AOpt2007(mAffineMax,"AffineMax","Affinate the local maxima",{eTA2007::HDV}) + << AOpt2007(mShowSteps,"ShowSteps","Show detail of computation steps by steps",{eTA2007::HDV}) + << AOpt2007(mZoomImL,"ZoomImL","Zoom for images of line",{eTA2007::HDV}) + << mPhProj.DPPointsMeasures().ArgDirInOpt("","Folder for ground truth measure") + ; +} + +std::vector cAppliExtractLine::Samples() const +{ + return { + "MMVII ExtractLine 'DSC_.*.JPG' ShowSteps=1 InOri=FB" + }; +} + +cPt2dr cAppliExtractLine::Redist(const cPt2dr & aP) const +{ + return mCalib ? mCalib->Redist(aP) : aP; +} + + +void cAppliExtractLine::DoOneImage(const std::string & aNameIm) +{ + tREAL8 aMulTeta = 1.0/M_PI; +// aMulTeta = 1.0; + bool mShow = true; + mCalib = nullptr; + if (mPhProj.DPOrient().DirInIsInit()) + mCalib = mPhProj.InternalCalibFromImage(aNameIm); + + cIm2D anIm = cIm2D::FromFile(aNameIm); + tREAL8 aTrhsCum = mRelThrsCum * Norm2(anIm.DIm().Sz()); + mExtrL = new cExtractLines (anIm); + + mExtrL->SetDericheGradAndMasq(2.0,10.0,2,mShow); // aAlphaDerich,aRayMaxLoc,aBorder + mExtrL->SetHough(cPt2dr(aMulTeta,1.0),0.1,mCalib,mAffineMax,mShow); + + std::vector aVMaxLoc = mExtrL->Hough().ExtractLocalMax(10,4.0,10.0,0.1); + for (const auto & aPMax : aVMaxLoc) + { + cHoughPS * aPS = mExtrL->Hough().PtToLine(aPMax); + if (aPS->Cumul() > aTrhsCum) + mVPS.push_back(aPS); + else + delete aPS; + } + // mVPS[0]->Test(*mVPS[1]); + // mVPS[1]->Test(*mVPS[0]); + + cHoughPS::SetMatch(mVPS,mLineIsWhite,mVParams.at(0),mVParams.at(1),mVParams.at(2)); + + if (mShowSteps) + MakeVisu(aNameIm); +} + +void cAppliExtractLine::MakeVisu(const std::string & aNameIm) +{ + const cDataIm2D& aDAccum =mExtrL->Hough().Accum().DIm(); + tREAL8 aVMax=0.0; + + // [0] Print some statistic to compared the compatcness of histogramm + { + std::vector aVCpt(mThreshCpt.size(),0); // count the number of point over a threshold + for (const auto & aPix : aDAccum) + { + tREAL8 aV = aDAccum.GetV(aPix); + UpdateMax(aVMax,aV); + for (size_t aK=0 ; aKmThreshCpt.at(aK)) + aVCpt.at(aK)++; + } + } + + // print max value + StdOut() << "VMAX=" << aVMax << std::endl; + for (size_t aK=0 ; aKMakeImageMaxLoc(mTransparencyCont); + aImV.ToJpgFileDeZoom(mPhProj.DirVisu() + "DetectL_"+ aNameTif,1); + } + + // [2] Visu module of gradient + { + std::string aNameGrad = mPhProj.DirVisu()+"Grad_" + aNameTif; + mExtrL->Grad().NormG().DIm().ToFile(aNameGrad); + // Convert_JPG(aNameGrad,true,90,"jpg"); + } + + + // [3] Visu the accum + local maximal + { + cRGBImage aVisAccum = RGBImFromGray(aDAccum,255.0/aVMax); + int aK=0; + for (const auto & aPS : mVPS) + { + //aVisAccum.SetRGBrectWithAlpha(ToI(aPS->IndTetaRho()) ,15,cRGBImage::Red,0.5); + cPt2dr aC=aPS->IndTetaRho(); + cPt2dr aSz(7,7); + aVisAccum.FillRectangle(cRGBImage::Red,ToI(aC-aSz),ToI(aC+aSz),cPt3dr(0.5,1.0,1.0)); + aK++; + } + aVisAccum.ToJpgFileDeZoom(mPhProj.DirVisu() + "Accum_" + aNameTif,1); + aDAccum.ToFile(mPhProj.DirVisu() + "RawAccum_" + aNameTif); + } + // [4] Visu of Image + Lines + { + cRGBImage aVisIm = cRGBImage::FromFile(aNameIm,mZoomImL); // Initialize with image + const auto & aDIm = aVisIm.ImR().DIm(); + for (size_t aKH=0 ; aKHMatched(); + // Compute colour, depend of ranking + bool isOk = (aHS2 != nullptr); + cPt3di aCoul = isOk ? cRGBImage::Red : cRGBImage::Blue ; + // if (aKH>=2) aCoul = cRGBImage::Green; + // if (aKH>=4) aCoul = cRGBImage::Blue; + if (isOk && (aHS1> aHS2)) + { + tREAL8 aDAng = aHS1->DistAnglAntiPar(*aHS2) * mExtrL->Hough().RhoMax() ; + tREAL8 aLarg = aHS1->DY(*aHS2) ; + tREAL8 aCumul = (aHS1->Cumul()+aHS2->Cumul())/2.0; + AddOneReportCSV(mNameReport,{aNameIm,ToStr(aDAng),ToStr(aLarg),ToStr(aCumul)}); + } + + // Compute Hough-Point -> line + cSegment aSeg = mVPS[aKH]->Seg(); + + tREAL8 aRay=5; + cPt2dr aC = Redist(aSeg.PMil() +VUnit(aSeg.V12()*cPt2dr(0,-1))* aRay); + aVisIm.DrawCircle(aCoul,aC ,aRay); + + // write point by point, in two direction + for (tREAL8 aSign : {-1.0,1.0}) + { + cPt2dr aPt = aSeg.PMil(); + cPt2dr aTgt = VUnit(aSeg.V12()) * aSign; // direction + cPt2dr aQ = Redist(aPt) ; // eventulay make invert distorsion correcion + while (aDIm.InsideBL(aQ)) + { + aVisIm.SetRGBPoint(aQ,aCoul); // print point + aPt = aPt+aTgt; // increment point + aQ = Redist(aPt); // eventulay make invert distorsion correcion + } + } + } + std::string aNameLine = mPhProj.DirVisu() + "Lines_" + aNameTif; + // convert dont handle well big file, so generate jpg only if zoom=1 + if (mZoomImL==1) + aVisIm.ToJpgFileDeZoom(mPhProj.DirVisu() + "Lines_" + aNameTif,1); + else + aVisIm.ToFile(aNameLine); + + } +} + + + +int cAppliExtractLine::Exe() +{ + mPhProj.FinishInit(); + InitReport(mNameReport,"csv",true); + + if (RunMultiSet(0,0)) + { + AddOneReportCSV(mNameReport,{"NameIm","Paral","Larg","Cumul"}); + return ResultMultiSet(); + } + DoOneImage(UniqueStr(0)); + return EXIT_SUCCESS; +} + + +tMMVII_UnikPApli Alloc_AppliExtractLine(const std::vector & aVArgs,const cSpecMMVII_Appli & aSpec) +{ + return tMMVII_UnikPApli(new cAppliExtractLine(aVArgs,aSpec)); +} + + +cSpecMMVII_Appli TheSpecAppliExtractLine +( + "ExtractLine", + Alloc_AppliExtractLine, + "Extraction of lines", + {eApF::Ori}, + {eApDT::Ori,eApDT::GCP}, + {eApDT::Console}, + __FILE__ +); + +}; diff --git a/MMVII/src/ImagesInfoExtract/cExtractLines.cpp b/MMVII/src/ImagesInfoExtract/cExtractLines.cpp new file mode 100755 index 0000000000..a0af97c758 --- /dev/null +++ b/MMVII/src/ImagesInfoExtract/cExtractLines.cpp @@ -0,0 +1,166 @@ +#include "MMVII_PCSens.h" +#include "MMVII_ImageInfoExtract.h" +#include "MMVII_ExtractLines.h" + + +namespace MMVII +{ + +/* ************************************************************************ */ +/* */ +/* cExtractLines */ +/* */ +/* ************************************************************************ */ + + + +template cExtractLines::cExtractLines(tIm anIm) : + mSz (anIm.DIm().Sz()), + mIm (anIm), + mImMasqCont (mSz,nullptr,eModeInitImage::eMIA_Null), + mGrad (nullptr), + mHough (nullptr), + mCalib (nullptr) +{ +} + +template cExtractLines::~cExtractLines() +{ + delete mGrad; + delete mHough; +} + +template void cExtractLines::SetHough + ( + const cPt2dr & aMulTetaRho, + tREAL8 aSigmTeta, + cPerspCamIntrCalib * aCalib, + bool AffineMax, + bool Show + ) +{ + // Memorize calbration and initialize hough structure + mCalib = aCalib; + mHough = new cHoughTransform(ToR(mSz),aMulTetaRho,aSigmTeta,aCalib); + + // compute average for weighting + tREAL8 aAvgIm=0; + for (const auto & aPix : mImMasqCont.DIm()) + aAvgIm += mGrad->NormG().DIm().GetV(aPix); + aAvgIm /= mGrad->NormG().DIm().NbElem() ; + + // Three measure of correction (tuning) + tREAL8 aSomCorAff=0; // sums the distance between point and its correction by refinement + tREAL8 aSomCorDist=0; // sums the distance between point and its correction by distorsion + tREAL8 aSomCorTeta=0; // sums the angulare distance of distorsion + int aNbDone=0; + + for (const auto & aPix : mImMasqCont.DIm()) + { + if ( mImMasqCont.DIm().GetV(aPix)) + { + if (Show && ((aNbDone%200000)==0) ) + StdOut() << "Remain to do " << mNbPtsCont-aNbDone << "\n"; + aNbDone++; + + // compute eventually refined position + cPt2dr aRPix0 = ToR(aPix); + cPt2dr aRPix = aRPix0; + if (AffineMax) + aRPix = mGrad->RefinePos(aRPix0); + aSomCorAff += Norm2(aRPix0-aRPix); + + // compute Teta of grad, for limiting the computation time + tREAL8 aTeta = Teta(mGrad->Grad(aPix)); + + // if calibration, correct position an teta + if (aCalib) + { + cPt2dr aCor = aCalib->Undist(aRPix); // point correction + cPt2dr aCor2 = aCalib->Undist(aRPix+ FromPolar(0.1,aTeta) ); // point in gradient direction + tREAL8 aTetaCorr = Teta(aCor2-aCor); // new teta + + aSomCorDist += Norm2(aCor-aRPix); // sum correction on point due to distorsion + aSomCorTeta += std::abs(aTetaCorr-aTeta); // sum correction on angle due to distorsion + aRPix = aCor; // set position to corrected + aTeta = aTetaCorr; // set teta to corrected + } + + // compute a weighting, growing with norm, but limited to 1.0 + tREAL8 aNorm = mGrad->NormG().DIm().GetV(aPix); + tREAL8 aW = aNorm / (aNorm + (2.0*aAvgIm)); + + // finnaly add the point to hough-accumulator + if (mImMasqCont.DIm().InsideBL(aRPix)) + mHough->AccumulatePtAndDir(aRPix,aTeta,aW); + } + } + // make some filter, not sure usefull + ExpFilterOfStdDev(mHough->Accum().DIm(),4,1.0); + + if (Show) + { + StdOut() + << " , Aff=" << (aSomCorAff/aNbDone) + << " , Dist-Pt=" << (aSomCorDist/aNbDone) + << " , Dist-Teta=" << (aSomCorTeta/aNbDone) + << std::endl; + } +} + +template void cExtractLines::SetDericheGradAndMasq(tREAL8 aAlpha,tREAL8 aRay,int aBorder,bool Show) +{ + // Create the data for storing gradient & init gradient + mGrad = new cImGradWithN(mIm.DIm(),aAlpha); + + cRect2 aRect(mImMasqCont.DIm().Dilate(-aBorder)); // rect interior + std::vector aVecNeigh = cImGradWithN::NeighborsForMaxLoc(aRay); // neigbours for compute max + + // count pts & pts of contour for stat + mNbPtsCont = 0; + int aNbPt =0; + // Parse all points to set the masq if is local maxima in direction of gradient + for (const auto & aPix : aRect) + { + aNbPt++; + if (mGrad->IsMaxLocDirGrad(aPix,aVecNeigh,1.0)) // aPix,Neigbours,aRatioXY + { + mImMasqCont.DIm().SetV(aPix,255); + mNbPtsCont++; + } + } + + if (Show) + StdOut()<< " Prop Contour = " << mNbPtsCont / double(aNbPt) << "\n"; +} + +/* Generate a RGB-image : + * - background is initial image + * - point of contour are set to red with alpha transparency + */ +template cRGBImage cExtractLines::MakeImageMaxLoc(tREAL8 aAlpha) +{ + cRGBImage aImV(mIm.DIm().Sz()); // init RGB with size + for (const auto & aPix : mImMasqCont.DIm()) + { + aImV.SetGrayPix(aPix,mIm.DIm().GetV(aPix)); // transfer image + // set contour + if (mImMasqCont.DIm().GetV(aPix)) + { + aImV.SetRGBPixWithAlpha(aPix,cRGBImage::Red,cPt3dr(aAlpha,aAlpha,aAlpha)); + } + } + return aImV; +} + + + +template cHoughTransform & cExtractLines::Hough() {return *mHough;} +template cImGradWithN & cExtractLines::Grad() {return *mGrad;} + + +// ========================= INSTANCIATION =============== + +template class cExtractLines; + +}; diff --git a/MMVII/src/ImagesInfoExtract/cHoughTransform.cpp b/MMVII/src/ImagesInfoExtract/cHoughTransform.cpp new file mode 100755 index 0000000000..25dcd48d8f --- /dev/null +++ b/MMVII/src/ImagesInfoExtract/cHoughTransform.cpp @@ -0,0 +1,191 @@ +#include "MMVII_Tpl_Images.h" +#include "MMVII_PCSens.h" +#include "MMVII_ImageInfoExtract.h" +#include "MMVII_ExtractLines.h" + + +namespace MMVII +{ + +/* ************************************************************************ */ +/* */ +/* cHoughTransform */ +/* */ +/* ************************************************************************ */ + + +cHoughTransform::cHoughTransform +( + const cPt2dr & aSzIn, + const cPt2dr & aMulTetaRho, + const tREAL8 & aSigmTeta, + cPerspCamIntrCalib * aCalib +) : + mMiddle (aSzIn / 2.0), + mRhoMax (Norm2(mMiddle)), + mMulTeta (aMulTetaRho.x()), + mMulRho (aMulTetaRho.y()), + mSigmTeta (aSigmTeta), + mCalib (aCalib), + mNbTeta (round_up(2*M_PI*mMulTeta * mRhoMax)), + mFactI2T ((2.0*M_PI)/mNbTeta), + mNbRho (2+round_up(2*mMulRho * mRhoMax)), + mTabSin (mNbTeta), + mDTabSin (mTabSin.DIm()), + mTabCos (mNbTeta), + mDTabCos (mTabCos.DIm()), + mAccum (cPt2di(mNbTeta,mNbRho),nullptr,eModeInitImage::eMIA_Null), + mDAccum (mAccum.DIm()) +{ + // Tabulate "cos" and "sin" + for (int aKTeta=0 ; aKTeta cHoughTransform::Accum() const {return mAccum;} +const tREAL8 & cHoughTransform::RhoMax() const {return mRhoMax;} + +tREAL8 cHoughTransform::AvgS2() const +{ + return std::sqrt(DotProduct(mDAccum,mDAccum) / mDAccum.NbElem()); +} + +cHoughPS* cHoughTransform::PtToLine(const cPt3dr & aPt) const +{ + // (x,y) in L <=> x cos(T) + y Sin(T) = R + tREAL8 aTeta = RInd2Teta(aPt.x()); + tREAL8 aRho = RInd2Rho(aPt.y()); + cPt2dr aTgt(-sin(aTeta),cos(aTeta)); + cPt2dr aP0 = mMiddle + cPt2dr(aRho*cos(aTeta),aRho*sin(aTeta)); + + return new cHoughPS(this,cPt2dr(aTeta,aRho),aPt.z(),aP0-aTgt,aP0+aTgt); +} + + +void cHoughTransform::AccumulatePtAndDir(const cPt2dr & aPt,tREAL8 aTetaC,tREAL8 aWeight) +{ + // compute the index-interval, centered on aTetaC, of size mSigmTeta + int iTeta0 = round_down(Teta2RInd(aTetaC-mSigmTeta)); + int iTeta1 = round_up(Teta2RInd(aTetaC+mSigmTeta)); + + // Angle are defined %2PI, make interval > 0 + if (iTeta0<0) + { + iTeta0 += mNbTeta; + iTeta1 += mNbTeta; + aTetaC += 2 * M_PI; + } + // Polar coordinate Rho-Teta, are defined arround point mMiddle + tREAL8 aX = aPt.x() - mMiddle.x(); + tREAL8 aY = aPt.y() - mMiddle.y(); + + for (int iTeta=iTeta0 ; iTeta<=iTeta1 ; iTeta++) + { + tREAL8 aTeta = Ind2Teta(iTeta); + // Compute a weighting function : 1 at aTetaC, 0 if dist > mSigmTeta + tREAL8 aWTot = aWeight * ( 1 - std::abs(aTeta-aTetaC) /mSigmTeta); + if (aWTot>0) // (weitgh can be <0 because of rounding) + { + int aITetaOK = iTeta%mNbTeta; // % 2PI => put indexe in interval + // equation : (x,y) in L <=> x cos(T) + y Sin(T) = R + tREAL8 aRho = aX * mDTabCos.GetV(aITetaOK) + aY*mDTabSin.GetV(aITetaOK) ; + tREAL8 aRIndRho = Rho2RInd(aRho); // convert Real rho in its index inside accumulator + // Make bi-linear repartition of contribution as "aRIndRho" is real + int aIRho0 = round_down(aRIndRho); + int aIRho1 = aIRho0 +1; + tREAL8 aW0 = (aIRho1-aRIndRho); + + // Finally accumulate + mDAccum.AddVal(cPt2di(aITetaOK,aIRho0), aW0 *aWTot); + mDAccum.AddVal(cPt2di(aITetaOK,aIRho1),(1-aW0)*aWTot); + } + } +} + +tREAL8 cHoughTransform::GetValueBlob(cPt2di aPt,int aMaxNeigh) const +{ + tREAL8 aRes = mDAccum.GetV(aPt); + for (const int aSign : {-1,1}) + { + int aNb =0; + tREAL8 aLastV = mDAccum.GetV(aPt); + cPt2di aPNext = aPt + cPt2di(0,aSign); + + while (mDAccum.Inside(aPNext) && (aNb cHoughTransform::ExtractLocalMax(size_t aNbMax,tREAL8 aDist,tREAL8 aThrAvg,tREAL8 aThrMax) const +{ + // [0] Make a duplication as we will modify the accum + cIm2D anAccum = mAccum.Dup(); + cDataIm2D& aDAccum = anAccum.DIm(); + + // [1] Compute average , max and threshold + // [1.A] : max & avg + tREAL8 aAvg = 0; + tREAL8 aVMax = 0; + for (const auto & aPix : aDAccum) + { + tREAL8 aVal = aDAccum.GetV(aPix); + aAvg += aVal; + UpdateMax(aVMax,aVal); + } + aAvg /= aDAccum.NbElem(); + + // [1.B] Threshlod is the stricter of both + tREAL8 aThrHold = std::max(aAvg*aThrAvg,aVMax*aThrMax); + + + // [2] Set to 0 point < aThrHold (to limitate number of pre-sel & accelerate IKthVal + for (const auto & aPix : aDAccum) + { + if (aDAccum.GetV(aPix) aRes; + cAffineExtremum aAffin(mAccum.DIm(),1.5); + for (const auto aPt : aExtr.mPtsMax) + { + cPt2dr aPAff = aAffin.OneIter(ToR(aPt)); + if ( mDAccum.InsideBL(aPAff)) + aRes.push_back(cPt3dr(aPAff.x(),aPAff.y(),GetValueBlob(aPt,7))); + } + + // [5] Sort with highest value first, then select NbMax + SortOnCriteria(aRes,[](const auto & aP) {return -aP.z();}); + while (aRes.size() > aNbMax) + aRes.pop_back(); + + + return aRes; +} + +}; diff --git a/MMVII/src/ImagesInfoExtract/cImGradWithN.cpp b/MMVII/src/ImagesInfoExtract/cImGradWithN.cpp new file mode 100755 index 0000000000..fbad28132b --- /dev/null +++ b/MMVII/src/ImagesInfoExtract/cImGradWithN.cpp @@ -0,0 +1,143 @@ +#include "MMVII_Linear2DFiltering.h" +#include "MMVII_Geom2D.h" + +#include "MMVII_ExtractLines.h" + +namespace MMVII +{ + +/* ************************************************************************ */ +/* */ +/* cImGradWithN */ +/* */ +/* ************************************************************************ */ + + +template + cImGradWithN::cImGradWithN(const cPt2di & aSz) : + cImGrad (aSz), + mNormG (aSz), + mDataNG (mNormG.DIm()) +{ +} + +template + cImGradWithN::cImGradWithN(const cDataIm2D & aImIn,Type aAlphaDeriche) : + cImGradWithN(aImIn.Sz()) +{ + ComputeDericheAndNorm(*this,aImIn,aAlphaDeriche); +} + + +/* The test is more complicated than traditionnal local maxima : + * + * - 1- On a given edge if we take all the 8-neighboors, one the neighbors in the edge will + * have higher value, that why we take the neigbours that are +- in the direction of grad + * See "NeighborsForMaxLoc()" + * + * - 2- with thin line, the oposite contour may have value that may be higer than the contour itself, + * that's why we consider only the point which are orientd in the same direction (test "Scal > 0"), + * note this work for a dark or light line on a average back-ground + */ +template bool cImGradWithN::IsMaxLocDirGrad(const cPt2di& aPix,const std::vector & aVP,tREAL8 aRatioXY) const +{ + // [1] Compute unitary vector + tREAL8 aN = mDataNG.GetV(aPix); + if (aN==0) return false; // avoid div by 0 + cPt2dr aDirGrad = ToR(this->Grad(aPix)) * (1.0/ aN); + + //[2] A Basic test to reject point on integer neighbourhood in the two direction + { + cPt2di aIDirGrad = ToI(aDirGrad); + for (int aSign : {-1,1}) + { + cPt2di aNeigh = aPix + aIDirGrad * aSign; + if ( (mDataNG.DefGetV(aNeigh,-1)>aN) && (Scal(aDirGrad,ToR(this->Grad(aNeigh)))>0) ) + { + return false; + } + } + } + + // [3] + + cPt2dr aConjDG = conj(aDirGrad); + for (const auto & aDeltaNeigh : aVP) + { + // cPt2dr aDirLoc = ToR(aDeltaNeigh) / aDirGrad; + cPt2dr aDirLoc = ToR(aDeltaNeigh) * aConjDG; + if (std::abs(aDirLoc.x()) >= std::abs(aDirLoc.y()*aRatioXY)) + { + cPt2di aNeigh = aPix + aDeltaNeigh; + + if ( (mDataNG.DefGetV(aNeigh,-1)>aN) && (Scal(aDirGrad,ToR(this->Grad(aNeigh))) >0)) + return false; + } + /*cPt2dr aNeigh = ToR(aPix) + ToR(aDeltaNeigh) * aDirGrad; + if ( (mDataNG.DefGetVBL(aNeigh,-1)>aN) && (Scal(aDirGrad,ToR(this->GradBL(aNeigh))) >0)) + return false; + */ + } + + return true; +} + +// return simply all the point except (0,0) +template std::vector cImGradWithN::NeighborsForMaxLoc(tREAL8 aRay) +{ + return SortedVectOfRadius(0.5,aRay); +} + +/* Make one iteration of position computation : + * + * -1- compute the value of gradient norm befor and after the point in gradient direction + * -2- extract the maximal by parabol fiting + */ + +template cPt2dr cImGradWithN::OneRefinePos(const cPt2dr & aP1) const +{ + if (! mDataNG.InsideBL(aP1)) + return aP1; + + cPt2dr aGr = VUnit(ToR(this->GradBL(aP1))); + + // extract point "before" + cPt2dr aP0 = aP1- aGr; + if (! mDataNG.InsideBL(aP0)) + return aP1; + + // extract point "after" + cPt2dr aP2 = aP1+ aGr; + if (! mDataNG.InsideBL(aP2)) + return aP1; + + // extract abscisse of maxim of parabol + tREAL8 aAbs = StableInterpoleExtr(mDataNG.GetVBL(aP0),mDataNG.GetVBL(aP1),mDataNG.GetVBL(aP2)); + + return aP1 + aGr * aAbs; +} + +// Very basic just max two iterations +template cPt2dr cImGradWithN::RefinePos(const cPt2dr & aP1) const +{ + return OneRefinePos(OneRefinePos(aP1)); +} + + /* ************************************************ */ + +template void ComputeDericheAndNorm(cImGradWithN & aResGrad,const cDataIm2D & aImIn,double aAlpha) +{ + ComputeDeriche(aResGrad,aImIn,aAlpha); + + auto & aDN = aResGrad.NormG().DIm(); + for (const auto & aPix : aDN) + { + aDN.SetV(aPix,Norm2(aResGrad.Grad(aPix))); + } +} + +template class cImGradWithN; +template void ComputeDericheAndNorm(cImGradWithN & aResGrad,const cDataIm2D & aImIn,double aAlpha) ; + +}; + diff --git a/MMVII/src/Instrumental/cBlockCamInit.cpp b/MMVII/src/Instrumental/cBlockCamInit.cpp index edc2134a67..048c3d2098 100644 --- a/MMVII/src/Instrumental/cBlockCamInit.cpp +++ b/MMVII/src/Instrumental/cBlockCamInit.cpp @@ -576,6 +576,8 @@ class cAppli_BlockCamInit : public cMMVII_Appli cCollecSpecArg2007 & ArgObl(cCollecSpecArg2007 & anArgObl) override; cCollecSpecArg2007 & ArgOpt(cCollecSpecArg2007 & anArgOpt) override; + std::vector Samples() const override; + private : std::string mSpecImIn; /// Pattern or xml file cPhotogrammetricProject mPhProj; @@ -590,6 +592,12 @@ class cAppli_BlockCamInit : public cMMVII_Appli bool mTestNoDel; ///< Do force an error on memory management to illustrate the }; +std::vector cAppli_BlockCamInit::Samples() const +{ + return { + "MMVII BlockCamInit SetFiltered_GCP_OK_Resec.xml BA_311_B '(.*)_(.*).JPG' [1,2] Rig_311_B" + }; +} cAppli_BlockCamInit::cAppli_BlockCamInit ( const std::vector & aVArgs, diff --git a/MMVII/src/Instrumental/cClinoInit.cpp b/MMVII/src/Instrumental/cClinoInit.cpp index 7bb7bc55c4..704ff4a95b 100644 --- a/MMVII/src/Instrumental/cClinoInit.cpp +++ b/MMVII/src/Instrumental/cClinoInit.cpp @@ -5,10 +5,9 @@ /** - \file GCPQuality.cpp - - MMVII ClinoInit ClinoMeasures.txt [949_,.JPG] MMVII-PhgrProj/Ori/Resec_311_All/ Nbit0=24 AmplSim=[0.1,0.1,1.5,0.1] SeedRand=22 + \file cClinoInit.cpp + Compute initial value of clinometers boresight relatively to a camera. */ @@ -18,41 +17,52 @@ namespace MMVII typedef cRotation3D tRot; +/** + */ + class cClinoCalMes1Cam { public : - cClinoCalMes1Cam(cSensorCamPC * aCam,const std::vector & aVAngles); + cClinoCalMes1Cam + ( + cSensorCamPC * aCam, // camera-> only the pose is usefuul + const std::vector & aVAngles, // vector of angles of all inclinometer + const cPt3dr & aVerticAbs = {0,0,-1} // position of vertical in current "absolut" system + ); void SetDirSimul(int aK,const cRotation3D &aR) ; - cSensorCamPC * mCam; - std::vector mVDir; - cPt3dr mVertInLoc; /// Vertical in camera coordinates - /// Theoretical Pos of needle, indicating vertical, in plane IJ, for Clino K, given a boresight aCam2Clino - cPt2dr PosNeedle(int aKClino,const cRotation3D &aCam2Clino) const; + cPt2dr PosNeedle(const cRotation3D &aCam2Clino) const; /// Vectorial difference, in plane IJ, between measure an theoreticall value for needle cPt2dr EcarVectRot(int aKClino,const cRotation3D &aCam2Clino) const; - /// Idem but with WPK - cPt2dr EcarVectWPK(int aKClino,const cPt3dr &aWPK) const; + /* => Use of WPK is currently deprecated + /// Idem but with WPK + cPt2dr EcarVectWPK(int aKClino,const cPt3dr &aWPK) const; + /// Idem but with WPK + tREAL8 ScoreWPK(int aKClino,const cPt3dr &aWPK) const; + */ /// Score (to minimize) between obs and theor for given rot tREAL8 ScoreRot(int aKClino,const cRotation3D &aR) const; - /// Idem but with WPK - tREAL8 ScoreWPK(int aKClino,const cPt3dr &aWPK) const; - /// Gradient / wpk of Needle's in x an y in plane IJ - std::pair Grad(int aKClino,const tRot & aRot,tREAL8 aEpsilon=1e-3) const; + /// Gradient / wpk of Needle's "Ecart vectorial" (in x an y) in plane IJ, computed using finite difference (now used for + std::pair GradEVR(int aKClino,const tRot & aRot,tREAL8 aEpsilon=1e-3) const; + private : + cSensorCamPC * mCam; ///< camera , memorization seems useless + std::vector mVDir; ///< measured position of needle, computed from angles + cPt3dr mVertInLoc; ///< Vertical in camera system, this is the only information usefull of camera orientation }; -cClinoCalMes1Cam::cClinoCalMes1Cam(cSensorCamPC * aCam,const std::vector & aVAngles) : +cClinoCalMes1Cam::cClinoCalMes1Cam(cSensorCamPC * aCam,const std::vector & aVAngles,const cPt3dr & aVertAbs) : mCam (aCam), mVertInLoc (mCam->Vec_W2L(cPt3dr(0,0,-1))) { + // tranformate angles in vector position of the needle in plane I,J for (auto & aTeta : aVAngles) { mVDir.push_back(FromPolar(1.0,aTeta)); @@ -60,18 +70,21 @@ cClinoCalMes1Cam::cClinoCalMes1Cam(cSensorCamPC * aCam,const std::vector } -std::pair cClinoCalMes1Cam::Grad(int aKClino,const tRot & aR0,tREAL8 aEps) const +std::pair cClinoCalMes1Cam::GradEVR(int aKClino,const tRot & aR0,tREAL8 aEps) const { cPt3dr aResX; cPt3dr aResY; + // Parse 3 angles for (int aK=0 ; aK<3 ; aK++) { - cPt3dr aPEps(0,0,0); - aPEps[aK] = aEps; + // compute Epsilon vector in 1 direction + cPt3dr aPEps = cPt3dr::P1Coord(aK,aEps); + + // small rotation in direction tRot aRPlus = tRot::RotFromWPK(aPEps); tRot aRMinus = tRot::RotFromWPK(-aPEps); - + // compute difference in needle position cPt2dr aDelta = EcarVectRot(aKClino,aR0*aRPlus) - EcarVectRot(aKClino,aR0*aRMinus); aResX[aK] = aDelta.x()/(2*aEps); aResY[aK] = aDelta.y()/(2*aEps); @@ -79,37 +92,46 @@ std::pair cClinoCalMes1Cam::Grad(int aKClino,const tRot & aR0,tR return std::pair(aResX,aResY); } -cPt2dr cClinoCalMes1Cam::PosNeedle(int aKClino,const cRotation3D &aCam2Clino) const + +cPt2dr cClinoCalMes1Cam::PosNeedle(const cRotation3D &aCam2Clino) const { - cPt3dr aVClin = aCam2Clino.Value(mVertInLoc); // Direction of Vert in clino repair + // "aVClin"= Direction of Vert in clino repair, mVertInLoc being in camera system + // we only need to know the boresight "aCam2Clino" to tranfer it in the clino system + cPt3dr aVClin = aCam2Clino.Value(mVertInLoc); + // VClin =(x,y,z) , the projection of VClin in plane "I,J" is given by Proj (x,y,z) -> (x,y) + // then we make a unitary vector to maintain the size of "needle" return VUnit(Proj(aVClin)); // Projection in plane I,J } cPt2dr cClinoCalMes1Cam::EcarVectRot(int aKClino,const cRotation3D &aCam2Clino) const { - return PosNeedle(aKClino,aCam2Clino) - mVDir[aKClino]; + // Theoretical Pos - Measured Pos + return PosNeedle(aCam2Clino) - mVDir[aKClino]; } tREAL8 cClinoCalMes1Cam::ScoreRot(int aKClino,const cRotation3D& aCam2Clino) const { + // score of a given rotation return SqN2(EcarVectRot(aKClino,aCam2Clino)); } +void cClinoCalMes1Cam::SetDirSimul(int aKClino,const cRotation3D &aR) +{ + mVDir[aKClino] = PosNeedle(aR); +} + +/* tREAL8 cClinoCalMes1Cam::ScoreWPK(int aKClino,const cPt3dr & aWPK) const { return ScoreRot(aKClino,cRotation3D::RotFromWPK(aWPK)); } - cPt2dr cClinoCalMes1Cam::EcarVectWPK(int aKClino,const cPt3dr &aWPK) const { return EcarVectRot(aKClino,cRotation3D::RotFromWPK(aWPK)); } +*/ -void cClinoCalMes1Cam::SetDirSimul(int aKClino,const cRotation3D &aR) -{ - mVDir[aKClino] = PosNeedle(aKClino,aR); -} @@ -139,28 +161,32 @@ class cAppli_ClinoInit : public cMMVII_Appli /// Make one iteration of research arround a current solution cWhichMin OneIter(const tRot & aR0,tREAL8 aStep, int aNbStep); - /// Make a first iteration using sampling by quaternions - cWhichMin IterInit(int aNbStep); + /// Make a first computation parsing "all" rotation at a given step, using sampling by quaternions + cWhichMin IterInit(int aNbStep) const; /// Compute the conditionning of least-square matrix - void ComputeCond(const tRot & aWPK,tREAL8 aEpsilon=1e-3); + cPt2dr ComputeCond(const tRot & aWPK,tREAL8 aEpsilon=1e-3); - cPhotogrammetricProject mPhProj; - std::string mNameClino; /// Pattern of xml file - std::vector mPrePost; /// Pattern of xml file - int mNbStep0; - int mNbStep; - int mNbIter; - std::vector mASim; + cPhotogrammetricProject mPhProj; ///< Classical structure for photogrammetric project + std::string mNameClino; ///< Pattern of xml file + std::vector mPrePost; ///< Pattern of xml file + int mNbStep0; ///< Number of step in initial first round + int mNbStepIter; ///< Number of step in each iteration + int mNbIter; ///< Number of iteration + std::vector mASim; ///< Amplitude of simulation parameters - std::vector mVMeasures; - std::vector mVKClino; - std::vector mComputedClino; /// Clino can be computed globally or independantly + std::vector mVMeasures; ///< store vectors of measure after file parsing + std::vector mVKClino; ///< Vecor of indexes of selected clino + /** If several Clino, they can be computed globally or independantly, this vector of bool allow to have + * the 2 computation, it's mainly for tuning and observing the effect on residual and conditionning + */ + std::vector mComputedClino; - std::string mNameRel12; - std::vector mOriRelClin; + std::string mNameRel12; ///< relative orientation of Clino1/Clino2 stored as "i-kj" + std::vector mOriRelClin; ///< vector of relative orientations Clino[0]/Clino[K] + bool mShowAll; ///< Do we show all the msg relative to residuals }; @@ -172,9 +198,10 @@ cAppli_ClinoInit::cAppli_ClinoInit cMMVII_Appli (aVArgs,aSpec), mPhProj (*this), mNbStep0 (65), - mNbStep (10), + mNbStepIter (10), mNbIter (50), - mNameRel12 ("i-kj") + mNameRel12 ("i-kj"), + mShowAll (false) { } @@ -194,7 +221,7 @@ cCollecSpecArg2007 & cAppli_ClinoInit::ArgOpt(cCollecSpecArg2007 & anArgOpt) return anArgOpt << AOpt2007(mNbStep0,"NbStep0","Number of step at first iteration",{eTA2007::HDV}) - << AOpt2007(mNbStep,"NbStep","Number of step at current iteration",{eTA2007::HDV}) + << AOpt2007(mNbStepIter,"NbStep","Number of step at current iteration",{eTA2007::HDV}) << AOpt2007(mNbIter,"NbIter","Number of iteration",{eTA2007::HDV}) << AOpt2007(mASim,"AmplSim","Amplitude of rotation is simul [W,P,K,BS]",{{eTA2007::ISizeV,"[5,5]"}}) << AOpt2007(mNameRel12,"Rel12","orientation relative 2 to 1, if several clino",{eTA2007::HDV}) @@ -206,16 +233,22 @@ std::vector cAppli_ClinoInit::VectCostRot(const tRot & aRot0) const { std::vector aRes; + // Parse all clino currently selected for (size_t aIndC=0 ; aIndC cAppli_ClinoInit::IterInit(int aNbStep) +cWhichMin cAppli_ClinoInit::IterInit(int aNbStep) const { - cWhichMin aWMin(tRot::Identity(),1e10); + cWhichMin aWMin(tRot::Identity(),1e10); // initialize with a "big" residual + // structure for parsing all rotation using quaternion , with aNbStep division in each direction cSampleQuat aSQ(aNbStep,true); - for (size_t aKQ=0 ; aKQ aMat = Quat2MatrRot(aSQ.KthQuat(aKQ)); + cDenseMatrix aMat = Quat2MatrRot(aSQ.KthQuat(aKQ)); // quaternion -> rotation tRot aRot(aMat,false); - tREAL8 aCost = CostRot(aRot); // cost of tested rot + tREAL8 aCost = CostRot(aRot); // compute the cost of tested rot aWMin.Add(aRot,aCost); // update lowest cost solution } - StdOut() << "MIN0= " << std::sqrt(aWMin.ValExtre()) << std::endl; + if (mShowAll) + StdOut() << "MIN0= " << std::sqrt(aWMin.ValExtre()) << std::endl; return aWMin; } @@ -269,12 +304,13 @@ cWhichMin cAppli_ClinoInit::OneIter(const tRot & aR0,tREAL8 aStep, } } - StdOut() << "MIN " << std::sqrt(aWMin.ValExtre()) << " Step=" << aStep << std::endl; + if (mShowAll) + StdOut() << "MIN " << std::sqrt(aWMin.ValExtre()) << " Step=" << aStep << std::endl; return aWMin; } -void cAppli_ClinoInit::ComputeCond(const tRot & aRot0,tREAL8 aEpsilon) +cPt2dr cAppli_ClinoInit::ComputeCond(const tRot & aRot0,tREAL8 aEpsilon) { cStrStat2 aStat(3); // covariant matrix @@ -286,7 +322,7 @@ void cAppli_ClinoInit::ComputeCond(const tRot & aRot0,tREAL8 aEpsilon) tRot aRot = OriOfClino(aIndC,aRot0); for (const auto & aMes : mVMeasures) { - auto [aGx,aGy] = aMes.Grad(mVKClino[aIndC],aRot,aEpsilon); + auto [aGx,aGy] = aMes.GradEVR(mVKClino[aIndC],aRot,aEpsilon); aStat.Add(aGx.ToVect()); aStat.Add(aGy.ToVect()); } @@ -296,7 +332,9 @@ void cAppli_ClinoInit::ComputeCond(const tRot & aRot0,tREAL8 aEpsilon) // print ratio with highest value cDenseVect aDV = aStat.DoEigen().EigenValues(); - StdOut() << "R2=" << std::sqrt(aDV(0)/aDV(2)) << " R1=" << std::sqrt(aDV(1)/aDV(2)) << std::endl; + + return cPt2dr ( std::sqrt(aDV(0)/aDV(2)) , std::sqrt(aDV(1)/aDV(2)) ); + // StdOut() << "R2=" << std::sqrt(aDV(0)/aDV(2)) << " R1=" << std::sqrt(aDV(1)/aDV(2)) << std::endl; } int cAppli_ClinoInit::Exe() @@ -305,29 +343,19 @@ int cAppli_ClinoInit::Exe() mComputedClino = std::vector(mVKClino.size(),true); // initially all user clino are active - std::vector> aVNames; // for reading names of camera - std::vector> aVAngles; // for reading angles of clinometers - std::vector aVFakePts; // not used, required by ReadFilesStruct - - std::string mFormat = "NFFFF"; - // read angles and camera - ReadFilesStruct - ( - mNameClino, - mFormat, - 0,-1, - '#', - aVNames, - aVFakePts,aVFakePts, - aVAngles, - false - ); + // std::vector> aVNames; // for reading names of camera + // std::vector> aVAngles; // for reading angles of clinometers + // std::vector aVFakePts; // not used, required by ReadFilesStruct + + std::string mFormat = "ISFSF"; + cReadFilesStruct aRFS(mNameClino,mFormat,0,-1,'#'); + aRFS.Read(); tRot aRSim = tRot::Identity(); // Rotation for simulation mOriRelClin.push_back(tRot::Identity()); mOriRelClin.push_back(tRot::RotFromCanonicalAxes(mNameRel12)); - size_t aNbMeasures = aVNames.size(); + size_t aNbMeasures = aRFS.NbRead(); // if we do simulation, generate if (IsInit(&mASim)) { @@ -357,9 +385,9 @@ int cAppli_ClinoInit::Exe() } else { - std::string aNameIm = mPrePost[0] + aVNames[aKLine][0] + mPrePost[1]; + std::string aNameIm = mPrePost[0] + aRFS.VNameIm().at(aKLine) + mPrePost[1]; cSensorCamPC * aCam = mPhProj.ReadCamPC(aNameIm,true); - mVMeasures.push_back(cClinoCalMes1Cam(aCam,aVAngles[aKLine])); + mVMeasures.push_back(cClinoCalMes1Cam(aCam,aRFS.VNums().at(aKLine))); } } @@ -370,19 +398,23 @@ int cAppli_ClinoInit::Exe() for (int aK=0 ; aK double MoyAbs(cIm2D aImIn) return aSom[0] / aSom[1]; } - -template cImGrad Deriche(const cDataIm2D & aImIn,double aAlpha) +template void ComputeDeriche(cImGrad & aResGrad,const cDataIm2D & aImIn,double aAlpha) { auto aV1In = cMMV1_Conv::ImToMMV1(aImIn); - - cImGrad aResGrad(aImIn.Sz()); auto aV1Gx = cMMV1_Conv::ImToMMV1(aResGrad.mGx.DIm()); auto aV1Gy = cMMV1_Conv::ImToMMV1(aResGrad.mGy.DIm()); @@ -58,7 +55,12 @@ template cImGrad Deriche(const cDataIm2D & aImIn,double deriche(aV1In.in_proj(),aAlpha,10), Virgule(aV1Gx.out(),aV1Gy.out()) ); +} +template cImGrad Deriche(const cDataIm2D & aImIn,double aAlpha) +{ + cImGrad aResGrad(aImIn.Sz()); + ComputeDeriche(aResGrad,aImIn,aAlpha); return aResGrad; } diff --git a/MMVII/src/MMV1/Numerics.cpp b/MMVII/src/MMV1/Numerics.cpp index 29fd96bd94..072cfdf883 100755 --- a/MMVII/src/MMV1/Numerics.cpp +++ b/MMVII/src/MMV1/Numerics.cpp @@ -19,6 +19,13 @@ double NC_KthVal(std::vector & aV, double aProportion) return ::KthValProp(aV,aProportion); } +// template TVal KthVal(std::vector & aV,int aKth) +double IKthVal(std::vector & aV, int aK) +{ + return ::KthVal(aV,aK); +} + + double Cst_KthVal(const std::vector & aV, double aProportion) { diff --git a/MMVII/src/Sensors/cCentralPerspCam.cpp b/MMVII/src/Sensors/cCentralPerspCam.cpp index 8797fd0547..54855eacb4 100644 --- a/MMVII/src/Sensors/cCentralPerspCam.cpp +++ b/MMVII/src/Sensors/cCentralPerspCam.cpp @@ -478,9 +478,8 @@ double cPerspCamIntrCalib::DegreeVisibility(const cPt3dr & aP) const const std::vector & cPerspCamIntrCalib::DirBundles(tVecIn & aV3 ,const tVecOut & aV0 ) const { - UpdateLSQDistIfRequired(); - - CheckBeforeInverse(aV0); + UpdateLSQDistIfRequired(); // Updtate inverse by least square + CheckBeforeInverse(aV0); // Check that V0 is inversible, nothing done in release static tVecOut aV1,aV2; mMapIm2PProj.Values(aV1,aV0); @@ -500,15 +499,37 @@ cPt3dr cPerspCamIntrCalib::DirBundle(const tPtOut & aPt) const return aRes; } + +/* Undist : compute the bundle then reproject on a perfect stenope camera. + + For a stenope camera, generate un-necessary computation , P->Bundle->P/z which + turn to identity. But this the price to pay for being generik. + Maybe maybe do a specialize version later + + cPt2dr cPerspCamIntrCalib::Undist_Quick(const tPtOut & aP0) const; + + +PP,*F Dist-1 InProj x,y/z -PP,/F + PIm ------> PHgr --------> PUd -------> Bundle ----------> P1 ----------> PUndist +*/ + cPt2dr cPerspCamIntrCalib::Undist(const tPtOut & aP0) const { cPt3dr aPt = DirBundle(aP0); - cPt2dr aP1 = Proj(aPt) / aPt.z(); - return mMapPProj2Im.Value(aP1); } +cPt2dr cPerspCamIntrCalib::Redist(const tPtOut & aP0) const +{ + cPt2dr aP1 = mMapIm2PProj.Value(aP0); + cPt3dr aP2(aP1.x(),aP1.y(),1.0); + + return Value(aP2); +} + + + + tREAL8 cPerspCamIntrCalib::InvProjIsDef(const tPtOut & aPix ) const { return mDefProj->P2DIsDef(mDist_DirInvertible->Inverse(mMapIm2PProj.Value(aPix))); diff --git a/MMVII/src/Sensors/cPhotogrammetricProject.cpp b/MMVII/src/Sensors/cPhotogrammetricProject.cpp index e1a75825a7..f757bf3cd3 100644 --- a/MMVII/src/Sensors/cPhotogrammetricProject.cpp +++ b/MMVII/src/Sensors/cPhotogrammetricProject.cpp @@ -614,12 +614,17 @@ cSensorImage* cPhotogrammetricProject::ReadSensorFromFolder(const std::string & cPerspCamIntrCalib * cPhotogrammetricProject::InternalCalibFromImage(const std::string & aNameIm) const { - // 4 now, pretty basic allox sensor, extract internal, destroy - // later will have to handle : - // * case where calib exist but not pose - // * case where nor calib nor pose exist, and must be created from xif + // allox sensor and if exist, extract internal, destroy + // else try to extract calib from standard name + // * case where nor calib nor pose exist, and must be created from xif still to implemant mDPOrient.AssertDirInIsInit(); - cSensorCamPC * aPC = ReadCamPC(aNameIm,false); + + cSensorCamPC * aPC = ReadCamPC(aNameIm,false,SVP::Yes); + if (aPC==nullptr) + { + return InternalCalibFromStdName(aNameIm); + } + cPerspCamIntrCalib * aCalib = aPC->InternalCalib(); delete aPC; diff --git a/MMVII/src/Serial/ElemStrToVal.cpp b/MMVII/src/Serial/ElemStrToVal.cpp index f005800928..fbf2398abb 100644 --- a/MMVII/src/Serial/ElemStrToVal.cpp +++ b/MMVII/src/Serial/ElemStrToVal.cpp @@ -208,6 +208,11 @@ template<> cE2Str::tMapE2Str cE2Str::mE2S {eTA2007::FileImage,"Im"}, {eTA2007::FileCloud,"Cloud"}, {eTA2007::File3DRegion,"3DReg"}, + {eTA2007::FileTagged,"FileTagged"}, + {eTA2007::FileTxt,"FileText"}, + {eTA2007::FileAny,"FileAny"}, + {eTA2007::FolderAny,"FolderAny"}, + {eTA2007::MPatFile,"MPF"}, {eTA2007::Orient,"Ori"}, {eTA2007::RadiomData,"RadData"}, diff --git a/MMVII/src/Utils/uti_set_sel.cpp b/MMVII/src/Utils/uti_set_sel.cpp index 9049766669..d1cd1e9a2e 100755 --- a/MMVII/src/Utils/uti_set_sel.cpp +++ b/MMVII/src/Utils/uti_set_sel.cpp @@ -1,5 +1,7 @@ #include "MMVII_MMV1Compat.h" #include "MMVII_2Include_Serial_Tpl.h" +#include "MMVII_util_tpl.h" + #include #include @@ -417,6 +419,52 @@ template cSelector Str2Interv(const std::string & aStr) return aRes; } +/******************************************************************/ +/* */ +/* cIdTransformator */ +/* */ +/******************************************************************/ + +template Type cIdTransformator::Transfo(const Type & aVal) const +{ + return aVal; +} + +/******************************************************************/ +/* */ +/* cPatternTransfo */ +/* */ +/******************************************************************/ + +// replace a pattern : yy(.*)zz , A$1 , yytotozz => Atoto +//std::string ReplacePattern(const std::string & aPattern,const std::string & aSubst,const std::string & aString); + + +cPatternTransfo::cPatternTransfo(const std::string & aPat,const std::string & aSubst) : + mPat (aPat), + mSubst (aSubst) +{ +} + +cPatternTransfo::cPatternTransfo(const std::vector & aVec) : + cPatternTransfo + ( + GetDef(aVec,0,std::string(".*")), + GetDef(aVec,1,std::string("$0")) + ) +{ + if ( (aVec.size()!=0) && (aVec.size() != 2)) + MMVII_UnclasseUsEr("cPatternTransfo bad size"); +} + + +std::string cPatternTransfo::Transfo(const std::string & aValue) const +{ + return ReplacePattern(mPat,mSubst,aValue); +} + + + #define MACRO_INSTANTIATE_SELECTOR(Type)\ template cSelector GenIntervalSelector(const std::optional & aV1, const std::optional & aV2, bool aInclLow,bool InclUp);\ @@ -435,6 +483,7 @@ template class cDataNegSelector;\ template class cDataAndSelector;\ template class cDataOrSelector;\ template class cDataCsteSelector;\ +template class cIdTransformator;\ MACRO_INSTANTIATE_SELECTOR(int) @@ -724,17 +773,26 @@ template bool cExtSet::Equal(const cExtSet & aSet) con return IncludedIn(aSet) && aSet.IncludedIn(*this); } -template void cExtSet::Filter(const cSelector & aSel) +template void cExtSet::Filter(const cSelector & aSel,const cTransformator & aTranfo) { std::vector aV; PutInVect(aV,false); for (const auto & ptr : aV) { - if (! aSel.Match(*ptr)) + if (! aSel.Match(aTranfo.Transfo(*ptr))) Suppress(*ptr); } } + +template void cExtSet::Filter(const cSelector & aSel) +{ + cIdTransformator aTransfoId; + + Filter(aSel,aTransfoId); +} + + template void cExtSet::OpAff(const eOpAff & anOp,const cExtSet aS) { if (anOp==eOpAff::ePlusEq) (*this) += aS; diff --git a/src/tiff/tiff_args_opt.cpp b/src/tiff/tiff_args_opt.cpp index ab1e513d2c..c6cdac4cd0 100755 --- a/src/tiff/tiff_args_opt.cpp +++ b/src/tiff/tiff_args_opt.cpp @@ -105,7 +105,7 @@ D_Tiff_ifd_Arg_opt::D_Tiff_ifd_Arg_opt() : mSzFileTile = Pt2di(Tiff_Im::DefTileFile(),Tiff_Im::DefTileFile()); - std::cout << "mSzFileTilemSzFileTile " <