codebase / stimlib

Browse Code »

Commit 7f1ab3c201bf09a8ff02f8bc66109d5be3a5252d

Authored by Pavel Govyadinov 2018-02-14 16:45:40 -0600

1 parent b7f8c759

fixed problems with dynamic network class

Showing 11 changed files with 2327 additions and 2782 deletions Show diff stats

Inline Side-by-side

stim/biomodels/centerline.h

Wrap text Show/Hide comments View file @7f1ab3c

1		-/*
2		-Copyright <2017> <David Mayerich>
3		-
4		-Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
5		-
6		-The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
7		-
8		-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
9		-*/
10		-
11		-#ifndef STIM_CENTERLINE_H
12		-#define STIM_CENTERLINE_H
13		-
14		-#include <vector>
15		-#include <stim/math/vec3.h>
16		-#include <stim/structures/kdtree.cuh>
17		-
18		-namespace stim{
19		-
20		-/** This class stores information about a single fiber represented as a set of geometric points
21		- * between two branch or end points. This class is used as a fundamental component of the stim::network
22		- * class to describe an interconnected (often biological) network.
23		- */
24		-template<typename T>
25		-class centerline : public std::vector< stim::vec3<T> >{
26		-
27		-protected:
28		-
29		- std::vector<T> L; //stores the integrated length along the fiber (used for parameterization)
30		-
31		- ///Return the normalized direction vector at point i (average of the incoming and outgoing directions)
32		- vec3<T> d(size_t i) {
33		- if (size() <= 1) return vec3<T>(0, 0, 0); //if there is insufficient information to calculate the direction, return a null vector
34		- 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
35		- if (i == 0) return (at(1) - at(0)).norm(); //the first direction vector is oriented towards the first line segment
36		- if (i == size() - 1) return (at(size() - 1) - at(size() - 2)).norm(); //the last direction vector is oriented towards the last line segment
37		-
38		- //all other direction vectors are the average direction of the two joined line segments
39		- vec3<T> a = at(i) - at(i - 1);
40		- vec3<T> b = at(i + 1) - at(i);
41		- vec3<T> ab = a.norm() + b.norm();
42		- return ab.norm();
43		- }
44		-
45		- //initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
46		- void update_L(size_t start = 0) {
47		- L.resize(size()); //allocate space for the L array
48		- if (start == 0) {
49		- L[0] = 0; //initialize the length value for the first point to zero (0)
50		- start++;
51		- }
52		-
53		- stim::vec3<T> d;
54		- for (size_t i = start; i < size(); i++) { //for each line segment in the centerline
55		- d = at(i) - at(i - 1);
56		- L[i] = L[i - 1] + d.len(); //calculate the running length total
57		- }
58		- }
59		-
60		- void init() {
61		- if (size() == 0) return; //return if there aren't any points
62		- update_L();
63		- }
64		-
65		- /// Returns a stim::vec representing the point at index i
66		-
67		- /// @param i is an index of the desired centerline point
68		- stim::vec<T> get_vec(unsigned i){
69		- return std::vector< stim::vec3<T> >::at(i);
70		- }
71		-
72		- ///finds the index of the point closest to the length l on the lower bound.
73		- ///binary search.
74		- size_t findIdx(T l) {
75		- for (size_t i = 1; i < L.size(); i++) { //for each point in the centerline
76		- if (L[i] > l) return i - 1; //if we have passed the desired length value, return i-1
77		- }
78		- return L.size() - 1;
79		- /*size_t i = L.size() / 2;
80		- size_t max = L.size() - 1;
81		- size_t min = 0;
82		- while (i < L.size() - 1){
83		- if (l < L[i]) {
84		- max = i;
85		- i = min + (max - min) / 2;
86		- }
87		- else if (L[i] <= l && L[i + 1] >= l) {
88		- break;
89		- }
90		- else {
91		- min = i;
92		- i = min + (max - min) / 2;
93		- }
94		- }
95		- return i;*/
96		- }
97		-
98		- ///Returns a position vector at the given length into the fiber (based on the pvalue).
99		- ///Interpolates the radius along the line.
100		- ///@param l: the location of the in the cylinder.
101		- ///@param idx: integer location of the point closest to l but prior to it.
102		- stim::vec3<T> p(T l, int idx) {
103		- T rat = (l - L[idx]) / (L[idx + 1] - L[idx]);
104		- stim::vec3<T> v1 = at(idx);
105		- stim::vec3<T> v2 = at(idx + 1);
106		- return(v1 + (v2 - v1)*rat);
107		- }
108		-
109		-
110		-public:
111		-
112		- using std::vector< stim::vec3<T> >::at;
113		- using std::vector< stim::vec3<T> >::size;
114		-
115		- centerline() : std::vector< stim::vec3<T> >() {
116		- init();
117		- }
118		- centerline(size_t n) : std::vector< stim::vec3<T> >(n){
119		- init();
120		- }
121		- centerline(std::vector<stim::vec3<T> > pos) :
122		- std::vector<stim::vec3<T> > (pos)
123		- {
124		- init();
125		- }
126		-
127		- //overload the push_back function to update the length vector
128		- void push_back(stim::vec3<T> p) {
129		- std::vector< stim::vec3<T> >::push_back(p);
130		- update_L(size() - 1);
131		- }
132		-
133		- ///Returns a position vector at the given p-value (p value ranges from 0 to 1).
134		- ///interpolates the position along the line.
135		- ///@param pvalue: the location of the in the cylinder, from 0 (beginning to 1).
136		- stim::vec3<T> p(T pvalue) {
137		- if (pvalue <= 0.0) return at(0); //return the first element
138		- if (pvalue >= 1.0) return back(); //return the last element
139		-
140		- T l = pvalue*L[L.size() - 1];
141		- int idx = findIdx(l);
142		- return p(l, idx);
143		- }
144		-
145		- ///Update centerline internal parameters (currently the L vector)
146		- void update() {
147		- init();
148		- }
149		- ///Return the length of the entire centerline
150		- T length() {
151		- return L.back();
152		- }
153		-
154		-
155		- /// stitch two centerlines
156		- ///@param c1, c2: two centerlines
157		- ///@param sigma: sample rate
158		- static std::vector< stim::centerline<T> > stitch(stim::centerline<T> c1, stim::centerline<T> c2 = stim::centerline<T>()) {
159		-
160		- std::vector< stim::centerline<T> > result;
161		- stim::centerline<T> new_centerline;
162		- stim::vec3<T> new_vertex;
163		-
164		- // if only one centerline, stitch itself!
165		- if (c2.size() == 0) {
166		- size_t num = c1.size();
167		- size_t id = 100000; // store the downsample start position
168		- T threshold;
169		- if (num < 4) { // if the number of vertex is less than 4, do nothing
170		- result.push_back(c1);
171		- return result;
172		- }
173		- else {
174		- // test geometry start vertex
175		- stim::vec3<T> v1 = c1[1] - c1[0]; // vector from c1[0] to c1[1]
176		- for (size_t p = 2; p < num; p++) { // 90° standard???
177		- stim::vec3<T> v2 = c1[p] - c1[0];
178		- float cosine = v2.dot(v1);
179		- if (cosine < 0) {
180		- id = p;
181		- threshold = v2.len();
182		- break;
183		- }
184		- }
185		- if (id != 100000) { // find a downsample position on the centerline
186		- T* c;
187		- c = (T)malloc(sizeof(T) (num - id) * 3);
188		- for (size_t p = id; p < num; p++) {
189		- for (size_t d = 0; d < 3; d++) {
190		- c[p * 3 + d] = c1[p][d];
191		- }
192		- }
193		- stim::kdtree<T, 3> kdt;
194		- kdt.create(c, num - id, 5); // create tree
195		-
196		- T* query = (T)malloc(sizeof(T) 3);
197		- for (size_t d = 0; d < 3; d++)
198		- query[d] = c1[0][d];
199		- size_t index;
200		- T dist;
201		-
202		- kdt.search(query, 1, &index, &dist);
203		-
204		- free(query);
205		- free(c);
206		-
207		- if (dist > threshold) {
208		- result.push_back(c1);
209		- }
210		- else {
211		- // the loop part
212		- new_vertex = c1[index];
213		- new_centerline.push_back(new_vertex);
214		- for (size_t p = 0; p < index + 1; p++) {
215		- new_vertex = c1[p];
216		- new_centerline.push_back(new_vertex);
217		- }
218		- result.push_back(new_centerline);
219		- new_centerline.clear();
220		-
221		- // the tail part
222		- for (size_t p = index; p < num; p++) {
223		- new_vertex = c1[p];
224		- new_centerline.push_back(new_vertex);
225		- }
226		- result.push_back(new_centerline);
227		- }
228		- }
229		- else { // there is one potential problem that two positions have to be stitched
230		- // test geometry end vertex
231		- stim::vec3<T> v1 = c1[num - 2] - c1[num - 1];
232		- for (size_t p = num - 2; p > 0; p--) { // 90° standard
233		- stim::vec3<T> v2 = c1[p - 1] - c1[num - 1];
234		- float cosine = v2.dot(v1);
235		- if (cosine < 0) {
236		- id = p;
237		- threshold = v2.len();
238		- break;
239		- }
240		- }
241		- if (id != 100000) { // find a downsample position
242		- T* c;
243		- c = (T)malloc(sizeof(T) (id + 1) * 3);
244		- for (size_t p = 0; p < id + 1; p++) {
245		- for (size_t d = 0; d < 3; d++) {
246		- c[p * 3 + d] = c1[p][d];
247		- }
248		- }
249		- stim::kdtree<T, 3> kdt;
250		- kdt.create(c, id + 1, 5); // create tree
251		-
252		- T* query = (T)malloc(sizeof(T) 1 * 3);
253		- for (size_t d = 0; d < 3; d++)
254		- query[d] = c1[num - 1][d];
255		- size_t index;
256		- T dist;
257		-
258		- kdt.search(query, 1, &index, &dist);
259		-
260		- free(query);
261		- free(c);
262		-
263		- if (dist > threshold) {
264		- result.push_back(c1);
265		- }
266		- else {
267		- // the tail part
268		- for (size_t p = 0; p < index + 1; p++) {
269		- new_vertex = c1[p];
270		- new_centerline.push_back(new_vertex);
271		- }
272		- result.push_back(new_centerline);
273		- new_centerline.clear();
274		-
275		- // the loop part
276		- for (size_t p = index; p < num; p++) {
277		- new_vertex = c1[p];
278		- new_centerline.push_back(new_vertex);
279		- }
280		- new_vertex = c1[index];
281		- new_centerline.push_back(new_vertex);
282		- result.push_back(new_centerline);
283		- }
284		- }
285		- else { // no stitch position
286		- result.push_back(c1);
287		- }
288		- }
289		- }
290		- }
291		-
292		-
293		- // two centerlines
294		- else {
295		- // find stitch position based on nearest neighbors
296		- size_t num1 = c1.size();
297		- T* c = (T)malloc(sizeof(T) num1 * 3); // c1 as reference point
298		- for (size_t p = 0; p < num1; p++) // centerline to array
299		- for (size_t d = 0; d < 3; d++) // because right now my kdtree code is a relatively close code, it has its own structure
300		- c[p * 3 + d] = c1[p][d]; // I will merge it into stimlib totally in the near future
301		-
302		- stim::kdtree<T, 3> kdt; // kdtree object
303		- kdt.create(c, num1, 5); // create tree
304		-
305		- size_t num2 = c2.size();
306		- T* query = (T)malloc(sizeof(T) num2 * 3); // c2 as query point
307		- for (size_t p = 0; p < num2; p++) {
308		- for (size_t d = 0; d < 3; d++) {
309		- query[p * 3 + d] = c2[p][d];
310		- }
311		- }
312		- std::vector<size_t> index(num2);
313		- std::vector<T> dist(num2);
314		-
315		- kdt.search(query, num2, &index[0], &dist[0]); // find the nearest neighbors in c1 for c2
316		-
317		- // clear up
318		- free(query);
319		- free(c);
320		-
321		- // find the average vertex distance of one centerline
322		- T sigma1 = 0;
323		- T sigma2 = 0;
324		- for (size_t p = 0; p < num1 - 1; p++)
325		- sigma1 += (c1[p] - c1[p + 1]).len();
326		- for (size_t p = 0; p < num2 - 1; p++)
327		- sigma2 += (c2[p] - c2[p + 1]).len();
328		- sigma1 /= (num1 - 1);
329		- sigma2 /= (num2 - 1);
330		- float threshold = 4 * (sigma1 + sigma2) / 2; // better way to do this?
331		-
332		- T min_d = *std::min_element(dist.begin(), dist.end()); // find the minimum distance between c1 and c2
333		-
334		- if (min_d > threshold) { // if the minimum distance is too large
335		- result.push_back(c1);
336		- result.push_back(c2);
337		-
338		-#ifdef DEBUG
339		- std::cout << "The distance between these two centerlines is too large" << std::endl;
340		-#endif
341		- }
342		- else {
343		- // auto smallest = std::min_element(dist.begin(), dist.end());
344		- unsigned int smallest = std::min_element(dist.begin(), dist.end());
345		- // auto i = std::distance(dist.begin(), smallest); // find the index of min-distance in distance list
346		- unsigned int i = std::distance(dist.begin(), smallest); // find the index of min-distance in distance list
347		-
348		- // stitch position in c1 and c2
349		- int id1 = index[i];
350		- int id2 = i;
351		-
352		- // actually there are two cases
353		- // first one inacceptable
354		- // second one acceptable
355		- if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) { // only stitch one end vertex to another centerline
356		- result.push_back(c1);
357		- result.push_back(c2);
358		- }
359		- else {
360		- if (id1 == 0 \|\| id1 == num1 - 1) { // if the stitch vertex is the first or last vertex of c1
361		- // for c2, consider two cases(one degenerate case)
362		- if (id2 == 0 \|\| id2 == num2 - 1) { // case 1, if stitch position is also on the end of c2
363		- // we have to decide which centerline get a new vertex, based on direction
364		- // for c1, computer the direction change angle
365		- stim::vec3<T> v1, v2;
366		- float alpha1, alpha2; // direction change angle
367		- if (id1 == 0)
368		- v1 = (c1[1] - c1[0]).norm();
369		- else
370		- v1 = (c1[num1 - 2] - c1[num1 - 1]).norm();
371		- v2 = (c2[id2] - c1[id1]).norm();
372		- alpha1 = v1.dot(v2);
373		- if (id2 == 0)
374		- v1 = (c2[1] - c2[0]).norm();
375		- else
376		- v1 = (c2[num2 - 2] - c2[num2 - 1]).norm();
377		- v2 = (c1[id1] - c2[id2]).norm();
378		- alpha2 = v1.dot(v2);
379		- if (abs(alpha1) > abs(alpha2)) { // add the vertex to c1 in order to get smooth connection
380		- // push back c1
381		- if (id1 == 0) { // keep geometry information
382		- new_vertex = c2[id2];
383		- new_centerline.push_back(new_vertex);
384		- for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
385		- new_vertex = c1[p];
386		- new_centerline.push_back(new_vertex);
387		- }
388		- }
389		- else {
390		- for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
391		- new_vertex = c1[p];
392		- new_centerline.push_back(new_vertex);
393		- }
394		- new_vertex = c2[id2];
395		- new_centerline.push_back(new_vertex);
396		- }
397		- result.push_back(new_centerline);
398		- new_centerline.clear();
399		-
400		- // push back c2
401		- for (size_t p = 0; p < num2; p++) {
402		- new_vertex = c2[p];
403		- new_centerline.push_back(new_vertex);
404		- }
405		- result.push_back(new_centerline);
406		- }
407		- else { // add the vertex to c2 in order to get smooth connection
408		- // push back c1
409		- for (size_t p = 0; p < num1; p++) {
410		- new_vertex = c1[p];
411		- new_centerline.push_back(new_vertex);
412		- }
413		- result.push_back(new_centerline);
414		- new_centerline.clear();
415		-
416		- // push back c2
417		- if (id2 == 0) { // keep geometry information
418		- new_vertex = c1[id1];
419		- new_centerline.push_back(new_vertex);
420		- for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
421		- new_vertex = c2[p];
422		- new_centerline.push_back(new_vertex);
423		- }
424		- }
425		- else {
426		- for (size_t p = 0; p < num2; p++) { // stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
427		- new_vertex = c2[p];
428		- new_centerline.push_back(new_vertex);
429		- }
430		- new_vertex = c1[id1];
431		- new_centerline.push_back(new_vertex);
432		- }
433		- result.push_back(new_centerline);
434		- }
435		- }
436		- else { // case 2, the stitch position is on c2
437		- // push back c1
438		- if (id1 == 0) { // keep geometry information
439		- new_vertex = c2[id2];
440		- new_centerline.push_back(new_vertex);
441		- for (size_t p = 0; p < num1; p++) { // stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
442		- new_vertex = c1[p];
443		- new_centerline.push_back(new_vertex);
444		- }
445		- }
446		- else {
447		- for (size_t p = 0; p < num1; p++) { // geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2
448		- new_vertex = c1[p];
449		- new_centerline.push_back(new_vertex);
450		- }
451		- new_vertex = c2[id2];
452		- new_centerline.push_back(new_vertex);
453		- }
454		- result.push_back(new_centerline);
455		- new_centerline.clear();
456		-
457		- // push back c2
458		- for (size_t p = 0; p < id2 + 1; p++) { // first part
459		- new_vertex = c2[p];
460		- new_centerline.push_back(new_vertex);
461		- }
462		- result.push_back(new_centerline);
463		- new_centerline.clear();
464		-
465		- for (size_t p = id2; p < num2; p++) { // second part
466		- new_vertex = c2[p];
467		- new_centerline.push_back(new_vertex);
468		- }
469		- result.push_back(new_centerline);
470		- }
471		- }
472		- else { // if the stitch vertex is the first or last vertex of c2
473		- // push back c2
474		- if (id2 == 0) { // keep geometry information
475		- new_vertex = c1[id1];
476		- new_centerline.push_back(new_vertex);
477		- for (size_t p = 0; p < num2; p++) { // stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2
478		- new_vertex = c2[p];
479		- new_centerline.push_back(new_vertex);
480		- }
481		- }
482		- else {
483		- for (size_t p = 0; p < num2; p++) { // geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1
484		- new_vertex = c2[p];
485		- new_centerline.push_back(new_vertex);
486		- }
487		- new_vertex = c1[id1];
488		- new_centerline.push_back(new_vertex);
489		- result.push_back(new_centerline);
490		- new_centerline.clear();
491		-
492		- // push back c1
493		- for (size_t p = 0; p < id1 + 1; p++) { // first part
494		- new_vertex = c1[p];
495		- new_centerline.push_back(new_vertex);
496		- }
497		- result.push_back(new_centerline);
498		- new_centerline.clear();
499		-
500		- for (size_t p = id1; p < num1; p++) { // second part
501		- new_vertex = c1[p];
502		- new_centerline.push_back(new_vertex);
503		- }
504		- result.push_back(new_centerline);
505		- }
506		- }
507		- }
508		- }
509		- }
510		- return result;
511		- }
512		-
513		- /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
514		- std::vector< stim::centerline<T> > split(unsigned int idx){
515		-
516		- std::vector< stim::centerline<T> > fl; //create an array to store up to two fibers
517		- size_t N = size();
518		-
519		- //if the index is an end point, only the existing fiber is returned
520		- if(idx == 0 \|\| idx == N-1){
521		- fl.resize(1); //set the size of the fiber to 1
522		- fl[0] = *this; //copy the current fiber
523		- }
524		-
525		- //if the index is not an end point
526		- else{
527		-
528		- unsigned int N1 = idx + 1; //calculate the size of both fibers
529		- unsigned int N2 = N - idx;
530		-
531		- fl.resize(2); //set the array size to 2
532		-
533		- fl[0] = stim::centerline<T>(N1); //set the size of each fiber
534		- fl[1] = stim::centerline<T>(N2);
535		-
536		- //copy both halves of the fiber
537		- unsigned int i;
538		-
539		- //first half
540		- for(i = 0; i < N1; i++) //for each centerline point
541		- fl[0][i] = std::vector< stim::vec3<T> >::at(i);
542		- fl[0].init(); //initialize the length vector
543		-
544		- //second half
545		- for(i = 0; i < N2; i++)
546		- fl[1][i] = std::vector< stim::vec3<T> >::at(idx+i);
547		- fl[1].init(); //initialize the length vector
548		- }
549		-
550		- return fl; //return the array
551		-
552		- }
553		-
554		- /// Outputs the fiber as a string
555		- std::string str(){
556		- std::stringstream ss;
557		- size_t N = std::vector< stim::vec3<T> >::size();
558		- ss << "---------[" << N << "]---------" << std::endl;
559		- for (size_t i = 0; i < N; i++)
560		- ss << std::vector< stim::vec3<T> >::at(i) << std::endl;
561		- ss << "--------------------" << std::endl;
562		-
563		- return ss.str();
564		- }
565		-
566		- /// Back method returns the last point in the fiber
567		- stim::vec3<T> back(){
568		- return std::vector< stim::vec3<T> >::back();
569		- }
570		-
571		- ////resample a fiber in the network
572		- stim::centerline<T> resample(T spacing)
573		- {
574		- //std::cout<<"fiber::resample()"<<std::endl;
575		-
576		- stim::vec3<T> v; //v-direction vector of the segment
577		- stim::vec3<T> p; //- intermediate point to be added
578		- stim::vec3<T> p1; // p1 - starting point of an segment on the fiber,
579		- stim::vec3<T> p2; // p2 - ending point,
580		- //double sum=0; //distance summation
581		-
582		- size_t N = size();
583		-
584		- centerline<T> new_c; // initialize list of new resampled points on the fiber
585		- // for each point on the centerline (skip if it is the last point on centerline)
586		- for(unsigned int f=0; f< N-1; f++)
587		- {
588		- p1 = at(f);
589		- p2 = at(f+1);
590		- v = p2 - p1;
591		-
592		- T lengthSegment = v.len(); //find Length of the segment as distance between the starting and ending points of the segment
593		-
594		- if(lengthSegment >= spacing){ // if length of the segment is greater than standard deviation resample
595		-
596		- // repeat resampling until accumulated stepsize is equsl to length of the segment
597		- for(T step=0.0; step<lengthSegment; step+=spacing){
598		- // calculate the resampled point by travelling step size in the direction of normalized gradient vector
599		- p = p1 + v * (step / lengthSegment);
600		-
601		- // add this resampled points to the new fiber list
602		- new_c.push_back(p);
603		- }
604		- }
605		- 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
606		- new_c.push_back(at(f));
607		- }
608		- new_c.push_back(at(N-1)); //add the last point on the fiber to the new fiber list
609		- //centerline newFiber(newPointList);
610		- return new_c;
611		- }
612		-
613		-};
614		-
615		-
616		-
617		-} //end namespace stim
618		-
619		-
620		-
621		-#endif
	1	+#ifndef JACK_CENTERLINE_H
	2	+#define JACK_CENTERLINE_H
	3	+
	4	+#include <stdlib.h>
	5	+#include <vector>
	6	+#include <stim/math/vec3.h>
	7	+#include <stim/structures/kdtree.cuh>
	8	+
	9	+namespace stim {
	10	+
	11	+ /// we always assume that one centerline has a flow direction even it actually does not have. Also, we allow loop centerline
	12	+ /// NOTE: centerline is not derived from std::vector<stim::vec3<T>> class!!!
	13	+ template<typename T>
	14	+ class centerline {
	15	+
	16	+ private:
	17	+ size_t n; // number of points on the centerline, can be zero if NULL
	18	+
	19	+ // update length information at each point (distance from starting point) starting from index "start"
	20	+ void update_L(size_t start = 0) {
	21	+ L.resize(n);
	22	+
	23	+ if (start == 0) {
	24	+ L[0] = 0.0f;
	25	+ start++;
	26	+ }
	27	+
	28	+ stim::vec3<T> dir; // temp direction vector for calculating length between two points
	29	+ for (size_t i = start; i < n; i++) {
	30	+ dir = C[i] - C[i - 1]; // calculate the certerline extending direction
	31	+ L[i] = L[i - 1] + dir.len(); // addition
	32	+ }
	33	+ }
	34	+
	35	+ protected:
	36	+ std::vector<stim::vec3<T> > C; // points on the centerline
	37	+ std::vector<T> L; // stores the integrated length along the fiber
	38	+
	39	+ public:
	40	+ /// constructors
	41	+ // empty constructor
	42	+ centerline() {
	43	+ n = 0;
	44	+ }
	45	+
	46	+ // constructor that allocate memory
	47	+ centerline(size_t s) {
	48	+ n = s;
	49	+ C.resize(s); // allocate memory for points
	50	+ L.resize(s); // allocate memory for lengths
	51	+
	52	+ update_L();
	53	+ }
	54	+
	55	+ // constructor that constructs a centerline based on a list of points
	56	+ centerline(std::vector<stim::vec3<T> > rhs) {
	57	+ n = rhs.size(); // get the number of points
	58	+ C.resize(n);
	59	+ for (size_t i = 0; i < n; i++)
	60	+ C[i] = rhs[i]; // copy data
	61	+ update_L();
	62	+ }
	63	+
	64	+
	65	+ /// vector operations
	66	+ // add a new point to current centerline
	67	+ void push_back(stim::vec3<T> p) {
	68	+ C.push_back(p);
	69	+ n++; // increase the number of points
	70	+ update_L(n - 1);
	71	+ }
	72	+
	73	+ // insert a new point at specific location to current centerline
	74	+ void insert(typename std::vector<stim::vec3<T> >::iterator pos, stim::vec3<T> p) {
	75	+ C.insert(pos, p); // insert a new point
	76	+ n++;
	77	+ size_t d = std::distance(C.begin(), pos); // get the index
	78	+ update_L(d);
	79	+ }
	80	+ // insert a new point at C[idx]
	81	+ void insert(size_t idx, stim::vec3<T> p) {
	82	+ n++;
	83	+ C.resize(n);
	84	+ for (size_t i = n - 1; i > idx; i--) // move point location
	85	+ C[i] = C[i - 1];
	86	+ C[idx] = p;
	87	+ update_L(idx);
	88	+ }
	89	+
	90	+ // assign a point at specific location to current centerline
	91	+ void assign(size_t idx, stim::vec3<T> p) {
	92	+ C[idx] = p;
	93	+ update_L(idx);
	94	+ }
	95	+
	96	+ // erase a point at specific location on current centerline
	97	+ void erase(typename std::vector<stim::vec3<T> >::iterator pos) {
	98	+ C.erase(pos); // erase a point
	99	+ n--;
	100	+ size_t d = std::distance(C.begin(), pos); // get the index
	101	+ update_L(d);
	102	+ }
	103	+ // erase a point at C[idx]
	104	+ void erase(size_t idx) {
	105	+ n--;
	106	+ for (size_t i = idx; i < n; i++)
	107	+ C[i] = C[i + 1];
	108	+ C.resize(n);
	109	+ update_L(idx);
	110	+ }
	111	+
	112	+ // clear up all the points
	113	+ void clear() {
	114	+ C.clear(); // clear list
	115	+ L.clear(); // clear length information
	116	+ n = 0; // set number to zero
	117	+ }
	118	+
	119	+ // reverse current centerline in terms of points order
	120	+ centerline<T> reverse() {
	121	+ centerline<T> result = *this;
	122	+
	123	+ std::reverse(result.C.begin(), result.C.end());
	124	+ result.update_L();
	125	+
	126	+ return result;
	127	+ }
	128	+
	129	+ /// functions for reading centerline information
	130	+ // return the number of points on current centerline
	131	+ size_t size() {
	132	+ return n;
	133	+ }
	134	+
	135	+ // return the length
	136	+ T length() {
	137	+ return L.back();
	138	+ }
	139	+
	140	+ // finds the index of the point closest to the length "l" on the lower bound
	141	+ size_t findIdx(T l) {
	142	+ for (size_t i = 1; i < L.size(); i++) {
	143	+ if (L[i] > l)
	144	+ return i - 1;
	145	+ }
	146	+
	147	+ return L.size() - 1;
	148	+ }
	149	+
	150	+ // get a position vector at the given length into the fiber (based on the pvalue), interpolate
	151	+ stim::vec3<T> p(T l, size_t idx) {
	152	+ T rate = (l - L[idx]) / (L[idx + 1] - L[idx]);
	153	+ stim::vec3<T> v1 = C[idx];
	154	+ stim::vec3<T> v2 = C[idx + 1];
	155	+
	156	+ return (v1 + (v2 - v1) * rate);
	157	+ }
	158	+ // get a position vector at the given pvalue(pvalue[0.0f, 1.0f])
	159	+ stim::vec3<T> p(T pvalue) {
	160	+ // degenerated cases
	161	+ if (pvalue <= 0.0f) return C[0];
	162	+ if (pvalue >= 1.0f) return C.back();
	163	+
	164	+ T l = pvalue * L.back(); // get the length based on the given pvalue
	165	+ size_t idx = findIdx(l);
	166	+
	167	+ return p(l, idx); // interpolation
	168	+ }
	169	+
	170	+ // get the normalized direction vector at point idx (average of the incoming and outgoing directions)
	171	+ stim::vec3<T> d(size_t idx) {
	172	+ if (n <= 1) return stim::vec3<T>(0.0f, 0.0f, 0.0f); // if there is insufficient information to calculate the direction, return null
	173	+ if (n == 2) return (C[1] - C[0]).norm(); // if there are only two points, the direction vector at both is the direction of the line segment
	174	+
	175	+ // degenerate cases at two ends
	176	+ if (idx == 0) return (C[1] - C[0]).norm(); // the first direction vector is oriented towards the first line segment
	177	+ if (idx == n - 1) return (C[n - 1] - C[n - 2]).norm(); // the last direction vector is oriented towards the last line segment
	178	+
	179	+ // all other direction vectors are the average direction of the two joined line segments
	180	+ stim::vec3<T> a = C[idx] - C[idx - 1];
	181	+ stim::vec3<T> b = C[idx + 1] - C[idx];
	182	+ stim::vec3<T> ab = a.norm() + b.norm();
	183	+ return ab.norm();
	184	+ }
	185	+
	186	+
	187	+ /// arithmetic operations
	188	+ // '=' operation
	189	+ centerline<T> & operator=(centerline<T> rhs) {
	190	+ n = rhs.n;
	191	+ C = rhs.C;
	192	+ L = rhs.L;
	193	+
	194	+ return *this;
	195	+ }
	196	+
	197	+ // "[]" operation
	198	+ stim::vec3<T> & operator[](size_t idx) {
	199	+ return C[idx];
	200	+ }
	201	+
	202	+ // "==" operation
	203	+ bool operator==(centerline<T> rhs) const {
	204	+
	205	+ if (n != rhs.size())
	206	+ return false;
	207	+ else {
	208	+ size_t num = rhs.size();
	209	+ stim::vec3<T> tmp; // weird situation that I can only use tmp instead of C itself in comparison
	210	+ for (size_t i = 0; i < num; i++) {
	211	+ stim::vec3<T> tmp = C[i];
	212	+ if (tmp[0] != rhs[i][0] \|\| tmp[1] != rhs[i][1] \|\| tmp[2] != rhs[i][2])
	213	+ return false;
	214	+ }
	215	+ return true;
	216	+ }
	217	+ }
	218	+
	219	+ // "+" operation
	220	+ centerline<T> operator+(stim::vec3<T> p) const {
	221	+ centerline<T> result(*this);
	222	+ result.C.push_back(p);
	223	+ result.n++;
	224	+ result.update_L(n - 1);
	225	+
	226	+ return result;
	227	+ }
	228	+ centerline<T> operator+(centerline<T> c) const {
	229	+ centerline<T> result(*this);
	230	+ size_t num1 = result.size();
	231	+ size_t num2 = c.size();
	232	+ for (size_t i = 0; i < num2; i++)
	233	+ result.push_back(c[i]);
	234	+ result.update_L(num1);
	235	+
	236	+ return result;
	237	+ }
	238	+
	239	+
	240	+ /// advanced operation
	241	+ // stitch two centerlines if possible (mutual-stitch and self-stitch)
	242	+ static std::vector<centerline<T> > stitch(centerline<T> c1, centerline<T> c2, size_t end = 0) {
	243	+ std::vector<centerline<T> > result;
	244	+ centerline<T> new_centerline;
	245	+ stim::vec3<T> new_vertex;
	246	+ // ******** for Pavel ********
	247	+ // ******** JACK thinks that ultimately we want it AUTOMATEDLY! ********
	248	+
	249	+ // check stitch case
	250	+ if (c1 == c2) { // self-stitch case
	251	+ // *** don't know how it works ***
	252	+ }
	253	+ else { // mutual-stitch case
	254	+ size_t num1 = c1.size(); // get the numbers of two centerlines
	255	+ size_t num2 = c2.size();
	256	+
	257	+ T* reference = (T)malloc(sizeof(T) num1 * 3); // c1 as reference set
	258	+ T* query = (T)malloc(sizeof(T) num2 * 3); // c2 as query set
	259	+ for (size_t p = 0; p < num1; p++) // read points
	260	+ for (size_t d = 0; d < 3; d++)
	261	+ reference[p * 3 + d] = c1[p][d]; // KDTREE is stilla close code, it has its own structure
	262	+ for (size_t p = 0; p < num2; p++) // read points
	263	+ for (size_t d = 0; d < 3; d++)
	264	+ query[p * 3 + d] = c2[p][d];
	265	+
	266	+ stim::kdtree<T, 3> kdt;
	267	+ kdt.create(reference, num1, 5); // create a tree based on reference set, max_tree_level is set to 5
	268	+
	269	+ std::vector<size_t> index(num2); // stores index results
	270	+ std::vector<float> dist(num2);
	271	+
	272	+ kdt.search(query, num2, &index[0], &dist[0]); // find the nearest neighbor in c1 for c2
	273	+
	274	+ // clear up
	275	+ free(reference);
	276	+ free(query);
	277	+
	278	+ std::vector<float>::iterator sm = std::min_element(dist.begin(), dist.end()); // find smallest distance
	279	+ size_t id = std::distance(dist.begin(), sm); // find the target index
	280	+
	281	+ size_t id1 = index[id];
	282	+ size_t id2 = id;
	283	+
	284	+ // for centerline c1
	285	+ bool flag = false; // flag indicates whether it has already added the bifurcation point to corresponding new centerline
	286	+ if (id1 == 0 \|\| id1 == num1 - 1) { // splitting bifurcation is on the end
	287	+ new_centerline = c1;
	288	+ new_vertex = c2[id2];
	289	+ if (id1 == 0) {
	290	+ new_centerline.insert(0, new_vertex);
	291	+ flag = true;
	292	+ }
	293	+ if (id1 == num1 - 1) {
	294	+ new_centerline.push_back(new_vertex);
	295	+ flag = true;
	296	+ }
	297	+
	298	+ result.push_back(new_centerline);
	299	+ }
	300	+ else { // splitting bifurcation is on the centerline
	301	+ std::vector<centerline<T> > tmp_centerline = c1.split(id1);
	302	+ result = tmp_centerline;
	303	+ }
	304	+
	305	+ // for centerline c2
	306	+ if (id2 == 0 \|\| id2 == num2 - 1) { // splitting bidurcation is on the end
	307	+ new_centerline = c2;
	308	+ if (flag)
	309	+ result.push_back(new_centerline);
	310	+ else { // add the bifurcation point to this centerline
	311	+ new_vertex = c1[id1];
	312	+ if (id2 == 0) {
	313	+ new_centerline.insert(0, new_vertex);
	314	+ flag = true;
	315	+ }
	316	+ if (id2 == num2 - 1) {
	317	+ new_centerline.push_back(new_vertex);
	318	+ flag = true;
	319	+ }
	320	+
	321	+ result.push_back(new_centerline);
	322	+ }
	323	+ }
	324	+ else { // splitting bifurcation is on the centerline
	325	+ std::vector<centerline<T> > tmp_centerline = c2.split(id2);
	326	+ result.push_back(tmp_centerline[0]);
	327	+ result.push_back(tmp_centerline[1]);
	328	+ }
	329	+ }
	330	+
	331	+ return result;
	332	+ }
	333	+
	334	+ // split current centerline at specific position
	335	+ std::vector<centerline<T> > split(size_t idx) {
	336	+ std::vector<centerline<T> > result;
	337	+
	338	+ // won't split
	339	+ if (idx <= 0 \|\| idx >= n - 1) {
	340	+ result.resize(1);
	341	+ result[0] = *this; // return current centerline
	342	+ }
	343	+ // do split
	344	+ else {
	345	+ size_t n1 = idx + 1; // vertex idx would appear twice
	346	+ size_t n2 = n - idx; // in total n + 1 points
	347	+
	348	+ centerline<T> tmp; // temp centerline
	349	+
	350	+ result.resize(2);
	351	+
	352	+ for (size_t i = 0; i < n1; i++) // first half
	353	+ tmp.push_back(C[i]);
	354	+ tmp.update_L();
	355	+ result[0] = tmp;
	356	+ tmp.clear(); // clear up for next computation
	357	+
	358	+ for (size_t i = 0; i < n2; i++) // second half
	359	+ tmp.push_back(C[i + idx]);
	360	+ tmp.update_L();
	361	+ result[1] = tmp;
	362	+ }
	363	+
	364	+ return result;
	365	+ }
	366	+
	367	+ // resample current centerline
	368	+ centerline<T> resample(T spacing) {
	369	+
	370	+ stim::vec3<T> dir; // direction vector
	371	+ stim::vec3<T> tmp; // intermiate point to be added
	372	+ stim::vec3<T> p1; // starting point
	373	+ stim::vec3<T> p2; // ending point
	374	+
	375	+ centerline<T> result;
	376	+
	377	+ for (size_t i = 0; i < n - 1; i++) {
	378	+ p1 = C[i];
	379	+ p2 = C[i + 1];
	380	+
	381	+ dir = p2 - p1; // compute the direction of current segment
	382	+ T seg_len = dir.len();
	383	+
	384	+ if (seg_len > spacing) { // current segment can be sampled
	385	+ for (T step = 0.0f; step < seg_len; step += spacing) {
	386	+ tmp = p1 + dir * (step / seg_len); // add new point
	387	+ result.push_back(tmp);
	388	+ }
	389	+ }
	390	+ else
	391	+ result.push_back(p1); // push back starting point
	392	+ }
	393	+ result.push_back(p2); // push back ending point
	394	+
	395	+ return result;
	396	+ }
	397	+ };
	398	+}
	399	+
	400	+#endif
