Merge branch 'mpd'

MMVII/Doc/Communication/Texte-formation-2024.txt

@@ -1,14 +1,22 @@
 Session de programmation sous MicMac version 2 (MMVII)
 
 
-Date : 3 jours pleins du 12 au 14 mars
+Date : 3 jours pleins du 12 au 14 mars, 
 Sujet : implanter dans MMVII un ajustement de faisceau sur image satelitte
 Objectif : apprendre l'architecture de MMVII, notamment vis à vis de l'otimisation non linéaire
-Public : développeurs ayant de bonnes connaissances en photogrammétrie à l'aise en C++ 
-Contact : [email protected]
+Public : développeurs ayant de bonnes connaissances en photogrammétrie et à l'aise en C++.
+organisation : en présentiel (IGN) et en visio (IGN et extérieur).
+Contact : [email protected], [email protected]
 
 
 
-Contexte : MicMac est une solution de photogrammétrie open source complète développée à l'IGN depuis 2003.  Comme beaucoup de logiciel de recherche il s'est construit pas empilement successif 
+Contexte : MicMac est une solution de photogrammétrie open source complète développée à l'IGN depuis 2003. Une version 2 visant à être à faciliter les contributions externes et à être plus maintenable sur le long terme est en développement depuis 2020. L'équipe de développement comporte actuellement 6 chercheurs et ingénieurs IGN (pour environs 2 ETP), le projet est actuellement soutenu par l'IGN et le CERN. Afin de promouvoir  des contributions externes, l'IGN organise des séances de programmation interactive (live-coding); après une session en novembre 2023,  centrée sur la prise de compte de la contrainte de bloc rigide, une session est prévue en mars 2024 sur l'ajustement de faisceaux sur des capteurs satelittes.
 
 
+
+
+
+ Après une séance organisée 
+
+  Comme beaucoup de logiciel de recherche il s'est construit pas empilement successif de couches et
+

MMVII/Doc/Methods/PushBroomSensor-Theory.tex

@@ -129,6 +129,8 @@ \subsection{Satellite model}
 
 \subsection{Satellite error model}
 
+\label{PB:SatErMOd}
+
 Let $R(i)$ be the initial orientation , and $R'(i)$ be the "real orientation"
 we want to estimate by adjsument.
 We admit that there exist a rotation $R_\Theta(i)$, very close to identity and
@@ -156,7 +158,7 @@ \subsection{Satellite error model}
 
 $\Theta(i)$ will be the unknown we want to compute, but for now
 we suppose we know it.  Eventually, later we will add intrinsic calibration
-unknown, for example if we add a scaling to modelize the focal, we will 
+unknown, for example if we add an unknown scaling $\lambda$ to modelize the focal, we will 
 have $\Theta=(\omega,\phi,\kappa, \lambda)$
 
 %-----------------------------------------------------------------------
@@ -229,10 +231,43 @@ \subsection{Compute projection function with perturbation}
     =  - \mathbb{J}_q ^{-1} \mathbb{J}_\Theta \begin{bmatrix} \omega \\ \phi \\ \kappa  \end{bmatrix}
 \end{equation}
 
+
+We note :
+\begin{equation}
+	\mathcal{J}(i,j,z) = - \mathbb{J}_q ^{-1} \mathbb{J}_\Theta
+\end{equation}
+
 We now have a way to compute $\pi'$ (to first order) : 
 
 \begin{equation}
-	 \pi'(P) =  \pi(P)  - \mathbb{J}_q ^{-1} \mathbb{J}_\Theta \begin{bmatrix} \omega \\ \phi \\ \kappa  \end{bmatrix}
+	 \pi'(P) =  \pi(P)  + \mathcal{J} \begin{bmatrix} \omega \\ \phi \\ \kappa  \end{bmatrix}
+\end{equation}
+
+%-----------------------------------------------------------------------
+
+\subsection{Correction model}
+
+As decribed in~\ref{PB:SatErMOd}, we assume that the corrective term between
+initial localisation and \emph{real} localisation :
+
+\begin{enumerate}
+        \item  is a rotation arround current centers, 
+	\item  that this rotation depends only of $i$,
+	\item  that it  is very close to identity 
+	\item  that it  has very smooth variation.
+\end{enumerate}
+
+The small rotation being coded by $\omega,\phi,\kappa$, the classical solution
+is to  select a  familly of smooth function and to express them as a linear combinaison
+of these smooth function.  Classically we select the polynoms as basis (also other
+alternative may be studied later).  So let $N$ be the degree of polynoms,
+and $a_i,b_i,c_i$  the coefficient, the projection fonction is :
+
+All these consideration  drive to select for $R(i)$ a pa
+
+\begin{equation}
+	 \pi'(i,j,z) =  \pi(i,j,z)  + \sum _{k=0}^{N}  \mathcal{J}(i,j,z)  
+                                                \begin{bmatrix} a_k \omega^k \\ b_k \phi^k  \\ c_k \kappa^k  \end{bmatrix}
 \end{equation}
 %-----------------------------------------------------------------------
 \section{Vrac}
