fixed bugs in merging new network format, temporary version

Jiaming Guo
1 parent 4ea8daa8
Showing 3 changed files with 107 additions and 48 deletions Show diff stats
stim/biomodels/centerline.h
stim/biomodels/network.h
stim/visualization/cylinder.h
@@ -31,20 +31,6 @@ protected:
 		vec3<T> ab = a.norm() + b.norm();
 		return ab.norm();
 	}
-	//initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
-	void update_L(size_t start = 0) {
-		L.resize(size());									//allocate space for the L array
-		if (start == 0) {
-			L[0] = 0;											//initialize the length value for the first point to zero (0)
-			start++;
-		}
-
-		stim::vec3<T> d;
-		for (size_t i = start; i < size(); i++) {		//for each line segment in the centerline
-			d = at(i) - at(i - 1);
-			L[i] = L[i - 1] + d.len();				//calculate the running length total
-		}
-	}
 	void init() {
 		if (size() == 0) return;								//return if there aren't any points
@@ -135,6 +121,21 @@ public:
 		return L.back();
 	}
+	//initializes the integrated length vector to make parameterization easier, starting with index idx (all previous indices are assumed to be correct)
+	void update_L(size_t start = 0) {
+		L.resize(size());									//allocate space for the L array
+		if (start == 0) {
+			L[0] = 0;											//initialize the length value for the first point to zero (0)
+			start++;
+		}
+
+		stim::vec3<T> d;
+		for (size_t i = start; i < size(); i++) {		//for each line segment in the centerline
+			d = at(i) - at(i - 1);
+			L[i] = L[i - 1] + d.len();				//calculate the running length total
+		}
+	}
+
 	/// stitch two centerlines
 	///@param c1, c2: two centerlines
 	///@param sigma: sample rate
