codebase / stimlib

#include "objJedi.h"
#include <list>
#include <vector>
#include <limits.h>
#include <fstream>

//itk includes
#include "../../VesselAnalysis/VesselAnalysis/VOL_to_ITK.h"
#include "itkDanielssonDistanceMapImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkExpandImageFilter.h"
#include "itkDiscreteGaussianImageFilter.h"
#include "itkCastImageFilter.h"


/* This code uses the octree library designed by Simon Perreault (http://nomis80.org/code/octree.html)*/
#include "octree/octree.h"

#define OCTREE_SIZE		1024
#define BUFFER_ZONE		0


using namespace std;



//definitions for simple structures
typedef vector<point3D<float>> FilamentType;
typedef point3D<float> CoordType;
typedef vector<FilamentType> NetworkType;

struct EdgeType
{
	unsigned int v0;
	unsigned int v1;
	float avg_radius;
	float min_radius;
	float max_radius;
	float volume;
	float length;
	bool valid;
};

struct NodeType
{
	CoordType p;
	vector<unsigned int> edges;
	bool valid;
};

struct ConnectType				//this structures stores the info necessary to combine two fibers
{			
	unsigned int edge0;
	unsigned int edge1;
	unsigned int node0;
	unsigned int node1;
	list<FilamentType>::iterator fiber0;
	list<FilamentType>::iterator fiber1;
};

struct ConnectableType
{
	FilamentType fiber;
	bool front;
	bool back;
	unsigned int id;
};

struct AttributeType
{
	float TotalFiberLength;
	int NumFibers;
	int NumPoints;
	float TotalVolume;
	int NumDisconnectedFibers;
	int NumBoundaryFibers;
	int NumBifurcations;
};

struct BranchType
{
	float angle;
	unsigned int branch_id;
};


class rtsNetwork
{
private:
	int id;							//used for swc output
	ofstream outfile;				//also for swc
	NetworkType network;			//list of filaments that make up the network
	vector<NodeType> NodeList;		//vector of endpoints
	vector<EdgeType> EdgeList;		//vector of vessels as edges
	vector3D<unsigned int> boundary;
	vector3D<float> position;
	

	float calcFiberLength(unsigned int f);
	float calcFiberBend(unsigned int f);
	float calcTotalLength();
	bool isBoundaryNode(unsigned int node);
	vector3D<float> calcTrajectory(FilamentType f, bool end);
	bool compareTrajectories(CoordType p0, vector3D<float> v0, CoordType p1, vector3D<float> v1, float distance, float cosine_angle);
	void calcGraph();
	vector<unsigned int> getConnections(unsigned int node);
	CoordType toTissue(CoordType grid_coord);
	FilamentType combineFibers(FilamentType f0, bool endpoint0, FilamentType f1, bool endpoint1);
	bool isConnectable(ConnectableType c0, ConnectableType c1, bool& c0_end, bool& c1_end, float distance, float cosine_angle);
	void traverseGraph(int node, int parent_id);
	void outputFiberToSWC(int edge, int node_from, int parent_id);

	vector<BranchType> getFiberAngles(unsigned int fiber);
	void reset();
	
public:
	//attribute variables
	AttributeType Attributes;

	//control variables
	point3D<float> voxel_size;

	//cleaning
	int cleanDegenerateEdges(float length);
	int cleanRedundantEdges();
	int cleanBarbs(float length);

	//methods
	rtsOBJ LoadOBJModel(const char* filename);
	void AddOBJModel(const char* filename);
	void SaveSWCModel(const char *filename, float root_x, float root_y);
	void SaveSWCModel(const char *filename);
	void Clean(float degenerate_length, float barb_length);
	void CalculateAttributes();
	rtsOBJ Repair(float max_distance, float cos_angle);
	rtsOBJ CreateModel();
	void printFiberStats(unsigned int fiber);
	void ReMesh(float resample_length);
	void SaveAngles(const char* filename);
	void SaveBends(const char* filename);
	void SaveLengths(const char* filename);
	void SaveMathematicaGraph(const char* filename);
	void SaveMathematicaNodeDistances(const char* filename);
	void SetVoxelSize(float dx, float dy, float dz){voxel_size.x = dx; voxel_size.y = dy; voxel_size.z = dz;}
	void RemoveTerminalComponents();
	void RemoveBoundaryComponents();
	void GetRadiusFromVolume(const char* filename, unsigned int threshold, unsigned int range);
	void GetRadiusFromDistanceMap(const char *filename, int range);
	point3D<float> GetFiberRadius(unsigned int fiber);
	void ScaleNetwork(float x, float y, float z);
	void SetOutputPosition(float x, float y, float z);
	

};
//private
void rtsNetwork::reset()
{
	vector<EdgeType>::iterator e;
	for(e=EdgeList.begin(); e!=EdgeList.end(); e++)
		(*e).valid = true;

	vector<NodeType>::iterator n;
	for(n=NodeList.begin(); n!=NodeList.end(); n++)
		(*n).valid = true;
}

float rtsNetwork::calcFiberBend(unsigned int f)
{
	/*calculates the bend of a fiber.*/

	//find the line connecting the two endpoints
	point3D<float> p1 = toTissue(NodeList[EdgeList[f].v0].p);
	point3D<float> p2 = toTissue(NodeList[EdgeList[f].v1].p);

	//find the largest distance between the filament and the line between endpoints
	FilamentType::iterator i;
	float max_distance = 0;
	for(i = network[f].begin(); i!= network[f].end(); i++)
	{
		point3D<float> p0 = toTissue((*i));

		//find the distance between the current point and the line made by the endpoints
		float d = ((p0 - p1).X(p0-p2).Length())/(p2-p1).Length();
		if(d > max_distance)
			max_distance = d;
	}
	float fiber_length = calcFiberLength(f);
	if(max_distance/fiber_length > 0.5)
	{
		cout<<"Something's wrong."<<endl;
		cout<<max_distance<<"---"<<fiber_length<<endl;
	}

	return max_distance/calcFiberLength(f);	

}
int rtsNetwork::cleanRedundantEdges()
{
	
	//redundant edges are edges with nodes of valence-2

	int combined = 0;

	//combine all valence-2 fibers
	NetworkType newNetwork;
	unsigned int numNodes = NodeList.size();
	unsigned int n;
	unsigned int e0, e1;
	for(n=0; n<numNodes; n++)
	{
		//if the node is valence-2
		if(getConnections(n).size() == 2)
		{
			//get both edges
			e0 = NodeList[n].edges[0];
			e1 = NodeList[n].edges[1];

			//if both edges are valid
			if(EdgeList[e0].valid && EdgeList[e1].valid)
			{
				//combine them
				if(EdgeList[e0].v0 == n)
				{
					if(EdgeList[e1].v0 == n)
						newNetwork.push_back(combineFibers(network[e0], 0, network[e1], 0));
					else
						newNetwork.push_back(combineFibers(network[e0], 0, network[e1], 1));
				}
				else if(EdgeList[e0].v1 == n)
				{
					if(EdgeList[e1].v0 == n)
						newNetwork.push_back(combineFibers(network[e0], 1, network[e1], 0));
					else
						newNetwork.push_back(combineFibers(network[e0], 1, network[e1], 1));
				}
				EdgeList[e0].valid = false;
				EdgeList[e1].valid = false;
				combined++;
			}

		}
	}

	unsigned int numEdges = network.size();
	unsigned int e;
	for(e=0; e<numEdges; e++)
	{
		if(EdgeList[e].valid)
			newNetwork.push_back(network[e]);
	}

	network = newNetwork;
	calcGraph();
	return combined;


}
int rtsNetwork::cleanBarbs(float length)
{
	int removed = EdgeList.size();
	//removes fibers that are shorter than the specified length
	NetworkType newNetwork;

	unsigned int numEdges, e;
	numEdges = EdgeList.size();
	unsigned int numVertices, v;
	float e_length;

	for(e=0; e<numEdges; e++)
	{
		//if the fiber is an end fiber (one vertex has valence 1)
		if(getConnections(EdgeList[e].v0).size() ==1 || getConnections(EdgeList[e].v1).size() ==1)
		{
			e_length = calcFiberLength(e);
			//if the fiber is larger than the given length
			if(e_length > length)
			{
				newNetwork.push_back(network[e]);
				removed--;
			}
			
		}
		//if the fiber is not an end fiber
		else
		{
			newNetwork.push_back(network[e]);
			removed--;
		}

	}

	network.clear();
	network = newNetwork;
	calcGraph();

	cout<<"Barbs removed: "<<removed<<endl;

	return removed;

}

