/*
 * This file is part of the ESO SINFONI Pipeline
 * Copyright (C) 2004,2005 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */
/*----------------------------------------------------------------------------

   File name     :    resampling.c
   Author         :    Nicolas Devillard
   Created on    :    Jan 04, 1996
   Description    :    resampling routines

 ---------------------------------------------------------------------------*/

/*
    $Id: sinfo_resampling.c,v 1.6 2012-03-03 10:18:26 amodigli Exp $
    $Author: amodigli $
    $Date: 2012-03-03 10:18:26 $
    $Revision: 1.6 $
 */

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif
/*---------------------------------------------------------------------------
                                  Includes
 ---------------------------------------------------------------------------*/
#include <string.h>
#include "sinfo_resampling.h"
#include <math.h>
#include "sinfo_globals.h"
/*---------------------------------------------------------------------------
                              Private functions
 ---------------------------------------------------------------------------*/
static void reverse_tanh_kernel(double * data, int nn) ;
static double * sinfo_generate_tanh_kernel(double steep);

/**@{*/
/**
 * @addtogroup sinfo_rec_jitter Image resampling
 *
 * TBD
 */



/*---------------------------------------------------------------------------
                              Function codes
 ---------------------------------------------------------------------------*/

/**
  @name        sinfo_generate_interpolation_kernel
  @memo        Generate an interpolation kernel to use in this module.
  @param    kernel_type        Type of interpolation kernel.
  @return    1 newly allocated array of doubles.
  @doc

  Provide the name of the kernel you want to generate. Supported kernel
  types are:

  \begin{tabular}{ll}
  NULL            &    default kernel, currently "tanh" \\
  "default"        &    default kernel, currently "tanh" \\
  "tanh"        &    Hyperbolic tangent \\
  "sinc2"        &    Square sinc \\
  "lanczos"        &    Lanczos2 kernel \\
  "hamming"        &    Hamming kernel \\
  "hann"        &    Hann kernel
  \end{tabular}

  The returned array of doubles is ready of use in the various re-sampling
  functions in this module. It must be deallocated using cpl_free().
 */

double   *
sinfo_generate_interpolation_kernel(const char * kernel_type)
{
    double  *    tab ;
    int         i ;
    double      x ;
    double        alpha ;
    double        inv_norm ;
    int         samples = KERNEL_SAMPLES ;

    if (kernel_type==NULL) {
        tab = sinfo_generate_interpolation_kernel("tanh") ;
    } else if (!strcmp(kernel_type, "default")) {
        tab = sinfo_generate_interpolation_kernel("tanh") ;
    } else if (!strcmp(kernel_type, "sinc")) {
        tab = cpl_malloc(samples * sizeof(double)) ;
        tab[0] = 1.0 ;
        tab[samples-1] = 0.0 ;
        for (i=1 ; i<samples ; i++) {
            x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
            tab[i] = sinfo_sinc(x) ;
        }
    } else if (!strcmp(kernel_type, "sinc2")) {
        tab = cpl_malloc(samples * sizeof(double)) ;
        tab[0] = 1.0 ;
        tab[samples-1] = 0.0 ;
        for (i=1 ; i<samples ; i++) {
            x = 2.0 * (double)i/(double)(samples-1) ;
            tab[i] = sinfo_sinc(x) ;
            tab[i] *= tab[i] ;
        }
    } else if (!strcmp(kernel_type, "lanczos")) {
        tab = cpl_malloc(samples * sizeof(double)) ;
        for (i=0 ; i<samples ; i++) {
            x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
            if (fabs(x)<2) {
                tab[i] = sinfo_sinc(x) * sinfo_sinc(x/2) ;
            } else {
                tab[i] = 0.00 ;
            }
        }
    } else if (!strcmp(kernel_type, "hamming")) {
        tab = cpl_malloc(samples * sizeof(double)) ;
        alpha = 0.54 ;
        inv_norm  = 1.00 / (double)(samples - 1) ;
        for (i=0 ; i<samples ; i++) {
            x = (double)i ;
            if (i<(samples-1)/2) {
                tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
            } else {
                tab[i] = 0.0 ;
            }
        }
    } else if (!strcmp(kernel_type, "hann")) {
        tab = cpl_malloc(samples * sizeof(double)) ;
        alpha = 0.50 ;
        inv_norm  = 1.00 / (double)(samples - 1) ;
        for (i=0 ; i<samples ; i++) {
            x = (double)i ;
            if (i<(samples-1)/2) {
                tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
            } else {
                tab[i] = 0.0 ;
            }
        }
    } else if (!strcmp(kernel_type, "tanh")) {
        tab = sinfo_generate_tanh_kernel(TANH_STEEPNESS) ;
    } else {
        sinfo_msg_error("unrecognized kernel type [%s]: aborting generation",
                        kernel_type) ;
        return NULL ;
    }

    return tab ;
}

/**
  @name        sinfo_sinc
  @memo        Cardinal sine.
  @param    x    double value.
  @return    1 double.
  @doc

  Compute the value of the function sinfo_sinc(x)=sin(pi*x)/(pi*x) at the
  requested x.
 */

double
sinfo_sinc(double x)
{
    if (fabs(x)<1e-4)
        return (double)1.00 ;
    else
        return ((sin(x * (double)PI_NUMB)) / (x * (double)PI_NUMB)) ;
}

/**
  @name        sinfo_generate_tanh_kernel
  @memo        Generate a hyperbolic tangent kernel.
  @param    steep    Steepness of the hyperbolic tangent parts.
  @return    1 pointer to a newly allocated array of doubles.
  @doc

  The following function builds up a good approximation of a box filter. It
  is built from a product of hyperbolic tangents. It has the following
  properties:

  \begin{itemize}
  \item It converges very quickly towards +/- 1.
  \item The converging transition is very sharp.
  \item It is infinitely differentiable everywhere (i.e. smooth).
  \item The transition sharpness is scalable.
  \end{itemize}

  The returned array must be deallocated using cpl_free().
 */

