Skip to content

Resolve rounding issues in reduce #1872

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Dec 12, 2020
Merged

Resolve rounding issues in reduce #1872

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
49d8051
Revert "Remove round-to-nearest behaviour"
kleisauke Oct 16, 2020
fdc140e
Round sum values to the nearest integer in *_notab
kleisauke Oct 18, 2020
4533375
Fix centre convention
kleisauke Oct 19, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 16 additions & 9 deletions libvips/resample/bicubic.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,15 +79,16 @@ typedef VipsInterpolate VipsInterpolateBicubic;
typedef VipsInterpolateClass VipsInterpolateBicubicClass;

/* 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
* works out slower since for many resizes the thing will no longer
* fit in L1.
*/
static int vips_bicubic_matrixi[VIPS_TRANSFORM_SCALE][4];
static double vips_bicubic_matrixf[VIPS_TRANSFORM_SCALE][4];
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.
*/
Expand Down Expand Up @@ -497,13 +498,19 @@ static void
vips_interpolate_bicubic_interpolate( VipsInterpolate *interpolate,
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 sy = y * VIPS_TRANSFORM_SCALE;
const int sx = x * VIPS_TRANSFORM_SCALE * 2;
const int sy = y * VIPS_TRANSFORM_SCALE * 2;

const int tx = sx & (VIPS_TRANSFORM_SCALE - 1);
const int ty = sy & (VIPS_TRANSFORM_SCALE - 1);
const int six = sx & (VIPS_TRANSFORM_SCALE * 2 - 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.
*/
Expand Down Expand Up @@ -636,7 +643,7 @@ vips_interpolate_bicubic_class_init( VipsInterpolateBicubicClass *iclass )

/* 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],
(float) x / VIPS_TRANSFORM_SCALE );

Expand Down
27 changes: 16 additions & 11 deletions libvips/resample/reduceh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,6 @@
* 6/6/20 kleisauke
* - deprecate @centre option, it's now always on
* - fix pixel shift
* - remove unnecessary round-to-nearest behaviour
*/

/*
Expand Down Expand Up @@ -79,10 +78,11 @@ typedef struct _VipsReduceh {
double hoffset;

/* 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];
double *matrixf[VIPS_TRANSFORM_SCALE];
int *matrixi[VIPS_TRANSFORM_SCALE + 1];
double *matrixf[VIPS_TRANSFORM_SCALE + 1];

/* Deprecated.
*/
Expand Down Expand Up @@ -281,8 +281,12 @@ reduceh_notab( VipsReduceh *reduceh,

vips_reduce_make_mask( cx, reduceh->kernel, reduceh->hshrink, x );

for( int z = 0; z < bands; z++ )
out[z] = reduce_sum<T, double>( in + z, bands, cx, n );
for( int z = 0; z < bands; z++ ) {
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
Expand Down Expand Up @@ -346,8 +350,9 @@ vips_reduceh_gen( VipsRegion *out_region, void *seq,
for( int x = 0; x < r->width; x++ ) {
const int ix = (int) X;
VipsPel *p = p0 + ix * ps;
const int sx = X * VIPS_TRANSFORM_SCALE;
const int tx = sx & (VIPS_TRANSFORM_SCALE - 1);
const int sx = X * VIPS_TRANSFORM_SCALE * 2;
const int six = sx & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int tx = (six + 1) >> 1;
const int *cxi = reduceh->matrixi[tx];
const double *cxf = reduceh->matrixf[tx];

Expand Down Expand Up @@ -475,11 +480,11 @@ vips_reduceh_build( VipsObject *object )
* by half that distance so that we discard pixels equally
* 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.
*/
for( int x = 0; x < VIPS_TRANSFORM_SCALE; x++ ) {
for( int x = 0; x < VIPS_TRANSFORM_SCALE + 1; x++ ) {
reduceh->matrixf[x] =
VIPS_ARRAY( object, reduceh->n_point, double );
reduceh->matrixi[x] =
Expand Down Expand Up @@ -513,7 +518,7 @@ vips_reduceh_build( VipsObject *object )
/* Add new pixels around the input so we can interpolate at the edges.
*/
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,
"extend", VIPS_EXTEND_COPY,
(void *) NULL ) )
Expand Down
38 changes: 22 additions & 16 deletions libvips/resample/reducev.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,6 @@
* - deprecate @centre option, it's now always on
* - fix pixel shift
* - speed up the mask construction for uchar/ushort images
* - remove unnecessary round-to-nearest behaviour
*/

/*
Expand Down Expand Up @@ -118,14 +117,15 @@ typedef struct _VipsReducev {
double voffset;

/* 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];
double *matrixf[VIPS_TRANSFORM_SCALE];
int *matrixi[VIPS_TRANSFORM_SCALE + 1];
double *matrixf[VIPS_TRANSFORM_SCALE + 1];

/* 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.
*/
Expand Down Expand Up @@ -154,7 +154,7 @@ vips_reducev_finalize( GObject *gobject )
for( int i = 0; i < reducev->n_pass; i++ )
VIPS_FREEF( vips_vector_free, reducev->pass[i].vector );
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->matrixi[i] );
VIPS_FREE( reducev->matrixo[i] );
Expand Down Expand Up @@ -511,8 +511,12 @@ reducev_notab( VipsReducev *reducev,

vips_reduce_make_mask( cy, reducev->kernel, reducev->vshrink, y );

for( int z = 0; z < ne; z++ )
out[z] = reduce_sum<T, double>( in + z, l1, cy, n );
for( int z = 0; z < ne; z++ ) {
double sum;
sum = reduce_sum<T, double>( in + z, l1, cy, n );

out[z] = VIPS_ROUND_UINT( sum );
}
}

static int
Expand Down Expand Up @@ -555,8 +559,9 @@ vips_reducev_gen( VipsRegion *out_region, void *vseq,
VIPS_REGION_ADDR( out_region, r->left, r->top + y );
const int py = (int) Y;
VipsPel *p = VIPS_REGION_ADDR( ir, r->left, py );
const int sy = Y * VIPS_TRANSFORM_SCALE;
const int ty = sy & (VIPS_TRANSFORM_SCALE - 1);
const int sy = Y * VIPS_TRANSFORM_SCALE * 2;
const int siy = sy & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int ty = (siy + 1) >> 1;
const int *cyi = reducev->matrixi[ty];
const double *cyf = reducev->matrixf[ty];
const int lskip = VIPS_REGION_LSKIP( ir );
Expand Down Expand Up @@ -677,8 +682,9 @@ vips_reducev_vector_gen( VipsRegion *out_region, void *vseq,
VipsPel *q =
VIPS_REGION_ADDR( out_region, r->left, r->top + y );
const int py = (int) Y;
const int sy = Y * VIPS_TRANSFORM_SCALE;
const int ty = sy & (VIPS_TRANSFORM_SCALE - 1);
const int sy = Y * VIPS_TRANSFORM_SCALE * 2;
const int siy = sy & (VIPS_TRANSFORM_SCALE * 2 - 1);
const int ty = (siy + 1) >> 1;
const int *cyo = reducev->matrixo[ty];

#ifdef DEBUG_PIXELS
Expand Down Expand Up @@ -739,7 +745,7 @@ vips_reducev_raw( VipsReducev *reducev, VipsImage *in, VipsImage **out )
*/
if( in->BandFmt == VIPS_FORMAT_UCHAR &&
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] =
VIPS_ARRAY( NULL, reducev->n_point, int );
if( !reducev->matrixo[y] )
Expand Down Expand Up @@ -847,11 +853,11 @@ vips_reducev_build( VipsObject *object )
* by half that distance so that we discard pixels equally
* 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.
*/
for( int y = 0; y < VIPS_TRANSFORM_SCALE; y++ ) {
for( int y = 0; y < VIPS_TRANSFORM_SCALE + 1; y++ ) {
reducev->matrixf[y] =
VIPS_ARRAY( NULL, reducev->n_point, double );
reducev->matrixi[y] =
Expand Down Expand Up @@ -885,7 +891,7 @@ vips_reducev_build( VipsObject *object )
/* Add new pixels around the input so we can interpolate at the edges.
*/
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,
"extend", VIPS_EXTEND_COPY,
(void *) NULL ) )
Expand Down