third-party/image-resampler/resampler.cpp

   1 // resampler.cpp, Separable filtering image rescaler v2.21, Rich Geldreich - richgel99@gmail.com
   2 // See unlicense at the bottom of resampler.h, or at http://unlicense.org/
   3 //
   4 // Feb. 1996: Creation, losely based on a heavily bugfixed version of Schumacher's resampler in Graphics Gems 3.
   5 // Oct. 2000: Ported to C++, tweaks.
   6 // May 2001: Continous to discrete mapping, box filter tweaks.
   7 // March 9, 2002: Kaiser filter grabbed from Jonathan Blow's GD magazine mipmap sample code.
   8 // Sept. 8, 2002: Comments cleaned up a bit.
   9 // Dec. 31, 2008: v2.2: Bit more cleanup, released as public domain.
  10 // June 4, 2012: v2.21: Switched to unlicense.org, integrated GCC fixes supplied by Peter Nagy <petern@crytek.com>, Anteru at anteru.net, and clay@coge.net,
  11 // added Codeblocks project (for testing with MinGW and GCC), VS2008 static code analysis pass.
  12 #include <stdlib.h>
  13 #include <math.h>
  14 #include <float.h>
  15 #include <assert.h>
  16 #include <string.h>
  17
  18 // resampler is written using C style casts
  19 #pragma GCC diagnostic push
  20 #pragma GCC diagnostic ignored "-Wold-style-cast"
  21
  22 #include "resampler.h"
  23
  24 #define resampler_assert assert
  25
  26 static inline int resampler_range_check(int v, int h) { (void)h; resampler_assert((v >= 0) && (v < h)); return v; }
  27
  28 #ifndef max
  29    #define max(a,b) (((a) > (b)) ? (a) : (b))
  30 #endif
  31
  32 #ifndef min
  33    #define min(a,b) (((a) < (b)) ? (a) : (b))
  34 #endif
  35
  36 #ifndef TRUE
  37    #define TRUE (1)
  38 #endif
  39
  40 #ifndef FALSE
  41    #define FALSE (0)
  42 #endif
  43
  44 #define RESAMPLER_DEBUG 0
  45
  46 #define M_PI 3.14159265358979323846
  47
  48 // Float to int cast with truncation.
  49 static inline int cast_to_int(Resample_Real i)
  50 {
  51    return (int)i;
  52 }
  53
  54 // (x mod y) with special handling for negative x values.
  55 static inline int posmod(int x, int y)
  56 {
  57    if (x >= 0)
  58       return (x % y);
  59    else
  60    {
  61       int m = (-x) % y;
  62
  63       if (m != 0)
  64          m = y - m;
  65
  66       return (m);
  67    }
  68 }
  69
  70 // To add your own filter, insert the new function below and update the filter table.
  71 // There is no need to make the filter function particularly fast, because it's
  72 // only called during initializing to create the X and Y axis contributor tables.
  73
  74 #define BOX_FILTER_SUPPORT (0.5f)
  75 static Resample_Real box_filter(Resample_Real t)    /* pulse/Fourier window */
  76 {
  77    // make_clist() calls the filter function with t inverted (pos = left, neg = right)
  78    if ((t >= -0.5f) && (t < 0.5f))
  79       return 1.0f;
  80    else
  81       return 0.0f;
  82 }
  83
  84 #define TENT_FILTER_SUPPORT (1.0f)
  85 static Resample_Real tent_filter(Resample_Real t)   /* box (*) box, bilinear/triangle */
  86 {
  87    if (t < 0.0f)
  88       t = -t;
  89
  90    if (t < 1.0f)
  91       return 1.0f - t;
  92    else
  93       return 0.0f;
  94 }
  95
  96 #define BELL_SUPPORT (1.5f)
  97 static Resample_Real bell_filter(Resample_Real t)    /* box (*) box (*) box */
  98 {
  99    if (t < 0.0f)
 100       t = -t;
 101
 102    if (t < .5f)
 103       return (.75f - (t * t));
 104
 105    if (t < 1.5f)
 106    {
 107       t = (t - 1.5f);
 108       return (.5f * (t * t));
 109    }
 110
 111    return (0.0f);
 112 }
 113
 114 #define B_SPLINE_SUPPORT (2.0f)
 115 static Resample_Real B_spline_filter(Resample_Real t)  /* box (*) box (*) box (*) box */
 116 {
 117    Resample_Real tt;
 118
 119    if (t < 0.0f)
 120       t = -t;
 121
 122    if (t < 1.0f)
 123    {
 124       tt = t * t;
 125       return ((.5f * tt * t) - tt + (2.0f / 3.0f));
 126    }
 127    else if (t < 2.0f)
 128    {
 129       t = 2.0f - t;
 130       return ((1.0f / 6.0f) * (t * t * t));
 131    }
 132
 133    return (0.0f);
 134 }
 135
 136 // Dodgson, N., "Quadratic Interpolation for Image Resampling"
 137 #define QUADRATIC_SUPPORT 1.5f
 138 static Resample_Real quadratic(Resample_Real t, const Resample_Real R)
 139 {
 140    if (t < 0.0f)
 141       t = -t;
 142    if (t < QUADRATIC_SUPPORT)
 143    {
 144       Resample_Real tt = t * t;
 145       if (t <= .5f)
 146          return (-2.0f * R) * tt + .5f * (R + 1.0f);
 147       else
 148          return (R * tt) + (-2.0f * R - .5f) * t + (3.0f / 4.0f) * (R + 1.0f);
 149    }
 150    else
 151       return 0.0f;
 152 }
 153
 154 static Resample_Real quadratic_interp_filter(Resample_Real t)
 155 {
 156    return quadratic(t, 1.0f);
 157 }
 158
 159 static Resample_Real quadratic_approx_filter(Resample_Real t)
 160 {
 161    return quadratic(t, .5f);
 162 }
 163
 164 static Resample_Real quadratic_mix_filter(Resample_Real t)
 165 {
 166    return quadratic(t, .8f);
 167 }
 168
 169 // Mitchell, D. and A. Netravali, "Reconstruction Filters in Computer Graphics."
 170 // Computer Graphics, Vol. 22, No. 4, pp. 221-228.
 171 // (B, C)
 172 // (1/3, 1/3) - Defaults recommended by Mitchell and Netravali
 173 // (1, 0)     - Equivalent to the Cubic B-Spline
 174 // (0, 0.5)   - Equivalent to the Catmull-Rom Spline
 175 // (0, C)     - The family of Cardinal Cubic Splines
 176 // (B, 0)     - Duff's tensioned B-Splines.
 177 static Resample_Real mitchell(Resample_Real t, const Resample_Real B, const Resample_Real C)
 178 {
 179    Resample_Real tt;
 180
 181    tt = t * t;
 182
 183    if(t < 0.0f)
 184       t = -t;
 185
 186    if(t < 1.0f)
 187    {
 188       t = (((12.0f - 9.0f * B - 6.0f * C) * (t * tt))
 189          + ((-18.0f + 12.0f * B + 6.0f * C) * tt)
 190          + (6.0f - 2.0f * B));
 191
 192       return (t / 6.0f);
 193    }
 194    else if (t < 2.0f)
 195    {
 196       t = (((-1.0f * B - 6.0f * C) * (t * tt))
 197          + ((6.0f * B + 30.0f * C) * tt)
 198          + ((-12.0f * B - 48.0f * C) * t)
 199          + (8.0f * B + 24.0f * C));
 200
 201       return (t / 6.0f);
 202    }
 203
 204    return (0.0f);
 205 }
 206
 207 #define MITCHELL_SUPPORT (2.0f)
 208 static Resample_Real mitchell_filter(Resample_Real t)
 209 {
 210    return mitchell(t, 1.0f / 3.0f, 1.0f / 3.0f);
 211 }
 212
 213 #define CATMULL_ROM_SUPPORT (2.0f)
 214 static Resample_Real catmull_rom_filter(Resample_Real t)
 215 {
 216    return mitchell(t, 0.0f, .5f);
 217 }
 218
 219 static double sinc(double x)
 220 {
 221    x = (x * M_PI);
 222
 223    if ((x < 0.01f) && (x > -0.01f))
 224       return 1.0f + x*x*(-1.0f/6.0f + x*x*1.0f/120.0f);
 225
 226    return sin(x) / x;
 227 }
 228
 229 static Resample_Real clean(double t)
 230 {
 231    const Resample_Real EPSILON = .0000125f;
 232    if (fabs(t) < EPSILON)
 233       return 0.0f;
 234    return (Resample_Real)t;
 235 }
 236
 237 static double blackman_exact_window(double x)
 238 {
 239    return 0.42659071f + 0.49656062f * cos(M_PI*x) + 0.07684867f * cos(2.0f*M_PI*x);
 240 }
 241
 242 #define BLACKMAN_SUPPORT (3.0f)
 243 static Resample_Real blackman_filter(Resample_Real t)
 244 {
 245    if (t < 0.0f)
 246       t = -t;
 247
 248    if (t < 3.0f)
 249       //return clean(sinc(t) * blackman_window(t / 3.0f));
 250       return clean(sinc(t) * blackman_exact_window(t / 3.0f));
 251    else
 252       return (0.0f);
 253 }
 254
 255 #define GAUSSIAN_SUPPORT (1.25f)
 256 static Resample_Real gaussian_filter(Resample_Real t) // with blackman window
 257 {
 258    if (t < 0)
 259       t = -t;
 260    if (t < GAUSSIAN_SUPPORT)
 261       return clean(exp(-2.0f * t * t) * sqrt(2.0f / M_PI) * blackman_exact_window(t / GAUSSIAN_SUPPORT));
 262    else
 263       return 0.0f;
 264 }
 265
 266 // Windowed sinc -- see "Jimm Blinn's Corner: Dirty Pixels" pg. 26.
 267 #define LANCZOS3_SUPPORT (3.0f)
 268 static Resample_Real lanczos3_filter(Resample_Real t)
 269 {
 270    if (t < 0.0f)
 271       t = -t;
 272
 273    if (t < 3.0f)
 274       return clean(sinc(t) * sinc(t / 3.0f));
 275    else
 276       return (0.0f);
 277 }
 278
 279 #define LANCZOS4_SUPPORT (4.0f)
 280 static Resample_Real lanczos4_filter(Resample_Real t)
 281 {
 282    if (t < 0.0f)
 283       t = -t;
 284
 285    if (t < 4.0f)
 286       return clean(sinc(t) * sinc(t / 4.0f));
 287    else
 288       return (0.0f);
 289 }
 290
 291 #define LANCZOS6_SUPPORT (6.0f)
 292 static Resample_Real lanczos6_filter(Resample_Real t)
 293 {
 294    if (t < 0.0f)
 295       t = -t;
 296
 297    if (t < 6.0f)
 298       return clean(sinc(t) * sinc(t / 6.0f));
 299    else
 300       return (0.0f);
 301 }
 302
 303 #define LANCZOS12_SUPPORT (12.0f)
 304 static Resample_Real lanczos12_filter(Resample_Real t)
 305 {
 306    if (t < 0.0f)
 307       t = -t;
 308
 309    if (t < 12.0f)
 310       return clean(sinc(t) * sinc(t / 12.0f));
 311    else
 312       return (0.0f);
 313 }
 314
 315 static double bessel0(double x)
 316 {
 317    const double EPSILON_RATIO = 1E-16;
 318    double xh, sum, pow, ds;
 319    int k;
 320
 321    xh = 0.5 * x;
 322    sum = 1.0;
 323    pow = 1.0;
 324    k = 0;
 325    ds = 1.0;
 326    while (ds > sum * EPSILON_RATIO) // FIXME: Shouldn't this stop after X iterations for max. safety?
 327    {
 328       ++k;
 329       pow = pow * (xh / k);
 330       ds = pow * pow;
 331       sum = sum + ds;
 332    }
 333
 334    return sum;
 335 }
 336
 337 static const Resample_Real KAISER_ALPHA = 4.0;
 338 static double kaiser(double alpha, double half_width, double x)
 339 {
 340    const double ratio = (x / half_width);
 341    return bessel0(alpha * sqrt(1 - ratio * ratio)) / bessel0(alpha);
 342 }
 343
 344 #define KAISER_SUPPORT 3
 345 static Resample_Real kaiser_filter(Resample_Real t)
 346 {
 347    if (t < 0.0f)
 348       t = -t;
 349
 350    if (t < KAISER_SUPPORT)
 351    {
 352       // db atten
 353       const Resample_Real att = 40.0f;
 354       const Resample_Real alpha = (Resample_Real)(exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886 * (att - 20.96));
 355       //const Resample_Real alpha = KAISER_ALPHA;
 356       return (Resample_Real)clean(sinc(t) * kaiser(alpha, KAISER_SUPPORT, t));
 357    }
 358
 359    return 0.0f;
 360 }
 361
 362 // filters[] is a list of all the available filter functions.
 363 static struct
 364 {
 365    Resampler::Filter name;
 366    Resample_Real (*func)(Resample_Real t);
 367    Resample_Real support;
 368 } g_filters[] =
 369 {
 370    { Resampler::BOX,              box_filter,              BOX_FILTER_SUPPORT },
 371    { Resampler::TENT,             tent_filter,             TENT_FILTER_SUPPORT },
 372    { Resampler::BELL,             bell_filter,             BELL_SUPPORT },
 373    { Resampler::B_SPLINE,         B_spline_filter,         B_SPLINE_SUPPORT },
 374    { Resampler::MITCHELL,         mitchell_filter,         MITCHELL_SUPPORT },
 375    { Resampler::LANCZOS3,         lanczos3_filter,         LANCZOS3_SUPPORT },
 376    { Resampler::BLACKMAN,         blackman_filter,         BLACKMAN_SUPPORT },
 377    { Resampler::LANCZOS4,         lanczos4_filter,         LANCZOS4_SUPPORT },
 378    { Resampler::LANCZOS6,         lanczos6_filter,         LANCZOS6_SUPPORT },
 379    { Resampler::LANCZOS12,        lanczos12_filter,        LANCZOS12_SUPPORT },
 380    { Resampler::KAISER,           kaiser_filter,           KAISER_SUPPORT },
 381    { Resampler::GAUSSIAN,         gaussian_filter,         GAUSSIAN_SUPPORT },
 382    { Resampler::CATMULLROM,       catmull_rom_filter,      CATMULL_ROM_SUPPORT },
 383    { Resampler::QUADRATIC_INTERPOLATION, quadratic_interp_filter, QUADRATIC_SUPPORT },
 384    { Resampler::QUADRATIC_APPROXIMATION, quadratic_approx_filter, QUADRATIC_SUPPORT },
 385    { Resampler::QUADRATIC_MIX,    quadratic_mix_filter,    QUADRATIC_SUPPORT },
 386 };
 387
 388 static const int NUM_FILTERS = sizeof(g_filters) / sizeof(g_filters[0]);
 389
 390 /* Ensure that the contributing source sample is
 391 * within bounds. If not, reflect, clamp, or wrap.
 392 */
 393 int Resampler::reflect(const int j, const int src_x, const Boundary_Op boundary_op)
 394 {
 395    int n;
 396
 397    if (j < 0)
 398    {
 399       if (boundary_op == BOUNDARY_REFLECT)
 400       {
 401          n = -j;
 402
 403          if (n >= src_x)
 404             n = src_x - 1;
 405       }
 406       else if (boundary_op == BOUNDARY_WRAP)
 407          n = posmod(j, src_x);
 408       else
 409          n = 0;
 410    }
 411    else if (j >= src_x)
 412    {
 413       if (boundary_op == BOUNDARY_REFLECT)
 414       {
 415          n = (src_x - j) + (src_x - 1);
 416
 417          if (n < 0)
 418             n = 0;
 419       }
 420       else if (boundary_op == BOUNDARY_WRAP)
 421          n = posmod(j, src_x);
 422       else
 423          n = src_x - 1;
 424    }
 425    else
 426       n = j;
 427
 428    return n;
 429 }
 430
 431 // The make_clist() method generates, for all destination samples,
 432 // the list of all source samples with non-zero weighted contributions.
 433 Resampler::Contrib_List* Resampler::make_clist(
 434    int src_x, int dst_x, Boundary_Op boundary_op,
 435    Resample_Real (*Pfilter)(Resample_Real),
 436    Resample_Real filter_support,
 437    Resample_Real filter_scale,
 438    Resample_Real src_ofs)
 439 {
 440    typedef struct
 441    {
 442       // The center of the range in DISCRETE coordinates (pixel center = 0.0f).
 443       Resample_Real center;
 444       int left, right;
 445    } Contrib_Bounds;
 446
 447    int i, j, k, n, left, right;
 448    Resample_Real total_weight;
 449    Resample_Real xscale, center, half_width, weight;
 450    Contrib_List* Pcontrib;
 451    Contrib* Pcpool;
 452    Contrib* Pcpool_next;
 453    Contrib_Bounds* Pcontrib_bounds;
 454
 455    if ((Pcontrib = (Contrib_List*)calloc(dst_x, sizeof(Contrib_List))) == NULL)
 456       return NULL;
 457
 458    Pcontrib_bounds = (Contrib_Bounds*)calloc(dst_x, sizeof(Contrib_Bounds));
 459    if (!Pcontrib_bounds)
 460    {
 461       free(Pcontrib);
 462       return (NULL);
 463    }
 464
 465    const Resample_Real oo_filter_scale = 1.0f / filter_scale;
 466
 467    const Resample_Real NUDGE = 0.5f;
 468    xscale = dst_x / (Resample_Real)src_x;
 469
 470    if (xscale < 1.0f)
 471    {
 472       int total; (void)total;
 473
 474       /* Handle case when there are fewer destination
 475       * samples than source samples (downsampling/minification).
 476       */
 477
 478       // stretched half width of filter
 479       half_width = (filter_support / xscale) * filter_scale;
 480
 481       // Find the range of source sample(s) that will contribute to each destination sample.
 482
 483       for (i = 0, n = 0; i < dst_x; i++)
 484       {
 485          // Convert from discrete to continuous coordinates, scale, then convert back to discrete.
 486          center = ((Resample_Real)i + NUDGE) / xscale;
 487          center -= NUDGE;
 488          center += src_ofs;
 489
 490          left   = cast_to_int((Resample_Real)floor(center - half_width));
 491          right  = cast_to_int((Resample_Real)ceil(center + half_width));
 492
 493          Pcontrib_bounds[i].center = center;
 494          Pcontrib_bounds[i].left = left;
 495          Pcontrib_bounds[i].right = right;
 496
 497          n += (right - left + 1);
 498       }
 499
 500       /* Allocate memory for contributors. */
 501
 502       if ((n == 0) || ((Pcpool = (Contrib*)calloc(n, sizeof(Contrib))) == NULL))
 503       {
 504          free(Pcontrib);
 505          free(Pcontrib_bounds);
 506          return NULL;
 507       }
 508       total = n;
 509
 510       Pcpool_next = Pcpool;
 511
 512       /* Create the list of source samples which
 513       * contribute to each destination sample.
 514       */
 515
 516       for (i = 0; i < dst_x; i++)
 517       {
 518          int max_k = -1;
 519          Resample_Real max_w = -1e+20f;
 520
 521          center = Pcontrib_bounds[i].center;
 522          left   = Pcontrib_bounds[i].left;
 523          right  = Pcontrib_bounds[i].right;
 524
 525          Pcontrib[i].n = 0;
 526          Pcontrib[i].p = Pcpool_next;
 527          Pcpool_next += (right - left + 1);
 528          resampler_assert ((Pcpool_next - Pcpool) <= total);
 529
 530          total_weight = 0;
 531
 532          for (j = left; j <= right; j++)
 533             total_weight += (*Pfilter)((center - (Resample_Real)j) * xscale * oo_filter_scale);
 534          const Resample_Real norm = static_cast<Resample_Real>(1.0f / total_weight);
 535
 536          total_weight = 0;
 537
 538 #if RESAMPLER_DEBUG
 539          printf("%i: ", i);
 540 #endif
 541
 542          for (j = left; j <= right; j++)
 543          {
 544             weight = (*Pfilter)((center - (Resample_Real)j) * xscale * oo_filter_scale) * norm;
 545             if (weight == 0.0f)
 546                continue;
 547
 548             n = reflect(j, src_x, boundary_op);
 549
 550 #if RESAMPLER_DEBUG
 551             printf("%i(%f), ", n, weight);
 552 #endif
 553
 554             /* Increment the number of source
 555             * samples which contribute to the
 556             * current destination sample.
 557             */
 558
 559             k = Pcontrib[i].n++;
 560
 561             Pcontrib[i].p[k].pixel  = (unsigned short)(n);       /* store src sample number */
 562             Pcontrib[i].p[k].weight = weight; /* store src sample weight */
 563
 564             total_weight += weight;          /* total weight of all contributors */
 565
 566             if (weight > max_w)
 567             {
 568                max_w = weight;
 569                max_k = k;
 570             }
 571          }
 572
 573 #if RESAMPLER_DEBUG
 574          printf("\n\n");
 575 #endif
 576
 577          //resampler_assert(Pcontrib[i].n);
 578          //resampler_assert(max_k != -1);
 579          if ((max_k == -1) || (Pcontrib[i].n == 0))
 580          {
 581             free(Pcpool);
 582             free(Pcontrib);
 583             free(Pcontrib_bounds);
 584             return NULL;
 585          }
 586
 587          if (total_weight != 1.0f)
 588             Pcontrib[i].p[max_k].weight += 1.0f - total_weight;
 589       }
 590    }
 591    else
 592    {
 593       /* Handle case when there are more
 594       * destination samples than source
 595       * samples (upsampling).
 596       */
 597
 598       half_width = filter_support * filter_scale;
 599
 600       // Find the source sample(s) that contribute to each destination sample.
 601
 602       for (i = 0, n = 0; i < dst_x; i++)
 603       {
 604          // Convert from discrete to continuous coordinates, scale, then convert back to discrete.
 605          center = ((Resample_Real)i + NUDGE) / xscale;
 606          center -= NUDGE;
 607          center += src_ofs;
 608
 609          left   = cast_to_int((Resample_Real)floor(center - half_width));
 610          right  = cast_to_int((Resample_Real)ceil(center + half_width));
 611
 612          Pcontrib_bounds[i].center = center;
 613          Pcontrib_bounds[i].left = left;
 614          Pcontrib_bounds[i].right = right;
 615
 616          n += (right - left + 1);
 617       }
 618
 619       /* Allocate memory for contributors. */
 620
 621       int total = n;
 622       if ((total == 0) || ((Pcpool = (Contrib*)calloc(total, sizeof(Contrib))) == NULL))
 623       {
 624          free(Pcontrib);
 625          free(Pcontrib_bounds);
 626          return NULL;
 627       }
 628
 629       Pcpool_next = Pcpool;
 630
 631       /* Create the list of source samples which
 632       * contribute to each destination sample.
 633       */
 634
 635       for (i = 0; i < dst_x; i++)
 636       {
 637          int max_k = -1;
 638          Resample_Real max_w = -1e+20f;
 639
 640          center = Pcontrib_bounds[i].center;
 641          left   = Pcontrib_bounds[i].left;
 642          right  = Pcontrib_bounds[i].right;
 643
 644          Pcontrib[i].n = 0;
 645          Pcontrib[i].p = Pcpool_next;
 646          Pcpool_next += (right - left + 1);
 647          resampler_assert((Pcpool_next - Pcpool) <= total);
 648
 649          total_weight = 0;
 650          for (j = left; j <= right; j++)
 651             total_weight += (*Pfilter)((center - (Resample_Real)j) * oo_filter_scale);
 652
 653          const Resample_Real norm = static_cast<Resample_Real>(1.0f / total_weight);
 654
 655          total_weight = 0;
 656
 657 #if RESAMPLER_DEBUG
 658          printf("%i: ", i);
 659 #endif
 660
 661          for (j = left; j <= right; j++)
 662          {
 663             weight = (*Pfilter)((center - (Resample_Real)j) * oo_filter_scale) * norm;
 664             if (weight == 0.0f)
 665                continue;
 666
 667             n = reflect(j, src_x, boundary_op);
 668
 669 #if RESAMPLER_DEBUG
 670             printf("%i(%f), ", n, weight);
 671 #endif
 672
 673             /* Increment the number of source
 674             * samples which contribute to the
 675             * current destination sample.
 676             */
 677
 678             k = Pcontrib[i].n++;
 679
 680             Pcontrib[i].p[k].pixel  = (unsigned short)(n);       /* store src sample number */
 681             Pcontrib[i].p[k].weight = weight; /* store src sample weight */
 682
 683             total_weight += weight;          /* total weight of all contributors */
 684
 685             if (weight > max_w)
 686             {
 687                max_w = weight;
 688                max_k = k;
 689             }
 690          }
 691
 692 #if RESAMPLER_DEBUG
 693          printf("\n\n");
 694 #endif
 695
 696          //resampler_assert(Pcontrib[i].n);
 697          //resampler_assert(max_k != -1);
 698
 699          if ((max_k == -1) || (Pcontrib[i].n == 0))
 700          {
 701             free(Pcpool);
 702             free(Pcontrib);
 703             free(Pcontrib_bounds);
 704             return NULL;
 705          }
 706
 707          if (total_weight != 1.0f)
 708             Pcontrib[i].p[max_k].weight += 1.0f - total_weight;
 709       }
 710    }
 711
 712 #if RESAMPLER_DEBUG
 713    printf("*******\n");
 714 #endif
 715
 716    free(Pcontrib_bounds);
 717
 718    return Pcontrib;
 719 }
 720
 721 void Resampler::resample_x(Sample* Pdst, const Sample* Psrc)
 722 {
 723    resampler_assert(Pdst);
 724    resampler_assert(Psrc);
 725
 726    int i, j;
 727    Sample total;
 728    Contrib_List *Pclist = m_Pclist_x;
 729    Contrib *p;
 730
 731    for (i = m_resample_dst_x; i > 0; i--, Pclist++)
 732    {
 733 #if RESAMPLER_DEBUG_OPS
 734       total_ops += Pclist->n;
 735 #endif
 736
 737       for (j = Pclist->n, p = Pclist->p, total = 0; j > 0; j--, p++)
 738          total += Psrc[p->pixel] * p->weight;
 739
 740       *Pdst++ = total;
 741    }
 742 }
 743
 744 void Resampler::scale_y_mov(Sample* Ptmp, const Sample* Psrc, Resample_Real weight, int dst_x)
 745 {
 746    int i;
 747
 748 #if RESAMPLER_DEBUG_OPS
 749    total_ops += dst_x;
 750 #endif
 751
 752    // Not += because temp buf wasn't cleared.
 753    for (i = dst_x; i > 0; i--)
 754       *Ptmp++ = *Psrc++ * weight;
 755 }
 756
 757 void Resampler::scale_y_add(Sample* Ptmp, const Sample* Psrc, Resample_Real weight, int dst_x)
 758 {
 759 #if RESAMPLER_DEBUG_OPS
 760    total_ops += dst_x;
 761 #endif
 762
 763    for (int i = dst_x; i > 0; i--)
 764       (*Ptmp++) += *Psrc++ * weight;
 765 }
 766
 767 void Resampler::clamp(Sample* Pdst, int n)
 768 {
 769    while (n > 0)
 770    {
 771       *Pdst = clamp_sample(*Pdst);
 772       ++Pdst;
 773       n--;
 774    }
 775 }
 776
 777 void Resampler::resample_y(Sample* Pdst)
 778 {
 779    int i, j;
 780    Sample* Psrc;
 781    Contrib_List* Pclist = &m_Pclist_y[m_cur_dst_y];
 782
 783    Sample* Ptmp = m_delay_x_resample ? m_Ptmp_buf : Pdst;
 784    resampler_assert(Ptmp);
 785
 786    /* Process each contributor. */
 787
 788    for (i = 0; i < Pclist->n; i++)
 789    {
 790       /* locate the contributor's location in the scan
 791       * buffer -- the contributor must always be found!
 792       */
 793
 794       for (j = 0; j < MAX_SCAN_BUF_SIZE; j++)
 795          if (m_Pscan_buf->scan_buf_y[j] == Pclist->p[i].pixel)
 796             break;
 797
 798       resampler_assert(j < MAX_SCAN_BUF_SIZE);
 799
 800       Psrc = m_Pscan_buf->scan_buf_l[j];
 801
 802       if (!i)
 803          scale_y_mov(Ptmp, Psrc, Pclist->p[i].weight, m_intermediate_x);
 804       else
 805          scale_y_add(Ptmp, Psrc, Pclist->p[i].weight, m_intermediate_x);
 806
 807       /* If this source line doesn't contribute to any
 808       * more destination lines then mark the scanline buffer slot
 809       * which holds this source line as free.
 810       * (The max. number of slots used depends on the Y
 811       * axis sampling factor and the scaled filter width.)
 812       */
 813
 814       if (--m_Psrc_y_count[resampler_range_check(Pclist->p[i].pixel, m_resample_src_y)] == 0)
 815       {
 816          m_Psrc_y_flag[resampler_range_check(Pclist->p[i].pixel, m_resample_src_y)] = FALSE;
 817          m_Pscan_buf->scan_buf_y[j] = -1;
 818       }
 819    }
 820
 821    /* Now generate the destination line */
 822
 823    if (m_delay_x_resample) // Was X resampling delayed until after Y resampling?
 824    {
 825       resampler_assert(Pdst != Ptmp);
 826       resample_x(Pdst, Ptmp);
 827    }
 828    else
 829    {
 830       resampler_assert(Pdst == Ptmp);
 831    }
 832
 833    if (m_lo < m_hi)
 834       clamp(Pdst, m_resample_dst_x);
 835 }
 836
 837 bool Resampler::put_line(const Sample* Psrc)
 838 {
 839    int i;
 840
 841    if (m_cur_src_y >= m_resample_src_y)
 842       return false;
 843
 844    /* Does this source line contribute
 845    * to any destination line? if not,
 846    * exit now.
 847    */
 848
 849    if (!m_Psrc_y_count[resampler_range_check(m_cur_src_y, m_resample_src_y)])
 850    {
 851       m_cur_src_y++;
 852       return true;
 853    }
 854
 855    /* Find an empty slot in the scanline buffer. (FIXME: Perf. is terrible here with extreme scaling ratios.) */
 856
 857    for (i = 0; i < MAX_SCAN_BUF_SIZE; i++)
 858       if (m_Pscan_buf->scan_buf_y[i] == -1)
 859          break;
 860
 861    /* If the buffer is full, exit with an error. */
 862
 863    if (i == MAX_SCAN_BUF_SIZE)
 864    {
 865       m_status = STATUS_SCAN_BUFFER_FULL;
 866       return false;
 867    }
 868
 869    m_Psrc_y_flag[resampler_range_check(m_cur_src_y, m_resample_src_y)] = TRUE;
 870    m_Pscan_buf->scan_buf_y[i]  = m_cur_src_y;
 871
 872    /* Does this slot have any memory allocated to it? */
 873
 874    if (!m_Pscan_buf->scan_buf_l[i])
 875    {
 876       if ((m_Pscan_buf->scan_buf_l[i] = (Sample*)malloc(m_intermediate_x * sizeof(Sample))) == NULL)
 877       {
 878          m_status = STATUS_OUT_OF_MEMORY;
 879          return false;
 880       }
 881    }
 882
 883    // Resampling on the X axis first?
 884    if (m_delay_x_resample)
 885    {
 886       resampler_assert(m_intermediate_x == m_resample_src_x);
 887
 888       // Y-X resampling order
 889       memcpy(m_Pscan_buf->scan_buf_l[i], Psrc, m_intermediate_x * sizeof(Sample));
 890    }
 891    else
 892    {
 893       resampler_assert(m_intermediate_x == m_resample_dst_x);
 894
 895       // X-Y resampling order
 896       resample_x(m_Pscan_buf->scan_buf_l[i], Psrc);
 897    }
 898
 899    m_cur_src_y++;
 900
 901    return true;
 902 }
 903
 904 const Resampler::Sample* Resampler::get_line()
 905 {
 906    int i;
 907
 908    /* If all the destination lines have been
 909    * generated, then always return NULL.
 910    */
 911
 912    if (m_cur_dst_y == m_resample_dst_y)
 913       return NULL;
 914
 915    /* Check to see if all the required
 916    * contributors are present, if not,
 917    * return NULL.
 918    */
 919
 920    for (i = 0; i < m_Pclist_y[m_cur_dst_y].n; i++)
 921       if (!m_Psrc_y_flag[resampler_range_check(m_Pclist_y[m_cur_dst_y].p[i].pixel, m_resample_src_y)])
 922          return NULL;
 923
 924    resample_y(m_Pdst_buf);
 925
 926    m_cur_dst_y++;
 927
 928    return m_Pdst_buf;
 929 }
 930
 931 Resampler::~Resampler()
 932 {
 933    int i;
 934
 935 #if RESAMPLER_DEBUG_OPS
 936    printf("actual ops: %i\n", total_ops);
 937 #endif
 938
 939    free(m_Pdst_buf);
 940    m_Pdst_buf = NULL;
 941
 942    if (m_Ptmp_buf)
 943    {
 944       free(m_Ptmp_buf);
 945       m_Ptmp_buf = NULL;
 946    }
 947
 948    /* Don't deallocate a contibutor list
 949    * if the user passed us one of their own.
 950    */
 951
 952    if ((m_Pclist_x) && (!m_clist_x_forced))
 953    {
 954       free(m_Pclist_x->p);
 955       free(m_Pclist_x);
 956       m_Pclist_x = NULL;
 957    }
 958
 959    if ((m_Pclist_y) && (!m_clist_y_forced))
 960    {
 961       free(m_Pclist_y->p);
 962       free(m_Pclist_y);
 963       m_Pclist_y = NULL;
 964    }
 965
 966    free(m_Psrc_y_count);
 967    m_Psrc_y_count = NULL;
 968
 969    free(m_Psrc_y_flag);
 970    m_Psrc_y_flag = NULL;
 971
 972    if (m_Pscan_buf)
 973    {
 974       for (i = 0; i < MAX_SCAN_BUF_SIZE; i++)
 975          free(m_Pscan_buf->scan_buf_l[i]);
 976
 977       free(m_Pscan_buf);
 978       m_Pscan_buf = NULL;
 979    }
 980 }
 981
 982 void Resampler::restart()
 983 {
 984    if (STATUS_OKAY != m_status)
 985       return;
 986
 987    m_cur_src_y = m_cur_dst_y = 0;
 988
 989    int i, j;
 990    for (i = 0; i < m_resample_src_y; i++)
 991    {
 992       m_Psrc_y_count[i] = 0;
 993       m_Psrc_y_flag[i] = FALSE;
 994    }
 995
 996    for (i = 0; i < m_resample_dst_y; i++)
 997    {
 998       for (j = 0; j < m_Pclist_y[i].n; j++)
 999          m_Psrc_y_count[resampler_range_check(m_Pclist_y[i].p[j].pixel, m_resample_src_y)]++;
1000    }
1001
1002    for (i = 0; i < MAX_SCAN_BUF_SIZE; i++)
1003    {
1004       m_Pscan_buf->scan_buf_y[i] = -1;
1005
1006       free(m_Pscan_buf->scan_buf_l[i]);
1007       m_Pscan_buf->scan_buf_l[i] = NULL;
1008    }
1009 }
1010
1011 Resampler::Resampler(int src_x, int src_y,
1012                      int dst_x, int dst_y,
1013                      Boundary_Op boundary_op,
1014                      Resample_Real sample_low, Resample_Real sample_high,
1015                      Resampler::Filter filter,
1016                      Contrib_List* Pclist_x,
1017                      Contrib_List* Pclist_y,
1018                      Resample_Real filter_x_scale,
1019                      Resample_Real filter_y_scale,
1020                      Resample_Real src_x_ofs,
1021                      Resample_Real src_y_ofs)
1022 {
1023    int i, j;
1024    Resample_Real support, (*func)(Resample_Real);
1025
1026    resampler_assert(src_x > 0);
1027    resampler_assert(src_y > 0);
1028    resampler_assert(dst_x > 0);
1029    resampler_assert(dst_y > 0);
1030
1031 #if RESAMPLER_DEBUG_OPS
1032    total_ops = 0;
1033 #endif
1034
1035    m_lo = sample_low;
1036    m_hi = sample_high;
1037
1038    m_delay_x_resample = false;
1039    m_intermediate_x = 0;
1040    m_Pdst_buf = NULL;
1041    m_Ptmp_buf = NULL;
1042    m_clist_x_forced = false;
1043    m_Pclist_x = NULL;
1044    m_clist_y_forced = false;
1045    m_Pclist_y = NULL;
1046    m_Psrc_y_count = NULL;
1047    m_Psrc_y_flag = NULL;
1048    m_Pscan_buf = NULL;
1049    m_status = STATUS_OKAY;
1050
1051    m_resample_src_x = src_x;
1052    m_resample_src_y = src_y;
1053    m_resample_dst_x = dst_x;
1054    m_resample_dst_y = dst_y;
1055
1056    m_boundary_op = boundary_op;
1057
1058    if ((m_Pdst_buf = (Sample*)malloc(m_resample_dst_x * sizeof(Sample))) == NULL)
1059    {
1060       m_status = STATUS_OUT_OF_MEMORY;
1061       return;
1062    }
1063
1064    // Find the specified filter.
1065    for (i = 0; i < NUM_FILTERS; i++)
1066       if ( filter ==  g_filters[i].name )
1067          break;
1068
1069    if (i == NUM_FILTERS)
1070    {
1071       m_status = STATUS_BAD_FILTER_TYPE;
1072       return;
1073    }
1074
1075    func = g_filters[i].func;
1076    support = g_filters[i].support;
1077
1078    /* Create contributor lists, unless the user supplied custom lists. */
1079
1080    if (!Pclist_x)
1081    {
1082       m_Pclist_x = make_clist(m_resample_src_x, m_resample_dst_x, m_boundary_op, func, support, filter_x_scale, src_x_ofs);
1083       if (!m_Pclist_x)
1084       {
1085          m_status = STATUS_OUT_OF_MEMORY;
1086          return;
1087       }
1088    }
1089    else
1090    {
1091       m_Pclist_x = Pclist_x;
1092       m_clist_x_forced = true;
1093    }
1094
1095    if (!Pclist_y)
1096    {
1097       m_Pclist_y = make_clist(m_resample_src_y, m_resample_dst_y, m_boundary_op, func, support, filter_y_scale, src_y_ofs);
1098       if (!m_Pclist_y)
1099       {
1100          m_status = STATUS_OUT_OF_MEMORY;
1101          return;
1102       }
1103    }
1104    else
1105    {
1106       m_Pclist_y = Pclist_y;
1107       m_clist_y_forced = true;
1108    }
1109
1110    if ((m_Psrc_y_count = (int*)calloc(m_resample_src_y, sizeof(int))) == NULL)
1111    {
1112       m_status = STATUS_OUT_OF_MEMORY;
1113       return;
1114    }
1115
1116    if ((m_Psrc_y_flag = (unsigned char*)calloc(m_resample_src_y, sizeof(unsigned char))) == NULL)
1117    {
1118       m_status = STATUS_OUT_OF_MEMORY;
1119       return;
1120    }
1121
1122    /* Count how many times each source line
1123    * contributes to a destination line.
1124    */
1125
1126    for (i = 0; i < m_resample_dst_y; i++)
1127       for (j = 0; j < m_Pclist_y[i].n; j++)
1128          m_Psrc_y_count[resampler_range_check(m_Pclist_y[i].p[j].pixel, m_resample_src_y)]++;
1129
1130    if ((m_Pscan_buf = (Scan_Buf*)malloc(sizeof(Scan_Buf))) == NULL)
1131    {
1132       m_status = STATUS_OUT_OF_MEMORY;
1133       return;
1134    }
1135
1136    for (i = 0; i < MAX_SCAN_BUF_SIZE; i++)
1137    {
1138       m_Pscan_buf->scan_buf_y[i] = -1;
1139       m_Pscan_buf->scan_buf_l[i] = NULL;
1140    }
1141
1142    m_cur_src_y = m_cur_dst_y = 0;
1143    {
1144       // Determine which axis to resample first by comparing the number of multiplies required
1145       // for each possibility.
1146       int x_ops = count_ops(m_Pclist_x, m_resample_dst_x);
1147       int y_ops = count_ops(m_Pclist_y, m_resample_dst_y);
1148
1149       // Hack 10/2000: Weight Y axis ops a little more than X axis ops.
1150       // (Y axis ops use more cache resources.)
1151       int xy_ops = x_ops * m_resample_src_y +
1152          (4 * y_ops * m_resample_dst_x)/3;
1153
1154       int yx_ops = (4 * y_ops * m_resample_src_x)/3 +
1155          x_ops * m_resample_dst_y;
1156
1157 #if RESAMPLER_DEBUG_OPS
1158       printf("src: %i %i\n", m_resample_src_x, m_resample_src_y);
1159       printf("dst: %i %i\n", m_resample_dst_x, m_resample_dst_y);
1160       printf("x_ops: %i\n", x_ops);
1161       printf("y_ops: %i\n", y_ops);
1162       printf("xy_ops: %i\n", xy_ops);
1163       printf("yx_ops: %i\n", yx_ops);
1164 #endif
1165
1166       // Now check which resample order is better. In case of a tie, choose the order
1167       // which buffers the least amount of data.
1168       if ((xy_ops > yx_ops) ||
1169          ((xy_ops == yx_ops) && (m_resample_src_x < m_resample_dst_x))
1170          )
1171       {
1172          m_delay_x_resample = true;
1173          m_intermediate_x = m_resample_src_x;
1174       }
1175       else
1176       {
1177          m_delay_x_resample = false;
1178          m_intermediate_x = m_resample_dst_x;
1179       }
1180 #if RESAMPLER_DEBUG_OPS
1181       printf("delaying: %i\n", m_delay_x_resample);
1182 #endif
1183    }
1184
1185    if (m_delay_x_resample)
1186    {
1187       if ((m_Ptmp_buf = (Sample*)malloc(m_intermediate_x * sizeof(Sample))) == NULL)
1188       {
1189          m_status = STATUS_OUT_OF_MEMORY;
1190          return;
1191       }
1192    }
1193 }
1194
1195 void Resampler::get_clists(Contrib_List** ptr_clist_x, Contrib_List** ptr_clist_y)
1196 {
1197    if (ptr_clist_x)
1198       *ptr_clist_x = m_Pclist_x;
1199
1200    if (ptr_clist_y)
1201       *ptr_clist_y = m_Pclist_y;
1202 }
1203
1204 #pragma GCC diagnostic pop