#define hk_gen(x,s) (((tanh(s*(x+0.5))+1)/2)*((tanh(s*(-x+0.5))+1)/2))

static double * 
sinfo_generate_tanh_kernel(double steep)
{
    double  *   kernel ;
    double  *   x ;
    double      width ;
    double      inv_np ;
    double      ind ;
    int         i ;
    int         np ;
    int         samples ;

    width   = (double)TABSPERPIX / 2.0 ; 
    samples = KERNEL_SAMPLES ;
    np      = 32768 ; /* Hardcoded: should never be changed */
    inv_np  = 1.00 / (double)np ;

    /*
     * Generate the kernel expression in Fourier space
     * with a correct frequency ordering to allow standard FT
     */
    x = cpl_malloc((2*np+1)*sizeof(double)) ;
    for (i=0 ; i<np/2 ; i++) {
        ind      = (double)i * 2.0 * width * inv_np ;
        x[2*i]   = hk_gen(ind, steep) ;
        x[2*i+1] = 0.00 ;
    }
    for (i=np/2 ; i<np ; i++) {
        ind      = (double)(i-np) * 2.0 * width * inv_np ;
        x[2*i]   = hk_gen(ind, steep) ;
        x[2*i+1] = 0.00 ;
    }

    /* 
     * Reverse Fourier to come back to image space
     */
    reverse_tanh_kernel(x, np) ;

    /*
     * Allocate and fill in returned array
     */
    kernel = cpl_malloc(samples * sizeof(double)) ;
    for (i=0 ; i<samples ; i++) {
        kernel[i] = 2.0 * width * x[2*i] * inv_np ;
    }
    cpl_free(x) ;
    return kernel ;
}

/**
  @name        reverse_tanh_kernel
  @memo        Bring a hyperbolic tangent kernel from Fourier to normal space.
  @param    data    Kernel samples in Fourier space.
  @param    nn        Number of samples in the input kernel.
  @return    void
  @doc

  Bring back a hyperbolic tangent kernel from Fourier to normal space. Do
  not try to understand the implementation and DO NOT MODIFY THIS FUNCTION.
 */

#define KERNEL_SW(a,b) tempr=(a);(a)=(b);(b)=tempr
static void 
reverse_tanh_kernel(double * data, int nn)
{
    unsigned long   n,
    mmax,
    m,
    i, j;

    double  tempr,
    tempi;

    n = (unsigned long)nn << 1;
    j = 1;
    for (i=1 ; i<n ; i+=2) {
        if (j > i) {
            KERNEL_SW(data[j-1],data[i-1]);
            KERNEL_SW(data[j],data[i]);
        }
        m = n >> 1;
        while (m>=2 && j>m) {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    mmax = 2;
    while (n > mmax) {
        unsigned long istep = mmax << 1;
        double theta = 2 * PI_NUMB / mmax;
        double wtemp = sin(0.5 * theta);
        double wpr = -2.0 * wtemp * wtemp;
        double wpi = sin(theta);
        double wr  = 1.0;
        double wi  = 0.0;
        for (m=1 ; m<mmax ; m+=2) {
            for (i=m ; i<=n ; i+=istep) {
                j = i + mmax;
                tempr = wr * data[j-1] - wi * data[j];
                tempi = wr * data[j]   + wi * data[j-1];
                data[j-1] = data[i-1] - tempr;
                data[j]   = data[i]   - tempi;
                data[i-1] += tempr;
                data[i]   += tempi;
            }
            wr = (wtemp = wr) * wpr - wi * wpi + wr;
            wi = wi * wpr + wtemp * wpi + wi;
        }
        mmax = istep;
    }
}
#undef KERNEL_SW

/**
  @name        sinfo_invert_linear_transform
  @memo        Invert a linear transformation.
  @param    trans    Transformation to sinfo_invert.
  @return    1 newly allocated array of 6 doubles.
  @doc

  Given 6 parameters a, b, c, d, e, f defining a linear transform such as:
  \begin{verbatim}
  u = ax + by + c
  v = dx + ey + f
  \end{verbatim}

  The inverse transform is also linear, and is defined by:
  \begin{verbatim}
  x = Au + Bv + C
  y = Du + Ev + F
  \end{verbatim}

  where:

  \begin{verbatim}
  if G = (ae-bd)

  A =  e/G
  B = -b/G
  C = (bf-ce)/G
  D = -d/G
  E =  a/G
  F = (cd-af)/G
  \end{verbatim}

  Notice that if G=0 (ae=bd) the transform cannot be reversed.

  The returned array must be deallocated using cpl_free().
 */

double *
sinfo_invert_linear_transform(double *trans)
{
    double   *    i_trans ;
    double       det ;

    if (trans==NULL) return NULL ;
    det = (trans[0] * trans[4]) - (trans[1] * trans[3]) ;
    if (fabs(det) < 1e-6) {
        sinfo_msg_error("NULL determinant: cannot sinfo_invert transform") ;
        return NULL ;
    }

    i_trans = cpl_calloc(6, sizeof(double)) ;

    i_trans[0] =  trans[4] / det ; 
    i_trans[1] = -trans[1] / det ;
    i_trans[2] = ((trans[1] * trans[5]) - (trans[2] * trans[4])) / det ;
    i_trans[3] = -trans[3] / det ;
    i_trans[4] =  trans[0] / det ;
    i_trans[5] = ((trans[2] * trans[3]) - (trans[0] * trans[5])) / det ;

    return i_trans ;
}
/**@}*/