@@ -243,3 +278,6 @@ \section{Vrac}
 \end{itemize}
 
 
+
+
+

MMVII/include/MMVII_Bench.h

@@ -127,6 +127,7 @@ void BenchHamming(cParamExeBench & aParam);
 void BenchPolynome(cParamExeBench & aParam);
 
 
+void BenchStenopeSat();
 void BenchUnCalibResection();
 void BenchPoseEstim(cParamExeBench & aParam);
 

MMVII/include/MMVII_Geom3D.h

@@ -132,10 +132,21 @@ template <class Type> class cRotation3D
        // Extract Axes of a rotation and compute its angle 
        void ExtractAxe(tPt & anAxe,Type & aTeta);
 
-       // conversion to Omega Phi Kapa
+       /// conversion to Omega Phi Kapa
        static cRotation3D<Type>  RotFromWPK(const tPt & aWPK);
+       /// extrecat Omega Phi Kapa from rotation
        tPt                       ToWPK() const;
 
+       /// Rotation arround X
+       static cDenseMatrix<Type> RotOmega(const tREAL8 & aOmega);
+       /// Rotation arround Y
+       static cDenseMatrix<Type> RotPhi(const tREAL8 & aPhi);
+       /// Rotation arround Z
+       static cDenseMatrix<Type> RotKappa(const tREAL8 & aKappa);
+
+       ///  0-> Omega   1->Phi  2-> Kappa
+       static cDenseMatrix<Type> Rot1WPK(int aK,const tREAL8 & aOmega);
+
        // conversion to Yaw Pitch Roll
        static cRotation3D<Type>  RotFromYPR(const tPt & aWPK);
        tPt                       ToYPR() const;
@@ -144,6 +155,7 @@ template <class Type> class cRotation3D
        cDenseMatrix<Type>  mMat;
 };
 
+typedef cRotation3D<tREAL8> tRotR; 
 
 /**  Class for 3D "affine" rotation of vector
 

MMVII/include/MMVII_Images.h

@@ -380,6 +380,7 @@ template <class Type,const int Dim> class cDataTypedIm : public cDataGenUnTypedI
         tREAL16  MoyVal() const;
 
         void DupIn(cDataTypedIm<Type,Dim> &) const;  ///< Duplicate raw data
+        void ChSignIn(cDataTypedIm<Type,Dim> &) const;  ///< Duplicate with sign change in raw data
         void DupInVect(std::vector<Type> &) const;  ///< Duplicate raw data in a vect
 
         // Defaults values quitt slow but may be usefull

MMVII/include/MMVII_Ptxd.h

@@ -80,6 +80,8 @@ template <class Type,const int Dim> class cPtxd
                 aRes.mCoords[aK]= aVal;
            return aRes;
        }
+       /// All coord 0 exect aCoord that values aVal
+       static cPtxd<Type,Dim>  P1Coord(size_t aCoord,const Type & aVal) ;
        /// Initialisation with nan value (to detect error asap)
        static cPtxd<Type,Dim>  Dummy();
        /// Initialisation from name "i..."  "-j..."    valide are "ijkl" 
@@ -832,8 +834,12 @@ template <class Type,const int Dim> class cSegmentCompiled : public cSegment<Typ
     public :
        typedef cPtxd<Type,Dim> tPt;
        cSegmentCompiled(const tPt& aP1,const tPt& aP2);
+       cSegmentCompiled(const cSegment<Type,Dim>&);
        tPt  Proj(const tPt &) const;
        Type Dist(const tPt &) const;
+
+       const Type & N2 () const;
+       const tPt  & Tgt() const;
     protected :
        Type    mN2;
        tPt     mTgt;

MMVII/include/MMVII_Sensor.h

@@ -90,14 +90,24 @@ class cSensorImage  :   public cObj2DelAtEnd,
 
 	 // =================   Image <-->  Ground  mappings  ===========================
 
-         ///  The most fundamental method, theoretically should be sufficient
+         /*  The most fundamental method, theoretically should be sufficient, when meaning full (ie non orthocentric) 
+            orientation must go from sensor to the scene */
+
          virtual tSeg3dr  Image2Bundle(const cPt2dr &) const =0;
 	 /// Basic method  GroundCoordinate ->  image coordinate of projection
          virtual cPt2dr Ground2Image(const cPt3dr &) const = 0;
