#include "rtsImplicit3D.h"
#include <iostream>
using namespace std;

#define DIST_MAX			255
#define DIST_UNSIGNED		0
#define DIST_SIGNED			1

float ComputeSurfaceDistance(point3D<float> p0, point3D<float> p1, float val0, float val1, float s)
{
	/*This function computes the distance from p0 to the surface, given two points p0 and p1
	on either side of the surface (with values v0 and v1 respectively).  surface specifies
	the value at the surface.
	*/

	//compute the normalized position of the surface between p0 and p1
	float s_norm_pos = (s - val0) / (val1 - val0);
	//compute the actual position of the surface
	point3D<float> s_pos = p0 + s_norm_pos * (p1 - p0);
	//compute the distance from p0 to the surface
	float result = (s_pos - p0).Length();
	//cout<<"distance: "<<result<<endl;
	return result;
}

void CreateBoundaryConditions(rtsImplicit3D<unsigned char>* function, float threshold,
							  rtsImplicit3D<float>* &result, rtsImplicit3D<bool>* &mask)
{
	/*This function creates an initial signed distance function from a threshold image.
	All voxels adjacent to the surface specified by the threshold are initialized with a
	distance value.  Low values are inside, high values are outside.
	*/
	//current and neighboring voxel flags (false = inside, true = outside)
	bool c, x_p, x_n, y_p, y_n, z_p, z_n;
	float d_xp, d_xn, d_yp, d_yn, d_zp, d_zn;
	float in_out = 1;

	//boundary condition function and the mask
	result = new rtsImplicit3D<float>(function->DimX(), function->DimY(), function->DimZ());
	//get the parameterization
	point3D<float> min_domain= function->getMinDomain();
	point3D<float> max_domain = function->getMaxDomain();
	result->Parameterize(min_domain.x, max_domain.x, min_domain.y, max_domain.y, min_domain.z, max_domain.z);
	(*result) = DIST_MAX;
	//create a mask
	mask = new rtsImplicit3D<bool>(function->DimX(), function->DimY(), function->DimZ());
	(*mask) = false;

	cout<<"done making boundary condition function"<<endl;
	//for each voxel
	vector3D<int> size(function->DimX(), function->DimY(), function->DimZ());	//get the function size
	int x, y, z;
	for(x=0; x<size.x; x++)
		for(y=0; y<size.y; y++)
			for(z=0; z<size.z; z++)
			{
				//reset flags
				c=x_p=x_n=y_p=y_n=z_p=z_n=true;
				in_out = 1.0;
				//look at the current voxel
				if((*function)(x, y, z) < threshold)
					c=false;
				else c=true;
				//look at each neighboring voxel
				if(x-1 < 0) x_n = c;	//X
				else if((*function)(x-1, y, z) < threshold) x_n = false;
				if(x+1 >= size.x) x_p = c;
				else if((*function)(x+1, y, z) < threshold) x_p = false;
				if(y-1 < 0) y_n = c;	//Y
				else if((*function)(x, y-1, z) < threshold) y_n = false;
				if(y+1 >= size.y) y_p = c;
				else if((*function)(x, y+1, z) < threshold) y_p = false;
				if(z-1 < 0) z_n = c;	//Z
				else if((*function)(x, y, z-1) < threshold) z_n = false;
				if(z+1 >= size.z) z_p = c;
				else if((*function)(x, y, z+1) < threshold) z_p = false;

				//set the distance from the isosurface
				if(c == false)
					in_out = -1.0;
				if(x_n != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
										 ComputeSurfaceDistance((*function).getParameter(x, y, z), 
													  (*function).getParameter(x-1, y, z),
													  (*function)(x, y, z),
													  (*function)(x-1, y, z),
													  threshold) * in_out);
				if(x_p != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											ComputeSurfaceDistance((*function).getParameter(x, y, z), 
													  (*function).getParameter(x+1, y, z),
													  (*function)(x, y, z),
													  (*function)(x+1, y, z),
													  threshold) * in_out);
				if(y_n != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											ComputeSurfaceDistance((*function).getParameter(x, y, z), 
													  (*function).getParameter(x, y-1, z),
													  (*function)(x, y, z),
													  (*function)(x, y-1, z),
													  threshold) * in_out);
				if(y_p != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											ComputeSurfaceDistance((*function).getParameter(x, y, z), 
													  (*function).getParameter(x, y+1, z),
													  (*function)(x, y, z),
													  (*function)(x, y+1, z),
													  threshold) * in_out);
				if(z_n != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											ComputeSurfaceDistance((*function).getParameter(x, y, z-1), 
													  (*function).getParameter(x, y, z),
													  (*function)(x, y, z),
													  (*function)(x, y, z-1),
													  threshold) * in_out);
				if(z_p != c)
					(*result)(x, y, z) = min((*result)(x,y,z),
											ComputeSurfaceDistance((*function).getParameter(x, y, z), 
													  (*function).getParameter(x, y, z+1),
													  (*function)(x, y, z),
													  (*function)(x, y, z+1),
													  threshold) * in_out);

				//set the mask to 1 if the voxel is on an edge node
				if(x_n != c || x_p != c || y_n != c || y_p != c || z_n != c || z_p != c)
					(*mask)(x, y, z) = true;
			}
				

	//if a line between the two voxels crosses the surface
		//find the distance between the voxel center and the surface


			cout<<"done computing boundary conditions"<<endl;
}

