diff --git a/src/nvimage/Filter.cpp b/src/nvimage/Filter.cpp index 37e35e5..3d93a6f 100644 --- a/src/nvimage/Filter.cpp +++ b/src/nvimage/Filter.cpp @@ -100,6 +100,10 @@ Filter::Filter(float width) : m_width(width) { } +/*virtual*/ Filter::~Filter() +{ +} + float Filter::sample(float x, float scale, int samples) const { // return evaluate(x * scale); @@ -223,11 +227,9 @@ KaiserFilter::KaiserFilter(float w) : Filter(w) { setParameters(4.0f, 1.0f); } float KaiserFilter::evaluate(float x) const { 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; + const float t = x / m_width; + if ((1 - t * t) >= 0) return sinc_value * bessel0(alpha * sqrtf(1 - t * t)) / bessel0(alpha); + else return 0; } void KaiserFilter::setParameters(float alpha, float stretch) @@ -239,17 +241,31 @@ void KaiserFilter::setParameters(float alpha, float stretch) /// Ctor. -Kernel1::Kernel1(uint ws) : m_windowSize(ws) +Kernel1::Kernel1(const Filter & f, int iscale, int samples/*= 32*/) { + nvDebugCheck(iscale > 1); + nvDebugCheck(samples > 0); + + const float scale = 1.0f / iscale; + + m_width = f.width() * iscale; + m_windowSize = ceilf(2 * m_width); m_data = new float[m_windowSize]; -} - -/// Copy ctor. -Kernel1::Kernel1(const Kernel1 & k) : m_windowSize(k.m_windowSize) -{ - m_data = new float[m_windowSize]; - for(uint i = 0; i < m_windowSize; i++) { - m_data[i] = k.m_data[i]; + + const float offset = float(m_windowSize) / 2; + + float total = 0.0f; + for (int i = 0; i < m_windowSize; i++) + { + const float sample = f.sample(i - offset, scale, samples); + m_data[i] = sample; + total += sample; + } + + const float inv = 1.0f / total; + for (int i = 0; i < m_windowSize; i++) + { + m_data[i] *= inv; } } @@ -259,103 +275,10 @@ Kernel1::~Kernel1() delete m_data; } -/// Normalize the filter. -void Kernel1::normalize() -{ - float total = 0.0f; - for(uint i = 0; i < m_windowSize; i++) { - total += m_data[i]; - } - - float inv = 1.0f / total; - for(uint i = 0; i < m_windowSize; i++) { - m_data[i] *= inv; - } -} - -#if 0 -/// Init 1D filter. -void Kernel1::initFilter(Filter::Enum f, int samples /*= 1*/) -{ - nvDebugCheck(f < Filter::Num); - nvDebugCheck(samples >= 1); - - float (* filter_function)(float) = s_filter_array[f].function; - const float support = s_filter_array[f].support; - - const float halfWindowSize = float(m_windowSize) / 2.0f; - const float scale = support / halfWindowSize; - - for(uint i = 0; i < m_windowSize; i++) - { - m_data[i] = sampleFilter(filter_function, i - halfWindowSize, scale, samples); - } - - normalize(); -} - - -/// Init 1D sinc filter. -void Kernel1::initSinc(float stretch /*= 1*/) -{ - const float halfWindowSize = float(m_windowSize) / 2; - const float nudge = 0.5f; - - for(uint i = 0; i < m_windowSize; i++) { - const float x = (i - halfWindowSize) + nudge; - m_data[i] = sincf(PI * x * stretch); - } - - normalize(); -} - - -/// Init 1D Kaiser-windowed sinc filter. -void Kernel1::initKaiser(float alpha /*= 4*/, float stretch /*= 0.5*/, int samples/*= 1*/) -{ - const float halfWindowSize = float(m_windowSize) / 2; - - const float s_scale = 1.0f / float(samples); - const float x_scale = 1.0f / halfWindowSize; - - for(uint i = 0; i < m_windowSize; i++) - { - float sum = 0; - for(int s = 0; s < samples; s++) - { - float x = i - halfWindowSize + (s + 0.5f) * s_scale; - - const float sinc_value = sincf(PI * x * stretch); - const float window_value = filter_kaiser(x * x_scale, alpha); // @@ should the window be streched? I don't think so. - - sum += sinc_value * window_value; - } - m_data[i] = sum; - } - - normalize(); -} - - -/// Init 1D Mitchell filter. -void Kernel1::initMitchell(float b, float c) -{ - const float halfWindowSize = float(m_windowSize) / 2; - const float nudge = 0.5f; - - for (uint i = 0; i < m_windowSize; i++) { - const float x = (i - halfWindowSize) + nudge; - m_data[i] = filter_mitchell(x / halfWindowSize, b, c); - } - - normalize(); -} -#endif - /// Print the kernel for debugging purposes. void Kernel1::debugPrint() { - for (uint i = 0; i < m_windowSize; i++) { + for (int i = 0; i < m_windowSize; i++) { nvDebug("%d: %f\n", i, m_data[i]); } } @@ -590,41 +513,13 @@ void Kernel2::initBlendedSobel(const Vector4 & scale) } -static float frac(float f) +PolyphaseKernel::PolyphaseKernel(const Filter & f, uint srcLength, uint dstLength, int samples/*= 32*/) { - return f - floorf(f); -} - -static bool isMonoPhase(float w) -{ - return isZero(frac(w)); -} - - -/* -PolyphaseKernel::PolyphaseKernel(float w, uint l) : - m_width(w), - m_size(ceilf(w) + 1), - m_length(l) -{ - // size = width + (length - 1) * phase - m_data = new float[m_size * m_length]; -} - -PolyphaseKernel::PolyphaseKernel(const PolyphaseKernel & k) : - m_width(k.m_width), - m_size(k.m_size), - m_length(k.m_length) -{ - m_data = new float[m_size * m_length]; - memcpy(m_data, k.m_data, sizeof(float) * m_size * m_length); -} -*/ - -PolyphaseKernel::PolyphaseKernel(const Filter & f, uint srcLength, uint dstLength) -{ - float scale = float(dstLength) / float(srcLength); - float iscale = 1.0f / scale; + nvCheck(srcLength >= dstLength); // @@ Upsampling not implemented! + nvDebugCheck(samples > 0); + + const float scale = float(dstLength) / float(srcLength); + const float iscale = 1.0f / scale; m_length = dstLength; m_width = f.width() * iscale; @@ -637,19 +532,19 @@ PolyphaseKernel::PolyphaseKernel(const Filter & f, uint srcLength, uint dstLengt { const float center = (0.5f + i) * iscale; - int left = floor(center - m_width); - int right = ceil(center + m_width); - nvCheck(right - left <= (int)m_windowSize); - + const int left = floorf(center - m_width); + const int right = ceilf(center + m_width); + nvDebugCheck(right - left <= m_windowSize); + float total = 0.0f; for (int j = 0; j < m_windowSize; j++) { - float sample = f.sample(left + j - center, scale, 40); - + const float sample = f.sample(left + j - center, scale, samples); + m_data[i * m_windowSize + j] = sample; total += sample; } - + // normalize weights. for (int j = 0; j < m_windowSize; j++) { @@ -670,7 +565,7 @@ void PolyphaseKernel::debugPrint() const for (uint i = 0; i < m_length; i++) { nvDebug("%d: ", i); - for (uint j = 0; j < m_windowSize; j++) + for (int j = 0; j < m_windowSize; j++) { nvDebug(" %6.4f", m_data[i * m_windowSize + j]); } diff --git a/src/nvimage/Filter.h b/src/nvimage/Filter.h index 7ce1a7e..098d810 100644 --- a/src/nvimage/Filter.h +++ b/src/nvimage/Filter.h @@ -4,6 +4,7 @@ #define NV_IMAGE_FILTER_H #include +#include namespace nv { @@ -14,6 +15,7 @@ namespace nv { public: NVIMAGE_API Filter(float width); + NVIMAGE_API virtual ~Filter(); NVIMAGE_API float width() const { return m_width; } NVIMAGE_API float sample(float x, float scale, int samples) const; @@ -116,31 +118,29 @@ namespace nv /// A 1D kernel. Used to precompute filter weights. class Kernel1 { + NV_FORBID_COPY(Kernel1); public: - NVIMAGE_API Kernel1(uint windowSize); - NVIMAGE_API Kernel1(const Kernel1 & k); + NVIMAGE_API Kernel1(const Filter & f, int iscale, int samples = 32); NVIMAGE_API ~Kernel1(); - NVIMAGE_API void normalize(); - float valueAt(uint x) const { + nvDebugCheck(x < (uint)m_windowSize); return m_data[x]; } - uint windowSize() const { + int 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); - */ + + float width() const { + return m_width; + } NVIMAGE_API void debugPrint(); private: - const uint m_windowSize; + int m_windowSize; + float m_width; float * m_data; }; @@ -180,12 +180,12 @@ namespace nv /// A 1D polyphase kernel class PolyphaseKernel { - NV_FORBID_COPY(PolyphaseKernel) + NV_FORBID_COPY(PolyphaseKernel); public: - NVIMAGE_API PolyphaseKernel(const Filter & f, uint srcLength, uint dstLength); + NVIMAGE_API PolyphaseKernel(const Filter & f, uint srcLength, uint dstLength, int samples = 32); NVIMAGE_API ~PolyphaseKernel(); - uint windowSize() const { + int windowSize() const { return m_windowSize; } @@ -199,14 +199,14 @@ namespace nv float valueAt(uint column, uint x) const { nvDebugCheck(column < m_length); - nvDebugCheck(x < m_windowSize); + nvDebugCheck(x < (uint)m_windowSize); return m_data[column * m_windowSize + x]; } NVIMAGE_API void debugPrint() const; private: - uint m_windowSize; + int m_windowSize; uint m_length; float m_width; float * m_data; diff --git a/src/nvimage/FloatImage.cpp b/src/nvimage/FloatImage.cpp index 5dd812f..93ba781 100644 --- a/src/nvimage/FloatImage.cpp +++ b/src/nvimage/FloatImage.cpp @@ -599,17 +599,18 @@ FloatImage * FloatImage::downSample(const Filter & filter, WrapMode wm) const return downSample(filter, w, h, wm); } + /// Downsample applying a 1D kernel separately in each dimension. FloatImage * FloatImage::downSample(const Filter & filter, uint w, uint h, WrapMode wm) const { // @@ Use monophase filters when frac(m_width / w) == 0 - PolyphaseKernel xkernel(filter, m_width, w); - PolyphaseKernel ykernel(filter, m_height, h); - AutoPtr tmp_image( new FloatImage() ); AutoPtr dst_image( new FloatImage() ); - + + PolyphaseKernel xkernel(filter, m_width, w, 32); + PolyphaseKernel ykernel(filter, m_height, h, 32); + // @@ Select fastest filtering order: //if (w * m_height <= h * m_width) { @@ -748,116 +749,64 @@ 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, int x, int c, WrapMode wm, float * output) const { - uint length = k.length(); - float scale = float(length) / float(m_height); - float iscale = 1.0f / scale; + const uint length = k.length(); + const float scale = float(length) / float(m_height); + const float iscale = 1.0f / scale; - float width = k.width(); - float windowSize = k.windowSize(); + const float width = k.width(); + const int windowSize = k.windowSize(); const float * channel = this->channel(c); - for (uint i = 0; i < length; i++) - { + 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 left = floor(center - width); + const 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]; - } + + 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++) - { - float sum = 0.0f; - for (int i = 0; i < kernelWindow; 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]; - } - - output[y] = sum; - } - */ } /// Apply 1D horizontal kernel at the given coordinates and return result. void FloatImage::applyKernelHorizontal(const PolyphaseKernel & k, int y, int c, WrapMode wm, float * output) const { - uint length = k.length(); - float scale = float(length) / float(m_width); - float iscale = 1.0f / scale; + const uint length = k.length(); + const float scale = float(length) / float(m_width); + const float iscale = 1.0f / scale; - float width = k.width(); - float windowSize = k.windowSize(); + const float width = k.width(); + const int windowSize = k.windowSize(); const float * channel = this->channel(c); - for (uint i = 0; i < length; i++) - { + 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 left = floorf(center - width); + const int right = ceilf(center + width); + nvDebugCheck(right - left <= 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]; - } + + 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++) - { - float sum = 0.0f; - for (int e = 0; e < kernelWindow; 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]; - } - - output[x] = sum; - } - */ }