diff --git a/stim/cuda/branch_detection.cuh b/stim/cuda/branch_detection.cuh index abd9b44..11d19e9 100644 --- a/stim/cuda/branch_detection.cuh +++ b/stim/cuda/branch_detection.cuh @@ -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 > 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( gpuI, sigma, x, y, t.getTexture(), t.getArray() ); - stim::gpu2image(gpuI, "Blur.jpg", x,y , 0, 255); +// stim::gpu2image(gpuI, "Blur.jpg", x,y , 0, 255); // stim::gpu2image(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(gpuCenters, gpuVote, final_t, conn, x, y); - stim::gpu2image(gpuCenters, "Centers.jpg", x, y, 0, 1); +// stim::gpu2image(gpuCenters, "Centers.jpg", x, y, 0, 1); cudaMemcpy(cpuCenters, gpuCenters, bytes_ds, cudaMemcpyDeviceToHost); - stim::cpu2image(cpuCenters, "CentersXPU.jpg", x, y, 0, 1); - std::cout << pixels << " " << x << " " << y << std::endl; - for(int i = 0; i < pixels; i++) +// stim::cpu2image(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((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((x_v/(float)x*360), (y_v/(float)y), 0.0)); } - } + } */ + +// std::cout << "After " << output.size() << std::endl; t.UnmapCudaTexture(); diff --git a/stim/gl/gl_spider.h b/stim/gl/gl_spider.h index 09f5975..8fe05bd 100644 --- a/stim/gl/gl_spider.h +++ b/stim/gl/gl_spider.h @@ -229,24 +229,49 @@ class gl_spider : public virtual gl_texture DrawLongCylinder(n, l_template, l_square); stim::cylinder cyl(cL, cM); std::vector< stim::vec > result = find_branch(btexbufferID, GL_TEXTURE_2D, n*l_square, (cL.size()-1)*l_template); + stim::vec 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 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 v = cyl.surf(pval, result[i][0]); // std::cout << v[0] << " ," << v[1] << " ," << v[2] << std::endl; - stim::vec 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(-di[0]+v[0], -di[1]+v[1], -di[2]+v[2]).norm()); - setSeedMag(cyl.r(result[i][1])); + stim::vec 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 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 { 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 || 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 //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 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(); diff --git a/stim/visualization/cylinder.h b/stim/visualization/cylinder.h index c2f7e11..ef263bc 100644 --- a/stim/visualization/cylinder.h +++ b/stim/visualization/cylinder.h @@ -12,8 +12,7 @@ class cylinder { private: stim::circle s; //an arbitrary circle - std::vector< stim::vec > pos; //positions of the cylinder. - std::vector< stim::vec > mags; //radii at each position + std::vector > e; std::vector< T > L; //length of the cylinder at each position. ///default init @@ -27,18 +26,44 @@ class cylinder void init(std::vector > inP, std::vector > inM) { - pos = inP; - mags = inM; +// pos = inP; +// mags = inM; + stim::vec v1; + stim::vec 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 dr = (inP[1] - inP[0]).norm(); + s = stim::circle(inP[0], inM[0][0], dr, stim::vec(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 v1 = (pos[idx]-pos[idx-1]).norm(); - stim::vec v2 = (pos[idx+1]-pos[idx]).norm(); +// return (e[idx+1].P - e[idx].P).norm(); + stim::vec v1 = (e[idx].P-e[idx-1].P).norm(); + stim::vec 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 + 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 ps = p(l, idx); T m = r(l, idx); - stim::vec dr = d(idx); - s = stim::circle(ps, m, dr, stim::vec(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 > > getPoints(int sides) { - if(pos.size() < 2) + std::vector > > points; + points.resize(e.size()); + for(int i = 0; i < e.size(); i++) { - return; - } else { - std::vector > > points; - points.resize(pos.size()); - stim::vec dr = (pos[1] - pos[0]).norm(); - s = stim::circle(pos[0], mags[0][0], dr, stim::vec(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(pos[i], mags[i][0], dr, stim::vec(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]); + } }; -- libgit2 0.21.4