fixed the issue with the seed points, fixed the issue with the no-terminating st…

…atement, this commit is still using ivote due to instability with the LaG version

fixed the issue with the seed points, fixed the issue with the no-terminating st…
…atement, this commit is still using ivote due to instability with the LaG version
Pavel Govyadinov
1 parent 26aee9ed
Showing 3 changed files with 147 additions and 73 deletions Show diff stats
stim/cuda/branch_detection.cuh
stim/gl/gl_spider.h
stim/visualization/cylinder.h
@@ -101,7 +101,7 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 	unsigned int 	x_ds		= (x + (x % 1 == 0 ? 0:1));
 	unsigned int 	y_ds		= (y + (x % 1 == 0 ? 0:1));
 	unsigned int	bytes_ds	= sizeof(float) * x_ds * y_ds;
-	unsigned int	conn		= 5;
+	unsigned int	conn		= 10;
 	float		final_t		= 200.0;
 	float*		cpuTable	= (float*) malloc(bytes_table);
 	float*		cpuCenters	= (float*) malloc(bytes_ds);
@@ -111,7 +111,7 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
-	test(texbufferID, texType, x, y);
+//	test(texbufferID, texType, x, y);
 	std::vector<stim::vec<float> >  output;
 	initCuda(bytes_table, bytes_ds); 
@@ -124,7 +124,7 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 	stim::cuda::tex_gaussian_blur2<float>(
 		gpuI, sigma, x, y, t.getTexture(), t.getArray()
 		);
-	stim::gpu2image<float>(gpuI, "Blur.jpg", x,y , 0, 255);
+//	stim::gpu2image<float>(gpuI, "Blur.jpg", x,y , 0, 255);
 //	stim::gpu2image<float>(t.getArray(), "ORIGINAL.jpg", x,y , 0, 255);
 	cudaDeviceSynchronize();
@@ -149,24 +149,50 @@ find_branch(GLint texbufferID, GLenum texType, unsigned int x, unsigned int y)
 	cudaDeviceSynchronize();
 	stim::cuda::gpu_local_max<float>(gpuCenters, gpuVote, final_t, conn, x, y);
-	stim::gpu2image<float>(gpuCenters, "Centers.jpg", x, y, 0, 1);
+//	stim::gpu2image<float>(gpuCenters, "Centers.jpg", x, y, 0, 1);
 	cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost);
-	stim::cpu2image<float>(cpuCenters, "CentersXPU.jpg", x, y, 0, 1);
-	std::cout << pixels << " " << x << " " << y << std::endl;
-	for(int i = 0; i < pixels; i++)
+//	stim::cpu2image<float>(cpuCenters, "CentersXPU.jpg", x, y, 0, 1);
+//	std::cerr << pixels << " " << x << " " << y << std::endl;
+//	std::cerr << "y is " << y << ", x is " << x << std::endl; 
+
+//	std::cout << "Before " << output.size() << std::endl;
+
+	for(int i = 0; i < x; i++)
+	{
+		for(int j = 0; j < y; j++)
+		{
+			int idx = x*j+i;
+			if(cpuCenters[idx] != 0)
+			{
+				float x_v = (float) i;
+				float y_v = (float) j;
+				std::cout << x_v/x*360.0 << std::endl;
+				std::cout << y_v/y << std::endl;
+				output.push_back(stim::vec<float>((x_v/(float)x*360.0),
+								  (y_v), y_v/8));	
+			}
+
+		} 
+	}
+	
+/*	for(int i = 0; i < pixels; i++)
 	{
 		int ix = (i % x);
 		int iy = (i / x);
-		if((cpuCenters[i] == 1))
+		if(cpuCenters[i] != 0)
 		{
 			float x_v = (float) ix;
 			float y_v = (float) iy;
+			std::cout << x_v/x*360 << std::endl;
+			std::cout << y_v/y << std::endl;
 			output.push_back(stim::vec<float>((x_v/(float)x*360),
 							  (y_v/(float)y), 0.0));	
 		}
-	}
+	} */
+
+//	std::cout << "After " <<  output.size() << std::endl;
 	t.UnmapCudaTexture();
@@ -229,24 +229,49 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			DrawLongCylinder(n, l_template, l_square);
 			stim::cylinder<float> cyl(cL, cM);
 			std::vector< stim::vec<float> > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template);