double ManhattanDistance(rtsImplicit3D<float>* function,
						 point3D<indextype> point,
						 vector3D<float> voxelsize,
						 int type = DIST_UNSIGNED)
{
	/*This function updates the manhattan distance from a surface using the manhattan
	distance of its neighboring points.
	*/
	indextype x, y, z;
	x=point.x; y=point.y, z=point.z;
	int sign = 1;
	float result = DIST_MAX;
	float near_value;				//the value of the neighbor being considered
	float possible_value;			
	if(x!=0)
	{ 
		near_value = (*function)(x-1, y, z);
		if(type == DIST_UNSIGNED)
			result = min(result, near_value + voxelsize.x);
		else if(type == DIST_SIGNED)
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + voxelsize.x);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}

	}
	if(x!=function->DimX()-1)
	{
		near_value = (*function)(x+1, y, z);
		if(type == DIST_UNSIGNED)
			result = min(result, near_value + voxelsize.x);
		else if(type == DIST_SIGNED)
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + voxelsize.x);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(y!=0)
	{
		near_value = (*function)(x, y-1, z);
		if(type == DIST_UNSIGNED)
			result = min(result, near_value + voxelsize.y);
		else if(type == DIST_SIGNED)
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + voxelsize.y);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(y!=function->DimY()-1)
	{
		near_value = (*function)(x, y+1, z);
		if(type == DIST_UNSIGNED)
			result = min(result, near_value + voxelsize.y);
		else if(type == DIST_SIGNED)
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + voxelsize.y);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(z!=0)
	{
		near_value = (*function)(x, y, z-1);
		if(type == DIST_UNSIGNED)
			result = min(result, near_value + voxelsize.z);
		else if(type == DIST_SIGNED)
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + voxelsize.z);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	if(z!=function->DimZ()-1)
	{
		near_value = (*function)(x, y, z+1);
		if(type == DIST_UNSIGNED)
			result = min(result, near_value + voxelsize.z);
		else if(type == DIST_SIGNED)
		{
			if(near_value<0) sign = -1; else sign = 1;	//determine if the value is inside or outside
			possible_value = sign*(fabs(near_value) + voxelsize.z);
			if(fabs(possible_value) < fabs(result))
				result = possible_value;
		}
	}
	return result;

}
void Eikonal_Manhattan(rtsImplicit3D<float>* &function, 
					   rtsImplicit3D<bool>* mask, 
					   int type = DIST_UNSIGNED)
{
	/*This function estimates the Eikonal equation based on the manhattan distance.  The
	function constantly increases values outwards from the boundary conditions.
	*/

	//compute the distance between two voxels
	vector3D<float> voxel_size;
	voxel_size.x = fabs(function->getParameter(0, 0, 0).x - function->getParameter(1, 0, 0).x);
	voxel_size.y = fabs(function->getParameter(0, 0, 0).y - function->getParameter(0, 1, 0).y);
	voxel_size.z = fabs(function->getParameter(0, 0, 0).z - function->getParameter(0, 0, 1).z);
	
	//use fast sweeping to compute the manhattan distance
	//0:X  0:Y  0:Z
	cout<<"first iteration..."<<endl;
	int x,y,z;
	for(x=0; x<function->DimX(); x++)
		for(y=0; y<function->DimY(); y++)
			for(z=0; z<function->DimZ(); z++)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
	cout<<"second iteration..."<<endl;
	//0:X 0:Y Z:0
	for(x=0; x<function->DimX(); x++)
		for(y=0; y<function->DimY(); y++)
			for(z=function->DimZ()-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
	cout<<"third iteration..."<<endl;
	//0:X Y:0 0:Z
	for(x=0; x<function->DimX(); x++)
		for(y=function->DimY()-1; y>=0; y--)
			for(z=0; z<function->DimZ(); z++)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
	cout<<"fourth iteration..."<<endl;
	//0:X Y:0 Z:0
	for(x=0; x<function->DimX(); x++)
		for(y=function->DimY()-1; y>=0; y--)
			for(z=function->DimZ()-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
	cout<<"fifth iteration..."<<endl;
	//X:0 0:Y 0:Z
	for(x=function->DimX()-1; x>=0; x--)
		for(y=0; y<function->DimY(); y++)
			for(z=0; z<function->DimZ(); z++)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
	cout<<"sixth iteration..."<<endl;
	//X:0 0:Y Z:0
	for(x=function->DimX()-1; x>=0; x--)
		for(y=0; y<function->DimY(); y++)
			for(z=function->DimZ()-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
	cout<<"seventh iteration..."<<endl;
	//X:0 Y:0 0:Z
	for(x=function->DimX()-1; x>=0; x--)
		for(y=function->DimY()-1; y>=0; y--)
			for(z=0; z<function->DimZ(); z++)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
	cout<<"eighth iteration..."<<endl;
	//X:0 Y:0 Z:0
	for(x=function->DimX()-1; x>=0; x--)
		for(y=function->DimY()-1; y>=0; y--)
			for(z=function->DimZ()-1; z>=0; z--)
				//if the current point is not a boundary value
				if(!(*mask)(x, y, z))
					(*function)(x,y,z) = ManhattanDistance(function, point3D<indextype>(x, y, z), voxel_size, type);
	cout<<"done."<<endl;
}