int rtsNetwork::cleanDegenerateEdges(float length)
{
	/*This function cleans unnecessary edges with length less than <length>*/

	//form lists of degenerate and non-degenerate edges
	int removed = 0;

	list<ConnectableType> degenerate;
	list<ConnectableType> complete;

	//iterate through each fiber and put it in the appropriate list
	int numFibers = network.size();
	int f;
	ConnectableType c;
	for(f=0; f<numFibers; f++)
	{
		c.fiber = network[f];
		c.id = f;
		if(calcFiberLength(f) <= length)
			degenerate.push_back(c);
		else
			complete.push_back(c);
	}

	
	vector<unsigned int> connected_edges;
	vector<unsigned int> shared_vertices;
	int numEdges, numVerts, e, e0, e1, v;
	int connectedEdge;
	//for each degenerate edge
	list<ConnectableType>::iterator i;
	for(i=degenerate.begin(); i!=degenerate.end(); i++)
	{
/*DEGENERATE CASE 1
In this case, we remove edges that are parts of very small cycles (ex short edges that
bridge between a bifurcation, forming a small triangle).
*/
		//get all of the edges connected to that edge
		connected_edges.clear();
		//add edges connected to v0
		v = EdgeList[(*i).id].v0;
		numEdges = NodeList[v].edges.size();
		for(e=0; e<numEdges; e++)
		{
			connectedEdge = NodeList[v].edges[e];
			if(connectedEdge != (*i).id && EdgeList[connectedEdge].valid)
				connected_edges.push_back(connectedEdge);
		}
		//add edges connected to v1
		v = EdgeList[(*i).id].v1;
		numEdges = NodeList[v].edges.size();
		for(e=0; e<numEdges; e++)
		{
			connectedEdge = NodeList[v].edges[e];
			if(connectedEdge != (*i).id && EdgeList[connectedEdge].valid)
				connected_edges.push_back(connectedEdge);
		}

		//find all of the shared vertices
		shared_vertices.clear();
		numEdges = connected_edges.size();
		for(e0=0; e0<numEdges; e0++)
		{
			for(e1 = e0; e1<numEdges; e1++)
			{
				if(e0 != e1)
				{
					if(EdgeList[connected_edges[e0]].v0 == EdgeList[connected_edges[e1]].v0 ||
					   EdgeList[connected_edges[e0]].v0 == EdgeList[connected_edges[e1]].v1)
						shared_vertices.push_back(EdgeList[connected_edges[e0]].v0);
					else if(EdgeList[connected_edges[e0]].v1 == EdgeList[connected_edges[e1]].v0 ||
							EdgeList[connected_edges[e0]].v1 == EdgeList[connected_edges[e1]].v1)
						shared_vertices.push_back(EdgeList[connected_edges[e0]].v1);
				}
			}
		}

		//are any of these vertices not part of the original edge?
		numVerts = shared_vertices.size();
		bool delete_fiber = false;
		for(v = 0; v<numVerts; v++)
		{
			if(shared_vertices[v] != EdgeList[(*i).id].v0 &&
			   shared_vertices[v] != EdgeList[(*i).id].v1)
			{
				if(getConnections(EdgeList[(*i).id].v0).size() ==1 ||
				   getConnections(EdgeList[(*i).id].v0).size() ==1)
				   cout<<"error  "<<(*i).id<<endl;
				delete_fiber = true;
				removed++;
				EdgeList[(*i).id].valid = false;
			}
		}
/*DEGENERATE CASE 2
In this case, we remove edges that are equal to each other (ex both ends are connected to the
same outside fiber).
*/
		
		//find the number of branches in v0
		int v0 = EdgeList[(*i).id].v0;
		int v1 = EdgeList[(*i).id].v1;

		int numV0 = NodeList[v0].edges.size();
		int numV1 = NodeList[v1].edges.size();

		//go through each fiber connected to both nodes
		for(int vi=0; vi<numV0; vi++)
			for(int vj=0; vj<numV1; vj++)
			{
				//if there is a match and it isn't the current fiber, delete it
				if(NodeList[v0].edges[vi] == NodeList[v1].edges[vj])
					if(NodeList[v0].edges[vi] != (*i).id && EdgeList[NodeList[v0].edges[vi]].valid)
					{
						delete_fiber = true;
						removed++;
						EdgeList[(*i).id].valid = false;
					}
			}

/*DEGENERATE CASE 3
Delete fibers that are less than length and have degree-4 connections
*/
		if(calcFiberLength((*i).id) < length && 
			NodeList[EdgeList[(*i).id].v1].edges.size() >3 &&
			NodeList[EdgeList[(*i).id].v0].edges.size() >3)
		{
			delete_fiber = true;
			removed++;
			EdgeList[(*i).id].valid = false;
		}




		if(!delete_fiber)
			complete.push_back((*i));

	}

	//re-create the network from <complete>
	NetworkType newNetwork;
	for(i = complete.begin(); i!=complete.end(); i++)
		newNetwork.push_back((*i).fiber);

	network = newNetwork;
	calcGraph();

	return removed;




}
bool rtsNetwork::isBoundaryNode(unsigned int node)
{
	CoordType p;

	p = NodeList[node].p;
	if(NodeList[node].edges.size() > 1)
		return false;
	if(p.x > BUFFER_ZONE && p.y > BUFFER_ZONE && p.z > BUFFER_ZONE &&
	   p.x < boundary.x - BUFFER_ZONE && p.y < boundary.y - BUFFER_ZONE && p.z < boundary.z - BUFFER_ZONE)
		return false;
	else
		return true;
}