+
 	 ///  add the the depth (to see if have a default with bundle+Gr2Ima)
          virtual cPt3dr Ground2ImageAndDepth(const cPt3dr &) const = 0;
          /// Invert of Ground2ImageAndDepth
          virtual cPt3dr ImageAndDepth2Ground(const cPt3dr &) const = 0;
+
+	 ///  add the the Z
+         virtual cPt3dr Ground2ImageAndZ(const cPt3dr &) const ;
+         /// Invert of Ground2ImageAndZ
+         virtual cPt3dr ImageAndZ2Ground(const cPt3dr &) const ;
+
+
 	 /// Facility for calling ImageeAndDepth2Ground(const cPt3dr &)
          cPt3dr ImageAndDepth2Ground(const cPt2dr &,const double & ) const;
 
@@ -169,6 +179,14 @@ class cSensorImage  :   public cObj2DelAtEnd,
 	 /// If the camera has its own "obs/cste" (like curent rot for PC-Cam) that's the place to say it
 	 virtual  void PushOwnObsColinearity( std::vector<double> &) = 0;
 
+         //  method used in push-broom perturbation model
+
+          /// Extract the pose a line of sensor, meaningfull for push-broom , but can be used in other  context
+          tPoseR   GetPoseLineSensor(tREAL8 aXorY,bool LineIsX,int aNbSample,bool*IsOk=0,std::vector<double>*aResCD=0) const;
+          /**  Once computed the Orientation of line, compute the differential of proj relatively to WPK
+               EpsXYZ and EpsWPK are used for computing derivative with relative differences */
+           cDenseMatrix<tREAL8> CalcDiffProjRot(const cPt3dr & aPt,const tPoseR &,const cPt3dr & aEpsXYZ,const tREAL8 & aEpsWPK) const;
+
      private :
 	 /// Return the calculator, adapted to the type, for computing colinearity equation
          virtual cCalculator<double> * CreateEqColinearity(bool WithDerives,int aSzBuf,bool ReUse) = 0;

MMVII/include/MMVII_Tpl_Images.h

@@ -139,6 +139,20 @@ template<class T2>   cDenseVect<T2> operator -= (cDenseVect<T2> & aI2,const cDen
 	return aI2;
 }
 
+template<class T1>   cIm2D<T1> operator - (const cIm2D<T1> & aIm)
+{
+   cIm2D<T1> aRes (aIm.DIm().Sz());
+   aIm.DIm().ChSignIn(aRes.DIm());
+
+   return aRes;
+}
+
+template<class T1>   cDenseMatrix<T1> operator - (const cDenseMatrix<T1> & aMat)
+{
+    return  cDenseMatrix<T1> (-aMat.Im());
+}
+
+
        //===========   Addition ===========
 
 template<class T1,class T2,class T3,int Dim>  

MMVII/include/MMVII_memory.h

@@ -146,20 +146,13 @@ class cMemManager
 };
 
 
+/** cMemCountable to check that all object allocated are destroyed */
 
-
-/**
-    This classe redefine operators new and delate to check alloc / desalloc and (some) bad access.
-    Allocation and desallocation is delegated to  cMemManager
-*/
-
-class  cMemCheck
+class cMemCountable
 {
-      public :
-         void * operator new    (size_t sz);
-         void operator delete   (void * ptr) ;
+     public :
 #if (The_MMVII_DebugLevel >= The_MMVII_DebugLevel_InternalError_tiny)
-         cMemCheck()  
+         cMemCountable()  
          {
               mCpt = TheCptObj;
               mActiveNbObj=   cMemManager::IsActiveMemoryCount();
@@ -169,8 +162,8 @@ class  cMemCheck
               }
 	      TheCptObj++;
          }
