added function to resample stim vec of points on fiber

pranathivemuri
1 parent cc7e4705
Showing 1 changed file with 58 additions and 15 deletions Show diff stats
stim/biomodels/fiber.h
@@ -18,11 +18,11 @@ protected:
 	double **c;						//centerline (array of double pointers)
 	T* r;						// array of fiber radii
-// 								// fibers this fiber intersects.
 	ANNkd_tree* kdt;			//kd-tree stores all points in the fiber for fast searching
 	/// Initialize an empty fiber
-	void init(){
+	void init()
+	{
 		kdt = NULL;
 		c=NULL;
 		r=NULL;
@@ -30,7 +30,8 @@ protected:
 	}
 	/// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
-	void init(unsigned int n){
+	void init(unsigned int n)
+	{
 		N = n;												//set the number of points
 		kdt = NULL;
@@ -61,27 +62,25 @@ protected:
 		gen_kdtree();							//generate the kd tree for the new fiber
 	}
-	void gen_kdtree(){
-
-		//create an array of data points
-		int n_data = N;
-
+	/// generate a KD tree for points on fiber 
+	void gen_kdtree()
+	{
+		int n_data = N; //create an array of data points
 		ANNpointArray pts = (ANNpointArray)c;			//cast the centerline list to an ANNpointArray
-
 		kdt = new ANNkd_tree(pts, n_data, 3);			//build a KD tree
 	}
-
+	/// find distance between two points
 	double dist(double* p0, double* p1){
-		double sum = 0;
+		double sum = 0; // initialize variables
 		float v;
-		for(unsigned int d = 0; d < 3; d++){
+		for(unsigned int d = 0; d < 3; d++)
+		{
 			v = p1[d] - p0[d];
 			sum +=v * v;
 		}
-
 		return sqrt(sum);
 	}
@@ -112,8 +111,6 @@ protected:
 	}
-
-
 public:
 	fiber(){
@@ -490,6 +487,52 @@ public:
 	stim::vec<T> back(){
 		return get_vec(N-1);
 	}
+		////resample a fiber in the network
+	std::vector<stim::vec<T> > Resample(std::vector< stim::vec<T> > fiberPositions, float spacing=25.0)
+	{
+
+		std::vector<T> v(3);    //v-direction vector of the segment
+		stim::vec<T> p(3);      //- intermediate point to be added
+		std::vector<T> p1(3);   // p1 - starting point of an segment on the fiber, 
+		std::vector<T> p2(3);   // p2 - ending point, 
+		double sum=0;  //distance summation
+
+		std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
+		// for each point on the centerline (skip if it is the last point on centerline)
+		unsigned int N = fiberPositions.size(); // number of points on the fiber
+		for(unsigned int f=0; f< N-1; f++)
+		{
+			for(unsigned int d = 0; d < 3; d++)
+			{
+				p1[d] = fiberPositions[f][d]; // starting point of the segment on the centerline
+				p2[d] = fiberPositions[f + 1][d]; // ending point of the segment on the centerline
+				v[d] = p2[d] - p1[d];            //direction vector
+				sum +=v[d] * v[d];               //length of segment-distance between starting and ending point
+			}
+			newPointList.push_back(fiberPositions[f]); //always push the starting point
+
+			T lengthSegment = sqrt(sum);  //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
+						for(unsigned int i=0; i<3;i++)
+							{
+								p[i] = p1[i] + v[i]*(step/lengthSegment);
+							} //for each dimension
+						// add this resampled points to the new fiber list
+						newPointList.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
+				newPointList.push_back(fiberPositions[f+1]);
+			}
+		newPointList.push_back(fiberPositions[N-1]);   //add the last point on the fiber to the new fiber list
+		return newPointList;
+	}
 };