git.cworth.org Git - vogl/blob - src/voglcore/vogl_resample_filters.cpp

   1 /**************************************************************************
   2  *
   3  * Copyright 2013-2014 RAD Game Tools and Valve Software
   4  * Copyright 2010-2014 Rich Geldreich and Tenacious Software LLC
   5  * All Rights Reserved.
   6  *
   7  * Permission is hereby granted, free of charge, to any person obtaining a copy
   8  * of this software and associated documentation files (the "Software"), to deal
   9  * in the Software without restriction, including without limitation the rights
  10  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  11  * copies of the Software, and to permit persons to whom the Software is
  12  * furnished to do so, subject to the following conditions:
  13  *
  14  * The above copyright notice and this permission notice shall be included in
  15  * all copies or substantial portions of the Software.
  16  *
  17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  18  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  20  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  21  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  22  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  23  * THE SOFTWARE.
  24  *
  25  **************************************************************************/
  26
  27 // File: vogl_resample_filters.cpp
  28 // RG: This is public domain code, originally derived from Graphics Gems 3, see: http://code.google.com/p/imageresampler/
  29 #include "vogl_core.h"
  30 #include "vogl_resample_filters.h"
  31 #include "vogl_strutils.h"
  32
  33 namespace vogl
  34 {
  35 #define M_PI 3.14159265358979323846
  36
  37 // To add your own filter, insert the new function below and update the filter table.
  38 // There is no need to make the filter function particularly fast, because it's
  39 // only called during initializing to create the X and Y axis contributor tables.
  40
  41 #define BOX_FILTER_SUPPORT (0.5f)
  42     static float box_filter(float t) /* pulse/Fourier window */
  43     {
  44         // make_clist() calls the filter function with t inverted (pos = left, neg = right)
  45         if ((t >= -0.5f) && (t < 0.5f))
  46             return 1.0f;
  47         else
  48             return 0.0f;
  49     }
  50
  51 #define TENT_FILTER_SUPPORT (1.0f)
  52     static float tent_filter(float t) /* box (*) box, bilinear/triangle */
  53     {
  54         if (t < 0.0f)
  55             t = -t;
  56
  57         if (t < 1.0f)
  58             return 1.0f - t;
  59         else
  60             return 0.0f;
  61     }
  62
  63 #define BELL_SUPPORT (1.5f)
  64     static float bell_filter(float t) /* box (*) box (*) box */
  65     {
  66         if (t < 0.0f)
  67             t = -t;
  68
  69         if (t < .5f)
  70             return (.75f - (t * t));
  71
  72         if (t < 1.5f)
  73         {
  74             t = (t - 1.5f);
  75             return (.5f * (t * t));
  76         }
  77
  78         return (0.0f);
  79     }
  80
  81 #define B_SPLINE_SUPPORT (2.0f)
  82     static float B_spline_filter(float t) /* box (*) box (*) box (*) box */
  83     {
  84         float tt;
  85
  86         if (t < 0.0f)
  87             t = -t;
  88
  89         if (t < 1.0f)
  90         {
  91             tt = t * t;
  92             return ((.5f * tt * t) - tt + (2.0f / 3.0f));
  93         }
  94         else if (t < 2.0f)
  95         {
  96             t = 2.0f - t;
  97             return ((1.0f / 6.0f) * (t * t * t));
  98         }
  99
 100         return (0.0f);
 101     }
 102
 103 // Dodgson, N., "Quadratic Interpolation for Image Resampling"
 104 #define QUADRATIC_SUPPORT 1.5f
 105     static float quadratic(float t, const float R)
 106     {
 107         if (t < 0.0f)
 108             t = -t;
 109         if (t < QUADRATIC_SUPPORT)
 110         {
 111             float tt = t * t;
 112             if (t <= .5f)
 113                 return (-2.0f * R) * tt + .5f * (R + 1.0f);
 114             else
 115                 return (R * tt) + (-2.0f * R - .5f) * t + (3.0f / 4.0f) * (R + 1.0f);
 116         }
 117         else
 118             return 0.0f;
 119     }
 120
 121     static float quadratic_interp_filter(float t)
 122     {
 123         return quadratic(t, 1.0f);
 124     }
 125
 126     static float quadratic_approx_filter(float t)
 127     {
 128         return quadratic(t, .5f);
 129     }
 130
 131     static float quadratic_mix_filter(float t)
 132     {
 133         return quadratic(t, .8f);
 134     }
 135
 136     // Mitchell, D. and A. Netravali, "Reconstruction Filters in Computer Graphics."
 137     // Computer Graphics, Vol. 22, No. 4, pp. 221-228.
 138     // (B, C)
 139     // (1/3, 1/3)  - Defaults recommended by Mitchell and Netravali
 140     // (1, 0)      - Equivalent to the Cubic B-Spline
 141     // (0, 0.5)         - Equivalent to the Catmull-Rom Spline
 142     // (0, C)           - The family of Cardinal Cubic Splines
 143     // (B, 0)           - Duff's tensioned B-Splines.
 144     static float mitchell(float t, const float B, const float C)
 145     {
 146         float tt;
 147
 148         tt = t * t;
 149
 150         if (t < 0.0f)
 151             t = -t;
 152
 153         if (t < 1.0f)
 154         {
 155             t = (((12.0f - 9.0f * B - 6.0f * C) * (t * tt)) + ((-18.0f + 12.0f * B + 6.0f * C) * tt) + (6.0f - 2.0f * B));
 156
 157             return (t / 6.0f);
 158         }
 159         else if (t < 2.0f)
 160         {
 161             t = (((-1.0f * B - 6.0f * C) * (t * tt)) + ((6.0f * B + 30.0f * C) * tt) + ((-12.0f * B - 48.0f * C) * t) + (8.0f * B + 24.0f * C));
 162
 163             return (t / 6.0f);
 164         }
 165
 166         return (0.0f);
 167     }
 168
 169 #define MITCHELL_SUPPORT (2.0f)
 170     static float mitchell_filter(float t)
 171     {
 172         return mitchell(t, 1.0f / 3.0f, 1.0f / 3.0f);
 173     }
 174
 175 #define CATMULL_ROM_SUPPORT (2.0f)
 176     static float catmull_rom_filter(float t)
 177     {
 178         return mitchell(t, 0.0f, .5f);
 179     }
 180
 181     static double sinc(double x)
 182     {
 183         x = (x * M_PI);
 184
 185         if ((x < 0.01f) && (x > -0.01f))
 186             return 1.0f + x * x * (-1.0f / 6.0f + x * x * 1.0f / 120.0f);
 187
 188         return sin(x) / x;
 189     }
 190
 191     static float clean(double t)
 192     {
 193         const float EPSILON = .0000125f;
 194         if (fabs(t) < EPSILON)
 195             return 0.0f;
 196         return (float)t;
 197     }
 198
 199     //static double blackman_window(double x)
 200     //{
 201     //  return .42f + .50f * cos(M_PI*x) + .08f * cos(2.0f*M_PI*x);
 202     //}
 203
 204     static double blackman_exact_window(double x)
 205     {
 206         return 0.42659071f + 0.49656062f * cos(M_PI * x) + 0.07684867f * cos(2.0f * M_PI * x);
 207     }
 208
 209 #define BLACKMAN_SUPPORT (3.0f)
 210     static float blackman_filter(float t)
 211     {
 212         if (t < 0.0f)
 213             t = -t;
 214
 215         if (t < 3.0f)
 216             //return clean(sinc(t) * blackman_window(t / 3.0f));
 217             return clean(sinc(t) * blackman_exact_window(t / 3.0f));
 218         else
 219             return (0.0f);
 220     }
 221
 222 #define GAUSSIAN_SUPPORT (1.25f)
 223     static float gaussian_filter(float t) // with blackman window
 224     {
 225         if (t < 0)
 226             t = -t;
 227         if (t < GAUSSIAN_SUPPORT)
 228             return clean(exp(-2.0f * t * t) * sqrt(2.0f / M_PI) * blackman_exact_window(t / GAUSSIAN_SUPPORT));
 229         else
 230             return 0.0f;
 231     }
 232
 233 // Windowed sinc -- see "Jimm Blinn's Corner: Dirty Pixels" pg. 26.
 234 #define LANCZOS3_SUPPORT (3.0f)
 235     static float lanczos3_filter(float t)
 236     {
 237         if (t < 0.0f)
 238             t = -t;
 239
 240         if (t < 3.0f)
 241             return clean(sinc(t) * sinc(t / 3.0f));
 242         else
 243             return (0.0f);
 244     }
 245
 246 #define LANCZOS4_SUPPORT (4.0f)
 247     static float lanczos4_filter(float t)
 248     {
 249         if (t < 0.0f)
 250             t = -t;
 251
 252         if (t < 4.0f)
 253             return clean(sinc(t) * sinc(t / 4.0f));
 254         else
 255             return (0.0f);
 256     }
 257
 258 #define LANCZOS6_SUPPORT (6.0f)
 259     static float lanczos6_filter(float t)
 260     {
 261         if (t < 0.0f)
 262             t = -t;
 263
 264         if (t < 6.0f)
 265             return clean(sinc(t) * sinc(t / 6.0f));
 266         else
 267             return (0.0f);
 268     }
 269
 270 #define LANCZOS12_SUPPORT (12.0f)
 271     static float lanczos12_filter(float t)
 272     {
 273         if (t < 0.0f)
 274             t = -t;
 275
 276         if (t < 12.0f)
 277             return clean(sinc(t) * sinc(t / 12.0f));
 278         else
 279             return (0.0f);
 280     }
 281
 282     static double bessel0(double x)
 283     {
 284         const double EPSILON_RATIO = 1E-16;
 285         double xh, sum, pow, ds;
 286         int k;
 287
 288         xh = 0.5 * x;
 289         sum = 1.0;
 290         pow = 1.0;
 291         k = 0;
 292         ds = 1.0;
 293         while (ds > sum * EPSILON_RATIO) // FIXME: Shouldn't this stop after X iterations for max. safety?
 294         {
 295             ++k;
 296             pow = pow * (xh / k);
 297             ds = pow * pow;
 298             sum = sum + ds;
 299         }
 300
 301         return sum;
 302     }
 303
 304     // static const float KAISER_ALPHA = 4.0;
 305     static double kaiser(double alpha, double half_width, double x)
 306     {
 307         const double ratio = (x / half_width);
 308         return bessel0(alpha * sqrt(1 - ratio * ratio)) / bessel0(alpha);
 309     }
 310
 311 #define KAISER_SUPPORT 3
 312     static float kaiser_filter(float t)
 313     {
 314         if (t < 0.0f)
 315             t = -t;
 316
 317         if (t < KAISER_SUPPORT)
 318         {
 319             // db atten
 320             const float att = 40.0f;
 321             const float alpha = (float)(exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886 * (att - 20.96));
 322             //const float alpha = KAISER_ALPHA;
 323             return (float)clean(sinc(t) * kaiser(alpha, KAISER_SUPPORT, t));
 324         }
 325
 326         return 0.0f;
 327     }
 328
 329     const resample_filter g_resample_filters[] =
 330         {
 331             { "box", box_filter, BOX_FILTER_SUPPORT },
 332             { "tent", tent_filter, TENT_FILTER_SUPPORT },
 333             { "bell", bell_filter, BELL_SUPPORT },
 334             { "b-spline", B_spline_filter, B_SPLINE_SUPPORT },
 335             { "mitchell", mitchell_filter, MITCHELL_SUPPORT },
 336             { "lanczos3", lanczos3_filter, LANCZOS3_SUPPORT },
 337             { "blackman", blackman_filter, BLACKMAN_SUPPORT },
 338             { "lanczos4", lanczos4_filter, LANCZOS4_SUPPORT },
 339             { "lanczos6", lanczos6_filter, LANCZOS6_SUPPORT },
 340             { "lanczos12", lanczos12_filter, LANCZOS12_SUPPORT },
 341             { "kaiser", kaiser_filter, KAISER_SUPPORT },
 342             { "gaussian", gaussian_filter, GAUSSIAN_SUPPORT },
 343             { "catmullrom", catmull_rom_filter, CATMULL_ROM_SUPPORT },
 344             { "quadratic_interp", quadratic_interp_filter, QUADRATIC_SUPPORT },
 345             { "quadratic_approx", quadratic_approx_filter, QUADRATIC_SUPPORT },
 346             { "quadratic_mix", quadratic_mix_filter, QUADRATIC_SUPPORT },
 347         };
 348
 349     const int g_num_resample_filters = sizeof(g_resample_filters) / sizeof(g_resample_filters[0]);
 350
 351     int find_resample_filter(const char *pName)
 352     {
 353         for (int i = 0; i < g_num_resample_filters; i++)
 354             if (vogl::vogl_stricmp(pName, g_resample_filters[i].name) == 0)
 355                 return i;
 356         return cInvalidIndex;
 357     }
 358
 359 } // namespace vogl