/*
* resample.c
* - a hacked version of
* ipow.c, poly2d.c, and resampling.c
* from version 3.6-0 of the Eclipse library (ESO)
* by Nicolas Devillard
*
* see http://www.eso.org/eclipse for further details
*/
#include "resample.h"
/*-------------------------------------------------------------------------*/
/**
@name ipow
@memo Same as pow(x,y) but for integer values of y only (faster).
@param x A double number.
@param p An integer power.
@return x to the power p.
@doc
This is much faster than the math function due to the integer. Some
compilers make this optimization already, some do not.
p can be positive, negative or null.
*/
/*--------------------------------------------------------------------------*/
double ipow(double x, int p)
{
double r, recip ;
/* Get rid of trivial cases */
switch (p) {
case 0:
return 1.00 ;
case 1:
return x ;
case 2:
return x*x ;
case 3:
return x*x*x ;
case -1:
return 1.00 / x ;
case -2:
return (1.00 / x) * (1.00 / x) ;
}
if (p>0) {
r = x ;
while (--p) r *= x ;
} else {
r = recip = 1.00 / x ;
while (++p) r *= recip ;
}
return r;
}
/*
* compute the value of a 2D polynomial at a point
* - it assumes that ncoeff is a small number, so there's
* little point in pre-calculating ipow(u,i)
*/
double
poly2d_compute( int ncoeff, double *c, double u, double *vpow ) {
double out;
int i, j, k;
out = 0.00;
k = 0;
for( j = 0; j < ncoeff; j++ ) {
for( i = 0; i < ncoeff; i++ ) {
out += c[k] * ipow( u, i ) * vpow[j];
k++;
}
}
return out;
}
/*-------------------------------------------------------------------------*/
/**
@name sinc
@memo Cardinal sine.
@param x double value.
@return 1 double.
@doc
Compute the value of the function sinc(x)=sin(pi*x)/(pi*x) at the
requested x.
*/
/*--------------------------------------------------------------------------*/
double
sinc(double x)
{
if (fabs(x)<1e-4)
return (double)1.00 ;
else
return ((sin(x * (double)PI_NUMB)) / (x * (double)PI_NUMB)) ;
} /* sinc() */
/**
@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,
istep ;
double wtemp,
wr,
wpr,
wpi,
wi,
theta;
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) {
istep = mmax << 1;
theta = 2 * M_PI / mmax;
wtemp = sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = sin(theta);
wr = 1.0;
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;
}
} /* reverse_tanh_kernel() */
#undef KERNEL_SW
/*-------------------------------------------------------------------------*/
/**
@name 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 free().
*/
/*--------------------------------------------------------------------------*/
#define hk_gen(x,s) (((tanh(s*(x+0.5))+1)/2)*((tanh(s*(-x+0.5))+1)/2))
double * 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 = (double *) 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 = (double *) malloc(samples * sizeof(double)) ;
for (i=0 ; i<samples ; i++) {
kernel[i] = 2.0 * width * x[2*i] * inv_np ;
}
free(x) ;
return kernel ;
} /* generate_tanh_kernel() */
/*-------------------------------------------------------------------------*/
/**
@name 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 free().
*/
/*--------------------------------------------------------------------------*/
double *
generate_interpolation_kernel(char * kernel_type)
{
double * tab ;
int i ;
double x ;
double alpha ;
double inv_norm ;
int samples = KERNEL_SAMPLES ;
if (kernel_type==NULL) {
tab = generate_interpolation_kernel("tanh") ;
} else if (!strcmp(kernel_type, "default")) {
tab = generate_interpolation_kernel("tanh") ;
} else if (!strcmp(kernel_type, "sinc")) {
tab = (double *) 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] = sinc(x) ;
}
} else if (!strcmp(kernel_type, "sinc2")) {
tab = (double *) 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] = sinc(x) ;
tab[i] *= tab[i] ;
}
} else if (!strcmp(kernel_type, "lanczos")) {
tab = (double *) 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] = sinc(x) * sinc(x/2) ;
} else {
tab[i] = 0.00 ;
}
}
} else if (!strcmp(kernel_type, "hamming")) {
tab = (double *) 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 = (double *) 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 = generate_tanh_kernel(TANH_STEEPNESS) ;
} else {
/**
** trapped at perl level, so should never reach here
e_error("unrecognized kernel type [%s]: aborting generation",
kernel_type) ;
**/
return NULL ;
}
return tab ;
} /* generate_interpolation_kernel() */
/***********
*** END ***
***********/