| /* |
| * Copyright (c) 2021, Alliance for Open Media. All rights reserved |
| * |
| * This source code is subject to the terms of the BSD 3-Clause Clear License |
| * and the Alliance for Open Media Patent License 1.0. If the BSD 3-Clause Clear |
| * License was not distributed with this source code in the LICENSE file, you |
| * can obtain it at aomedia.org/license/software-license/bsd-3-c-c/. If the |
| * Alliance for Open Media Patent License 1.0 was not distributed with this |
| * source code in the PATENTS file, you can obtain it at |
| * aomedia.org/license/patent-license/. |
| */ |
| |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <stdint.h> |
| #include <stdbool.h> |
| #include <string.h> |
| #include <assert.h> |
| |
| #include "tools/lanczos/lanczos_resample.h" |
| |
| /* Shift down with rounding for use when n >= 0, value >= 0 */ |
| #define ROUND_POWER_OF_TWO(value, n) (((value) + (((1 << (n)) >> 1))) >> (n)) |
| |
| /* Shift down with rounding for signed integers, for use when n >= 0 */ |
| #define ROUND_POWER_OF_TWO_SIGNED(value, n) \ |
| (((value) < 0) ? -ROUND_POWER_OF_TWO(-(value), (n)) \ |
| : ROUND_POWER_OF_TWO((value), (n))) |
| |
| #define MAX(a, b) ((a) < (b) ? (b) : (a)) |
| #define MIN(a, b) ((a) < (b) ? (a) : (b)) |
| |
| #ifndef M_PI |
| #define M_PI (3.14159265358979323846) |
| #endif |
| |
| double get_centered_x0(int p, int q) { return (double)(q - p) / (2 * p); } |
| |
| double get_cosited_chroma_x0(int p, int q) { return (double)(q - p) / (4 * p); } |
| |
| double get_inverse_x0_numeric(int p, int q, double x0) { return -x0 * p / q; } |
| |
| double get_inverse_x0(int p, int q, double x0, int subsampled) { |
| if (x0 == (double)('c')) |
| x0 = get_centered_x0(p, q); |
| else if (x0 == (double)('d')) |
| x0 = subsampled ? get_cosited_chroma_x0(p, q) : get_centered_x0(p, q); |
| return get_inverse_x0_numeric(p, q, x0); |
| } |
| |
| static inline int doclip(int x, int low, int high) { |
| return (x < low ? low : x > high ? high : x); |
| } |
| |
| void show_resample_filter(RationalResampleFilter *rf) { |
| printf("Resample factor: %d / %d\n", rf->p, rf->q); |
| printf("Extension type: %s\n", ext2str(rf->ext_type)); |
| printf("Start = %d\n", rf->start); |
| printf("Steps = "); |
| for (int i = 0; i < rf->p; ++i) { |
| printf("%d, ", rf->steps[i]); |
| } |
| printf("\n"); |
| printf("Phases = "); |
| for (int i = 0; i < rf->p; ++i) { |
| printf("%f, ", rf->phases[i]); |
| } |
| printf("\n"); |
| printf("Filters [length %d, bits %d]:\n", rf->length, rf->filter_bits); |
| for (int i = 0; i < rf->p; ++i) { |
| printf(" { "); |
| for (int j = 0; j < rf->length; ++j) printf("%d, ", rf->filter[i][j]); |
| printf(" }\n"); |
| } |
| printf("\n"); |
| } |
| |
| static double sinc(double x) { |
| if (fabs(x) < 1e-12) return 1.0; |
| return sin(M_PI * x) / (M_PI * x); |
| } |
| |
| static double mod_bessel_first(double x) { |
| const double t = 0.25 * x * x; |
| double fact = 1.0; |
| double tpow = 1.0; |
| double v = 1.0; |
| double dv; |
| int k = 1; |
| do { |
| fact *= k; |
| tpow *= t; |
| dv = tpow / (fact * fact); |
| v += dv; |
| k++; |
| } while (fabs(dv) > fabs(v) * 1e-8); |
| return v; |
| } |
| |
| // This is a window function assumed to be defined between [-1, 1] and |
| // with the value at y=0 being 1. |
| static double window(double y, WIN_TYPE win) { |
| switch (win) { |
| case WIN_LANCZOS: { |
| return sinc(y); |
| } |
| case WIN_LANCZOS_DIL: { |
| return sinc(y * 0.95); |
| } |
| case WIN_GAUSSIAN: { |
| const double sigma = 0.66; |
| const double sigma2 = sigma * sigma; |
| return exp(-y * y / sigma2); |
| } |
| case WIN_GENGAUSSIAN: { |
| const double alpha = 4; |
| const double sigma = 0.78; |
| return exp(-pow(fabs(y / sigma), alpha)); |
| } |
| case WIN_COSINE: { |
| return cos(M_PI * y / 2); |
| } |
| case WIN_HAMMING: { |
| const double a0 = 25.0 / 46.0; |
| const double a1 = 1.0 - a0; |
| return (a0 + a1 * cos(M_PI * y)); |
| } |
| case WIN_BLACKMAN: { |
| const double a0 = 0.42659; |
| const double a1 = 0.49656; |
| const double a2 = 1.0 - a0 - a1; |
| return a0 + a1 * cos(M_PI * y) + a2 * cos(2 * M_PI * y); |
| } |
| case WIN_KAISER: { |
| const double alpha = 1.32; |
| const double u = M_PI * alpha; |
| const double v = M_PI * alpha * sqrt(1 - y * y); |
| return mod_bessel_first(v) / mod_bessel_first(u); |
| } |
| default: { |
| assert(0 && "Unknown window type"); |
| return 0; |
| } |
| } |
| } |
| |
| static double kernel(double x, int a, WIN_TYPE win_type) { |
| const double absx = fabs(x); |
| if (absx < (double)a) { |
| return sinc(x) * window(x / a, win_type); |
| } else { |
| return 0.0; |
| } |
| } |
| |
| static int get_lanczos_downsampler_filter_length(int p, int q, int a) { |
| assert(p < q); |
| return 2 * ((a * q + p - 1) / p); |
| } |
| |
| static int get_lanczos_upsampler_filter_length(int p, int q, int a) { |
| (void)p; |
| (void)q; |
| assert(p >= q); |
| return 2 * a; |
| } |
| |
| static void integerize_array(double *x, int len, int bits, int16_t *y) { |
| int sumy = 0; |
| for (int i = 0; i < len; ++i) { |
| y[i] = (int16_t)rint(x[i] * (1 << bits)); |
| sumy += y[i]; |
| } |
| while (sumy > (1 << bits)) { |
| double mx = -65536.0; |
| int imx = -1; |
| for (int i = 0; i < len; ++i) { |
| const double v = (double)y[i] - (x[i] * (1 << bits)); |
| if (v > mx) { |
| mx = v; |
| imx = i; |
| } |
| } |
| y[imx] -= 1; |
| sumy -= 1; |
| } |
| while (sumy < (1 << bits)) { |
| double mx = 65536.0; |
| int imx = -1; |
| for (int i = 0; i < len; ++i) { |
| const double v = (double)y[i] - (x[i] * (1 << bits)); |
| if (v < mx) { |
| mx = v; |
| imx = i; |
| } |
| } |
| y[imx] += 1; |
| sumy += 1; |
| } |
| sumy = 0; |
| for (int i = 0; i < len; ++i) { |
| sumy += y[i]; |
| } |
| assert(sumy == (1 << bits)); |
| } |
| |
| static void get_lanczos_downsampler(double x, int p, int q, int a, int bits, |
| WIN_TYPE win_type, int16_t *ifilter) { |
| double filter[MAX_FILTER_LEN] = { 0.0 }; |
| int tapsby2 = get_lanczos_downsampler_filter_length(p, q, a) / 2; |
| assert(tapsby2 * 2 <= MAX_FILTER_LEN); |
| double filter_sum = 0; |
| for (int i = -tapsby2 + 1; i <= tapsby2; ++i) { |
| const double tap = kernel((i - x) * p / q, a, win_type); |
| filter[i + tapsby2 - 1] = tap; |
| filter_sum += tap; |
| } |
| assert(filter_sum != 0.0); |
| for (int i = -tapsby2 + 1; i <= tapsby2; ++i) { |
| filter[i + tapsby2 - 1] /= filter_sum; |
| } |
| integerize_array(filter, 2 * tapsby2, bits, ifilter); |
| } |
| |
| static void get_lanczos_upsampler(double x, int p, int q, int a, int bits, |
| WIN_TYPE win_type, int16_t *ifilter) { |
| double filter[MAX_FILTER_LEN] = { 0.0 }; |
| int tapsby2 = get_lanczos_upsampler_filter_length(p, q, a) / 2; |
| assert(tapsby2 * 2 <= MAX_FILTER_LEN); |
| double filter_sum = 0; |
| for (int i = -tapsby2 + 1; i <= tapsby2; ++i) { |
| const double tap = kernel(i - x, a, win_type); |
| filter[i + tapsby2 - 1] = tap; |
| filter_sum += tap; |
| } |
| assert(filter_sum != 0.0); |
| for (int i = -tapsby2 + 1; i <= tapsby2; ++i) { |
| filter[i + tapsby2 - 1] /= filter_sum; |
| } |
| integerize_array(filter, 2 * tapsby2, bits, ifilter); |
| } |
| |
| static int gcd(int p, int q) { |
| int p1 = (p < q ? p : q); |
| int q1 = (p1 == p ? q : p); |
| while (p1) { |
| const int t = p1; |
| p1 = q1 % p1; |
| q1 = t; |
| } |
| return q1; |
| } |
| |
| const char *ext_names[] = { "Repeat", "Symmetric", "Reflect", "Gradient" }; |
| const char *ext2str(EXT_TYPE ext_type) { return ext_names[(int)ext_type]; } |
| |
| int get_resample_filter(int p, int q, int a, double x0, EXT_TYPE ext_type, |
| WIN_TYPE win_type, int subsampled, int bits, |
| RationalResampleFilter *rf) { |
| double offset[MAX_RATIONAL_FACTOR + 1]; |
| int intpel[MAX_RATIONAL_FACTOR]; |
| if (p <= 0 || q <= 0) { |
| fprintf(stderr, "Resampling numerator or denominator must be positive\n"); |
| return 0; |
| } |
| const int g = gcd(p, q); |
| assert(g > 0); |
| rf->p = p / g; |
| rf->q = q / g; |
| if (rf->p <= 0 || rf->p > MAX_RATIONAL_FACTOR) { |
| fprintf(stderr, "Resampling numerator %d ratio exceeds maximum allowed\n", |
| rf->p); |
| return 0; |
| } |
| if (rf->q <= 0 || rf->q > MAX_RATIONAL_FACTOR) { |
| fprintf(stderr, "Resampling denominator %d ratio exceeds maximum allowed\n", |
| rf->q); |
| return 0; |
| } |
| rf->ext_type = ext_type; |
| rf->win_type = win_type; |
| if (x0 == (double)('c')) |
| x0 = get_centered_x0(rf->p, rf->q); |
| else if (x0 == (double)('d')) |
| x0 = subsampled ? get_cosited_chroma_x0(rf->p, rf->q) |
| : get_centered_x0(rf->p, rf->q); |
| rf->filter_bits = bits; |
| for (int i = 0; i < rf->p; ++i) { |
| offset[i] = (double)rf->q / (double)rf->p * i + x0; |
| intpel[i] = (int)floor(offset[i]); |
| rf->phases[i] = offset[i] - intpel[i]; |
| } |
| offset[rf->p] = rf->q + x0; |
| intpel[rf->p] = (int)floor(offset[rf->p]); |
| |
| rf->start = intpel[0]; |
| for (int i = 0; i < rf->p; ++i) rf->steps[i] = intpel[i + 1] - intpel[i]; |
| if (rf->p < rf->q) { // downsampling |
| rf->length = get_lanczos_downsampler_filter_length(rf->p, rf->q, a); |
| if (rf->length > MAX_FILTER_LEN) { |
| fprintf(stderr, "Filter length %d ratio exceeds maximum allowed\n", |
| rf->length); |
| return 0; |
| } |
| for (int i = 0; i < rf->p; ++i) { |
| get_lanczos_downsampler(rf->phases[i], rf->p, rf->q, a, rf->filter_bits, |
| rf->win_type, rf->filter[i]); |
| } |
| } else if (rf->p >= rf->q) { // upsampling |
| rf->length = get_lanczos_upsampler_filter_length(rf->p, rf->q, a); |
| if (rf->length > MAX_FILTER_LEN) { |
| fprintf(stderr, "Filter length %d ratio exceeds maximum allowed\n", |
| rf->length); |
| return 0; |
| } |
| for (int i = 0; i < rf->p; ++i) { |
| get_lanczos_upsampler(rf->phases[i], rf->p, rf->q, a, rf->filter_bits, |
| rf->win_type, rf->filter[i]); |
| } |
| } |
| return 1; |
| } |
| |
| int is_resampler_noop(RationalResampleFilter *rf) { |
| return (rf->p == 1 && rf->q == 1 && rf->phases[0] == 0.0); |
| } |
| |
| int get_resample_filter_inv(int p, int q, int a, double x0, EXT_TYPE ext_type, |
| WIN_TYPE win_type, int subsampled, int bits, |
| RationalResampleFilter *rf) { |
| double y0 = get_inverse_x0(p, q, x0, subsampled); |
| return get_resample_filter(q, p, a, y0, ext_type, win_type, subsampled, bits, |
| rf); |
| } |
| |
| // Assume x buffer is already extended on both sides with x pointing to the |
| // leftmost pixel, and the extension values are already filled up. |
| static void resample_1d_core(const int16_t *x, int inlen, |
| RationalResampleFilter *rf, int downshift, |
| ClipProfile *clip, int16_t *y, int outlen) { |
| (void)inlen; |
| const int tapsby2 = rf->length / 2; |
| const int16_t *xext = x; |
| xext += rf->start; |
| for (int i = 0, p = 0; i < outlen; ++i, p = (p + 1) % rf->p) { |
| int64_t sum = 0; |
| for (int j = -tapsby2 + 1; j <= tapsby2; ++j) { |
| sum += (int)rf->filter[p][j + tapsby2 - 1] * (int)xext[j]; |
| } |
| sum = ROUND_POWER_OF_TWO_SIGNED(sum, downshift); |
| if (clip) { |
| y[i] = (int16_t)(clip->issigned |
| ? doclip((int)sum, -(1 << (clip->bits - 1)), |
| (1 << (clip->bits - 1)) - 1) |
| : doclip((int)sum, 0, (1 << clip->bits) - 1)); |
| } else { |
| y[i] = (int16_t)doclip((int)sum, -(1 << 15), (1 << 15) - 1); |
| } |
| xext += rf->steps[p]; |
| } |
| } |
| |
| static void extend_border(int16_t *x, int inlen, EXT_TYPE ext_type, |
| int border) { |
| switch (ext_type) { |
| case EXT_REPEAT: |
| for (int i = -border; i < 0; ++i) x[i] = x[0]; |
| for (int i = 0; i < border; ++i) x[i + inlen] = x[inlen - 1]; |
| break; |
| case EXT_SYMMETRIC: |
| if (inlen >= border) { |
| for (int i = -border; i < 0; ++i) x[i] = x[-i - 1]; |
| for (int i = 0; i < border; ++i) x[i + inlen] = x[inlen - 1 - i]; |
| } else { |
| for (int i = -border; i < 0; ++i) |
| x[i] = x[(-i - 1 > inlen - 1 ? inlen - 1 : -i - 1)]; |
| for (int i = 0; i < border; ++i) |
| x[i + inlen] = x[(inlen - 1 - i < 0 ? 0 : inlen - 1 - i)]; |
| } |
| break; |
| case EXT_REFLECT: |
| if (inlen > border) { |
| for (int i = -border; i < 0; ++i) x[i] = x[-i]; |
| for (int i = 0; i < border; ++i) x[i + inlen] = x[inlen - 2 - i]; |
| } else { |
| for (int i = -border; i < 0; ++i) |
| x[i] = x[(-i > inlen - 1 ? inlen - 1 : -i)]; |
| for (int i = 0; i < border; ++i) |
| x[i + inlen] = x[(inlen - 2 - i < 0 ? 0 : inlen - 2 - i)]; |
| } |
| break; |
| case EXT_GRADIENT: |
| if (inlen > border) { |
| for (int i = -border; i < 0; ++i) x[i] = 2 * x[0] - x[-i]; |
| for (int i = 0; i < border; ++i) |
| x[i + inlen] = 2 * x[inlen - 1] - x[inlen - 2 - i]; |
| } else { |
| for (int i = -border; i < 0; ++i) |
| x[i] = 2 * x[0] - x[(-i > inlen - 1 ? inlen - 1 : -i)]; |
| for (int i = 0; i < border; ++i) |
| x[i + inlen] = |
| 2 * x[inlen - 1] - x[(inlen - 2 - i < 0 ? 0 : inlen - 2 - i)]; |
| } |
| break; |
| } |
| } |
| |
| // Assume x buffer is already extended on both sides with x pointing to the |
| // leftmost pixel, but the extension values are not filled up. |
| static void resample_1d_xt(int16_t *x, int inlen, RationalResampleFilter *rf, |
| int downshift, ClipProfile *clip, int16_t *y, |
| int outlen) { |
| extend_border(x, inlen, rf->ext_type, rf->length / 2); |
| resample_1d_core(x, inlen, rf, downshift, clip, y, outlen); |
| } |
| |
| // Assume a scratch buffer xext of size inlen + rf->length is provided |
| static void resample_1d_xc(const int16_t *x, int inlen, |
| RationalResampleFilter *rf, int downshift, |
| ClipProfile *clip, int16_t *y, int outlen, |
| int16_t *xext) { |
| memcpy(xext, x, sizeof(*x) * inlen); |
| |
| resample_1d_xt(xext, inlen, rf, downshift, clip, y, outlen); |
| } |
| |
| static void fill_col_to_arr(const int16_t *img, int stride, int len, |
| int16_t *arr) { |
| int i; |
| const int16_t *iptr = img; |
| int16_t *aptr = arr; |
| for (i = 0; i < len; ++i, iptr += stride) { |
| *aptr++ = *iptr; |
| } |
| } |
| |
| static void fill_arr_to_col(int16_t *img, int stride, int len, |
| const int16_t *arr) { |
| int i; |
| int16_t *iptr = img; |
| const int16_t *aptr = arr; |
| for (i = 0; i < len; ++i, iptr += stride) { |
| *iptr = *aptr++; |
| } |
| } |
| |
| void resample_1d(const int16_t *x, int inlen, RationalResampleFilter *rf, |
| int downshift, ClipProfile *clip, int16_t *y, int outlen) { |
| const int tapsby2 = rf->length / 2; |
| int16_t *xext_ = (int16_t *)malloc((inlen + rf->length) * sizeof(*x)); |
| int16_t *xext = xext_ + tapsby2; |
| |
| resample_1d_xc(xext, inlen, rf, downshift, clip, y, outlen, xext); |
| |
| free(xext_); |
| } |
| |
| void resample_2d(const int16_t *x, int inwidth, int inheight, int instride, |
| RationalResampleFilter *rfh, RationalResampleFilter *rfv, |
| int int_extra_bits, ClipProfile *clip, int16_t *y, |
| int outwidth, int outheight, int outstride) { |
| if (rfv == NULL || is_resampler_noop(rfv)) { |
| resample_horz(x, inwidth, inheight, instride, rfh, clip, y, outwidth, |
| outstride); |
| return; |
| } |
| if (rfh == NULL || is_resampler_noop(rfh)) { |
| resample_vert(x, inwidth, inheight, instride, rfv, clip, y, outheight, |
| outstride); |
| return; |
| } |
| int16_t *tmpbuf = (int16_t *)malloc(sizeof(int16_t) * outwidth * inheight); |
| const int arrsize = |
| outheight + ((inheight + rfv->length > inwidth + rfh->length) |
| ? (inheight + rfv->length) |
| : (inwidth + rfh->length)); |
| int16_t *tmparr_ = (int16_t *)calloc(arrsize, sizeof(int16_t)); |
| int16_t *tmparrh = tmparr_ + outheight + rfh->length / 2; |
| int16_t *tmparrv = tmparr_ + outheight + rfv->length / 2; |
| int16_t *tmparro = tmparr_; |
| int tmpstride = outwidth; |
| // intermediate data is stored in 16 bit buffers, so limit int_extra_bits |
| int_extra_bits = MIN(int_extra_bits, 14 - clip->bits); |
| const int downshifth = rfh->filter_bits - int_extra_bits; |
| const int downshiftv = rfh->filter_bits + int_extra_bits; |
| for (int i = 0; i < inheight; ++i) { |
| resample_1d_xc(x + instride * i, inwidth, rfh, downshifth, NULL, |
| tmpbuf + i * tmpstride, outwidth, tmparrh); |
| } |
| for (int i = 0; i < outwidth; ++i) { |
| fill_col_to_arr(tmpbuf + i, outwidth, inheight, tmparrv); |
| resample_1d_xt(tmparrv, inheight, rfv, downshiftv, clip, tmparro, |
| outheight); |
| fill_arr_to_col(y + i, outstride, outheight, tmparro); |
| } |
| free(tmpbuf); |
| free(tmparr_); |
| } |
| |
| void resample_horz(const int16_t *x, int inwidth, int inheight, int instride, |
| RationalResampleFilter *rfh, ClipProfile *clip, int16_t *y, |
| int outwidth, int outstride) { |
| const int arrsize = inwidth + rfh->length; |
| int16_t *tmparr_ = (int16_t *)calloc(arrsize, sizeof(int16_t)); |
| int16_t *tmparrh = tmparr_ + rfh->length / 2; |
| for (int i = 0; i < inheight; ++i) { |
| resample_1d_xc(x + instride * i, inwidth, rfh, rfh->filter_bits, clip, |
| y + i * outstride, outwidth, tmparrh); |
| } |
| free(tmparr_); |
| } |
| |
| void resample_vert(const int16_t *x, int inwidth, int inheight, int instride, |
| RationalResampleFilter *rfv, ClipProfile *clip, int16_t *y, |
| int outheight, int outstride) { |
| const int arrsize = outheight + inheight + rfv->length; |
| int16_t *tmparr_ = (int16_t *)calloc(arrsize, sizeof(int16_t)); |
| int16_t *tmparrv = tmparr_ + outheight + rfv->length / 2; |
| int16_t *tmparro = tmparr_; |
| for (int i = 0; i < inwidth; ++i) { |
| fill_col_to_arr(x + i, instride, inheight, tmparrv); |
| resample_1d_xt(tmparrv, inheight, rfv, rfv->filter_bits, clip, tmparro, |
| outheight); |
| fill_arr_to_col(y + i, outstride, outheight, tmparro); |
| } |
| free(tmparr_); |
| } |
| |
| // Assume x buffer is already extended on both sides with x pointing to the |
| // leftmost pixel, and the extension values are already filled up. |
| static void resample_1d_core_in8b(const uint8_t *x, int inlen, |
| RationalResampleFilter *rf, int downshift, |
| ClipProfile *clip, int16_t *y, int outlen) { |
| (void)inlen; |
| const int tapsby2 = rf->length / 2; |
| const uint8_t *xext = x; |
| xext += rf->start; |
| for (int i = 0, p = 0; i < outlen; ++i, p = (p + 1) % rf->p) { |
| int64_t sum = 0; |
| for (int j = -tapsby2 + 1; j <= tapsby2; ++j) { |
| sum += (int)rf->filter[p][j + tapsby2 - 1] * (int)xext[j]; |
| } |
| sum = ROUND_POWER_OF_TWO_SIGNED(sum, downshift); |
| if (clip) { |
| y[i] = (int16_t)(clip->issigned |
| ? doclip((int)sum, -(1 << (clip->bits - 1)), |
| (1 << (clip->bits - 1)) - 1) |
| : doclip((int)sum, 0, (1 << clip->bits) - 1)); |
| } else { |
| y[i] = (int16_t)doclip((int)sum, -(1 << 15), (1 << 15) - 1); |
| } |
| xext += rf->steps[p]; |
| } |
| } |
| |
| // Assume x buffer is already extended on both sides with x pointing to the |
| // leftmost pixel, and the extension values are already filled up. |
| static void resample_1d_core_8b(const uint8_t *x, int inlen, |
| RationalResampleFilter *rf, int downshift, |
| ClipProfile *clip, uint8_t *y, int outlen) { |
| (void)inlen; |
| const int tapsby2 = rf->length / 2; |
| const uint8_t *xext = x; |
| xext += rf->start; |
| for (int i = 0, p = 0; i < outlen; ++i, p = (p + 1) % rf->p) { |
| int64_t sum = 0; |
| for (int j = -tapsby2 + 1; j <= tapsby2; ++j) { |
| sum += (int)rf->filter[p][j + tapsby2 - 1] * (int)xext[j]; |
| } |
| sum = ROUND_POWER_OF_TWO_SIGNED(sum, downshift); |
| if (clip) { |
| y[i] = (uint8_t)(clip->issigned |
| ? doclip((int)sum, -(1 << (clip->bits - 1)), |
| (1 << (clip->bits - 1)) - 1) |
| : doclip((int)sum, 0, (1 << clip->bits) - 1)); |
| } else { |
| y[i] = (uint8_t)sum; |
| } |
| xext += rf->steps[p]; |
| } |
| } |
| |
| static void extend_border_8b(uint8_t *x, int inlen, EXT_TYPE ext_type, |
| int border) { |
| switch (ext_type) { |
| case EXT_REPEAT: |
| for (int i = -border; i < 0; ++i) x[i] = x[0]; |
| for (int i = 0; i < border; ++i) x[i + inlen] = x[inlen - 1]; |
| break; |
| case EXT_SYMMETRIC: |
| if (inlen >= border) { |
| for (int i = -border; i < 0; ++i) x[i] = x[-i - 1]; |
| for (int i = 0; i < border; ++i) x[i + inlen] = x[inlen - 1 - i]; |
| } else { |
| for (int i = -border; i < 0; ++i) |
| x[i] = x[(-i - 1 > inlen - 1 ? inlen - 1 : -i - 1)]; |
| for (int i = 0; i < border; ++i) |
| x[i + inlen] = x[(inlen - 1 - i < 0 ? 0 : inlen - 1 - i)]; |
| } |
| break; |
| case EXT_REFLECT: |
| if (inlen > border) { |
| for (int i = -border; i < 0; ++i) x[i] = x[-i]; |
| for (int i = 0; i < border; ++i) x[i + inlen] = x[inlen - 2 - i]; |
| } else { |
| for (int i = -border; i < 0; ++i) |
| x[i] = x[(-i > inlen - 1 ? inlen - 1 : -i)]; |
| for (int i = 0; i < border; ++i) |
| x[i + inlen] = x[(inlen - 2 - i < 0 ? 0 : inlen - 2 - i)]; |
| } |
| break; |
| case EXT_GRADIENT: |
| if (inlen > border) { |
| for (int i = -border; i < 0; ++i) { |
| const int t = 2 * x[0] - x[-i]; |
| x[i] = (uint8_t)doclip(t, 0, 255); |
| } |
| for (int i = 0; i < border; ++i) { |
| const int t = 2 * x[inlen - 1] - x[inlen - 2 - i]; |
| x[i + inlen] = (uint8_t)doclip(t, 0, 255); |
| } |
| } else { |
| for (int i = -border; i < 0; ++i) { |
| const int t = 2 * x[0] - x[(-i > inlen - 1 ? inlen - 1 : -i)]; |
| x[i] = (uint8_t)doclip(t, 0, 255); |
| } |
| for (int i = 0; i < border; ++i) { |
| const int t = |
| 2 * x[inlen - 1] - x[(inlen - 2 - i < 0 ? 0 : inlen - 2 - i)]; |
| x[i + inlen] = (uint8_t)doclip(t, 0, 255); |
| } |
| } |
| break; |
| } |
| } |
| |
| static void resample_1d_xt_8b(uint8_t *x, int inlen, RationalResampleFilter *rf, |
| int downshift, ClipProfile *clip, uint8_t *y, |
| int outlen) { |
| extend_border_8b(x, inlen, rf->ext_type, rf->length / 2); |
| resample_1d_core_8b(x, inlen, rf, downshift, clip, y, outlen); |
| } |
| |
| static void resample_1d_xc_8b(const uint8_t *x, int inlen, |
| RationalResampleFilter *rf, int downshift, |
| ClipProfile *clip, uint8_t *y, int outlen, |
| uint8_t *xext) { |
| memcpy(xext, x, inlen * sizeof(*x)); |
| |
| resample_1d_xt_8b(xext, inlen, rf, downshift, clip, y, outlen); |
| } |
| |
| static void resample_1d_xt_in8b(uint8_t *x, int inlen, |
| RationalResampleFilter *rf, int downshift, |
| ClipProfile *clip, int16_t *y, int outlen) { |
| extend_border_8b(x, inlen, rf->ext_type, rf->length / 2); |
| resample_1d_core_in8b(x, inlen, rf, downshift, clip, y, outlen); |
| } |
| |
| static void resample_1d_xc_in8b(const uint8_t *x, int inlen, |
| RationalResampleFilter *rf, int downshift, |
| ClipProfile *clip, int16_t *y, int outlen, |
| uint8_t *xext) { |
| memcpy(xext, x, inlen * sizeof(*x)); |
| |
| resample_1d_xt_in8b(xext, inlen, rf, downshift, clip, y, outlen); |
| } |
| |
| static void fill_col_to_arr_in8b(const uint8_t *img, int stride, int len, |
| int16_t *arr) { |
| int i; |
| const uint8_t *iptr = img; |
| int16_t *aptr = arr; |
| for (i = 0; i < len; ++i, iptr += stride) { |
| *aptr++ = (int16_t)(*iptr); |
| } |
| } |
| |
| static void fill_arr_to_col_out8b(uint8_t *img, int stride, int len, |
| const int16_t *arr) { |
| int i; |
| uint8_t *iptr = img; |
| const int16_t *aptr = arr; |
| for (i = 0; i < len; ++i, iptr += stride) { |
| *iptr = (uint8_t)*aptr++; |
| } |
| } |
| |
| void resample_1d_8b(const uint8_t *x, int inlen, RationalResampleFilter *rf, |
| int downshift, ClipProfile *clip, uint8_t *y, int outlen) { |
| const int tapsby2 = rf->length / 2; |
| uint8_t *xext_ = (uint8_t *)malloc((inlen + rf->length) * sizeof(*x)); |
| uint8_t *xext = xext_ + tapsby2; |
| |
| resample_1d_xc_8b(x, inlen, rf, downshift, clip, y, outlen, xext); |
| |
| free(xext_); |
| } |
| |
| void resample_2d_8b(const uint8_t *x, int inwidth, int inheight, int instride, |
| RationalResampleFilter *rfh, RationalResampleFilter *rfv, |
| int int_extra_bits, ClipProfile *clip, uint8_t *y, |
| int outwidth, int outheight, int outstride) { |
| if (rfv == NULL || is_resampler_noop(rfv)) { |
| resample_horz_8b(x, inwidth, inheight, instride, rfh, clip, y, outwidth, |
| outstride); |
| return; |
| } |
| if (rfh == NULL || is_resampler_noop(rfh)) { |
| resample_vert_8b(x, inwidth, inheight, instride, rfv, clip, y, outheight, |
| outstride); |
| return; |
| } |
| int16_t *tmpbuf = (int16_t *)malloc(sizeof(int16_t) * outwidth * inheight); |
| const int arrsize = |
| outheight + ((inheight + rfv->length > inwidth + rfh->length) |
| ? (inheight + rfv->length) |
| : (inwidth + rfh->length)); |
| int16_t *tmparr_ = (int16_t *)calloc(arrsize, sizeof(int16_t)); |
| int16_t *tmparrh = tmparr_ + outheight + rfh->length / 2; |
| int16_t *tmparrv = tmparr_ + outheight + rfv->length / 2; |
| int16_t *tmparro = tmparr_; |
| int tmpstride = outwidth; |
| // intermediate data is stored in 16 bit buffers, so limit int_extra_bits |
| int_extra_bits = MIN(int_extra_bits, 14 - clip->bits); |
| const int downshifth = rfh->filter_bits - int_extra_bits; |
| const int downshiftv = rfh->filter_bits + int_extra_bits; |
| for (int i = 0; i < inheight; ++i) { |
| resample_1d_xc_in8b(x + instride * i, inwidth, rfh, downshifth, NULL, |
| tmpbuf + i * tmpstride, outwidth, (uint8_t *)tmparrh); |
| } |
| for (int i = 0; i < outwidth; ++i) { |
| fill_col_to_arr(tmpbuf + i, outwidth, inheight, tmparrv); |
| resample_1d_xt(tmparrv, inheight, rfv, downshiftv, clip, tmparro, |
| outheight); |
| fill_arr_to_col_out8b(y + i, outstride, outheight, tmparro); |
| } |
| free(tmpbuf); |
| free(tmparr_); |
| } |
| |
| void resample_horz_8b(const uint8_t *x, int inwidth, int inheight, int instride, |
| RationalResampleFilter *rfh, ClipProfile *clip, |
| uint8_t *y, int outwidth, int outstride) { |
| const int arrsize = inwidth + rfh->length; |
| uint8_t *tmparr_ = (uint8_t *)calloc(arrsize, sizeof(*tmparr_)); |
| uint8_t *tmparrh = tmparr_ + rfh->length / 2; |
| for (int i = 0; i < inheight; ++i) { |
| resample_1d_xc_8b(x + instride * i, inwidth, rfh, rfh->filter_bits, clip, |
| y + i * outstride, outwidth, tmparrh); |
| } |
| free(tmparr_); |
| } |
| |
| void resample_vert_8b(const uint8_t *x, int inwidth, int inheight, int instride, |
| RationalResampleFilter *rfv, ClipProfile *clip, |
| uint8_t *y, int outheight, int outstride) { |
| const int arrsize = outheight + inheight + rfv->length; |
| int16_t *tmparr_ = (int16_t *)calloc(arrsize, sizeof(int16_t)); |
| int16_t *tmparrv = tmparr_ + outheight + rfv->length / 2; |
| int16_t *tmparro = tmparr_; |
| for (int i = 0; i < inwidth; ++i) { |
| fill_col_to_arr_in8b(x + i, instride, inheight, tmparrv); |
| resample_1d_xt(tmparrv, inheight, rfv, rfv->filter_bits, clip, tmparro, |
| outheight); |
| fill_arr_to_col_out8b(y + i, outstride, outheight, tmparro); |
| } |
| free(tmparr_); |
| } |
| |
| int get_resampled_output_length(int inlen, int p, int q, int force_even) { |
| if (!force_even) { |
| // round |
| return (inlen * p + q / 2) / q; |
| } |
| int outlen_floor = inlen * p / q; |
| // choose floor or ceil depending on which one is even |
| if ((outlen_floor % 2) == 1) |
| return outlen_floor + 1; |
| else |
| return outlen_floor; |
| } |