+			stim::vec<float> size(S[0]*R[0], S[1]*R[1], S[2]*R[2]);
+			float pval;
 //			std::cerr << "the number of points is " << result.size() << std::endl;
 			if(!result.empty())
 			{
 				for(int i = 0; i < result.size(); i++)
 				{
-					std::cout << "Testing "<< i << ": " << result[i][0] << ", " << result[i][1] << std::endl;
-					stim::vec<float> v = cyl.surf(result[i][1], result[i][0]);
+					int id = result[i][2];
+					if(fmod(result[i][2], id) != 0 && id != 0)
+					{
+						
+						pval = ((cyl.getl(id+1)-cyl.getl(id))*
+							(fmod(result[i][2], id))+cyl.getl(id))/cyl.getl(cL.size()-1);
+//						std::cout << id << " " <<  cyl.getl(id) << " " << pval << " " << cyl.getl(cL.size()-1) << fmod(result[i][2], id) << std::endl;
+					}
+					else if(id == 0)
+					{
+						pval = (cyl.getl(id+1)*result[i][2])/cyl.getl(cL.size()-1);
+					}
+					else
+					{
+						pval = (cyl.getl(id)/cyl.getl(cL.size()-1));
+					}
+//					std::cout << "Testing "<< i << ": " << result[i][0] << ", " << result[i][1] << ", " <<  result[i][2] <<  std::endl;
+//					std::cout << "Testing " << pval << std::endl;
+					stim::vec<float> v = cyl.surf(pval, result[i][0]);
 //					std::cout <<  v[0] << " ," << v[1] << " ," << v[2] << std::endl;
-					stim::vec<float> di = cyl.p(result[i][1]);
-					std::cout <<  di[0] <<  " ," << di[1] << " ," << di[2] << std::endl;
-					float rad = cyl.r(result[i][1]);
-//					std::cout << rad << std::endl;
-					setSeed(v);
-					setSeedVec(stim::vec<float>(-di[0]+v[0], -di[1]+v[1], -di[2]+v[2]).norm());
-					setSeedMag(cyl.r(result[i][1]));
+					stim::vec<float> di = cyl.p(pval);
+//					std::cout <<  di[0] <<  " ," << di[1] << " ," << di[2] << std::endl;
+					float rad = cyl.r(pval);
+					std::cout << rad << std::endl;
+					if(
+					 !(v[0] > size[0] || v[1] > size[1]
+					 || v[2] > size[2] || v[0] < 0
+					 || v[1] < 0 || v[2] < 0))
+					{
+						setSeed(v);
+						setSeedVec((v-di).norm());
+						setSeedMag(rad);
+					}
 				}
 			}
-			std::cout << "I ran the new branch detection" << std::endl;
+//			std::cout << "I ran the new branch detection" << std::endl;
 			}
 		}
@@ -1264,6 +1289,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			 glEndList();
 		}
