#include <math.h>
#include <stdio.h>
#include <omp.h>

#include "project_volumes.h"
#include "project_screen.h"
#include "math_vector3D.h"

Cvolume::Cvolume(int width, int height, int depth)
{
	this->width = width;
	this->height = height;
	this->depth = depth;

	this->iso_depth = default_iso_depth;

	this->data = new unsigned char[width * height * depth];
	this->pal = new unsigned int[256];
}

Cvolume::~Cvolume()
{
}

int Cvolume::GetWidth()
{
	return width;
}

int Cvolume::GetHeight()
{
	return height;
}

int Cvolume::GetDepth()
{
	return depth;
}

int Cvolume::GetIsoDepth()
{
	return iso_depth;
}

void Cvolume::SetIsoDepth(int iso_depth)
{
	this->iso_depth = iso_depth;
}

unsigned char* Cvolume::GetDataPointer()
{
	return data;
}

unsigned int* Cvolume::GetPalPointer()
{
	return pal;
}

unsigned char Cvolume::GetSample(Cvector pos)
{
	float l = (1.0f / width) * 1.025f;

	if (pos.x < l) pos.x = l;
	if (pos.x > 1.0f-l) pos.x = 1.0f-l;
	if (pos.y < l) pos.y = l;
	if (pos.y > 1.0f-l) pos.y = 1.0f-l;
	if (pos.z < l) pos.z = l;
	if (pos.z > 1.0f-l) pos.z = 1.0f-l;

	int ipx = (int)(pos.x * (float)width);
	int ipy = (int)(pos.y * (float)height);
	int ipz = (int)(pos.z * (float)depth);

	int v000 = ipx + ipy * width + ipz * width * height;

	return data[v000];
}

unsigned char Cvolume::GetSampleFiltered(Cvector pos)
{
	float l = (1.0f / width) * 1.025f;

	if (pos.x < l) pos.x = l;
	if (pos.x > 1.0f - l) pos.x = 1.0f - l;
	if (pos.y < l) pos.y = l;
	if (pos.y > 1.0f - l) pos.y = 1.0f - l;
	if (pos.z < l) pos.z = l;
	if (pos.z > 1.0f - l) pos.z = 1.0f - l;

	float vpx = pos.x * (float)width;
	float vpy = pos.y * (float)height;
	float vpz = pos.z * (float)depth;

	int ipx = (int)vpx;
	int ipy = (int)vpy;
	int ipz = (int)vpz;

	int v000 = ipx + ipy * width + ipz * width * height;
	int v001 = v000 + 1;
	int v010 = v000 + width;
	int v011 = v010 + 1;
	int v100 = v000 + width * height;
	int v101 = v100 + 1;
	int v110 = v100 + width;
	int v111 = v110 + 1;

	float pdx = vpx - (float)ipx;
	float pdy = vpy - (float)ipy;
	float pdz = vpz - (float)ipz;

	float c_00x = (float)((float)data[v000] * (1.0f - pdx) + (float)data[v001] * pdx);
	float c_01x = (float)((float)data[v010] * (1.0f - pdx) + (float)data[v011] * pdx);
	float c_0yx = (float)(c_00x * (1.0f - pdy) + c_01x * pdy);

	float c_10x = (float)((float)data[v100] * (1.0f - pdx) + (float)data[v101] * pdx);
	float c_11x = (float)((float)data[v110] * (1.0f - pdx) + (float)data[v111] * pdx);
	float c_1yx = (float)(c_10x * (1.0f - pdy) + c_11x * pdy);

	unsigned char c = (unsigned char)(c_0yx * (1.0f - pdz) + c_1yx * pdz);

	return c;
}

// Volume Generation Algorithms
// ----------------------------

#define fsin_size 65536
bool pvolume_init = false;
int fsin1[fsin_size], fsin2[fsin_size], fsin3[fsin_size];

void init_pvolume()
{
	for (int i=0; i<fsin_size; i++)
	{
		fsin1[i] = (int)(sin(i / 45.0) * 123.0);
		fsin2[i] = (int)(sin(i / 63.0) * 97.0);
		fsin3[i] = (int)(sin(i / 200.0) * 248.0);
	}

	pvolume_init = true;
}

int EvalSines(int x, int y, int z, int octave)
{
	return ((fsin1[x << octave] >> octave) + (fsin2[y << octave] >> octave) + (fsin3[z << octave] >> octave) + (fsin3[(x + y + 64) << octave] >> octave));
}

