From 8fa3e337ddf976f669d7c9a35489ea399dbf2a69 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Sat, 26 Oct 2024 13:48:42 +0200 Subject: [PATCH 1/2] Pb sysco????? --- MMVII/include/MMVII_SysCo.h | 1 + 1 file changed, 1 insertion(+) diff --git a/MMVII/include/MMVII_SysCo.h b/MMVII/include/MMVII_SysCo.h index 4c1fea176b..edad6a240f 100644 --- a/MMVII/include/MMVII_SysCo.h +++ b/MMVII/include/MMVII_SysCo.h @@ -1,6 +1,7 @@ #ifndef _MMVII_SYSCO_H_ #define _MMVII_SYSCO_H_ + #include "MMVII_DeclareCste.h" #include "MMVII_Mappings.h" #include "MMVII_AllClassDeclare.h" From 00a5164cd93badfdcecbf815c23049e1c86dd561 Mon Sep 17 00:00:00 2001 From: deseilligny Date: Sun, 15 Dec 2024 22:07:42 +0100 Subject: [PATCH 2/2] Supress dependance to MMV1 in ImageFilterMMV1.cpp --- MMVII/include/MMVII_Bench.h | 2 + MMVII/include/MMVII_ImageInfoExtract.h | 2 +- MMVII/include/MMVII_Images.h | 1 + MMVII/include/MMVII_NonLinear2DFiltering.h | 10 +- MMVII/include/MMVII_nums.h | 6 + MMVII/include/V1VII.h | 2 +- MMVII/src/Bench/BenchImages.cpp | 1 + MMVII/src/ImagesBase/BaseImage.cpp | 14 + MMVII/src/ImagesFiltrLinear/Deriche.cpp | 240 ++++++++++++ MMVII/src/ImagesFiltrLinear/LocalFilters.cpp | 343 ++++++++++++++++++ .../src/ImagesInfoExtract/ExtracExtremum.cpp | 28 ++ MMVII/src/MMV1/ImageFilterMMV1.cpp | 322 ++++++++-------- MMVII/src/MorphoMat/ConnectedComponents.cpp | 11 + MMVII/src/UtiMaths/uti_rand.cpp | 5 + 14 files changed, 837 insertions(+), 150 deletions(-) create mode 100755 MMVII/src/ImagesFiltrLinear/Deriche.cpp create mode 100755 MMVII/src/ImagesFiltrLinear/LocalFilters.cpp diff --git a/MMVII/include/MMVII_Bench.h b/MMVII/include/MMVII_Bench.h index 03cf166ddb..41ba712bb6 100755 --- a/MMVII/include/MMVII_Bench.h +++ b/MMVII/include/MMVII_Bench.h @@ -162,8 +162,10 @@ void BenchBaseImage(); ///< Global bench on image void BenchImNDim(); void BenchIm3D(); ///< Bench on fulll 3D Images + "Layer" images + void BenchGlobImage2d(); ///< Global bench on image void BenchFileImage(); ///< Global bench on image +void BenchImFilterV1V2(); ///< Temporary, check that new implementation are close enough to V1's void TestTimeV1V2(); ///< Not a formal Bench, require visual inspection diff --git a/MMVII/include/MMVII_ImageInfoExtract.h b/MMVII/include/MMVII_ImageInfoExtract.h index 35bed3ecd1..f62d5fbd1e 100755 --- a/MMVII/include/MMVII_ImageInfoExtract.h +++ b/MMVII/include/MMVII_ImageInfoExtract.h @@ -31,7 +31,7 @@ double ValueStableInterpoleExtr(double V1,double V2,double V3); template double MoyAbs(cIm2D aImIn); ///< Compute average of Abs of Image -template cPt2dr ValExtre(cIm2D aImIn); ///< X -> Min, Y -> Max +template std::pair ValExtre(cIm2D aImIn); ///< X -> Min, Y -> Max /// Class to store results of extremum struct cResultExtremum diff --git a/MMVII/include/MMVII_Images.h b/MMVII/include/MMVII_Images.h index e736ff8fdd..c613015fb1 100755 --- a/MMVII/include/MMVII_Images.h +++ b/MMVII/include/MMVII_Images.h @@ -379,6 +379,7 @@ template class cDataTypedIm : public cDataGenUnTypedI double L2Dist(const cDataTypedIm & aV,bool Avg=true) const; ///< Dist som square double SqL2Dist(const cDataTypedIm & aV,bool Avg=true) const; ///< Square L2Dist double LInfDist(const cDataTypedIm & aV) const; ///< Dist max + double SafeMaxRelDif(const cDataTypedIm & aV,tREAL8 aEps) const; ///< Maximal safe relative dif double L1Norm(bool Avg=true) const; ///< Norm som abs double L2Norm(bool Avg=true) const; ///< Norm square double LInfNorm() const; ///< Nomr max diff --git a/MMVII/include/MMVII_NonLinear2DFiltering.h b/MMVII/include/MMVII_NonLinear2DFiltering.h index 4b735d08b3..b054345cd6 100755 --- a/MMVII/include/MMVII_NonLinear2DFiltering.h +++ b/MMVII/include/MMVII_NonLinear2DFiltering.h @@ -27,10 +27,16 @@ template void SelfLabMaj(cIm2D aImIn,const cBox2di &); /// Basic laplacien -cIm2D Lapl(cIm2D aImIn); // Well linear ... +template cIm2D Lapl(cIm2D aImIn); // Well linear ... + +/// return the label majoritar in a neighbourhood +template cIm2D LabMaj(cIm2D aImIn,const cBox2di &); + + /// Extincion function = dist of image (V!=0) to its complementar (V==0) -void MakeImageDist(cIm2D aImIn,const std::string & aNameChamfer="32"); +// No longer implemented, used MMV1 +// void MakeImageDist(cIm2D aImIn,const std::string & aNameChamfer="32"); diff --git a/MMVII/include/MMVII_nums.h b/MMVII/include/MMVII_nums.h index e68148faa6..e641345193 100755 --- a/MMVII/include/MMVII_nums.h +++ b/MMVII/include/MMVII_nums.h @@ -73,6 +73,8 @@ double RandInInterval(double a,double b); ///< Uniform distribution in [a,b] double RandInInterval(const cPt2dr &interval); ///< Uniform distribution in [interval.x,interval.y] double RandInInterval_C(const cPt2dr &interval); ///< Uniform distribution in [-interval.y,-interval.x]U[interval.x,interval.y] +int RandUnif_M_N(int aM,int aN); ///< Uniform disrtibution in [M,N] + /** Class for mapping object R->R */ class cFctrRR { @@ -324,6 +326,8 @@ template <> class tElemNumTrait : public tBaseNumTrait { public : static tINT2 DummyVal() {MMVII_INTERNAL_ERROR("No DummyVal for type");return 0;} + static tINT2 MaxVal() {return std::numeric_limits::max();} + static tINT2 MinVal() {return std::numeric_limits::min();} static bool Signed() {return true;} static eTyNums TyNum() {return eTyNums::eTN_INT2;} typedef tREAL4 tFloatAssoc; @@ -353,6 +357,8 @@ template <> class tElemNumTrait : public tBaseNumTrait template <> class tElemNumTrait : public tBaseNumTrait { public : + static tREAL4 MaxVal() {return std::numeric_limits::max();} + static tREAL4 MinVal() {return std::numeric_limits::min();} static tREAL4 DummyVal() {return std::nanf("");} static tREAL4 Accuracy() {return 1e-2f;} static bool Signed() {return true;} ///< Not usefull but have same interface diff --git a/MMVII/include/V1VII.h b/MMVII/include/V1VII.h index 42a48fe22a..5ed50ff13b 100755 --- a/MMVII/include/V1VII.h +++ b/MMVII/include/V1VII.h @@ -33,7 +33,7 @@ template std::string ToStrComMMV1(const cTplBox & aBox) void ExportHomMMV1(const std::string & aIm1,const std::string & aIm2,const std::string & SH,const std::vector & aVP); void ExportHomMMV1(const std::string & aIm1,const std::string & aIm2,const std::string & SH,const std::vector & aVP); -void MakeStdIm8BIts(cIm2D aImIn,const std::string& aName); +void MakeStdIm8BIts(const cIm2D &aImIn,const std::string& aName); diff --git a/MMVII/src/Bench/BenchImages.cpp b/MMVII/src/Bench/BenchImages.cpp index 55a3eb264e..aee124b2c4 100755 --- a/MMVII/src/Bench/BenchImages.cpp +++ b/MMVII/src/Bench/BenchImages.cpp @@ -585,6 +585,7 @@ void BenchGlobImage(cParamExeBench & aParam) BenchRectObj(); BenchBaseImage(); BenchGlobImage2d(); + BenchImFilterV1V2(); aParam.EndBench(); } diff --git a/MMVII/src/ImagesBase/BaseImage.cpp b/MMVII/src/ImagesBase/BaseImage.cpp index 52d7d442d9..a28fe1203c 100755 --- a/MMVII/src/ImagesBase/BaseImage.cpp +++ b/MMVII/src/ImagesBase/BaseImage.cpp @@ -126,6 +126,20 @@ template } +template + double cDataTypedIm::SafeMaxRelDif(const cDataTypedIm & aI2,tREAL8 aEps) const +{ + tPB::AssertSameArea(aI2); + double aRes = 0; + for (int aK=0 ; aK double cDataTypedIm::LInfDist(const cDataTypedIm & aI2) const diff --git a/MMVII/src/ImagesFiltrLinear/Deriche.cpp b/MMVII/src/ImagesFiltrLinear/Deriche.cpp new file mode 100755 index 0000000000..fa4962181b --- /dev/null +++ b/MMVII/src/ImagesFiltrLinear/Deriche.cpp @@ -0,0 +1,240 @@ +#include "MMVII_Linear2DFiltering.h" + +//=================================================== +//=================================================== +//=================================================== + +namespace MMVII +{ + + +class Canny_Deriche_Param +{ + public : + Canny_Deriche_Param(tREAL8 alpha); + tREAL8 a,a1,a2,a3,a4; + tREAL8 b1,b2; + + tREAL8 mAmpl; +}; + +Canny_Deriche_Param::Canny_Deriche_Param(tREAL8 alpha) +{ + tREAL8 exp_alpha = (tREAL8) exp( (double) - alpha); + tREAL8 exp_2alpha = exp_alpha * exp_alpha; + + tREAL8 kd = - (1-exp_alpha)*(1-exp_alpha)/exp_alpha; + tREAL8 ks = (1-exp_alpha)*(1-exp_alpha)/(1+ 2*alpha*exp_alpha - exp_2alpha); + + a = kd * exp_alpha; + b1 = 2 * exp_alpha; + b2 = exp_2alpha; + + a1 = ks; + a2= ks * exp_alpha * (alpha-1) ; + a3 = ks * exp_alpha * (alpha+1); + a4 = ks * exp_2alpha; + + + mAmpl = (2*(b1-a-2*b2))/(1-b1+b2); +} + +/* Merci a Canny, Rachid Deriche et Tuan Dang */ + +template + void deriche_ll // Deriche Low Level interface + ( + TypeOut * IGX, + TypeOut * IGY, + TypeIn * IE, + tINT8 DX, + tINT8 DY, + tREAL8 alpha + ) +{ + // TypeOut *YP, *YM; + //int i, j; + + cPt2di aSz(DX,DY); + + cIm2D aImYP(aSz); + cIm2D aImYM(aSz); + + TypeOut * YP = aImYP.DIm().RawDataLin(); + TypeOut * YM = aImYM.DIm().RawDataLin(); + + + const Canny_Deriche_Param p (alpha); + + /***********************************************/ + /***** Derivation selon LIGNE : on fixe y *****/ + /***********************************************/ + /**** Filtrage de Gauche a Droite ****/ + for (int i=0; i< DY; i++) { + YM[i*DX+0] = 0; + YM[i*DX+1] = 0; + for (int j=2; j< DX; j++) { + YM[i*DX+j] = p.a * (tREAL8)IE[i*DX+j-1] + + p.b1 * YM[i*DX + j-1] + - p.b2 * YM[i*DX + j-2] ; + } + } + /**** Filtrage de Droite a Gauche ****/ + for (int i=0; i< DY; i++) { + YP[i*DX+DX-1] = 0; + YP[i*DX+DX-2] = 0; + for (int j= DX-3; j >= 0; j--) { + YP[i*DX+j] = -p.a *(tREAL8)IE[i*DX+j+1] + + p.b1 * YP[i*DX + j+1] + - p.b2 * YP[i*DX + j+2] ; + } + } + /**** Derivee selon LIGNE ****/ + for (int i=0; i< DY; i++) { + for (int j=0; j< DX; j++) { + IGX[i*DX + j] = YP[i*DX + j] + YM[i*DX + j] ; + + } + } + /*******************************************************************/ + /***** Lissage suivant COLONNE --> composante en x du gradient *****/ + /*******************************************************************/ + /**** Filtrage de Gauche a Droite ****/ + for (int j=0; j< DX; j++) { + YM[0*DX+j]= 0; + YM[1*DX+j]= 0; + for (int i= 2; i< DY; i++) { + YM[i*DX+j] = p.a1 * IGX[i*DX+j] + p.a2 * IGX[(i-1)*DX + j] + + p.b1 * YM[(i-1)*DX + j] - p.b2 * YM[(i-2)*DX + j]; + } + } + /**** Filtrage de Droite a Gauche ****/ + for (int j=0; j< DX; j++) { + YP[(DY-1)*DX+j] =0; + YP[(DY-2)*DX+j] =0; + for (int i= DY-3; i>= 0; i--) { + YP[i*DX+j] = p.a3 * IGX[(i+1)*DX+j] - p.a4 * IGX[(i+2)*DX+j] + + p.b1 * YP[(i+1)*DX+j] - p.b2 * YP[(i+2)*DX+j]; + } + } + /**** Lissage suivant COLONNE ****/ + for (int i= 0; i< DY; i++) { + for (int j= 0; j< DX; j++) { + IGX[i*DX + j] = YP[i*DX + j] + YM[i*DX + j]; + } + } + + /************************************************/ + /***** Derivation selon COLONNE : on fixe x *****/ + /************************************************/ + /**** Filtrage de Gauche a Droite ****/ + for (int j=0; j< DX; j++) { + YM[0*DX+j] = 0; + YM[1*DX+j] = 0; + for (int i=2; i< DY; i++) { + YM[i*DX+j] = p.a * (tREAL8)IE[(i-1)*DX+j] + + p.b1 * YM[(i-1)*DX +j] + - p.b2 * YM[(i-2)*DX + j] ; + } + } + /**** Filtrage de Droite a Gauche ****/ + for (int j=0; j< DX; j++) { + YP[(DY-1)*DX+j] = 0; + YP[(DY-2)*DX+j] = 0; + for (int i= DY-3; i >= 0; i--) { + YP[i*DX+j] = -p.a * (tREAL8)IE[(i+1)*DX+j] + + p.b1 * YP[(i+1)*DX + j] + - p.b2 * YP[(i+2)*DX+j] ; + } + } + /**** Derivee selon COLONNE ****/ + for (int i=0; i< DY; i++) { + for (int j=0; j< DX; j++) { + IGY[i*DX + j] = YP[i*DX + j] + YM[i*DX + j]; + } + } + /******************************************************************/ + /* Lissage suivant LIGNE --> composante en y du gradient *****/ + /******************************************************************/ + /**** Filtrage de Gauche a Droite ****/ + + for (int i=0; i< DY; i++) { + YM[i*DX+0]= 0; + YM[i*DX+1]= 0; + for (int j= 2; j< DX; j++) { + YM[i*DX+j] = p.a1 * IGY[i*DX+j] + p.a2 * IGY[i*DX + j-1] + + p.b1 * YM[i*DX + j-1] - p.b2 * YM[i*DX + j-2]; + } + } + + /**** Filtrage de Droite a Gauche ****/ + for (int i=0; i< DY; i++) { + YP[i*DX+ DX-1] =0; + YP[i*DX+ DX-2] =0; + for (int j= DX-3; j>= 0; j--) { + YP[i*DX+j] = p.a3 * IGY[i*DX+j+1] - p.a4 * IGY[i*DX+j+2] + + p.b1 * YP[i*DX+j+1] - p.b2 * YP[i*DX+j+2]; + } + } + /**** Lissage suivant LIGNE ****/ + for (int i= 0; i< DY; i++) { + for (int j= 0; j< DX; j++) { + IGY[i*DX + j] = YP[i*DX + j] + YM[i*DX + j]; + } + } + // Not in initial code , but in mmv1 + + for (int i= 0; i< DY; i++) { + for (int j= 0; j< DX; j++) { + IGX[i*DX + j] /= p.mAmpl; + IGY[i*DX + j] /= p.mAmpl; + } + } + +} + +template void ComputeDeriche(cImGrad & aResGrad,const cDataIm2D & aImIn,double aAlpha) +{ + aImIn.AssertSameArea(aResGrad.mGx.DIm()); + aImIn.AssertSameArea(aResGrad.mGy.DIm()); + + deriche_ll + ( + aResGrad.mGx.DIm().RawDataLin(), + aResGrad.mGy.DIm().RawDataLin(), + aImIn.RawDataLin(), + aImIn.Sz().x(), + aImIn.Sz().y(), + aAlpha + ); + +/* + auto aV1In = cMMV1_Conv::ImToMMV1(aImIn); + auto aV1Gx = cMMV1_Conv::ImToMMV1(aResGrad.mGx.DIm()); + auto aV1Gy = cMMV1_Conv::ImToMMV1(aResGrad.mGy.DIm()); + + ELISE_COPY + ( + aV1In.all_pts(), + 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; +} + + +template void deriche_ll(tREAL4 *IGX,tREAL4 * IGY,tREAL4 * IE,tINT8 DX,tINT8 DY,tREAL8 alpha); +template cImGrad Deriche(const cDataIm2D & aImIn,double aAlpha); +template void ComputeDeriche(cImGrad & aResGrad,const cDataIm2D & aImIn,double aAlpha); + +}; + + diff --git a/MMVII/src/ImagesFiltrLinear/LocalFilters.cpp b/MMVII/src/ImagesFiltrLinear/LocalFilters.cpp new file mode 100755 index 0000000000..f876617130 --- /dev/null +++ b/MMVII/src/ImagesFiltrLinear/LocalFilters.cpp @@ -0,0 +1,343 @@ +#include "MMVII_Matrix.h" +#include "MMVII_Linear2DFiltering.h" +#include "MMVII_ImageInfoExtract.h" + + +/** \file : contain implemenation for filter that are : + - invariant by translation + - defined on bounded neighboorhoud + +*/ + + + +namespace MMVII +{ + +void MakeStdIm8BIts(const cIm2D& aImIn,const std::string& aName) +{ + cComputeStdDev aStdDev; + const cDataIm2D & aDIm = aImIn.DIm(); + + for (const auto & aPt : aDIm) + aStdDev.Add(aDIm.GetV(aPt)); + + aStdDev = aStdDev.Normalize(); + + cIm2D aImINorm(aDIm.Sz()); + for (const auto & aPt : aDIm) + { + tREAL8 aNormV = aStdDev.NormalizedVal(aDIm.GetV(aPt)); + aNormV = 128 + 128 * std::erf(aNormV); + aImINorm.DIm().SetVTrunc(aPt,round_ni(aNormV)); + } + aImINorm.DIm().ToFile(aName); +} + + + +static const cPt2di TheP00(0,0); +static const cPt2di ThePXp(1,0); +static const cPt2di ThePXm(-1,0); +static const cPt2di ThePYp(0,1); +static const cPt2di ThePYm(0,-1); + +static const cPt2di ThePXpYp( 1, 1); +static const cPt2di ThePXmYm(-1,-1); +static const cPt2di ThePXpYm( 1,-1); +static const cPt2di ThePXmYp(-1, 1); + +/* ********************************************************** * +* * +* Functers class for image access * +* * +************************************************************** */ + +/* To limit code duplication , the tem +*/ + +/// Idem cImAccesStd but add a central point P0 +template class cImAccesStd_CenteredPt +{ + public : + cImAccesStd_CenteredPt(const cDataIm2D & aIm,const cPt2di & aP0) : mIm (aIm), mP0(aP0) {} + inline Type operator() (const cPt2di & aP) const {return mIm.GetV(aP+mP0);} + typedef Type tVal; + typedef typename tBaseNumTrait::tBase tBase; + + private : + const cDataIm2D& mIm; + cPt2di mP0; +}; + + + +/// Idem cImAccesProj but add a central point P0 +template class cImAccesProj_CenteredPt +{ + public : + cImAccesProj_CenteredPt(const cDataIm2D & aIm,const cPt2di& aP0) : mIm (aIm) , mP0 (aP0) {} + inline Type operator() (const cPt2di & aP) const {return mIm.GetV(mIm.Proj(aP+mP0));} + typedef Type tVal; + typedef typename tBaseNumTrait::tBase tBase; + private : + const cDataIm2D& mIm; + cPt2di mP0; +}; + +template + cIm2D GenLocal_Filter + ( + const cIm2D& anImIn, + const tFctrInside & aFctrI, + const tFctrBorder & aFctBrd, + int aSzBrd + ) +{ + const cDataIm2D & aDImIn(anImIn.DIm()); + cIm2D aRes(aDImIn.Sz()); + + for (const auto & aPt : aDImIn.Interior(aSzBrd)) + { + typename tFctrInside::tFuncter aImAccessStd (aDImIn,aPt); + aRes.DIm().SetV(aPt,aFctrI(aImAccessStd)); + } + + if (aSzBrd) + { + for (const auto & aPt : aDImIn.Border(aSzBrd) ) + { + typename tFctrBorder::tFuncter aImAccessBrd (aDImIn,aPt); + aRes.DIm().SetV(aPt,aFctBrd(aImAccessBrd)); + } + } + + return aRes; +} + + +/* Not used for now +/// Make a functor for an image, giving standard access to its value +template class cImAccesStd +{ + public : + cImAccesStd(const cDataIm2D & aIm) : mIm (aIm) {} + inline Type operator() (const cPt2di & aP) const {return mIm.GetV(aP);} + typedef Type tVal; + typedef typename tBaseNumTrait::tBase tBase; + + private : + const cDataIm2D& mIm; +}; +/// Make a "safe" functor for an image, giving access to its value of projected point +template class cImAccesProj +{ + public : + cImAccesProj(const cDataIm2D & aIm) : mIm (aIm) {} + inline Type operator() (const cPt2di & aP) const {return mIm.GetV(mIm.Proj(aP));} + typedef Type tVal; + typedef typename tBaseNumTrait::tBase tBase; + private : + const cDataIm2D& mIm; +}; +*/ + + +/* ****************************************************** */ +/* */ +/* CourbTgt */ +/* */ +/* ****************************************************** */ + +/* +*/ + +/** Fct for computing the CourbTgt filter */ + +template class cCourbTgtOfFctr +{ + public : + typedef typename tFunc::tBase tBase; + typedef typename tFunc::tVal tVal ; + typedef tFunc tFuncter; + + cCourbTgtOfFctr(tREAL8 aExp) : mExp(aExp) {} + + tREAL8 operator() (const tFunc & aFunc) const + { + // compute values that will be used several times + tVal a2V00 = aFunc(TheP00) * 2; + tVal aVxP1 = aFunc(ThePXp); + tVal aVxM1 = aFunc(ThePXm); + tVal aVyP1 = aFunc(ThePYp); + tVal aVyM1 = aFunc(ThePYm); + + // compute gradient and its norm + tBase aGx = (aVxP1 - aVxM1) / 2; + tBase aGy = (aVyP1 - aVyM1) / 2; + tREAL8 aG2 = Square(aGx)+Square(aGy); + + if (aG2 == 0) + { + return 0; + } + + // compute hessian + tVal aD2xx = (aVxP1+aVxM1- a2V00); + tVal aD2yy = (aVyP1+aVyM1- a2V00); + tVal aD2xy = ( aFunc(ThePXpYp) + aFunc(ThePXmYm) - aFunc(ThePXmYp)-aFunc(ThePXpYm)) / 4; + + // apply the hessian to direction orthonal to gradient (= direction of the level curve) + return ( (aD2xx * aGy) * aGy - 2 * (aD2xy * aGx) * aGy + (aD2yy * aGx) * aGx) / std::pow(aG2,mExp); + } + private : + tREAL8 mExp; +}; + + + +template cIm2D CourbTgt(cIm2D aImIn) +{ + return GenLocal_Filter + ( + aImIn, + cCourbTgtOfFctr>(0.5), + cCourbTgtOfFctr>(0.5), + 1 + ); +} +template void SelfCourbTgt(cIm2D anIm) +{ + cIm2D aRes = CourbTgt(anIm); + aRes.DIm().DupIn(anIm.DIm()); +} + +/* ****************************************************** */ +/* */ +/* Lapl */ +/* */ +/* ****************************************************** */ + +template class cLaplOfFctr +{ + public : + typedef typename tFunc::tBase tBase; + typedef typename tFunc::tVal tVal ; + typedef tFunc tFuncter; + + tREAL8 operator() (const tFunc & aFunc) const + { + return 4 * aFunc(TheP00) - (aFunc(ThePXp) + aFunc(ThePXm) + aFunc(ThePYp) + aFunc(ThePYm)); + } +}; + +template cIm2D Lapl(cIm2D aImIn) +{ + return GenLocal_Filter + ( + aImIn, + cLaplOfFctr>(), + cLaplOfFctr>(), + 1 + ); +} + +/* ****************************************************** */ +/* */ +/* MajLab */ +/* */ +/* ****************************************************** */ + +template class cMajLabOfFctr +{ + public : + typedef typename tFunc::tBase tBase; + typedef typename tFunc::tVal tVal ; + typedef tFunc tFuncter; + + cMajLabOfFctr(const cBox2di & aBox,int aVMin,int aVMax) : + mBox (aBox), + mVMin (aVMin), + mVMax (aVMax), + mHisto (1+mVMax-mVMin) + { + } + + tREAL8 operator() (const tFunc & aFunc) const + { + // reset histogramm to 0 + for (auto & aH : mHisto) + aH=0; + // increment histogramm for each label met + for (const auto & aP : cRect2(mBox)) // cRect2::BoxWindow(mSz)) + mHisto.at(aFunc(aP)-mVMin) ++; + + cWhichMax aWMax; + for (size_t aK=0 ; aK mHisto; +}; + +template cIm2D LabMaj(cIm2D aImIn,const cBox2di & aBox) +{ + auto [aMin,aMax] = ValExtre(aImIn); + + return GenLocal_Filter + ( + aImIn, + cMajLabOfFctr>(aBox,aMin,aMax), + cMajLabOfFctr>(aBox,aMin,aMax), + std::max(NormInf(aBox.P0()),NormInf(aBox.P1())) + ); +} +template void SelfLabMaj(cIm2D anImIn,const cBox2di & aBox) +{ + cIm2D aRes = LabMaj(anImIn,aBox); + aRes.DIm().DupIn(anImIn.DIm()); +} + + +template cIm2D LabMaj(cIm2D aImIn,const cBox2di &); +template void SelfLabMaj(cIm2D anImIn,const cBox2di & aBox); + + +/* +*/ + + //================= Global Solution ================== + + +/* Apparently no longer used, maintain it it if some regret .... +cIm2D Lapl(cIm2D aImIn) +{ + cIm2D aRes(aImIn.DIm().Sz()); + + Im2D aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); + Im2D aV1Res = cMMV1_Conv::ImToMMV1( aRes.DIm()); + + ELISE-COPY(aV1In.all_pts(),Laplacien(aV1In.in_proj()),aV1Res.out()); + + return aRes; +} +*/ + + + +#define INSTANTIATE_TRAIT_AIME(TYPE)\ +template cIm2D CourbTgt(cIm2D aImIn);\ +template cIm2D Lapl(cIm2D aImIn);\ +template void SelfCourbTgt(cIm2D aImIn); + +INSTANTIATE_TRAIT_AIME(tREAL4) +INSTANTIATE_TRAIT_AIME(tINT2) + + + +}; diff --git a/MMVII/src/ImagesInfoExtract/ExtracExtremum.cpp b/MMVII/src/ImagesInfoExtract/ExtracExtremum.cpp index 283024c5c4..89b59375c2 100755 --- a/MMVII/src/ImagesInfoExtract/ExtracExtremum.cpp +++ b/MMVII/src/ImagesInfoExtract/ExtracExtremum.cpp @@ -802,5 +802,33 @@ void BenchExtre(cParamExeBench & aParam) aParam.EndBench(); } +template std::pair ValExtre(cIm2D aImIn) +{ + Type aVMin = tElemNumTrait::MaxVal(); + Type aVMax = tElemNumTrait::MinVal(); + const cDataIm2D & aDIm = aImIn.DIm(); + + for (const auto & aPt : aDIm) + UpdateMinMax(aVMin,aVMax,aDIm.GetV(aPt)); + + return std::pair(aVMin,aVMax); +} +template double MoyAbs(cIm2D aImIn) +{ + cWeightAv aAvg; + const cDataIm2D & aDIm = aImIn.DIm(); + for (const auto & aPt : aDIm) + aAvg.Add(1.0,std::fabs(aDIm.GetV(aPt))); + + return aAvg.Average(); +} + +template double MoyAbs(cIm2D aImIn); +template double MoyAbs(cIm2D aImIn); + +template std::pair ValExtre(cIm2D aImIn); +template std::pair ValExtre(cIm2D aImIn); + + }; diff --git a/MMVII/src/MMV1/ImageFilterMMV1.cpp b/MMVII/src/MMV1/ImageFilterMMV1.cpp index b40ceb1b6a..94302bf1a7 100755 --- a/MMVII/src/MMV1/ImageFilterMMV1.cpp +++ b/MMVII/src/MMV1/ImageFilterMMV1.cpp @@ -1,6 +1,12 @@ +#define WITH_MMV1_BENCH false + +#if (WITH_MMV1_BENCH) #include "V1VII.h" +#endif + #include "MMVII_Matrix.h" #include "MMVII_Linear2DFiltering.h" +#include "MMVII_NonLinear2DFiltering.h" // FIXME CM->MPD: Mail a MPD ! Qu'est ce qu'on fait ici pour remplacer ELISE_COPY ? @@ -13,93 +19,89 @@ exactly what I want, it's much faster to implement them with MMV1. */ -// Test git namespace MMVII { -/* - * - REAL** i0 = in[0]; // Image - REAL** p0 = in[1]; // power - REAL* o0 = out[0]; - for (INT x=arg.x0() ; x aILab(aSz); + cDataIm2D& aDILab = aILab.DIm(); - */ + for (const auto & aPt : aDILab) + aDILab.SetV(aPt,RandUnif_M_N(aLabMin,aLabMax)); -/// Make a functor for an image, giving standard access to its value -template class cImAccesStd -{ - public : - cImAccesStd(const cDataIm2D & aIm) : mIm (aIm) {} - Type operator() (const cPt2di & aP) {return mIm.GetV(aP);} - typedef Type tVal; - private : - const cDataIm2D& mIm; -}; -/// Make a "safe" functor for an image, giving standard access to its value -template class cImAccesProj -{ - public : - cImAccesProj(const cDataIm2D & aIm) : mIm (aIm) {} - Type operator() (const cPt2di & aP) {return mIm.GetV(mIm.Proj(aP));} - typedef Type tVal; - private : - const cDataIm2D& mIm; -}; + cIm2D aV2Maj = LabMaj(aILab,aBox); + cDataIm2D& aDV2Maj = aV2Maj.DIm(); + for (const auto & aPt : aDV2Maj) + { + std::vector aHisto(aLabMax-aLabMin+1,0); + for (const auto & aV : cRect2(aBox)) + { + aHisto.at(aDILab.GetV(aDILab.Proj(aPt+aV))-aLabMin)++; + } + int aNb = aHisto.at(aDV2Maj.GetV(aPt)-aLabMin); + for (const auto & aElem : aHisto) + MMVII_INTERNAL_ASSERT_bench( aElem<=aNb,"Bench_LabMaj"); + } -/* -template typename tFunc::tVal NewCourbTgt(const tFunc & aFunc,const cPt2di & aPt) +} +void Bench_LabMaj() { + for (int aK=0 ; aK<100 ; aK++) + { + cPt2di aSz(RandUnif_M_N(20,30),RandUnif_M_N(20,30)); + int aVMin = RandUnif_M_N(-5,2); + int aVMax = aVMin+1+RandUnif_N(7); + cPt2di aP0(-RandUnif_N(3),-RandUnif_N(3)); + cPt2di aP1 = aP0 + cPt2di(1+RandUnif_N(6),1+RandUnif_N(6)); + Bench_LabMaj(aSz,cBox2di(aP0,aP1),aVMin,aVMax); + } } -*/ +#if (WITH_MMV1_BENCH) -template cIm2D MMVII_CourbTgt(const cIm2D & aImIn) +// Not implemanted/used in MMVII 4 now +void MMV1_MakeImageDist(cIm2D aImIn,const std::string & aNameChamfer) { - const cDataIm2D & aDImIn(aImIn.DIm()); - - cImAccesStd aIAStd (aDImIn); - cImAccesProj aIAProj(aDImIn); + const Chamfer & aChamf = Chamfer::ChamferFromName(aNameChamfer); + aChamf.im_dist(cMMV1_Conv::ImToMMV1(aImIn.DIm())) ; +} +void ExportHomMMV1(const std::string & aIm1,const std::string & aIm2,const std::string & SH,const std::vector & aVP) +{ + ElPackHomologue aPackH; - cIm2D aRes(aDImIn.Sz()); - for (const auto & aPt : aDImIn.Interior(1) ) + for (int aK=0 ; aKAssoc1To2(aKeyH,aIm1,aIm2,true); - return aRes; -} + aPackH.StdPutInFile(aNameF); -template cIm2D CourbTgt(cIm2D aImIn) +} +void ExportHomMMV1(const std::string & aIm1,const std::string & aIm2,const std::string & SH,const std::vector & aVI) { -MMVII_CourbTgt(aImIn); + std::vector aVR; + for (const auto & aPI: aVI) + aVR.push_back(ToR(aPI)); + ExportHomMMV1(aIm1,aIm2,SH,aVR); +} + + //================= V1 Solution ================== + +template cIm2D MMV1_CourbTgt(cIm2D aImIn) +{ cIm2D aRes(aImIn.DIm().Sz()); auto aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); @@ -109,15 +111,48 @@ MMVII_CourbTgt(aImIn); return aRes; } +template void MMV1_SelfCourbTgt(cIm2D aImIn) +{ + auto aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); + ELISE_COPY(aV1In.all_pts(),courb_tgt(aV1In.in_proj(),0.5),aV1In.out()); +} + +cIm2D MMV1_Lapl(cIm2D aImIn) +{ + cIm2D aRes(aImIn.DIm().Sz()); + + Im2D aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); + Im2D aV1Res = cMMV1_Conv::ImToMMV1( aRes.DIm()); + + ELISE_COPY(aV1In.all_pts(),Laplacien(aV1In.in_proj()),aV1Res.out()); + + return aRes; +} -template void SelfCourbTgt(cIm2D aImIn) + //================= Global Solution ================== + + +template void BenchImFilterV1V2(cIm2D anI1,cIm2D anI2,tREAL8 aEps) +{ + + tREAL8 aRD = anI1.DIm().SafeMaxRelDif(anI2.DIm(),1e-2); + if ( aRD > aEps ) + { + MMVII_INTERNAL_ASSERT_bench( false,"BenchImFilterV1V2 RD="+ToStr(aRD)); + } + +} + +template std::pair MMV1_ValExtre(cIm2D aImIn) { auto aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); - ELISE_COPY(aV1In.all_pts(),courb_tgt(aV1In.in_proj(),0.5),aV1In.out()); + double aVMin,aVMax; + ELISE_COPY(aV1In.all_pts(),aV1In.in(),VMin(aVMin)|VMax(aVMax)); + return std::pair((Type)aVMin,(Type)aVMax); } -template double MoyAbs(cIm2D aImIn) +template double MMV1_MoyAbs(cIm2D aImIn) { auto aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); double aSom[2]; @@ -125,7 +160,7 @@ template double MoyAbs(cIm2D aImIn) return aSom[0] / aSom[1]; } -template void ComputeDeriche(cImGrad & aResGrad,const cDataIm2D & aImIn,double aAlpha) +template void MMV1_ComputeDeriche(cImGrad & aResGrad,const cDataIm2D & aImIn,double aAlpha) { auto aV1In = cMMV1_Conv::ImToMMV1(aImIn); auto aV1Gx = cMMV1_Conv::ImToMMV1(aResGrad.mGx.DIm()); @@ -138,38 +173,14 @@ template void ComputeDeriche(cImGrad & aResGrad,const cDataIm2 Virgule(aV1Gx.out(),aV1Gy.out()) ); } - -template cImGrad Deriche(const cDataIm2D & aImIn,double aAlpha) +template cImGrad MMV1_Deriche(const cDataIm2D & aImIn,double aAlpha) { cImGrad aResGrad(aImIn.Sz()); - ComputeDeriche(aResGrad,aImIn,aAlpha); + MMV1_ComputeDeriche(aResGrad,aImIn,aAlpha); return aResGrad; } - -#define INSTANTIATE_TRAIT_AIME(TYPE)\ -template cImGrad Deriche(const cDataIm2D &aImIn,double aAlpha);\ -template double MoyAbs(cIm2D aImIn);\ -template cIm2D CourbTgt(cIm2D aImIn);\ -template void SelfCourbTgt(cIm2D aImIn); - -INSTANTIATE_TRAIT_AIME(tREAL4) -INSTANTIATE_TRAIT_AIME(tINT2) - - -cIm2D Lapl(cIm2D aImIn) -{ - cIm2D aRes(aImIn.DIm().Sz()); - - Im2D aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); - Im2D aV1Res = cMMV1_Conv::ImToMMV1( aRes.DIm()); - - ELISE_COPY(aV1In.all_pts(),Laplacien(aV1In.in_proj()),aV1Res.out()); - - return aRes; -} - -void MakeStdIm8BIts(cIm2D aImIn,const std::string& aName) +void MMV1_MakeStdIm8BIts(cIm2D aImIn,const std::string& aName) { Im2D aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); double aS0,aS1,aS2; @@ -194,74 +205,93 @@ void MakeStdIm8BIts(cIm2D aImIn,const std::string& aName) -template cPt2dr ValExtre(cIm2D aImIn) +template void Tpl_BenchImFilterV1V2(const cPt2di& aSz) { - // cDataIm2D aDIm(aImIn.DIm()); - // Im2D aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); - auto aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); - double aVMin,aVMax; - ELISE_COPY(aV1In.all_pts(),aV1In.in(),VMin(aVMin)|VMax(aVMax)); - return cPt2dr(aVMin,aVMax); -} + cIm2D aI1(aSz,nullptr,eModeInitImage::eMIA_Rand); + tREAL8 aSigma = 5; + aI1.DIm().InitDirac(aSz/2, Square(aSigma) * 10000); + aI1=aI1.GaussFilter(aSigma,20); // filter image, because corner on a pure random is too noisy + BenchImFilterV1V2(MMV1_Lapl(aI1),Lapl(aI1),2e-3); -template void SelfLabMaj(cIm2D aImIn,const cBox2di & aBox) -{ - cPt2dr aVE = ValExtre(aImIn); - int aMin = round_ni(aVE.x()); - int aMax = round_ni(aVE.y()); - auto aV1In = cMMV1_Conv::ImToMMV1(aImIn.DIm()); - ELISE_COPY - ( - aV1In.all_pts(), - aMin+label_maj(aV1In.in_proj()-aMin,1+aMax-aMin,ToMMV1(aBox)), - aV1In.out() - ); + cIm2D aCV1 = MMV1_CourbTgt(aI1); + SelfCourbTgt(aI1); + BenchImFilterV1V2(aCV1,aI1,1e-3); } - -template cPt2dr ValExtre(cIm2D aImIn); -template void SelfLabMaj(cIm2D aImIn,const cBox2di & aBox); - -//============================== - -void MakeImageDist(cIm2D aImIn,const std::string & aNameChamfer) +template void Tpl_BenchDericheV1V2(const cPt2di& aSz,tREAL8 aAlpha) { - const Chamfer & aChamf = Chamfer::ChamferFromName(aNameChamfer); - aChamf.im_dist(cMMV1_Conv::ImToMMV1(aImIn.DIm())) ; -} + cIm2D aI1(aSz,nullptr,eModeInitImage::eMIA_Rand); + + // we dont make check on diff V1/V2, because it's difficult and as its purely a visual tool + // not important + if (0) + { + MakeStdIm8BIts(aI1,"Grad-Norm-V2.tif"); + MMV1_MakeStdIm8BIts(aI1,"Grad-NormV1.tif"); + } + tREAL8 aSigma = 5; + aI1.DIm().InitDirac(aSz/2, Square(aSigma) * 10000); + aI1=aI1.GaussFilter(aSigma,20); // filter image, because corner on a pure random is too noisy + + cImGrad aDerV2 = Deriche(aI1.DIm(),aAlpha); + cImGrad aDerV1 = MMV1_Deriche(aI1.DIm(),aAlpha); + + MMVII_INTERNAL_ASSERT_bench(aDerV2.mDGx->LInfDist(*(aDerV1.mDGx))<1e-2,"BenchDericheVBenchDericheV"); + MMVII_INTERNAL_ASSERT_bench(aDerV2.mDGy->LInfDist(*(aDerV1.mDGy))<1e-2,"BenchDericheVBenchDericheV"); + + // a Test to justify the IGX[i*DX + j] /= p.mAmpl in deriche_ll + if (0) + { + aDerV2.mDGx->ToFile("GradX-V2.tif"); + aDerV2.mDGy->ToFile("GradY-V2.tif"); + aDerV1.mDGx->ToFile("GradX-V1.tif"); + aDerV1.mDGy->ToFile("GradY-V1.tif"); + aI1.DIm().ToFile("GradOri.tif"); + + cPt2di aP= aSz/2 - cPt2di(1*aSigma,0); + tREAL8 aDif = (aI1.DIm().GetV(aP+cPt2di(1,0)) -aI1.DIm().GetV(aP+cPt2di(-1,0))) / 2.0; + tREAL8 aGx = aDerV2.mDGx->GetV(aP); + StdOut() << " GGGGG " << aDif << " " << aGx << "\n"; + } +} +template void Tpl_BenchInfoExtr_V1V2(const cPt2di& aSz) +{ + cIm2D aI1(aSz,nullptr,eModeInitImage::eMIA_Rand); + auto [aV1Min,aV1Max] = MMV1_ValExtre(aI1); + auto [aV2Min,aV2Max] = ValExtre(aI1); + MMVII_INTERNAL_ASSERT_bench(aV1Min == aV2Min,"Tpl_BenchInfoExtr VMin"); + MMVII_INTERNAL_ASSERT_bench(aV1Max == aV2Max,"Tpl_BenchInfoExtr VMax"); -//============================== + MMVII_INTERNAL_ASSERT_bench(std::fabs(MMV1_MoyAbs(aI1) - MoyAbs(aI1))<=1e-8 ,"Tpl_BenchInfoExtr VMax"); +} -void ExportHomMMV1(const std::string & aIm1,const std::string & aIm2,const std::string & SH,const std::vector & aVP) +void BenchImFilterV1V2() { - ElPackHomologue aPackH; + Tpl_BenchImFilterV1V2(cPt2di(100,110)); - for (int aK=0 ; aK(cPt2di(10,8)); + Tpl_BenchInfoExtr_V1V2(cPt2di(10,12)); - std::string aKeyH = "NKS-Assoc-CplIm2Hom@"+ SH + "@dat"; - cInterfChantierNameManipulateur* anICNM = cInterfChantierNameManipulateur::BasicAlloc("./"); - - std::string aNameF = anICNM->Assoc1To2(aKeyH,aIm1,aIm2,true); - - aPackH.StdPutInFile(aNameF); + Bench_LabMaj(); + //Tpl_BenchDericheV1V2(cPt2di(100,120),2.0); + //Tpl_BenchDericheV1V2(cPt2di(140,120),1.0); + Tpl_BenchDericheV1V2(cPt2di(240,220),2.0); } -void ExportHomMMV1(const std::string & aIm1,const std::string & aIm2,const std::string & SH,const std::vector & aVI) + +#else +void BenchImFilterV1V2() { - std::vector aVR; - for (const auto & aPI: aVI) - aVR.push_back(ToR(aPI)); - ExportHomMMV1(aIm1,aIm2,SH,aVR); + Bench_LabMaj(); } - +#endif + + +//============================== + }; diff --git a/MMVII/src/MorphoMat/ConnectedComponents.cpp b/MMVII/src/MorphoMat/ConnectedComponents.cpp index e3558dbe65..07b8a38306 100644 --- a/MMVII/src/MorphoMat/ConnectedComponents.cpp +++ b/MMVII/src/MorphoMat/ConnectedComponents.cpp @@ -61,4 +61,15 @@ void ConnectedComponent ConnectedComponent(aVPts,aDIm,aNeighbourhood,aVSeed,aMarqInit,aNewMarq,Restore); } +void MakeImageDist(cIm2D aImIn,const std::string & aNameChamfer) +{ + // Apparently Distance image are not used for now, will write it if needed + MMVII_INTERNAL_ERROR("No MakeImageDist"); + // const Chamfer & aChamf = Chamfer::ChamferFromName(aNameChamfer); + // aChamf.im_dist(cMMV1_Conv::ImToMMV1(aImIn.DIm())) ; +} +/* +*/ + + }; diff --git a/MMVII/src/UtiMaths/uti_rand.cpp b/MMVII/src/UtiMaths/uti_rand.cpp index d7eada98ca..fd77c3a54b 100755 --- a/MMVII/src/UtiMaths/uti_rand.cpp +++ b/MMVII/src/UtiMaths/uti_rand.cpp @@ -233,6 +233,11 @@ double RandUnif_N(int aN) return cRandGenerator::TheOne()->Unif_N(aN); } +int RandUnif_M_N(int aM,int aN) +{ + return RandUnif_N(1+aN-aM) + aM; +} + bool HeadOrTail() { return RandUnif_0_1() > 0.5;