bool rtsNetwork::isConnectable(ConnectableType c0, ConnectableType c1, bool& c0_end, bool& c1_end, float distance, float cosine_angle)
{
	//tests to see if two ConnectableType fibers can be connected.  If so, the function
	//returns true as well as the optimal connected endpoints as <c0_end> and <c1_end>.
	//For these parameters, 0 = the first point on the fiber and 1 = the last point.

	if(c0.id == 1381 && c1.id == 1407)
		cout<<"here"<<endl;


	vector3D<float> t0;
	vector3D<float> t1;
	if(c0.front == 0)
	{
		t0 = calcTrajectory(c0.fiber, 0);
		if(c1.front == 0)
		{
			t1 = calcTrajectory(c1.fiber, 0);
			if(compareTrajectories(c0.fiber.front(), t0, c1.fiber.front(), t1, distance, cosine_angle))
			{
				c0_end = 0;
				c1_end = 0;
				return true;
			}
		}
		if(c1.back == 0)
		{
			t1 = calcTrajectory(c1.fiber, 1);
			if(compareTrajectories(c0.fiber.front(), t0, c1.fiber.back(), t1, distance, cosine_angle))
			{
				c0_end = 0;
				c1_end = 1;
				return true;
			}
		}
	}
	if(c0.back == 0)
	{
		t0 = calcTrajectory(c0.fiber, 1);
		if(c1.front == 0)
		{
			t1 = calcTrajectory(c1.fiber, 0);
			if(compareTrajectories(c0.fiber.back(), t0, c1.fiber.front(), t1, distance, cosine_angle))
			{
				c0_end = 1;
				c1_end = 0;
				return true;
			}
		}
		if(c1.back == 0)
		{
			t1 = calcTrajectory(c1.fiber, 1);
			if(compareTrajectories(c0.fiber.back(), t0, c1.fiber.back(), t1, distance, cosine_angle))
			{
				c0_end = 1;
				c1_end = 1;
				return true;
			}
		}
	}


	return false;

}
FilamentType rtsNetwork::combineFibers(FilamentType f0, bool endpoint0, FilamentType f1, bool endpoint1)
{

	//create the new filament
	FilamentType newFilament;
	//insert the first filament
	//if the valence-1 vertex is at the beginning
	if(endpoint0 == 0)
	{
		//add the edge backwards
		for(int v=f0.size() - 1; v>=0; v--)
			newFilament.push_back(f0[v]);
	}
	else
	{
		//otherwise we can just copy the filament over
		newFilament = f0;
	}

	//insert the second filament
	//if the valence-1 vertex is at the beginning
	if(endpoint1 == 0)
	{
		//add the edge in forwards
		for(int v=0; v<f1.size(); v++)
			newFilament.push_back(f1[v]);
	}
	else
	{
		for(int v=f1.size() - 1; v>=0; v--)
			newFilament.push_back(f1[v]);
	}

	return newFilament;
}
vector3D<float> rtsNetwork::calcTrajectory(FilamentType f, bool end)
{
	//calculates the trajectory of a fiber at the given endpoint (0 = first, 1 = last)

	CoordType pEndpoint;
	CoordType pPrevpoint;
	if(end == 0)
	{
		pEndpoint = f.front();
		pPrevpoint = f[1];
		//pPrevpoint = f[ceil(f.size()/4.0)];
	}
	else
	{
		pEndpoint = f.back();
		pPrevpoint = f[f.size() - 2];
		//pPrevpoint = f[floor(f.size()*3.0/4.0)];
	}

	//pPrevpoint = f[f.size()/2];
	

	vector3D<float> result = toTissue(pEndpoint) - toTissue(pPrevpoint);
	result.Normalize();
	return result;
}
bool rtsNetwork::compareTrajectories(CoordType p0, vector3D<float> v0, CoordType p1, vector3D<float> v1, float distance, float cosine_angle)
{
	//determines of two points/trajectory pairs are compatible for connection
	vector3D<float> difference = p1 - p0;
	//test distance
	if(difference.Length() <= distance)
	{
		difference.Normalize();
		if(v0*difference >= cosine_angle && v1*difference <= -cosine_angle)
			return true;
	}
	return false;
}
vector<unsigned int> rtsNetwork::getConnections(unsigned int node)
{
	return NodeList[node].edges;

}
point3D<float> rtsNetwork::toTissue(CoordType grid_coord)
{
	//converts a coordinate from the original voxel grid coordinates to tissue-space coordinates
	return CoordType(grid_coord.x * voxel_size.x,
					 grid_coord.y * voxel_size.y,
				     grid_coord.z * voxel_size.z);
					 
	//return grid_coord;
}

float rtsNetwork::calcFiberLength(unsigned int f)
{
	//This function calculates the total length of the given fiber

	float result = 0.0;

	//find the length of each fiber
	int vertNum, v;

	vertNum = network[f].size();
	for(v=1; v<vertNum; v++)
	{
		result += (toTissue(network[f][v-1]) - toTissue(network[f][v])).Length();
	}

	return result;

}

float rtsNetwork::calcTotalLength()
{
	float result = 0.0;

	//find the length of each fiber
	int numFibers = network.size();
	int f;

	for(f = 0; f < numFibers; f++)
	{
		result += calcFiberLength(f);		
	}
	return result;

}

rtsOBJ rtsNetwork::LoadOBJModel(const char *filename)
{
	rtsOBJ model;
	model.LoadFile(filename);

	//empty the current network if it isn't already
	network.clear();

	//get the number of filaments
	int lineNum = model.getNumLines();

	int vertexNum, v, vIndex;		//for keeping track of the vertices in each line

	for(int l=0; l<lineNum; l++)
	{
		FilamentType fiber;							//create a new fiber
		vertexNum = model.getNumLineVertices(l);	//get the number of vertices in the line

		for(v = 0; v<vertexNum; v++)				//convert the line to a fiber
		{
			vIndex = model.getLineVertex(l, v);
			fiber.push_back(model.getVertex3d(vIndex));
		}
		network.push_back(fiber);					//store the converted fiber in the network
	}

	//calculate the graph components
	calcGraph();

	//set the boundary
	AABB bound = model.getBoundingBox();
	boundary.x = ceil(bound.max.x);
	boundary.y = ceil(bound.max.y);
	boundary.z = ceil(bound.max.z);

	return model;
	
}

void rtsNetwork::AddOBJModel(const char *filename)
{
	rtsOBJ model;
	model.LoadFile(filename);

	//get the number of filaments
	int lineNum = model.getNumLines();

	int vertexNum, v, vIndex;		//for keeping track of the vertices in each line

	for(int l=0; l<lineNum; l++)
	{
		FilamentType fiber;							//create a new fiber
		vertexNum = model.getNumLineVertices(l);	//get the number of vertices in the line

		for(v = 0; v<vertexNum; v++)				//convert the line to a fiber
		{
			vIndex = model.getLineVertex(l, v);
			fiber.push_back(model.getVertex3d(vIndex));
		}
		network.push_back(fiber);					//store the converted fiber in the network
	}

	//calculate the graph components
	calcGraph();

	/*//set the boundary
	AABB bound = model.getBoundingBox();
	boundary.x = ceil(bound.max.x);
	boundary.y = ceil(bound.max.y);
	boundary.z = ceil(bound.max.z);
	*/

}

