#include "math.h"
#include "resample.h"
#include "strings.h"
/* periodic symmetric extension of src */
#define REFOLD(J,L)   { if(J<0) J =-J; J = J % (2*L-1); if(J>=L) J=(2*L-1)-J;}
/* border extension of src */
#define BORD(J,L)   { if((J)<0) (J)=0; else if ((J)>(L)) (J)=(L); }
/* zero extension of src */
#define ZERO(J,L)   { if(J<0 || J>L) continue; }

/* enforce boundary condition <-> array bounds */
//#define BC(J,L) REFOLD(J,L)



#undef  MAX
#define MAX(a, b)  (((a) > (b)) ? (a) : (b))

#undef  MIN
#define MIN(a, b)  (((a) < (b)) ? (a) : (b))

#undef  ABS
#define ABS(a)     (((a) < 0) ? -(a) : (a))

#undef  CLAMP
#define CLAMP(x, low, high)  (((x) > (high)) ? (high) : (((x) < (low)) ? (low) : (x)))


#include "stdio.h"
#include "math.h"


/* these lines should help if you are using  braindead MS C compiler */
#ifdef MSVC 
#define inline 
#define static 
#warning static and inline are disabled 
#endif



/* void (*resample_array_inv_choices[10])(const double *F,  */
/*    const PIXEL_TYPE *src, int s_len, int s_stride,  */
/*    PIXEL_TYPE *dst, int d_len, int d_stride); */

typedef  void (* resample_t)(const double *F, 
   const PIXEL_TYPE *src, int s_len, int s_stride, 
   PIXEL_TYPE *dst, int d_len, int d_stride);



/*****************************************************************/

// WARNING this is not ANTIALIASED in any way but is very simple and fast

static inline void resample_array_inv_near_neighbor
  (const double *F, 
   const PIXEL_TYPE *src, int s_len, int s_stride, 
         PIXEL_TYPE *dst, int d_len, int d_stride)
{
  int i,p=0,j;
  for(i=0;i<d_len;i++) {
    j=F[i];
    BORD(j,s_len-1)
    dst[p]=src[j* s_stride];
    p+= d_stride;
  }
}


/*********** antialiased using bilinear interpixel**************
 * works well when the scale is not reduced and/or when there are not
 * very fine grains; it is much faster; it is approximate in many points
 * For animations it is better to use antialiased
 */

static void resample_array_inv_bilinear
  (const double *F, 
   const PIXEL_TYPE *src, int s_len, int s_stride, 
         PIXEL_TYPE *dst, int d_len, int d_stride)
{
  int i,p=0,j,nj;
  double v,x,dj;
  for(i=0;i<d_len;i++) {
    x=F[i];
    BORD(x,s_len-1);
    j=floor(x);
    dj=x-j;
    nj=j+1;
    if(nj>=s_len) 
      v=src[j* s_stride];
    else
      v=src[j* s_stride]*(1-dj) + src[nj* s_stride]*dj;    
    dst[p]=v;
    p+= d_stride;
  }
}

#ifdef SEEMS_UNNEEDED
      int c;
      /* I have tested this improvement to 
	 the above, and this seems to introduce ripples;
	 I don't want to improve it, use AA if you really want AA*/
      if( oj>=j+2) {
	c=oj-j+1;
	if(oj>=s_len) oj=s_len-1;
	while (oj>=j) { v+=src[oj* s_stride]; oj--;  }
	v/=c;
      }else 
	if ( oj+2<=j) {
	  c=j-oj+1;
	  if(oj<0) oj=0;
	  while (oj<=j) { v+=src[oj* s_stride]; oj++;  }
	  v/=c;
	} else 
#endif

/*********** antialiased using convolutional kernels***************/

#if 256 >= (PIXEL_MAX-PIXEL_MIN)
#include "sinc_256.h"
#else
#include "sinc_1024.h"
#endif


static double sinc(const double x)
{
  if (x<0.0001 && x > -0.0001)
    return 1;
  else
    {
      double tmp = M_PI * x  ;
      return sin(tmp) / tmp;
    }
}

static double sinc_by_table(const double x)
{
  if (x < -4 || x > 4) {
    double tmp = M_PI * x  ;
    return sin(tmp) / tmp;
  }  else
    return sinc_table [(int)((ABS(x))*SINC_TABLE_UNIT )];
}


static double lanczos(const double x)
{
  if (x < -2 || x > 2)
    return 0;
  else
    return sinc(x);
}

static double lanczos4(const double x)
{
  if (x < -4 || x > 4)
    return 0;
  else
    return sinc(x);
}

static double triangle(const double x)
{
  if (x<-1) return 0;
  else if(x<0) return x+1;
  else if(x<1) return 1-x;
  else return 0;
}



// this generates the function name
#define FUN(A) XFUN(A)
#define XFUN(A) resample_array_inv_ ## A


/*************  choice of antialiasing kernel **************/
/* optimizations  for the case of any kernel  based on sinc
   suggested by lvalero, oberger 05/05/2004 */
#define KERNEL_sinc_fast


/* if the above is undefined, then this will be used instead 
 this also is used in creating the function name*/
#define KERNEL lanczos

/* this is the half of the width of the kernel 
 for lanczos4, it must be 4, otherwise 2*/
#define KERNEL_WIDTH 2

/* this creates the function */
#include "resample_snippet.h"

/* then another function */
#undef KERNEL
#define KERNEL lanczos4
#undef KERNEL_WIDTH
#define KERNEL_WIDTH 4
#include "resample_snippet.h"
/***********end  choice of antialiasing kernel **************/









//#struct { fun resample_t; const char *name;]
resample_t resample_choices[10] =
  { FUN(near_neighbor) ,     FUN(bilinear),
    FUN(lanczos),FUN(lanczos4), 
  NULL};

