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