void rtsNetwork::outputFiberToSWC(int e, int node_from, int parent_id)
{
	EdgeList[e].valid = false;
	int p;
	point3D<float> outPoint;
	//if we are traversing the fiber front-to-back
	if(EdgeList[e].v0 == node_from)
	{
		if(parent_id == -1)
		{
			p=0;
			outPoint = network[e][p] + position;
			outfile<<id<<" "<<2<<" "<<outPoint.x<<" "<<outPoint.y<<" "<<outPoint.z<<" "<<1.0<<" "<<parent_id<<endl;
		}
		else
		{
			p=1;
			outPoint = network[e][p] + position;
			outfile<<id<<" "<<2<<" "<<outPoint.x<<" "<<outPoint.y<<" "<<outPoint.z<<" "<<1.0<<" "<<parent_id<<endl;
		}
		p++;
		id++;
		for(; p<network[e].size(); p++)
		{
			outPoint = network[e][p] + position;
			outfile<<id<<" "<<2<<" "<<outPoint.x<<" "<<outPoint.y<<" "<<outPoint.z<<" "<<1.0<<" "<<id-1<<endl;
			id++;
		}
	}
	//if we are traversing the fiber back-to-front
	else
	{
		if(parent_id == -1)
		{
			p=network[e].size()-1;
			outPoint = network[e][p] + position;
			outfile<<id<<" "<<2<<" "<<outPoint.x<<" "<<outPoint.y<<" "<<outPoint.z<<" "<<1.0<<" "<<parent_id<<endl;
		}
		else
		{
			p=network[e].size()-2;
			outPoint = network[e][p] + position;
			outfile<<id<<" "<<2<<" "<<outPoint.x<<" "<<outPoint.y<<" "<<outPoint.z<<" "<<1.0<<" "<<parent_id<<endl;
		}
		p--;
		id++;
		for(; p>=0; p--)
		{
			outPoint = network[e][p] + position;
			outfile<<id<<" "<<2<<" "<<outPoint.x<<" "<<outPoint.y<<" "<<outPoint.z<<" "<<1.0<<" "<<id-1<<endl;
			id++;
		}
	}

}
void rtsNetwork::traverseGraph(int node, int parent_id)
{
	//if the node is valid
	if(NodeList[node].valid)
	{
		NodeList[node].valid = false;
		int num_edges = NodeList[node].edges.size();
		for(int e = 0; e<num_edges; e++)
		{
			if(EdgeList[NodeList[node].edges[e]].valid)
			{
				outputFiberToSWC(NodeList[node].edges[e], node, parent_id);
				int newNode = EdgeList[NodeList[node].edges[e]].v0;
				if(newNode != node)
					traverseGraph(newNode, id-1);
				else
					traverseGraph(EdgeList[NodeList[node].edges[e]].v1, id-1);
			}
		}
	}


}

void rtsNetwork::SetOutputPosition(float x, float y, float z)
{
	position.x = x;
	position.y = y;
	position.z = z;
}
void rtsNetwork::SaveSWCModel(const char *filename, float root_x, float root_y)
{
	cout<<filename<<endl;
	outfile.open(filename);
	//recursively iterate through the network saving the elements to an SWC file
	int node;
	id = 1;

	//set the root seed point
	point3D<float> rootSeed(root_x, root_y, 0);
	float distance = 9999999;
	int closest_node;

	//find the closest node to the root seed
	for(node = 0; node != NodeList.size(); node++)
	{
		float length = (NodeList[node].p - rootSeed).Length();
		if(length < distance)
		{
			closest_node = node;
			distance = length;
			cout<<length<<endl;
		}
	}
	traverseGraph(closest_node, -1);

	outfile.close();
	
}

void rtsNetwork::SaveSWCModel(const char* filename)
{
	cout<<filename<<endl;
	outfile.open(filename);
	//recursively iterate through the network saving the elements to an SWC file
	int node;
	id = 1;

	//set the root seed point
	
	//find the closest node to the root seed
	for(node = 0; node != NodeList.size(); node++)
	{
		if(NodeList[node].valid)
			traverseGraph(node, -1);
	}
	

	outfile.close();
}

void rtsNetwork::ScaleNetwork(float x, float y, float z)
{

	NetworkType::iterator f;
	FilamentType::iterator p;

	for(f = network.begin(); f!=network.end(); f++)
		for(p = (*f).begin(); p!= (*f).end(); p++)
		{
			(*p).x = (*p).x * x;
			(*p).y = (*p).y * y;
			(*p).z = (*p).z * z;
		}

	calcGraph();
}
rtsOBJ rtsNetwork::CreateModel()
{
	rtsOBJ model;
	/*
	NetworkType::iterator f;
	unsigned int numVerts, v;
	for(f=network.begin(); f!=network.end(); f++)
	{
		model.objBegin(OBJ_LINE_STRIP);
		
		numVerts = (*f).size();
		for(v=0; v<numVerts; v++)
		{
			model.objVertex3f((*f)[v].x, (*f)[v].y, (*f)[v].z);
		}
		
		model.objEnd();
	}
	*/

	//add all of the vertices to the OBJ
	int num_vertices = NodeList.size();
	int v;
	for(v=0; v<num_vertices; v++)
	{
		model.insertVertexPosition(NodeList[v].p.x, NodeList[v].p.y, NodeList[v].p.z);
	}

	//for each edge
	int num_edges = EdgeList.size();
	int e;
	unsigned int* buffer = new unsigned int[10000];
	
	for(e=1; e<num_edges-1; e++)
	{	
		buffer[0] = EdgeList[e].v0;
		
		//walk along the fiber
		num_vertices = network[e].size()-1;
		for(v=1; v<num_vertices; v++)
		{
			buffer[v] = model.getNumVertices();
			model.insertVertexPosition(network[e][v].x, network[e][v].y, network[e][v].z);
		}
		
		buffer[v] = EdgeList[e].v1;

		model.insertLine(v+1, buffer, NULL, NULL);
	}
	

	return model;

}



void rtsNetwork::calcGraph()
{
	//make sure that all graph information is cleared
	EdgeList.clear();
	NodeList.clear();

	//find the maximum extents of the fiber network
	int numFibers = network.size();
	int f;
	float Max = 0;
	float temp;
	for(f=0; f<numFibers; f++)
	{
		temp = max(network[f].front().x,
			   max(network[f].front().y,
			   max(network[f].front().z,
			   max(network[f].back().x,
			   max(network[f].back().y, 
			   network[f].back().z)))));
		if(temp > Max) Max = temp;
	}

	//place each point ID inside the octree
	unsigned int id = 0;
	Octree<unsigned int> tree(OCTREE_SIZE);
	tree.setEmptyValue(UINT_MAX);
	unsigned int v0, v1;


	for(f=0; f<numFibers; f++)
	{
		//get the appropriate id from the tree
		NodeType v0pos, v1pos;
		v0pos.p = network[f].front();
		v1pos.p = network[f].back();
		v0pos.valid = 1;
		v1pos.valid = 1;
		v0 = tree(v0pos.p.x, v0pos.p.y, v0pos.p.z);
		v1 = tree(v1pos.p.x, v1pos.p.y, v1pos.p.z);

		//if the id doesn't exist, add it
		if(v0 == UINT_MAX)
		{
			v0 = id;								//set the id of the current vertex
			tree(v0pos.p.x, v0pos.p.y, v0pos.p.z) = id;	//add the id to the octree
			NodeList.push_back(v0pos);	//add the node position to the node list

			id++;									//increment the id
		}
		if(v1 == UINT_MAX)
		{
			v1 = id;								//set the id of the current vertex
			tree(v1pos.p.x, v1pos.p.y, v1pos.p.z) = id;	//add the id to the octree
			NodeList.push_back(v1pos);	//add the node position to the node list

			id++;									//increment the id
		}

		//add the current edge to the edge list
		EdgeType e;
		e.v0 = v0;
		e.v1 = v1;
		e.valid = 1;
		e.length = calcFiberLength(f);
		EdgeList.push_back(e);
		//add this edge to each point's edge list
		NodeList[e.v0].edges.push_back(EdgeList.size() - 1);
		NodeList[e.v1].edges.push_back(EdgeList.size() - 1);

		//cout<<"Edge: "<<v0<<"--->"<<v1<<"  "<<"w = "<<e.w<<endl;
	}

	/*Computes a connectivity matrix based on the edge list and vertices in Attributes*/

	/*
	//create the matrix and zero it out
	Array2D<unsigned int> matrix(NodeList.size(), NodeList.size());
	unsigned int numNodes = NodeList.size();
	int x, y;
	for(x=0; x<numNodes; x++)
		for(y=0; y<numNodes; y++)
			matrix(x, y) = 0;

	//iterate through each edge, assigning values to the connectivity matrix
	vector<EdgeType>::iterator e;
	for(e = EdgeList.begin(); e!= EdgeList.end(); e++)
	{
		matrix((*e).v0, (*e).v1) = 1;
		matrix((*e).v1, (*e).v0) = 1;
	}

	ConnectivityMatrix = matrix;
	*/
}