-         cMemCheck(const cMemCheck &)  : cMemCheck () {}
-         ~cMemCheck()
+         cMemCountable(const cMemCountable &)  : cMemCountable () {}
+         ~cMemCountable()
          {
             if (mActiveNbObj)
             {
@@ -184,13 +177,28 @@ class  cMemCheck
       private :
          static int    TheNbObjLive;
          static int    TheCptObj;
+};
+
+
+/**
+    This classe redefine operators new and delate to check alloc / desalloc and (some) bad access.
+    Allocation and desallocation is delegated to  cMemManager
+*/
+
+class  cMemCheck : public cMemCountable
+{
+      public :
+         void * operator new    (size_t sz);
+         void operator delete   (void * ptr) ;
+         cMemCheck () : cMemCountable() {}
+         cMemCheck(const cMemCheck &)  : cMemCheck () {}
        // to avoid use 
 };
 
 /**  some object have to live untill the end of the process, just before the verif is done in main
  *   appli destructor.
  */
-class cObj2DelAtEnd
+class cObj2DelAtEnd  : public cMemCountable
 {
         public :
                 virtual ~cObj2DelAtEnd();

MMVII/src/Appli/cSpecMMVII_Appli.cpp

@@ -86,12 +86,12 @@ int cSpecMMVII_Appli::AllocExecuteDestruct(const std::vector<std::string> & aVAr
         }
     }
     cMemManager::CheckRestoration(aMemoState);
-    MMVII_INTERNAL_ASSERT_always(cMemCheck::NbObjLive()==aNbObjLive,"Mem check obj not killed");
+    MMVII_INTERNAL_ASSERT_always(cMemCountable::NbObjLive()==aNbObjLive,"Mem check obj not killed");
     aCptCallIntern--;
     // This was the initial test, stricter, maintain it when call by main
     if (aCptCallIntern==0)
     {
-         MMVII_INTERNAL_ASSERT_always(cMemCheck::NbObjLive()==0,"Mem check obj not killed");
+         MMVII_INTERNAL_ASSERT_always(cMemCountable::NbObjLive()==0,"Mem check obj not killed");
     }
     return aRes;
 }

MMVII/src/Geom3D/cRotation3D.cpp

@@ -386,19 +386,45 @@ template <class Type> cPtxd<Type,3>  cRotation3D<Type>::ToYPR() const
 	return tPt(-aWPK.z(),-aWPK.y(),-aWPK.x());
 }
 
-template <class Type> cRotation3D<Type>  cRotation3D<Type>::RotFromWPK(const tPt & aWPK)
+template <class Type> cDenseMatrix<Type>  cRotation3D<Type>::RotOmega(const tREAL8 & aOmega)
+{
+   Type aCx = std::cos(aOmega);
+   Type aSx = std::sin(aOmega);
+   return M3x3FromLines(tPt(1,0,0),tPt(0,aCx,-aSx),tPt(0,aSx,aCx));
+}
+
+template <class Type> cDenseMatrix<Type>  cRotation3D<Type>::RotPhi(const tREAL8 & aPhi)
 {
-   Type aCx = std::cos(aWPK.x());
-   Type aSx = std::sin(aWPK.x());
-   auto aRx = M3x3FromLines(tPt(1,0,0),tPt(0,aCx,-aSx),tPt(0,aSx,aCx));
+   Type aCy = std::cos(aPhi);
+   Type aSy = std::sin(aPhi);
+   return M3x3FromLines(tPt(aCy,0,aSy),tPt(0,1,0),tPt(-aSy,0,aCy));
+}
 
-   Type aCy = std::cos(aWPK.y());
-   Type aSy = std::sin(aWPK.y());
-   auto aRy = M3x3FromLines(tPt(aCy,0,aSy),tPt(0,1,0),tPt(-aSy,0,aCy));
+template <class Type> cDenseMatrix<Type>  cRotation3D<Type>::RotKappa(const tREAL8 & aKappa)
+{
+   Type aCz = std::cos(aKappa);
+   Type aSz = std::sin(aKappa);
+   return M3x3FromLines(tPt(aCz,-aSz,0),tPt(aSz,aCz,0),tPt(0,0,1));
+}
 
-   Type aCz = std::cos(aWPK.z());
-   Type aSz = std::sin(aWPK.z());
-   auto aRz = M3x3FromLines(tPt(aCz,-aSz,0),tPt(aSz,aCz,0),tPt(0,0,1));
+template <class Type> cDenseMatrix<Type>  cRotation3D<Type>::Rot1WPK(int aK,const tREAL8 & aTeta)
+{
+    switch (aK)
+    {
+        case 0 : return RotOmega(aTeta);
+        case 1 : return RotPhi(aTeta);
+        case 2 : return RotKappa(aTeta);
+    }
+    MMVII_INTERNAL_ERROR("Bad value of K in Rot1WPK : " + ToStr(aK));
+    return  cDenseMatrix<Type>(3);
+}
+
+
+template <class Type> cRotation3D<Type>  cRotation3D<Type>::RotFromWPK(const tPt & aWPK)
+{
+   auto aRx = RotOmega(aWPK.x());
+   auto aRy = RotPhi(aWPK.y());
+   auto aRz = RotKappa(aWPK.z());
 
    return cRotation3D<Type>(aRx*aRy*aRz,false);
 /*