interpolation precision hacking

This commit is contained in:
John Cupitt 2010-05-28 21:16:28 +00:00
parent 53cc08881d
commit d0573100b3
5 changed files with 151 additions and 79 deletions

88
TODO
View File

@ -1,4 +1,92 @@
mask selection seems reasonable:
vips_interpolate_bicubic_interpolate: 2 2
left=1, top=1, width=4, height=4
maskx=0, masky=0
vips_interpolate_bicubic_interpolate: 2.1 2
left=1, top=1, width=4, height=4
maskx=3, masky=0
vips_interpolate_bicubic_interpolate: 2.2 2
left=1, top=1, width=4, height=4
maskx=6, masky=0
vips_interpolate_bicubic_interpolate: 2.3 2
left=1, top=1, width=4, height=4
maskx=9, masky=0
vips_interpolate_bicubic_interpolate: 2.4 2
left=1, top=1, width=4, height=4
maskx=12, masky=0
vips_interpolate_bicubic_interpolate: 2.5 2
left=1, top=1, width=4, height=4
maskx=16, masky=0
vips_interpolate_bicubic_interpolate: 2.6 2
left=1, top=1, width=4, height=4
maskx=19, masky=0
vips_interpolate_bicubic_interpolate: 2.7 2
left=1, top=1, width=4, height=4
maskx=22, masky=0
vips_interpolate_bicubic_interpolate: 2.8 2
left=1, top=1, width=4, height=4
maskx=25, masky=0
vips_interpolate_bicubic_interpolate: 2.9 2
left=1, top=1, width=4, height=4
maskx=28, masky=0
vips_interpolate_bicubic_interpolate: 3 2
left=2, top=1, width=4, height=4
maskx=0, masky=0
so for (eg.) 2.4 we pick mask 12, which is (12 / 32), or 0.375
however, mask13 would be closer: 13 / 32) = 0.4076
calculation is
2.4 * 32 == 76.8
FLOOR(76.8) == 76
76 & 31 == 12
vips_interpolate_bicubic_class_init:
mask = 0, calculate_coefficients_catmull: 0
mask = 1, calculate_coefficients_catmull: 0.03125
mask = 2, calculate_coefficients_catmull: 0.0625
mask = 3, calculate_coefficients_catmull: 0.09375
mask = 4, calculate_coefficients_catmull: 0.125
mask = 5, calculate_coefficients_catmull: 0.15625
mask = 6, calculate_coefficients_catmull: 0.1875
mask = 7, calculate_coefficients_catmull: 0.21875
mask = 8, calculate_coefficients_catmull: 0.25
mask = 9, calculate_coefficients_catmull: 0.28125
mask = 10, calculate_coefficients_catmull: 0.3125
mask = 11, calculate_coefficients_catmull: 0.34375
mask = 12, calculate_coefficients_catmull: 0.375
mask = 13, calculate_coefficients_catmull: 0.40625
mask = 14, calculate_coefficients_catmull: 0.4375
mask = 15, calculate_coefficients_catmull: 0.46875
mask = 16, calculate_coefficients_catmull: 0.5
mask = 17, calculate_coefficients_catmull: 0.53125
mask = 18, calculate_coefficients_catmull: 0.5625
mask = 19, calculate_coefficients_catmull: 0.59375
mask = 20, calculate_coefficients_catmull: 0.625
mask = 21, calculate_coefficients_catmull: 0.65625
mask = 22, calculate_coefficients_catmull: 0.6875
mask = 23, calculate_coefficients_catmull: 0.71875
mask = 24, calculate_coefficients_catmull: 0.75
mask = 25, calculate_coefficients_catmull: 0.78125
mask = 26, calculate_coefficients_catmull: 0.8125
mask = 27, calculate_coefficients_catmull: 0.84375
mask = 28, calculate_coefficients_catmull: 0.875
mask = 29, calculate_coefficients_catmull: 0.90625
mask = 30, calculate_coefficients_catmull: 0.9375
mask = 31, calculate_coefficients_catmull: 0.96875
mask = 32, calculate_coefficients_catmull: 1
- tools subdirs are now pretty stupid :-( just have a single dir
- int bicubic needs more precision? try a 10x upscale and compare to the

View File

@ -96,7 +96,7 @@ int vips_interpolate_get_window_offset( VipsInterpolate *interpolate );
/* How many bits of precision we keep for transformations, ie. how many
* pre-computed matricies we have.
*/
#define VIPS_TRANSFORM_SHIFT (5)
#define VIPS_TRANSFORM_SHIFT (2)
#define VIPS_TRANSFORM_SCALE (1 << VIPS_TRANSFORM_SHIFT)
/* How many bits of precision we keep for interpolation, ie. where the decimal

View File

@ -68,26 +68,21 @@
(G_TYPE_INSTANCE_GET_CLASS( (obj), \
VIPS_TYPE_INTERPOLATE_BICUBIC, VipsInterpolateBicubicClass ))
typedef struct _VipsInterpolateBicubic {
VipsInterpolate parent_object;
typedef VipsInterpolate VipsInterpolateBicubic;
} VipsInterpolateBicubic;
typedef VipsInterpolateClass VipsInterpolateBicubicClass;
typedef struct _VipsInterpolateBicubicClass {
VipsInterpolateClass parent_class;
/* Precalculated interpolation matrices. int (used for pel
* sizes up to short), and double (for all others). We go to
* scale + 1 so we can round-to-nearest safely.
*/
/* Precalculated interpolation matrices. int (used for pel
* sizes up to short), and double (for all others). We go to
* scale + 1, so we can round-to-nearest safely.
*/
/* We could keep a large set of 2d 4x4 matricies, but this actually
* works out slower, since for many resizes the thing will no longer
* fit in L1.
*/
int matrixi[VIPS_TRANSFORM_SCALE + 1][4];
double matrixf[VIPS_TRANSFORM_SCALE + 1][4];
} VipsInterpolateBicubicClass;
/* We could keep a large set of 2d 4x4 matricies, but this actually
* works out slower since for many resizes the thing will no longer
* fit in L1.
*/
static int vips_bicubic_matrixi[VIPS_TRANSFORM_SCALE + 1][4];
static double vips_bicubic_matrixf[VIPS_TRANSFORM_SCALE + 1][4];
/* We need C linkage for this.
*/
@ -307,15 +302,15 @@ static void
vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateBicubicClass *bicubic_class =
VIPS_INTERPOLATE_BICUBIC_GET_CLASS( interpolate );
/* Scaled int.
*/
const double sx = x * VIPS_TRANSFORM_SCALE;
const double sy = y * VIPS_TRANSFORM_SCALE;
const int sxi = FAST_PSEUDO_FLOOR( sx );
const int syi = FAST_PSEUDO_FLOOR( sy );
/* We know sx/sy are always positive, so we can just (int) them.
*/
const int sxi = (int) sx;
const int syi = (int) sy;
/* Get index into interpolation table and unscaled integer
* position.
@ -327,10 +322,10 @@ vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
/* Look up the tables we need.
*/
const int *cxi = bicubic_class->matrixi[tx];
const int *cyi = bicubic_class->matrixi[ty];
const double *cxf = bicubic_class->matrixf[tx];
const double *cyf = bicubic_class->matrixf[ty];
const int *cxi = vips_bicubic_matrixi[tx];
const int *cyi = vips_bicubic_matrixi[ty];
const double *cxf = vips_bicubic_matrixf[tx];
const double *cyf = vips_bicubic_matrixf[ty];
/* Back and up one to get the top-left of the 4x4.
*/
@ -345,6 +340,7 @@ vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
printf( "vips_interpolate_bicubic_interpolate: %g %g\n", x, y );
printf( "\tleft=%d, top=%d, width=%d, height=%d\n",
ix - 1, iy - 1, 4, 4 );
printf( "\tmaskx=%d, masky=%d\n", tx, ty );
#endif /*DEBUG*/
switch( in->im->BandFmt ) {
@ -442,11 +438,12 @@ vips_interpolate_bicubic_class_init( VipsInterpolateBicubicClass *iclass )
for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
calculate_coefficients_catmull(
(float) x / VIPS_TRANSFORM_SCALE,
iclass->matrixf[x] );
vips_bicubic_matrixf[x] );
for( int i = 0; i < 4; i++ )
iclass->matrixi[x][i] =
iclass->matrixf[x][i] * VIPS_INTERPOLATE_SCALE;
vips_bicubic_matrixi[x][i] =
vips_bicubic_matrixf[x][i] *
VIPS_INTERPOLATE_SCALE;
}
}