vector<BranchType> rtsNetwork::getFiberAngles(unsigned int fiber)
{
	vector<BranchType> result;

	int b;
	vector3D<float> currentTrajectory;
	vector3D<float> branchTrajectory;
	unsigned int f;

	unsigned int node = EdgeList[fiber].v0;

	int numBranches = NodeList[node].edges.size();	
	currentTrajectory = calcTrajectory(network[fiber], 0).Normalize();
	//currentTrajectory.print();
	for(b=0; b<numBranches; b++)
	{
		f = NodeList[node].edges[b];
		if(f != fiber && EdgeList[f].valid)
		{
			if(EdgeList[f].v0 == node)
				branchTrajectory = calcTrajectory(network[f], 0).Normalize();
			else
				branchTrajectory = calcTrajectory(network[f], 1).Normalize();
			//branchTrajectory.print();
			BranchType newBranch;
			newBranch.angle = acos((branchTrajectory*(-1))*currentTrajectory)*(180.0/3.14159);
			newBranch.branch_id = f;
			result.push_back(newBranch);
		}
		
	}

	node = EdgeList[fiber].v1;

	numBranches = NodeList[node].edges.size();	
	currentTrajectory = calcTrajectory(network[fiber], 1).Normalize();
	//currentTrajectory.print();
	for(b=0; b<numBranches; b++)
	{
		f = NodeList[node].edges[b];
		if(f != fiber && EdgeList[f].valid)
		{
			if(EdgeList[f].v0 == node)
				branchTrajectory = calcTrajectory(network[f], 0).Normalize();
			else
				branchTrajectory = calcTrajectory(network[f], 1).Normalize();
			//branchTrajectory.print();
			BranchType newBranch;
			newBranch.angle = acos((branchTrajectory*(-1))*currentTrajectory)*(180.0/3.14159);
			newBranch.branch_id = f;
			result.push_back(newBranch);
		}
		
	}
	
	//mark the fiber as traversed
	EdgeList[fiber].valid = false;

	return result;
}
//public
void rtsNetwork::CalculateAttributes()
{
	Attributes.TotalFiberLength = calcTotalLength();
	Attributes.NumFibers = network.size();
	
	Attributes.NumPoints = 0;
	Attributes.TotalVolume = 0;
	Attributes.NumDisconnectedFibers = 0;
	Attributes.NumBoundaryFibers = 0;
	int numFibers = network.size();
	for(int f=0; f<numFibers; f++)
	{
		Attributes.NumPoints += network[f].size();
		Attributes.TotalVolume += EdgeList[f].volume;
		if(NodeList[EdgeList[f].v0].edges.size() == 1 ||
			NodeList[EdgeList[f].v1].edges.size() ==1)
			if(isBoundaryNode(EdgeList[f].v0) || isBoundaryNode(EdgeList[f].v1))
				Attributes.NumBoundaryFibers++;
			else
				Attributes.NumDisconnectedFibers++;

	}

	Attributes.NumBifurcations = 0;
	int numNodes = NodeList.size();
	for(int n=0; n<numNodes; n++)
		if(NodeList[n].edges.size() >= 2)
			Attributes.NumBifurcations++;


}