void Cvolume::ClearBorders()
{
	int wh = width * height;

	for (int y=0; y<height; y++)
	{
		for (int x=0; x<width; x++)
		{
			data[y * width + x] = 0;
			data[(depth-1) * wh + y * width + x] = 0;
		}
	}

	for (int z=0; z<depth; z++)
	{
		for (int y=0; y<height; y++)
		{
			data[z * wh + y * width] = 0;
			data[z * wh + y * width + (width-1)] = 0;
		}
	}

	for (int z=0; z<depth; z++)
	{
		for (int x=0; x<width; x++)
		{
			data[z * wh + x] = 0;
			data[z * wh + (height-1) * width + x] = 0;
		}
	}
}

void Cvolume::GenVolumePlasma()
{
	if (!pvolume_init) init_pvolume();

	int i = 0;
	for (int z=0; z<depth; z++)
	{
		float zz = (float)(depth - z) / (float)depth;
		for (int y=0; y<height; y++)
		{
			for (int x=0; x<width; x++)
			{
				int ii = 0;
				for (int j=1; j<=6; j++)
				{
					ii +=  EvalSines(x,y,z,j);
				}
				int cc = ii & 511;

				if (x==0 || x==width-1 || y==0 || y==height-1 || z==0 || z==depth-1) cc = 0;
				if (cc > 255) cc = 511 - cc;
				int aa = 255;
				if (cc < 192) aa = 0;

				data[i++] = cc;
			}
		}
		printf("%d / %d\r",z,depth);
	}
}

void Cvolume::GenVolumeSphere()
{
	printf("\n");

	int i = 0;
	for (int z=0; z<depth; z++)
	{
		float z_center = z - depth / 2.0f;
		for (int y=0; y<height; y++)
		{
			float y_center = y - height / 2.0f;
			for (int x=0; x<width; x++)
			{
				float x_center = x - width / 2.0f;
				float c = 2.0f * (width / 2.0f - sqrt(x_center * x_center + y_center * y_center + z_center * z_center));
				if (c < 0) c = 0;
				if (c > width) c = (float)width;

				float cc = c / (float)(width);

				unsigned char ci = (unsigned char)(255 * cc);

				data[i++] = ci;
			}
		}
		printf("%d / %d\r",z,depth);
	}
}

void Cvolume::FadeByDistanceFromCenter(float radius)
{
	printf("\n");

	int i = 0;
	for (int z=0; z<depth; z++)
	{
		float z_center = z - depth / 2.0f;
		for (int y=0; y<height; y++)
		{
			float y_center = y - height / 2.0f;
			for (int x=0; x<width; x++)
			{
				float x_center = x - width / 2.0f;
				float c = radius * (width / 2.0f - sqrt(x_center * x_center + y_center * y_center + z_center * z_center));
				if (c < 0) c = 0;
				if (c > width) c = (float)width;

				float cc = c / (float)width;

				unsigned char ccc = (unsigned char)((float)data[i] * cc);

				data[i++] = ccc;
			}
		}
		printf("%d / %d\r",z,depth);
	}
}

void Cvolume::ZeroMemory()
{
	for (int i=0; i<width * height * depth; i++)
	{
		data[i] = 0;
	}
}

void Cvolume::GenVolume(int way)
{
	switch(way)
	{
		case VGEN_EMPTY:
			ZeroMemory();
		break;

		case VGEN_PLASMA:
			GenVolumePlasma();
			FadeByDistanceFromCenter(16);
		break;

		case VGEN_SPHERE:
			GenVolumeSphere();
		break;

		default:
			ZeroMemory();
		break;
	}
}

// Volume Loading Functions
// ------------------------

void Cvolume::LoadVolume(char *filename)
{
	// load file

	unsigned short* fdata;
	FILE *fvol;
	fvol = fopen(filename,"rt");
	int count = 0;

	fseek(fvol, 0, SEEK_END);
    count = ftell(fvol) / 2;
    rewind(fvol);

	fdata = new unsigned short[count];
	fread(fdata, sizeof(char), count, fvol);
	fclose(fvol);


	// file loaded

	for (int i=0; i<count; i++)
	{
		int c = fdata[i] >> 6;
		int a = c >> 8;
		//if (a < 128) a = 0;
		unsigned char cc = c;
		data[i] = cc;
	}

	delete [] fdata;
}

