segmentation / netmets

#ifndef STIM_CENTERLINE_H
#define STIM_CENTERLINE_H

#include <vector>
#include <stim/math/vec3.h>
#include <stim/structures/kdtree.cuh>

namespace stim{

/**	This class stores information about a single fiber represented as a set of geometric points
 *	between two branch or end points. This class is used as a fundamental component of the stim::network
 *	class to describe an interconnected (often biological) network.
 */
template<typename T>
class centerline : public std::vector< stim::vec3<T> >{

protected:

	std::vector<T> L;										//stores the integrated length along the fiber (used for parameterization)

	///Return the normalized direction vector at point i (average of the incoming and outgoing directions)
	vec3<T> d(size_t i) {
		if (size() <= 1) return vec3<T>(0, 0, 0);						//if there is insufficient information to calculate the direction, return a null vector
		if (size() == 2) return (at(1) - at(0)).norm();					//if there are only two points, the direction vector at both is the direction of the line segment
		if (i == 0) return (at(1) - at(0)).norm();						//the first direction vector is oriented towards the first line segment
		if (i == size() - 1) return (at(size() - 1) - at(size() - 2)).norm();	//the last direction vector is oriented towards the last line segment

		//all other direction vectors are the average direction of the two joined line segments
		vec3<T> a = at(i) - at(i - 1);
		vec3<T> b = at(i + 1) - at(i);
		vec3<T> ab = a.norm() + b.norm();
		return ab.norm();
	}

	//initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
	void update_L(size_t start = 0) {
		L.resize(size());									//allocate space for the L array
		if (start == 0) {
			L[0] = 0;											//initialize the length value for the first point to zero (0)
			start++;
		}

		stim::vec3<T> d;
		for (size_t i = start; i < size(); i++) {		//for each line segment in the centerline
			d = at(i) - at(i - 1);
			L[i] = L[i - 1] + d.len();				//calculate the running length total
		}
	}

	void init() {
		if (size() == 0) return;								//return if there aren't any points
		update_L();
	}

	/// Returns a stim::vec representing the point at index i

	/// @param i is an index of the desired centerline point
	stim::vec<T> get_vec(unsigned i){
		return std::vector< stim::vec3<T> >::at(i);
	}

	///finds the index of the point closest to the length l on the lower bound.
	///binary search.
	size_t findIdx(T l) {
		for (size_t i = 1; i < L.size(); i++) {				//for each point in the centerline
			if (L[i] > l) return i - 1;						//if we have passed the desired length value, return i-1
		}
		return L.size() - 1;
		/*size_t i = L.size() / 2;
		size_t max = L.size() - 1;
		size_t min = 0;
		while (i < L.size() - 1){
			if (l < L[i]) {
				max = i;
				i = min + (max - min) / 2;
			}
			else if (L[i] <= l && L[i + 1] >= l) {
				break;
			}
			else {
				min = i;
				i = min + (max - min) / 2;
			}
		}
		return i;*/
	}

	///Returns a position vector at the given length into the fiber (based on the pvalue).
	///Interpolates the radius along the line.
	///@param l: the location of the in the cylinder.
	///@param idx: integer location of the point closest to l but prior to it.
	stim::vec3<T> p(T l, int idx) {
		T rat = (l - L[idx]) / (L[idx + 1] - L[idx]);
		stim::vec3<T> v1 = at(idx);
		stim::vec3<T> v2 = at(idx + 1);
		return(v1 + (v2 - v1)*rat);
	}


public:

	using std::vector< stim::vec3<T> >::at;
	using std::vector< stim::vec3<T> >::size;

	centerline() : std::vector< stim::vec3<T> >() {
		init();
	}
	centerline(size_t n) : std::vector< stim::vec3<T> >(n){
		init();
	}
	centerline(std::vector<stim::vec3<T> > pos) :
		std::vector<stim::vec3<T> > (pos)
	{
		init();
	}
	