...	...

stim/biomodels/network.h

Wrap text Show/Hide comments View file @7f1ab3c

1		-/*
2		-Copyright <2017> <David Mayerich>
3		-
4		-Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
5		-
6		-The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
7		-
8		-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
9		-*/
10		-#ifndef STIM_NETWORK_H
11		-#define STIM_NETWORK_H
12		-
13		-#include <stdlib.h>
14		-#include <assert.h>
15		-#include <sstream>
16		-#include <fstream>
17		-#include <algorithm>
18		-#include <string.h>
19		-#include <math.h>
20		-#include <stim/math/vec3.h>
21		-#include <stim/visualization/obj.h>
22		-#include <stim/visualization/swc.h>
23		-#include <stim/visualization/cylinder.h>
24		-#include <stim/cuda/cudatools/timer.h>
25		-#include <stim/cuda/cudatools/callable.h>
26		-#include <stim/structures/kdtree.cuh>
27		-//******************help function******************
28		-// gaussian_function
29		-CUDA_CALLABLE float gaussianFunction(float x, float std = 25) { return exp(-x / (2 * std*std)); } // std default sigma value is 25
30		-
31		-// compute metric in parallel
32		-#ifdef __CUDACC__
33		-template <typename T>
34		-__global__ void find_metric_parallel(T* M, size_t n, T* D, float sigma){
35		- size_t x = blockDim.x * blockIdx.x + threadIdx.x;
36		- if(x >= n) return;
37		- M[x] = 1.0f - gaussianFunction(D[x], sigma);
38		-}
39		-
40		-//find the corresponding edge index from array index
41		-__global__ void find_edge_index_parallel(size_t* I, size_t n, unsigned* R, size_t* E, size_t ne){
42		- size_t x = blockDim.x * blockIdx.x + threadIdx.x;
43		- if(x >= n) return;
44		- unsigned i = 0;
45		- size_t N = 0;
46		- for(unsigned e = 0; e < ne; e++){
47		- N += E[e];
48		- if(I[x] < N){
49		- R[x] = i;
50		- break;
51		- }
52		- i++;
53		- }
54		-}
55		-#endif
56		-
57		-//hard-coded factor
58		-int threshold_fac;
59		-
60		-namespace stim{
61		-/** This is the a class that interfaces with gl_spider in order to store the currently
62		- * segmented network. The following data is stored and can be extracted:
63		- * 1)Network geometry and centerline.
64		- * 2)Network connectivity (a graph of nodes and edges), reconstructed using kdtree.
65		-*/
66		-
67		-template<typename T>
68		-class network{
69		-
70		- ///Each edge is a fiber with two nodes.
71		- ///Each node is an in index to the endpoint of the fiber in the nodes array.
72		- class edge : public cylinder<T>
73		- {
74		- public:
75		-
76		- unsigned int v[2]; //unique id's designating the starting and ending
77		- // default constructor
78		- edge() : cylinder<T>() {
79		- v[1] = (unsigned)(-1); v[0] = (unsigned)(-1);
80		- }
81		- /// Constructor - creates an edge from a list of points by calling the stim::fiber constructor
82		-/*
83		- ///@param v0, the starting index.
84		- ///@param v1, the ending index.
85		- ///@param sz, the number of point in the fiber.
86		- edge(unsigned int v0, unsigned int v1, unsigned int sz) : cylinder<T>(
87		- {
88		-
89		- }
90		-*/
91		- edge(std::vector<stim::vec3<T> > p, std::vector<T> s)
92		- : cylinder<T>(p,s)
93		- {
94		- }
95		- ///@param p is an array of positions in space
96		- edge(stim::centerline<T> p) : cylinder<T>(p){}
97		-
98		- /// Copy constructor creates an edge from a cylinder
99		- edge(stim::cylinder<T> f) : cylinder<T>(f) {}
100		-
101		- /// Resamples an edge by calling the fiber resampling function
102		- edge resample(T spacing){
103		- edge e(cylinder<T>::resample(spacing)); //call the fiber->edge constructor
104		- e.v[0] = v[0]; //copy the vertex data
105		- e.v[1] = v[1];
106		-
107		- return e; //return the new edge
108		- }
109		-
110		- /// Output the edge information as a string
111		- std::string str(){
112		- std::stringstream ss;
113		- ss<<"("<<cylinder<T>::size()<<")\tl = "<<this->length()<<"\t"<<v[0]<<"----"<<v[1];
114		- return ss.str();
115		- }
116		-
117		- std::vector<edge> split(unsigned int idx){
118		-
119		- std::vector< stim::cylinder<T> > C;
120		- C.resize(2);
121		- C = (*this).cylinder<T>::split(idx);
122		- std::vector<edge> E(C.size());
123		-
124		- for(unsigned e = 0; e < E.size(); e++){
125		- E[e] = C[e];
126		- }
127		- return E;
128		- }
129		-
130		- /// operator for writing the edge information into a binary .nwt file.
131		- friend std::ofstream& operator<<(std::ofstream& out, const edge& e)
132		- {
133		- out.write(reinterpret_cast<const char*>(&e.v[0]), sizeof(unsigned int)); ///write the starting point.
134		- out.write(reinterpret_cast<const char*>(&e.v[1]), sizeof(unsigned int)); ///write the ending point.
135		- unsigned int sz = e.size(); ///write the number of point in the edge.
136		- out.write(reinterpret_cast<const char*>(&sz), sizeof(unsigned int));
137		- for(int i = 0; i < sz; i++) ///write each point
138		- {
139		- stim::vec3<T> point = e[i];
140		- out.write(reinterpret_cast<const char>(&point[0]), 3sizeof(T));
141		- // for(int j = 0; j < nmags(); j++) //future code for multiple mags
142		- // {
143		- out.write(reinterpret_cast<const char*>(&e.R[i]), sizeof(T)); ///write the radius
144		- //std::cout << point.str() << " " << e.R[i] << std::endl;
145		- // }
146		- }
147		- return out; //return stream
148		- }
149		-
150		- /// operator for reading an edge from a binary .nwt file.
151		- friend std::ifstream& operator>>(std::ifstream& in, edge& e)
152		- {
153		- unsigned int v0, v1, sz;
154		- in.read(reinterpret_cast<char*>(&v0), sizeof(unsigned int)); //read the staring point.
155		- in.read(reinterpret_cast<char*>(&v1), sizeof(unsigned int)); //read the ending point
156		- in.read(reinterpret_cast<char*>(&sz), sizeof(unsigned int)); //read the number of points in the edge
157		-// stim::centerline<T> temp = stim::centerline<T>(sz); //allocate the new edge
158		-// e = edge(temp);
159		- std::vector<stim::vec3<T> > p(sz);
160		- std::vector<T> r(sz);
161		- for(int i = 0; i < sz; i++) //set the points and radii to the newly read values
162		- {
163		- stim::vec3<T> point;
164		- in.read(reinterpret_cast<char>(&point[0]), 3sizeof(T));
165		- p[i] = point;
166		- T mag;
167		- // for(int j = 0; j < nmags(); j++) ///future code for mags
168		- // {
169		- in.read(reinterpret_cast<char*>(&mag), sizeof(T));
170		- r[i] = mag;
171		- //std::cout << point.str() << " " << mag << std::endl;
172		- // }
173		- }
174		- e = edge(p,r);
175		- e.v[0] = v0; e.v[1] = v1;
176		- return in;
177		- }
178		- };
179		-
180		- ///Node class that stores the physical position of the node as well as the edges it is connected to (edges that connect to it), As well as any additional data necessary.
181		- class vertex : public stim::vec3<T>
182		- {
183		- public:
184		- //std::vector<unsigned int> edges; //indices of edges connected to this node.
185		- std::vector<unsigned int> e[2]; //indices of edges going out (e[0]) and coming in (e[1])
186		- //stim::vec3<T> p; //position of this node in physical space.
187		- //default constructor
188		- vertex() : stim::vec3<T>()
189		- {
190		- }
191		- //constructor takes a stim::vec
192		- vertex(stim::vec3<T> p) : stim::vec3<T>(p){}
193		-
194		- /// Output the vertex information as a string
195		- std::string
196		- str(){
197		- std::stringstream ss;
198		- ss<<"\t(x, y, z) = "<<stim::vec3<T>::str();
199		-
200		- if(e[0].size() > 0){
201		- ss<<"\t> ";
202		- for(unsigned int o = 0; o < e[0].size(); o++)
203		- ss<<e[0][o]<<" ";
204		- }
205		- if(e[1].size() > 0){
206		- ss<<"\t< ";
207		- for(unsigned int i = 0; i < e[1].size(); i++)
208		- ss<<e[1][i]<<" ";
209		- }
210		-
211		- return ss.str();
212		- }
213		- ///operator for writing the vector into the stream;
214		- friend std::ofstream& operator<<(std::ofstream& out, const vertex& v)
215		- {
216		- unsigned int s0, s1;
217		- s0 = v.e[0].size();
218		- s1 = v.e[1].size();
219		- out.write(reinterpret_cast<const char>(&v.ptr[0]), 3sizeof(T)); ///write physical vertex location
220		- out.write(reinterpret_cast<const char*>(&s0), sizeof(unsigned int)); ///write the number of "outgoing edges"
221		- out.write(reinterpret_cast<const char*>(&s1), sizeof(unsigned int)); ///write the number of "incoming edges"
222		- if (s0 != 0)
223		- out.write(reinterpret_cast<const char>(&v.e[0][0]), sizeof(unsigned int)v.e[0].size()); ///write the "outgoing edges"
224		- if (s1 != 0)
225		- out.write(reinterpret_cast<const char>(&v.e[1][0]), sizeof(unsigned int)v.e[1].size()); ///write the "incoming edges"
226		- return out;
227		- }
228		-
229		- ///operator for reading the vector out of the stream;
230		- friend std::ifstream& operator>>(std::ifstream& in, vertex& v)
231		- {
232		- in.read(reinterpret_cast<char>(&v[0]), 3sizeof(T)); ///read the physical position
233		- unsigned int s[2];
234		- in.read(reinterpret_cast<char>(&s[0]), 2sizeof(unsigned int)); ///read the sizes of incoming and outgoing edge arrays
235		-
236		- std::vector<unsigned int> one(s[0]);
237		- std::vector<unsigned int> two(s[1]);
238		- v.e[0] = one;
239		- v.e[1] = two;
240		- if (one.size() != 0)
241		- in.read(reinterpret_cast<char>(&v.e[0][0]), s[0] sizeof(unsigned int)); ///read the arrays of "outgoing edges"
242		- if (two.size() != 0)
243		- in.read(reinterpret_cast<char>(&v.e[1][0]), s[1] sizeof(unsigned int)); ///read the arrays of "incoming edges"
244		- return in;
245		- }
246		-
247		- };
248		-
249		-protected:
250		-
251		- std::vector<edge> E; //list of edges
252		- std::vector<vertex> V; //list of vertices.
253		-
254		-public:
255		-
256		- ///default constructor
257		- network()
258		- {
259		-
260		- }
261		-
262		- ///constructor with a file to load.
263		- network(std::string fileLocation)
264		- {
265		- load_obj(fileLocation);
266		- }
267		-
268		- ///Returns the number of edges in the network.
269		- unsigned int edges(){
270		- return E.size();
271		- }
272		-
273		- ///Returns the number of nodes in the network.
274		- unsigned int vertices(){
275		- return V.size();
276		- }
277		-
278		- ///Returns the radius at specific point in the edge
279		- T get_r(unsigned e, unsigned i) {
280		- return E[e].r(i);
281		- }
282		-
283		- ///Returns the average radius of specific edge
284		- T get_average_r(unsigned e) {
285		- T result = 0.0;
286		- unsigned n = E[e].size();
287		- for (unsigned p = 0; p < n; p++)
288		- result += E[e].r(p);
289		-
290		- return (T)result / n;
291		- }
292		-
293		- ///Returns the length of current edge
294		- T get_l(unsigned e) {
295		- return E[e].length();
296		- }
297		-
298		- ///Returns the start vertex of current edge
299		- size_t get_start_vertex(unsigned e) {
300		- return E[e].v[0];
301		- }
302		-
303		- ///Returns the end vertex of current edge
304		- size_t get_end_vertex(unsigned e) {
305		- return E[e].v[1];
306		- }
307		-
308		- ///Returns one vertex
309		- stim::vec3<T> get_vertex(unsigned i) {
310		- return V[i];
311		- }
312		-
313		- ///Returns the boundary vertices' indices
314		- std::vector<unsigned> get_boundary_vertex() {
315		- std::vector<unsigned> result;
316		-
317		- for (unsigned v = 0; v < V.size(); v++) {
318		- if (V[v].e[0].size() + V[v].e[1].size() == 1) { // boundary vertex
319		- result.push_back(v);
320		- }
321		- }
322		-
323		- return result;
324		- }
325		-
326		- ///Set radius
327		- void set_r(unsigned e, std::vector<T> radius) {
328		- E[e].cylinder<T>::copy_r(radius);
329		- }
330		-
331		- void set_r(unsigned e, T radius) {
332		- for (size_t i = 0; i < E[e].size(); i++)
333		- E[e].cylinder<T>::set_r(i, radius);
334		- }
335		- //scale the network by some constant value
336		- // I don't think these work??????
337		- /std::vector<vertex> operator(T s){
338		- for (unsigned i=0; i< vertices; i ++ ){
339		- V[i] = V[i] * s;
340		- }
341		- return V;
342		- }
343		-
344		- std::vector<vertex> operator*(vec<T> s){
345		- for (unsigned i=0; i< vertices; i ++ ){
346		- for (unsigned dim = 0 ; dim< 3; dim ++){
347		- V[i][dim] = V[i][dim] * s[dim];
348		- }
349		- }
350		- return V;
351		- }*/
352		-
353		- // Returns an average of branching index in the network
354		- double BranchingIndex(){
355		- double B=0;
356		- for(unsigned v=0; v < V.size(); v ++){
357		- B += ((V[v].e[0].size()) + (V[v].e[1].size()));
358		- }
359		- B = B / V.size();
360		- return B;
361		-
362		- }
363		-
364		- // Returns number of branch points in thenetwork
365		- unsigned int BranchP(){
366		- unsigned int B=0;
367		- unsigned int c;
368		- for(unsigned v=0; v < V.size(); v ++){
369		- c = ((V[v].e[0].size()) + (V[v].e[1].size()));
370		- if (c > 2){
371		- B += 1;}
372		- }
373		- return B;
374		-
375		- }
376		-
377		- // Returns number of end points (tips) in thenetwork
378		- unsigned int EndP(){
379		- unsigned int B=0;
380		- unsigned int c;
381		- for(unsigned v=0; v < V.size(); v ++){
382		- c = ((V[v].e[0].size()) + (V[v].e[1].size()));
383		- if (c == 1){
384		- B += 1;}
385		- }
386		- return B;
387		-
388		- }
389		-
390		- //// Returns a dictionary with the key as the vertex
391		- //std::map<std::vector<vertex>,unsigned int> DegreeDict(){
392		- // std::map<std::vector<vertex>,unsigned int> dd;
393		- // unsigned int c = 0;
394		- // for(unsigned v=0; v < V.size(); v ++){
395		- // c = ((V[v].e[0].size()) + (V[v].e[1].size()));
396		- // dd[V[v]] = c;
397		- // }
398		- // return dd;
399		- //}
400		-
401		- //// Return number of branching stems
402		- //unsigned int Stems(){
403		- // unsigned int s = 0;
404		- // std::map<std::vector<vertex>,unsigned int> dd;
405		- // dd = DegreeDict();
406		- // //for(unsigned v=0; v < V.size(); v ++){
407		- // // V[v].e[0].
408		- // return s;
409		- //}
410		-
411		- //Calculate Metrics---------------------------------------------------
412		- // Returns an average of fiber/edge lengths in the network
413		- double Lengths(){
414		- stim::vec<T> L;
415		- double sumLength = 0;
416		- for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
417		- L.push_back(E[e].length()); //append the edge length
418		- sumLength = sumLength + E[e].length();
419		- }
420		- double avg = sumLength / E.size();
421		- return avg;
422		- }
423		-
424		-
425		- // Returns an average of tortuosities in the network
426		- double Tortuosities(){
427		- stim::vec<T> t;
428		- stim::vec<T> id1, id2; // starting and ending vertices of the edge
429		- double distance;double tortuosity;double sumTortuosity = 0;
430		- for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
431		- id1 = E[e][0]; //get the edge starting point
432		- id2 = E[e][E[e].size() - 1]; //get the edge ending point
433		- distance = (id1 - id2).len(); //displacement between the starting and ending points
434		- if(distance > 0){
435		- tortuosity = E[e].length()/ distance ; // tortuoisty = edge length / edge displacement
436		- }
437		- else{
438		- tortuosity = 0;}
439		- t.push_back(tortuosity);
440		- sumTortuosity += tortuosity;
441		- }
442		- double avg = sumTortuosity / E.size();
443		- return avg;
444		- }
445		-
446		- // Returns average contraction of the network
447		- double Contractions(){
448		- stim::vec<T> t;
449		- stim::vec<T> id1, id2; // starting and ending vertices of the edge
450		- double distance;double contraction;double sumContraction = 0;
451		- for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
452		- id1 = E[e][0]; //get the edge starting point
453		- id2 = E[e][E[e].size() - 1]; //get the edge ending point
454		- distance = (id1 - id2).len(); //displacement between the starting and ending points
455		- contraction = distance / E[e].length(); // tortuoisty = edge length / edge displacement
456		- t.push_back(contraction);
457		- sumContraction += contraction;
458		- }
459		- double avg = sumContraction / E.size();
460		- return avg;
461		- }
462		-
463		- // returns average fractal dimension of the branches of the network
464		- double FractalDimensions(){
465		- stim::vec<T> t;
466		- stim::vec<T> id1, id2; // starting and ending vertices of the edge
467		- double distance;double fract;double sumFractDim = 0;
468		- for(unsigned e = 0; e < E.size(); e++){ //for each edge in the network
469		- id1 = E[e][0]; //get the edge starting point
470		- id2 = E[e][E[e].size() - 1]; //get the edge ending point
471		- distance = (id1 - id2).len(); //displacement between the starting and ending points
472		- fract = std::log(distance) / std::log(E[e].length()); // tortuoisty = edge length / edge displacement
473		- t.push_back(sumFractDim);
474		- sumFractDim += fract;
475		- }
476		- double avg = sumFractDim / E.size();
477		- return avg;
478		- }
479		-
480		- //returns a cylinder represented a given fiber (based on edge index)
481		- stim::cylinder<T> get_cylinder(unsigned e){
482		- return E[e]; //return the specified edge (casting it to a fiber)
483		- }
484		-
485		- //load a network from an OBJ file
486		- void load_obj(std::string filename){
487		-
488		- stim::obj<T> O; //create an OBJ object
489		- O.load(filename); //load the OBJ file as an object
490		-
491		- std::vector<unsigned> id2vert; //this list stores the OBJ vertex ID associated with each network vertex
492		-
493		- unsigned i[2]; //temporary, IDs associated with the first and last points in an OBJ line
494		-
495		- //for each line in the OBJ object
496		- for(unsigned int l = 1; l <= O.numL(); l++){
497		-
498		- std::vector< stim::vec<T> > c; //allocate an array of points for the vessel centerline
499		- O.getLine(l, c); //get the fiber centerline
500		-
501		- stim::centerline<T> c3(c.size());
502		- for(size_t j = 0; j < c.size(); j++)
503		- c3[j] = c[j];
504		- c3.update();
505		-
506		- // edge new_edge = c3; ///This is dangerous.
507		- edge new_edge(c3);
508		-
509		- //create an edge from the given centerline
510		- unsigned int I = new_edge.size(); //calculate the number of points on the centerline
511		-
512		- //get the first and last vertex IDs for the line
513		- std::vector< unsigned > id; //create an array to store the centerline point IDs
514		- O.getLinei(l, id); //get the list of point IDs for the line
515		- i[0] = id.front(); //get the OBJ ID for the first element of the line
516		- i[1] = id.back(); //get the OBJ ID for the last element of the line
517		-
518		- std::vector<unsigned>::iterator it; //create an iterator for searching the id2vert array
519		- unsigned it_idx; //create an integer for the id2vert entry index
520		-
521		- //find out if the nodes for this fiber have already been created
522		- it = find(id2vert.begin(), id2vert.end(), i[0]); //look for the first node
523		- if(it == id2vert.end()){ //if i[0] hasn't already been used
524		- vertex new_vertex = new_edge[0]; //create a new vertex, assign it a position
525		- bool flag = false;
526		- unsigned j = 0;
527		- for (; j < V.size(); j++) { // check whether current vertex is already exist
528		- if (new_vertex == V[j]) {
529		- flag = true;
530		- break;
531		- }
532		- }
533		- if (!flag) { // unique one
534		- new_vertex.e[0].push_back(E.size()); //add the current edge as outgoing
535		- new_edge.v[0] = V.size(); //add the new edge to the edge
536		- V.push_back(new_vertex); //add the new vertex to the vertex list
537		- id2vert.push_back(i[0]); //add the ID to the ID->vertex conversion list
538		- }
539		- else {
540		- V[j].e[0].push_back(E.size());
541		- new_edge.v[0] = j;
542		- }
543		- }
544		- else{ //if the vertex already exists
545		- it_idx = std::distance(id2vert.begin(), it);
546		- V[it_idx].e[0].push_back(E.size()); //add the current edge as outgoing
547		- new_edge.v[0] = it_idx;
548		- }
549		-
550		- it = find(id2vert.begin(), id2vert.end(), i[1]); //look for the second ID
551		- if(it == id2vert.end()){ //if i[1] hasn't already been used
552		- vertex new_vertex = new_edge[I-1]; //create a new vertex, assign it a position
553		- bool flag = false;
554		- unsigned j = 0;
555		- for (; j < V.size(); j++) { // check whether current vertex is already exist
556		- if (new_vertex == V[j]) {
557		- flag = true;
558		- break;
559		- }
560		- }
561		- if (!flag) {
562		- new_vertex.e[1].push_back(E.size()); //add the current edge as incoming
563		- new_edge.v[1] = V.size(); //add the new vertex to the edge
564		- V.push_back(new_vertex); //add the new vertex to the vertex list
565		- id2vert.push_back(i[1]); //add the ID to the ID->vertex conversion list
566		- }
567		- else {
568		- V[j].e[1].push_back(E.size());
569		- new_edge.v[1] = j;
570		- }
571		- }
572		- else{ //if the vertex already exists
573		- it_idx = std::distance(id2vert.begin(), it);
574		- V[it_idx].e[1].push_back(E.size()); //add the current edge as incoming
575		- new_edge.v[1] = it_idx;
576		- }
577		-
578		- E.push_back(new_edge); //push the edge to the list
579		-
580		- }
581		-
582		- // copy the radii information from OBJ
583		- /*if (O.numVT()) {
584		- unsigned k = 0;
585		- for (unsigned i = 0; i < E.size(); i++) {
586		- for (unsigned j = 0; j < E[i].size(); j++) {
587		- E[i].cylinder<T>::set_r(j, O.getVT(k)[0] / 2);
588		- k++;
589		- }
590		- }
591		- }*/
592		- // OBJ class assumes that in L the two values are equal
593		- if (O.numVT()) {
594		- std::vector< unsigned > id; //create an array to store the centerline point IDs
595		- for (unsigned i = 0; i < O.numL(); i++) {
596		- id.clear();
597		- O.getLinei(i + 1, id); //get the list of point IDs for the line
598		- for (unsigned j = 0; j < id.size(); j++)
599		- E[i].cylinder<T>::set_r(j, O.getVT(id[j] - 1)[0] / 2);
600		- }
601		- }
602		- }
603		-
604		- ///loads a .nwt file. Reads the header and loads the data into the network according to the header.
605		- void
606		- loadNwt(std::string filename)
607		- {
608		- int dims[2]; ///number of vertex, number of edges
609		- readHeader(filename, &dims[0]); //read header
610		- std::ifstream file;
611		- file.open(filename.c_str(), std::ios::in \| std::ios::binary); ///skip header information.
612		- file.seekg(14+58+4+4, file.beg);
613		- vertex v;
614		- for(int i = 0; i < dims[0]; i++) ///for every vertex, read vertex, add to network.
615		- {
616		- file >> v;
617		- V.push_back(v);
618		-// std::cout << i << " " << v.str() << std::endl;
619		- }
620		-
621		- std::cout << std::endl;
622		- for(int i = 0; i < dims[1]; i++) ///for every edge, read edge, add to network.
623		- {
624		- edge e;
625		- file >> e;
626		- E.push_back(e);
627		- //std::cout << i << " " << E[i].str() << std::endl; // not necessary?
628		- }
629		- file.close();
630		- }
631		-
632		- ///saves a .nwt file. Writes the header in raw text format, then saves the network as a binary file.
633		- void
634		- saveNwt(std::string filename)
635		- {
636		- writeHeader(filename);
637		- std::ofstream file;
638		- file.open(filename.c_str(), std::ios::out \| std::ios::binary \| std::ios::app); ///since we have written the header we are not appending.
639		- for(int i = 0; i < V.size(); i++) ///look through the Vertices and write each one.
640		- {
641		-// std::cout << i << " " << V[i].str() << std::endl;
642		- file << V[i];
643		- }
644		- for(int i = 0; i < E.size(); i++) ///loop through the Edges and write each one.
645		- {
646		- //std::cout << i << " " << E[i].str() << std::endl; // not necesarry?
647		- file << E[i];
648		- }
649		- file.close();
650		- }
651		-
652		-
653		- ///Writes the header information to a .nwt file.
654		- void
655		- writeHeader(std::string filename)
656		- {
657		- std::string magicString = "nwtFileFormat "; ///identifier for the file.
658		- std::string desc = "fileid(14B), desc(58B), #vertices(4B), #edges(4B): bindata";
659		- int hNumVertices = V.size(); ///int byte header storing the number of vertices in the file
660		- int hNumEdges = E.size(); ///int byte header storing the number of edges.
661		- std::ofstream file;
662		- file.open(filename.c_str(), std::ios::out \| std::ios::binary);
663		- std::cout << hNumVertices << " " << hNumEdges << std::endl;
664		- file.write(reinterpret_cast<const char*>(&magicString.c_str()[0]), 14); //write the file id
665		- file.write(reinterpret_cast<const char*>(&desc.c_str()[0]), 58); //write the description
666		- file.write(reinterpret_cast<const char*>(&hNumVertices), sizeof(int)); //write #vert.
667		- file.write(reinterpret_cast<const char*>(&hNumEdges), sizeof(int)); //write #edges
668		-// file << magicString.c_str() << desc.c_str() << hNumVertices << hNumEdges;
669		- file.close();
670		-
671		- }
672		-
673		- ///Reads the header information from a .nwt file.
674		- void
675		- readHeader(std::string filename, int *dims)
676		- {
677		- char magicString[14]; ///id
678		- char desc[58]; ///description
679		- int hNumVertices; ///#vert
680		- int hNumEdges; ///#edges
681		- std::ifstream file; ////create stream
682		- file.open(filename.c_str(), std::ios::in \| std::ios::binary);
683		- file.read(reinterpret_cast<char*>(&magicString[0]), 14); ///read the file id.
684		- file.read(reinterpret_cast<char*>(&desc[0]), 58); ///read the description
685		- file.read(reinterpret_cast<char*>(&hNumVertices), sizeof(int)); ///read the number of vertices
686		- file.read(reinterpret_cast<char*>(&hNumEdges), sizeof(int)); ///read the number of edges
687		-// std::cout << magicString << desc << hNumVertices << " " << hNumEdges << std::endl;
688		- file.close(); ///close the file.
689		- dims[0] = hNumVertices; ///fill the returned reference.
690		- dims[1] = hNumEdges;
691		- }
692		-
693		- //load a network from an SWC file
694		- void load_swc(std::string filename) {
695		- stim::swc<T> S; // create swc variable
696		- S.load(filename); // load the node information
697		- S.create_tree(); // link those node according to their linking relationships as a tree
698		- S.resample();
699		-
700		- //NT.push_back(S.node[0].type); // set the neuronal_type value to the first vertex in the network
701		- std::vector<unsigned> id2vert; // this list stores the SWC vertex ID associated with each network vertex
702		- unsigned i[2]; // temporary, IDs associated with the first and last points
703		-
704		- for (unsigned int l = 0; l < S.numE(); l++) { // for every edge
705		- //NT.push_back(S.node[l].type);
706		-
707		- std::vector< stim::vec3<T> > c;
708		- S.get_points(l, c);
709		-
710		- stim::centerline<T> c3(c.size()); // new fiber
711		-
712		- for (unsigned j = 0; j < c.size(); j++)
713		- c3[j] = c[j]; // copy the points
714		-
715		- c3.update(); // update the L information
716		-
717		- stim::cylinder<T> C3(c3); // create a new cylinder in order to copy the origin radius information
718		- // upadate the R information
719		- std::vector<T> radius;
720		- S.get_radius(l, radius);
721		-
722		- C3.copy_r(radius);
723		-
724		- edge new_edge(C3); // new edge
725		-
726		- //create an edge from the given centerline
727		- unsigned int I = new_edge.size(); //calculate the number of points on the centerline
728