char * resample_array_inv_names[10] =
  { "near_neighbor" , //choose nearest pixel: fastest, looks bad
    "bilinear", //bilinear: same as with the old libmorph warping code
    "lanczos", //Lanczos:  much better quality, a must for animations
    //and/or fine grained images; it is though  slower,
    "lanczos4",// even better than before, but no noticeable difference on 
    // most images
  NULL};


void (*resample_array_inv)(const double *F, 
   const PIXEL_TYPE *src, int s_len, int s_stride, 
   PIXEL_TYPE *dst, int d_len, int d_stride)= FUN(lanczos);


void
mesh_resample_choose_aa(int f)
{

  //  fprintf(stderr,"\n%s:%d:  choice '%s' for kernel!\n",__FILE__,__LINE__,resample_array_inv_names[f] );
  resample_array_inv=resample_choices[f];
}

#include <string.h> //strcmp

void
mesh_resample_choose_aa_by_name(char * s)
{
  int f=0;
  while(resample_choices[f]) {
    if (0==strcmp (s,resample_array_inv_names[f]))
      {
	resample_array_inv=resample_choices[f];
	return;
      }
    f++;
  }
  
  fprintf(stderr,"\n%s:%d: no choice '%s' for kernel!\n",__FILE__,__LINE__,s);
}






//EOF

/************************ extra stuff   *************************************/

#ifdef ARTIFACTS

/* A Mennucc: this routine that follows was provided by lvalero and oberger;
it uses  integers for every double variable,
and directly accesses the sinc table; so it is very very fast;
unfortunately if x and dx are integers, then warping is not precise
to the subpixel , and there are visible artifacts

I have rewritten resample_array_inv_conv (above)
so that it directly accesses the sinc_table , and it uses integers,
but so that it is precise to subpixel

in my tests I could achieve these results

old routine     9.55 sec
new routine     3.50 sec
lvalero,oberger 2.86 sec 

note that this routine needs the old style sinc table which is not provided
any more 
*/
#include "sinc_table.h"




/* this test was provided by lvalero, oberger 05/05/2004 :*/

/*With one image 2048*2048 in 16 bits, Athlon 1.4Ghz + parhelia 128 Mo + 785 Mo RAM
 we have :
for resample_array_inv_conv (optimised) 3s
for resample_array_inv_bilin 2.5s
for resample_array_inv_noaa_ 1.4s*/     



//VERSION OPTIMISEE MODIF lvalero, oberger 05/05/2004
static void resample_array_inv_conv (const double *F, const PIXEL_TYPE *src, int s_len, int s_stride, PIXEL_TYPE *dst, int d_len, int d_stride)
{
  int i, p = 0, j, lastj ;
  int x, nx, px, dx ;
  double c, s, v ;
  int firstj,index_stride,index,increment ;
  
  for( i = 0 ; i < d_len ; ++i ) 
  {
   v=0; c=0;
   x=(int)F[i];
   nx=(int)F[i+1];
   px=(int)F[i-1];
   //dx=ABS(nx-px)/2.;
   dx=ABS(nx-px) >> 1 ;

   if(dx<1) dx=1;

   //lastj=ceil(x+dx+dx);
   lastj = x+dx+dx ;
   
   // optimisation pour v+=s* (double)src[j*s_stride]; => on economise une multiplication
   // par defaut, g met floor pour passer en int, mais ce peut etre mauvais
   //firstj=floor(x-dx-dx);
   firstj = x-dx-dx;

   //index_stride = floor(firstj*s_stride) ;
   index_stride = firstj*s_stride ;
   // en deroulant la boucle avec j = x-dx-dx on a :
   //
   // (x-j)/dx = (x-x+dx+dx)/dx = (2*dx)/dx = 2  <- etape 0
   // (x-(j+1))/dx = (x-(x+1)+dx+dx)/dx  = 2-(1/dx)  <- etape 1
   // (x-(j+2))/dx = (x-(x+2)+dx+dx)/dx  = 2-(2/dx)  <- etape 2
   // (x-(j+3))/dx = (x-(x+3)+dx+dx)/dx  = 2-(3/dx)  <- etape 3
   // ...
   // (x-lastj)/dx = (x-x-dx-dx)/dx  = -2  <- etape fin
   //
   // Donc on a pour l'etape n :
   //
   // ((2 - (n/dx)) + 4)*100 = (2+4)*100 - (n/dx)*100
   //           = (2+4)*100 - n*(1/dx)*100
   //           = (2+4)*100 - n*(100/dx)
   
   // avec le commentaire precedent, on s'apercoit que la premiere valeur est :
   //index = (2+4)*100 ;
   index = 600 ;
   // et l'increment devient :
   //increment = floor(100./dx) ;//  correspond a (1/dx)*100, mais on racle une multiplication
   increment = 100/dx ;
   
   for ( j = firstj ; j <= lastj  ; ++j )
   {
    // avec cette optimisation il convient de changer le test : laisser les operations
    // comme elle sont ecrites, le compilo fera le calcul
    //if ( index < ((-2+4)*100) /*|| index > ((2+4)*100)*/ )
    if ( index < (200) /*|| index > ((2+4)*100)*/ )
     s = 0 ;    
    else
     s=sinc_table [index];
    if ( j>=0 && j<s_len)
     v+=s* (double)src[index_stride];
     c+=s;
    
    index -= increment ;
    ++index_stride ;
   }
   if (!(c<0.0001 && c > -0.0001)) v /=c ;
   dst[p]=(PIXEL_TYPE)CLAMP(v,0,255);
   p+= d_stride;    
  }
}



#endif