sergiotarxz/libvips
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge pull request #1872 from kleisauke/fix-rounding-reduce

Browse Source 
Resolve rounding issues in reduce
This commit is contained in:
John Cupitt 2020-12-12 14:51:15 +00:00 committed by GitHub
commit a54cec93ac
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
3 changed files with 54 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
libvips/resample/bicubic.cpp
View File

@ -79,15 +79,16 @@ typedef VipsInterpolate VipsInterpolateBicubic;
typedef VipsInterpolateClass VipsInterpolateBicubicClass; typedef VipsInterpolateClass VipsInterpolateBicubicClass;
/* Precalculated interpolation matrices. int (used for pel /* Precalculated interpolation matrices. int (used for pel
* sizes up to short), and double (for all others). * 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 /* 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 * works out slower since for many resizes the thing will no longer
* fit in L1. * fit in L1.
*/ */
static int vips_bicubic_matrixi[VIPS_TRANSFORM_SCALE][4]; static int vips_bicubic_matrixi[VIPS_TRANSFORM_SCALE + 1][4];
static double vips_bicubic_matrixf[VIPS_TRANSFORM_SCALE][4]; static double vips_bicubic_matrixf[VIPS_TRANSFORM_SCALE + 1][4];
/* We need C linkage for this. /* We need C linkage for this.
*/ */
@ -497,13 +498,19 @@ static void
vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate, vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
void *out, VipsRegion *in, double x, double y ) void *out, VipsRegion *in, double x, double y )
{ {
/* Find the mask index. /* Find the mask index. We round-to-nearest, so we need to generate
* indexes in 0 to VIPS_TRANSFORM_SCALE, 2^n + 1 values. We multiply
* by 2 more than we need to, add one, mask, then shift down again to
* get the extra range.
*/ */
const int sx = x * VIPS_TRANSFORM_SCALE; const int sx = x * VIPS_TRANSFORM_SCALE * 2;
const int sy = y * VIPS_TRANSFORM_SCALE; const int sy = y * VIPS_TRANSFORM_SCALE * 2;
const int tx = sx & (VIPS_TRANSFORM_SCALE - 1); const int six = sx & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int ty = sy & (VIPS_TRANSFORM_SCALE - 1); const int siy = sy & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int tx = (six + 1) >> 1;
const int ty = (siy + 1) >> 1;
/* We know x/y are always positive, so we can just (int) them. /* We know x/y are always positive, so we can just (int) them.
*/ */
@ -636,7 +643,7 @@ vips_interpolate_bicubic_class_init( VipsInterpolateBicubicClass *iclass )
/* Build the tables of pre-computed coefficients. /* Build the tables of pre-computed coefficients.
*/ */
for( int x = 0; x < VIPS_TRANSFORM_SCALE; x++ ) { for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
calculate_coefficients_catmull( vips_bicubic_matrixf[x], calculate_coefficients_catmull( vips_bicubic_matrixf[x],
(float) x / VIPS_TRANSFORM_SCALE ); (float) x / VIPS_TRANSFORM_SCALE );

27
libvips/resample/reduceh.cpp
View File

@ -11,7 +11,6 @@
* 6/6/20 kleisauke * 6/6/20 kleisauke
* - deprecate @centre option, it's now always on * - deprecate @centre option, it's now always on
* - fix pixel shift * - fix pixel shift
* - remove unnecessary round-to-nearest behaviour
*/ */
/* /*
@ -79,10 +78,11 @@ typedef struct _VipsReduceh {
double hoffset; double hoffset;
/* Precalculated interpolation matrices. int (used for pel /* Precalculated interpolation matrices. int (used for pel
* sizes up to short), and double (for all others). * sizes up to short), and double (for all others). We go to
* scale + 1 so we can round-to-nearest safely.
*/ */
int *matrixi[VIPS_TRANSFORM_SCALE]; int *matrixi[VIPS_TRANSFORM_SCALE + 1];
double *matrixf[VIPS_TRANSFORM_SCALE]; double *matrixf[VIPS_TRANSFORM_SCALE + 1];
/* Deprecated. /* Deprecated.
*/ */
@ -281,8 +281,12 @@ reduceh_notab( VipsReduceh *reduceh,
vips_reduce_make_mask( cx, reduceh->kernel, reduceh->hshrink, x ); vips_reduce_make_mask( cx, reduceh->kernel, reduceh->hshrink, x );
for( int z = 0; z < bands; z++ ) for( int z = 0; z < bands; z++ ) {
out[z] = reduce_sum<T, double>( in + z, bands, cx, n ); double sum;
sum = reduce_sum<T, double>( in + z, bands, cx, n );
out[z] = VIPS_ROUND_UINT( sum );
}
} }
/* Tried a vector path (see reducev) but it was slower. The vectors for /* Tried a vector path (see reducev) but it was slower. The vectors for
@ -346,8 +350,9 @@ vips_reduceh_gen( VipsRegion *out_region, void *seq,
for( int x = 0; x < r->width; x++ ) { for( int x = 0; x < r->width; x++ ) {
const int ix = (int) X; const int ix = (int) X;
VipsPel *p = p0 + ix * ps; VipsPel *p = p0 + ix * ps;
const int sx = X * VIPS_TRANSFORM_SCALE; const int sx = X * VIPS_TRANSFORM_SCALE * 2;
const int tx = sx & (VIPS_TRANSFORM_SCALE - 1); const int six = sx & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int tx = (six + 1) >> 1;
const int *cxi = reduceh->matrixi[tx]; const int *cxi = reduceh->matrixi[tx];
const double *cxf = reduceh->matrixf[tx]; const double *cxf = reduceh->matrixf[tx];
@ -475,11 +480,11 @@ vips_reduceh_build( VipsObject *object )
* by half that distance so that we discard pixels equally * by half that distance so that we discard pixels equally
* from left and right. * from left and right.
*/ */
reduceh->hoffset = (1 + extra_pixels) / 2.0; reduceh->hoffset = (1 + extra_pixels) / 2.0 - 1;
/* Build the tables of pre-computed coefficients. /* Build the tables of pre-computed coefficients.
*/ */
for( int x = 0; x < VIPS_TRANSFORM_SCALE; x++ ) { for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
reduceh->matrixf[x] = reduceh->matrixf[x] =
VIPS_ARRAY( object, reduceh->n_point, double ); VIPS_ARRAY( object, reduceh->n_point, double );
reduceh->matrixi[x] = reduceh->matrixi[x] =
@ -513,7 +518,7 @@ vips_reduceh_build( VipsObject *object )
/* Add new pixels around the input so we can interpolate at the edges. /* Add new pixels around the input so we can interpolate at the edges.
*/ */
if( vips_embed( in, &t[1], if( vips_embed( in, &t[1],
reduceh->n_point / 2 - 1, 0, VIPS_CEIL( reduceh->n_point / 2.0 ) - 1, 0,
in->Xsize + reduceh->n_point, in->Ysize, in->Xsize + reduceh->n_point, in->Ysize,
"extend", VIPS_EXTEND_COPY, "extend", VIPS_EXTEND_COPY,
(void *) NULL ) ) (void *) NULL ) )

38
libvips/resample/reducev.cpp
View File

@ -21,7 +21,6 @@
* - deprecate @centre option, it's now always on * - deprecate @centre option, it's now always on
* - fix pixel shift * - fix pixel shift
* - speed up the mask construction for uchar/ushort images * - speed up the mask construction for uchar/ushort images
* - remove unnecessary round-to-nearest behaviour
*/ */
/* /*
@ -118,14 +117,15 @@ typedef struct _VipsReducev {
double voffset; double voffset;
/* Precalculated interpolation matrices. int (used for pel /* Precalculated interpolation matrices. int (used for pel
* sizes up to short), and double (for all others). * sizes up to short), and double (for all others). We go to
* scale + 1 so we can round-to-nearest safely.
*/ */
int *matrixi[VIPS_TRANSFORM_SCALE]; int *matrixi[VIPS_TRANSFORM_SCALE + 1];
double *matrixf[VIPS_TRANSFORM_SCALE]; double *matrixf[VIPS_TRANSFORM_SCALE + 1];
/* And another set for orc: we want 2.6 precision. /* And another set for orc: we want 2.6 precision.
*/ */
int *matrixo[VIPS_TRANSFORM_SCALE]; int *matrixo[VIPS_TRANSFORM_SCALE + 1];
/* The passes we generate for this mask. /* The passes we generate for this mask.
*/ */
@ -154,7 +154,7 @@ vips_reducev_finalize( GObject *gobject )
for( int i = 0; i < reducev->n_pass; i++ ) for( int i = 0; i < reducev->n_pass; i++ )
VIPS_FREEF( vips_vector_free, reducev->pass[i].vector ); VIPS_FREEF( vips_vector_free, reducev->pass[i].vector );
reducev->n_pass = 0; reducev->n_pass = 0;
for( int i = 0; i < VIPS_TRANSFORM_SCALE; i++ ) { for( int i = 0; i < VIPS_TRANSFORM_SCALE + 1; i++ ) {
VIPS_FREE( reducev->matrixf[i] ); VIPS_FREE( reducev->matrixf[i] );
VIPS_FREE( reducev->matrixi[i] ); VIPS_FREE( reducev->matrixi[i] );
VIPS_FREE( reducev->matrixo[i] ); VIPS_FREE( reducev->matrixo[i] );
@ -511,8 +511,12 @@ reducev_notab( VipsReducev *reducev,
vips_reduce_make_mask( cy, reducev->kernel, reducev->vshrink, y ); vips_reduce_make_mask( cy, reducev->kernel, reducev->vshrink, y );
for( int z = 0; z < ne; z++ ) for( int z = 0; z < ne; z++ ) {
out[z] = reduce_sum<T, double>( in + z, l1, cy, n ); double sum;
sum = reduce_sum<T, double>( in + z, l1, cy, n );
out[z] = VIPS_ROUND_UINT( sum );
}
} }
static int static int
@ -555,8 +559,9 @@ vips_reducev_gen( VipsRegion *out_region, void *vseq,
VIPS_REGION_ADDR( out_region, r->left, r->top + y ); VIPS_REGION_ADDR( out_region, r->left, r->top + y );
const int py = (int) Y; const int py = (int) Y;
VipsPel *p = VIPS_REGION_ADDR( ir, r->left, py ); VipsPel *p = VIPS_REGION_ADDR( ir, r->left, py );
const int sy = Y * VIPS_TRANSFORM_SCALE; const int sy = Y * VIPS_TRANSFORM_SCALE * 2;
const int ty = sy & (VIPS_TRANSFORM_SCALE - 1); const int siy = sy & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int ty = (siy + 1) >> 1;
const int *cyi = reducev->matrixi[ty]; const int *cyi = reducev->matrixi[ty];
const double *cyf = reducev->matrixf[ty]; const double *cyf = reducev->matrixf[ty];
const int lskip = VIPS_REGION_LSKIP( ir ); const int lskip = VIPS_REGION_LSKIP( ir );
@ -677,8 +682,9 @@ vips_reducev_vector_gen( VipsRegion *out_region, void *vseq,
VipsPel *q = VipsPel *q =
VIPS_REGION_ADDR( out_region, r->left, r->top + y ); VIPS_REGION_ADDR( out_region, r->left, r->top + y );
const int py = (int) Y; const int py = (int) Y;
const int sy = Y * VIPS_TRANSFORM_SCALE; const int sy = Y * VIPS_TRANSFORM_SCALE * 2;
const int ty = sy & (VIPS_TRANSFORM_SCALE - 1); const int siy = sy & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int ty = (siy + 1) >> 1;
const int *cyo = reducev->matrixo[ty]; const int *cyo = reducev->matrixo[ty];
#ifdef DEBUG_PIXELS #ifdef DEBUG_PIXELS
@ -739,7 +745,7 @@ vips_reducev_raw( VipsReducev *reducev, VipsImage *in, VipsImage **out )
*/ */
if( in->BandFmt == VIPS_FORMAT_UCHAR && if( in->BandFmt == VIPS_FORMAT_UCHAR &&
vips_vector_isenabled() ) vips_vector_isenabled() )
for( int y = 0; y < VIPS_TRANSFORM_SCALE; y++ ) { for( int y = 0; y < VIPS_TRANSFORM_SCALE + 1; y++ ) {
reducev->matrixo[y] = reducev->matrixo[y] =
VIPS_ARRAY( NULL, reducev->n_point, int ); VIPS_ARRAY( NULL, reducev->n_point, int );
if( !reducev->matrixo[y] ) if( !reducev->matrixo[y] )
@ -847,11 +853,11 @@ vips_reducev_build( VipsObject *object )
* by half that distance so that we discard pixels equally * by half that distance so that we discard pixels equally
* from left and right. * from left and right.
*/ */
reducev->voffset = (1 + extra_pixels) / 2.0; reducev->voffset = (1 + extra_pixels) / 2.0 - 1;
/* Build the tables of pre-computed coefficients. /* Build the tables of pre-computed coefficients.
*/ */
for( int y = 0; y < VIPS_TRANSFORM_SCALE; y++ ) { for( int y = 0; y < VIPS_TRANSFORM_SCALE + 1; y++ ) {
reducev->matrixf[y] = reducev->matrixf[y] =
VIPS_ARRAY( NULL, reducev->n_point, double ); VIPS_ARRAY( NULL, reducev->n_point, double );
reducev->matrixi[y] = reducev->matrixi[y] =
@ -885,7 +891,7 @@ vips_reducev_build( VipsObject *object )
/* Add new pixels around the input so we can interpolate at the edges. /* Add new pixels around the input so we can interpolate at the edges.
*/ */
if( vips_embed( in, &t[1], if( vips_embed( in, &t[1],
0, reducev->n_point / 2 - 1, 0, VIPS_CEIL( reducev->n_point / 2.0 ) - 1,
in->Xsize, in->Ysize + reducev->n_point, in->Xsize, in->Ysize + reducev->n_point,
"extend", VIPS_EXTEND_COPY, "extend", VIPS_EXTEND_COPY,
(void *) NULL ) ) (void *) NULL ) )