Correct polyphase filters.

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");
 	}

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;
 	};
 

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<FloatImage> 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<FloatImage> tmp_image( new FloatImage() );
 	AutoPtr<FloatImage> dst_image( new FloatImage() );	
-	dst_image->allocate(m_componentNum, w, h);
-	
-	Array<float> 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<float> 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<float> 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;
 	}
+	*/
 }
 

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);

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.

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<nv::FloatImage> fresult(fimage.downSample(k, image.width() * scale, image.height() * scale, nv::FloatImage::WrapMode_Clamp));
 	
-	nv::AutoPtr<nv::FloatImage> fresult(fimage.downSample(image.width() * scale, image.height() * scale, nv::FloatImage::WrapMode_Mirror));
+	nv::BoxFilter filter;
+	nv::AutoPtr<nv::FloatImage> fresult(fimage.downSample(filter, image.width() * scale, image.height() * scale, nv::FloatImage::WrapMode_Mirror));
 	
 	
 	nv::AutoPtr<nv::Image> result(fresult->createImageGammaCorrect(1.0));