From c01958adbad61ab26da7cafe2bd52749322819fc Mon Sep 17 00:00:00 2001 From: pgovyadi Date: Mon, 20 Feb 2017 18:20:07 -0600 Subject: [PATCH] merged the two stim versions to work for volume-spider, and ONLY volume spider. circle.h, cylinder.h, centerline.h and some other have been changed. Other compilation errors on linux also fixxed in spharmonic.h and gl_spharmonic.h --- stim/biomodels/centerline.h | 389 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------ stim/cuda/filter.cuh | 1 + stim/cuda/ivote/local_max.cuh | 25 ++++++++++++------------- stim/gl/gl_spider.h | 188 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------- stim/image/bmp.h | 6 +++--- stim/image/image.h | 2 +- stim/math/circle.h | 142 ++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------- stim/math/random.h | 2 +- stim/math/spharmonics.h | 52 +++++++++++++++++++++++++++++++++++++++++++++++++--- stim/visualization/cylinder.h | 327 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- stim/visualization/gl_spharmonics.h | 13 ++++++++++--- 11 files changed, 551 insertions(+), 596 deletions(-) diff --git a/stim/biomodels/centerline.h b/stim/biomodels/centerline.h index 46d9369..40d56eb 100644 --- a/stim/biomodels/centerline.h +++ b/stim/biomodels/centerline.h @@ -3,6 +3,7 @@ #include #include +//#include namespace stim{ @@ -11,134 +12,195 @@ namespace stim{ * class to describe an interconnected (often biological) network. */ template -class centerline : public std::vector< stim::vec3 >{ +class centerline{ protected: + unsigned int N; //number of points in the fiber + double **c; //centerline (array of double pointers) +// ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching + + /// Initialize an empty fiber + void init() + { + N=0; + c=NULL; +// kdt = NULL; + } - std::vector 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 d(size_t i) { - if (size() <= 1) return vec3(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 a = at(i) - at(i - 1); - vec3 b = at(i + 1) - at(i); - vec3 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++; - } + /// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0) + void init(unsigned int n) + { + + N = n; //set the number of points +// kdt = NULL; + c = (double**) malloc(sizeof(double*) * N); //allocate the array pointer + + for(unsigned int i = 0; i < N; i++) //allocate space for each point + c[i] = (double*) malloc(sizeof(double) * 3); + } + + /// Copies an existing fiber to the current fiber + + /// @param cpy stores the new copy of the fiber + void copy( const stim::centerline& cpy, bool kd = 0){ + + ///allocate space for the new fiber + init(cpy.N); - stim::vec3 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 + ///copy the points + for(unsigned int i = 0; i < N; i++){ + for(unsigned int d = 0; d < 3; d++) //for each dimension + c[i][d] = cpy.c[i][d]; //copy the coordinate } +// if(kd) +// gen_kdtree(); //generate the kd tree for the new fiber } - void init() { - if (size() == 0) return; //return if there aren't any points - update_L(); + /// 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; // initialize variables + float v; + for(unsigned int d = 0; d < 3; d++) + { + v = p1[d] - p0[d]; + sum +=v * v; + } + return sqrt(sum); } + /// This function retreives the index for the fiber point closest to q + + /// @param q is a reference point used to find the closest point on the fiber center line +// unsigned int ann( stim::vec q ){ +// ANNidxArray idx = new ANNidx[1]; //variable used to hold the nearest point +// ANNdistArray sq_dist = new ANNdist[1]; //variable used to hold the squared distance to the nearest point +// kdt->annkSearch(q.data(), 1, idx, sq_dist); //search the KD tree for the nearest neighbor +// return *idx; +// } + /// Returns a stim::vec representing the point at index i /// @param i is an index of the desired centerline point stim::vec get_vec(unsigned i){ - return std::vector< stim::vec3 >::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;*/ - } + stim::vec3 r; + r.resize(3); + r[0] = c[i][0]; + r[1] = c[i][1]; + r[2] = c[i][2]; - ///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 p(T l, int idx) { - T rat = (l - L[idx]) / (L[idx + 1] - L[idx]); - stim::vec3 v1 = at(idx); - stim::vec3 v2 = at(idx + 1); - return(v1 + (v2 - v1)*rat); + return r; } public: - using std::vector< stim::vec3 >::at; - using std::vector< stim::vec3 >::size; - - centerline() : std::vector< stim::vec3 >() { + centerline(){ init(); } - centerline(size_t n) : std::vector< stim::vec3 >(n){ - init(); + + /// Copy constructor + centerline(const stim::centerline &obj){ + copy(obj); } - //overload the push_back function to update the length vector - void push_back(stim::vec3 p) { - std::vector< stim::vec3 >::push_back(p); - update_L(size() - 1); + //temp constructor for graph visualization + centerline(int n) + { + init(n); } - ///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 p(T pvalue) { - if (pvalue <= 0.0) return at(0); //return the first element - if (pvalue >= 1.0) return back(); //return the last element + /// Constructor takes a list of stim::vec points, the radius at each point is set to zero + centerline(std::vector< stim::vec > p, bool kd = 0){ + init(p.size()); //initialize the fiber + + //for each point, set the centerline position and radius + for(unsigned int i = 0; i < N; i++){ - T l = pvalue*L[L.size() - 1]; - int idx = findIdx(l); - return p(l, idx); + //set the centerline position + for(unsigned int d = 0; d < 3; d++) + c[i][d] = (double) p[i][d]; + + //set the radius + } + //generate a kd tree +// if(kd) +// gen_kdtree(); } - ///Update centerline internal parameters (currently the L vector) - void update() { - init(); + /// constructor takes a list of points + centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){ + init(pos.size()); //initialize the fiber + + //for each point, set the centerline position and radius + for(unsigned int i = 0; i < N; i++){ + //set the centerline position + for(unsigned int d = 0; d < 3; d++) + c[i][d] = (double) pos[i][d]; + //set the radius + } + + //generate a kd tree + //if(kd) + // gen_kdtree(); + } + + /// Assignment operation + centerline& operator=(const centerline &rhs){ + if(this == &rhs) return *this; //test for and handle self-assignment + copy(rhs); + return *this; } - ///Return the length of the entire centerline - T length() { - return L.back(); + + + /// Return the point on the fiber closest to q + /// @param q is the query point used to locate the nearest point on the fiber centerline +// stim::vec nearest(stim::vec q){ +// +// stim::vec temp( (double) q[0], (double) q[1], (double) q[2]); +// +// unsigned int idx = ann(temp); //determine the index of the nearest neighbor +// +// return stim::vec((T) c[idx][0], (T) c[idx][1], (T) c[idx][2]); //return the nearest centerline point +// } + + /// Return the point index on the fiber closest to q + /// @param q is the query point used to locate the nearest point on the fiber centerline +// unsigned int nearest_idx(stim::vec q){ +// +// stim::vec temp((double) q[0], (double) q[1], (double) q[2]); +// +// unsigned int idx = ann(temp); //determine the index of the nearest neighbor +// +// return idx; //return the nearest centerline point index +// } + + /// Returns the fiber centerline as an array of stim::vec points + std::vector< stim::vec > get_centerline(){ + + //create an array of stim vectors + std::vector< stim::vec3 > pts(N); + + //cast each point to a stim::vec, keeping only the position information + for(unsigned int i = 0; i < N; i++) + pts[i] = stim::vec3((T) c[i][0], (T) c[i][1], (T) c[i][2]); + + //return the centerline array + return pts; } /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned std::vector< stim::centerline > split(unsigned int idx){ - std::vector< stim::centerline > fl; //create an array to store up to two fibers - size_t N = size(); + std::vector< stim::centerline > fl; //create an array to store up to two fibers //if the index is an end point, only the existing fiber is returned if(idx == 0 || idx == N-1){ @@ -154,84 +216,123 @@ public: fl.resize(2); //set the array size to 2 - fl[0] = stim::centerline(N1); //set the size of each fiber - fl[1] = stim::centerline(N2); + fl[0].init(N1); //set the size of each fiber + fl[1].init(N2); //copy both halves of the fiber - unsigned int i; + unsigned int i, d; //first half - for(i = 0; i < N1; i++) //for each centerline point - fl[0][i] = std::vector< stim::vec3 >::at(i); - fl[0].init(); //initialize the length vector + for(i = 0; i < N1; i++){ //for each centerline point + for(d = 0; d < 3; d++) + fl[0].c[i][d] = c[i][d]; //copy each coordinate + } //second half - for(i = 0; i < N2; i++) - fl[1][i] = std::vector< stim::vec3 >::at(idx+i); - fl[1].init(); //initialize the length vector + for(i = 0; i < N2; i++){ + for(d = 0; d < 3; d++) + fl[1].c[i][d] = c[idx + i][d]; + + } } return fl; //return the array } + /// Calculates the set of fibers resulting from a connection between the current fiber and a fiber f + + /// @param f is the fiber that will be connected to the current fiber +/* std::vector< stim::centerline > connect( stim::centerline &f, double dist){ + + double min_dist; + unsigned int idx0, idx1; + + //go through each point in the query fiber, looking for the indices for the closest points + for(unsigned int i = 0; i < f.n_pts(); i++){ + //Run through all points and find the index with the closest point, then partition the fiber and return two fibers. + + } + + + + } +*/ /// Outputs the fiber as a string std::string str(){ std::stringstream ss; - size_t N = std::vector< stim::vec3 >::size(); - ss << "---------[" << N << "]---------" << std::endl; - for (size_t i = 0; i < N; i++) - ss << std::vector< stim::vec3 >::at(i) << std::endl; - ss << "--------------------" << std::endl; + + //create an iterator for the point list + //typename std::list< point >::iterator i; + for(unsigned int i = 0; i < N; i++){ + ss<<" [ "; + for(unsigned int d = 0; d < 3; d++){ + ss< back(){ - return std::vector< stim::vec3 >::back(); + /// Returns the number of centerline points in the fiber + unsigned int size(){ + return N; } - ////resample a fiber in the network - stim::centerline resample(T spacing) - { - //std::cout<<"fiber::resample()"< v; //v-direction vector of the segment - stim::vec3 p; //- intermediate point to be added - stim::vec3 p1; // p1 - starting point of an segment on the fiber, - stim::vec3 p2; // p2 - ending point, - //double sum=0; //distance summation + /// Bracket operator returns the element at index i - size_t N = size(); + /// @param i is the index of the element to be returned as a stim::vec + stim::vec operator[](unsigned i){ + return get_vec(i); + } - centerline new_c; // initialize list of new resampled points on the fiber + /// Back method returns the last point in the fiber + stim::vec back(){ + return get_vec(N-1); + } + ////resample a fiber in the network + stim::centerline resample(T spacing) + { + std::cout<<"fiber::resample()"< v(3); //v-direction vector of the segment + stim::vec p(3); //- intermediate point to be added + stim::vec p1(3); // p1 - starting point of an segment on the fiber, + stim::vec p2(3); // p2 - ending point, + double sum=0; //distance summation + std::vector > fiberPositions = centerline(); + std::vector > 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++) - { - 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= 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(centers, res, conn, DIM_X, DIM_Y); stim::cuda::gpu_local_max(centers, res, threshold, conn, DIM_X, DIM_Y); cudaDeviceSynchronize(); diff --git a/stim/cuda/ivote/local_max.cuh b/stim/cuda/ivote/local_max.cuh index 38b7096..b342a6f 100644 --- a/stim/cuda/ivote/local_max.cuh +++ b/stim/cuda/ivote/local_max.cuh @@ -10,24 +10,24 @@ namespace stim{ // this kernel calculates the local maximum for finding the cell centers template - __global__ void cuda_local_max(T* gpuCenters, T* gpuVote, int conn, size_t x, size_t y){ + __global__ void cuda_local_max(T* gpuCenters, T* gpuVote, T final_t, int conn, int x, int y){ // calculate the 2D coordinates for this current thread. - size_t xi = blockIdx.x * blockDim.x + threadIdx.x; - size_t yi = blockIdx.y * blockDim.y + threadIdx.y; + int xi = blockIdx.x * blockDim.x + threadIdx.x; + int yi = blockIdx.y * blockDim.y + threadIdx.y; if(xi >= x || yi >= y) return; // convert 2D coordinates to 1D - size_t i = yi * x + xi; + int i = yi * x + xi; gpuCenters[i] = 0; //initialize the value at this location to zero T val = gpuVote[i]; //compare to the threshold - //if(val < final_t) return; + if(val < final_t) return; //define an array to store indices with same vote value /*int * IdxEq; @@ -56,25 +56,24 @@ namespace stim{ return; } } */ - //gpuCenters[i] = 1; - gpuCenters[i] = gpuVote[i]; + gpuCenters[i] = 1; } template - void gpu_local_max(T* gpuCenters, T* gpuVote, unsigned int conn, size_t x, size_t y){ + void gpu_local_max(T* gpuCenters, T* gpuVote, T final_t, unsigned int conn, unsigned int x, unsigned int y){ unsigned int max_threads = stim::maxThreadsPerBlock(); /*dim3 threads(max_threads, 1); dim3 blocks(x/threads.x + (x %threads.x == 0 ? 0:1) , y);*/ - dim3 threads((unsigned int)sqrt(max_threads), (unsigned int)sqrt(max_threads) ); - dim3 blocks((unsigned int)x/threads.x + 1, (unsigned int)y/threads.y + 1); + dim3 threads( sqrt(max_threads), sqrt(max_threads) ); + dim3 blocks(x/threads.x + 1, y/threads.y + 1); //call the kernel to find the local maximum. - cuda_local_max <<< blocks, threads >>>(gpuCenters, gpuVote, conn, x, y); + cuda_local_max <<< blocks, threads >>>(gpuCenters, gpuVote, final_t, conn, x, y); } template - void cpu_local_max(T* cpuCenters, T* cpuVote, unsigned int conn, unsigned int x, unsigned int y){ + void cpu_local_max(T* cpuCenters, T* cpuVote, T final_t, unsigned int conn, unsigned int x, unsigned int y){ //calculate the number of bytes in the array unsigned int bytes = x * y * sizeof(T); @@ -91,7 +90,7 @@ namespace stim{ HANDLE_ERROR(cudaMemcpy(gpuVote, cpuVote, bytes, cudaMemcpyHostToDevice)); //call the GPU version of the local max function - gpu_local_max(gpuCenters, gpuVote, conn, x, y); + gpu_local_max(gpuCenters, gpuVote, final_t, conn, x, y); //copy the cell centers data to the CPU cudaMemcpy(cpuCenters, gpuCenters, bytes, cudaMemcpyDeviceToHost) ; diff --git a/stim/gl/gl_spider.h b/stim/gl/gl_spider.h index e2ae0ea..7dfd7e7 100644 --- a/stim/gl/gl_spider.h +++ b/stim/gl/gl_spider.h @@ -2,33 +2,47 @@ #define STIM_GL_SPIDER_H //#include +///basic GL and Cuda libs #include #include #include #include -#include + +///stim .*h for gl tracking and visualization #include -#include #include +#include +#include +#include +#include +#include + +///stim *.h for CUDA +#include +#include + +///stim *.h for branch detection +#include +#include +#include "../../../volume-spider/glnetwork.h" + +/// *.h for math +#include #include #include #include -#include #include -#include -#include #include -#include -#include -#include -#include +#include + + +///c++ libs for everything #include #include -#include -#include "../../../volume-spider/glnetwork.h" -#include #include #include + +///Timing and debugging #ifdef TIMING #include #include @@ -39,6 +53,7 @@ #include #endif +// #include #ifdef DEBUG #include #endif @@ -61,6 +76,7 @@ class gl_spider // : public virtual gl_texture double cost_time;// = 0; double network_time;// = 0; double hit_time;// = 0; + int iter_t; #endif stim::vec3 p; //vector designating the position of the spider. @@ -72,7 +88,7 @@ class gl_spider // : public virtual gl_texture std::vector mV; //A list of all the size vectors. std::vector lV; //A list of all the size vectors. - stim::matrix cT; //current Transformation matrix (tissue)->(texture) + stim::matrix_sq cT; //current Transformation matrix (tissue)->(texture) GLuint texID; //OpenGL ID for the texture to be traced stim::vec3 S; //Size of a voxel in the volume. stim::vec3 R; //Dimensions of the volume. @@ -339,31 +355,6 @@ class gl_spider // : public virtual gl_texture #endif } - - float uniformRandom() - { - return ( (float)(rand()))/( (float)(RAND_MAX)); ///generates a random number between 0 and 1 using the uniform distribution. - } - - float normalRandom() - { - float u1 = uniformRandom(); - float u2 = uniformRandom(); - return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1)); ///generate a random number using the normal distribution between 0 and 1. - } - - stim::vec3 uniformRandVector() - { - stim::vec3 r(uniformRandom(), uniformRandom(), 1.0); ///generate a random vector using the uniform distribution between 0 and 1. - return r; - } - - stim::vec3 normalRandVector() - { - stim::vec3 r(normalRandom(), normalRandom(), 1.0); ///generate a random vector using the normal distribution between 0 and 1. - return r; - } - //--------------------------------------------------------------------------// //---------------------TEMPLATE CREATION METHODS----------------------------// @@ -400,8 +391,8 @@ class gl_spider // : public virtual gl_texture glNewList(dList, GL_COMPILE); ///create a display list of all the direction templates. for(int i = 0; i < numSamples; i++) ///for each sample { - z = uniformRandom()*range + Z[1]; ///generate a z coordinate - theta = uniformRandom()*stim::TAU; ///generate a theta coordinate + z = stim::Random::uniformRandom()*range + Z[1]; ///generate a z coordinate + theta = stim::Random::uniformRandom()*stim::TAU; ///generate a theta coordinate phi = acos(z); ///generate a phi from the z. stim::vec3 sph(1, theta, phi); ///combine into a vector in spherical coordinates. stim::vec3 cart = sph.sph2cart();///convert to cartesian. @@ -460,7 +451,7 @@ class gl_spider // : public virtual gl_texture UpdateBuffer(0.0, 0.0+0*n_pixels); for(int i = 1; i < numSamplesPos; i++) ///for the number of position samples { - stim::vec3 temp = uniformRandVector(); ///generate a random point on a plane. + stim::vec3 temp = stim::Random::uniformRandVector(); ///generate a random point on a plane. temp[0] = temp[0]*delta; temp[1] = temp[1]*2*stim::PI; @@ -953,6 +944,7 @@ class gl_spider // : public virtual gl_texture cost_time = 0; network_time = 0; hit_time = 0; + iter_t = 0; #endif #ifdef DEBUG iter = 0; @@ -960,7 +952,7 @@ class gl_spider // : public virtual gl_texture iter_dir = 0; iter_siz = 0; #endif - stepsize = 6.0; + stepsize = 3.0; n_pixels = 16.0; srand(100); @@ -1274,12 +1266,19 @@ class gl_spider // : public virtual gl_texture { std::vector ret; ret.resize(7); + // ret[0] = branch_time/(float)iter; ret[0] = branch_time; + // ret[1] = direction_time/(float)iter; ret[1] = direction_time; + // ret[2] = position_time/(float)iter; ret[2] = position_time; + // ret[3] = size_time/(float)iter; ret[3] = size_time; + // ret[4] = cost_time/(float)iter; ret[4] = cost_time; + // ret[5] = network_time/(float)iter; ret[5] = network_time; + // ret[6] = hit_time/(float)iter; ret[6] = hit_time; return ret; @@ -1311,6 +1310,7 @@ class gl_spider // : public virtual gl_texture setSeed(x, y, z); setSeedVec(u, v, w); setSeedMag(m); + std::cout << getLastSeed() << getLastSeedVec() << std::cout; } myfile.close(); } else { @@ -1320,7 +1320,7 @@ class gl_spider // : public virtual gl_texture ///Saves the network to a file. void - saveNetwork(std::string name) + saveNetwork(std::string name, int xoffset = 0, int yoffset = 0, int zoffset = 0) { stim::glObj sk1; for(int i = 0; i < nt.sizeE(); i++) @@ -1331,7 +1331,7 @@ class gl_spider // : public virtual gl_texture for(int j = 0; j < ce.size(); j++) { sk1.TexCoord(cm[j]); - sk1.Vertex(ce[j][0], ce[j][1], ce[j][2]); + sk1.Vertex(ce[j][0]+xoffset, ce[j][1]+yoffset, ce[j][2]+zoffset); } sk1.End(); } @@ -1408,6 +1408,11 @@ class gl_spider // : public virtual gl_texture #ifdef DEBUG std::cerr << std::endl; #endif + + #ifdef TIMING + iter_t++; + #endif + return current_cost; } @@ -1422,36 +1427,6 @@ class gl_spider // : public virtual gl_texture //--------------------------------------------------------------------------// //-----------------------------EXPERIMENTAL METHODS-------------------------// //--------------------------------------------------------------------------// - - void - MonteCarloDirectionVectors(int nSamples, float solidAngle = stim::TAU) - { -// float PHI[2];//, Z[2];//, range; -// PHI[0] = asin(solidAngle/2); -// PHI[1] = asin(0); - -// Z[0] = cos(PHI[0]); -// Z[1] = cos(PHI[1]); - -// range = Z[0] - Z[1]; - - std::vector > vecsUni; - for(int i = 0; i < numSamplesPos; i++) - { - stim::vec3 a(uniformRandom()*0.8, uniformRandom()*0.8, 0.0); - a[0] = a[0]-0.4; - a[1] = a[1]-0.4; - vecsUni.push_back(a); - } - - stringstream name; - for(int i = 0; i < numSamplesPos; i++) - name << vecsUni[i].str() << std::endl; - - std::ofstream outFile; - outFile.open("New_Pos_Vectors.txt"); - outFile << name.str().c_str(); - } /* void DrawCylinder() @@ -1570,9 +1545,6 @@ class gl_spider // : public virtual gl_texture { //Define the varibles and turn on Selection Mode - #ifdef TIMING - gpuStartTimer(); - #endif GLuint selectBuf[2048]; ///size of the selection buffer in bytes. GLint hits; ///hit fibers @@ -1602,6 +1574,9 @@ class gl_spider // : public virtual gl_texture glPushMatrix(); glLoadIdentity(); + #ifdef TIMING + gpuStartTimer(); + #endif CHECK_OPENGL_ERROR gluLookAt(ps[0], ps[1], ps[2], ds[0], ds[1], ds[2], @@ -1849,6 +1824,65 @@ class gl_spider // : public virtual gl_texture std::cout << "I broke out!" << std::endl; #endif } +//--------------------------------------------------------------------------// +//--------------------------------DEPRECIATED-------------------------------// +//--------------------------------------------------------------------------// +/* + float uniformRandom() + { + return ( (float)(rand()))/( (float)(RAND_MAX)); ///generates a random number between 0 and 1 using the uniform distribution. + } + + float normalRandom() + { + float u1 = uniformRandom(); + float u2 = uniformRandom(); + return cos(2.0*atan(1.0)*u2)*sqrt(-1.0*log(u1)); ///generate a random number using the normal distribution between 0 and 1. + } + + stim::vec3 uniformRandVector() + { + stim::vec3 r(uniformRandom(), uniformRandom(), 1.0); ///generate a random vector using the uniform distribution between 0 and 1. + return r; + } + + stim::vec3 normalRandVector() + { + stim::vec3 r(normalRandom(), normalRandom(), 1.0); ///generate a random vector using the normal distribution between 0 and 1. + return r; + } + + void + MonteCarloDirectionVectors(int nSamples, float solidAngle = stim::TAU) + { +// float PHI[2];//, Z[2];//, range; +// PHI[0] = asin(solidAngle/2); +// PHI[1] = asin(0); + +// Z[0] = cos(PHI[0]); +// Z[1] = cos(PHI[1]); + +// range = Z[0] - Z[1]; + + std::vector > vecsUni; + for(int i = 0; i < numSamplesPos; i++) + { + stim::vec3 a(uniformRandom()*0.8, uniformRandom()*0.8, 0.0); + a[0] = a[0]-0.4; + a[1] = a[1]-0.4; + vecsUni.push_back(a); + } + + stringstream name; + for(int i = 0; i < numSamplesPos; i++) + name << vecsUni[i].str() << std::endl; + + std::ofstream outFile; + outFile.open("New_Pos_Vectors.txt"); + outFile << name.str().c_str(); + } + +*/ }; } #endif diff --git a/stim/image/bmp.h b/stim/image/bmp.h index 855d438..defd6dc 100644 --- a/stim/image/bmp.h +++ b/stim/image/bmp.h @@ -187,7 +187,7 @@ namespace stim { } bool open(std::string filename) { //open the bitmap file and read the header data - file.open(filename, std::ifstream::binary); + file.open(filename.c_str(), std::ifstream::binary); if (!file) { std::cout << "stim::bmp ERROR: error opening file: " << filename.c_str() << std::endl; return false; @@ -246,7 +246,7 @@ namespace stim { info_header.bcSize = sizeof(tagBITMAPCOREHEADER); info_header.bcPlanes = 1; - std::ofstream outfile(filename, std::ios::binary); //open the output file for binary writing + std::ofstream outfile(filename.c_str(), std::ios::binary); //open the output file for binary writing outfile.write((char*)&file_header, sizeof(tagBITMAPFILEHEADER)); //write the file header outfile.write((char*)&info_header, sizeof(tagBITMAPCOREHEADER)); //write the information header @@ -259,4 +259,4 @@ namespace stim { free(pad); return true; } -} \ No newline at end of file +} diff --git a/stim/image/image.h b/stim/image/image.h index a479eb5..561c18b 100644 --- a/stim/image/image.h +++ b/stim/image/image.h @@ -606,7 +606,7 @@ public: //crop regions given by an array of 1D index values std::vector< image > crop_idx(size_t w, size_t h, std::vector idx) { - std::vector> result(idx.size()); //create an array of image files to return + std::vector > result(idx.size()); //create an array of image files to return for (size_t i = 0; i < idx.size(); i++) { //for each specified index point size_t y = idx[i] / X(); //calculate the y coordinate from the 1D index (center of ROI) size_t x = idx[i] - y * X(); //calculate the x coordinate (center of ROI) diff --git a/stim/math/circle.h b/stim/math/circle.h index 23b5c1a..423a9d3 100644 --- a/stim/math/circle.h +++ b/stim/math/circle.h @@ -18,14 +18,13 @@ class circle : plane private: - //stim::vec3 Y; - T R; //radius of the circle + stim::vec3 Y; - /*CUDA_CALLABLE void + CUDA_CALLABLE void init() { Y = U.cross(N).norm(); - }*/ + } public: using stim::plane::n; @@ -35,8 +34,6 @@ public: using stim::plane::rotate; using stim::plane::setU; - using stim::plane::init; - ///base constructor ///@param th value of the angle of the starting point from 0 to 360. CUDA_CALLABLE @@ -45,28 +42,26 @@ public: init(); } - ///create a circle given a size and position in Z space. + ///create a rectangle given a size and position in Z space. ///@param size: size of the rectangle in ND space. ///@param z_pos z coordinate of the rectangle. CUDA_CALLABLE - circle(T radius, T z_pos = (T)0) : plane(z_pos) + circle(T size, T z_pos = (T)0) : plane() { - //center(stim::vec3(0, 0, z_pos)); - //scale(size); - //init(); - R = radius; + center(stim::vec3(0,0,z_pos)); + scale(size); + init(); } ///create a rectangle from a center point, normal ///@param c: x,y,z location of the center. ///@param n: x,y,z direction of the normal. CUDA_CALLABLE - circle(vec3 c, vec3 n = vec3(0,0,1)) : plane(n, c) + circle(vec3 c, vec3 n = vec3(0,0,1)) : plane() { - //center(c); - //normal(n); - //init(); - R = (T)1; + center(c); + normal(n); + init(); } /*///create a rectangle from a center point, normal, and size @@ -89,18 +84,14 @@ public: ///@param n: x,y,z direction of the normal. ///@param u: x,y,z direction for the zero vector (from where the rotation starts) CUDA_CALLABLE - circle(vec3 c, T r, vec3 n = vec3(0,0,1), vec3 u = vec3(1, 0, 0)) : plane() + circle(vec3 c, T s, vec3 n = vec3(0,0,1), vec3 u = vec3(1, 0, 0)) : plane() { - P = c; - N = n; + init(); setU(u); - R = r; - //init(); - //setU(u); // U = u; - //center(c); - //normal(n); - //scale(s); + center(c); + normal(n); + scale(s); } ///scales the circle by a certain factor @@ -108,111 +99,86 @@ public: CUDA_CALLABLE void scale(T factor) { - //U *= factor; - //Y *= factor; - R *= factor; - } - - ///set the radius of circle to a certain value - ///@param value: the new radius of the circle - CUDA_CALLABLE - void set_R(T value) - { - R = value; + U *= factor; + Y *= factor; } ///sets the normal for the cirlce ///@param n: x,y,z direction of the normal. CUDA_CALLABLE void - normal(vec3 n){ - rotate(n); + normal(vec3 n) + { + rotate(n, Y); } ///sets the center of the circle. ///@param n: x,y,z location of the center. CUDA_CALLABLE void center(vec3 p){ - P = p; + this->P = p; } ///boolean comparison bool operator==(const circle & rhs) { - if(P == rhs.P && U == rhs.U) + if(P == rhs.P && U == rhs.U && Y == rhs.Y) return true; else return false; } - //returns the point in world space corresponding to the polar coordinate (r, theta) - CUDA_CALLABLE stim::vec3 - p(T r, T theta) { - T u = r * cos(theta); //calculate the coordinates in the planar space defined by the circle - T v = r * sin(theta); - - vec3 V = U.cross(N); //calculate the orthogonal vector V - return P + U * u + V * v; //calculate the cartesian coordinate of the specified point - } - - //returns the point in world space corresponding to the value theta at radius R - CUDA_CALLABLE stim::vec3 - p(T theta) { - return p(R, theta); - } - - //get the world space value given the polar coordinates (r, theta) - ///get the world space value given the planar coordinates a, b in [0, 1] - /*CUDA_CALLABLE stim::vec3 p(T a, T b) + CUDA_CALLABLE stim::vec3 p(T a, T b) { stim::vec3 result; vec3 A = this->P - this->U * (T)0.5 - Y * (T)0.5; result = A + this->U * a + Y * b; return result; - }*/ + } ///parenthesis operator returns the world space given rectangular coordinates a and b in [0 1] - CUDA_CALLABLE stim::vec3 operator()(T r, T theta) + CUDA_CALLABLE stim::vec3 operator()(T a, T b) { - return p(r, theta); - } - - //parenthesis operator returns the world space coordinate at the edge of the circle given theta - CUDA_CALLABLE stim::vec3 operator()(T theta) { - return p(theta); + return p(a,b); } ///returns a vector with the points on the initialized circle. ///connecting the points results in a circle. ///@param n: integer for the number of points representing the circle. - std::vector< stim::vec3 > points(unsigned n) { - std::vector< stim::vec3 > result(n); //initialize a vector of n points - - float dt = stim::TAU / n; - for (unsigned i = 0; i < n; i++) - result[i] = p(i * dt); //calculate a point on the edge of the circle + std::vector > + getPoints(int n) + { + std::vector > result; + stim::vec3 point; + T x,y; + float step = 360.0/(float) n; + for(float j = 0; j <= 360.0; j += step) + { + y = 0.5*cos(j*stim::TAU/360.0)+0.5; + x = 0.5*sin(j*stim::TAU/360.0)+0.5; + result.push_back(p(x,y)); + } return result; } - - ///returns a vector with the points on the initialized circle - ///connecting the points results in a circle - ///@param n: integer for the number of points representing the circle - ///the only difference between points and glpoints is that the first point appears twice in the returning lists - std::vector< stim::vec3 > glpoints(unsigned n) { - std::vector< stim::vec3 > result(n + 1); - float dt = stim::TAU / n; - for (unsigned i = 0; i < n; i++) - result[i] = p(i * dt); - result[n] = p(0); //close the circle! - return result; - } + ///returns a vector with the points on the initialized circle. + ///connecting the points results in a circle. + ///@param n: integer for the number of points representing the circle. + CUDA_CALLABLE stim::vec3 + p(T theta) + { + T x,y; + y = 0.5*cos(theta*STIM_TAU/360.0)+0.5; + x = 0.5*sin(theta*STIM_TAU/360.0)+0.5; + return p(x,y); + } + std::string str() const { std::stringstream ss; - ss << "r = "< > sph_pts, unsigned int l, int m, stim::vec3 c = stim::vec3(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector w = std::vector()) + void pdf(std::vector > sph_pts, unsigned int l, int m, stim::vec3 c = stim::vec3(0, 0, 0), unsigned int n = 1000, bool norm = true, std::vector w = std::vector()) { std::vector weights; ///the weight at each point on the surface of the sphere. // weights.resize(n); @@ -223,7 +223,7 @@ namespace stim { } } mcEnd(); - }*/ + } std::string str() { @@ -313,6 +313,52 @@ namespace stim { T operator()(T theta, T phi) { return p(theta, phi); } + + //overload arithmetic operations + + spharmonics operator*(T rhs) const { + + spharmonics result(C.size()); //create a new spherical harmonics object + + for (size_t c = 0; c < C.size(); c++) //for each coefficient + + result.C[c] = C[c] * rhs; //calculate the factor and store the result in the new spharmonics object + + return result; + + } + + + + spharmonics operator+(spharmonics rhs) { + + size_t low = std::min(C.size(), rhs.C.size()); //store the number of coefficients in the lowest object + size_t high = std::max(C.size(), rhs.C.size()); //store the number of coefficients in the result + bool rhs_lowest = false; //true if rhs has the lowest number of coefficients + if (rhs.C.size() < C.size()) rhs_lowest = true; //if rhs has a lower number of coefficients, set the flag + + + + spharmonics result(high); //create a new object + + size_t c; + for (c = 0; c < low; c++) //perform the first batch of additions + result.C[c] = C[c] + rhs.C[c]; //perform the addition + + for (c = low; c < high; c++) { + if (rhs_lowest) + result.C[c] = C[c]; + else + result.C[c] = rhs.C[c]; + } + return result; + } + + + + spharmonics operator-(spharmonics rhs) { + return (*this) + (rhs * (T)(-1)); + } /// Fill an NxN grid with the spherical function for theta = [0 2pi] and phi = [0 pi] void get_func(T* data, size_t X, size_t Y) { T dt = stim::TAU / (T)X; //calculate the step size in each direction @@ -351,4 +397,4 @@ namespace stim { } -#endif \ No newline at end of file +#endif diff --git a/stim/visualization/cylinder.h b/stim/visualization/cylinder.h index 9c218cf..ff3fb82 100644 --- a/stim/visualization/cylinder.h +++ b/stim/visualization/cylinder.h @@ -3,268 +3,56 @@ #include #include #include -#include +/* + +*/ namespace stim { template -class cylinder : public centerline { -protected: - - using stim::centerline::d; +class cylinder + : public centerline +{ + private: + stim::circle s; //an arbitrary circle + std::vector > e; //an array of circles that store the centerline - std::vector< stim::vec3 > U; //stores the array of U vectors defining the Frenet frame - std::vector< T > R; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius) + std::vector > norms; + std::vector > Us; + std::vector > mags; //stores a list of magnitudes for each point in the centerline (assuming mags[0] is the radius) + std::vector< T > L; //length of the cylinder at each position (pre-integration) - using stim::centerline::findIdx; - - //calculates the U values for each point to initialize the frenet frame - // this function assumes that the centerline has already been set - void init() { - U.resize(size()); //allocate space for the frenet frame vectors - R.resize(size()); - - stim::circle c; //create a circle - stim::vec3 d0, d1; - for (size_t i = 0; i < size() - 1; i++) { //for each line segment in the centerline - c.rotate(d(i)); //rotate the circle to match that normal - U[i] = c.U; //save the U vector from the circle - } - U[size() - 1] = c.U; //for the last point, duplicate the final frenet frame vector - } - -public: - - using stim::centerline::size; - using stim::centerline::at; - using stim::centerline::L; - using stim::centerline::length; - - cylinder() : centerline(){} - - cylinder(centerlinec) : centerline(c) { - init(); - } - - //cylinder(centerlinec, T r) : centerline(c) { - // init(); - // //add_mag(r); - //} - - //initialize a cylinder with a list of points and magnitude values - //cylinder(centerlinec, std::vector r) : centerline(c){ - // init(); - // //add_mag(r); - //} - - //copy the original radius - void copy_r(std::vector radius) { - for (unsigned i = 0; i < radius.size(); i++) - R[i] = radius[i]; - } - - ///Returns magnitude i at the given length into the fiber (based on the pvalue). - ///Interpolates the position 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. - T r(T l, int idx) { - T a = (l - L[idx]) / (L[idx + 1] - L[idx]); - T v2 = (R[idx] + (R[idx + 1] - R[idx])*a); - - return v2; - } - - ///Returns the ith magnitude at the given p-value (p value ranges from 0 to 1). - ///interpolates the radius along the line. - ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1). - T rl(T pvalue) { - if (pvalue <= 0.0) return R[0]; - if (pvalue >= 1.0) return R[size() - 1]; - - T l = pvalue*L[L.size() - 1]; - int idx = findIdx(l); - return r(l, idx); - } - - /// Returns the magnitude at the given index - /// @param i is the index of the desired point - /// @param r is the index of the magnitude value - T r(unsigned i) { - return R[i]; - } - - ///adds a magnitude to each point in the cylinder - /*void add_mag(V val = 0) { - if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline - for (size_t i = 0; i < size(); i++) //for each point - R[i].push_back(val); //add this value to the magnitude vector at each point - }*/ - - //adds a magnitude based on a list of magnitudes for each point - /*void add_mag(std::vector val) { - if (M.size() == 0) M.resize(size()); //if the magnitude vector isn't initialized, resize it to match the centerline - for (size_t i = 0; i < size(); i++) //for each point - R[i].push_back(val[i]); //add this value to the magnitude vector at each point - }*/ - - //sets the value of magnitude m at point i - void set_r(size_t i, T r) { - R[i] = r; - } - - /*size_t nmags() { - if (M.size() == 0) return 0; - else return R[0].size(); - }*/ - - ///Returns a circle representing the cylinder cross section at point i - stim::circle circ(size_t i) { - return stim::circle(at(i), R[i], d(i), U[i]); - } - - ///Returns an OBJ object representing the cylinder with a radial tesselation value of N using magnitude m - stim::obj OBJ(size_t N) { - stim::obj out; //create an OBJ object - stim::circle c0, c1; - std::vector< stim::vec3 > p0, p1; //allocate space for the point sets representing the circles bounding each cylinder segment - T u0, u1, v0, v1; //allocate variables to store running texture coordinates - for (size_t i = 1; i < size(); i++) { //for each line segment in the cylinder - c0 = circ(i - 1); //get the two circles bounding the line segment - c1 = circ(i); - - p0 = c0.points(N); //get t points for each of the end caps - p1 = c1.points(N); - - u0 = L[i - 1] / length(); //calculate the texture coordinate (u, v) where u runs along the cylinder length - u1 = L[i] / length(); - - for (size_t n = 1; n < N; n++) { //for each point in the circle - v0 = (double)(n-1) / (double)(N - 1); //v texture coordinate runs around the cylinder - v1 = (double)(n) / (double)(N - 1); - out.Begin(OBJ_FACE); //start a face (quad) - out.TexCoord(u0, v0); - out.Vertex(p0[n - 1][0], p0[n - 1][1], p0[n - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment - out.TexCoord(u1, v0); - out.Vertex(p1[n - 1][0], p1[n - 1][1], p1[n - 1][2]); - - out.TexCoord(u0, v1); - out.Vertex(p1[n][0], p1[n][1], p1[n][2]); - out.TexCoord(u1, v1); - out.Vertex(p0[n][0], p0[n][1], p0[n][2]); - out.End(); - } - v0 = (double)(N - 2) / (double)(N - 1); //v texture coordinate runs around the cylinder - v1 = 1.0; - out.Begin(OBJ_FACE); - out.TexCoord(u0, v0); - out.Vertex(p0[N - 1][0], p0[N - 1][1], p0[N - 1][2]); //output the points composing a strip of quads wrapping the cylinder segment - out.TexCoord(u1, v0); - out.Vertex(p1[N - 1][0], p1[N - 1][1], p1[N - 1][2]); - - out.TexCoord(u0, v1); - out.Vertex(p1[0][0], p1[0][1], p1[0][2]); - out.TexCoord(u1, v1); - out.Vertex(p0[0][0], p0[0][1], p0[0][2]); - out.End(); - } - return out; - } - - std::string str() { - std::stringstream ss; - size_t N = std::vector< stim::vec3 >::size(); - ss << "---------[" << N << "]---------" << std::endl; - for (size_t i = 0; i < N; i++) - ss << std::vector< stim::vec3 >::at(i) << " r = " << R[i] << " u = " << U[i] << std::endl; - ss << "--------------------" << std::endl; - - return ss.str(); - } - - /// Integrates a magnitude value along the cylinder. - /// @param m is the magnitude value to be integrated (this is usually the radius) - T integrate() { - T sum = 0; //initialize the integral to zero - if (L.size() == 1) - return sum; - else { - T m0, m1; //allocate space for both magnitudes in a single segment - m0 = R[0]; //initialize the first point and magnitude to the first point in the cylinder - T len = L[1]; //allocate space for the segment length + using stim::centerline::c; + using stim::centerline::N; + + ///default init + void + init() + { - for (unsigned p = 1; p < size(); p++) { //for every consecutive point in the cylinder - m1 = R[p]; - if (p > 1) len = (L[p] - L[p - 1]); //calculate the segment length using the L array - sum += (m0 + m1) / (T)2.0 * len; //add the average magnitude, weighted by the segment length - m0 = m1; //move to the next segment by shifting points - } - return sum; //return the integral - } - } - - /// Resamples the cylinder to provide a maximum distance of "spacing" between centerline points. All current - /// centerline points are guaranteed to exist in the new cylinder - /// @param spacing is the maximum spacing allowed between sample points - cylinder resample(T spacing) { - cylinder c = stim::centerline::resample(spacing); //resample the centerline and use it to create a new cylinder - - //size_t nm = nmags(); //get the number of magnitude values in the current cylinder - //if (nm > 0) { //if there are magnitude values - // std::vector magvec(nm, 0); //create a magnitude vector for a single point - // c.M.resize(c.size(), magvec); //allocate space for a magnitude vector at each point of the new cylinder - //} - - T l, t; - for (size_t i = 0; i < c.size(); i++) { //for each point in the new cylinder - l = c.L[i]; //get the length along the new cylinder - t = l / length(); //calculate the parameter value along the new cylinder - //for (size_t mag = 0; mag < nm; mag++) { //for each magnitude value - c.R[i] = r(t); //retrieve the interpolated magnitude from the current cylinder and store it in the new one - //} - } - return c; - } - - std::vector< stim::cylinder > split(unsigned int idx) { - - unsigned N = size(); - std::vector< stim::centerline > LL; - LL.resize(2); - LL = (*this).centerline::split(idx); - std::vector< stim::cylinder > C(LL.size()); - unsigned i = 0; - C[0] = LL[0]; - //C[0].R.resize(idx); - for (; i < idx + 1; i++) { - //for(unsigned d = 0; d < 3; d++) - //C[0][i][d] = LL[0].c[i][d]; - C[0].R[i] = R[i]; - //C[0].R[i].resize(1); - } - if (C.size() == 2) { - C[1] = LL[1]; - i--; - //C[1].M.resize(N - idx); - for (; i < N; i++) { - //for(unsigned d = 0; d < 3; d++) - //C[1][i - idx][d] = LL[1].c[i - idx][d]; - C[1].R[i - idx] = R[i]; - //C[1].M[i - idx].resize(1); - } } - return C; - } - + ///Calculate the physical length of the fiber at each point in the fiber. + void + calculateL() + { +/* L.resize(inP.size()); + T temp = (T)0; + L[0] = 0; + for(size_t i = 1; i < L.size(); i++) + { + temp += (inP[i-1] - inP[i]).len(); + L[i] = temp; + } +*/ } - /* - ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and magnitudes (inM) + ///inits the cylinder from a list of points (std::vector of stim::vec3 --inP) and radii (inM) void - init(centerline inP, std::vector< std::vector > inM) { - M = inM; //the magnitude vector can be copied directly - (*this) = inP; //the centerline can be copied to this class directly + init(std::vector > inP, std::vector > inM) + { + mags = inM; stim::vec3 v1; stim::vec3 v2; e.resize(inP.size()); @@ -286,7 +74,7 @@ public: } stim::vec3 dr = (inP[1] - inP[0]).norm(); - s = stim::circle(inP[0], inR[0][0], dr, stim::vec3(1,0,0)); + s = stim::circle(inP[0], inM[0][0], dr, stim::vec3(1,0,0)); norms[0] = s.N; e[0] = s; Us[0] = e[0].U; @@ -300,7 +88,7 @@ public: norms[i] = s.N; - s.scale(inR[i][0]/inR[i-1][0]); + s.scale(inM[i][0]/inM[i-1][0]); e[i] = s; Us[i] = e[i].U; } @@ -312,7 +100,7 @@ public: norms[j] = s.N; - s.scale(inR[j][0]/inR[j-1][0]); + s.scale(inM[j][0]/inM[j-1][0]); e[j] = s; Us[j] = e[j].U; } @@ -439,7 +227,7 @@ public: stim::vec zero(0.0,0.0); temp.resize(inM.size(), zero); for(int i = 0; i < inM.size(); i++) - temp[i][0] = inR[i]; + temp[i][0] = inM[i]; init(inP, temp); } @@ -458,18 +246,13 @@ public: init(inP, inM); } - //assignment operator creates a cylinder from a centerline (default radius is zero) - cylinder& operator=(stim::centerline c) { - init(c); - } - ///Returns the number of points on the cylinder centerline unsigned int size(){ return N; } - + ///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). @@ -483,7 +266,22 @@ public: T l = pvalue*L[L.size()-1]; int idx = findIdx(l); return (p(l,idx)); - } +/* T rat = (l-L[idx])/(L[idx+1]-L[idx]); + stim::vec3 v1( + c[idx][0], + c[idx][1], + c[idx][2] + ); + + stim::vec3 v2( + c[idx+1][0], + c[idx+1][1], + c[idx+1][2] + ); + +// return( e[idx].P + (e[idx+1].P-e[idx].P)*rat); + return( v1 + (v2 - v1)*rat); +*/ } ///Returns a position vector at the given length into the fiber (based on the pvalue). ///Interpolates the radius along the line. @@ -523,7 +321,10 @@ public: T l = pvalue*L[L.size()-1]; int idx = findIdx(l); return (r(l,idx)); - } +/* T rat = (l-L[idx])/(L[idx+1]-L[idx]); +// return (mags[idx][0] + (mags[idx+1][0]-mags[idx][0])*rat); + return (e[idx].U.len() + (e[idx+1].U.len() - e[idx].U.len())*rat); +*/ } ///Returns a radius at the given length into the fiber (based on the pvalue). ///Interpolates the position along the line. @@ -705,10 +506,10 @@ public: return cylinder(result); - }*/ + } }; } -#endif \ No newline at end of file +#endif diff --git a/stim/visualization/gl_spharmonics.h b/stim/visualization/gl_spharmonics.h index 9882325..098f51f 100644 --- a/stim/visualization/gl_spharmonics.h +++ b/stim/visualization/gl_spharmonics.h @@ -51,6 +51,13 @@ namespace stim { magnitude = true; } + gl_spharmonics(stim::spharmonics disp, stim::spharmonics color, size_t slices) : + gl_spharmonics(slices) + { + Sc = color; + Sd = disp; + } + ~gl_spharmonics() { if (dlist) glDeleteLists(dlist, 1); //delete the display list when the object is destroyed } @@ -166,13 +173,13 @@ namespace stim { } /// Project a set of samples onto the basis - void project(std::vector>& vlist, size_t nc) { + void project(std::vector >& vlist, size_t nc) { Sd.project(vlist, nc); Sc = Sd; } /// Calculate a density function from a list of points in spherical coordinates - void pdf(std::vector>& vlist, size_t nc) { + void pdf(std::vector >& vlist, size_t nc) { Sd.pdf(vlist, nc); Sc = Sd; } @@ -196,4 +203,4 @@ namespace stim { -#endif \ No newline at end of file +#endif -- libgit2 0.21.4