Commit eb119a4f8520f37c648a94435f159c025a4b2743
1 parent
9f5c0d4a
renamed fiber into centerline and removed all of the radius information from fiber.
Showing
1 changed file
with
0 additions
and
483 deletions
Show diff stats
stim/biomodels/fiber.h deleted
| 1 | -#ifndef STIM_FIBER_H | ||
| 2 | -#define STIM_FIBER_H | ||
| 3 | - | ||
| 4 | -#include <vector> | ||
| 5 | -#include <ANN/ANN.h> | ||
| 6 | - | ||
| 7 | -namespace stim{ | ||
| 8 | - | ||
| 9 | -/** This class stores information about a single fiber represented as a set of geometric points | ||
| 10 | - * between two branch or end points. This class is used as a fundamental component of the stim::network | ||
| 11 | - * class to describe an interconnected (often biological) network. | ||
| 12 | - */ | ||
| 13 | -template<typename T> | ||
| 14 | -class fiber{ | ||
| 15 | - | ||
| 16 | -protected: | ||
| 17 | - unsigned int N; //number of points in the fiber | ||
| 18 | - double **c; //centerline (array of double pointers) | ||
| 19 | - | ||
| 20 | - T* r; // array of fiber radii | ||
| 21 | - ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching | ||
| 22 | - | ||
| 23 | - /// Initialize an empty fiber | ||
| 24 | - void init() | ||
| 25 | - { | ||
| 26 | - kdt = NULL; | ||
| 27 | - c=NULL; | ||
| 28 | - r=NULL; | ||
| 29 | - N=0; | ||
| 30 | - } | ||
| 31 | - | ||
| 32 | - /// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0) | ||
| 33 | - void init(unsigned int n) | ||
| 34 | - { | ||
| 35 | - | ||
| 36 | - N = n; //set the number of points | ||
| 37 | - kdt = NULL; | ||
| 38 | - c = (double**) malloc(sizeof(double*) * N); //allocate the array pointer | ||
| 39 | - | ||
| 40 | - for(unsigned int i = 0; i < N; i++) //allocate space for each point | ||
| 41 | - c[i] = (double*) malloc(sizeof(double) * 3); | ||
| 42 | - | ||
| 43 | - r = (T*) malloc(sizeof(T) * N); //allocate space for the radii | ||
| 44 | - } | ||
| 45 | - | ||
| 46 | - /// Copies an existing fiber to the current fiber | ||
| 47 | - | ||
| 48 | - /// @param cpy stores the new copy of the fiber | ||
| 49 | - void copy( const stim::fiber<T>& cpy ){ | ||
| 50 | - | ||
| 51 | - ///allocate space for the new fiber | ||
| 52 | - init(cpy.N); | ||
| 53 | - | ||
| 54 | - ///copy the points | ||
| 55 | - for(unsigned int i = 0; i < N; i++){ | ||
| 56 | - for(unsigned int d = 0; d < 3; d++) //for each dimension | ||
| 57 | - c[i][d] = cpy.c[i][d]; //copy the coordinate | ||
| 58 | - | ||
| 59 | - r[i] = cpy.r[i]; //copy the radius | ||
| 60 | - } | ||
| 61 | - | ||
| 62 | - gen_kdtree(); //generate the kd tree for the new fiber | ||
| 63 | - } | ||
| 64 | - | ||
| 65 | - /// generate a KD tree for points on fiber | ||
| 66 | - void gen_kdtree() | ||
| 67 | - { | ||
| 68 | - int n_data = N; //create an array of data points | ||
| 69 | - ANNpointArray pts = (ANNpointArray)c; //cast the centerline list to an ANNpointArray | ||
| 70 | - kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree | ||
| 71 | - } | ||
| 72 | - | ||
| 73 | - /// find distance between two points | ||
| 74 | - double dist(double* p0, double* p1){ | ||
| 75 | - | ||
| 76 | - double sum = 0; // initialize variables | ||
| 77 | - float v; | ||
| 78 | - for(unsigned int d = 0; d < 3; d++) | ||
| 79 | - { | ||
| 80 | - v = p1[d] - p0[d]; | ||
| 81 | - sum +=v * v; | ||
| 82 | - | ||
| 83 | - } | ||
| 84 | - return sqrt(sum); | ||
| 85 | - } | ||
| 86 | - | ||
| 87 | - /// This function retreives the index for the fiber point closest to q | ||
| 88 | - | ||
| 89 | - /// @param q is a reference point used to find the closest point on the fiber center line | ||
| 90 | - unsigned int ann( stim::vec<double> q ){ | ||
| 91 | - | ||
| 92 | - ANNidxArray idx = new ANNidx[1]; //variable used to hold the nearest point | ||
| 93 | - ANNdistArray sq_dist = new ANNdist[1]; //variable used to hold the squared distance to the nearest point | ||
| 94 | - | ||
| 95 | - kdt->annkSearch(q.data(), 1, idx, sq_dist); //search the KD tree for the nearest neighbor | ||
| 96 | - | ||
| 97 | - return *idx; | ||
| 98 | - } | ||
| 99 | - | ||
| 100 | - /// Returns a stim::vec representing the point at index i | ||
| 101 | - | ||
| 102 | - /// @param i is an index of the desired centerline point | ||
| 103 | - stim::vec<T> get_vec(unsigned i){ | ||
| 104 | - stim::vec<T> r; | ||
| 105 | - r.resize(3); | ||
| 106 | - r[0] = c[i][0]; | ||
| 107 | - r[1] = c[i][1]; | ||
| 108 | - r[2] = c[i][2]; | ||
| 109 | - | ||
| 110 | - return r; | ||
| 111 | - } | ||
| 112 | - | ||
| 113 | - | ||
| 114 | -public: | ||
| 115 | - | ||
| 116 | - fiber(){ | ||
| 117 | - init(); | ||
| 118 | - } | ||
| 119 | - | ||
| 120 | - /// Copy constructor | ||
| 121 | - fiber(const stim::fiber<T> &obj){ | ||
| 122 | - | ||
| 123 | - copy(obj); | ||
| 124 | - | ||
| 125 | - } | ||
| 126 | - | ||
| 127 | - //temp constructor for graph visualization | ||
| 128 | - fiber(int n) | ||
| 129 | - { | ||
| 130 | - init(n); | ||
| 131 | - } | ||
| 132 | - | ||
| 133 | - /// Constructor takes a list of stim::vec points, the radius at each point is set to zero | ||
| 134 | - fiber(std::vector< stim::vec<T> > p){ | ||
| 135 | - init(p.size()); //initialize the fiber | ||
| 136 | - | ||
| 137 | - //for each point, set the centerline position and radius | ||
| 138 | - for(unsigned int i = 0; i < N; i++){ | ||
| 139 | - | ||
| 140 | - //set the centerline position | ||
| 141 | - for(unsigned int d = 0; d < 3; d++) | ||
| 142 | - c[i][d] = (double) p[i][d]; | ||
| 143 | - | ||
| 144 | - //set the radius | ||
| 145 | - r[i] = 0; | ||
| 146 | - } | ||
| 147 | - | ||
| 148 | - //generate a kd tree | ||
| 149 | - gen_kdtree(); | ||
| 150 | - } | ||
| 151 | - | ||
| 152 | - /// constructor takes a list of points and radii | ||
| 153 | - fiber(std::vector< stim::vec< T > > pos, std::vector< T > radii){ | ||
| 154 | - init(pos.size()); //initialize the fiber | ||
| 155 | - | ||
| 156 | - //for each point, set the centerline position and radius | ||
| 157 | - for(unsigned int i = 0; i < N; i++){ | ||
| 158 | - | ||
| 159 | - //set the centerline position | ||
| 160 | - for(unsigned int d = 0; d < 3; d++) | ||
| 161 | - c[i][d] = (double) pos[i][d]; | ||
| 162 | - | ||
| 163 | - //set the radius | ||
| 164 | - r[i] = radii[i]; | ||
| 165 | - } | ||
| 166 | - | ||
| 167 | - //generate a kd tree | ||
| 168 | - gen_kdtree(); | ||
| 169 | - } | ||
| 170 | - | ||
| 171 | - /// constructor takes an array of points and radii | ||
| 172 | - // this function is used when the radii are represented as a stim::vec, | ||
| 173 | - // since this may be easier when importing OBJs | ||
| 174 | - fiber(std::vector< stim::vec<T> > pos, std::vector< stim::vec<T> > radii){ | ||
| 175 | - | ||
| 176 | - init(pos.size()); | ||
| 177 | - | ||
| 178 | - //for each point, set the position and radius | ||
| 179 | - for(unsigned int i = 0; i < N; i++){ | ||
| 180 | - //at(i) = (double*)malloc(sizeof(double) * 3); | ||
| 181 | - for(unsigned int d = 0; d < 3; d++) | ||
| 182 | - c[i][d] = (double) pos[i][d]; | ||
| 183 | - | ||
| 184 | - r[i] = radii[i][(unsigned int)0]; | ||
| 185 | - } | ||
| 186 | - | ||
| 187 | - gen_kdtree(); | ||
| 188 | - } | ||
| 189 | - | ||
| 190 | - /// Assignment operation | ||
| 191 | - fiber& operator=(const fiber &rhs){ | ||
| 192 | - | ||
| 193 | - if(this == &rhs) return *this; //test for and handle self-assignment | ||
| 194 | - | ||
| 195 | - copy(rhs); | ||
| 196 | - } | ||
| 197 | - | ||
| 198 | - /// Calculate the length of the fiber and return it. | ||
| 199 | - double length(){ | ||
| 200 | - | ||
| 201 | - double* p0; | ||
| 202 | - double *p1; | ||
| 203 | - double l = 0; //initialize the length to zero | ||
| 204 | - | ||
| 205 | - //for each point | ||
| 206 | - //typename std::list< point<T> >::iterator i; //create a point iterator | ||
| 207 | - for(unsigned int i = 0; i < N; i++){ //for each point in the fiber | ||
| 208 | - | ||
| 209 | - if(i == 0) //if this is the first point, just store it | ||
| 210 | - p1 = c[0]; | ||
| 211 | - else{ //if this is any other point | ||
| 212 | - p0 = p1; //shift p1->p0 | ||
| 213 | - p1 = c[i]; //set p1 to the new point | ||
| 214 | - l += dist(p0, p1); //add the length of p1 - p0 to the running sum | ||
| 215 | - } | ||
| 216 | - } | ||
| 217 | - | ||
| 218 | - return l; //return the length | ||
| 219 | - } | ||
| 220 | - | ||
| 221 | - /// Calculates the length and average radius of the fiber | ||
| 222 | - | ||
| 223 | - /// @param length is filled with the fiber length | ||
| 224 | - T radius(T& length){ | ||
| 225 | - | ||
| 226 | - double* p0; //temporary variables to store point positions | ||
| 227 | - double* p1; | ||
| 228 | - T r0, r1; //temporary variables to store radii at points | ||
| 229 | - double l; | ||
| 230 | - T r_mean; //temporary variable to store the length and average radius of a fiber segment | ||
| 231 | - double length_sum = 0; //initialize the length to zero | ||
| 232 | - T radius_sum = 0; //initialize the radius sum to zero | ||
| 233 | - | ||
| 234 | - //for each point | ||
| 235 | - //typename std::list< point<T> >::iterator i; //create a point iterator | ||
| 236 | - for(unsigned int i = 0; i < N; i++){ //for each point in the fiber | ||
| 237 | - | ||
| 238 | - if(i == 0){ //if this is the first point, just store it | ||
| 239 | - p1 = c[0]; | ||
| 240 | - r1 = r[0]; | ||
| 241 | - } | ||
| 242 | - else{ //if this is any other point | ||
| 243 | - p0 = p1; //shift p1->p0 and r1->r0 | ||
| 244 | - r0 = r1; | ||
| 245 | - p1 = c[i]; //set p1 to the new point | ||
| 246 | - r1 = r[i]; | ||
| 247 | - | ||
| 248 | - l = dist(p0, p1); //calculate the length of the p0-p1 segment | ||
| 249 | - r_mean = (r0 + r1) / 2; //calculate the average radius of the segment | ||
| 250 | - | ||
| 251 | - radius_sum += r_mean * (T) l; //add the radius scaled by the length to a running sum | ||
| 252 | - length_sum += l; //add the length of p1 - p0 to the running sum | ||
| 253 | - } | ||
| 254 | - } | ||
| 255 | - | ||
| 256 | - length = length_sum; //store the total length | ||
| 257 | - | ||
| 258 | - //if the total length is zero, store a radius of zero | ||
| 259 | - if(length == 0) | ||
| 260 | - return 0; | ||
| 261 | - else | ||
| 262 | - return (T)(radius_sum / length); //return the average radius of the fiber | ||
| 263 | - } | ||
| 264 | - T average_radius() | ||
| 265 | - { | ||
| 266 | - T r_sum = 0.; | ||
| 267 | - for(unsigned int i = 0; i < N; i++) | ||
| 268 | - { | ||
| 269 | - r_sum = r_sum + r[i]; | ||
| 270 | - } | ||
| 271 | - return r_sum/((T) N); | ||
| 272 | - } | ||
| 273 | - | ||
| 274 | - /// Calculates the average radius of the fiber | ||
| 275 | - T radius(){ | ||
| 276 | - T length; | ||
| 277 | - return radius(length); | ||
| 278 | - } | ||
| 279 | - | ||
| 280 | - /// Returns the radius at index idx. | ||
| 281 | - T radius(int idx){ | ||
| 282 | - return r[idx]; | ||
| 283 | - } | ||
| 284 | - | ||
| 285 | - /// Return the point on the fiber closest to q | ||
| 286 | - /// @param q is the query point used to locate the nearest point on the fiber centerline | ||
| 287 | - stim::vec<T> nearest(stim::vec<T> q){ | ||
| 288 | - | ||
| 289 | - stim::vec<double> temp( (double) q[0], (double) q[1], (double) q[2]); | ||
| 290 | - | ||
| 291 | - unsigned int idx = ann(temp); //determine the index of the nearest neighbor | ||
| 292 | - | ||
| 293 | - return stim::vec<T>((T) c[idx][0], (T) c[idx][1], (T) c[idx][2]); //return the nearest centerline point | ||
| 294 | - } | ||
| 295 | - | ||
| 296 | - /// Return the point index on the fiber closest to q | ||
| 297 | - /// @param q is the query point used to locate the nearest point on the fiber centerline | ||
| 298 | - unsigned int nearest_idx(stim::vec<T> q){ | ||
| 299 | - | ||
| 300 | - stim::vec<double> temp((double) q[0], (double) q[1], (double) q[2]); | ||
| 301 | - | ||
| 302 | - unsigned int idx = ann(temp); //determine the index of the nearest neighbor | ||
| 303 | - | ||
| 304 | - return idx; //return the nearest centerline point index | ||
| 305 | - } | ||
| 306 | - | ||
| 307 | - /// Returns the fiber centerline as an array of stim::vec points | ||
| 308 | - std::vector< stim::vec<T> > centerline(){ | ||
| 309 | - | ||
| 310 | - //create an array of stim vectors | ||
| 311 | - std::vector< stim::vec<T> > pts(N); | ||
| 312 | - | ||
| 313 | - //cast each point to a stim::vec, keeping only the position information | ||
| 314 | - for(unsigned int i = 0; i < N; i++) | ||
| 315 | - pts[i] = stim::vec<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]); | ||
| 316 | - | ||
| 317 | - //return the centerline array | ||
| 318 | - return pts; | ||
| 319 | - } | ||
| 320 | - | ||
| 321 | - /// Returns the fiber centerline magnitudes as an array of stim::vec points | ||
| 322 | - std::vector< stim::vec<T> > centerlinemag(){ | ||
| 323 | - | ||
| 324 | - //create an array of stim vectors | ||
| 325 | - std::vector< stim::vec<T> > pts(N); | ||
| 326 | - | ||
| 327 | - //cast each point to a stim::vec, keeping only the position information | ||
| 328 | - for(unsigned int i = 0; i < N; i++) | ||
| 329 | - pts[i] = stim::vec<T>(r[i], r[i]);; | ||
| 330 | - | ||
| 331 | - //return the centerline array | ||
| 332 | - return pts; | ||
| 333 | - } | ||
| 334 | - | ||
| 335 | - /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned | ||
| 336 | - std::vector< stim::fiber<T> > split(unsigned int idx){ | ||
| 337 | - | ||
| 338 | - std::vector< stim::fiber<T> > fl; //create an array to store up to two fibers | ||
| 339 | - | ||
| 340 | - //if the index is an end point, only the existing fiber is returned | ||
| 341 | - if(idx == 0 || idx == N-1){ | ||
| 342 | - fl.resize(1); //set the size of the fiber to 1 | ||
| 343 | - fl[0] = *this; //copy the current fiber | ||
| 344 | - } | ||
| 345 | - | ||
| 346 | - //if the index is not an end point | ||
| 347 | - else{ | ||
| 348 | - | ||
| 349 | - unsigned int N1 = idx + 1; //calculate the size of both fibers | ||
| 350 | - unsigned int N2 = N - idx; | ||
| 351 | - | ||
| 352 | - fl.resize(2); //set the array size to 2 | ||
| 353 | - | ||
| 354 | - fl[0].init(N1); //set the size of each fiber | ||
| 355 | - fl[1].init(N2); | ||
| 356 | - | ||
| 357 | - //copy both halves of the fiber | ||
| 358 | - unsigned int i, d; | ||
| 359 | - | ||
| 360 | - //first half | ||
| 361 | - for(i = 0; i < N1; i++){ //for each centerline point | ||
| 362 | - for(d = 0; d < 3; d++) | ||
| 363 | - fl[0].c[i][d] = c[i][d]; //copy each coordinate | ||
| 364 | - fl[0].r[i] = r[i]; //copy the corresponding radius | ||
| 365 | - } | ||
| 366 | - | ||
| 367 | - //second half | ||
| 368 | - for(i = 0; i < N2; i++){ | ||
| 369 | - for(d = 0; d < 3; d++) | ||
| 370 | - fl[1].c[i][d] = c[idx + i][d]; | ||
| 371 | - fl[1].r[i] = r[idx + i]; | ||
| 372 | - } | ||
| 373 | - } | ||
| 374 | - | ||
| 375 | - return fl; //return the array | ||
| 376 | - | ||
| 377 | - } | ||
| 378 | - | ||
| 379 | - /// Calculates the set of fibers resulting from a connection between the current fiber and a fiber f | ||
| 380 | - | ||
| 381 | - /// @param f is the fiber that will be connected to the current fiber | ||
| 382 | - std::vector< stim::fiber<T> > connect( stim::fiber<T> &f, double dist){ | ||
| 383 | - | ||
| 384 | - double min_dist; | ||
| 385 | - unsigned int idx0, idx1; | ||
| 386 | - | ||
| 387 | - //go through each point in the query fiber, looking for the indices for the closest points | ||
| 388 | - for(unsigned int i = 0; i < f.n_pts(); i++){ | ||
| 389 | - //Run through all points and find the index with the closest point, then partition the fiber and return two fibers. | ||
| 390 | - | ||
| 391 | - } | ||
| 392 | - | ||
| 393 | - | ||
| 394 | - | ||
| 395 | - } | ||
| 396 | - | ||
| 397 | - /// Outputs the fiber as a string | ||
| 398 | - std::string str(){ | ||
| 399 | - std::stringstream ss; | ||
| 400 | - | ||
| 401 | - //create an iterator for the point list | ||
| 402 | - //typename std::list< point<T> >::iterator i; | ||
| 403 | - for(unsigned int i = 0; i < N; i++){ | ||
| 404 | - ss<<" [ "; | ||
| 405 | - for(unsigned int d = 0; d < 3; d++){ | ||
| 406 | - ss<<c[i][d]<<" "; | ||
| 407 | - } | ||
| 408 | - ss<<"] r = "<<r[i]<<std::endl; | ||
| 409 | - } | ||
| 410 | - | ||
| 411 | - return ss.str(); | ||
| 412 | - } | ||
| 413 | - /// Returns the number of centerline points in the fiber | ||
| 414 | - unsigned int size(){ | ||
| 415 | - return N; | ||
| 416 | - } | ||
| 417 | - | ||
| 418 | - | ||
| 419 | - /// Bracket operator returns the element at index i | ||
| 420 | - | ||
| 421 | - /// @param i is the index of the element to be returned as a stim::vec | ||
| 422 | - stim::vec<T> operator[](unsigned i){ | ||
| 423 | - return get_vec(i); | ||
| 424 | - } | ||
| 425 | - | ||
| 426 | - /// Back method returns the last point in the fiber | ||
| 427 | - stim::vec<T> back(){ | ||
| 428 | - return get_vec(N-1); | ||
| 429 | - } | ||
| 430 | - ////resample a fiber in the network | ||
| 431 | - stim::fiber<T> resample(T spacing) | ||
| 432 | - { | ||
| 433 | - std::cout<<"fiber::resample()"<<std::endl; | ||
| 434 | - | ||
| 435 | - std::vector<T> v(3); //v-direction vector of the segment | ||
| 436 | - stim::vec<T> p(3); //- intermediate point to be added | ||
| 437 | - stim::vec<T> p1(3); // p1 - starting point of an segment on the fiber, | ||
| 438 | - stim::vec<T> p2(3); // p2 - ending point, | ||
| 439 | - double sum=0; //distance summation | ||
| 440 | - std::vector<stim::vec<T> > fiberPositions = centerline(); | ||
| 441 | - std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber | ||
| 442 | - // for each point on the centerline (skip if it is the last point on centerline) | ||
| 443 | - //unsigned int N = fiberPositions.size(); // number of points on the fiber | ||
| 444 | - for(unsigned int f=0; f< N-1; f++) | ||
| 445 | - { | ||
| 446 | - | ||
| 447 | - p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1; | ||
| 448 | - for(unsigned int d = 0; d < 3; d++){ | ||
| 449 | - sum +=v[d] * v[d];} //length of segment-distance between starting and ending point | ||
| 450 | - | ||
| 451 | - T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment | ||
| 452 | - | ||
| 453 | - if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample | ||
| 454 | - { | ||
| 455 | - // repeat resampling until accumulated stepsize is equsl to length of the segment | ||
| 456 | - for(T step=0.0; step<lengthSegment; step+=spacing) | ||
| 457 | - { | ||
| 458 | - // calculate the resampled point by travelling step size in the direction of normalized gradient vector | ||
| 459 | - for(unsigned int i=0; i<3;i++) | ||
| 460 | - { | ||
| 461 | - p[i] = p1[i] + v[i]*(step/lengthSegment); | ||
| 462 | - } //for each dimension | ||
| 463 | - // add this resampled points to the new fiber list | ||
| 464 | - newPointList.push_back(p); | ||
| 465 | - } | ||
| 466 | - } | ||
| 467 | - 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 | ||
| 468 | - newPointList.push_back(fiberPositions[f+1]); | ||
| 469 | - } | ||
| 470 | - newPointList.push_back(fiberPositions[N-1]); //add the last point on the fiber to the new fiber list | ||
| 471 | - fiber newFiber(newPointList); | ||
| 472 | - return newFiber; | ||
| 473 | - } | ||
| 474 | - | ||
| 475 | -}; | ||
| 476 | - | ||
| 477 | - | ||
| 478 | - | ||
| 479 | -} //end namespace stim | ||
| 480 | - | ||
| 481 | - | ||
| 482 | - | ||
| 483 | -#endif |