	//overload the push_back function to update the length vector
	void push_back(stim::vec3<T> p) {
		std::vector< stim::vec3<T> >::push_back(p);
		update_L(size() - 1);
	}

	///Returns a position vector at the given p-value (p value ranges from 0 to 1).
	///interpolates the position along the line.
	///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
	stim::vec3<T> p(T pvalue) {
		if (pvalue <= 0.0) return at(0);			//return the first element
		if (pvalue >= 1.0) return back();			//return the last element

		T l = pvalue*L[L.size() - 1];
		int idx = findIdx(l);
		return p(l, idx);
	}

	///Update centerline internal parameters (currently the L vector)
	void update() {
		init();
	}
	///Return the length of the entire centerline
	T length() {
		return L.back();
	}


	/// stitch two centerlines
	///@param c1, c2: two centerlines
	///@param sigma: sample rate
	static std::vector< stim::centerline<T> > stitch(stim::centerline<T> c1, stim::centerline<T> c2 = stim::centerline<T>()) {
		
		std::vector< stim::centerline<T> > result;
		stim::centerline<T> new_centerline;
		stim::vec3<T> new_vertex;

		// if only one centerline, stitch itself!
		if (c2.size() == 0) {
			size_t num = c1.size();
			size_t id = 100000;							// store the downsample start position
			T threshold;
			if (num < 4) {								// if the number of vertex is less than 4, do nothing
				result.push_back(c1);
				return result;
			}
			else {
				// test geometry start vertex
				stim::vec3<T> v1 = c1[1] - c1[0];		// vector from c1[0] to c1[1]
				for (size_t p = 2; p < num; p++) {		// 90° standard???
					stim::vec3<T> v2 = c1[p] - c1[0];
					float cosine = v2.dot(v1);
					if (cosine < 0) {
						id = p;
						threshold = v2.len();
						break;
					}
				}
				if (id != 100000) {						// find a downsample position on the centerline
					T* c;
					c = (T*)malloc(sizeof(T) * (num - id) * 3);
					for (size_t p = id; p < num; p++) {
						for (size_t d = 0; d < 3; d++) {
							c[p * 3 + d] = c1[p][d];
						}
					}
					stim::kdtree<T, 3> kdt;
					kdt.create(c, num - id, 5);			// create tree

					T* query = (T*)malloc(sizeof(T) * 3);
					for (size_t d = 0; d < 3; d++)
						query[d] = c1[0][d];
					size_t index;
					T dist;

					kdt.search(query, 1, &index, &dist);

					free(query);
					free(c);

					if (dist > threshold) {
						result.push_back(c1);
					}
					else {
						// the loop part
						new_vertex = c1[index];
						new_centerline.push_back(new_vertex);
						for (size_t p = 0; p < index + 1; p++) {
							new_vertex = c1[p];
							new_centerline.push_back(new_vertex);
						}
						result.push_back(new_centerline);
						new_centerline.clear();

						// the tail part
						for (size_t p = index; p < num; p++) {
							new_vertex = c1[p];
							new_centerline.push_back(new_vertex);
						}
						result.push_back(new_centerline);
					}
				}
				else {	// there is one potential problem that two positions have to be stitched
						// test geometry end vertex
					stim::vec3<T> v1 = c1[num - 2] - c1[num - 1];
					for (size_t p = num - 2; p > 0; p--) {		// 90° standard
						stim::vec3<T> v2 = c1[p - 1] - c1[num - 1];
						float cosine = v2.dot(v1);
						if (cosine < 0) {
							id = p;
							threshold = v2.len();
							break;
						}
					}
					if (id != 100000) {						// find a downsample position
						T* c;
						c = (T*)malloc(sizeof(T) * (id + 1) * 3);
						for (size_t p = 0; p < id + 1; p++) {
							for (size_t d = 0; d < 3; d++) {
								c[p * 3 + d] = c1[p][d];
							}
						}
						stim::kdtree<T, 3> kdt;
						kdt.create(c, id + 1, 5);				// create tree

						T* query = (T*)malloc(sizeof(T) * 1 * 3);
						for (size_t d = 0; d < 3; d++)
							query[d] = c1[num - 1][d];
						size_t index;
						T dist;

						kdt.search(query, 1, &index, &dist);

						free(query);
						free(c);

						if (dist > threshold) {
							result.push_back(c1);
						}
						else {
							// the tail part
							for (size_t p = 0; p < index + 1; p++) {
								new_vertex = c1[p];
								new_centerline.push_back(new_vertex);
							}
							result.push_back(new_centerline);
							new_centerline.clear();

							// the loop part
							for (size_t p = index; p < num; p++) {
								new_vertex = c1[p];
								new_centerline.push_back(new_vertex);
							}
							new_vertex = c1[index];
							new_centerline.push_back(new_vertex);
							result.push_back(new_centerline);
						}
					}
					else {	// no stitch position
						result.push_back(c1);
					}
				}
			}
		}


		// two centerlines
		else {
			// find stitch position based on nearest neighbors												
			size_t num1 = c1.size();
			T* c = (T*)malloc(sizeof(T) * num1 * 3);		// c1 as reference point
			for (size_t p = 0; p < num1; p++)				// centerline to array
				for (size_t d = 0; d < 3; d++)				// because right now my kdtree code is a relatively close code, it has its own structure
					c[p * 3 + d] = c1[p][d];				// I will merge it into stimlib totally in the near future

			stim::kdtree<T, 3> kdt;							// kdtree object
			kdt.create(c, num1, 5);							// create tree

			size_t num2 = c2.size();						
			T* query = (T*)malloc(sizeof(T) * num2 * 3);	// c2 as query point
			for (size_t p = 0; p < num2; p++) {
				for (size_t d = 0; d < 3; d++) {
					query[p * 3 + d] = c2[p][d];
				}
			}
			std::vector<size_t> index(num2);
			std::vector<T> dist(num2);

			kdt.search(query, num2, &index[0], &dist[0]);	// find the nearest neighbors in c1 for c2

			// clear up
			free(query);
			free(c);

			// find the average vertex distance of one centerline
			T sigma1 = 0;
			T sigma2 = 0;
			for (size_t p = 0; p < num1 - 1; p++) 
				sigma1 += (c1[p] - c1[p + 1]).len();
			for (size_t p = 0; p < num2 - 1; p++)
				sigma2 += (c2[p] - c2[p + 1]).len();
			sigma1 /= (num1 - 1);
			sigma2 /= (num2 - 1);
			float threshold = 4 * (sigma1 + sigma2) / 2;			// better way to do this?

			T min_d = *std::min_element(dist.begin(), dist.end());	// find the minimum distance between c1 and c2

			if (min_d > threshold) {								// if the minimum distance is too large
				result.push_back(c1);
				result.push_back(c2);

#ifdef DEBUG
				std::cout << "The distance between these two centerlines is too large" << std::endl;
#endif
			}
			else {
			//	auto smallest = std::min_element(dist.begin(), dist.end());
				unsigned int smallest = std::min_element(dist.begin(), dist.end());
			//	auto i = std::distance(dist.begin(), smallest);		// find the index of min-distance in distance list
				unsigned int i = std::distance(dist.begin(), smallest);		// find the index of min-distance in distance list

				// stitch position in c1 and c2
				int id1 = index[i];
				int id2 = i;

				// actually there are two cases
				// first one inacceptable
				// second one acceptable
				if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) {		// only stitch one end vertex to another centerline
					result.push_back(c1);
					result.push_back(c2);
				}
				else {
					if (id1 == 0 || id1 == num1 - 1) {			// if the stitch vertex is the first or last vertex of c1
						// for c2, consider two cases(one degenerate case)
						if (id2 == 0 || id2 == num2 - 1) {		// case 1, if stitch position is also on the end of c2
							// we have to decide which centerline get a new vertex, based on direction
							// for c1, computer the direction change angle
							stim::vec3<T> v1, v2;
							float alpha1, alpha2;				// direction change angle
							if (id1 == 0)
								v1 = (c1[1] - c1[0]).norm();
							else
								v1 = (c1[num1 - 2] - c1[num1 - 1]).norm();
							v2 = (c2[id2] - c1[id1]).norm();
							alpha1 = v1.dot(v2);
							if (id2 == 0)
								v1 = (c2[1] - c2[0]).norm();
							else
								v1 = (c2[num2 - 2] - c2[num2 - 1]).norm();
							v2 = (c1[id1] - c2[id2]).norm();
							alpha2 = v1.dot(v2);
							if (abs(alpha1) > abs(alpha2)) {					// add the vertex to c1 in order to get smooth connection
								// push back c1
								if (id1 == 0) {									// keep geometry information
									new_vertex = c2[id2];
									new_centerline.push_back(new_vertex);
									for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
										new_vertex = c1[p];
										new_centerline.push_back(new_vertex);
									}
								}
								else {
									for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
										new_vertex = c1[p];
										new_centerline.push_back(new_vertex);
									}
									new_vertex = c2[id2];
									new_centerline.push_back(new_vertex);
								}
								result.push_back(new_centerline);
								new_centerline.clear();

								// push back c2
								for (size_t p = 0; p < num2; p++) {
									new_vertex = c2[p];
									new_centerline.push_back(new_vertex);
								}
								result.push_back(new_centerline);
							}
							else {												// add the vertex to c2 in order to get smooth connection
								// push back c1
								for (size_t p = 0; p < num1; p++) {
									new_vertex = c1[p];
									new_centerline.push_back(new_vertex);
								}
								result.push_back(new_centerline);
								new_centerline.clear();

								// push back c2
								if (id2 == 0) {									// keep geometry information
									new_vertex = c1[id1];
									new_centerline.push_back(new_vertex);
									for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
										new_vertex = c2[p];
										new_centerline.push_back(new_vertex);
									}
								}
								else {
									for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
										new_vertex = c2[p];
										new_centerline.push_back(new_vertex);
									}
									new_vertex = c1[id1];
									new_centerline.push_back(new_vertex);
								}
								result.push_back(new_centerline);
							}
						}
						else {												// case 2, the stitch position is on c2
							// push back c1
							if (id1 == 0) {									// keep geometry information
								new_vertex = c2[id2];
								new_centerline.push_back(new_vertex);
								for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
									new_vertex = c1[p];
									new_centerline.push_back(new_vertex);
								}
							}
							else {
								for (size_t p = 0; p < num1; p++) {			// geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2
									new_vertex = c1[p];
									new_centerline.push_back(new_vertex);
								}
								new_vertex = c2[id2];
								new_centerline.push_back(new_vertex);
							}
							result.push_back(new_centerline);
							new_centerline.clear();