// Volume Rendering Functions
// --------------------------

bool bbox_check(Cvector origin, Cvector dir, float &far_d, float &near_d)
{
	float left = 0.0f;
	float right = 1.0f;
	float down = 0.0f;
	float up = 1.0f;
	float back = 0.0f;
	float front = 1.0f;

	float tNearX, tFarX;
	float tNearY, tFarY;
	float tNearZ, tFarZ;

	// Slab X
	if (dir.x >= 0.0f)
	{
		tNearX = (left - origin.x) / dir.x;
		tFarX = (right - origin.x) / dir.x;
	}
	else if (dir.x < 0.0f)
	{
		tNearX = (right - origin.x) / dir.x;
		tFarX = (left - origin.x) / dir.x;
	}

	// Slab Y
	if (dir.y >= 0.0f)
	{
		tNearY = (down - origin.y) / dir.y;
		tFarY = (up - origin.y) / dir.y;
	}
	else if (dir.y < 0.0f)
	{
		tNearY = (up - origin.y) / dir.y;
		tFarY = (down - origin.y) / dir.y;
	}

	// Slab Z
	if (dir.z >= 0.0f)
	{
		tNearZ = (back - origin.z) / dir.z;
		tFarZ = (front - origin.z) / dir.z;
	}
	else if (dir.z < 0.0f)
	{
		tNearZ = (front - origin.z) / dir.z;
		tFarZ = (back - origin.z) / dir.z;
	}

	// We have to find now the the largest of the Near and the the smallest of the Far

	float tNearMax, tFarMin;

	if (tNearX > tNearY)
		tNearMax = tNearX;
	else
		tNearMax = tNearY;
	if (tNearZ > tNearMax)
		tNearMax = tNearZ;

	if (tFarX < tFarY)
		tFarMin = tFarX;
	else
		tFarMin = tFarY;
	if (tFarZ < tFarMin)
		tFarMin = tFarZ;

	far_d = tFarMin;
	near_d = tNearMax;

	// Final test
	return (tNearMax < tFarMin);
}

unsigned char Cvolume::back_step_sample(Cvector pos, Cvector dir, int num)
{
	unsigned char c;
	for (int i=0; i<num; i++)
	{
		pos = pos + dir;
		c = GetSampleFiltered(pos);

		if (c <= iso_depth)
		{
			c = GetSampleFiltered(pos - dir);
			break;
		}
	}
	return c;
}

void Cvolume::RenderVolume(Cvector pos, Cvector dir, Cvector up, Cvector right)
{
	unsigned int *vram = (unsigned int*)screen->pixels;

	const float jump_size = 4.0f;

	float small_step = 1.0f / this->GetWidth();
	float step = small_step * jump_size;

	float xtd = 1.0f / (float)ScreenWidth;
	float ytd = 1.0f / (float)ScreenHeight;

	#pragma omp parallel for num_threads(4)
	for (int y=0; y<ScreenHeight; y++)
	{
		for (int x=0; x<ScreenWidth; x++)
		{
			float yt = -0.5f + ytd * y;
			float xt = -0.5f + xtd * x;
			Cvector pscreen = right * xt + up * yt + dir;
			Cvector vdir = pscreen.normalized();
			Cvector vdirstep = vdir * step;

			Cvector vpos = pos + Cvector(0.5f, 0.5f, 0.5f);

			float far_d = 0.0f;
			float near_d = 0.0f;
			bool bchk = bbox_check(vpos, vdirstep, far_d, near_d);

			int col = 0;

			if (bchk)
			{
				if (near_d > 0.0)
				{
					Cvector mdir = vdirstep * near_d;
					vpos = vpos + mdir;
				}
				else near_d = 0;

				float f_num_steps = far_d - near_d;

				float ns = 0.0;
				do	
				{
					unsigned char c = GetSampleFiltered(vpos);

					if (c > iso_depth)
					{
						Cvector small_backstep = vdir.negate() * small_step;
						c = back_step_sample(vpos, small_backstep, int(jump_size));
						col = (int)(vpos.z * (float)c);
						break;
					}
					vpos = vpos + vdirstep;

				}while(ns++ < f_num_steps);
			}

			//*(vram + x + y * ScreenWidth) = (col << 16) | (col << 8) | col;
			*(vram + x + y * ScreenWidth) = pal[col];
		}
	}
}