#include "ToIntel.h"

using namespace std;

static bool SameCuboid(const IpprCuboid& a, const IpprCuboid& b)
{
	return (a.x == b.x && a.y == b.y && a.z == b.z && a.width == b.width && a.height == b.height && a.depth == b.depth);
}

static unsigned short absDiff(unsigned short a, unsigned short b)
{
	int diff = abs(static_cast<int>(a) - static_cast<int>(b));
	return static_cast<unsigned short>(diff);
}

static float Trilinear(float x, float y, float z)
{
	assert(x >= 0.0f && x <= 1.0f && y >= 0.0f && y <= 1.0f && z >= 0.0f && z <= 1.0f);

	return (1.0f - x) * (1.0f - y) * (1.0f - z);
}

//We assume the center of voxel with integer indices (i,j,k) is (i+0.5f,j+0.5f,k+0.5f).
//The extreme corners of the volume have voxel coordinates (0,0,0) and (xdim,ydim,zdim).
static unsigned short Interpolate(
	unsigned short* pVol,
	int xdim,
	int ydim,
	int zdim,
	float xf,
	float yf,
	float zf
)
{
	if (xf < 0.0f || xf > static_cast<float>(xdim) || yf < 0.0f || yf > static_cast<float>(ydim) || zf < 0.0f || zf > static_cast<float>(zdim)) {
		return 0; //Points outside the volume get 0
	}

	const int i = static_cast<int>(floorf(xf - 0.5f));
	const int j = static_cast<int>(floorf(yf - 0.5f));
	const int k = static_cast<int>(floorf(zf - 0.5f));

	//(xmin,ymin,zmin) is center of a voxel
	const float xmin = (i + 0.5f);
	const float ymin = (j + 0.5f);
	const float zmin = (k + 0.5f);

	//(xmax,ymax,zmax) is center of a voxel (if inside the volume)
	const float xmax = xmin + 1.0f;
	const float ymax = ymin + 1.0f;
	const float zmax = zmin + 1.0f;

	const float zz[8] = { zf - zmin, zf - zmin, zf - zmin, zf - zmin, zmax - zf, zmax - zf, zmax - zf, zmax - zf };
	const float yy[8] = { yf - ymin, yf - ymin, ymax - yf, ymax - yf, yf - ymin, yf - ymin, ymax - yf, ymax - yf };
	const float xx[8] = { xf - xmin, xmax - xf, xf - xmin, xmax - xf, xf - xmin, xmax - xf, xf - xmin, xmax - xf };

	const int kk[8] = { k, k,   k,   k,   k + 1, k + 1, k + 1, k + 1 };
	const int jj[8] = { j, j,   j + 1, j + 1, j,   j,   j + 1, j + 1 };
	const int ii[8] = { i, i + 1, i,   i + 1, i,   i + 1, i,   i + 1 };

	float val, coef, sum = 0.0f, prod = 0.0f;

	for (int n = 0; n < 8; n++) {
		if (ii[n] < 0 || ii[n] >= xdim || jj[n] < 0 || jj[n] >= ydim || kk[n] < 0 || kk[n] >= zdim) continue; //Don't use voxels outside the volume

		coef = Trilinear(xx[n], yy[n], zz[n]);

		val = static_cast<float>(pVol[kk[n] * ydim * xdim + jj[n] * xdim + ii[n]]);

		sum += coef;

		prod += coef * val;
	}

	if (fabsf(sum) < FLT_EPSILON) {
		assert(false);
		return 0;
	}

	val = prod / sum;

	int n = static_cast<int>(val + 0.5f);
	n = (n < 0 ? 0 : (n > USHRT_MAX ? USHRT_MAX : n)); //n should already be in range, but for safety

	return static_cast<unsigned short>(n);
}

enum class ValType { Random, Fixed, Zero };

//make a volume of size (dim x dim x dim)
static vector<unsigned short> makeVolume(int dim, ValType type)
{
	vector<unsigned short> vol(dim * dim * dim);

	if (type == ValType::Random) {
		default_random_engine generator;
		generator.seed(12345);
		uniform_int_distribution<unsigned short> distribution(100, 1000);

		for (unsigned short& v : vol) {
			v = distribution(generator);
		}
	}
	else if (type == ValType::Fixed) {
		int i = 0, j, z;
		for (z = 0; z < dim; z++) {
			unsigned short val = 100 * (z + 1);
			for (j = 0; j < dim * dim; j++, i++) vol[i] = val;
		}
	}
	else {
		assert(type == ValType::Zero);
		fill(vol.begin(), vol.end(), 0);
	}

	return vol;
} //end makeVolume