rtsOBJ rtsNetwork::Repair(float max_distance, float cos_angle)
{
	/*This function repairs the network by serching a solid angle around each endpoint
	for a suitable candidate endpoint to connect to.*/

	cout<<"repairing...."<<endl;

	rtsOBJ preview;

	//create the complete and incomplete connectable lists
	list<ConnectableType> complete;
	list<ConnectableType> incomplete;

	cout<<"Finding Valence-1 Fibers..."<<endl;
	//go through each fiber, adding ones with valence-1 endpoints to <incomplete>
	int numFibers = network.size();
	int f;
	unsigned int v0;
	unsigned int v1;
	int numVerts, v;
	for(f=0; f<numFibers; f++)
	{
		ConnectableType c;
		v0 = EdgeList[f].v0;
		v1 = EdgeList[f].v1;
		//test valence and boundary
		if(NodeList[v0].edges.size() == 1 && !isBoundaryNode(v0))
		{
			c.front = false;
		}
		else c.front = true;

		if(NodeList[v1].edges.size() == 1 && !isBoundaryNode(v1))
		{
			c.back = false;
		}
		else c.back = true;

		if(c.front == false || c.back == false)
		{
			c.fiber = network[f];
			c.id = f;
			incomplete.push_back(c);
			numVerts = c.fiber.size();
			//preview.objBegin(OBJ_LINE_STRIP);
			//for(v=0; v<numVerts; v++)
			//	preview.objVertex3f(c.fiber[v].x, c.fiber[v].y, c.fiber[v].z);
			//preview.objEnd();
		}
		else
		{
			c.fiber = network[f];
			c.id = f;
			complete.push_back(c);
		}
		
	}
	cout<<"Valence-1 Fibers: "<<incomplete.size()<<endl;

	//find all of the possible connections
	list<ConnectableType>::iterator i;
	list<ConnectableType>::iterator j;
	list<ConnectableType>::iterator temp;
	FilamentType newFiber;
	bool endpoint0, endpoint1;
	
	int iIndex, jIndex;

	iIndex = 0;
	for(i=incomplete.begin(); i!=incomplete.end(); i++)
	{
		jIndex = iIndex;
		j = i;
		while(j != incomplete.end() && i != incomplete.end())
		{
			if(i!=j)
			{
				if(isConnectable((*i), (*j), endpoint0, endpoint1, max_distance, cos_angle))
				{
					ConnectableType newConnectable;
					newConnectable.fiber = combineFibers((*i).fiber, endpoint0, (*j).fiber, endpoint1);
					//determine the new connectivity
					if(endpoint0 == 0)
						newConnectable.front = (*i).back;
					else
						newConnectable.front = (*i).front;

					if(endpoint1 == 0)
						newConnectable.back = (*j).back;
					else
						newConnectable.back = (*j).front;

					//remove the previous connectables
					temp = i;
					i++;
					iIndex++;
					incomplete.erase(temp);
					temp = j;
					//this is just in case the fibers are next to each otehr
					if(j == i)
					{
						i++;
						iIndex++;
						j = i;
						jIndex = iIndex;
					}
					else
					{
						j = i;
						jIndex = iIndex;
					}
					incomplete.erase(temp);
					//insert the new connectable in the proper list
					if(newConnectable.front == 0 || newConnectable.back == 0)
						incomplete.push_back(newConnectable);
					else
						complete.push_back(newConnectable);
					
					

					//numVerts = newConnectable.fiber.size();
					//preview.objBegin(OBJ_LINE_STRIP);
					//for(v=0; v<numVerts; v++)
					//	preview.objVertex3f(newConnectable.fiber[v].x, newConnectable.fiber[v].y, newConnectable.fiber[v].z);
					//preview.objEnd();
				}
				else
				{
					j++;
					jIndex++;
				}
				
			}
			else
			{
				j++;
				jIndex++;
			}
			//cout<<"i: "<<iIndex<<"j: "<<jIndex<<endl;
			//cout<<"size of incomplete: "<<incomplete.size()<<endl;
		}
		iIndex++;
	}

	//reconstruct the fiber list
	vector<FilamentType> newNetwork;
	for(i=incomplete.begin(); i!=incomplete.end(); i++)
		newNetwork.push_back((*i).fiber);
	for(i=complete.begin(); i!=complete.end(); i++)
		newNetwork.push_back((*i).fiber);

	network = newNetwork;
	calcGraph();

	
	//build the model to return
	preview = CreateModel();
	
	return preview;

}
void rtsNetwork::printFiberStats(unsigned int fiber)
{
	/*
	cout<<"-----------------------------------------------------"<<endl;
	cout<<"Fiber nodes:"<<endl;
	int numVerts = network[fiber].size();
	for(int v=0; v<numVerts; v++)
	{
		cout<<network[fiber][v].x<<","<<network[fiber][v].y<<","<<network[fiber][v].z<<endl;
	}
	*/
	cout<<"------------------------------------------------------"<<endl;
	cout<<"Selected node: "<<fiber<<endl;


	int b;
	vector3D<float> currentTrajectory;
	vector3D<float> branchTrajectory;
	unsigned int f;

	//get the list of angles
	vector<BranchType> angles = getFiberAngles(fiber);

	unsigned int node = EdgeList[fiber].v0;
	int numBranches = NodeList[node].edges.size()-1;	
	cout<<"V0 Valence: "<<getConnections(node).size()<<endl; 

	for(b=0; b<numBranches; b++)
	{
		cout<<"     ["<<angles[b].branch_id<<"] "<<angles[b].angle<<(char)248<<endl;
	}


	node = EdgeList[fiber].v1;	
	cout<<"V1 Valence: "<<getConnections(node).size()<<endl;
	numBranches = angles.size();
	for(; b<numBranches; b++)
	{
		cout<<"     ["<<angles[b].branch_id<<"] "<<angles[b].angle<<(char)248<<endl;
	}

	reset();

	cout<<"End Points: ["<<NodeList[EdgeList[fiber].v0].p.x<<","
						<<NodeList[EdgeList[fiber].v0].p.y<<","
						<<NodeList[EdgeList[fiber].v0].p.z<<"]"
						<<"["<<NodeList[EdgeList[fiber].v1].p.x<<","
						<<NodeList[EdgeList[fiber].v1].p.y<<","
						<<NodeList[EdgeList[fiber].v1].p.z<<"]"<<endl;
	cout<<"Fiber Length: "<<calcFiberLength(fiber)<<"um"<<endl;
	cout<<"Average Radius: "<<EdgeList[fiber].avg_radius<<"um"<<endl;
	cout<<"Min Radius: "<<EdgeList[fiber].min_radius<<"um"<<endl;
	cout<<"Max Radius: "<<EdgeList[fiber].max_radius<<"um"<<endl;
	cout<<"Fiber Volume: "<<EdgeList[fiber].volume<<"um^3"<<endl;
	cout<<"Maximum Bend: "<<calcFiberBend(fiber)<<endl;
	cout<<"------------------------------------------------------"<<endl;

	int numVertices = network[fiber].size();
	//for(int v=0; v<numVertices; v++)
	//	cout<<network[fiber][v].x<<","<<network[fiber][v].y<<","<<network[fiber][v].z<<endl;


}
void rtsNetwork::Clean(float degenerate_length = 5, float barb_length = 30)
{
	int changes;
	do
	{
		changes = 0;
		int degen = cleanDegenerateEdges(degenerate_length);
		cout<<"Degen: "<<degen<<endl;
		changes += degen;
		
		int barbs = cleanBarbs(barb_length);
		cout<<"Barbs: "<<barbs<<endl;
		changes += barbs;
		
		int redundant = cleanRedundantEdges();
		cout<<"Redundant: "<<redundant<<endl;
		changes += redundant;

	}while(changes != 0);


}
void rtsNetwork::ReMesh(float resample_length)
{
	//Resample the skeleton at the given frequency
	NetworkType newNetwork;

	int numFibers = network.size();
	int f;
	int numVerts, v;
	FilamentType newFiber;
	vector3D<float> distanceVector;
	point3D<float> prevPoint;
	for(f=0; f<numFibers; f++)
	{
		newFiber.clear();
		numVerts = network[f].size();
		newFiber.push_back(network[f][0]);
		prevPoint = network[f][0];
		for(v=1; v<numVerts-1; v++)
		{
			distanceVector = toTissue(network[f][v]) - toTissue(prevPoint);
			if(distanceVector.Length() >= resample_length)
			{
				newFiber.push_back(network[f][v]);
				prevPoint = network[f][v];
			}
		}
		
		//deal with the last vertex
		//newFiber.push_back(network[f][v]);
		
		distanceVector = toTissue(network[f][v]) - toTissue(prevPoint);
		if(distanceVector.Length() >= resample_length || newFiber.size() == 1)
		{
			newFiber.push_back(network[f][v]);
		}
		else
		{
			//change the last point
			newFiber[newFiber.size() - 1] = network[f][v];
		}
		
		
		//cout<<"old: "<<network[f].size()<<endl;
		//cout<<"new: "<<newFiber.size()<<endl;
		newNetwork.push_back(newFiber);
	}

	//replace the network
	//cout<<"old network: "<<network.size()<<endl;
	//cout<<"new network: "<<newNetwork.size()<<endl;
	network = newNetwork;
	calcGraph();

}
void rtsNetwork::SaveAngles(const char *filename)
{
	ofstream outfile;
	outfile.open(filename);

	int errors = 0;

	//get the list of angles for each fiber
	unsigned int f;
	int numFibers = EdgeList.size();
	for(f=0; f<numFibers; f++)
	{
		vector<BranchType> angles = getFiberAngles(f);
		for(int a=0; a<angles.size(); a++)
		{
			outfile<<angles[a].angle<<endl;

			//print out errors
			if(angles[a].angle >179 && angles[a].angle < 181)
			{
				//cout<<"ERROR"<<endl;
				//printFiberStats(f);
				errors++;
			}
			if(angles[a].angle >90 && angles[a].angle < 90.01)
			{
				//cout<<"ERROR"<<endl;
				//printFiberStats(f);
				errors++;
			}
			if(angles[a].angle < 0)
			{
				//cout<<"ERROR"<<endl;
				//printFiberStats(f);
				errors++;
			}
		}
	}

	cout<<"ERRORS: "<<errors<<endl;

	outfile.close();

	reset();

}
void rtsNetwork::SaveBends(const char *filename)
{
	ofstream outfile;
	outfile.open(filename);

	int numFibers = EdgeList.size();
	for(int f = 0; f<numFibers; f++)
		if(network[f].size() >2)
			outfile<<calcFiberBend(f)<<endl;

	outfile.close();

}

