Radial-Functions-Interpolation/interp.hpp

#ifndef INTERP_HPP
#define INTERP_HPP
template<class real> inline void quadrati(real *vnodes,
					  real& first, real& second)
{

  first  =  0.375*vnodes[0] + 0.75*vnodes[1] - 0.125*vnodes[2];
  second = -0.125*vnodes[0] + 0.75*vnodes[1] + 0.375*vnodes[2];


}
template<class real> inline void interp1(real *temp, real *image,int nx,
					 int i,int j)
{
  auto indix= [nx](int i, int j){return i*nx+j;};
  int I=i-2,J=j-2;
  //1) ligne horizontale en  haut :
  temp[0]=image[indix(I-2,J-2)];
  temp[1]=image[indix(I-2,J)];
  temp[2]=image[indix(I-2,J+2)];
  quadrati(temp,image[indix(I-2,J-1)],image[indix(I-2,J+1)]);
	
  //2) colonne verticale gauche :
  temp[0]=image[indix(I-2,J-2)];
  temp[1]=image[indix(I,J-2)];
  temp[2]=image[indix(I+2,J-2)];
  quadrati(temp,image[indix(I-1,J-2)],image[indix(I+1,J-2)]);

  //3) colonne verticale centrale :
  temp[0]=image[indix(I-2,J)];
  temp[1]=image[indix(I,J)];
  temp[2]=image[indix(I+2,J)];
  quadrati(temp,image[indix(I-1,J)],image[indix(I+1,J)]);

  //4)  colonne verticale à droite :
  temp[0]=image[indix(I-2,J+2)];
  temp[1]=image[indix(I,J+2)];
  temp[2]=image[indix(I+2,J+2)];
  quadrati(temp,image[indix(I-1,J+2)],image[indix(I+1,J+2)]);

  //5) première ligne:
  temp[0]=image[indix(I-2,J-2)];
  temp[1]=image[indix(I-2,J)];
  temp[2]=image[indix(I-2,J+2)]; //
  quadrati(temp,image[indix(I-2,J-1)],image[indix(I-2,J+1)]);
	
  //6) deuxième ligne:
  temp[0]=image[indix(I-1,J-2)];
  temp[1]=image[indix(I-1,J)];
  temp[2]=image[indix(I-1,J+2)]; //
  quadrati(temp,image[indix(I-1,J-1)],image[indix(I-1,J+1)]);

  //7) troisième ligne
  temp[0]=image[indix(I,J-2)];
  temp[1]=image[indix(I,J)];
  temp[2]=image[indix(I,J+2)]; //
  quadrati(temp,image[indix(I,J-1)],image[indix(I,J+1)]);
	
  //8) quatrième  ligne:
  temp[0]=image[indix(I+1,J-2)];
  temp[1]=image[indix(I+1,J)];
  temp[2]=image[indix(I+1,J+2)];
  quadrati(temp,image[indix(I+1,J-1)],image[indix(I+1,J+1)]);
}
// Un point sur deux étant calculé  dans image, interpoler qudratiquement
// les autres points.
// image : l'image
// nx, ny: les dimensions de l'image
template<class real> void interpquad(real* image,int nx,int ny)
{
  
  // tester si nx et ny sont divisibles par 4
  // à faire!
  real temp[3];
  
  for(int i=4;i<nx;i+=4)
    for(int j=4;j<ny;j+=4)
      interp1(temp,image,nx,i,j);

  for(int j=4;j<ny;j+=4)
      interp1(temp,image,nx,nx-1,j);
  for(int i=4;i<nx;i+=4)
    interp1(temp,image,nx,i,ny-1);
 
       
}
#endif