static unsigned short testResize(int dimIn, int dimOut)
{
	vector<unsigned short> src = makeVolume(dimIn, ValType::Fixed);
	vector<unsigned short> dstIntel = makeVolume(dimOut, ValType::Zero);
	vector<unsigned short> dstUS = makeVolume(dimOut, ValType::Zero);

	IpprVolume srcVol = { dimIn, dimIn, dimIn };
	IpprCuboid srcVOI = { 0, 0, 0, dimIn, dimIn, dimIn }, dstVOI = { 0, 0, 0, dimOut, dimOut, dimOut }, dstVOItmp;

	int srcStep = dimIn * sizeof(unsigned short);
	int srcPlaneStep = dimIn * dimIn * sizeof(unsigned short);

	int dstStep = dimOut * sizeof(unsigned short);
	int dstPlaneStep = dimOut * dimOut * sizeof(unsigned short);

	double scale = static_cast<double>(dimOut) / static_cast<double>(dimIn);
	double shift = 0.0;
	
	IppStatus status;

	assert(IPPI_INTER_LINEAR == ippLinear);

	status = ipprGetResizeCuboid(srcVOI, &dstVOItmp, scale, scale, scale, shift, shift, shift, ippLinear);
	assert(SameCuboid(dstVOI, dstVOItmp));

	int bufSize{ 0 };
	ipprResizeGetBufSize(srcVOI, dstVOI, 1, ippLinear, &bufSize);
	bufSize = std::max(bufSize, 10); //To always have positive size

	Ipp8u* pBuffer = ippsMalloc_8u(bufSize);
	assert(pBuffer != nullptr);

	status = ipprResize_16u_C1V(src.data(), srcVol, srcStep, srcPlaneStep, srcVOI, dstIntel.data(), dstStep, dstPlaneStep, dstVOI, scale, scale, scale, shift, shift, shift, ippLinear, pBuffer);
	assert(status == ippStsOk);

	ippsFree(pBuffer);

	//For information only, no voxels should have the value 0
	int numZero = count(dstIntel.begin(), dstIntel.end(), static_cast<unsigned short>(0));

	unsigned short maxDiff = 0; //maximum difference between Intel interpolation and our interpolation
	tuple<int, int, int> maxLocation{ 0,0,0 }; //location where the difference is maximum

	int i = 0, x, y, z;
	for (z = 0; z < dimOut; z++) {
		float zOut = static_cast<float>(z + 0.5f); //Move to center of voxel
		float zIn = static_cast<float>((zOut - shift) / scale);

		for (y = 0; y < dimOut; y++) {
			float yOut = static_cast<float>(y + 0.5f); //Move to center of voxel
			float yIn = static_cast<float>((yOut - shift) / scale);

			for (x = 0; x < dimOut; x++, i++) {
				float xOut = static_cast<float>(x + 0.5f); //Move to center of voxel
				float xIn = static_cast<float>((xOut - shift) / scale);

				dstUS[i] = Interpolate(src.data(), dimIn, dimIn, dimIn, xIn, yIn, zIn);
				assert(dstUS[i] != 0);

				unsigned short valIntel = dstIntel[i];
				unsigned short valUS = dstUS[i];

				unsigned short diff = absDiff(valIntel, valUS);

				if (diff > maxDiff) {
					maxDiff = diff;
					maxLocation = { x,y,z };
				}
			}
		}
	}

	return maxDiff;
} //End testResize

bool testResize()
{
	ippInit();

	const unsigned short maxAllowedDiff = 2; //Could make a percentage error, but here just use absolute error of 2

	const int dimMin = 2;
	const int dimMax = 50;
	const int dimStep = 1;

	int numTrials = 0, numErrors = 0;

	//2 to 50 in both dimIn and dimOut leads to 49*48 = 2352 trials, of which 186 fail
	ofstream f("C:\\ipprResize_failures.txt");
	if (!f.is_open()) return false;
	f << "Dim In, Dim Out, Max Intensity Difference" << endl;

	for (int dimIn = dimMin; dimIn <= dimMax; dimIn += dimStep) {
		for (int dimOut = dimMin; dimOut <= dimMax; dimOut += dimStep) {
			if (dimOut == dimIn) continue;

			//cout << dimIn << " , " << dimOut << endl;

			unsigned short maxDiff = testResize(dimIn, dimOut);

			if (maxDiff > maxAllowedDiff) {
				f << dimIn << " , " << dimOut << " , " << maxDiff << endl;
				numErrors++;
			}

			numTrials++;
		}
	}

	f << "Finished: Num Trials " << numTrials << " Num Errors " << numErrors << endl;
	f.close();

	return (numErrors == 0);
}