+///need to return the cylinder.
 		void
 		DrawLongCylinder(int n = 8, int l_template = 8,int l_square = 8)
 		{
@@ -1590,6 +1616,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 			{
 				int cost = Step();
 				if (cost > min_cost){
+					std::cout << "Cost Limit" << std::endl;
 					running = false;
 					sk.End();
 		branchDetection2();
@@ -1604,6 +1631,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					 || pos[2] > size[2] || pos[0] < 0
 					 || pos[1] < 0 || pos[2] < 0)
 					{
+						std::cout << "Edge Limit" << std::endl;
 //					       std::cout << "Found Edge" << std::endl;
 						running = false;
 						sk.End();
@@ -1624,7 +1652,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 					//Has the template size gotten unreasonable?
 					mag = getMagnitude();
 					if(mag[0] > 75 || mag[0] < 1){
-//						std::cout << "Magnitude Limit" << std::endl;
+						std::cout << "Magnitude Limit" << std::endl;
 						running = false;
 						sk.End();
 		branchDetection2();
@@ -1638,6 +1666,7 @@ class gl_spider : public virtual gl_texture&lt;T&gt;
 						h = selectObject(p, getDirection(), m[0]);
 						//Have we hit something previously traced?
 						if(h != -1){
+							std::cout << "Hit Limit" << std::endl;
 							running = false;
 							sk.End();
 		branchDetection2();
@@ -12,8 +12,7 @@ class cylinder
 {
 	private:
 		stim::circle<T> s;			//an arbitrary circle
-		std::vector< stim::vec<T> > pos;	//positions of the cylinder.
-		std::vector< stim::vec<T> > mags;	//radii at each position
+		std::vector<stim::circle<T> > e;
 		std::vector< T > L;			//length of the cylinder at each position.
 		///default init	
@@ -27,18 +26,44 @@ class cylinder
 		void
 		init(std::vector<stim::vec<T> > inP, std::vector<stim::vec<T> > inM)
 		{
-			pos = inP;
-			mags = inM;
+//			pos = inP;
+//			mags = inM;
+			stim::vec<float> v1;
+			stim::vec<float> v2;
+			e.resize(inP.size());
+			if(inP.size() < 2)
+				return;
 			//calculate each L.
-			L.resize(pos.size());
+			L.resize(inP.size());
 			T temp = (T)0;
 			L[0] = 0;
 			for(int i = 1; i < L.size(); i++)
 			{
-				temp += (pos[i-1] - pos[i]).len();
+				temp += (inP[i-1] - inP[i]).len();
 				L[i] = temp;
 			}
+
+			stim::vec<T> dr = (inP[1] - inP[0]).norm();
+			s = stim::circle<T>(inP[0], inM[0][0], dr, stim::vec<T>(1,0,0));
+			e[0] = s;
+			for(int i = 1; i < inP.size()-1; i++)
+			{
+				s.center(inP[i]);
+				v1 = (inP[i] - inP[i-1]).norm();
+				v2 = (inP[i+1] - inP[i]).norm();
+				dr = (v1+v2).norm();
+				s.normal(dr);
+				s.scale(inM[i][0]/inM[i-1][0]);
+				e[i] = s;
+			}
+			
+			int j = inP.size()-1;
+			s.center(inP[j]);
+			dr = (inP[j] - inP[j-1]).norm();
+			s.normal(dr);
+			s.scale(inM[j][0]/inM[j-1][0]);
+			e[j] = s;
 		}
 		///returns the direction vector at point idx.
@@ -47,33 +72,38 @@ class cylinder
 		{
 			if(idx == 0)
 			{
-				return (pos[idx+1] - pos[idx]).norm();
+				return (e[idx+1].P - e[idx].P).norm();
 			}
-			else if(idx == pos.size()-1)
+			else if(idx == e.size()-1)
 			{
-				return (pos[idx] - pos[idx-1]).norm();
+				return (e[idx].P - e[idx-1].P).norm();
 			}
 			else
 			{
-				stim::vec<float> v1 = (pos[idx]-pos[idx-1]).norm();
-				stim::vec<float> v2 = (pos[idx+1]-pos[idx]).norm();
+//				return (e[idx+1].P - e[idx].P).norm();
+				stim::vec<float> v1 = (e[idx].P-e[idx-1].P).norm();
+				stim::vec<float> v2 = (e[idx+1].P-e[idx].P).norm();
 				return (v1+v2).norm();			
-			}
-			
+			} 
+	//		return e[idx].N;	
+
 		}
-		///returns the total length of the line at index j.
-		T
-		getl(int j)
+		stim::vec<T>
+		d(T l, int idx)
 		{
-			T temp = (T) 0;
-			for(int i = 0; i < j; ++i)
+			if(idx == 0 || idx == e.size()-1)
 			{
-				temp += (pos[i] - pos[i+1]).len();
-				L[i] = temp;
+				return e[idx].N;
 			}
+			else
+			{
+				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
+				return(	e[idx].N + (e[idx+1].N - e[idx].N)*rat);
+			} 
 		}
+
 		///finds the index of the point closest to the length l on the lower bound.
 		///binary search.
 		int
@@ -82,26 +112,22 @@ class cylinder
 			unsigned int i = L.size()/2;
 			unsigned int max = L.size()-1;
 			unsigned int min = 0;
-//			std::cerr << "Index initially: " << i << std::endl;
-//			std::cerr << "l initially: " << l << std::endl;
 			while(i > 0 && i < L.size()-1)
 			{
-//				std::cerr << "L[i] =  " << L[i] << " i= "  << i << " maxval= " << L.size()-1 << std::endl;
+//				std::cerr << "Trying " << i << std::endl;
+//				std::cerr << "l is " << l << ", L[" << i << "]" << L[i] << std::endl;
 				if(l < L[i])
 				{
 					max = i;
-//				std::cerr << l << " < " << L[i] << "so max is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl;
 					i = min+(max-min)/2;
 				}
 				else if(L[i] <= l && L[i+1] >= l)
 				{
-//				std::cerr << L[i] << " < " << l << " < " << L[i+1] << std::endl;
 					break;
 				}
 				else
 				{
 					min = i;
-//				std::cerr << l << " > " << L[i] << "so min is set to " << i << " and the new i is " << min+(max-min)/2 << std::endl;
 					i = min+(max-min)/2;
 				}
 			}
@@ -137,7 +163,7 @@ class cylinder
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
 				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-			return(	pos[idx] + (pos[idx+1]-pos[idx])*rat);
+			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
 		}
 		///Returns a position vector at the given length into the fiber (based on the pvalue).
@@ -148,7 +174,7 @@ class cylinder
 		p(T l, int idx)
 		{
 				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-			return(	pos[idx] + (pos[idx+1]-pos[idx])*rat);
+			return(	e[idx].P + (e[idx+1].P-e[idx].P)*rat);
 //			return(
 //			return (pos[idx] + (pos[idx+1]-pos[idx])*((l-L[idx])/(L[idx+1]- L[idx])));
 		}
@@ -163,7 +189,7 @@ class cylinder
 				return;
 			T l = pvalue*L[L.size()-1];
 			int idx = findIdx(l);
-			return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0];
+			return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*((l-L[idx])/(L[idx+1]- L[idx])));
 		}
 		///Returns a radius at the given length into the fiber (based on the pvalue).
@@ -174,7 +200,7 @@ class cylinder
 		r(T l, int idx)
 		{
 				T rat = (l-L[idx])/(L[idx+1]-L[idx]);
-			return(	mags[idx] + (mags[idx+1]-mags[idx])*rat)[0];
+			return(	e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat);
 //			return (mags[idx] + (mags[idx+1]-mags[idx])*((l-L[idx])/(L[idx+1]- L[idx])))[0];
 		}
@@ -197,8 +223,10 @@ class cylinder
 //			std::cerr << "Index: " << idx << std::endl;
 			stim::vec<T> ps = p(l, idx); 
 			T m = r(l, idx);
-			stim::vec<T> dr = d(idx);
-			s = stim::circle<T>(ps, m, dr, stim::vec<T>(1,0,0));
+			s = e[idx];
+			s.center(ps);
+			s.normal(d(l, idx));
+			s.scale(m/e[idx].U.len());
 			return(s.p(theta));
 			}
 		}
