Merge branch 'master' of git.stim.ee.uh.edu:codebase/stimlib

David Mayerich
2 parents fa8763a6 72175dcf
Showing 3 changed files with 363 additions and 6 deletions Show diff stats
stim/biomodels/centerline.h
stim/envi/bip.h
stim/envi/hsi.h
@@ -3,6 +3,7 @@
  
 #include <vector>
 #include <stim/math/vec3.h>
+#include <stim/structures/kdtree.cuh>
  
 namespace stim{
  
@@ -134,6 +135,362 @@ public:
 		return L.back();
 	}
  
+	/// stitch two centerlines
+	///@param c1, c2: two centerlines
+	///@param sigma: sample rate
+	static std::vector< stim::centerline<T> > stitch(stim::centerline<T> c1, stim::centerline<T> c2 = stim::centerline<T>()) {
+		
+		std::vector< stim::centerline<T> > result;
+		stim::centerline<T> new_centerline;
+		stim::vec3<T> new_vertex;
+
+		// if only one centerline, stitch itself!
+		if (c2.size() == 0) {
+			size_t num = c1.size();
+			size_t id = 100000;							// store the downsample start position
+			T threshold;
+			if (num < 4) {								// if the number of vertex is less than 4, do nothing
+				result.push_back(c1);
+				return result;
+			}
+			else {
+				// test geometry start vertex
+				stim::vec3<T> v1 = c1[1] - c1[0];		// vector from c1[0] to c1[1]
+				for (size_t p = 2; p < num; p++) {		// 90° standard???
+					stim::vec3<T> v2 = c1[p] - c1[0];
+					float cosine = v2.dot(v1);
+					if (cosine < 0) {
+						id = p;
+						threshold = v2.len();
+						break;
+					}
+				}
+				if (id != 100000) {						// find a downsample position on the centerline
+					T* c;
+					c = (T*)malloc(sizeof(T) * (num - id) * 3);
+					for (size_t p = id; p < num; p++) {
+						for (size_t d = 0; d < 3; d++) {
+							c[p * 3 + d] = c1[p][d];
+						}
+					}
+					stim::kdtree<T, 3> kdt;
+					kdt.create(c, num - id, 5);			// create tree
+
+					T* query = (T*)malloc(sizeof(T) * 3);
+					for (size_t d = 0; d < 3; d++)
+						query[d] = c1[0][d];
+					size_t index;
+					T dist;
+
+					kdt.search(query, 1, &index, &dist);
+
+					free(query);
+					free(c);
+
+					if (dist > threshold) {
+						result.push_back(c1);
+					}
+					else {
+						// the loop part
+						new_vertex = c1[index];
+						new_centerline.push_back(new_vertex);
+						for (size_t p = 0; p < index + 1; p++) {
+							new_vertex = c1[p];
+							new_centerline.push_back(new_vertex);
+						}
+						result.push_back(new_centerline);
+						new_centerline.clear();
+
+						// the tail part
+						for (size_t p = index; p < num; p++) {
+							new_vertex = c1[p];
+							new_centerline.push_back(new_vertex);
+						}
+						result.push_back(new_centerline);
+					}
+				}
+				else {	// there is one potential problem that two positions have to be stitched
+						// test geometry end vertex
+					stim::vec3<T> v1 = c1[num - 2] - c1[num - 1];
+					for (size_t p = num - 2; p > 0; p--) {		// 90° standard
+						stim::vec3<T> v2 = c1[p - 1] - c1[num - 1];
+						float cosine = v2.dot(v1);
+						if (cosine < 0) {
+							id = p;
+							threshold = v2.len();
+							break;
+						}
+					}
+					if (id != 100000) {						// find a downsample position
+						T* c;
+						c = (T*)malloc(sizeof(T) * (id + 1) * 3);
+						for (size_t p = 0; p < id + 1; p++) {
+							for (size_t d = 0; d < 3; d++) {
+								c[p * 3 + d] = c1[p][d];
+							}
+						}
+						stim::kdtree<T, 3> kdt;
+						kdt.create(c, id + 1, 5);				// create tree
+
+						T* query = (T*)malloc(sizeof(T) * 1 * 3);
+						for (size_t d = 0; d < 3; d++)
+							query[d] = c1[num - 1][d];
+						size_t index;
+						T dist;
+
+						kdt.search(query, 1, &index, &dist);
+
+						free(query);
+						free(c);
+
+						if (dist > threshold) {
+							result.push_back(c1);
+						}
+						else {
+							// the tail part
+							for (size_t p = 0; p < index + 1; p++) {
+								new_vertex = c1[p];
+								new_centerline.push_back(new_vertex);
+							}
+							result.push_back(new_centerline);
+							new_centerline.clear();
+
+							// the loop part
+							for (size_t p = index; p < num; p++) {
+								new_vertex = c1[p];
+								new_centerline.push_back(new_vertex);
+							}
+							new_vertex = c1[index];
+							new_centerline.push_back(new_vertex);
+							result.push_back(new_centerline);
+						}
+					}
+					else {	// no stitch position
+						result.push_back(c1);
+					}
+				}
+			}
+		}
+
+
+		// two centerlines
+		else {
+			// find stitch position based on nearest neighbors												
+			size_t num1 = c1.size();
+			T* c = (T*)malloc(sizeof(T) * num1 * 3);		// c1 as reference point
+			for (size_t p = 0; p < num1; p++)				// centerline to array
+				for (size_t d = 0; d < 3; d++)				// because right now my kdtree code is a relatively close code, it has its own structure
+					c[p * 3 + d] = c1[p][d];				// I will merge it into stimlib totally in the near future
+
+			stim::kdtree<T, 3> kdt;							// kdtree object
+			kdt.create(c, num1, 5);							// create tree
+
+			size_t num2 = c2.size();						
+			T* query = (T*)malloc(sizeof(T) * num2 * 3);	// c2 as query point
+			for (size_t p = 0; p < num2; p++) {
+				for (size_t d = 0; d < 3; d++) {
+					query[p * 3 + d] = c2[p][d];
+				}
+			}
+			std::vector<size_t> index(num2);
+			std::vector<T> dist(num2);
+
+			kdt.search(query, num2, &index[0], &dist[0]);	// find the nearest neighbors in c1 for c2
+
+			// clear up
+			free(query);
+			free(c);
+
+			// find the average vertex distance of one centerline
+			T sigma1 = 0;
+			T sigma2 = 0;
+			for (size_t p = 0; p < num1 - 1; p++) 
+				sigma1 += (c1[p] - c1[p + 1]).len();
+			for (size_t p = 0; p < num2 - 1; p++)
+				sigma2 += (c2[p] - c2[p + 1]).len();
+			sigma1 /= (num1 - 1);
+			sigma2 /= (num2 - 1);
+			float threshold = 4 * (sigma1 + sigma2) / 2;			// better way to do this?
+
+			T min_d = *std::min_element(dist.begin(), dist.end());	// find the minimum distance between c1 and c2
+
+			if (min_d > threshold) {								// if the minimum distance is too large
+				result.push_back(c1);
+				result.push_back(c2);
+
+#ifdef DEBUG
+				std::cout << "The distance between these two centerlines is too large" << std::endl;
+#endif
+			}
+			else {
+				auto smallest = std::min_element(dist.begin(), dist.end());
+				auto i = std::distance(dist.begin(), smallest);		// find the index of min-distance in distance list
+
+				// stitch position in c1 and c2
+				int id1 = index[i];
+				int id2 = i;
+
+				// actually there are two cases
+				// first one inacceptable
+				// second one acceptable
+				if (id1 != 0 && id1 != num1 - 1 && id2 != 0 && id2 != num2 - 1) {		// only stitch one end vertex to another centerline
+					result.push_back(c1);
+					result.push_back(c2);
+				}
+				else {
+					if (id1 == 0 || id1 == num1 - 1) {			// if the stitch vertex is the first or last vertex of c1
+						// for c2, consider two cases(one degenerate case)
+						if (id2 == 0 || id2 == num2 - 1) {		// case 1, if stitch position is also on the end of c2
+							// we have to decide which centerline get a new vertex, based on direction
+							// for c1, computer the direction change angle
+							stim::vec3<T> v1, v2;
+							float alpha1, alpha2;				// direction change angle
+							if (id1 == 0)
+								v1 = (c1[1] - c1[0]).norm();
+							else
+								v1 = (c1[num1 - 2] - c1[num1 - 1]).norm();
+							v2 = (c2[id2] - c1[id1]).norm();
+							alpha1 = v1.dot(v2);
+							if (id2 == 0)
+								v1 = (c2[1] - c2[0]).norm();
+							else
+								v1 = (c2[num2 - 2] - c2[num2 - 1]).norm();
+							v2 = (c1[id1] - c2[id2]).norm();
+							alpha2 = v1.dot(v2);
+							if (abs(alpha1) > abs(alpha2)) {					// add the vertex to c1 in order to get smooth connection
+								// push back c1
+								if (id1 == 0) {									// keep geometry information
+									new_vertex = c2[id2];
+									new_centerline.push_back(new_vertex);
+									for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
+										new_vertex = c1[p];
+										new_centerline.push_back(new_vertex);
+									}
+								}
+								else {
+									for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
+										new_vertex = c1[p];
+										new_centerline.push_back(new_vertex);
+									}
+									new_vertex = c2[id2];
+									new_centerline.push_back(new_vertex);
+								}
+								result.push_back(new_centerline);
+								new_centerline.clear();
+
+								// push back c2
+								for (size_t p = 0; p < num2; p++) {
+									new_vertex = c2[p];
+									new_centerline.push_back(new_vertex);
+								}
+								result.push_back(new_centerline);
+							}
+							else {												// add the vertex to c2 in order to get smooth connection
+								// push back c1
+								for (size_t p = 0; p < num1; p++) {
+									new_vertex = c1[p];
+									new_centerline.push_back(new_vertex);
+								}
+								result.push_back(new_centerline);
+								new_centerline.clear();
+
+								// push back c2
+								if (id2 == 0) {									// keep geometry information
+									new_vertex = c1[id1];
+									new_centerline.push_back(new_vertex);
+									for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
+										new_vertex = c2[p];
+										new_centerline.push_back(new_vertex);
+									}
+								}
+								else {
+									for (size_t p = 0; p < num2; p++) {			// stitch vertex on c2 -> geometry end vertex on c1 -> geometry start vertex on c1
+										new_vertex = c2[p];
+										new_centerline.push_back(new_vertex);
+									}
+									new_vertex = c1[id1];
+									new_centerline.push_back(new_vertex);
+								}
+								result.push_back(new_centerline);
+							}
+						}
+						else {												// case 2, the stitch position is on c2
+							// push back c1
+							if (id1 == 0) {									// keep geometry information
+								new_vertex = c2[id2];
+								new_centerline.push_back(new_vertex);
+								for (size_t p = 0; p < num1; p++) {			// stitch vertex on c2 -> geometry start vertex on c1 -> geometry end vertex on c1
+									new_vertex = c1[p];
+									new_centerline.push_back(new_vertex);
+								}
+							}
+							else {
+								for (size_t p = 0; p < num1; p++) {			// geometry end vertex on c1 -> geometry start vertex on c1 -> stitch vertex on c2
+									new_vertex = c1[p];
+									new_centerline.push_back(new_vertex);
+								}
+								new_vertex = c2[id2];
+								new_centerline.push_back(new_vertex);
+							}
+							result.push_back(new_centerline);
+							new_centerline.clear();
+
+							// push back c2
+							for (size_t p = 0; p < id2 + 1; p++) {			// first part
+								new_vertex = c2[p];
+								new_centerline.push_back(new_vertex);
+							}
+							result.push_back(new_centerline);
+							new_centerline.clear();
+
+							for (size_t p = id2; p < num2; p++) {			// second part
+								new_vertex = c2[p];
+								new_centerline.push_back(new_vertex);
+							}
+							result.push_back(new_centerline);
+						}
+					}
+					else {							// if the stitch vertex is the first or last vertex of c2
+						// push back c2
+						if (id2 == 0) {										// keep geometry information
+							new_vertex = c1[id1];
+							new_centerline.push_back(new_vertex);
+							for (size_t p = 0; p < num2; p++) {				// stitch vertex on c1 -> geometry start vertex on c2 -> geometry end vertex on c2
+								new_vertex = c2[p];
+								new_centerline.push_back(new_vertex);
+							}
+						}
+						else {
+							for (size_t p = 0; p < num2; p++) {				// geometry end vertex on c2 -> geometry start vertex on c2 -> stitch vertex on c1
+								new_vertex = c2[p];
+								new_centerline.push_back(new_vertex);
+							}
+							new_vertex = c1[id1];
+							new_centerline.push_back(new_vertex);
+							result.push_back(new_centerline);
+							new_centerline.clear();
+
+							// push back c1
+							for (size_t p = 0; p < id1 + 1; p++) {			// first part
+								new_vertex = c1[p];
+								new_centerline.push_back(new_vertex);
+							}
+							result.push_back(new_centerline);
+							new_centerline.clear();
+
+							for (size_t p = id1; p < num1; p++) {			// second part
+								new_vertex = c1[p];
+								new_centerline.push_back(new_vertex);
+							}
+							result.push_back(new_centerline);
+						}
+					}
+				}
+			}
+		}
+		return result;
+	}
+
 	/// Split the fiber at the specified index. If the index is an end point, only one fiber is returned
 	std::vector< stim::centerline<T> > split(unsigned int idx){
  
@@ -419,10 +419,10 @@ public:
 				file.read((char*)in, in_bytes);						//read an input spectrum
 				for (size_t b = 0; b < Bout; b++) {					//for each band
 					wc = bandlist[b];									//set the desired wavelength
-					band_bounds(wc, b0, b1);								//get the surrounding bands
+					hsi<T>::band_bounds(wc, b0, b1);								//get the surrounding bands
 					w0 = w[b0];											//get the wavelength for the lower band
 					w1 = w[b1];											//get the wavelength for the higher band
-					out[b] = lerp(wc, in[b0], w0, in[b1], w1);			//interpolate the spectral values to get the desired output value
+					out[b] = hsi<T>::lerp(wc, in[b0], w0, in[b1], w1);			//interpolate the spectral values to get the desired output value
 				}
 			}
 			else
@@ -63,8 +63,8 @@ protected:
 	}
  
 	/// Gets the two band indices surrounding a given wavelength
-	void band_bounds(double wavelength, unsigned long long& low, unsigned long long& high){
-		unsigned long long B = Z();
+	void band_bounds(double wavelength, size_t& low, size_t& high){
+		size_t B = Z();
 		for(high = 0; high < B; high++){
 			if(w[high] > wavelength) break;
 		}
@@ -92,7 +92,7 @@ protected:
 	/// @param s is the spectrum in main memory of length Z()
 	/// @param wavelength is the wavelength value to interpolate out
 	T interp_spectrum(T* s, double wavelength){
-		unsigned long long low, high;								//indices for the bands surrounding wavelength
+		size_t low, high;								//indices for the bands surrounding wavelength
 		band_bounds(wavelength, low, high);							//get the surrounding band indices
  
 		if(high == w.size()) return s[w.size()-1];					//if the high band is above the wavelength range, return the highest wavelength value
@@ -225,4 +225,4 @@ public:
  
 }		//end namespace STIM
  
-#endif
 \ No newline at end of file
+#endif