Ver Fonte

commentaires

Thierry Dumont há 2 anos atrás
pai
commit
4c460c956a
3 ficheiros alterados com 39 adições e 16 exclusões
  1. 21 10
      RBFInterp.hpp
  2. 7 1
      RBF_functions.hpp
  3. 11 5
      main.cpp

+ 21 - 10
RBFInterp.hpp

@@ -13,15 +13,23 @@ using namespace lapack_c;
 // Interpolation rbf.
 //-------------------------------------------------------------------------
 template<typename real,class Radial>  class RBFInterp{
-  double * matrix , *b;
-  real * image;
+  // on met tout dans une classe paramétrée par le type de flottants,
+  // et par Radial qui est le type (cllase) de fonction d'interpolation.
+  double * matrix , *b; // pour le système linaéire
+  real * image; // l'"image"
   real eps;
-  int imgsize_x,imgsize_y;
+  int imgsize_x,imgsize_y; // tailles
+  // work sert à lapack (tableau de travail).
   double rcond,innorm; double *work; bool rcondok;
+  
 public:
+  // Ce qui définit les points de mesure
   // (x,y) -> valeur
-  // Il ne faut pas d'ex-aequeo pour (x,y), d'où l'utilisation d'une map
+  // NB:
+  // Il ne faut pas d'ex-aequo pour (x,y), d'où l'utilisation d'une map
   // pour décrire les noeuds de l'interpolation.
+  //Clé de la map: paire d'entiers (ordonnée lexicographiquement)
+  // Les clés sont uniques.
   typedef std::map<pair<int,int>,real> PointValue;
 
   // consructeur:
@@ -33,7 +41,7 @@ public:
   // comp_cond : on calcule ou non le conditionnement L\infinity de
   // la matrice.
   {
-    Radial  R(eps);
+    Radial  R(eps);  // fontion d'interpolation, à eps fixé.
     int np= values.size();
    
     matrix = new double[np*np];
@@ -43,18 +51,22 @@ public:
     // matrice).
 
     // 1) matrice:
-    auto  indix= [np](int i, int j){return i*np+j;};
+    auto  indix= [np](int i, int j){return i*np+j;};//lambda function.
     int i=0;
     for(typename PointValue::const_iterator I=values.begin();
-	I!=values.end();I++)
+	I!=values.end();I++) //parcourir la map. Boucle sur les fonctions
       {
-	auto v= I->first; //la fonction est centrée en ce point.
+	auto v= I->first; //la fonction est centrée en ce point (paire i,j)
 	int j=i;
 	for(typename PointValue::const_iterator J=I;J!=values.end();++J)
+	  // boucle sur les points
 	  {
+	    //calculer la distance r entre le point I et le point J
+	    // et y applique la fonction radiale ->s.
+	   
 	    auto s = R(r(v,J->first));
 	    matrix[indix(i,j)]=s;
-	    if(i!=j)
+	    if(i!=j)//matrice symétrique
 	      matrix[indix(j,i)] =s;
 	    ++j;
 	  }
@@ -127,7 +139,6 @@ public:
   }
   real Interpolator(PointValue& values,int nx,int i,int j,Radial R)
   {
-    //auto  indix= [nx](int i, int j){return i*nx+j;};
     real s=0.0;
     int k=0;
     for(typename PointValue::const_iterator K=values.begin();

+ 7 - 1
RBF_functions.hpp

@@ -1,6 +1,7 @@
 #ifndef RBF_functions_hpp
 #define RBF_functions_hpp
 // Les RBF
+// classe de base
 template<typename real> class RBF{
   // Radial Basis Functions: base.
   protected:
@@ -14,11 +15,16 @@ template<typename real> class MQI: public RBF<real> {
   // Multiquadratique inverse.
  
 public:
-  using RBF<real>::set_epsilon;
+  using RBF<real>::set_epsilon; //utiliser le set-epsilon de la classe
+  // de base (C++ 11 !).
+
+  // constructeur: on appelle de constructeur de la base.
   MQI(real eps): RBF<real>(eps){}
+  // la veur en x (distance au noeud)
   inline real operator()(real x){return 1./sqrt(1.+
 						RBF<real>::epsilon*x*x);}
   };
+///
 template<typename real> class TPS: public RBF<real> {
   // plaque mince.
 public:

+ 11 - 5
main.cpp

@@ -19,8 +19,9 @@ int main()
 {
 
   typedef float real;
-  //
+  //  voir RBFInterp qui définit ceci:
   typedef typename RBFInterp<real,Gauss<real>>::PointValue PointValue;
+  // qui est une map (i,j) -> valeur flottante.
 
   // taille de l'"image":
   int nx=1500,ny=1500;
@@ -35,11 +36,13 @@ int main()
   int npts= 50;
   for(int i=0;i<npts;i++)
     {
+      // la position des points est choisie au hasard:
       auto vrand=float(std::rand())/RAND_MAX;
       int xx= static_cast<int>(nx*vrand);
       vrand=float(std::rand())/RAND_MAX;
       int yy= static_cast<int>(ny*vrand);
       auto k=make_pair(xx,yy);
+      // valeur de l'image en ce point.
       Values[k]=finit(xx,yy);
       if(xx<0 || xx>=nx||yy<0 || yy>=ny)
 	// vérifier que les points sont dans l'image.
@@ -66,12 +69,13 @@ int main()
     " nb. points de test : "<<Test.size()<<endl;
   
   // interpolant:
+  // définition de la fontion radiale utilisée (MQI semble mailleure)
   typedef MQI<real> typRad;
-  //  typedef TPS<real> typRad;
+  //typedef TPS<real> typRad;
   //typedef TPSD<real> typRad;
   //typedef Gauss<real> typRad;
 
-  real eps= 0.01;
+  real eps= 0.01; //valeur de départ de eps
   RBFInterp<real,typRad> A(eps);
 
   //
@@ -81,7 +85,7 @@ int main()
 
 
   //--1--- Détermination d'un epsilon "optimal" pour les RBF.
-  ofstream f; f.open("eps-err");
+  ofstream f; f.open("eps-err");//pour faire un graphique precision(epsilon)
   start=get_time();
   real epsm= eps;
   real testm=1.e10;
@@ -118,6 +122,7 @@ int main()
       f<<eps<<" "<<test<<endl;
       if(test<testm)
 	{
+	  // la valeur epsm choisie pour eps.
 	  testm=test;
 	  epsm=eps;
 	}
@@ -144,7 +149,8 @@ int main()
   cout<<"interpolants ok"<<endl;
   
   start=get_time();
-  int increment=2;
+  int increment=2; //pour utiliser l'interp. quad.
+  // (un point sur 2 dans caque direction)
   A.Interp(Values,image,nx,ny,increment);
   stop=get_time();
   cout <<"--- (durée)  interpolation RBF : "<<stop-start<<endl;