void rtsNetwork::SaveLengths(const char *filename)
{
	ofstream outfile;
	outfile.open(filename);

	int numFibers = EdgeList.size();
	for(int f = 0; f<numFibers; f++)
		if(network[f].size() >2)
			outfile<<calcFiberLength(f)<<endl;

	outfile.close();
}

void rtsNetwork::SaveMathematicaGraph(const char* filename)
{
	ofstream outfile;
	outfile.open(filename);

	outfile<<"<< Combinatorica`"<<endl;
	outfile<<"<< GraphUtilities`"<<endl;
	outfile<<"e = {";

	int numEdges = EdgeList.size();
	int e;

	for(e=0; e<numEdges-1; e++)
	{
		outfile<<""<<EdgeList[e].v0+1<<"->"<<EdgeList[e].v1+1<<",";
	}
	outfile<<""<<EdgeList[e].v0+1<<"->"<<EdgeList[e].v1+1<<"};";

	//assign weights
	outfile<<endl<<"w = {";
	for(e=0; e<numEdges-1; e++)
	{
		outfile<<calcFiberLength(e)<<", ";
	}
	outfile<<calcFiberLength(e)<<"};";

	int numVertices = NodeList.size();
	int v;
	outfile<<endl<<"v = {";
	for(v=0; v<numVertices-1; v++)
	{
		//outfile<<"{{"<<NodeList[v].p.x<<","<<NodeList[v].p.y<<"}}, ";
		outfile<<"{"<<NodeList[v].p.x<<","<<NodeList[v].p.y<<","<<NodeList[v].p.z<<"}, ";
	}
	//outfile<<"{{"<<NodeList[v].p.x<<","<<NodeList[v].p.y<<"}}};";
	outfile<<"{"<<NodeList[v].p.x<<","<<NodeList[v].p.y<<","<<NodeList[v].p.z<<"}};";


	outfile.close();
}

void rtsNetwork::SaveMathematicaNodeDistances(const char* filename)
{
	ofstream outfile;
	outfile.open(filename);

	//for each node, save the distance between all other nodes
	int numNodes = NodeList.size();
	int i, j;
	vector3D<float> distance;
	outfile<<"d = {";
	for(i=0; i<numNodes-1; i++)
	{
		outfile<<"{";
		for(j=0; j<numNodes-1; j++)
		{
			distance = toTissue(NodeList[i].p) - toTissue(NodeList[j].p);
			outfile<<distance.Length()<<",";
		}
		distance = toTissue(NodeList[i].p) - toTissue(NodeList[j].p);
		outfile<<distance.Length()<<"},";
	}
	outfile<<"{";
	for(j=0; j<numNodes-1; j++)
	{
		distance = toTissue(NodeList[i].p) - toTissue(NodeList[j].p);
		outfile<<distance.Length()<<",";
	}
	distance = toTissue(NodeList[i].p) - toTissue(NodeList[j].p);
	outfile<<distance.Length()<<"}};";

	outfile.close();
		

}
void rtsNetwork::RemoveTerminalComponents()
{
	int removed;
	do{
		removed = EdgeList.size();
		//removes fibers that are shorter than the specified length
		NetworkType newNetwork;

		unsigned int numEdges, e;
		numEdges = EdgeList.size();
		unsigned int numVertices, v;

		for(e=0; e<numEdges; e++)
		{
			//if the fiber is not an end fiber (one vertex has valence 1)
			if(getConnections(EdgeList[e].v0).size() != 1 && getConnections(EdgeList[e].v1).size() != 1)
			{
					newNetwork.push_back(network[e]);
					removed--;
				
			}

		}

		network.clear();
		network = newNetwork;

		Clean();
		calcGraph();
	}while(removed != 0);

}

