/*
 				interpolate.c

*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*
*	Part of:	SWarp
*
*	Author:		E.BERTIN (IAP)
*
*	Contents:	Function related to vignet manipulations.
*
*	Last modify:	14/02/2001
*
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
*/

#include	<math.h>
#include	<stdio.h>
#include	<stdlib.h>

#include	"define.h"
#include	"types.h"
#include	"globals.h"
#include	"interpolate.h"

int		interp_kernwidth[5]={1,2,4,6,8};

PIXTYPE		*kernel_buffer, *wkernel_buffer;
interpenum	kernel_interptype[INTERP_MAXDIM];
int		kernel_width[INTERP_MAXDIM],
		kernel_nlines;

/****** interpolate_pix *******************************************************
PROTO	int interpolate_pix(fieldstruct *field, fieldstruct *wfield,
		double *pos, PIXTYPE *pixout, PIXTYPE *wpixout)
PURPOSE	Interpolate pixel data through sinc interpolation.
INPUT	Field structure pointer,
	Weight field structure pointer,
	Position vector,
	Pointer to the output pixel,
	Pointer to the output weight.
OUTPUT	RETURN_OK if pixel falls within the input frame, RETURN_ERROR
	otherwise.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	14/02/2001
 ***/
int	interpolate_pix(fieldstruct *field, fieldstruct *wfield,
		double *pos, PIXTYPE *outpix, PIXTYPE *woutpix)

  {
   PIXTYPE		*pixin,*pixout;
   static double	dpos[INTERP_MAXDIM],
			kernel_vector[INTERP_MAXKERNELWIDTH];
   double		*kvector,
			max, val;
   static long		step[INTERP_MAXDIM];
   long			start, fac;
   static int		ipos[INTERP_MAXDIM], linecount[INTERP_MAXDIM];
   int			*naxisn,
			i,j,n, ival, naxis, nlines, kwidth,width;

  naxis = field->tab->naxis;
  naxisn = field->tab->naxisn;
  width = field->width;
  start = 0;
  fac = 1;
  for (n=0; n<naxis; n++)
    {
    val = *(pos++);
    width = *(naxisn++);
/*-- Get the integer part of the current coordinate or nearest neighbour */
    ival = (kernel_interptype[n]==INTERP_NEARESTNEIGHBOUR)?
					(int)(val-0.50001):(int)val;
/*-- Store the fractional part of the current coordinate */
    dpos[n] = val - ival;
/*-- Check if interpolation start/end exceed image boudary... */
    kwidth = kernel_width[n];
    ipos[n] = (ival-=kwidth/2);
    if (ival<0 || ival+kwidth>width)
      {
      *outpix = 0.0;
      if (woutpix)
        *woutpix = BIG;
      return RETURN_ERROR;
      }
/*-- Update starting pointer */
    start += ival*fac;
/*-- Update step between interpolated regions */
    step[n] = fac*(width-kwidth);
    fac *= width;
    }

/* Update Interpolation kernel vectors */
  make_kernel(*dpos, kernel_vector, *kernel_interptype);
  kwidth = *kernel_width;
  nlines = kernel_nlines;
/* First step: interpolate along NAXIS1 from the data themselves */
  pixin = field->pix+start;
  pixout = kernel_buffer;
  for (j=nlines; j--;)
    {
    val = 0.0;
    kvector = kernel_vector;
    for (i=kwidth; i--;)
      val += *(kvector++)**(pixin++);
    *(pixout++) = val;
    for (n=1; n<naxis; n++)
      {
      pixin+=step[n-1];
      if (++linecount[n]<kernel_width[n])
        break;
      else
        linecount[n] = 0;	/* No need to initialize it to 0! */
      }
    }

/* Now the weight (variance, in fact) map */
  if (wfield)
    {
    pixin = wfield->pix+start;
    pixout = wkernel_buffer;
    for (j=nlines; j--;)
      {
      max = 0.0;
      for (i=kwidth; i--;)
        if ((val = *(pixin++))>max)
          max = val;
      *(pixout++) = max;
      for (n=1; n<naxis; n++)
        {
        pixin+=step[n-1];
        if (++linecount[n]<kernel_width[n])
          break;
        else
          linecount[n] = 0;
        }
      }
    }

/* Second step: interpolate along other axes from the interpolation buffer */
  for (n=1; n<naxis; n++)
    {
    make_kernel(dpos[n], kernel_vector, kernel_interptype[n]);
    kwidth = kernel_width[n];
    pixout = pixin = kernel_buffer;
    for (j = (nlines/=kwidth); j--;)
      {
      val = 0.0;
      kvector = kernel_vector;
      for (i=kwidth; i--;)
        val += *(kvector++)**(pixin++);
      *(pixout++) = val;
     }
/*-- Now the weight (variance, in fact) map */
    if (wfield)
      {
      pixout = pixin = wkernel_buffer;
      for (j = nlines; j--;)
        {
        max = 0.0;
        for (i=kwidth; i--;)
          if ((val = *(pixin++))>max)
            max = val;
        *(pixout++) = max;
        }
      }
    }

/* Finally, fill the output pointer(s) */
  *outpix = *kernel_buffer;
  if (woutpix)
    *woutpix = wfield? *wkernel_buffer : field->backsig*field->backsig;

  return RETURN_OK;
  }


/****** make_kernel **********************************************************
PROTO	void make_kernel(double pos, double *kernel, interpenum interptype)
PURPOSE	Conpute interpolation-kernel data
INPUT	Position,
	Pointer to the output kernel data,
	Interpolation method.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	26/08/2000
 ***/
