/*
* 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 <string.h>
#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.
*/
/*--------------------------------------------------------------------------*/
long double ipow(long double x, int p)
{
long 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)
*/
long double
poly2d_compute( int ncoeff, long double *c, long double u, long double *vpow ) {
long 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(long double * data, int nn)
{
unsigned long n,
mmax,
m,
i, j,
istep ;
long double wtemp,
wr,
wpr,
wpi,
wi,
theta;
long 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.
@param samples Number of samples
@param kernel Block of (samples) * long double
@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}
*/
/*--------------------------------------------------------------------------*/
#define hk_gen(x,s) (((tanh(s*(x+0.5))+1)/2)*((tanh(s*(-x+0.5))+1)/2))
void generate_tanh_kernel(long double steep, int samples, long double *kernel)
{
long double * x ;
long double width ;
long double inv_np ;
long double ind ;
int i ;
int np ;
width = (long double)TABSPERPIX / 2.0 ;
np = 32768 ; /* Hardcoded: should never be changed */
inv_np = 1.00 / (long double)np ;
/*
* Generate the kernel expression in Fourier space
* with a correct frequency ordering to allow standard FT
*/
x = (long double *) malloc((2*np+1)*sizeof(long double)) ;
for (i=0 ; i<np/2 ; i++) {
ind = (long 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 = (long 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) ;
/*
* fill in passed-in array
*/
for (i=0 ; i<samples ; i++) {
kernel[i] = 2.0 * width * x[2*i] * inv_np ;
}
free(x) ;
} /* generate_tanh_kernel() */
/*-------------------------------------------------------------------------*/
/**
@name generate_interpolation_kernel
@memo Generate an interpolation kernel to use in this module.
@param kernel_type Type of interpolation kernel.
@param samples Number of samples
@param kernel Block of (samples) * long double
@return true on success
@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 filled-in array of long doubles is ready to use in the various re-sampling
functions in this module.
*/
/*--------------------------------------------------------------------------*/
char
generate_interpolation_kernel(char *kernel_type, int samples, long double *tab)
{
int i ;
long double x ;
long double alpha ;
long double inv_norm ;
if (kernel_type==NULL || !strcmp(kernel_type, "default") || !strcmp(kernel_type, "tanh")) {
generate_tanh_kernel(TANH_STEEPNESS, samples, tab) ;
} else if (!strcmp(kernel_type, "sinc")) {
tab[0] = 1.0 ;
tab[samples-1] = 0.0 ;
for (i=1 ; i<samples ; i++) {
x = (long double)KERNEL_WIDTH * (long double)i/(long double)(samples-1) ;
tab[i] = sinc(x) ;
}
} else if (!strcmp(kernel_type, "sinc2")) {
tab[0] = 1.0 ;
tab[samples-1] = 0.0 ;
for (i=1 ; i<samples ; i++) {
x = 2.0 * (long double)i/(long double)(samples-1) ;
tab[i] = sinc(x) ;
tab[i] *= tab[i] ;
}
} else if (!strcmp(kernel_type, "lanczos")) {
for (i=0 ; i<samples ; i++) {
x = (long double)KERNEL_WIDTH * (long double)i/(long double)(samples-1) ;
if (fabsl(x)<2) {
tab[i] = sinc(x) * sinc(x/2) ;
} else {
tab[i] = 0.00 ;
}
}
} else if (!strcmp(kernel_type, "hamming")) {
alpha = 0.54 ;
inv_norm = 1.00 / (long double)(samples - 1) ;
for (i=0 ; i<samples ; i++) {
x = (long 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")) {
alpha = 0.50 ;
inv_norm = 1.00 / (long double)(samples - 1) ;
for (i=0 ; i<samples ; i++) {
x = (long 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 {
/**
** trapped at perl level, so should never reach here
e_error("unrecognized kernel type [%s]: aborting generation",
kernel_type) ;
**/
return 0;
}
return 1;
} /* generate_interpolation_kernel() */