diff --git a/libvips/resample/presample.h b/libvips/resample/presample.h index 18cb4add00..b2799beaf7 100644 --- a/libvips/resample/presample.h +++ b/libvips/resample/presample.h @@ -70,8 +70,6 @@ GType vips_resample_get_type(void); #define MAX_POINT (2000) int vips_reduce_get_points(VipsKernel kernel, double shrink); -void vips_reduce_make_mask(double *c, - VipsKernel kernel, double shrink, double x); void vips_reduceh_uchar_hwy(VipsPel *pout, VipsPel *pin, int n, int width, int bands, diff --git a/libvips/resample/reduceh.cpp b/libvips/resample/reduceh.cpp index 8322f99c53..afbe91b3f0 100644 --- a/libvips/resample/reduceh.cpp +++ b/libvips/resample/reduceh.cpp @@ -120,8 +120,6 @@ vips_reduce_get_points(VipsKernel kernel, double shrink) return 2 * rint(2 * shrink) + 1; case VIPS_KERNEL_LANCZOS2: - /* Needs to be in sync with calculate_coefficients_lanczos(). - */ return 2 * rint(2 * shrink) + 1; case VIPS_KERNEL_LANCZOS3: @@ -133,45 +131,6 @@ vips_reduce_get_points(VipsKernel kernel, double shrink) } } -/* Calculate a mask element. - */ -void -vips_reduce_make_mask(double *c, VipsKernel kernel, double shrink, double x) -{ - switch (kernel) { - case VIPS_KERNEL_NEAREST: - c[0] = 1.0; - break; - - case VIPS_KERNEL_LINEAR: - calculate_coefficients_triangle(c, shrink, x); - break; - - case VIPS_KERNEL_CUBIC: - /* Catmull-Rom. - */ - calculate_coefficients_cubic(c, shrink, x, 0.0, 0.5); - break; - - case VIPS_KERNEL_MITCHELL: - calculate_coefficients_cubic(c, shrink, x, - 1.0 / 3.0, 1.0 / 3.0); - break; - - case VIPS_KERNEL_LANCZOS2: - calculate_coefficients_lanczos(c, 2, shrink, x); - break; - - case VIPS_KERNEL_LANCZOS3: - calculate_coefficients_lanczos(c, 3, shrink, x); - break; - - default: - g_assert_not_reached(); - break; - } -} - template static void inline reduceh_unsigned_int_tab(VipsReduceh *reduceh, VipsPel *pout, const VipsPel *pin, @@ -275,7 +234,8 @@ static void inline reduceh_notab(VipsReduceh *reduceh, double cx[MAX_POINT]; - vips_reduce_make_mask(cx, reduceh->kernel, reduceh->hshrink, x); + vips_reduce_make_mask(cx, reduceh->kernel, reduceh->n_point, + reduceh->hshrink, x); for (int z = 0; z < bands; z++) { double sum; @@ -555,8 +515,8 @@ vips_reduceh_build(VipsObject *object) !reduceh->matrixs[x]) return -1; - vips_reduce_make_mask(reduceh->matrixf[x], - reduceh->kernel, reduceh->hshrink, + vips_reduce_make_mask(reduceh->matrixf[x], reduceh->kernel, + reduceh->n_point, reduceh->hshrink, (float) x / VIPS_TRANSFORM_SCALE); for (int i = 0; i < reduceh->n_point; i++) diff --git a/libvips/resample/reducev.cpp b/libvips/resample/reducev.cpp index 190a6efd9c..88164547fd 100644 --- a/libvips/resample/reducev.cpp +++ b/libvips/resample/reducev.cpp @@ -511,7 +511,8 @@ static void inline reducev_notab(VipsReducev *reducev, double cy[MAX_POINT]; - vips_reduce_make_mask(cy, reducev->kernel, reducev->vshrink, y); + vips_reduce_make_mask(cy, reducev->kernel, reducev->n_point, + reducev->vshrink, y); for (int z = 0; z < ne; z++) { double sum; @@ -956,8 +957,8 @@ vips_reducev_build(VipsObject *object) !reducev->matrixs[y]) return -1; - vips_reduce_make_mask(reducev->matrixf[y], - reducev->kernel, reducev->vshrink, + vips_reduce_make_mask(reducev->matrixf[y], reducev->kernel, + reducev->n_point, reducev->vshrink, (float) y / VIPS_TRANSFORM_SCALE); for (int i = 0; i < reducev->n_point; i++) diff --git a/libvips/resample/templates.h b/libvips/resample/templates.h index 66fdcb17a6..c0b54c5825 100644 --- a/libvips/resample/templates.h +++ b/libvips/resample/templates.h @@ -28,8 +28,6 @@ */ -#include - /* * Various casts which assume that the data is already in range. (That * is, they are to be used with monotone samplers.) @@ -312,124 +310,117 @@ static void inline calculate_coefficients_catmull(double c[4], const double x) 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. +/* 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_triangle(double *c, - const double shrink, const double x) +static double inline cubic_filter(double x, double B, double C) { - /* 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; + const double axp = VIPS_FABS(x); + const double axp2 = axp * axp; + const double axp3 = axp2 * axp; + + if (axp <= 1) + return ((12 - 9 * B - 6 * C) * axp3 + + (-18 + 12 * B + 6 * C) * axp2 + + (6 - 2 * B)) / + 6; + + if (axp <= 2) + return ((-B - 6 * C) * axp3 + + (6 * B + 30 * C) * axp2 + + (-12 * B - 48 * C) * axp + + (8 * B + 24 * C)) / + 6; + + return 0.0; +} - int i; - double sum; +static double inline sinc_filter(double x) +{ + if (x == 0.0) + return 1.0; - sum = 0; - for (i = 0; i < n_points; i++) { - const double xp = (i - half) / shrink; + x = x * VIPS_PI; - double l; + return sin(x) / x; +} - l = 1.0 - VIPS_FABS(xp); - l = VIPS_MAX(0.0, l); +using VipsFilterFn = double (*)(double); - c[i] = l; - sum += l; - } +template +static double inline filter(double x); - for (i = 0; i < n_points; i++) - c[i] /= sum; +template <> +double inline filter(double x) +{ + if (x < 0.0) + x = -x; + + if (x < 1.0) + return 1.0 - x; + + return 0.0; } -/* 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 +/* Catmull-Rom. */ -static void inline calculate_coefficients_cubic(double *c, - const double shrink, const double x, double B, double C) +template <> +double inline filter(double x) { - /* 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; + return cubic_filter(x, 0.0, 0.5); +} - int i; - double sum; +template <> +double inline filter(double x) +{ + return cubic_filter(x, 1.0 / 3.0, 1.0 / 3.0); +} - 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; +template <> +double inline filter(double x) +{ + if (x >= -2 && x <= 2) + return sinc_filter(x) * sinc_filter(x / 2); - c[i] = l; - sum += l; - } + return 0.0; +} - for (i = 0; i < n_points; i++) - c[i] /= sum; +template <> +double inline filter(double x) +{ + if (x >= -3 && x <= 3) + return sinc_filter(x) * sinc_filter(x / 3); + + return 0.0; } /* 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 + * calculate c0 .. c(@n_points), the 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. + * @shrink is the reduction factor, so 1 for interpolation, 2 for a + * x2 reduction, etc. */ -static void inline calculate_coefficients_lanczos(double *c, - const int a, const double shrink, const double x) +template +static void +calculate_coefficients(T *c, const int n_points, + VipsFilterFn filter_fn, 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; + T sum; - sum = 0; + sum = 0.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); + double l = filter_fn(xp); c[i] = l; sum += l; @@ -439,41 +430,61 @@ static void inline calculate_coefficients_lanczos(double *c, c[i] /= sum; } -/* Simplified version of std::enable_if::type +/* Calculate a mask element. */ -template -using Requires = typename std::enable_if::type; /* C++11 */ -// using Requires = std::enable_if_t; /* C++14 */ - -/* Our inner loop for resampling with a convolution. Operate on elements of - * type T, gather results in an intermediate of type IT. - */ -template ::value> = true> -static IT -reduce_sum(const T *restrict in, int stride, const short *restrict c, int n) +template +static void +vips_reduce_make_mask(T *c, VipsKernel kernel, const int n_points, + const double shrink, const double x) { - IT sum; - - sum = 0; - for (int i = 0; i < n; i++) { - sum += (IT) c[i] * in[0]; - in += stride; + switch (kernel) { + case VIPS_KERNEL_NEAREST: + c[0] = 1.0; + break; + + case VIPS_KERNEL_LINEAR: + calculate_coefficients(c, n_points, + filter, shrink, x); + break; + + case VIPS_KERNEL_CUBIC: + calculate_coefficients(c, n_points, + filter, shrink, x); + break; + + case VIPS_KERNEL_MITCHELL: + calculate_coefficients(c, n_points, + filter, shrink, x); + break; + + case VIPS_KERNEL_LANCZOS2: + calculate_coefficients(c, n_points, + filter, shrink, x); + break; + + case VIPS_KERNEL_LANCZOS3: + calculate_coefficients(c, n_points, + filter, shrink, x); + break; + + default: + g_assert_not_reached(); + break; } - - return sum; } -/* Same as above, but specialized for floating point types. +/* Our inner loop for resampling with a convolution of type CT. Operate on + * elements of type T, gather results in an intermediate of type IT. */ -template ::value> = true> -static IT -reduce_sum(const T *restrict in, int stride, const double *restrict c, int n) +template +static IT inline reduce_sum(const T *restrict in, int stride, + const CT *restrict c, int n) { IT sum; sum = 0; for (int i = 0; i < n; i++) { - sum += c[i] * in[0]; + sum += (IT) c[i] * in[0]; in += stride; }