Some progress with polyphase filters.

src/nvimage/Filter.cpp

@@ -109,7 +109,7 @@ inline float sincf( const float x )
 inline static float filter_lanczos3(float x)
 {
 	if( x < 0.0f ) x = -x;
-	if( x < 3.0f ) return(sincf(x) * sincf(x / 3.0f));
+	if( x < 3.0f ) return sincf(x) * sincf(x / 3.0f);
 	return 0.0f;
 }
 
@@ -165,7 +165,7 @@ static float bessel0(float x)
 	return sum;
 }
 
-// Alternative bessel function from Paul Heckbert.
+/*// Alternative bessel function from Paul Heckbert.
 static float _bessel0(float x)
 {
 	const float EPSILON_RATIO = 1E-6;
@@ -177,7 +177,7 @@ static float _bessel0(float x)
 		t *= y / float(i * i);
     }
     return sum;
-}
+}*/
 
 // support = 1.0
 inline static float filter_kaiser(float x, float alpha)
@@ -244,21 +244,29 @@ void Kernel1::normalize()
 
 
 /// Init 1D Box filter.
-void Kernel1::initFilter(Filter::Enum f)
+void Kernel1::initFilter(Filter::Enum f, int samples /*= 1*/)
 {
-	nvCheck((w & 1) == 0);
-	nvCheck(f < Filter::Num);
+	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 half_width = float(w / 2);
+	const float half_width = float(w) / 2;
 	const float offset = -half_width;
-	const float nudge = 0.5f;
+
+	const float s_scale = 1.0f / float(samples);
+	const float x_scale = support / half_width;
 	
-	for(uint i = 0; i < w; i++) {
-		const float x = (i + offset) + nudge;
-		data[i] = filter_function(x * support / half_width);
+	for(uint i = 0; i < w; 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;
 	}
 	
 	normalize();
@@ -268,9 +276,7 @@ void Kernel1::initFilter(Filter::Enum f)
 /// Init 1D sinc filter.
 void Kernel1::initSinc(float stretch /*= 1*/)
 {
-	nvCheck((w & 1) == 0);
-	
-	const float half_width = float(w / 2);
+	const float half_width = float(w) / 2;
 	const float offset = -half_width;
 	const float nudge = 0.5f;
 	
@@ -283,21 +289,28 @@ void Kernel1::initSinc(float stretch /*= 1*/)
 }
 
 
-/// Init 1D windowed Kaiser filter.
-void Kernel1::initKaiser(float alpha, float stretch /*= 1*/)
+/// Init 1D Kaiser-windowed sinc filter.
+void Kernel1::initKaiser(float alpha /*= 4*/, float stretch /*= 0.5*/, int samples/*= 1*/)
 {
-	nvCheck((w & 1) == 0);
-	
-	const float half_width = float(w / 2);
+	const float half_width = float(w) / 2;
 	const float offset = -half_width;
-	const float nudge = 0.5f;
+
+	const float s_scale = 1.0f / float(samples);
+	const float x_scale = 1.0f / half_width;
 	
-	for(uint i = 0; i < w; i++) {
-		const float x = (i + offset) + nudge;
-		const float sinc_value = sincf(PI * x * stretch);
-		const float window_value = filter_kaiser(x / half_width, alpha);
-		
-		data[i] = sinc_value * window_value;	// @@ sinc windowed by kaiser
+	for(uint i = 0; i < w; i++)
+	{
+		float sum = 0;
+		for(int s = 0; s < samples; s++)
+		{
+			float x = i + offset + (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;
 	}
 
 	normalize();
@@ -307,9 +320,7 @@ void Kernel1::initKaiser(float alpha, float stretch /*= 1*/)
 /// Init 1D Mitchell filter.
 void Kernel1::initMitchell(float b, float c)
 {
-	nvCheck((w & 1) == 0);
-	
-	const float half_width = float(w / 2);
+	const float half_width = float(w) / 2;
 	const float offset = -half_width;
 	const float nudge = 0.5f;
 	
@@ -496,7 +507,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale)
 	nvCheck(w == 9);
 
 	{
-		float elements[] = {
+		const float elements[] = {
             -1, -2, -3, -4, 0, 4, 3, 2, 1,
             -2, -3, -4, -5, 0, 5, 4, 3, 2,
             -3, -4, -5, -6, 0, 6, 5, 4, 3,
@@ -513,7 +524,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale)
 		}
 	}
 	{
-		float elements[] = {
+		const float elements[] = {
             -1, -2, -3, 0, 3, 2, 1,
             -2, -3, -4, 0, 4, 3, 2,
             -3, -4, -5, 0, 5, 4, 3,
@@ -530,7 +541,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale)
 		}
 	}
 	{
-		float elements[] = {
+		const float elements[] = {
             -1, -2, 0, 2, 1,
             -2, -3, 0, 3, 2,
             -3, -4, 0, 4, 3,
@@ -545,7 +556,7 @@ void Kernel2::initBlendedSobel(const Vector4 & scale)
 		}
 	}
 	{
-		float elements[] = {
+		const float elements[] = {
             -1, 0, 1,
             -2, 0, 2,
             -1, 0, 1,
@@ -560,13 +571,100 @@ void Kernel2::initBlendedSobel(const Vector4 & scale)
 }
 
 
-/*PI_DECLARE_TEST(BesselTest) {
+static float frac(float f)
+{
+	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()
+{
+	delete [] m_data;
+}
 
-	for(int i = 0; i < 8; i++) {
-		nvDebug("bessel0(%i) %f =? %f\n", i, bessel0(i), _bessel0(i));
-		PI_TEST(equalf(bessel0(i), _bessel0(i)));
+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;
+	
+	for (uint j = 0; j < m_length; j++)
+	{
+		const float phase = frac(m_width * 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;
+		}
+		
+		for (uint i = 0; i < m_size; i++)
+		{
+			m_data[j * m_size + i] /= total;
+			//m_data[j * m_size + i] /= samples;
+		}
 	}
+}
 
-	return PiTestUnit::Succeed;
-}*/
+void PolyphaseKernel::initKaiser(float alpha /*= 4.0f*/, float stretch /*= 1.0f*/)
+{
+	
+}
+
+/// Print the kernel for debugging purposes.
+void PolyphaseKernel::debugPrint()
+{
+	for (uint j = 0; j < m_length; j++)
+	{
+		nvDebug("%d: ", j);
+		for (uint i = 0; i < m_size; i++)
+		{
+			nvDebug(" %6.4f", m_data[j * m_size + i]);
+		}
+		nvDebug("\n");
+	}
+}
 

src/nvimage/Filter.h

@@ -25,7 +25,7 @@ namespace nv
 			Kaiser,		// Kaiser-windowed sinc filter
 			Num
 		};
-
+		
 		float (*function)(float x);
 		float support;
 	};
@@ -49,9 +49,9 @@ namespace nv
 			return w;
 		}
 		
-		NVIMAGE_API void initFilter(Filter::Enum filter);
+		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);
+		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();
@@ -85,18 +85,49 @@ namespace nv
 		NVIMAGE_API void initEdgeDetection();
 		NVIMAGE_API void initSobel();
 		NVIMAGE_API void initPrewitt();
-
+		
 		NVIMAGE_API void initBlendedSobel(const Vector4 & scale);
-
+		
 	private:
 		const uint w;
 		float * data;
 	};
 
-
-	// @@ Implement non linear filters:
-	// Kuwahara filter
-	// Median filter
+	/// A 1D polyphase kernel
+	class PolyphaseKernel
+	{
+	public:
+		NVIMAGE_API PolyphaseKernel(float width, uint lineLength);
+		NVIMAGE_API PolyphaseKernel(const PolyphaseKernel & k);
+		NVIMAGE_API ~PolyphaseKernel();
+		
+		float valueAt(uint column, uint x) const {
+			return m_data[column * m_size + x];
+		}
+		
+		float width() const {
+			return m_width;
+		}
+		
+		uint size() const {
+			return m_size;
+		}
+		
+		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();
+		
+	private:
+		const float m_width;
+		const uint m_size;
+		const uint m_length;
+		float * m_data;
+	};
 
 } // nv namespace
 

src/nvimage/FloatImage.cpp

@@ -82,18 +82,13 @@ Image * FloatImage::createImage(uint base_component/*= 0*/, uint num/*= 4*/) con
 	for(uint i = 0; i < size; i++) {
 		
 		uint c;
-		uint8 rgba[4];
+		uint8 rgba[4]= {0, 0, 0, 0xff};
 
 		for(c = 0; c < num; c++) {
 			float f = m_mem[size * (base_component + c) + i];
 			rgba[c] = nv::clamp(int(255.0f * f), 0, 255);
 		}
 
-		// Fill the rest with 0xff000000;
-		for(; c < 4; c++) {
-			rgba[c] = c != 3 ? 0 : 0xff;
-		}
-		
 		img->pixel(i) = Color32(rgba[0], rgba[1], rgba[2], rgba[3]);
 	}
 	
@@ -595,6 +590,62 @@ FloatImage * FloatImage::downSample(const Kernel1 & kernel, uint w, uint h, Wrap
 }
 
 
+/// Downsample applying a 1D kernel separately in each dimension.
+FloatImage * FloatImage::downSample(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();	
+	
+	AutoPtr<FloatImage> tmp_image( new FloatImage() );
+	tmp_image->allocate(m_componentNum, w, m_height);
+	
+	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++)
+	{
+		float * tmp_channel = tmp_image->channel(c);
+		
+		for(uint y = 0; y < m_height; y++) {
+			this->applyKernelHorizontal(&xkernel, xkernel.size(), y, c, wm, tmp_channel + y * w);
+		}
+		
+		float * dst_channel = dst_image->channel(c);
+		
+		for(uint x = 0; x < w; x++) {
+			this->applyKernelVertical(&ykernel, x * xscale, yscale, c, wm, tmp_column.unsecureBuffer());
+			
+			for(uint y = 0; y < h; y++) {
+				dst_channel[y * w + x] = tmp_column[y];
+			}
+		}
+	}
+	
+	//return tmp_image.release();
+	return dst_image.release();
+}
+
+
+
 /// Apply 2D kernel at the given coordinates and return result.
 float FloatImage::applyKernel(const Kernel2 * k, int x, int y, int c, WrapMode wm) const
 {
@@ -669,171 +720,61 @@ float FloatImage::applyKernelHorizontal(const Kernel1 * k, int x, int y, int c,
 }
 
 
-
-#if 0
-
-Vec3d bilinear(double u, double v) const
-{
-	u = mod(u*(W-1),W);
-	v = mod(v*(H-1),H);
-
-	Vec3d v1,v2,v3,v4;
-
-	int x_small	= (int)floor(u);
-	int x_big	  = x_small + 1;
-	int y_small = (int)floor(v);
-	int y_big   = y_small + 1;
-
-	if (x_small < 0)
-		x_small = W-1;
-	else if (x_big >= W)
-		x_big = 0;
-	if (y_small < 0)
-		y_small = H-1;
-	else if (y_big >= H)
-		y_big = 0;
-	
-	double fractional_X = u - x_small;
-	double fractional_Y = v - y_small;
-
-	if (nchan == 3)
-	{
-		v1 = Vec3d(pixel(x_small, y_small)[0], pixel(x_small, y_small)[1], pixel(x_small, y_small)[2]);
-		v2 = Vec3d(pixel(x_big, y_small)[0], pixel(x_big, y_small)[1], pixel(x_big, y_small)[2]);
-		v3 = Vec3d(pixel(x_small, y_big)[0], pixel(x_small, y_big)[1], pixel(x_small, y_big)[2]);
-		v4 = Vec3d(pixel(x_big, y_big)[0], pixel(x_big, y_big)[1], pixel(x_big, y_big)[2]);
-	}
-			
-	Vec3d i1 = lerp(v1, v2, fractional_X);
-	Vec3d i2 = lerp(v3, v4, fractional_X);
-
-	return lerp(i1, i2, fractional_Y);
-}
-
-Vec3d bicubic(double u, double v) const
+/// Apply 1D vertical kernel at the given coordinates and return result.
+void FloatImage::applyKernelVertical(const PolyphaseKernel * k, int x, float yscale, int c, WrapMode wm, float * output) const
 {
-	u = mod(u*(W-1),W);
-	v = mod(v*(H-1),H);
-
-	int x_small1	= (int)floor(u),
-		x_small2	= x_small1 - 1,
-		x_big1		= x_small1 + 1,				
-		x_big2		= x_small1 + 2;
-
-	int y_small1	= (int)floor(v),
-		y_small2	= y_small1 - 1,
-		y_big1		= y_small1 + 1,
-		y_big2		= y_small1 + 2;
-
-	x_small1	= (int)mod(x_small1,W);
-	x_small2	= (int)mod(x_small2,W);
-	x_big1		= (int)mod(x_big1,W);
-	x_big2		= (int)mod(x_big2,W);
+	nvDebugCheck(k != NULL);
 	
-	y_small1	= (int)mod(y_small1,H);
-	y_small2	= (int)mod(y_small2,H);
-	y_big1		= (int)mod(y_big1,H);
-	y_big2		= (int)mod(y_big2,H);
+	const float kernelWidth = k->width();
+	const float kernelOffset = kernelWidth / 2;
+	const int kernelLength = k->length();
+	const int kernelWindow = k->size();
 	
-	double fractional_X = u - x_small1;
-	double fractional_Y = v - y_small1;
+	const float * channel = this->channel(c);
 
-	if (nchan == 3)
+	for (int y = 0; y < kernelLength; y++)
 	{
-		// the interpolations across the rows
-		Vec3d row1 = cubic(Vec3d(pixel(x_small2, y_small2)[0], pixel(x_small2, y_small2)[1], pixel(x_small2, y_small2)[2]),
-							Vec3d(pixel(x_small1, y_small2)[0], pixel(x_small1, y_small2)[1], pixel(x_small1, y_small2)[2]),
-							Vec3d(pixel(x_big1, y_small2)[0], pixel(x_big1, y_small2)[1], pixel(x_big1, y_small2)[2]),
-							Vec3d(pixel(x_big2, y_small2)[0], pixel(x_big2, y_small2)[1], pixel(x_big2, y_small2)[2]),
-							fractional_X);
-
-		Vec3d row2 = cubic(Vec3d(pixel(x_small2, y_small1)[0], pixel(x_small2, y_small1)[1], pixel(x_small2, y_small1)[2]),
-							Vec3d(pixel(x_small1, y_small1)[0], pixel(x_small1, y_small1)[1], pixel(x_small1, y_small1)[2]),
-							Vec3d(pixel(x_big1, y_small1)[0], pixel(x_big1, y_small1)[1], pixel(x_big1, y_small1)[2]),
-							Vec3d(pixel(x_big2, y_small1)[0], pixel(x_big2, y_small1)[1], pixel(x_big2, y_small1)[2]),
-							fractional_X);
-
-		Vec3d row3 = cubic(Vec3d(pixel(x_small2, y_big1)[0], pixel(x_small2, y_big1)[1], pixel(x_small2, y_big1)[2]),
-							Vec3d(pixel(x_small1, y_big1)[0], pixel(x_small1, y_big1)[1], pixel(x_small1, y_big1)[2]),
-							Vec3d(pixel(x_big1, y_big1)[0], pixel(x_big1, y_big1)[1], pixel(x_big1, y_big1)[2]),
-							Vec3d(pixel(x_big2, y_big1)[0], pixel(x_big2, y_big1)[1], pixel(x_big2, y_big1)[2]),
-							fractional_X);
-
-		Vec3d row4 = cubic(Vec3d(pixel(x_small2, y_big2)[0], pixel(x_small2, y_big2)[1], pixel(x_small2, y_big2)[2]),
-							Vec3d(pixel(x_small1, y_big2)[0], pixel(x_small1, y_big2)[1], pixel(x_small1, y_big2)[2]),
-							Vec3d(pixel(x_big1, y_big2)[0], pixel(x_big1, y_big2)[1], pixel(x_big1, y_big2)[2]),
-							Vec3d(pixel(x_big2, y_big2)[0], pixel(x_big2, y_big2)[1], pixel(x_big2, y_big2)[2]),
-							fractional_X);
-
-		// now interpolate across the interpolated rows (the columns)
-
-		return cubic(row1,row2,row3,row4,fractional_Y);
+		const float kernelPhase = (kernelWidth * y) - floor(kernelWidth * y);
+		
+		float sum = 0.0f;
+		for (int i = 0; i < kernelWindow; i++)
+		{
+			const int src_y = int(y * kernelWidth - kernelOffset - kernelPhase) + i;
+			const int idx = this->index(x, src_y, wm);
+			
+			sum += k->valueAt(y, i) * channel[idx];
+		}
+		
+		output[y] = sum;
 	}
-	else
-		return Vec3d(0.0);
 }
 
-Vec3d bicubic2(double u, double v) const
+/// Apply 1D horizontal kernel at the given coordinates and return result.
+void FloatImage::applyKernelHorizontal(const PolyphaseKernel * k, float xscale, int y, int c, WrapMode wm, float * output) const
 {
-	u = mod(u*(W-1),W);
-	v = mod(v*(H-1),H);
-
-	int x_small1	= floorf(u),
-		x_small2	= x_small1 - 1,
-		x_big1		= int(x_small1 + 1),
-		x_big2		= int(x_small1 + 2);
-
-	int y_small1	= floorf(v),
-		y_small2	= y_small1 - 1,
-		y_big1		= y_small1 + 1,
-		y_big2		= y_small1 + 2;
-
-	x_small1	= (int)mod(x_small1,W);
-	x_small2	= (int)mod(x_small2,W);
-	x_big1		= (int)mod(x_big1,W);
-	x_big2		= (int)mod(x_big2,W);
+	nvDebugCheck(k != NULL);
 	
-	y_small1	= (int)mod(y_small1,H);
-	y_small2	= (int)mod(y_small2,H);
-	y_big1		= (int)mod(y_big1,H);
-	y_big2		= (int)mod(y_big2,H);
+	const float kernelWidth = k->width();
+	const float kernelOffset = kernelWidth / 2;
+	const int kernelLength = k->length();
+	const int kernelWindow = k->size();
 	
-	double fractional_X = u - x_small1;
-	double fractional_Y = v - y_small1;
+	const float * channel = this->channel(c);
 
-	if (nchan == 3)
+	for (int x = 0; x < kernelLength; x++)
 	{
-		// the interpolations across the rows
-		Vec3d row1 = cubic2(Vec3d(pixel(x_small2, y_small2)[0], pixel(x_small2, y_small2)[1], pixel(x_small2, y_small2)[2]),
-							Vec3d(pixel(x_small1, y_small2)[0], pixel(x_small1, y_small2)[1], pixel(x_small1, y_small2)[2]),
-							Vec3d(pixel(x_big1, y_small2)[0], pixel(x_big1, y_small2)[1], pixel(x_big1, y_small2)[2]),
-							Vec3d(pixel(x_big2, y_small2)[0], pixel(x_big2, y_small2)[1], pixel(x_big2, y_small2)[2]),
-							fractional_X);
-
-		Vec3d row2 = cubic2(Vec3d(pixel(x_small2, y_small1)[0], pixel(x_small2, y_small1)[1], pixel(x_small2, y_small1)[2]),
-							Vec3d(pixel(x_small1, y_small1)[0], pixel(x_small1, y_small1)[1], pixel(x_small1, y_small1)[2]),
-							Vec3d(pixel(x_big1, y_small1)[0], pixel(x_big1, y_small1)[1], pixel(x_big1, y_small1)[2]),
-							Vec3d(pixel(x_big2, y_small1)[0], pixel(x_big2, y_small1)[1], pixel(x_big2, y_small1)[2]),
-							fractional_X);
-
-		Vec3d row3 = cubic2(Vec3d(pixel(x_small2, y_big1)[0], pixel(x_small2, y_big1)[1], pixel(x_small2, y_big1)[2]),
-							Vec3d(pixel(x_small1, y_big1)[0], pixel(x_small1, y_big1)[1], pixel(x_small1, y_big1)[2]),
-							Vec3d(pixel(x_big1, y_big1)[0], pixel(x_big1, y_big1)[1], pixel(x_big1, y_big1)[2]),
-							Vec3d(pixel(x_big2, y_big1)[0], pixel(x_big2, y_big1)[1], pixel(x_big2, y_big1)[2]),
-							fractional_X);
-
-		Vec3d row4 = cubic2(Vec3d(pixel(x_small2, y_big2)[0], pixel(x_small2, y_big2)[1], pixel(x_small2, y_big2)[2]),
-							Vec3d(pixel(x_small1, y_big2)[0], pixel(x_small1, y_big2)[1], pixel(x_small1, y_big2)[2]),
-							Vec3d(pixel(x_big1, y_big2)[0], pixel(x_big1, y_big2)[1], pixel(x_big1, y_big2)[2]),
-							Vec3d(pixel(x_big2, y_big2)[0], pixel(x_big2, y_big2)[1], pixel(x_big2, y_big2)[2]),
-							fractional_X);
-
-		// now interpolate across the interpolated rows (the columns)
-
-		return cubic2(row1,row2,row3,row4,fractional_Y);
+		const float kernelPhase = (kernelWidth * x) - floor(kernelWidth * x);
+		
+		float sum = 0.0f;
+		for (int e = 0; e < kernelWindow; e++)
+		{
+			const int src_x = int(x * kernelWidth - kernelOffset - kernelPhase) + e;
+			const int idx = this->index(src_x, y, wm);
+			
+			sum += k->valueAt(x, e) * channel[idx];
+		}
+		
+		output[x] = sum;
 	}
-	else
-		return Vec3d(0.0);
 }
 
-#endif

src/nvimage/FloatImage.h

@@ -12,7 +12,7 @@ namespace nv
 class Image;
 class Kernel1;
 class Kernel2;
-
+class PolyphaseKernel;
 
 /// Multicomponent floating point image class.
 class FloatImage
@@ -63,11 +63,17 @@ 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 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, int x, float yscale, int c, WrapMode wm, float * output) const;
+	NVIMAGE_API void applyKernelHorizontal(const PolyphaseKernel * k, float xscale, int y, int c, WrapMode wm, float * output) const;
+	
 	
 	uint width() const { return m_width; }
 	uint height() const { return m_height; }