gistfile1.cpp

#include "Rcpp.h"

using namespace std;
using namespace Rcpp;

NumericMatrix lowPass6xDecX(NumericMatrix src)
{
	const int ws = src.ncol(), h = src.nrow();
	const int h2 = h * 2, h3 = h * 3, h4 = h * 4, h5 = h * 5;
	int wr = ws / 2;
	if (wr == 0) wr = 1;

	NumericMatrix result(wr,h);
	NumericMatrix::iterator rptr = result.begin();
	NumericMatrix::const_iterator sptr = src.begin();

	if (ws <= 1)
		result = src;
	else if (ws == 2)
		for (int y = 0; y < h; ++y)
		{
			// use kernel [1 1] / 2
			*rptr++ = (sptr[0] + sptr[h]) / 2.0;
			++sptr;
		}
	else if (ws == 3)
		for (int y = 0; y < h; ++y)
		{
			// use kernel [1 2 1] / 4
			*rptr++ = (sptr[0] + sptr[h] * 2.0 + sptr[h2]) / 4.0;
			++sptr;
		}
	else // general case for ws >= 4
	{
		// left most point - use kernel [10 10 5 1] / 26
		for (int y = 0; y < h; ++y)
		{
			*rptr++ = ((sptr[0] + sptr[h]) * 10.0 +
					sptr[h2] * 5.0 + sptr[h3]) / 26.0;
			++sptr;
		}
		sptr -= h;

		// general case
		int x;
		for (x = 0; x < (ws - 5); x += 2)
		{
			for (int y = 0; y < h; ++y)
			{
				// use kernel [1 5 10 10 5 1] / 32
				*rptr++ = ((sptr[h]  + sptr[h4])  *  5.0 +
						(sptr[h2] + sptr[h3])  * 10.0 +
						(sptr[0]  + sptr[h5])) / 32.0;
				++sptr;
			}
			sptr += h;
		}

		// find out how to treat the right most point
		if (x == (ws - 5))
			for (int y = 0; y < h; ++y)
			{
				// use kernel [1 5 10 10 5] / 31
				*rptr++ = ((sptr[h]  + sptr[h4])  *  5.0 +
						(sptr[h2] + sptr[h3])  * 10.0 +
						sptr[0]) / 31.0;
				++sptr;
			}
		else
			for (int y = 0; y < h; ++y)
			{
				// use kernel [1 5 10 10] / 26
				*rptr++ = ( sptr[0]  + sptr[h]   * 5.0 +
						(sptr[h2] + sptr[h3]) * 10.0) / 26.0;
				++sptr;
			}
	}
	return result;
}

NumericMatrix lowPass6yDecY(NumericMatrix src)
{
	const int w = src.ncol(), hs = src.nrow();
	int hr = hs / 2;
	if (hr == 0) hr = 1;

	NumericMatrix result(w,hr);
	NumericMatrix::iterator rptr = result.begin();
	NumericMatrix::const_iterator sptr = src.begin();

	if (hs <= 1)
		result = src;
	else if (hs == 2)
		for (int x = 0; x < w; ++x)
		{
			// use kernel [1 1]^T / 2
			*rptr++ = (sptr[0] + sptr[1]) / 2.0;
			sptr += 2;
		}
	else if (hs == 3)
		for (int x = 0; x < w; ++x)
		{
			// use kernel [1 2 1]^T / 4
			*rptr++ = (sptr[0] + sptr[1] * 2.0 + sptr[2]) / 4.0;
			sptr += 3;
		}
	else // general case with hs >= 4
		for (int x = 0; x < w; ++x)
		{
			// top most point - use kernel [10 10 5 1]^T / 26
			*rptr++ = ((sptr[0] + sptr[1]) * 10.0 +
					sptr[2] * 5.0 + sptr[3]) / 26.0;
			//++sptr;

			// general case
			int y;
			for (y = 0; y < (hs - 5); y += 2)
			{
				// use kernel [1 5 10 10 5 1]^T / 32
				*rptr++ = ((sptr[1] + sptr[4])  *  5.0 +
						(sptr[2] + sptr[3])  * 10.0 +
						(sptr[0] + sptr[5])) / 32.0;
				sptr += 2;
			}

			// find out how to treat the bottom most point
			if (y == (hs - 5))
			{
				// use kernel [1 5 10 10 5]^T / 31
				*rptr++ = ((sptr[1] + sptr[4])  *  5.0 +
						(sptr[2] + sptr[3])  * 10.0 +
						sptr[0])            / 31.0;
				sptr += 5;
			}
			else
			{
				// use kernel [1 5 10 10]^T / 26
				*rptr++ = ( sptr[0] + sptr[1]  *  5.0 +
						(sptr[2] + sptr[3]) * 10.0) / 26.0;
				sptr += 4;
			}
		}

	return result;
}

//' Gaussian Subsample
//' 
//' Convolves image with a 6x6 separable Gaussian kernel 
//' and subsamples by a factor of two all in one integrated operation. 
//' 
//' @param img a numeric matrix
//' 
//' @return a numeric matrix
//' 
// [[Rcpp::export("GaussianSubsample")]]
NumericMatrix GaussianSubsample(NumericMatrix img)
{
	return(lowPass6yDecY(lowPass6xDecX(img)));
}