diff --git a/src/nvimage/Filter.cpp b/src/nvimage/Filter.cpp index 583e47f..37e35e5 100644 --- a/src/nvimage/Filter.cpp +++ b/src/nvimage/Filter.cpp @@ -9,6 +9,7 @@ * * References from Thacher Ulrich: * See _Graphics Gems III_ "General Filtered Image Rescaling", Dale A. Schumacher + * http://tog.acm.org/GraphicsGems/gemsiii/filter.c * * References from Paul Heckbert: * A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975 @@ -27,6 +28,9 @@ * Reconstruction Filters in Computer Graphics * http://www.mentallandscape.com/Papers_siggraph88.pdf * + * More references: + * http://www.worldserver.com/turk/computergraphics/ResamplingFilters.pdf + * http://www.dspguide.com/ch16.htm */ @@ -39,28 +43,107 @@ using namespace nv; namespace { + // Sinc function. + inline static float sincf(const float x) + { + if (fabs(x) < NV_EPSILON) { + //return 1.0; + return 1.0f + x*x*(-1.0f/6.0f + x*x*1.0f/120.0f); + } + else { + return sin(x) / x; + } + } + + // Bessel function of the first kind from Jon Blow's article. + // http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html + // http://en.wikipedia.org/wiki/Bessel_function + inline static float bessel0(float x) + { + const float EPSILON_RATIO = 1E-6; + float xh, sum, pow, ds; + int k; + + xh = 0.5 * x; + sum = 1.0; + pow = 1.0; + k = 0; + ds = 1.0; + while (ds > sum * EPSILON_RATIO) { + ++k; + pow = pow * (xh / k); + ds = pow * pow; + sum = sum + ds; + } + + return sum; + } + + /*// Alternative bessel function from Paul Heckbert. + static float _bessel0(float x) + { + const float EPSILON_RATIO = 1E-6; + float sum = 1.0f; + float y = x * x / 4.0f; + float t = y; + for(int i = 2; t > EPSILON_RATIO; i++) { + sum += t; + t *= y / float(i * i); + } + return sum; + }*/ -// support = 0.5 -inline static float filter_box(float x) +} // namespace + + +Filter::Filter(float width) : m_width(width) { - if( x < -0.5f ) return 0.0f; - if( x <= 0.5 ) return 1.0f; - return 0.0f; } -// support = 1.0 -inline static float filter_triangle(float x) +float Filter::sample(float x, float scale, int samples) const { - if( x < -1.0f ) return 0.0f; - if( x < 0.0f ) return 1.0f + x; - if( x < 1.0f ) return 1.0f - x; +// return evaluate(x * scale); + + float sum = 0; + float isamples = 1.0f / float(samples); + + for(int s = 0; s < samples; s++) + { + float p = (x + (float(s) + 0.5f) * isamples) * scale; + float value = evaluate(p); + sum += value; + } + + return sum * isamples; +} + + +BoxFilter::BoxFilter() : Filter(0.5f) {} +BoxFilter::BoxFilter(float width) : Filter(width) {} + +float BoxFilter::evaluate(float x) const +{ + if (fabs(x) <= m_width) return 1.0f; + else return 0.0f; +} + + +TriangleFilter::TriangleFilter() : Filter(1.0f) {} +TriangleFilter::TriangleFilter(float width) : Filter(width) {} + +float TriangleFilter::evaluate(float x) const +{ + x = fabs(x); + if( x < m_width ) return m_width - x; return 0.0f; } -// support = 1.5 -inline static float filter_quadratic(float x) + +QuadraticFilter::QuadraticFilter() : Filter(1.5f) {} + +float QuadraticFilter::evaluate(float x) const { - if( x < 0.0f ) x = -x; + x = fabs(x); if( x < 0.5f ) return 0.75f - x * x; if( x < 1.5f ) { float t = x - 1.5f; @@ -69,22 +152,23 @@ inline static float filter_quadratic(float x) return 0.0f; } -// @@ Filter from tulrich. -// support 1.0 -inline static float filter_cubic(float x) + +CubicFilter::CubicFilter() : Filter(1.0f) {} + +float CubicFilter::evaluate(float x) const { // f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1 - if( x < 0.0f ) x = -x; + x = fabs(x); if( x < 1.0f ) return((2.0f * x - 3.0f) * x * x + 1.0f); return 0.0f; } -// @@ Paul Heckbert calls this cubic instead of spline. -// support = 2.0 -inline static float filter_spline(float x) +BSplineFilter::BSplineFilter() : Filter(2.0f) {} + +float BSplineFilter::evaluate(float x) const { - if( x < 0.0f ) x = -x; + x = fabs(x); if( x < 1.0f ) return (4.0f + x * x * (-6.0f + x * 3.0f)) / 6.0f; if( x < 2.0f ) { float t = 2.0f - x; @@ -93,132 +177,64 @@ inline static float filter_spline(float x) return 0.0f; } -/// Sinc function. -inline float sincf( const float x ) -{ - if( fabs(x) < NV_EPSILON ) { - return 1.0 ; - //return 1.0f + x*x*(-1.0f/6.0f + x*x*1.0f/120.0f); - } - else { - return sin(x) / x; - } -} - -// support = 3.0 -inline static float filter_lanczos3(float x) -{ - if( x < 0.0f ) x = -x; - if( x < 3.0f ) return sincf(x) * sincf(x / 3.0f); - return 0.0f; -} - +MitchellFilter::MitchellFilter() : Filter(2.0f) { setParameters(1.0f/3.0f, 1.0f/3.0f); } -// Mitchell & Netravali's two-param cubic -// see "Reconstruction Filters in Computer Graphics", SIGGRAPH 88 -// support = 2.0 -inline static float filter_mitchell(float x, float b, float c) +float MitchellFilter::evaluate(float x) const { - // @@ Coefficients could be precomputed. - // @@ if b and c are fixed, these are constants. - const float p0 = (6.0f - 2.0f * b) / 6.0f; - const float p2 = (-18.0f + 12.0f * b + 6.0f * c) / 6.0f; - const float p3 = (12.0f - 9.0f * b - 6.0f * c) / 6.0f; - const float q0 = (8.0f * b + 24.0f * c) / 6.0f; - const float q1 = (-12.0f * b - 48.0f * c) / 6.0f; - const float q2 = (6.0f * b + 30.0f * c) / 6.0f; - const float q3 = (-b - 6.0f * c) / 6.0f; - - if( x < 0.0f ) x = -x; + x = fabs(x); if( x < 1.0f ) return p0 + x * x * (p2 + x * p3); if( x < 2.0f ) return q0 + x * (q1 + x * (q2 + x * q3)); return 0.0f; } -inline static float filter_mitchell(float x) +void MitchellFilter::setParameters(float b, float c) { - return filter_mitchell(x, 1.0f/3.0f, 1.0f/3.0f); + p0 = (6.0f - 2.0f * b) / 6.0f; + p2 = (-18.0f + 12.0f * b + 6.0f * c) / 6.0f; + p3 = (12.0f - 9.0f * b - 6.0f * c) / 6.0f; + q0 = (8.0f * b + 24.0f * c) / 6.0f; + q1 = (-12.0f * b - 48.0f * c) / 6.0f; + q2 = (6.0f * b + 30.0f * c) / 6.0f; + q3 = (-b - 6.0f * c) / 6.0f; } -// Bessel function of the first kind from Jon Blow's article. -// http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html -// http://en.wikipedia.org/wiki/Bessel_function -static float bessel0(float x) -{ - const float EPSILON_RATIO = 1E-6; - float xh, sum, pow, ds; - int k; - xh = 0.5 * x; - sum = 1.0; - pow = 1.0; - k = 0; - ds = 1.0; - while (ds > sum * EPSILON_RATIO) { - ++k; - pow = pow * (xh / k); - ds = pow * pow; - sum = sum + ds; - } +LanczosFilter::LanczosFilter() : Filter(3.0f) {} - return sum; +float LanczosFilter::evaluate(float x) const +{ + x = fabs(x); + if( x < 3.0f ) return sincf(PI * x) * sincf(PI * x / 3.0f); + return 0.0f; } -/*// Alternative bessel function from Paul Heckbert. -static float _bessel0(float x) -{ - const float EPSILON_RATIO = 1E-6; - float sum = 1.0f; - float y = x * x / 4.0f; - float t = y; - for(int i = 2; t > EPSILON_RATIO; i++) { - sum += t; - t *= y / float(i * i); - } - return sum; -}*/ -// support = 1.0 -inline static float filter_kaiser(float x, float alpha) -{ - return bessel0(alpha * sqrtf(1 - x * x)) / bessel0(alpha); -} +SincFilter::SincFilter(float w) : Filter(w) {} -inline static float filter_kaiser(float x) +float SincFilter::evaluate(float x) const { - return filter_kaiser(x, 4.0f); + return 0.0f; } -// Array of filters. -static Filter s_filter_array[] = { - {filter_box, 0.5f}, // Box - {filter_triangle, 1.0f}, // Triangle - {filter_quadratic, 1.5f}, // Quadratic - {filter_cubic, 1.0f}, // Cubic - {filter_spline, 2.0f}, // Spline - {filter_lanczos3, 3.0f}, // Lanczos - {filter_mitchell, 1.0f}, // Mitchell - {filter_kaiser, 1.0f}, // Kaiser -}; +KaiserFilter::KaiserFilter(float w) : Filter(w) { setParameters(4.0f, 1.0f); } - -inline static float sampleFilter(Filter::Function func, float x, float scale, int samples) +float KaiserFilter::evaluate(float x) const { - float sum = 0; - - for(int s = 0; s < samples; s++) - { - sum += func((x + (float(s) + 0.5f) * (1.0f / float(samples))) * scale); - } - - return sum; + const float sinc_value = sincf(PI * x * stretch); + float t = x / m_width; + if (t * t <= 1.0f) + return sinc_value * bessel0(alpha * sqrtf(1 - t * t)) / bessel0(alpha); + else + return 0; } - - -} // namespace +void KaiserFilter::setParameters(float alpha, float stretch) +{ + this->alpha = alpha; + this->stretch = stretch; +} @@ -257,7 +273,7 @@ void Kernel1::normalize() } } - +#if 0 /// Init 1D filter. void Kernel1::initFilter(Filter::Enum f, int samples /*= 1*/) { @@ -334,7 +350,7 @@ void Kernel1::initMitchell(float b, float c) normalize(); } - +#endif /// Print the kernel for debugging purposes. void Kernel1::debugPrint() @@ -585,7 +601,7 @@ static bool isMonoPhase(float w) } - +/* PolyphaseKernel::PolyphaseKernel(float w, uint l) : m_width(w), m_size(ceilf(w) + 1), @@ -603,98 +619,60 @@ PolyphaseKernel::PolyphaseKernel(const PolyphaseKernel & k) : m_data = new float[m_size * m_length]; memcpy(m_data, k.m_data, sizeof(float) * m_size * m_length); } +*/ -PolyphaseKernel::~PolyphaseKernel() +PolyphaseKernel::PolyphaseKernel(const Filter & f, uint srcLength, uint dstLength) { - delete [] m_data; -} + float scale = float(dstLength) / float(srcLength); + float iscale = 1.0f / scale; + m_length = dstLength; + m_width = f.width() * iscale; + m_windowSize = ceilf(m_width * 2) + 1; -/* @@ Should we precompute left & right? - - // scale factor - double dScale = double(uDstSize) / double(uSrcSize); - - if(dScale < 1.0) { - // minification - dWidth = dFilterWidth / dScale; - dFScale = dScale; - } else { - // magnification - dWidth= dFilterWidth; - } - - // window size is the number of sampled pixels - m_WindowSize = 2 * (int)ceil(dWidth) + 1; - m_LineLength = uDstSize; - - // offset for discrete to continuous coordinate conversion - double dOffset = (0.5 / dScale) - 0.5; - - for(u = 0; u < m_LineLength; u++) { - // scan through line of contributions - double dCenter = (double)u / dScale + dOffset; // reverse mapping - // find the significant edge points that affect the pixel - int iLeft = MAX (0, (int)floor (dCenter - dWidth)); - int iRight = MIN ((int)ceil (dCenter + dWidth), int(uSrcSize) - 1); - - ... - } -*/ + m_data = new float[m_windowSize * m_length]; + memset(m_data, 0, sizeof(float) * m_windowSize * m_length); -void PolyphaseKernel::initFilter(Filter::Enum f, int samples/*= 1*/) -{ - nvCheck(f < Filter::Num); - - float (* filter_function)(float) = s_filter_array[f].function; - const float support = s_filter_array[f].support; - - const float half_width = m_width / 2; - const float scale = support / half_width; - - for (uint j = 0; j < m_length; j++) + for (uint i = 0; i < m_length; i++) { - const float phase = frac(m_width * j); - const float offset = half_width + phase; + const float center = (0.5f + i) * iscale; - nvDebug("%d: ", j); + int left = floor(center - m_width); + int right = ceil(center + m_width); + nvCheck(right - left <= (int)m_windowSize); float total = 0.0f; - for (uint i = 0; i < m_size; i++) + for (int j = 0; j < m_windowSize; j++) { - float sample = sampleFilter(filter_function, i - offset, scale, samples); - - nvDebug("(%5.3f | %d) ", sample, j + i - m_size/2); + float sample = f.sample(left + j - center, scale, 40); - m_data[j * m_size + i] = sample; + m_data[i * m_windowSize + j] = sample; total += sample; } - nvDebug("\n"); - // normalize weights. - for (uint i = 0; i < m_size; i++) + for (int j = 0; j < m_windowSize; j++) { - m_data[j * m_size + i] /= total; - //m_data[j * m_size + i] /= samples; + m_data[i * m_windowSize + j] /= total; } } } -void PolyphaseKernel::initKaiser(float alpha /*= 4.0f*/, float stretch /*= 1.0f*/) +PolyphaseKernel::~PolyphaseKernel() { - + delete [] m_data; } + /// Print the kernel for debugging purposes. -void PolyphaseKernel::debugPrint() +void PolyphaseKernel::debugPrint() const { - for (uint j = 0; j < m_length; j++) + for (uint i = 0; i < m_length; i++) { - nvDebug("%d: ", j); - for (uint i = 0; i < m_size; i++) + nvDebug("%d: ", i); + for (uint j = 0; j < m_windowSize; j++) { - nvDebug(" %6.4f", m_data[j * m_size + i]); + nvDebug(" %6.4f", m_data[i * m_windowSize + j]); } nvDebug("\n"); } diff --git a/src/nvimage/Filter.h b/src/nvimage/Filter.h index 4e160ae..7ce1a7e 100644 --- a/src/nvimage/Filter.h +++ b/src/nvimage/Filter.h @@ -9,30 +9,110 @@ namespace nv { class Vector4; - /// A filter function. - struct Filter + /// Base filter class. + class Filter { - // Standard filters. - enum Enum - { - Box, - Triangle, - Quadratic, // Bell - Cubic, - Spline, - Lanczos, - Mitchell, - Kaiser, // Kaiser-windowed sinc filter - Num - }; - - typedef float (* Function)(float); + public: + NVIMAGE_API Filter(float width); + + NVIMAGE_API float width() const { return m_width; } + NVIMAGE_API float sample(float x, float scale, int samples) const; + + virtual float evaluate(float x) const = 0; + + protected: + const float m_width; + }; + + // Box filter. + class BoxFilter : public Filter + { + public: + NVIMAGE_API BoxFilter(); + NVIMAGE_API BoxFilter(float width); + NVIMAGE_API virtual float evaluate(float x) const; + }; + + // Triangle (bilinear/tent) filter. + class TriangleFilter : public Filter + { + public: + NVIMAGE_API TriangleFilter(); + NVIMAGE_API TriangleFilter(float width); + NVIMAGE_API virtual float evaluate(float x) const; + }; + + // Quadratic (bell) filter. + class QuadraticFilter : public Filter + { + public: + NVIMAGE_API QuadraticFilter(); + NVIMAGE_API virtual float evaluate(float x) const; + }; + + // Cubic filter from Thatcher Ulrich. + class CubicFilter : public Filter + { + public: + NVIMAGE_API CubicFilter(); + NVIMAGE_API virtual float evaluate(float x) const; + }; + + // Cubic b-spline filter from Paul Heckbert. + class BSplineFilter : public Filter + { + public: + NVIMAGE_API BSplineFilter(); + NVIMAGE_API virtual float evaluate(float x) const; + }; + + /// Mitchell & Netravali's two-param cubic + /// @see "Reconstruction Filters in Computer Graphics", SIGGRAPH 88 + class MitchellFilter : public Filter + { + public: + NVIMAGE_API MitchellFilter(); + NVIMAGE_API virtual float evaluate(float x) const; + + NVIMAGE_API void setParameters(float a, float b); + + private: + float p0, p2, p3; + float q0, q1, q2, q3; + }; + + // Lanczos3 filter. + class LanczosFilter : public Filter + { + public: + NVIMAGE_API LanczosFilter(); + NVIMAGE_API virtual float evaluate(float x) const; + }; + + // Sinc filter. + class SincFilter : public Filter + { + public: + NVIMAGE_API SincFilter(float w); + NVIMAGE_API virtual float evaluate(float x) const; + }; + + // Kaiser filter. + class KaiserFilter : public Filter + { + public: + NVIMAGE_API KaiserFilter(float w); + NVIMAGE_API virtual float evaluate(float x) const; + + NVIMAGE_API void setParameters(float a, float stretch); - Function function; - float support; + private: + float alpha; + float stretch; }; + /// A 1D kernel. Used to precompute filter weights. class Kernel1 { @@ -50,11 +130,12 @@ namespace nv uint windowSize() const { return m_windowSize; } - + /* NVIMAGE_API void initFilter(Filter::Enum filter, int samples = 1); NVIMAGE_API void initSinc(float stretch = 1); NVIMAGE_API void initKaiser(float alpha = 4.0f, float stretch = 1.0f, int sampes = 1); NVIMAGE_API void initMitchell(float b = 1.0f/3.0f, float c = 1.0f/3.0f); + */ NVIMAGE_API void debugPrint(); @@ -95,39 +176,39 @@ namespace nv float * m_data; }; + /// A 1D polyphase kernel class PolyphaseKernel { + NV_FORBID_COPY(PolyphaseKernel) public: - NVIMAGE_API PolyphaseKernel(float width, uint lineLength); - NVIMAGE_API PolyphaseKernel(const PolyphaseKernel & k); + NVIMAGE_API PolyphaseKernel(const Filter & f, uint srcLength, uint dstLength); NVIMAGE_API ~PolyphaseKernel(); - float valueAt(uint column, uint x) const { - return m_data[column * m_size + x]; - } - - float width() const { - return m_width; - } - uint windowSize() const { - return m_size; + return m_windowSize; } uint length() const { return m_length; } - - NVIMAGE_API void initFilter(Filter::Enum filter, int samples = 1); - NVIMAGE_API void initKaiser(float alpha = 4.0f, float stretch = 0.5f); - - NVIMAGE_API void debugPrint(); + + float width() const { + return m_width; + } + + float valueAt(uint column, uint x) const { + nvDebugCheck(column < m_length); + nvDebugCheck(x < m_windowSize); + return m_data[column * m_windowSize + x]; + } + + NVIMAGE_API void debugPrint() const; private: - const float m_width; - const uint m_size; - const uint m_length; + uint m_windowSize; + uint m_length; + float m_width; float * m_data; }; diff --git a/src/nvimage/FloatImage.cpp b/src/nvimage/FloatImage.cpp index d720c8d..5dd812f 100644 --- a/src/nvimage/FloatImage.cpp +++ b/src/nvimage/FloatImage.cpp @@ -535,7 +535,7 @@ FloatImage * FloatImage::fastDownSample() const return dst_image.release(); } - +/* /// Downsample applying a 1D kernel separately in each dimension. FloatImage * FloatImage::downSample(const Kernel1 & kernel, WrapMode wm) const { @@ -588,59 +588,83 @@ FloatImage * FloatImage::downSample(const Kernel1 & kernel, uint w, uint h, Wrap return dst_image.release(); } +*/ +/// Downsample applying a 1D kernel separately in each dimension. +FloatImage * FloatImage::downSample(const Filter & filter, WrapMode wm) const +{ + const uint w = max(1, m_width / 2); + const uint h = max(1, m_height / 2); + + return downSample(filter, w, h, wm); +} /// Downsample applying a 1D kernel separately in each dimension. -FloatImage * FloatImage::downSample(uint w, uint h, WrapMode wm) const +FloatImage * FloatImage::downSample(const Filter & filter, uint w, uint h, WrapMode wm) const { - // Build polyphase kernels. - const float xscale = float(m_width) / float(w); - const float yscale = float(m_height) / float(h); - - int kw = 1; - float xwidth = kw * xscale; - float ywidth = kw * yscale; - - PolyphaseKernel xkernel(xwidth, w); - PolyphaseKernel ykernel(ywidth, h); - - xkernel.initFilter(Filter::Box, 32); - ykernel.initFilter(Filter::Box, 32); -// xkernel.initKaiser(4, 1.0f / xscale); -// ykernel.initKaiser(4, 1.0f / yscale); - - xkernel.debugPrint(); - - // @@ Select fastest filtering order: - // w * m_height <= h * m_width -> XY, else -> YX + // @@ Use monophase filters when frac(m_width / w) == 0 - AutoPtr tmp_image( new FloatImage() ); - tmp_image->allocate(m_componentNum, w, m_height); + PolyphaseKernel xkernel(filter, m_width, w); + PolyphaseKernel ykernel(filter, m_height, h); + AutoPtr tmp_image( new FloatImage() ); AutoPtr dst_image( new FloatImage() ); - dst_image->allocate(m_componentNum, w, h); - - Array tmp_column(h); - tmp_column.resize(h); - - for (uint c = 0; c < m_componentNum; c++) + + // @@ Select fastest filtering order: + //if (w * m_height <= h * m_width) { - float * tmp_channel = tmp_image->channel(c); + tmp_image->allocate(m_componentNum, w, m_height); + dst_image->allocate(m_componentNum, w, h); - for(uint y = 0; y < m_height; y++) { - this->applyKernelHorizontal(&xkernel, xscale, y, c, wm, tmp_channel + y * w); - } + Array tmp_column(h); + tmp_column.resize(h); - float * dst_channel = dst_image->channel(c); - - for (uint x = 0; x < w; x++) { - tmp_image->applyKernelVertical(&ykernel, yscale, x, c, wm, tmp_column.unsecureBuffer()); + for (uint c = 0; c < m_componentNum; c++) + { + float * tmp_channel = tmp_image->channel(c); + + for (uint y = 0; y < m_height; y++) { + this->applyKernelHorizontal(xkernel, y, c, wm, tmp_channel + y * w); + } - for(uint y = 0; y < h; y++) { - dst_channel[y * w + x] = tmp_column[y]; + float * dst_channel = dst_image->channel(c); + + for (uint x = 0; x < w; x++) { + tmp_image->applyKernelVertical(ykernel, x, c, wm, tmp_column.unsecureBuffer()); + + for (uint y = 0; y < h; y++) { + dst_channel[y * w + x] = tmp_column[y]; + } } } } + /*else + { + tmp_image->allocate(m_componentNum, m_width, h); + dst_image->allocate(m_componentNum, w, h); + + Array tmp_column(h); + tmp_column.resize(h); + + for (uint c = 0; c < m_componentNum; c++) + { + float * tmp_channel = tmp_image->channel(c); + + for (uint x = 0; x < w; x++) { + tmp_image->applyKernelVertical(ykernel, x, c, wm, tmp_column.unsecureBuffer()); + + for (uint y = 0; y < h; y++) { + tmp_channel[y * w + x] = tmp_column[y]; + } + } + + float * dst_channel = dst_image->channel(c); + + for (uint y = 0; y < m_height; y++) { + this->applyKernelHorizontal(xkernel, y, c, wm, dst_channel + y * w); + } + } + }*/ return dst_image.release(); } @@ -722,15 +746,46 @@ float FloatImage::applyKernelHorizontal(const Kernel1 * k, int x, int y, int c, /// Apply 1D vertical kernel at the given coordinates and return result. -void FloatImage::applyKernelVertical(const PolyphaseKernel * k, float scale, int x, int c, WrapMode wm, float * output) const +void FloatImage::applyKernelVertical(const PolyphaseKernel & k, int x, int c, WrapMode wm, float * output) const { - nvDebugCheck(k != NULL); - + uint length = k.length(); + float scale = float(length) / float(m_height); + float iscale = 1.0f / scale; + + float width = k.width(); + float windowSize = k.windowSize(); + + const float * channel = this->channel(c); + + for (uint i = 0; i < length; i++) + { + const float center = (0.5f + i) * iscale; + + int left = floor(center - width); + int right = ceil(center + width); + + nvCheck(right - left <= windowSize); + + float sum = 0; + for (int j = 0; j < windowSize; ++j) + { + const int idx = this->index(x, j+left, wm); + + sum += k.valueAt(i, j) * channel[idx]; + } + + output[i] = sum; + } + + /* const float kernelWidth = k->width(); const float kernelOffset = kernelWidth * 0.5f; const int kernelLength = k->length(); const int kernelWindow = k->windowSize(); + //const float offset = 0.5f * scale * (1 - kw); + const float offset = (0.5f * scale) - kernelOffset; + const float * channel = this->channel(c); for (int y = 0; y < kernelLength; y++) @@ -738,7 +793,7 @@ void FloatImage::applyKernelVertical(const PolyphaseKernel * k, float scale, int float sum = 0.0f; for (int i = 0; i < kernelWindow; i++) { - const int src_y = int(y * scale) + i; + const int src_y = int(y * scale + offset) + i; const int idx = this->index(x, src_y, wm); sum += k->valueAt(y, i) * channel[idx]; @@ -746,18 +801,48 @@ void FloatImage::applyKernelVertical(const PolyphaseKernel * k, float scale, int output[y] = sum; } + */ } /// Apply 1D horizontal kernel at the given coordinates and return result. -void FloatImage::applyKernelHorizontal(const PolyphaseKernel * k, float scale, int y, int c, WrapMode wm, float * output) const +void FloatImage::applyKernelHorizontal(const PolyphaseKernel & k, int y, int c, WrapMode wm, float * output) const { - nvDebugCheck(k != NULL); - + uint length = k.length(); + float scale = float(length) / float(m_width); + float iscale = 1.0f / scale; + + float width = k.width(); + float windowSize = k.windowSize(); + + const float * channel = this->channel(c); + + for (uint i = 0; i < length; i++) + { + const float center = (0.5f + i) * iscale; + + int left = floor(center - width); + int right = ceil(center + width); + nvCheck(right - left <= (int)windowSize); + + float sum = 0; + for (int j = 0; j < windowSize; ++j) + { + const int idx = this->index(left + j, y, wm); + + sum += k.valueAt(i, j) * channel[idx]; + } + + output[i] = sum; + } + + /* const float kernelWidth = k->width(); const float kernelOffset = kernelWidth * 0.5f; const int kernelLength = k->length(); const int kernelWindow = k->windowSize(); - + + const float offset = (0.5f * scale) - kernelOffset; + const float * channel = this->channel(c); for (int x = 0; x < kernelLength; x++) @@ -765,7 +850,7 @@ void FloatImage::applyKernelHorizontal(const PolyphaseKernel * k, float scale, i float sum = 0.0f; for (int e = 0; e < kernelWindow; e++) { - const int src_x = int(x * scale) + e; + const int src_x = int(x * scale + offset) + e; const int idx = this->index(src_x, y, wm); sum += k->valueAt(x, e) * channel[idx]; @@ -773,5 +858,6 @@ void FloatImage::applyKernelHorizontal(const PolyphaseKernel * k, float scale, i output[x] = sum; } + */ } diff --git a/src/nvimage/FloatImage.h b/src/nvimage/FloatImage.h index bcc878d..83e74df 100644 --- a/src/nvimage/FloatImage.h +++ b/src/nvimage/FloatImage.h @@ -10,6 +10,7 @@ namespace nv { class Image; +class Filter; class Kernel1; class Kernel2; class PolyphaseKernel; @@ -61,18 +62,18 @@ public: NVIMAGE_API FloatImage * fastDownSample() const; - NVIMAGE_API FloatImage * downSample(const Kernel1 & filter, WrapMode wm) const; - NVIMAGE_API FloatImage * downSample(const Kernel1 & filter, uint w, uint h, WrapMode wm) const; - - // experimental polyphase filter: - NVIMAGE_API FloatImage * downSample(uint w, uint h, WrapMode wm) const; + NVIMAGE_API FloatImage * downSample(const Filter & filter, WrapMode wm) const; + NVIMAGE_API FloatImage * downSample(const Filter & filter, uint w, uint h, WrapMode wm) const; + + //NVIMAGE_API FloatImage * downSample(const Kernel1 & filter, WrapMode wm) const; + //NVIMAGE_API FloatImage * downSample(const Kernel1 & filter, uint w, uint h, WrapMode wm) const; //@} NVIMAGE_API float applyKernel(const Kernel2 * k, int x, int y, int c, WrapMode wm) const; NVIMAGE_API float applyKernelVertical(const Kernel1 * k, int x, int y, int c, WrapMode wm) const; NVIMAGE_API float applyKernelHorizontal(const Kernel1 * k, int x, int y, int c, WrapMode wm) const; - NVIMAGE_API void applyKernelVertical(const PolyphaseKernel * k, float scale, int x, int c, WrapMode wm, float * output) const; - NVIMAGE_API void applyKernelHorizontal(const PolyphaseKernel * k, float scale, int y, int c, WrapMode wm, float * output) const; + NVIMAGE_API void applyKernelVertical(const PolyphaseKernel & k, int x, int c, WrapMode wm, float * output) const; + NVIMAGE_API void applyKernelHorizontal(const PolyphaseKernel & k, int y, int c, WrapMode wm, float * output) const; uint width() const { return m_width; } @@ -225,17 +226,16 @@ inline uint FloatImage::indexRepeat(int x, int y) const return repeat_remainder(y, m_height) * m_width + repeat_remainder(x, m_width); } -// @@ This could be way more efficient. inline uint FloatImage::indexMirror(int x, int y) const { - while ((x < 0) || (x > (m_width - 1))) { - if (x < 0) x = -x; - if (x >= m_width) x = m_width + m_width - x - 1; + x = abs(x); + while (x >= m_width) { + x = m_width + m_width - x - 2; } - while ((y < 0) || (y > (m_height - 1))) { - if (y < 0) y = -y; - if (y >= m_height) y = m_height + m_height - y - 1; + y = abs(y); + while (y >= m_height) { + y = m_height + m_height - y - 2; } return index(x, y); diff --git a/src/nvtt/nvtt.cpp b/src/nvtt/nvtt.cpp index 74fde81..e47a18f 100644 --- a/src/nvtt/nvtt.cpp +++ b/src/nvtt/nvtt.cpp @@ -325,15 +325,15 @@ static FloatImage * createMipmap(const FloatImage * floatImage, const InputOptio } else if (inputOptions.mipmapFilter == MipmapFilter_Triangle) { - Kernel1 kernel(4); - kernel.initFilter(Filter::Triangle); - result = floatImage->downSample(kernel, (FloatImage::WrapMode)inputOptions.wrapMode); + TriangleFilter filter; + result = floatImage->downSample(filter, (FloatImage::WrapMode)inputOptions.wrapMode); } else /*if (inputOptions.mipmapFilter == MipmapFilter_Kaiser)*/ { - Kernel1 kernel(inputOptions.kaiserWidth); - kernel.initKaiser(inputOptions.kaiserAlpha, inputOptions.kaiserStretch); - result = floatImage->downSample(kernel, (FloatImage::WrapMode)inputOptions.wrapMode); + nvDebugCheck(inputOptions.mipmapFilter == MipmapFilter_Kaiser); + KaiserFilter filter(inputOptions.kaiserWidth); + filter.setParameters(inputOptions.kaiserAlpha, inputOptions.kaiserStretch); + result = floatImage->downSample(filter, (FloatImage::WrapMode)inputOptions.wrapMode); } // Normalize mipmap. diff --git a/src/nvtt/tools/resize.cpp b/src/nvtt/tools/resize.cpp index 89dbe96..8e817c9 100644 --- a/src/nvtt/tools/resize.cpp +++ b/src/nvtt/tools/resize.cpp @@ -68,7 +68,7 @@ static bool loadImage(nv::Image & image, const char * fileName) int main(int argc, char *argv[]) { - MyAssertHandler assertHandler; + //MyAssertHandler assertHandler; MyMessageHandler messageHandler; float scale = 0.5f; @@ -122,7 +122,8 @@ int main(int argc, char *argv[]) // k.initKaiser(4, scale, 20); // nv::AutoPtr fresult(fimage.downSample(k, image.width() * scale, image.height() * scale, nv::FloatImage::WrapMode_Clamp)); - nv::AutoPtr fresult(fimage.downSample(image.width() * scale, image.height() * scale, nv::FloatImage::WrapMode_Mirror)); + nv::BoxFilter filter; + nv::AutoPtr fresult(fimage.downSample(filter, image.width() * scale, image.height() * scale, nv::FloatImage::WrapMode_Mirror)); nv::AutoPtr result(fresult->createImageGammaCorrect(1.0));