From 5938966c037de1185d36c5e3c4e712943a7f2601 Mon Sep 17 00:00:00 2001
From: deseilligny <marc.pierrot-deseilligny@ensg.eu>
Date: Wed, 14 Mar 2018 21:00:07 +0100
Subject: [PATCH] Hom filter + Rolling sur FITS

---
 src/uti_image/NewRechPH/NH_InvarRad.cpp | 75 +++++++++++++++++++++----
 src/uti_image/NewRechPH/NewRechPH.cpp   |  8 ++-
 src/uti_image/NewRechPH/NewRechPH.h     |  7 ++-
 src/uti_phgrm/CPP_HomFilterMasq.cpp     |  2 +-
 4 files changed, 77 insertions(+), 15 deletions(-)

diff --git a/src/uti_image/NewRechPH/NH_InvarRad.cpp b/src/uti_image/NewRechPH/NH_InvarRad.cpp
index 460aaa2af7..68a11f5de9 100644
--- a/src/uti_image/NewRechPH/NH_InvarRad.cpp
+++ b/src/uti_image/NewRechPH/NH_InvarRad.cpp
@@ -40,9 +40,14 @@ Header-MicMac-eLiSe-25/06/2007*/
 
 #include "NewRechPH.h"
 
-Pt2di cAppli_NewRechPH::SzInvRad()
+Pt2di cAppli_NewRechPH::SzInvRadUse()
 {
-   return Pt2di(mNbSR, mNbTetaInv);
+   return Pt2di(mNbSR2Use, mNbTetaInv);
+}
+
+Pt2di cAppli_NewRechPH::SzInvRadCalc()
+{
+   return Pt2di(mNbSR2Calc, mNbTetaInv);
 }
 
 /*
@@ -169,6 +174,23 @@ class cRadInvStat
         double mS3;
 };
 
+   // return Pt2di(mNbSR2Use, mNbTetaInv);
+void Normalise(tImNRPH aImBuf,tImNRPH aImOut,int aX0In,int aX1In,int aSzXOut)
+{
+    int aSzY = aImBuf.sz().y;
+    double aS0,aS1,aS2;
+    ELISE_COPY
+    (
+         rectangle(Pt2di(aX0In,0),Pt2di(aX1In,aSzY)),
+         Virgule(1,aImBuf.in(),Square(aImBuf.in())),
+         Virgule(sigma(aS0),sigma(aS1),sigma(aS2))
+    );
+    aS1 /= aS0;
+    aS2 /= aS0;
+    aS2 -= ElSquare(aS1);
+    aS2 = sqrt(ElMax(1e-10,aS2));
+    ELISE_COPY(rectangle(Pt2di(aX0In,0),Pt2di(aX0In+aSzXOut,aSzY)),(aImBuf.in()-aS1)/aS2, aImBuf.out());
+}
 
 
 bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
@@ -182,7 +204,7 @@ bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
    }
 
    // Buf[KTeta][KRho]   pour KRho=0, duplication de la valeur centrale
-   tImNRPH aImBuf(SzInvRad().x,SzInvRad().y);
+   tImNRPH aImBuf(SzInvRadCalc().x,SzInvRadCalc().y);
    tTImNRPH aTBuf(aImBuf);
 
    std::vector<cOneScaleImRechPH *>  aVIm;
@@ -194,7 +216,7 @@ bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
 
    int aN0 = aPt.NivScale();
    // aVIm.push_back(mVI1.at(aN0));
-   for (int aKRho=0 ; aKRho <mNbSR ; aKRho++)
+   for (int aKRho=0 ; aKRho <mNbSR2Calc ; aKRho++)
    {
        aVIm.push_back(mVI1.at(aN0 + aKRho * mDeltaSR));
    }
@@ -238,6 +260,36 @@ bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
       }
    }
 
+   // Normalisation a priori, pour l'instant sans rolling
+   if (mRollNorm)
+   {
+       for (int aX= 0 ; aX<mNbSR2Use ; aX++)
+       {
+            Normalise(aImBuf,aImBuf,aX,aX+mNbSR2Use,1);
+       }
+   }
+   else
+   {
+      Normalise(aImBuf,aImBuf,0,mNbSR2Use,mNbSR2Use);
+   }
+
+/*
+   {
+      double aS0,aS1,aS2;
+      ELISE_COPY
+      (
+         rectangle(Pt2di(0,0),SzInvRadUse()),
+         Virgule(1,aImBuf.in(),Square(aImBuf.in())),
+         Virgule(sigma(aS0),sigma(aS1),sigma(aS2))
+      );
+      aS1 /= aS0;
+      aS2 /= aS0;
+      aS2 -= ElSquare(aS1);
+      aS2 = sqrt(ElMax(1e-10,aS2));
+      ELISE_COPY(aImBuf.all_pts(),(aImBuf.in()-aS1)/aS2, aImBuf.out());
+   }
+*/
+
    if (BUG)
    {
        std::cout << "PTBBUGGGGG 111\n";
@@ -247,7 +299,7 @@ bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
    int aKPS2 = mNbTetaInv /4 ;
    int aKPi  = mNbTetaInv /2 ;
    int aNbGrand = ((int) eTIR_NoLabel) / 3;
-   for (int aKRho=0 ; aKRho<int(aVIm.size()) ; aKRho++)
+   for (int aKRho=0 ; aKRho<mNbSR2Use ; aKRho++)
    {
       double aRealDTeta =  aVDeltaRad[aKRho] / aVDeltaTang[aKRho];
       int aDTeta = round_ni(aRealDTeta); // Delta correspondant a 1 rho
@@ -340,11 +392,11 @@ bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
        std::string aName= aDir + "InvRad" + aNamePt +  ".tif";
 
        L_Arg_Opt_Tiff aLarg;
-       aLarg = aLarg + Arg_Tiff(Tiff_Im::AStrip( arrondi_sup(SzInvRad().x,8)));
+       aLarg = aLarg + Arg_Tiff(Tiff_Im::AStrip( arrondi_sup(SzInvRadCalc().x,8)));
        Tiff_Im  aSaveBuf
        (
            aName.c_str(),
-           SzInvRad(),
+           SzInvRadCalc(),
            GenIm::real4,
            Tiff_Im::No_Compr,
            Tiff_Im::BlackIsZero,
@@ -357,10 +409,11 @@ bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
 
    // Pour l'export xml
    {
+/*
       double aS0,aS1,aS2;
       ELISE_COPY
       (
-         aImBuf.all_pts(),
+         rectangle(Pt2di(0,0),SzInvRadUse()),
          Virgule(1,aImBuf.in(),Square(aImBuf.in())),
          Virgule(sigma(aS0),sigma(aS1),sigma(aS2))
       );
@@ -369,8 +422,10 @@ bool  cAppli_NewRechPH::CalvInvariantRot(cOnePCarac & aPt)
       aS2 -= ElSquare(aS1);
       aS2 = sqrt(ElMax(1e-10,aS2));
       ELISE_COPY(aImBuf.all_pts(),(aImBuf.in()-aS1)/aS2, aImBuf.out());
-      aPt.ImLogPol() =  Im2D_INT1(aImBuf.sz().x,aImBuf.sz().y);
-      ELISE_COPY(aImBuf.all_pts(),Max(-128,Min(127,round_ni(aImBuf.in()*32))),aPt.ImLogPol().out());
+*/
+
+      aPt.ImLogPol() =  Im2D_INT1(SzInvRadUse().x,SzInvRadUse().y);
+      ELISE_COPY(aPt.ImLogPol().all_pts(),Max(-128,Min(127,round_ni(aImBuf.in()*32))),aPt.ImLogPol().out());
       aPt.VectRho() = aVRho;
       
       if (BUG)
diff --git a/src/uti_image/NewRechPH/NewRechPH.cpp b/src/uti_image/NewRechPH/NewRechPH.cpp
index 46d4770ed9..6704339502 100644
--- a/src/uti_image/NewRechPH/NewRechPH.cpp
+++ b/src/uti_image/NewRechPH/NewRechPH.cpp
@@ -124,9 +124,9 @@ cAppli_NewRechPH::cAppli_NewRechPH(int argc,char ** argv,bool ModeTest) :
     mNbS         (30),
     mISF         (-1,1e10),
     mStepSR      (1.0),
-    mNbSR        (10),
+    mNbSR2Use    (10),
     mDeltaSR     (1),
-    mMaxLevR     (mNbS - (mNbSR-1) * mDeltaSR),
+    mRollNorm    (true),
     mNbTetaIm    (16),
     mMulNbTetaInv (4),
     mNbTetaInv   (mNbTetaIm*mMulNbTetaInv),
@@ -202,6 +202,10 @@ cAppli_NewRechPH::cAppli_NewRechPH(int argc,char ** argv,bool ModeTest) :
                       << EAM(mISF, "ISF",true,"Interval scale forced")
    );
 
+    mNbSR2Calc = mRollNorm ?  (mNbSR2Use * 2 - 1) : mNbSR2Use;
+    mMaxLevR     = mNbS - (mNbSR2Calc-1) * mDeltaSR;
+
+
    if (! EAMIsInit(&mExtSave))
    {
       mExtSave  = mBasic ? "Basic" : "Std";
diff --git a/src/uti_image/NewRechPH/NewRechPH.h b/src/uti_image/NewRechPH/NewRechPH.h
index 6ddd87718d..f019a22163 100644
--- a/src/uti_image/NewRechPH/NewRechPH.h
+++ b/src/uti_image/NewRechPH/NewRechPH.h
@@ -162,7 +162,8 @@ class cOneScaleImRechPH
 class cAppli_NewRechPH
 {
     public :
-        Pt2di SzInvRad();
+        Pt2di SzInvRadUse();
+        Pt2di SzInvRadCalc();
 
         cAppli_NewRechPH(int argc,char ** argv,bool ModeTest);
 
@@ -221,7 +222,9 @@ class cAppli_NewRechPH
   Invariant Rotation
 */
         double      mStepSR;  // Pas entre les pixel de reechantillonage pour les desc
-        int         mNbSR;     // Nbre de niveau radial (entre 
+        bool        mRollNorm; // si true , normalisation par fenetre glissante, :
+        int         mNbSR2Use;     // Nbre de niveau pour utilisation
+        int         mNbSR2Calc;     // Nbre de niveau a calculer, different si rolling
         int         mDeltaSR;  // Delta entre deux niveau radiaux, genre 1 ou 2 ?
         int         mMaxLevR;  // Niv max permettant le calcul (calcule a partir des autres)
 
diff --git a/src/uti_phgrm/CPP_HomFilterMasq.cpp b/src/uti_phgrm/CPP_HomFilterMasq.cpp
index 7be4deb26a..905d4b1437 100644
--- a/src/uti_phgrm/CPP_HomFilterMasq.cpp
+++ b/src/uti_phgrm/CPP_HomFilterMasq.cpp
@@ -292,7 +292,7 @@ std::cout << aNameIm1  << " # " << aNameIm2 << "\n";
                               aNbInTer += OkTer;
                           }
 
-                          if (aDistHom >0 )
+                          if (Ok && (aDistHom >0 ))
                           {
                               Pt2dr aRP1 =  aVCam[aKN1]->Ter2Capteur(aPTer);
                               Pt2dr aRP2 =  aVCam[aKN2]->Ter2Capteur(aPTer);