From c772a00b8f57371f42b571d4cee518168d2244cc Mon Sep 17 00:00:00 2001 From: castano Date: Wed, 28 Nov 2007 10:34:40 +0000 Subject: [PATCH] More progress with polyphase filters. --- src/nvimage/Filter.cpp | 252 +++++++++++++++++++++---------------- src/nvimage/Filter.h | 28 +++-- src/nvimage/FloatImage.cpp | 38 +++--- 3 files changed, 177 insertions(+), 141 deletions(-) diff --git a/src/nvimage/Filter.cpp b/src/nvimage/Filter.cpp index a8582bc..583e47f 100644 --- a/src/nvimage/Filter.cpp +++ b/src/nvimage/Filter.cpp @@ -203,47 +203,62 @@ static Filter s_filter_array[] = { {filter_kaiser, 1.0f}, // Kaiser }; + +inline static float sampleFilter(Filter::Function func, float x, float scale, int samples) +{ + float sum = 0; + + for(int s = 0; s < samples; s++) + { + sum += func((x + (float(s) + 0.5f) * (1.0f / float(samples))) * scale); + } + + return sum; +} + + + } // namespace /// Ctor. -Kernel1::Kernel1(uint width) : w(width) +Kernel1::Kernel1(uint ws) : m_windowSize(ws) { - data = new float[w]; + m_data = new float[m_windowSize]; } /// Copy ctor. -Kernel1::Kernel1(const Kernel1 & k) : w(k.w) +Kernel1::Kernel1(const Kernel1 & k) : m_windowSize(k.m_windowSize) { - data = new float[w]; - for(uint i = 0; i < w; i++) { - data[i] = k.data[i]; + m_data = new float[m_windowSize]; + for(uint i = 0; i < m_windowSize; i++) { + m_data[i] = k.m_data[i]; } } /// Dtor. Kernel1::~Kernel1() { - delete data; + delete m_data; } /// Normalize the filter. void Kernel1::normalize() { float total = 0.0f; - for(uint i = 0; i < w; i++) { - total += data[i]; + for(uint i = 0; i < m_windowSize; i++) { + total += m_data[i]; } float inv = 1.0f / total; - for(uint i = 0; i < w; i++) { - data[i] *= inv; + for(uint i = 0; i < m_windowSize; i++) { + m_data[i] *= inv; } } -/// Init 1D Box filter. +/// Init 1D filter. void Kernel1::initFilter(Filter::Enum f, int samples /*= 1*/) { nvDebugCheck(f < Filter::Num); @@ -252,21 +267,12 @@ void Kernel1::initFilter(Filter::Enum f, int samples /*= 1*/) float (* filter_function)(float) = s_filter_array[f].function; const float support = s_filter_array[f].support; - const float half_width = float(w) / 2; - const float offset = -half_width; - - const float s_scale = 1.0f / float(samples); - const float x_scale = support / half_width; + const float halfWindowSize = float(m_windowSize) / 2.0f; + const float scale = support / halfWindowSize; - for(uint i = 0; i < w; i++) + for(uint i = 0; i < m_windowSize; i++) { - float sum = 0; - for(int s = 0; s < samples; s++) - { - float x = i + offset + (s + 0.5f) * s_scale; - sum += filter_function(x * x_scale); - } - data[i] = sum; + m_data[i] = sampleFilter(filter_function, i - halfWindowSize, scale, samples); } normalize(); @@ -276,13 +282,12 @@ void Kernel1::initFilter(Filter::Enum f, int samples /*= 1*/) /// Init 1D sinc filter. void Kernel1::initSinc(float stretch /*= 1*/) { - const float half_width = float(w) / 2; - const float offset = -half_width; + const float halfWindowSize = float(m_windowSize) / 2; const float nudge = 0.5f; - for(uint i = 0; i < w; i++) { - const float x = (i + offset) + nudge; - data[i] = sincf(PI * x * stretch); + for(uint i = 0; i < m_windowSize; i++) { + const float x = (i - halfWindowSize) + nudge; + m_data[i] = sincf(PI * x * stretch); } normalize(); @@ -292,25 +297,24 @@ void Kernel1::initSinc(float stretch /*= 1*/) /// Init 1D Kaiser-windowed sinc filter. void Kernel1::initKaiser(float alpha /*= 4*/, float stretch /*= 0.5*/, int samples/*= 1*/) { - const float half_width = float(w) / 2; - const float offset = -half_width; + const float halfWindowSize = float(m_windowSize) / 2; const float s_scale = 1.0f / float(samples); - const float x_scale = 1.0f / half_width; + const float x_scale = 1.0f / halfWindowSize; - for(uint i = 0; i < w; i++) + for(uint i = 0; i < m_windowSize; i++) { float sum = 0; for(int s = 0; s < samples; s++) { - float x = i + offset + (s + 0.5f) * s_scale; + 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; } - data[i] = sum; + m_data[i] = sum; } normalize(); @@ -320,13 +324,12 @@ void Kernel1::initKaiser(float alpha /*= 4*/, float stretch /*= 0.5*/, int sampl /// Init 1D Mitchell filter. void Kernel1::initMitchell(float b, float c) { - const float half_width = float(w) / 2; - const float offset = -half_width; + const float halfWindowSize = float(m_windowSize) / 2; const float nudge = 0.5f; - for(uint i = 0; i < w; i++) { - const float x = (i + offset) + nudge; - data[i] = filter_mitchell(x / half_width, b, c); + for (uint i = 0; i < m_windowSize; i++) { + const float x = (i - halfWindowSize) + nudge; + m_data[i] = filter_mitchell(x / halfWindowSize, b, c); } normalize(); @@ -336,25 +339,25 @@ void Kernel1::initMitchell(float b, float c) /// Print the kernel for debugging purposes. void Kernel1::debugPrint() { - for(uint i = 0; i < w; i++) { - nvDebug("%d: %f\n", i, data[i]); + for (uint i = 0; i < m_windowSize; i++) { + nvDebug("%d: %f\n", i, m_data[i]); } } /// Ctor. -Kernel2::Kernel2(uint width) : w(width) +Kernel2::Kernel2(uint ws) : m_windowSize(ws) { - data = new float[w*w]; + m_data = new float[m_windowSize * m_windowSize]; } /// Copy ctor. -Kernel2::Kernel2(const Kernel2 & k) : w(k.w) +Kernel2::Kernel2(const Kernel2 & k) : m_windowSize(k.m_windowSize) { - data = new float[w*w]; - for(uint i = 0; i < w*w; i++) { - data[i] = k.data[i]; + m_data = new float[m_windowSize * m_windowSize]; + for (uint i = 0; i < m_windowSize * m_windowSize; i++) { + m_data[i] = k.m_data[i]; } } @@ -362,29 +365,29 @@ Kernel2::Kernel2(const Kernel2 & k) : w(k.w) /// Dtor. Kernel2::~Kernel2() { - delete data; + delete m_data; } /// Normalize the filter. void Kernel2::normalize() { float total = 0.0f; - for(uint i = 0; i < w*w; i++) { - total += fabs(data[i]); + for(uint i = 0; i < m_windowSize*m_windowSize; i++) { + total += fabs(m_data[i]); } float inv = 1.0f / total; - for(uint i = 0; i < w*w; i++) { - data[i] *= inv; + for(uint i = 0; i < m_windowSize*m_windowSize; i++) { + m_data[i] *= inv; } } /// Transpose the kernel. void Kernel2::transpose() { - for(uint i = 0; i < w; i++) { - for(uint j = i+1; j < w; j++) { - swap(data[i*w + j], data[j*w + i]); + for(uint i = 0; i < m_windowSize; i++) { + for(uint j = i+1; j < m_windowSize; j++) { + swap(m_data[i*m_windowSize + j], m_data[j*m_windowSize + i]); } } } @@ -392,40 +395,40 @@ void Kernel2::transpose() /// Init laplacian filter, usually used for sharpening. void Kernel2::initLaplacian() { - nvDebugCheck(w == 3); -// data[0] = -1; data[1] = -1; data[2] = -1; -// data[3] = -1; data[4] = +8; data[5] = -1; -// data[6] = -1; data[7] = -1; data[8] = -1; + nvDebugCheck(m_windowSize == 3); +// m_data[0] = -1; m_data[1] = -1; m_data[2] = -1; +// m_data[3] = -1; m_data[4] = +8; m_data[5] = -1; +// m_data[6] = -1; m_data[7] = -1; m_data[8] = -1; - data[0] = +0; data[1] = -1; data[2] = +0; - data[3] = -1; data[4] = +4; data[5] = -1; - data[6] = +0; data[7] = -1; data[8] = +0; + m_data[0] = +0; m_data[1] = -1; m_data[2] = +0; + m_data[3] = -1; m_data[4] = +4; m_data[5] = -1; + m_data[6] = +0; m_data[7] = -1; m_data[8] = +0; -// data[0] = +1; data[1] = -2; data[2] = +1; -// data[3] = -2; data[4] = +4; data[5] = -2; -// data[6] = +1; data[7] = -2; data[8] = +1; +// m_data[0] = +1; m_data[1] = -2; m_data[2] = +1; +// m_data[3] = -2; m_data[4] = +4; m_data[5] = -2; +// m_data[6] = +1; m_data[7] = -2; m_data[8] = +1; } /// Init simple edge detection filter. void Kernel2::initEdgeDetection() { - nvCheck(w == 3); - data[0] = 0; data[1] = 0; data[2] = 0; - data[3] = -1; data[4] = 0; data[5] = 1; - data[6] = 0; data[7] = 0; data[8] = 0; + nvCheck(m_windowSize == 3); + m_data[0] = 0; m_data[1] = 0; m_data[2] = 0; + m_data[3] =-1; m_data[4] = 0; m_data[5] = 1; + m_data[6] = 0; m_data[7] = 0; m_data[8] = 0; } /// Init sobel filter. void Kernel2::initSobel() { - if (w == 3) + if (m_windowSize == 3) { - data[0] = -1; data[1] = 0; data[2] = 1; - data[3] = -2; data[4] = 0; data[5] = 2; - data[6] = -1; data[7] = 0; data[8] = 1; + m_data[0] = -1; m_data[1] = 0; m_data[2] = 1; + m_data[3] = -2; m_data[4] = 0; m_data[5] = 2; + m_data[6] = -1; m_data[7] = 0; m_data[8] = 1; } - else if (w == 5) + else if (m_windowSize == 5) { float elements[] = { -1, -2, 0, 2, 1, @@ -436,10 +439,10 @@ void Kernel2::initSobel() }; for (int i = 0; i < 5*5; i++) { - data[i] = elements[i]; + m_data[i] = elements[i]; } } - else if (w == 7) + else if (m_windowSize == 7) { float elements[] = { -1, -2, -3, 0, 3, 2, 1, @@ -452,10 +455,10 @@ void Kernel2::initSobel() }; for (int i = 0; i < 7*7; i++) { - data[i] = elements[i]; + m_data[i] = elements[i]; } } - else if (w == 9) + else if (m_windowSize == 9) { float elements[] = { -1, -2, -3, -4, 0, 4, 3, 2, 1, @@ -470,7 +473,7 @@ void Kernel2::initSobel() }; for (int i = 0; i < 9*9; i++) { - data[i] = elements[i]; + m_data[i] = elements[i]; } } } @@ -478,13 +481,13 @@ void Kernel2::initSobel() /// Init prewitt filter. void Kernel2::initPrewitt() { - if (w == 3) + if (m_windowSize == 3) { - data[0] = -1; data[1] = 0; data[2] = -1; - data[3] = -1; data[4] = 0; data[5] = -1; - data[6] = -1; data[7] = 0; data[8] = -1; + m_data[0] = -1; m_data[1] = 0; m_data[2] = -1; + m_data[3] = -1; m_data[4] = 0; m_data[5] = -1; + m_data[6] = -1; m_data[7] = 0; m_data[8] = -1; } - else if (w == 5) + else if (m_windowSize == 5) { // @@ Is this correct? float elements[] = { @@ -496,7 +499,7 @@ void Kernel2::initPrewitt() }; for (int i = 0; i < 5*5; i++) { - data[i] = elements[i]; + m_data[i] = elements[i]; } } } @@ -504,7 +507,7 @@ void Kernel2::initPrewitt() /// Init blended sobel filter. void Kernel2::initBlendedSobel(const Vector4 & scale) { - nvCheck(w == 9); + nvCheck(m_windowSize == 9); { const float elements[] = { @@ -520,7 +523,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale) }; for (int i = 0; i < 9*9; i++) { - data[i] = elements[i] * scale.w(); + m_data[i] = elements[i] * scale.w(); } } { @@ -536,7 +539,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale) for (int i = 0; i < 7; i++) { for (int e = 0; e < 7; e++) { - data[i * 9 + e + 1] += elements[i * 7 + e] * scale.z(); + m_data[i * 9 + e + 1] += elements[i * 7 + e] * scale.z(); } } } @@ -551,7 +554,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale) for (int i = 0; i < 5; i++) { for (int e = 0; e < 5; e++) { - data[i * 9 + e + 2] += elements[i * 5 + e] * scale.y(); + m_data[i * 9 + e + 2] += elements[i * 5 + e] * scale.y(); } } } @@ -564,7 +567,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale) for (int i = 0; i < 3; i++) { for (int e = 0; e < 3; e++) { - data[i * 9 + e + 3] += elements[i * 3 + e] * scale.x(); + m_data[i * 9 + e + 3] += elements[i * 3 + e] * scale.x(); } } } @@ -582,6 +585,7 @@ static bool isMonoPhase(float w) } + PolyphaseKernel::PolyphaseKernel(float w, uint l) : m_width(w), m_size(ceilf(w) + 1), @@ -605,42 +609,70 @@ PolyphaseKernel::~PolyphaseKernel() delete [] m_data; } + +/* @@ 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); + + ... + } +*/ + void PolyphaseKernel::initFilter(Filter::Enum f, int samples/*= 1*/) { nvCheck(f < Filter::Num); - nvDebug("Init Filter:\n"); - nvDebug(" width = %f\n", m_width); - 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 offset = -half_width; - - const float s_scale = 1.0f / float(samples); - const float x_scale = support / half_width; + const float scale = support / half_width; for (uint j = 0; j < m_length; j++) { const float phase = frac(m_width * j); + const float offset = half_width + phase; + nvDebug("%d: ", j); + float total = 0.0f; for (uint i = 0; i < m_size; i++) { - const float x = i + offset - phase; - - float sum = 0; - for(int s = 0; s < samples; s++) - { - const float xx = x + (s + 0.5f) * s_scale; - sum += filter_function(xx * x_scale); - } - m_data[j * m_size + i] = sum; - - total += sum; + float sample = sampleFilter(filter_function, i - offset, scale, samples); + + nvDebug("(%5.3f | %d) ", sample, j + i - m_size/2); + + m_data[j * m_size + i] = sample; + total += sample; } + + nvDebug("\n"); + // normalize weights. for (uint i = 0; i < m_size; i++) { m_data[j * m_size + i] /= total; diff --git a/src/nvimage/Filter.h b/src/nvimage/Filter.h index 31943ce..4e160ae 100644 --- a/src/nvimage/Filter.h +++ b/src/nvimage/Filter.h @@ -26,7 +26,9 @@ namespace nv Num }; - float (*function)(float x); + typedef float (* Function)(float); + + Function function; float support; }; @@ -35,18 +37,18 @@ namespace nv class Kernel1 { public: - NVIMAGE_API Kernel1(uint width); + NVIMAGE_API Kernel1(uint windowSize); NVIMAGE_API Kernel1(const Kernel1 & k); NVIMAGE_API ~Kernel1(); NVIMAGE_API void normalize(); float valueAt(uint x) const { - return data[x]; + return m_data[x]; } - uint width() const { - return w; + uint windowSize() const { + return m_windowSize; } NVIMAGE_API void initFilter(Filter::Enum filter, int samples = 1); @@ -57,8 +59,8 @@ namespace nv NVIMAGE_API void debugPrint(); private: - const uint w; - float * data; + const uint m_windowSize; + float * m_data; }; @@ -74,11 +76,11 @@ namespace nv NVIMAGE_API void transpose(); float valueAt(uint x, uint y) const { - return data[y * w + x]; + return m_data[y * m_windowSize + x]; } - uint width() const { - return w; + uint windowSize() const { + return m_windowSize; } NVIMAGE_API void initLaplacian(); @@ -89,8 +91,8 @@ namespace nv NVIMAGE_API void initBlendedSobel(const Vector4 & scale); private: - const uint w; - float * data; + const uint m_windowSize; + float * m_data; }; /// A 1D polyphase kernel @@ -109,7 +111,7 @@ namespace nv return m_width; } - uint size() const { + uint windowSize() const { return m_size; } diff --git a/src/nvimage/FloatImage.cpp b/src/nvimage/FloatImage.cpp index 0c3616f..d720c8d 100644 --- a/src/nvimage/FloatImage.cpp +++ b/src/nvimage/FloatImage.cpp @@ -549,7 +549,7 @@ FloatImage * FloatImage::downSample(const Kernel1 & kernel, WrapMode wm) const /// Downsample applying a 1D kernel separately in each dimension. FloatImage * FloatImage::downSample(const Kernel1 & kernel, uint w, uint h, WrapMode wm) const { - nvCheck(!(kernel.width() & 1)); // Make sure that kernel m_width is even. + nvCheck(!(kernel.windowSize() & 1)); // Make sure that kernel m_width is even. AutoPtr tmp_image( new FloatImage() ); tmp_image->allocate(m_componentNum, w, m_height); @@ -611,6 +611,9 @@ FloatImage * FloatImage::downSample(uint w, uint h, WrapMode wm) const xkernel.debugPrint(); + // @@ Select fastest filtering order: + // w * m_height <= h * m_width -> XY, else -> YX + AutoPtr tmp_image( new FloatImage() ); tmp_image->allocate(m_componentNum, w, m_height); @@ -620,7 +623,7 @@ FloatImage * FloatImage::downSample(uint w, uint h, WrapMode wm) const Array tmp_column(h); tmp_column.resize(h); - for(uint c = 0; c < m_componentNum; c++) + for (uint c = 0; c < m_componentNum; c++) { float * tmp_channel = tmp_image->channel(c); @@ -630,7 +633,7 @@ FloatImage * FloatImage::downSample(uint w, uint h, WrapMode wm) const float * dst_channel = dst_image->channel(c); - for(uint x = 0; x < w; x++) { + for (uint x = 0; x < w; x++) { tmp_image->applyKernelVertical(&ykernel, yscale, x, c, wm, tmp_column.unsecureBuffer()); for(uint y = 0; y < h; y++) { @@ -639,7 +642,6 @@ FloatImage * FloatImage::downSample(uint w, uint h, WrapMode wm) const } } - //return tmp_image.release(); return dst_image.release(); } @@ -650,17 +652,17 @@ float FloatImage::applyKernel(const Kernel2 * k, int x, int y, int c, WrapMode w { nvDebugCheck(k != NULL); - const uint kernelWidth = k->width(); - const int kernelOffset = int(kernelWidth / 2) - 1; + const uint kernelWindow = k->windowSize(); + const int kernelOffset = int(kernelWindow / 2) - 1; const float * channel = this->channel(c); float sum = 0.0f; - for(uint i = 0; i < kernelWidth; i++) + for (uint i = 0; i < kernelWindow; i++) { const int src_y = int(y + i) - kernelOffset; - for(uint e = 0; e < kernelWidth; e++) + for (uint e = 0; e < kernelWindow; e++) { const int src_x = int(x + e) - kernelOffset; @@ -679,13 +681,13 @@ float FloatImage::applyKernelVertical(const Kernel1 * k, int x, int y, int c, Wr { nvDebugCheck(k != NULL); - const uint kernelWidth = k->width(); - const int kernelOffset = int(kernelWidth / 2) - 1; + const uint kernelWindow = k->windowSize(); + const int kernelOffset = int(kernelWindow / 2) - 1; const float * channel = this->channel(c); float sum = 0.0f; - for(uint i = 0; i < kernelWidth; i++) + for (uint i = 0; i < kernelWindow; i++) { const int src_y = int(y + i) - kernelOffset; const int idx = this->index(x, src_y, wm); @@ -701,13 +703,13 @@ float FloatImage::applyKernelHorizontal(const Kernel1 * k, int x, int y, int c, { nvDebugCheck(k != NULL); - const uint kernelWidth = k->width(); - const int kernelOffset = int(kernelWidth / 2) - 1; + const uint kernelWindow = k->windowSize(); + const int kernelOffset = int(kernelWindow / 2) - 1; const float * channel = this->channel(c); float sum = 0.0f; - for(uint e = 0; e < kernelWidth; e++) + for (uint e = 0; e < kernelWindow; e++) { const int src_x = int(x + e) - kernelOffset; const int idx = this->index(src_x, y, wm); @@ -725,9 +727,9 @@ void FloatImage::applyKernelVertical(const PolyphaseKernel * k, float scale, int nvDebugCheck(k != NULL); const float kernelWidth = k->width(); - const float kernelOffset = scale * 0.5f; + const float kernelOffset = kernelWidth * 0.5f; const int kernelLength = k->length(); - const int kernelWindow = k->size(); + const int kernelWindow = k->windowSize(); const float * channel = this->channel(c); @@ -752,9 +754,9 @@ void FloatImage::applyKernelHorizontal(const PolyphaseKernel * k, float scale, i nvDebugCheck(k != NULL); const float kernelWidth = k->width(); - const float kernelOffset = scale * 0.5f; + const float kernelOffset = kernelWidth * 0.5f; const int kernelLength = k->length(); - const int kernelWindow = k->size(); + const int kernelWindow = k->windowSize(); const float * channel = this->channel(c);