459 lines
11 KiB
C++
459 lines
11 KiB
C++
/* various interpolation templates
|
|
*/
|
|
|
|
/*
|
|
|
|
This file is part of VIPS.
|
|
|
|
VIPS is free software; you can redistribute it and/or modify
|
|
it under the terms of the GNU Lesser 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 Lesser General Public License for more details.
|
|
|
|
You should have received a copy of the GNU Lesser General Public License
|
|
along with this program; if not, write to the Free Software
|
|
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
|
|
02110-1301 USA
|
|
|
|
*/
|
|
|
|
/*
|
|
|
|
These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
|
|
|
|
*/
|
|
|
|
/*
|
|
* Various casts which assume that the data is already in range. (That
|
|
* is, they are to be used with monotone samplers.)
|
|
*/
|
|
template <typename T> static T inline
|
|
to_fptypes( const double val )
|
|
{
|
|
const T newval = val;
|
|
|
|
return( newval );
|
|
}
|
|
|
|
template <typename T> static T inline
|
|
to_withsign( const double val )
|
|
{
|
|
const int sign_of_val = 2 * ( val >= 0. ) - 1;
|
|
const int rounded_abs_val = .5 + sign_of_val * val;
|
|
const T newval = sign_of_val * rounded_abs_val;
|
|
|
|
return( newval );
|
|
}
|
|
|
|
template <typename T> static T inline
|
|
to_nosign( const double val )
|
|
{
|
|
const T newval = .5 + val;
|
|
|
|
return( newval );
|
|
}
|
|
|
|
/*
|
|
* Various bilinear implementation templates. Note that no clampling
|
|
* is used: There is an assumption that the data is such that
|
|
* over/underflow is not an issue:
|
|
*/
|
|
|
|
/*
|
|
* Bilinear interpolation for float and double types. The first four
|
|
* inputs are weights, the last four are the corresponding pixel
|
|
* values:
|
|
*/
|
|
template <typename T> static T inline
|
|
bilinear_fptypes(
|
|
const double w_times_z,
|
|
const double x_times_z,
|
|
const double w_times_y,
|
|
const double x_times_y,
|
|
const double tre_thr,
|
|
const double tre_thrfou,
|
|
const double trequa_thr,
|
|
const double trequa_thrfou )
|
|
{
|
|
const T newval =
|
|
w_times_z * tre_thr +
|
|
x_times_z * tre_thrfou +
|
|
w_times_y * trequa_thr +
|
|
x_times_y * trequa_thrfou;
|
|
|
|
return( newval );
|
|
}
|
|
|
|
/*
|
|
* Bilinear interpolation for signed integer types:
|
|
*/
|
|
template <typename T> static T inline
|
|
bilinear_withsign(
|
|
const double w_times_z,
|
|
const double x_times_z,
|
|
const double w_times_y,
|
|
const double x_times_y,
|
|
const double tre_thr,
|
|
const double tre_thrfou,
|
|
const double trequa_thr,
|
|
const double trequa_thrfou )
|
|
{
|
|
const double val =
|
|
w_times_z * tre_thr +
|
|
x_times_z * tre_thrfou +
|
|
w_times_y * trequa_thr +
|
|
x_times_y * trequa_thrfou;
|
|
|
|
const int sign_of_val = 2 * ( val >= 0. ) - 1;
|
|
|
|
const int rounded_abs_val = .5 + sign_of_val * val;
|
|
|
|
const T newval = sign_of_val * rounded_abs_val;
|
|
|
|
return( newval );
|
|
}
|
|
|
|
/*
|
|
* Bilinear Interpolation for unsigned integer types:
|
|
*/
|
|
template <typename T> static T inline
|
|
bilinear_nosign(
|
|
const double w_times_z,
|
|
const double x_times_z,
|
|
const double w_times_y,
|
|
const double x_times_y,
|
|
const double tre_thr,
|
|
const double tre_thrfou,
|
|
const double trequa_thr,
|
|
const double trequa_thrfou )
|
|
{
|
|
const T newval =
|
|
w_times_z * tre_thr +
|
|
x_times_z * tre_thrfou +
|
|
w_times_y * trequa_thr +
|
|
x_times_y * trequa_thrfou +
|
|
0.5;
|
|
|
|
return( newval );
|
|
}
|
|
|
|
/*
|
|
* Bicubic (Catmull-Rom) interpolation templates:
|
|
*/
|
|
|
|
static int inline
|
|
unsigned_fixed_round( int v )
|
|
{
|
|
const int round_by = VIPS_INTERPOLATE_SCALE >> 1;
|
|
|
|
return( (v + round_by) >> VIPS_INTERPOLATE_SHIFT );
|
|
}
|
|
|
|
/* Fixed-point integer bicubic, used for 8-bit types.
|
|
*/
|
|
template <typename T> static int inline
|
|
bicubic_unsigned_int(
|
|
const T uno_one, const T uno_two, const T uno_thr, const T uno_fou,
|
|
const T dos_one, const T dos_two, const T dos_thr, const T dos_fou,
|
|
const T tre_one, const T tre_two, const T tre_thr, const T tre_fou,
|
|
const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
|
|
const int* restrict cx, const int* restrict cy )
|
|
{
|
|
const int c0 = cx[0];
|
|
const int c1 = cx[1];
|
|
const int c2 = cx[2];
|
|
const int c3 = cx[3];
|
|
|
|
const int r0 = unsigned_fixed_round(
|
|
c0 * uno_one +
|
|
c1 * uno_two +
|
|
c2 * uno_thr +
|
|
c3 * uno_fou );
|
|
const int r1 = unsigned_fixed_round(
|
|
c0 * dos_one +
|
|
c1 * dos_two +
|
|
c2 * dos_thr +
|
|
c3 * dos_fou );
|
|
const int r2 = unsigned_fixed_round(
|
|
c0 * tre_one +
|
|
c1 * tre_two +
|
|
c2 * tre_thr +
|
|
c3 * tre_fou );
|
|
const int r3 = unsigned_fixed_round(
|
|
c0 * qua_one +
|
|
c1 * qua_two +
|
|
c2 * qua_thr +
|
|
c3 * qua_fou );
|
|
|
|
return( unsigned_fixed_round(
|
|
cy[0] * r0 +
|
|
cy[1] * r1 +
|
|
cy[2] * r2 +
|
|
cy[3] * r3 ) );
|
|
}
|
|
|
|
static int inline
|
|
signed_fixed_round( int v )
|
|
{
|
|
const int sign_of_v = 2 * (v > 0) - 1;
|
|
const int round_by = sign_of_v * (VIPS_INTERPOLATE_SCALE >> 1);
|
|
|
|
return( (v + round_by) >> VIPS_INTERPOLATE_SHIFT );
|
|
}
|
|
|
|
/* Fixed-point integer bicubic, used for 8-bit types.
|
|
*/
|
|
template <typename T> static int inline
|
|
bicubic_signed_int(
|
|
const T uno_one, const T uno_two, const T uno_thr, const T uno_fou,
|
|
const T dos_one, const T dos_two, const T dos_thr, const T dos_fou,
|
|
const T tre_one, const T tre_two, const T tre_thr, const T tre_fou,
|
|
const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
|
|
const int* restrict cx, const int* restrict cy )
|
|
{
|
|
const int c0 = cx[0];
|
|
const int c1 = cx[1];
|
|
const int c2 = cx[2];
|
|
const int c3 = cx[3];
|
|
|
|
const int r0 = signed_fixed_round(
|
|
c0 * uno_one +
|
|
c1 * uno_two +
|
|
c2 * uno_thr +
|
|
c3 * uno_fou );
|
|
const int r1 = signed_fixed_round(
|
|
c0 * dos_one +
|
|
c1 * dos_two +
|
|
c2 * dos_thr +
|
|
c3 * dos_fou );
|
|
const int r2 = signed_fixed_round(
|
|
c0 * tre_one +
|
|
c1 * tre_two +
|
|
c2 * tre_thr +
|
|
c3 * tre_fou );
|
|
const int r3 = signed_fixed_round(
|
|
c0 * qua_one +
|
|
c1 * qua_two +
|
|
c2 * qua_thr +
|
|
c3 * qua_fou );
|
|
|
|
return( signed_fixed_round(
|
|
cy[0] * r0 +
|
|
cy[1] * r1 +
|
|
cy[2] * r2 +
|
|
cy[3] * r3 ) );
|
|
}
|
|
|
|
template <typename T> static T inline
|
|
cubic_float(
|
|
const T one, const T two, const T thr, const T fou,
|
|
const double* restrict cx )
|
|
{
|
|
return( cx[0] * one +
|
|
cx[1] * two +
|
|
cx[2] * thr +
|
|
cx[3] * fou );
|
|
}
|
|
|
|
/* Floating-point bicubic, used for int/float/double types.
|
|
*/
|
|
template <typename T> static T inline
|
|
bicubic_float(
|
|
const T uno_one, const T uno_two, const T uno_thr, const T uno_fou,
|
|
const T dos_one, const T dos_two, const T dos_thr, const T dos_fou,
|
|
const T tre_one, const T tre_two, const T tre_thr, const T tre_fou,
|
|
const T qua_one, const T qua_two, const T qua_thr, const T qua_fou,
|
|
const double* restrict cx, const double* restrict cy )
|
|
{
|
|
const double r0 = cubic_float<T>(
|
|
uno_one, uno_two, uno_thr, uno_fou, cx );
|
|
const double r1 = cubic_float<T>(
|
|
dos_one, dos_two, dos_thr, dos_fou, cx );
|
|
const double r2 = cubic_float<T>(
|
|
tre_one, tre_two, tre_thr, tre_fou, cx );
|
|
const double r3 = cubic_float<T>(
|
|
qua_one, qua_two, qua_thr, qua_fou, cx );
|
|
|
|
return( cubic_float<T>( r0, r1, r2, r3, cy ) );
|
|
}
|
|
|
|
/* Given an offset in [0,1] (we can have x == 1 when building tables),
|
|
* calculate c0, c1, c2, c3, the catmull-rom coefficients. This is called
|
|
* from the interpolator as well as from the table builder.
|
|
*/
|
|
static void inline
|
|
calculate_coefficients_catmull( double c[4], const double x )
|
|
{
|
|
/* Nicolas believes that the following is an hitherto unknown
|
|
* hyper-efficient method of computing Catmull-Rom coefficients. It
|
|
* only uses 4* & 1+ & 5- for a total of only 10 flops to compute
|
|
* four coefficients.
|
|
*/
|
|
const double cr1 = 1. - x;
|
|
const double cr2 = -.5 * x;
|
|
const double cr3 = cr1 * cr2;
|
|
const double cone = cr1 * cr3;
|
|
const double cfou = x * cr3;
|
|
const double cr4 = cfou - cone;
|
|
const double ctwo = cr1 - cone + cr4;
|
|
const double cthr = x - cfou - cr4;
|
|
|
|
g_assert( x >= 0. && x <= 1. );
|
|
|
|
c[0] = cone;
|
|
c[3] = cfou;
|
|
c[1] = ctwo;
|
|
c[2] = cthr;
|
|
}
|
|
|
|
/* Given an x in [0,1] (we can have x == 1 when building tables),
|
|
* calculate c0 .. c(@shrink + 1), the triangle coefficients. This is called
|
|
* from the interpolator as well as from the table builder.
|
|
*/
|
|
static void inline
|
|
calculate_coefficients_triangle( double *c,
|
|
const double shrink, const double x )
|
|
{
|
|
/* Needs to be in sync with vips_reduce_get_points().
|
|
*/
|
|
const int n_points = 2 * rint( shrink ) + 1;
|
|
const double half = x + n_points / 2.0 - 1;
|
|
|
|
int i;
|
|
double sum;
|
|
|
|
sum = 0;
|
|
for( i = 0; i < n_points; i++ ) {
|
|
const double xp = (i - half) / shrink;
|
|
|
|
double l;
|
|
|
|
l = 1.0 - VIPS_FABS( xp );
|
|
l = VIPS_MAX( 0.0, l );
|
|
|
|
c[i] = l;
|
|
sum += l;
|
|
}
|
|
|
|
for( i = 0; i < n_points; i++ )
|
|
c[i] /= sum;
|
|
}
|
|
|
|
/* Generate a cubic filter. See:
|
|
*
|
|
* Mitchell and Netravali, Reconstruction Filters in Computer Graphics
|
|
* Computer Graphics, Volume 22, Number 4, August 1988.
|
|
*
|
|
* B = 1, C = 0 - cubic B-spline
|
|
* B = 1/3, C = 1/3 - Mitchell
|
|
* B = 0, C = 1/2 - Catmull-Rom spline
|
|
*/
|
|
static void inline
|
|
calculate_coefficients_cubic( double *c,
|
|
const double shrink, const double x, double B, double C )
|
|
{
|
|
/* Needs to be in sync with vips_reduce_get_points().
|
|
*/
|
|
const int n_points = 2 * rint( 2 * shrink ) + 1;
|
|
const double half = x + n_points / 2.0 - 1;
|
|
|
|
int i;
|
|
double sum;
|
|
|
|
sum = 0;
|
|
for( i = 0; i < n_points; i++ ) {
|
|
const double xp = (i - half) / shrink;
|
|
const double axp = VIPS_FABS( xp );
|
|
const double axp2 = axp * axp;
|
|
const double axp3 = axp2 * axp;
|
|
|
|
double l;
|
|
|
|
if( axp <= 1 )
|
|
l = ((12 - 9 * B - 6 * C) * axp3 +
|
|
(-18 + 12 * B + 6 * C) * axp2 +
|
|
(6 - 2 * B)) / 6;
|
|
else if( axp <= 2 )
|
|
l = ((-B - 6 * C) * axp3 +
|
|
(6 * B + 30 * C) * axp2 +
|
|
(-12 * B - 48 * C) * axp +
|
|
(8 * B + 24 * C)) / 6;
|
|
else
|
|
l = 0.0;
|
|
|
|
c[i] = l;
|
|
sum += l;
|
|
}
|
|
|
|
for( i = 0; i < n_points; i++ )
|
|
c[i] /= sum;
|
|
}
|
|
|
|
/* Given an x in [0,1] (we can have x == 1 when building tables),
|
|
* calculate c0 .. c(@a * @shrink + 1), the lanczos coefficients. This is called
|
|
* from the interpolator as well as from the table builder.
|
|
*
|
|
* @a is the number of lobes, so usually 2 or 3. @shrink is the reduction
|
|
* factor, so 1 for interpolation, 2 for a x2 reduction, etc. We need more
|
|
* points for large decimations to avoid aliasing.
|
|
*/
|
|
static void inline
|
|
calculate_coefficients_lanczos( double *c,
|
|
const int a, const double shrink, const double x )
|
|
{
|
|
/* Needs to be in sync with vips_reduce_get_points().
|
|
*/
|
|
const int n_points = 2 * rint( a * shrink ) + 1;
|
|
const double half = x + n_points / 2.0 - 1;
|
|
|
|
int i;
|
|
double sum;
|
|
|
|
sum = 0;
|
|
for( i = 0; i < n_points; i++ ) {
|
|
const double xp = (i - half) / shrink;
|
|
|
|
double l;
|
|
|
|
if( xp == 0.0 )
|
|
l = 1.0;
|
|
else if( xp < -a )
|
|
l = 0.0;
|
|
else if( xp > a )
|
|
l = 0.0;
|
|
else
|
|
l = (double) a * sin( VIPS_PI * xp ) *
|
|
sin( VIPS_PI * xp / (double) a ) /
|
|
(VIPS_PI * VIPS_PI * xp * xp);
|
|
|
|
c[i] = l;
|
|
sum += l;
|
|
}
|
|
|
|
for( i = 0; i < n_points; i++ )
|
|
c[i] /= sum;
|
|
}
|
|
|
|
/* Our inner loop for resampling with a convolution. Operate on elements of
|
|
* type T, gather results in an intermediate of type IT.
|
|
*/
|
|
template <typename T, typename IT>
|
|
static IT
|
|
reduce_sum( const T * restrict in, int stride, const IT * restrict c, int n )
|
|
{
|
|
IT sum;
|
|
|
|
sum = 0;
|
|
for( int i = 0; i < n; i++ ) {
|
|
sum += c[i] * in[0];
|
|
in += stride;
|
|
}
|
|
|
|
return( sum );
|
|
}
|