gistfile1.cpp

#include "Rcpp.h"

using namespace std;
using namespace Rcpp;

/*
inline unsigned char saturate( float x )
{
    return x > 255.0f ? 255
         : x < 0.0f   ? 0
         :              unsigned char(x);
}
 */

//' Bicubic Resize
//'
//' @param img a numeric matrix
//'
//' @return a numeric matrix
//'
// [[Rcpp::export("GaussianSubsample")]]
NumericMatrix BicubicResize(NumericMatrix img, int new_width, int new_height)
{
	NumericMatrix out(new_height, new_width, 3);

	const double tx = (double)img.ncol() / (double)new_width;
	const double ty = (double)img.nrow() / (double)new_height;

	//const std::size_t row_stride = dest_width * channels;
	float C[5] = { 0 };

	for (unsigned i = 0; i < new_height; ++i)
	{
		for (unsigned j = 0; j < new_width; ++j)
		{
			const float x = float(tx * j);
			const float y = float(ty * i);
			const float dx = tx * j - x, dx2 = dx * dx, dx3 = dx2 * dx;
			const float dy = ty * i - y, dy2 = dy * dy, dy3 = dy2 * dy;

			for (int k = 0; k < 3; ++k)
			{
				for (int jj = 0; jj < 4; ++jj)
				{
					const int idx = y - 1 + jj;
					float a0 = img[idx, x, k];
					float d0 = img[idx, x - 1, k] - a0;
					float d2 = img[idx, x + 1, k] - a0;
					float d3 = img[idx, x + 2, k] - a0;
					float a1 = -(1.0f / 3.0f) * d0 + d2 - (1.0f / 6.0f) * d3;
					float a2 =          0.5f  * d0 +              0.5f *  d2;
					float a3 = -(1.0f / 6.0f) * d0 - 0.5f * d2 + (1.0f / 6.0f) * d3;
					C[jj] = a0 + a1 * dx + a2 * dx2 + a3 * dx3;

					d0 = C[0] - C[1];
					d2 = C[2] - C[1];
					d3 = C[3] - C[1];
					a0 = C[1];
					a1 = -(1.0f / 3.0f) * d0 + d2 -(1.0f / 6.0f) * d3;
					a2 =          0.5f  * d0 +             0.5f  * d2;
					a3 = -(1.0f / 6.0f) * d0 - 0.5f * d2 + (1.0f / 6.0f) * d3;
					out[i * row_stride + j * channels + k] = saturate( a0 + a1 * dy + a2 * dy2 + a3 * dy3 );
				}
			}
		}
	}

	return out;
}