View File

@ -305,22 +305,21 @@ vips_interpolate_nearest_static( void )
VIPS_TYPE_INTERPOLATE_BILINEAR, VipsInterpolateBilinearClass ))
typedef VipsInterpolate VipsInterpolateBilinear;
typedef struct _VipsInterpolateBilinearClass {
VipsInterpolateClass parent_class;
/* Precalculated interpolation matricies. int (used for pel sizes up
* to short), and float (for all others). We go to scale + 1, so
* we can round-to-nearest safely. Don't bother with double, since
* this is an approximation anyway.
*/
int matrixi[VIPS_TRANSFORM_SCALE + 1][VIPS_TRANSFORM_SCALE + 1][4];
float matrixd[VIPS_TRANSFORM_SCALE + 1][VIPS_TRANSFORM_SCALE + 1][4];
} VipsInterpolateBilinearClass;
typedef VipsInterpolateClass VipsInterpolateBilinearClass;
G_DEFINE_TYPE( VipsInterpolateBilinear, vips_interpolate_bilinear,
VIPS_TYPE_INTERPOLATE );
/* Precalculated interpolation matricies. int (used for pel sizes up
* to short), and float (for all others). We go to scale + 1 so
* we can round-to-nearest safely. Don't bother with double, since
* this is an approximation anyway.
*/
static int vips_bilinear_matrixi
[VIPS_TRANSFORM_SCALE + 1][VIPS_TRANSFORM_SCALE + 1][4];
static float vips_bilinear_matrixd
[VIPS_TRANSFORM_SCALE + 1][VIPS_TRANSFORM_SCALE + 1][4];
/* in this class, name vars in the 2x2 grid as eg.
* p1 p2
* p3 p4
@ -331,10 +330,10 @@ G_DEFINE_TYPE( VipsInterpolateBilinear, vips_interpolate_bilinear,
#define BILINEAR_INT( TYPE ) { \
TYPE *tq = (TYPE *) out; \
\
const int c1 = class->matrixi[xi][yi][0]; \
const int c2 = class->matrixi[xi][yi][1]; \
const int c3 = class->matrixi[xi][yi][2]; \
const int c4 = class->matrixi[xi][yi][3]; \
const int c1 = vips_bilinear_matrixi[xi][yi][0]; \
const int c2 = vips_bilinear_matrixi[xi][yi][1]; \
const int c3 = vips_bilinear_matrixi[xi][yi][2]; \
const int c4 = vips_bilinear_matrixi[xi][yi][3]; \
\
const TYPE *tp1 = (TYPE *) p1; \
const TYPE *tp2 = (TYPE *) p2; \
@ -351,10 +350,10 @@ G_DEFINE_TYPE( VipsInterpolateBilinear, vips_interpolate_bilinear,
#define BILINEAR_FLOAT( TYPE ) { \
TYPE *tq = (TYPE *) out; \
\
const double c1 = class->matrixd[xi][yi][0]; \
const double c2 = class->matrixd[xi][yi][1]; \
const double c3 = class->matrixd[xi][yi][2]; \
const double c4 = class->matrixd[xi][yi][3]; \
const double c1 = vips_bilinear_matrixd[xi][yi][0]; \
const double c2 = vips_bilinear_matrixd[xi][yi][1]; \
const double c3 = vips_bilinear_matrixd[xi][yi][2]; \
const double c4 = vips_bilinear_matrixd[xi][yi][3]; \
\
const TYPE *tp1 = (TYPE *) p1; \
const TYPE *tp2 = (TYPE *) p2; \
@ -388,9 +387,6 @@ static void
vips_interpolate_bilinear_interpolate( VipsInterpolate *interpolate,
PEL *out, REGION *in, double x, double y )
{
VipsInterpolateBilinearClass *class =
VIPS_INTERPOLATE_BILINEAR_GET_CLASS( interpolate );
/* Pel size and line size.
*/
const int ps = IM_IMAGE_SIZEOF_PEL( in->im );
@ -401,8 +397,11 @@ vips_interpolate_bilinear_interpolate( VipsInterpolate *interpolate,
*/
const double sx = x * VIPS_TRANSFORM_SCALE;
const double sy = y * VIPS_TRANSFORM_SCALE;
const int sxi = FAST_PSEUDO_FLOOR( sx );
const int syi = FAST_PSEUDO_FLOOR( sy );
/* We know sx/sy are always positive so we can just (int) them.
*/
const int sxi = (int) sx;
const int syi = (int) sy;
/* Get index into interpolation table and unscaled integer
* position.
@ -459,15 +458,19 @@ vips_interpolate_bilinear_class_init( VipsInterpolateBilinearClass *class )
c3 = Xd * Y;
c4 = X * Y;
class->matrixd[x][y][0] = c1;
class->matrixd[x][y][1] = c2;
class->matrixd[x][y][2] = c3;
class->matrixd[x][y][3] = c4;
vips_bilinear_matrixd[x][y][0] = c1;
vips_bilinear_matrixd[x][y][1] = c2;
vips_bilinear_matrixd[x][y][2] = c3;
vips_bilinear_matrixd[x][y][3] = c4;
class->matrixi[x][y][0] = c1 * VIPS_INTERPOLATE_SCALE;
class->matrixi[x][y][1] = c2 * VIPS_INTERPOLATE_SCALE;
class->matrixi[x][y][2] = c3 * VIPS_INTERPOLATE_SCALE;
class->matrixi[x][y][3] = c4 * VIPS_INTERPOLATE_SCALE;
vips_bilinear_matrixi[x][y][0] =
c1 * VIPS_INTERPOLATE_SCALE;
vips_bilinear_matrixi[x][y][1] =
c2 * VIPS_INTERPOLATE_SCALE;
vips_bilinear_matrixi[x][y][2] =
c3 * VIPS_INTERPOLATE_SCALE;
vips_bilinear_matrixi[x][y][3] =
c4 * VIPS_INTERPOLATE_SCALE;
}
}

View File

@ -60,22 +60,6 @@
#define FAST_MINMOD(a,b,ab,abminusaa) \
( (ab)>=0. ? ( (abminusaa)>=0. ? (a) : (b) ) : 0. )
/*
* Comment from Nicolas: I don't understand why the following restrict
* defs cannot be offloaded to config files.
*/
#ifndef restrict
#ifdef __restrict
#define restrict __restrict
#else
#ifdef __restrict__
#define restrict __restrict__
#else
#define restrict
#endif
#endif
#endif
/*
* Various casts which assume that the data is already in range. (That
* is, they are to be used with monotone samplers.)