void	make_kernel(double pos, double *kernel, interpenum interptype)
  {
   double	x, val;

  switch(interptype)
    {
    case INTERP_NEARESTNEIGHBOUR:
      *kernel = 1;
      break;
    case INTERP_BILINEAR:
      *(kernel++) = 1.0-pos;
      *kernel = pos;
      break;
    case INTERP_LANCZOS2:
      if (fabs(pos)<1e-5)
        {
        *(kernel++) = 0.0;
        *(kernel++) = 1.0;
        *(kernel++) = 0.0;
        *kernel = 0.0;
        }
      else
        {
        val = 0.0;
        x = PI*(-pos-1.0);
        val += (*(kernel++) = 2.0*sin(x/2.0)*sin(x)/(x*x));
        x = PI*(-pos);
        val += (*(kernel++) = 2.0*sin(x/2.0)*sin(x)/(x*x));
        x = PI*(1.0-pos);
        val += (*(kernel++) = 2.0*sin(x/2.0)*sin(x)/(x*x));
        x = PI*(2.0-pos);
        val += (*kernel = 2.0*sin(x/2.0)*sin(x)/(x*x));
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *kernel /= val;
        }
      break;
    case INTERP_LANCZOS3:
      if (fabs(pos)<1e-5)
        {
        *(kernel++) = 0.0;
        *(kernel++) = 0.0;
        *(kernel++) = 1.0;
        *(kernel++) = 0.0;
        *(kernel++) = 0.0;
        *kernel = 0.0;
        }
      else
        {
        val = 0.0;
        x = PI*(-pos-2.0);
        val += (*(kernel++) = 3.0*sin(x/3.0)*sin(x)/(x*x));
        x = PI*(-pos-1.0);
        val += (*(kernel++) = 3.0*sin(x/3.0)*sin(x)/(x*x));
        x = PI*(-pos);
        val += (*(kernel++) = 3.0*sin(x/3.0)*sin(x)/(x*x));
        x = PI*(1.0-pos);
        val += (*(kernel++) = 3.0*sin(x/3.0)*sin(x)/(x*x));
        x = PI*(2.0-pos);
        val += (*(kernel++) = 3.0*sin(x/3.0)*sin(x)/(x*x));
        x = PI*(3.0-pos);
        val += (*kernel = 3.0*sin(x/3.0)*sin(x)/(x*x));
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *kernel /= val;
        }
      break;
    case INTERP_LANCZOS4:
      if (fabs(pos)<1e-5)
        {
        *(kernel++) = 0.0;
        *(kernel++) = 0.0;
        *(kernel++) = 0.0;
        *(kernel++) = 1.0;
        *(kernel++) = 0.0;
        *(kernel++) = 0.0;
        *(kernel++) = 0.0;
        *kernel = 0.0;
        }
      else
        {
        val = 0.0;
        x = PI*(-pos-3.0);
        val += (*(kernel++) = 4.0*sin(x/4.0)*sin(x)/(x*x));
        x = PI*(-pos-2.0);
        val +=(*(kernel++) = 4.0*sin(x/4.0)*sin(x)/(x*x));
        x = PI*(-pos-1.0);
        val += (*(kernel++) = 4.0*sin(x/4.0)*sin(x)/(x*x));
        x = PI*(-pos);
        val += (*(kernel++) = 4.0*sin(x/4.0)*sin(x)/(x*x));
        x = PI*(1.0-pos);
        val += (*(kernel++) = 4.0*sin(x/4.0)*sin(x)/(x*x));
        x = PI*(2.0-pos);
        val += (*(kernel++) = 4.0*sin(x/4.0)*sin(x)/(x*x));
        x = PI*(3.0-pos);
        val += (*(kernel++) = 4.0*sin(x/4.0)*sin(x)/(x*x));
        x = PI*(4.0-pos);
        val += (*kernel = 4.0*sin(x/4.0)*sin(x)/(x*x));
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *(kernel--) /= val;
        *kernel /= val;
        }
      break;
    default:
      error(EXIT_FAILURE, "*Internal Error*: Unknown interpolation type in ",
		"make_kernel()");
    }

  return;
  }


/****** init_kernel ***********************************************************
PROTO	void init_kernel(interpenum *interptype, int naxis)
PURPOSE	Prepare interpolation operations.
INPUT	Interpolation type for each axis,
	Number of axes.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	14/01/2000
 ***/
void	init_kernel(interpenum *interptype, int naxis)
  {
   int	n;

  kernel_nlines = 1;
  for (n=0; n<naxis; n++)
    {
    kernel_nlines *= (kernel_width[n] = interp_kernwidth[(int)interptype[n]]);
    kernel_interptype[n] = interptype[n];
    }
  kernel_nlines /= kernel_width[0];
  QMALLOC(kernel_buffer, PIXTYPE, kernel_nlines);
  QMALLOC(wkernel_buffer, PIXTYPE, kernel_nlines);

  return;
  }


/****** free_kernel ***********************************************************
PROTO	void free_kernel(void)
PURPOSE	Free buffers allocated for interpolation operations.
INPUT	-.
OUTPUT	-.
NOTES	-.
AUTHOR	E. Bertin (IAP)
VERSION	14/01/2000
 ***/
void	free_kernel(void)
  {
  free(kernel_buffer);
  free(wkernel_buffer);
  kernel_nlines = 0;

  return;
  }