codebase / stimlib

Browse Code »

Commit a9ed210f0018dafe10c4bfe1e07d62535b92a489

Authored by Pavel Govyadinov 2016-08-26 13:35:45 -0500

1 parent 0311ab81

added centerline

Showing 1 changed file with 346 additions and 0 deletions Show diff stats

Inline Side-by-side

stim/biomodels/centerline.h 0 → 100644

Wrap text Show/Hide comments View file @a9ed210

	1	+#ifndef STIM_CENTERLINE_H
	2	+#define STIM_CENTERLINE_H
	3	+
	4	+#include <vector>
	5	+#include <stim/math/vec3.h>
	6	+#include <ANN/ANN.h>
	7	+
	8	+namespace stim{
	9	+
	10	+/** This class stores information about a single fiber represented as a set of geometric points
	11	+ * between two branch or end points. This class is used as a fundamental component of the stim::network
	12	+ * class to describe an interconnected (often biological) network.
	13	+ */
	14	+template<typename T>
	15	+class centerline{
	16	+
	17	+protected:
	18	+ unsigned int N; //number of points in the fiber
	19	+ double **c; //centerline (array of double pointers)
	20	+ ANNkd_tree* kdt; //kd-tree stores all points in the fiber for fast searching
	21	+
	22	+ /// Initialize an empty fiber
	23	+ void init()
	24	+ {
	25	+ N=0;
	26	+ c=NULL;
	27	+ kdt = NULL;
	28	+ }
	29	+
	30	+ /// Initialize a fiber with N centerline points (all located at [0, 0, 0] with radius 0)
	31	+ void init(unsigned int n)
	32	+ {
	33	+
	34	+ N = n; //set the number of points
	35	+ kdt = NULL;
	36	+ c = (double*) malloc(sizeof(double) * N); //allocate the array pointer
	37	+
	38	+ for(unsigned int i = 0; i < N; i++) //allocate space for each point
	39	+ c[i] = (double) malloc(sizeof(double) 3);
	40	+ }
	41	+
	42	+ /// Copies an existing fiber to the current fiber
	43	+
	44	+ /// @param cpy stores the new copy of the fiber
	45	+ void copy( const stim::centerline<T>& cpy, bool kd = 0){
	46	+
	47	+ ///allocate space for the new fiber
	48	+ init(cpy.N);
	49	+
	50	+ ///copy the points
	51	+ for(unsigned int i = 0; i < N; i++){
	52	+ for(unsigned int d = 0; d < 3; d++) //for each dimension
	53	+ c[i][d] = cpy.c[i][d]; //copy the coordinate
	54	+ }
	55	+ if(kd)
	56	+ gen_kdtree(); //generate the kd tree for the new fiber
	57	+ }
	58	+
	59	+ /// generate a KD tree for points on fiber
	60	+ void gen_kdtree()
	61	+ {
	62	+ int n_data = N; //create an array of data points
	63	+ ANNpointArray pts = (ANNpointArray)c; //cast the centerline list to an ANNpointArray
	64	+ kdt = new ANNkd_tree(pts, n_data, 3); //build a KD tree
	65	+ }
	66	+
	67	+ /// find distance between two points
	68	+ double dist(double* p0, double* p1){
	69	+
	70	+ double sum = 0; // initialize variables
	71	+ float v;
	72	+ for(unsigned int d = 0; d < 3; d++)
	73	+ {
	74	+ v = p1[d] - p0[d];
	75	+ sum +=v * v;
	76	+ }
	77	+ return sqrt(sum);
	78	+ }
	79	+
	80	+ /// This function retreives the index for the fiber point closest to q
	81	+
	82	+ /// @param q is a reference point used to find the closest point on the fiber center line
	83	+ unsigned int ann( stim::vec<double> q ){
	84	+ ANNidxArray idx = new ANNidx[1]; //variable used to hold the nearest point
	85	+ ANNdistArray sq_dist = new ANNdist[1]; //variable used to hold the squared distance to the nearest point
	86	+ kdt->annkSearch(q.data(), 1, idx, sq_dist); //search the KD tree for the nearest neighbor
	87	+ return *idx;
	88	+ }
	89	+
	90	+ /// Returns a stim::vec representing the point at index i
	91	+
	92	+ /// @param i is an index of the desired centerline point
	93	+ stim::vec<T> get_vec(unsigned i){
	94	+ stim::vec3<T> r;
	95	+ r.resize(3);
	96	+ r[0] = c[i][0];
	97	+ r[1] = c[i][1];
	98	+ r[2] = c[i][2];
	99	+
	100	+ return r;
	101	+ }
	102	+
	103	+
	104	+public:
	105	+
	106	+ centerline(){
	107	+ init();
	108	+ }
	109	+
	110	+ /// Copy constructor
	111	+ centerline(const stim::centerline<T> &obj){
	112	+ copy(obj);
	113	+ }
	114	+
	115	+ //temp constructor for graph visualization
	116	+ centerline(int n)
	117	+ {
	118	+ init(n);
	119	+ }
	120	+
	121	+ /// Constructor takes a list of stim::vec points, the radius at each point is set to zero
	122	+ centerline(std::vector< stim::vec<T> > p, bool kd = 0){
	123	+ init(p.size()); //initialize the fiber
	124	+
	125	+ //for each point, set the centerline position and radius
	126	+ for(unsigned int i = 0; i < N; i++){
	127	+
	128	+ //set the centerline position
	129	+ for(unsigned int d = 0; d < 3; d++)
	130	+ c[i][d] = (double) p[i][d];
	131	+
	132	+ //set the radius
	133	+ }
	134	+ //generate a kd tree
	135	+ if(kd)
	136	+ gen_kdtree();
	137	+ }
	138	+
	139	+ /// constructor takes a list of points
	140	+ centerline(std::vector< stim::vec3< T > > pos, bool kd = 0){
	141	+ init(pos.size()); //initialize the fiber
	142	+
	143	+ //for each point, set the centerline position and radius
	144	+ for(unsigned int i = 0; i < N; i++){
	145	+ //set the centerline position
	146	+ for(unsigned int d = 0; d < 3; d++)
	147	+ c[i][d] = (double) pos[i][d];
	148	+ //set the radius
	149	+ }
	150	+
	151	+ //generate a kd tree
	152	+ if(kd)
	153	+ gen_kdtree();
	154	+ }
	155	+
	156	+ /// Assignment operation
	157	+ centerline& operator=(const centerline &rhs){
	158	+ if(this == &rhs) return *this; //test for and handle self-assignment
	159	+ copy(rhs);
	160	+ return *this;
	161	+ }
	162	+
	163	+
	164	+ /// Return the point on the fiber closest to q
	165	+ /// @param q is the query point used to locate the nearest point on the fiber centerline
	166	+ stim::vec<T> nearest(stim::vec<T> q){
	167	+
	168	+ stim::vec<double> temp( (double) q[0], (double) q[1], (double) q[2]);
	169	+
	170	+ unsigned int idx = ann(temp); //determine the index of the nearest neighbor
	171	+
	172	+ return stim::vec<T>((T) c[idx][0], (T) c[idx][1], (T) c[idx][2]); //return the nearest centerline point
	173	+ }
	174	+
	175	+ /// Return the point index on the fiber closest to q
	176	+ /// @param q is the query point used to locate the nearest point on the fiber centerline
	177	+ unsigned int nearest_idx(stim::vec<T> q){
	178	+
	179	+ stim::vec<double> temp((double) q[0], (double) q[1], (double) q[2]);
	180	+
	181	+ unsigned int idx = ann(temp); //determine the index of the nearest neighbor
	182	+
	183	+ return idx; //return the nearest centerline point index
	184	+ }
	185	+
	186	+ /// Returns the fiber centerline as an array of stim::vec points
	187	+ std::vector< stim::vec<T> > get_centerline(){
	188	+
	189	+ //create an array of stim vectors
	190	+ std::vector< stim::vec3<T> > pts(N);
	191	+
	192	+ //cast each point to a stim::vec, keeping only the position information
	193	+ for(unsigned int i = 0; i < N; i++)
	194	+ pts[i] = stim::vec3<T>((T) c[i][0], (T) c[i][1], (T) c[i][2]);
	195	+
	196	+ //return the centerline array
	197	+ return pts;
	198	+ }
	199	+
	200	+ /// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
	201	+ std::vector< stim::centerline<T> > split(unsigned int idx){
	202	+
	203	+ std::vector< stim::centerline<T> > fl; //create an array to store up to two fibers
	204	+
	205	+ //if the index is an end point, only the existing fiber is returned
	206	+ if(idx == 0 \|\| idx == N-1){
	207	+ fl.resize(1); //set the size of the fiber to 1
	208	+ fl[0] = *this; //copy the current fiber
	209	+ }
	210	+
	211	+ //if the index is not an end point
	212	+ else{
	213	+
	214	+ unsigned int N1 = idx + 1; //calculate the size of both fibers
	215	+ unsigned int N2 = N - idx;
	216	+
	217	+ fl.resize(2); //set the array size to 2
	218	+
	219	+ fl[0].init(N1); //set the size of each fiber
	220	+ fl[1].init(N2);
	221	+
	222	+ //copy both halves of the fiber
	223	+ unsigned int i, d;
	224	+
	225	+ //first half
	226	+ for(i = 0; i < N1; i++){ //for each centerline point
	227	+ for(d = 0; d < 3; d++)
	228	+ fl[0].c[i][d] = c[i][d]; //copy each coordinate
	229	+ }
	230	+
	231	+ //second half
	232	+ for(i = 0; i < N2; i++){
	233	+ for(d = 0; d < 3; d++)
	234	+ fl[1].c[i][d] = c[idx + i][d];
	235	+
	236	+ }
	237	+ }
	238	+
	239	+ return fl; //return the array
	240	+
	241	+ }
	242	+
	243	+ /// Calculates the set of fibers resulting from a connection between the current fiber and a fiber f
	244	+
	245	+ /// @param f is the fiber that will be connected to the current fiber
	246	+/* std::vector< stim::centerline<T> > connect( stim::centerline<T> &f, double dist){
	247	+
	248	+ double min_dist;
	249	+ unsigned int idx0, idx1;
	250	+
	251	+ //go through each point in the query fiber, looking for the indices for the closest points
	252	+ for(unsigned int i = 0; i < f.n_pts(); i++){
	253	+ //Run through all points and find the index with the closest point, then partition the fiber and return two fibers.
	254	+
	255	+ }
	256	+
	257	+
	258	+
	259	+ }
	260	+*/
	261	+ /// Outputs the fiber as a string
	262	+ std::string str(){
	263	+ std::stringstream ss;
	264	+
	265	+ //create an iterator for the point list
	266	+ //typename std::list< point<T> >::iterator i;
	267	+ for(unsigned int i = 0; i < N; i++){
	268	+ ss<<" [ ";
	269	+ for(unsigned int d = 0; d < 3; d++){
	270	+ ss<<c[i][d]<<" ";
	271	+ }
	272	+ }
	273	+
	274	+ return ss.str();
	275	+ }
	276	+ /// Returns the number of centerline points in the fiber
	277	+ unsigned int size(){
	278	+ return N;
	279	+ }
	280	+
	281	+
	282	+ /// Bracket operator returns the element at index i
	283	+
	284	+ /// @param i is the index of the element to be returned as a stim::vec
	285	+ stim::vec<T> operator[](unsigned i){
	286	+ return get_vec(i);
	287	+ }
	288	+
	289	+ /// Back method returns the last point in the fiber
	290	+ stim::vec<T> back(){
	291	+ return get_vec(N-1);
	292	+ }
	293	+ ////resample a fiber in the network
	294	+ stim::centerline<T> resample(T spacing)
	295	+ {
	296	+ std::cout<<"fiber::resample()"<<std::endl;
	297	+
	298	+ std::vector<T> v(3); //v-direction vector of the segment
	299	+ stim::vec<T> p(3); //- intermediate point to be added
	300	+ stim::vec<T> p1(3); // p1 - starting point of an segment on the fiber,
	301	+ stim::vec<T> p2(3); // p2 - ending point,
	302	+ double sum=0; //distance summation
	303	+ std::vector<stim::vec<T> > fiberPositions = centerline();
	304	+ std::vector<stim::vec<T> > newPointList; // initialize list of new resampled points on the fiber
	305	+ // for each point on the centerline (skip if it is the last point on centerline)
	306	+ //unsigned int N = fiberPositions.size(); // number of points on the fiber
	307	+ for(unsigned int f=0; f< N-1; f++)
	308	+ {
	309	+
	310	+ p1 = fiberPositions[f]; p2 = fiberPositions[f + 1]; v = p2 - p1;
	311	+ for(unsigned int d = 0; d < 3; d++){
	312	+ sum +=v[d] * v[d];} //length of segment-distance between starting and ending point
	313	+
	314	+ T lengthSegment = sqrt(sum); //find Length of the segment as distance between the starting and ending points of the segment
	315	+
	316	+ if(lengthSegment >= spacing) // if length of the segment is greater than standard deviation resample
	317	+ {
	318	+ // repeat resampling until accumulated stepsize is equsl to length of the segment
	319	+ for(T step=0.0; step<lengthSegment; step+=spacing)
	320	+ {
	321	+ // calculate the resampled point by travelling step size in the direction of normalized gradient vector
	322	+ for(unsigned int i=0; i<3;i++)
	323	+ {
	324	+ p[i] = p1[i] + v[i]*(step/lengthSegment);
	325	+ } //for each dimension
	326	+ // add this resampled points to the new fiber list
	327	+ newPointList.push_back(p);
	328	+ }
	329	+ }
	330	+ 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
	331	+ newPointList.push_back(fiberPositions[f+1]);
	332	+ }
	333	+ newPointList.push_back(fiberPositions[N-1]); //add the last point on the fiber to the new fiber list
	334	+ centerline newFiber(newPointList);
	335	+ return newFiber;
	336	+ }
	337	+
	338	+};
	339	+
	340	+
	341	+
	342	+} //end namespace stim
	343	+
	344	+
	345	+
	346	+#endif
...	...