void rtsNetwork::RemoveBoundaryComponents()
{
	int removed;
	do{
		removed = EdgeList.size();
		//removes fibers that are shorter than the specified length
		NetworkType newNetwork;

		unsigned int numEdges, e;
		numEdges = EdgeList.size();
		unsigned int numVertices, v;

		for(e=0; e<numEdges; e++)
		{
			//if the fiber is not an end fiber (one vertex has valence 1)
			if(getConnections(EdgeList[e].v0).size() != 1 && getConnections(EdgeList[e].v1).size() != 1)
			{
					newNetwork.push_back(network[e]);
					removed--;
				
			}
			//if the fiber is an end fiber
			else if(!isBoundaryNode(EdgeList[e].v0) && !isBoundaryNode(EdgeList[e].v1))
			{
				newNetwork.push_back(network[e]);
				removed--;
			}


		}

		network.clear();
		network = newNetwork;

		Clean();
		calcGraph();
	}while(removed != 0);

}
void rtsNetwork::GetRadiusFromVolume(const char *filename, unsigned int threshold, unsigned int range = 2)
{
	//load the volume image
	VOLType::Pointer volume = LoadVOL(filename);
	VOLType::SpacingType spacing;
	// Note: measurement units (e.g., mm, inches, etc.) are defined by the application.
	spacing[0] = 0.86; // spacing along X
	spacing[1] = 1.0; // spacing along Y
	spacing[2] = 1.4; // spacing along Z
	volume->SetSpacing(spacing);
	SaveSlice(volume, 20, "original.bmp");
	//volume->ReleaseDataFlagOn();

/*
	cout<<"Expanding Image..."<<endl;
	typedef itk::ExpandImageFilter<VOLType, VOLType> ResizeFilterType;
	ResizeFilterType::Pointer resizeFilter = ResizeFilterType::New();
	resizeFilter->SetExpandFactors(resample);
	resizeFilter->SetInput(volume);
	resizeFilter->Update();
	VOLType::Pointer resizedImage = resizeFilter->GetOutput();
	SaveSlice(resizedImage, 20, "resized.bmp");
	cout<<"done."<<endl;
*/	

	typedef itk::DiscreteGaussianImageFilter<VOLType, FloatType> BlurFilterType;
	BlurFilterType::Pointer blurFilter = BlurFilterType::New();
	blurFilter->SetUseImageSpacingOn();
	float variance[3] = {0.61, 1.5, 0.61};
	blurFilter->SetVariance(variance);
	blurFilter->SetMaximumKernelWidth(8);
	blurFilter->SetInput(volume);
	try
	{
		blurFilter->Update();
	}
	catch(itk::ExceptionObject & exp)
	{
		std::cout<<exp<<std::endl;
	}
	FloatType::Pointer blurredImage = blurFilter->GetOutput();

	typedef itk::CastImageFilter<FloatType, VOLType> CastingFilterType;
	CastingFilterType::Pointer castFilter = CastingFilterType::New();
	castFilter->SetInput(blurredImage);
	castFilter->Update();
	VOLType::Pointer resizedImage = castFilter->GetOutput();

	SaveSlice(resizedImage, 20, "resized.bmp");


	//threshold the input image
	cout<<"Thresholding Volume..."<<endl;
	typedef itk::BinaryThresholdImageFilter<VOLType, VOLType> ThresholdFilterType;
	ThresholdFilterType::Pointer thresholdFilter = ThresholdFilterType::New();
	thresholdFilter->SetInsideValue(255);
	thresholdFilter->SetOutsideValue(0);
	thresholdFilter->SetLowerThreshold(threshold);
	thresholdFilter->SetUpperThreshold(255);
	thresholdFilter->SetInput(resizedImage);
	try
	{
		thresholdFilter->Update();
	}
	catch(itk::ExceptionObject & exp)
	{
		std::cout<<exp<<std::endl;
	}
	VOLType::Pointer thresholdImage = thresholdFilter->GetOutput();
	//volume->ReleaseData();
	//thresholdFilter->ReleaseDataFlagOn();

	SaveSlice(thresholdImage, 20, "threshold.bmp");

	cout<<"done."<<endl;

	
	cout<<"Computing Distance Field..."<<endl;
	//create the distance transform filter
	typedef itk::DanielssonDistanceMapImageFilter<VOLType, FloatType> DistanceFilterType;
	DistanceFilterType::Pointer distanceFilter = DistanceFilterType::New();
	distanceFilter->SetInput(thresholdImage);
	distanceFilter->UseImageSpacingOn();

	try
	{
		distanceFilter->Update();
	}
	catch(itk::ExceptionObject & exp)
	{
		std::cout<<exp<<std::endl;
	}
	//SaveSlice(distanceFilter->GetOutput(), 20, "distance.bmp");
	FloatType::Pointer distanceMap = distanceFilter->GetOutput();
	//SaveFloatRAW(distanceMap, "distanceMap.raw");
	cout<<"done."<<endl;

	cout<<"Computing Radii..."<<endl;

	

	//create a neighborhood iterator
	typedef itk::ConstNeighborhoodIterator<FloatType> NeighborhoodIteratorType;
	NeighborhoodIteratorType::RadiusType radius;
	radius.Fill(range);
	NeighborhoodIteratorType i(radius, distanceMap, distanceMap->GetBufferedRegion());
	vector<NeighborhoodIteratorType::OffsetType> OffsetVector;
	OffsetVector.resize((range*2+1)*(range*2+1)*(range*2+1));
	int x, y, z, j;
	j=0;
	for(x=-range; x<=range; x++)
		for(y=-range; y<=range; y++)
			for(z=-range; z<=range; z++)
			{
				OffsetVector[j][0]=x;
				OffsetVector[j][1]=y;
				OffsetVector[j][2]=z;
				//cout<<"offset: "<<OffsetVector[j]<<endl;
				j++;
			}


	int numFibers = network.size();
	int f, numVertices, v;
	point3D<float> position;
	float total_radius, near_radius, total_volume, max_radius, min_radius;
	FloatType::IndexType pixelIndex;
	float prev_radius;
	point3D<float> p0, p1;
	float height;
	for(f=0; f<numFibers; f++)
	{
		total_radius = 0;
		total_volume = 0;
		max_radius = 0;
		min_radius = 999;
		numVertices = network[f].size();
		for(v=0; v<numVertices; v++)
		{
			position = network[f][v];
			//pixelIndex[0] = position.x*resample;
			//pixelIndex[1] = position.y*resample;
			//pixelIndex[2] = position.z*resample;
			pixelIndex[0] = position.x;
			pixelIndex[1] = position.y;
			pixelIndex[2] = position.z;
			i.SetLocation(pixelIndex);
			near_radius = 0;
			for(j=0; j<OffsetVector.size(); j++)
			{
				near_radius = max(near_radius, (float)i.GetPixel(OffsetVector[j]));
				//cout<<(float)i.GetPixel(OffsetVector[j])<<endl;
				//cout<<"j = "<<j<<"   "<<i.GetPixel(OffsetVector[j])<<endl;
			}
			//cout<<"v = "<<v<<"   "<<near_radius<<endl;
			if(near_radius > max_radius)
				max_radius = near_radius;
			if(near_radius < min_radius)
				min_radius = near_radius;
			total_radius += near_radius;
			if(v>0)
			{
				//compute the volume of the segment
				p0 = toTissue(network[f][v-1]);
				p1 = toTissue(network[f][v]);
				height = (p1 - p0).Length();
				total_volume += ((3.14159*height)/3.0)*(prev_radius * prev_radius + near_radius*near_radius + prev_radius*near_radius);

			}
			prev_radius = near_radius;
		}
		EdgeList[f].avg_radius = total_radius/numVertices;
		EdgeList[f].min_radius = min_radius;
		EdgeList[f].max_radius = max_radius;
		EdgeList[f].volume = total_volume;
	}

	cout<<"done."<<endl;

}
void rtsNetwork::GetRadiusFromDistanceMap(const char *filename, int range = 2)
{
	//load the volume image
	FloatType::Pointer distanceMap = LoadRAWFloat(filename, 256, 256, 256);

	//create a neighborhood iterator
	typedef itk::ConstNeighborhoodIterator<FloatType> NeighborhoodIteratorType;
	NeighborhoodIteratorType::RadiusType radius;
	radius.Fill(range);
	NeighborhoodIteratorType i(radius, distanceMap, distanceMap->GetRequestedRegion());
	vector<NeighborhoodIteratorType::OffsetType> OffsetVector;
	OffsetVector.resize((range*2+1)*(range*2+1)*(range*2+1));
	int x, y, z, j;
	j=0;
	for(x=-range; x<=range; x++)
		for(y=-range; y<=range; y++)
			for(z=-range; z<=range; z++)
			{
				OffsetVector[j][0]=x;
				OffsetVector[j][1]=y;
				OffsetVector[j][2]=z;
				//cout<<"offset: "<<OffsetVector[j]<<endl;
				j++;
			}

	

	cout<<"Computing Radii..."<<endl;
	int numFibers = network.size();
	int f, numVertices, v;
	point3D<float> position;
	float total_radius, max_radius;
	FloatType::IndexType pixelIndex;
	for(f=0; f<numFibers; f++)
	{
		total_radius = 0;
		numVertices = network[f].size();
		for(v=0; v<numVertices; v++)
		{
			position = network[f][v];
			pixelIndex[0] = position.x;
			pixelIndex[1] = position.y;
			pixelIndex[2] = position.z;
			i.SetLocation(pixelIndex);
			max_radius = 0;
			for(j=0; j<OffsetVector.size(); j++)
			{
				max_radius = max(max_radius, (float)i.GetPixel(OffsetVector[j]));
				//cout<<"j = "<<j<<"   "<<i.GetPixel(OffsetVector[j])<<endl;
			}

			total_radius += max_radius;
		}
		EdgeList[f].avg_radius = total_radius/numVertices;
	}

	cout<<"done."<<endl;

}

point3D<float> rtsNetwork::GetFiberRadius(unsigned int fiber)
{
	point3D<float> result;
	result.x = EdgeList[fiber].min_radius;
	result.y = EdgeList[fiber].max_radius;
	result.z = EdgeList[fiber].avg_radius;
	return result;
}