@@ -208,29 +236,13 @@ class cylinder
 		std::vector<std::vector<vec<T> > >
 		getPoints(int sides)
 		{
-			if(pos.size() < 2)
+			std::vector<std::vector <vec<T> > > points;
+			points.resize(e.size());
+			for(int i = 0; i < e.size(); i++)
 			{
-				return;
-			} else {
-				std::vector<std::vector <vec<T> > > points;
-				points.resize(pos.size());
-				stim::vec<T> dr = (pos[1] - pos[0]).norm();
-				s = stim::circle<T>(pos[0], mags[0][0], dr, stim::vec<T>(1,0,0));
-				points[0] = s.getPoints(sides);
-				for(int i = 1; i < pos.size(); i++)
-				{
-					dr = d(i);
-//					dr = (pos[i-1] - pos[i]).norm();
-					
-//					s = stim::circle<T>(pos[i], mags[i][0], dr, stim::vec<T>(1,0,0));
-					s.center(pos[i]);
-					s.normal(dr);
-//					s.scale(mags[i][0]/mags[i-1][0], mags[i][0]/mags[i-1][0]);
-					s.scale(mags[i][0]/mags[i-1][0]);
-					points[i] = s.getPoints(sides);
-				}
-				return points;
+				points[i] = e[i].getPoints(sides);
 			}
+			return points;
 		}
 		void
@@ -238,6 +250,13 @@ class cylinder
 		{
 			std::cout << d(idx) << std::endl;
 		}
+
+		///returns the total length of the line at index j.
+		T
+		getl(int j)
+		{
+			return (L[j]);
+		}
 };