#include <vector>
using namespace std;

#include "rtsPoint3d.h"
#include "rts_itkVolume.h"
typedef vector<point3D<double>> filament;

#define BACKGROUND_NODE	0
#define SKELETON_NODE	255
#define ENDPOINT_NODE	128
#define VISITED_NODE	64

class rtsSkeleton
{
	rts_itkVolume<unsigned char> implicit_skeleton;
	vector<filament> explicit_filaments;
	vector<point3D<unsigned int>> explicit_endpoints;

	void ComputeEndPoints();

public:
	void ConstructExplicitSkeleton();
	void SmoothExplicitSkeleton(int iterations);
	void SaveExplicitSkeleton(const char* filename);
}

void rtsSkeleton::ComputeEndPoints()
{
	/*This function finds all of the filament end points in the implicit skeleton.
	The end points are stored in the explicit_endpoints vector.
	An end point is defined as a node with (n == 1) or (n > 2) connected components.
	*/

	unsigned int max_x = implicit_skeleton.DimX();
	unsigned int max_y = implicit_skeleton.DimY();
	unsigned int max_z = implicit_skeleton.DimZ();
	unsigned int conn;

	indextype x, y, z;
	for(x=0; x<max_x; x++)		//go through each voxel in the data set
		for(y=0; y<max_y; y++)
			for(z=0; z<max_z; z++)
			{
				if(implicit_skeleton.GetPixel(x, y, z) == SKELETON_NODE)	//if the voxel is part of the skeleton
				{
					conn = Neighbors26(&implicit_skeleton, x, y, z, 1);	//test the number of connected components
					if(conn == 1 || conn > 2)	//if the number of connected components meets the
												//criteria, store the node in the explicit_endpoints vector.
					{
						explicit_endpoints.push_back(point3D<unsigned int>(x, y, z));
					}
				}

			}

	//tag the endpoints in the implicit function
	vector<point3D<unsigned int>>::iterator i;
	for(i=explicit_endpoints.begin(); i!=explicit_endpoints.end(); i++)
		implicit_skeleton.xyz(i->x, i->y, i->z) = ENDPOINT_NODE;

			//cout<<"Number of end points:  "<<explicit_endpoints.size()<<endl;
}

void rtsSkeleton::ConstructExplicitSkeleton()
{
	/*The implicit skeleton is constructed as follows:
	1)  Label all branch nodes.
	2)  For each branch, find the number of fibers connected to that branch.
		a)  Create each fiber as a seperate structure and store them in a queue.
		b)  Iteratively track each fiber.
	*/
	ComputeEndPoints();	//compute the endpoints for the explicit skeleton and
						//label them in the implicit skeleton
	
	//iterate through each end point
	vector<point3D<unsigned int>>::iterator i;
	rtsFunction3D<unsigned char> local;
	unsigned int x, y, z;
	for(i = explicit_endpoints.begin(); i != explicit_endpoints.end(); i++)
	{
		//cout<<"Tracing Filament..."<<endl;
		//i->print();
		local = getLocalRegion(i->x, i->y, i->z);
		local(1, 1, 1) = 0;
		//local.toConsole();
		for(x=0; x<3; x++)
			for(y=0; y<3; y++)
				for(z=0; z<3; z++)
					if(local.xyz(x, y, z) > VISITED_NODE)
					{
						//cout<<(int)local.xyz(x, y, z)<<endl;
						//cout<<(int)implicit_skeleton.xyz(i->x - 1 + x,
						//							  i->y - 1 + y,
						//							  i->z - 1 + z)<<endl;
						TraceFilament((*i), point3D<unsigned int>(i->x - 1 + x,
																  i->y - 1 + y,
																  i->z - 1 + z));
					}
		//since we have tracked all filaments attached to an endpoint
		//set the endpoint as visited
		implicit_skeleton.xyz(i->x, i->y, i->z) = VISITED_NODE;
	}

	//cout<<"Tracing complete"<<endl;


}