/*
* 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() */