@@ -63,6 +63,7 @@ class network{
 	class edge : public cylinder<T>
 	{
 		public:
+	
 		unsigned int v[2];		//unique id's designating the starting and ending
 		// default constructor
 		edge() : cylinder<T>() {
@@ -127,7 +128,7 @@ class network{
 			//	for(int j = 0; j < nmags(); j++)	//future code for multiple mags
 			//	{
 				out.write(reinterpret_cast<const char*>(&e.R[i]), sizeof(T));	///write the radius
-				std::cout << point.str() << " " << e.R[i] << std::endl;
+				//std::cout << point.str() << " " << e.R[i] << std::endl;
 			//	}
 			}
 			return out;	//return stream
@@ -141,7 +142,7 @@ class network{
 			in.read(reinterpret_cast<char*>(&v1), sizeof(unsigned int));	//read the ending point
 			in.read(reinterpret_cast<char*>(&sz), sizeof(unsigned int));	//read the number of points in the edge
 			stim::centerline<T> temp = stim::centerline<T>(sz);		//allocate the new edge
-			e = edge(temp);		
+			e = edge(temp);
 			e.v[0] = v0; e.v[1] = v1;
 			for(int i = 0; i < sz; i++)		//set the points and radii to the newly read values
 			{
@@ -152,10 +153,12 @@ class network{
 		//				for(int j = 0; j < nmags(); j++)		///future code for mags
 		//				{
 				in.read(reinterpret_cast<char*>(&mag), sizeof(T));	
-				e.set_r(0, mag);
-				std::cout << point.str() << " " << mag << std::endl;
+				e.set_r(i, mag);
+				//std::cout << point.str() << " " << mag << std::endl;
 		//				}
 			}
+			e.init();
+			e.update();
 			return in;
 		}
 	};
@@ -193,21 +196,23 @@ class network{
 				return ss.str();
 			}
-			///operator for writing the edge into the stream;
+			///operator for writing the vector into the stream;
 			friend std::ofstream& operator<<(std::ofstream& out, const vertex& v)
 			{
 				unsigned int s0, s1;
 				s0 = v.e[0].size();
 				s1 = v.e[1].size();
 				out.write(reinterpret_cast<const char*>(&v.ptr[0]), 3*sizeof(T));	///write physical vertex location
-				out.write(reinterpret_cast<const char*>(&s0), sizeof(unsigned int));	///write the number of "incoming edges"
-				out.write(reinterpret_cast<const char*>(&s1), sizeof(unsigned int));	///write the number of "outgoing edges"	
-				out.write(reinterpret_cast<const char*>(&v.e[0][0]), sizeof(unsigned int)*v.e[0].size());	///write the "incoming edges"
-				out.write(reinterpret_cast<const char*>(&v.e[1][0]), sizeof(unsigned int)*v.e[1].size());	///write the "outgoing edges"
+				out.write(reinterpret_cast<const char*>(&s0), sizeof(unsigned int));	///write the number of "outgoing edges"
+				out.write(reinterpret_cast<const char*>(&s1), sizeof(unsigned int));	///write the number of "incoming edges"	
+				if (s0 != 0)
+					out.write(reinterpret_cast<const char*>(&v.e[0][0]), sizeof(unsigned int)*v.e[0].size());	///write the "outgoing edges"
+				if (s1 != 0)
+					out.write(reinterpret_cast<const char*>(&v.e[1][0]), sizeof(unsigned int)*v.e[1].size());	///write the "incoming edges"
 				return out;
 			}
-			///operator for reading the edge out of the stream;
+			///operator for reading the vector out of the stream;
 			friend std::ifstream& operator>>(std::ifstream& in, vertex& v)
 			{
 				in.read(reinterpret_cast<char*>(&v[0]), 3*sizeof(T));	///read the physical position
@@ -218,8 +223,10 @@ class network{
 				std::vector<unsigned int> two(s[1]);
 				v.e[0] = one;
 				v.e[1] = two;
-				in.read(reinterpret_cast<char*>(&v.e[0][0]),s[0]*sizeof(unsigned int));		///read the arrays of "incoming edges"
-				in.read(reinterpret_cast<char*>(&v.e[1][0]),s[1]*sizeof(unsigned int));		///read the arrays of "outgoing edges"
+				if (one.size() != 0)
+					in.read(reinterpret_cast<char*>(&v.e[0][0]), s[0] * sizeof(unsigned int));		///read the arrays of "outgoing edges"
+				if (two.size() != 0)
+					in.read(reinterpret_cast<char*>(&v.e[1][0]), s[1] * sizeof(unsigned int));		///read the arrays of "incoming edges"
 				return in;
 			}
@@ -501,10 +508,24 @@ public:
 			it = find(id2vert.begin(), id2vert.end(), i[0]);	//look for the first node
 			if(it == id2vert.end()){							//if i[0] hasn't already been used
 				vertex new_vertex = new_edge[0];				//create a new vertex, assign it a position
-				new_vertex.e[0].push_back(E.size());				//add the current edge as outgoing
-				new_edge.v[0] = V.size();					//add the new edge to the edge
-				V.push_back(new_vertex);					//add the new vertex to the vertex list
-				id2vert.push_back(i[0]);					//add the ID to the ID->vertex conversion list
+				bool flag = false;
+				unsigned j = 0;
+				for (; j < V.size(); j++) {						// check whether current vertex is already exist
+					if (new_vertex == V[j]) {
+						flag = true;
+						break;
+					}
+				}
+				if (!flag) {									// unique one
+					new_vertex.e[0].push_back(E.size());				//add the current edge as outgoing
+					new_edge.v[0] = V.size();					//add the new edge to the edge
+					V.push_back(new_vertex);					//add the new vertex to the vertex list
+					id2vert.push_back(i[0]);					//add the ID to the ID->vertex conversion list
+				}
+				else {
+					V[j].e[0].push_back(E.size());
+					new_edge.v[0] = j;
+				}
 			}
 			else{									//if the vertex already exists
 				it_idx = std::distance(id2vert.begin(), it);
@@ -515,10 +536,24 @@ public:
 			it = find(id2vert.begin(), id2vert.end(), i[1]);			//look for the second ID
 			if(it == id2vert.end()){						//if i[1] hasn't already been used
 				vertex new_vertex = new_edge[I-1];				//create a new vertex, assign it a position
-				new_vertex.e[1].push_back(E.size());				//add the current edge as incoming
-				new_edge.v[1] = V.size();                                  	//add the new vertex to the edge
-				V.push_back(new_vertex);					//add the new vertex to the vertex list
-				id2vert.push_back(i[1]);					//add the ID to the ID->vertex conversion list
+				bool flag = false;
+				unsigned j = 0;
+				for (; j < V.size(); j++) {					// check whether current vertex is already exist
+					if (new_vertex == V[j]) {
+						flag = true;
+						break;
+					}
+				}
+				if (!flag) {
+					new_vertex.e[1].push_back(E.size());				//add the current edge as incoming
+					new_edge.v[1] = V.size();                                  	//add the new vertex to the edge
+					V.push_back(new_vertex);					//add the new vertex to the vertex list
+					id2vert.push_back(i[1]);					//add the ID to the ID->vertex conversion list
+				}
+				else {
+					V[j].e[1].push_back(E.size());
+					new_edge.v[1] = j;
+				}
 			}
 			else{									//if the vertex already exists
 				it_idx = std::distance(id2vert.begin(), it);
@@ -529,6 +564,27 @@ public:
 			E.push_back(new_edge);							//push the edge to the list
 		}
+
+		// copy the radii information from OBJ
+		/*if (O.numVT()) {
+			unsigned k = 0;
+			for (unsigned i = 0; i < E.size(); i++) {
+				for (unsigned j = 0; j < E[i].size(); j++) {
+					E[i].cylinder<T>::set_r(j, O.getVT(k)[0] / 2);
+					k++;
+				}
+			}
+		}*/
+		// OBJ class assumes that in L the two values are equal
+		if (O.numVT()) {
+			std::vector< unsigned > id;						//create an array to store the centerline point IDs
+			for (unsigned i = 0; i < O.numL(); i++) {
+				id.clear();
+				O.getLinei(i + 1, id);							//get the list of point IDs for the line
+				for (unsigned j = 0; j < id.size(); j++)
+					E[i].cylinder<T>::set_r(j, O.getVT(id[j] - 1)[0] / 2);
+			}
+		}
 	}
 	///loads a .nwt file. Reads the header and loads the data into the network according to the header.
@@ -554,7 +610,7 @@ public:
 			edge e;
 			file >> e;
 			E.push_back(e);
-			std::cout << i << " " << E[i].str() << std::endl;
+			//std::cout << i << " " << E[i].str() << std::endl;		// not necessary?
 		}
 		file.close();
 	}
@@ -573,7 +629,7 @@ public:
 		}
 		for(int i = 0; i < E.size(); i++)	///loop through the Edges and write each one.
 		{
-			std::cout << i << " " << E[i].str() << std::endl;
+			//std::cout << i << " " << E[i].str() << std::endl;		// not necesarry?
 			file << E[i];
 		}
 		file.close();
@@ -21,18 +21,6 @@ protected:
 	//calculates the U values for each point to initialize the frenet frame
 	//	this function assumes that the centerline has already been set
-	void init() {
-		U.resize(size());								//allocate space for the frenet frame vectors
-		R.resize(size());
-	
-		stim::circle<T> c;								//create a circle
-		stim::vec3<T> d0, d1;
-		for (size_t i = 0; i < size() - 1; i++) {		//for each line segment in the centerline
-			c.rotate(d(i));								//rotate the circle to match that normal
-			U[i] = c.U;									//save the U vector from the circle
-		}
-		U[size() - 1] = c.U;							//for the last point, duplicate the final frenet frame vector
-	}
 public:
@@ -94,6 +82,20 @@ public:
 		return R[i];
 	}
+	void init() {
+		U.resize(size());			//allocate space for the frenet frame vectors
+//		if (R.size() != 0)
+		R.resize(size());
+
+		stim::circle<T> c;								//create a circle
+		stim::vec3<T> d0, d1;
+		for (size_t i = 0; i < size() - 1; i++) {		//for each line segment in the centerline
+			c.rotate(d(i));								//rotate the circle to match that normal
+			U[i] = c.U;									//save the U vector from the circle
+		}
+		U[size() - 1] = c.U;							//for the last point, duplicate the final frenet frame vector
+	}
+
 	///adds a magnitude to each point in the cylinder
 	/*void add_mag(V val = 0) {
 		if (M.size() == 0) M.resize(size());	//if the magnitude vector isn't initialized, resize it to match the centerline