1 /**************************************************************************
3 * Copyright 2013-2014 RAD Game Tools and Valve Software
4 * Copyright 2010-2014 Rich Geldreich and Tenacious Software LLC
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:
14 * The above copyright notice and this permission notice shall be included in
15 * all copies or substantial portions of the Software.
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
25 **************************************************************************/
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"
35 #define M_PI 3.14159265358979323846
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.
41 #define BOX_FILTER_SUPPORT (0.5f)
42 static float box_filter(float t) /* pulse/Fourier window */
44 // make_clist() calls the filter function with t inverted (pos = left, neg = right)
45 if ((t >= -0.5f) && (t < 0.5f))
51 #define TENT_FILTER_SUPPORT (1.0f)
52 static float tent_filter(float t) /* box (*) box, bilinear/triangle */
63 #define BELL_SUPPORT (1.5f)
64 static float bell_filter(float t) /* box (*) box (*) box */
70 return (.75f - (t * t));
75 return (.5f * (t * t));
81 #define B_SPLINE_SUPPORT (2.0f)
82 static float B_spline_filter(float t) /* box (*) box (*) box (*) box */
92 return ((.5f * tt * t) - tt + (2.0f / 3.0f));
97 return ((1.0f / 6.0f) * (t * t * t));
103 // Dodgson, N., "Quadratic Interpolation for Image Resampling"
104 #define QUADRATIC_SUPPORT 1.5f
105 static float quadratic(float t, const float R)
109 if (t < QUADRATIC_SUPPORT)
113 return (-2.0f * R) * tt + .5f * (R + 1.0f);
115 return (R * tt) + (-2.0f * R - .5f) * t + (3.0f / 4.0f) * (R + 1.0f);
121 static float quadratic_interp_filter(float t)
123 return quadratic(t, 1.0f);
126 static float quadratic_approx_filter(float t)
128 return quadratic(t, .5f);
131 static float quadratic_mix_filter(float t)
133 return quadratic(t, .8f);
136 // Mitchell, D. and A. Netravali, "Reconstruction Filters in Computer Graphics."
137 // Computer Graphics, Vol. 22, No. 4, pp. 221-228.
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)
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));
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));
169 #define MITCHELL_SUPPORT (2.0f)
170 static float mitchell_filter(float t)
172 return mitchell(t, 1.0f / 3.0f, 1.0f / 3.0f);
175 #define CATMULL_ROM_SUPPORT (2.0f)
176 static float catmull_rom_filter(float t)
178 return mitchell(t, 0.0f, .5f);
181 static double sinc(double x)
185 if ((x < 0.01f) && (x > -0.01f))
186 return 1.0f + x * x * (-1.0f / 6.0f + x * x * 1.0f / 120.0f);
191 static float clean(double t)
193 const float EPSILON = .0000125f;
194 if (fabs(t) < EPSILON)
199 //static double blackman_window(double x)
201 // return .42f + .50f * cos(M_PI*x) + .08f * cos(2.0f*M_PI*x);
204 static double blackman_exact_window(double x)
206 return 0.42659071f + 0.49656062f * cos(M_PI * x) + 0.07684867f * cos(2.0f * M_PI * x);
209 #define BLACKMAN_SUPPORT (3.0f)
210 static float blackman_filter(float t)
216 //return clean(sinc(t) * blackman_window(t / 3.0f));
217 return clean(sinc(t) * blackman_exact_window(t / 3.0f));
222 #define GAUSSIAN_SUPPORT (1.25f)
223 static float gaussian_filter(float t) // with blackman window
227 if (t < GAUSSIAN_SUPPORT)
228 return clean(exp(-2.0f * t * t) * sqrt(2.0f / M_PI) * blackman_exact_window(t / GAUSSIAN_SUPPORT));
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)
241 return clean(sinc(t) * sinc(t / 3.0f));
246 #define LANCZOS4_SUPPORT (4.0f)
247 static float lanczos4_filter(float t)
253 return clean(sinc(t) * sinc(t / 4.0f));
258 #define LANCZOS6_SUPPORT (6.0f)
259 static float lanczos6_filter(float t)
265 return clean(sinc(t) * sinc(t / 6.0f));
270 #define LANCZOS12_SUPPORT (12.0f)
271 static float lanczos12_filter(float t)
277 return clean(sinc(t) * sinc(t / 12.0f));
282 static double bessel0(double x)
284 const double EPSILON_RATIO = 1E-16;
285 double xh, sum, pow, ds;
293 while (ds > sum * EPSILON_RATIO) // FIXME: Shouldn't this stop after X iterations for max. safety?
296 pow = pow * (xh / k);
304 // static const float KAISER_ALPHA = 4.0;
305 static double kaiser(double alpha, double half_width, double x)
307 const double ratio = (x / half_width);
308 return bessel0(alpha * sqrt(1 - ratio * ratio)) / bessel0(alpha);
311 #define KAISER_SUPPORT 3
312 static float kaiser_filter(float t)
317 if (t < KAISER_SUPPORT)
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));
329 const resample_filter g_resample_filters[] =
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 },
349 const int g_num_resample_filters = sizeof(g_resample_filters) / sizeof(g_resample_filters[0]);
351 int find_resample_filter(const char *pName)
353 for (int i = 0; i < g_num_resample_filters; i++)
354 if (vogl::vogl_stricmp(pName, g_resample_filters[i].name) == 0)
356 return cInvalidIndex;