							// push back c2
							for (size_t p = 0; p < id2 + 1; p++) {			// first part
								new_vertex = c2[p];
								new_centerline.push_back(new_vertex);
							}
							result.push_back(new_centerline);
							new_centerline.clear();

							for (size_t p = id2; p < num2; p++) {			// second part
								new_vertex = c2[p];
								new_centerline.push_back(new_vertex);
							}
							result.push_back(new_centerline);
						}
					}
					else {							// if the stitch vertex is the first or last vertex of c2
						// push back c2
						if (id2 == 0) {										// keep geometry information
							new_vertex = c1[id1];
							new_centerline.push_back(new_vertex);
							for (size_t p = 0; p < num2; p++) {				// stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2
								new_vertex = c2[p];
								new_centerline.push_back(new_vertex);
							}
						}
						else {
							for (size_t p = 0; p < num2; p++) {				// geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1
								new_vertex = c2[p];
								new_centerline.push_back(new_vertex);
							}
							new_vertex = c1[id1];
							new_centerline.push_back(new_vertex);
							result.push_back(new_centerline);
							new_centerline.clear();

							// push back c1
							for (size_t p = 0; p < id1 + 1; p++) {			// first part
								new_vertex = c1[p];
								new_centerline.push_back(new_vertex);
							}
							result.push_back(new_centerline);
							new_centerline.clear();

							for (size_t p = id1; p < num1; p++) {			// second part
								new_vertex = c1[p];
								new_centerline.push_back(new_vertex);
							}
							result.push_back(new_centerline);
						}
					}
				}
			}
		}
		return result;
	}

	/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
	std::vector< stim::centerline<T> > split(unsigned int idx){

		std::vector< stim::centerline<T> > fl;				//create an array to store up to two fibers
		size_t N = size();

		//if the index is an end point, only the existing fiber is returned
		if(idx == 0 || idx == N-1){
			fl.resize(1);							//set the size of the fiber to 1
			fl[0] = *this;							//copy the current fiber
		}

		//if the index is not an end point
		else{

			unsigned int N1 = idx + 1;					//calculate the size of both fibers
			unsigned int N2 = N - idx;

			fl.resize(2);								//set the array size to 2

			fl[0] = stim::centerline<T>(N1);			//set the size of each fiber
			fl[1] = stim::centerline<T>(N2);

			//copy both halves of the fiber
			unsigned int i;

			//first half
			for(i = 0; i < N1; i++)					//for each centerline point
				fl[0][i] = std::vector< stim::vec3<T> >::at(i);
			fl[0].init();							//initialize the length vector

			//second half
			for(i = 0; i < N2; i++)
				fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
			fl[1].init();							//initialize the length vector
		}

		return fl;										//return the array

	}

	/// Outputs the fiber as a string
	std::string str(){
		std::stringstream ss;
		size_t N = std::vector< stim::vec3<T> >::size();
		ss << "---------[" << N << "]---------" << std::endl;
		for (size_t i = 0; i < N; i++)
			ss << std::vector< stim::vec3<T> >::at(i) << std::endl;
		ss << "--------------------" << std::endl;

		return ss.str();
	}

	/// Back method returns the last point in the fiber
	stim::vec3<T> back(){
		return std::vector< stim::vec3<T> >::back();
	}

		////resample a fiber in the network
	stim::centerline<T> resample(T spacing)
	{	
		//std::cout<<"fiber::resample()"<<std::endl;

		stim::vec3<T> v;    //v-direction vector of the segment
		stim::vec3<T> p;      //- intermediate point to be added
		stim::vec3<T> p1;   // p1 - starting point of an segment on the fiber,
		stim::vec3<T> p2;   // p2 - ending point,
		//double sum=0;  //distance summation

		size_t N = size();

		centerline<T> new_c; // initialize list of new resampled points on the fiber
		// for each point on the centerline (skip if it is the last point on centerline)
		for(unsigned int f=0; f< N-1; f++)
		{			
			p1 = at(f); 
			p2 = at(f+1);
			v = p2 - p1;
			
			T lengthSegment = v.len();			//find Length of the segment as distance between the starting and ending points of the segment

			if(lengthSegment >= spacing){ // if length of the segment is greater than standard deviation resample
				
				// repeat resampling until accumulated stepsize is equsl to length of the segment
				for(T step=0.0; step<lengthSegment; step+=spacing){
					// calculate the resampled point by travelling step size in the direction of normalized gradient vector
					p = p1 + v * (step / lengthSegment);
					
					// add this resampled points to the new fiber list
					new_c.push_back(p);
				}
			}
			else       // length of the segment is now less than standard deviation, push the ending point of the segment and proceed to the next point in the fiber
				new_c.push_back(at(f));
		}
		new_c.push_back(at(N-1));   //add the last point on the fiber to the new fiber list
		//centerline newFiber(newPointList);
		return new_c;
	}

};



}	//end namespace stim



#endif