Jiaming Guo · David Mayerich · David Mayerich · David Mayerich · David Mayerich · Jiaming Guo
Showing 9 changed files Show diff stats
python/classify.py
python/digitalstain.py
python/envi.py
python/stim_spectral.py → python/hyperspectral.py
stim/biomodels/network.h
stim/grids/image_stack.h
stim/image/image.h
stim/visualization/cylinder.h
stim/visualization/gl_network.h
@@ -12,7 +12,7 @@ import sklearn.metrics
 import scipy
 import scipy.misc
 import envi
-import spectral
+import hyperspectral
 import random
 import progressbar
 import matplotlib.pyplot as plt
@@ -177,7 +177,7 @@ def envi_batch_predict(E, C, batch=10000):
         else:
             Tv = numpy.concatenate((Tv, C.predict(Fv.transpose()).transpose()), 0)
         tempmask = E.batchmask()
-        Lv = spectral.unsift2(Tv, tempmask)
+        Lv = hyperspectral.unsift2(Tv, tempmask)
         Cv = label2class(Lv.squeeze(), background=0)
         RGB = class2color(Cv)
         plt.imshow(RGB)
@@ -5,7 +5,7 @@ Created on Tue Jul 25 16:28:37 2017
 @author: david
 """
  
-import spectral
+import hyperspectral
 import envi
 import classify
 import numpy
@@ -18,22 +18,36 @@ import glob
 import matplotlib.pyplot as plt
 import random
  
-def generate_stain(envifile, stainfile, N=5000, batch_size=10000, validate=True):
-    E = envi.envi(envifile)
+def generate_stain(envifile, stainfile, maskfile="", trainmask="", N=5000, batch_size=10000, validate=True):
+    if trainmask == "":
+        E = envi.envi(envifile)
+    else:
+        mask = scipy.misc.imread(trainmask, flatten=True)
+        E = envi.envi(envifile, mask=mask)
+        
     mask = classify.random_mask(E.mask, N)
+    scipy.misc.imsave("random.bmp", mask)
  
     Ft = E.loadmask(mask).transpose()
  
     stain = numpy.rollaxis(scipy.misc.imread(stainfile), 2)
-    Tt = spectral.sift2(stain, mask).transpose()
+    Tt = hyperspectral.sift2(stain, mask).transpose()
  
+    print("Training MLPRegressor...")
     CLASS = sklearn.neural_network.MLPRegressor(solver='lbfgs', alpha=1e-5, hidden_layer_sizes=(), random_state=1, verbose=True)
     CLASS.fit(Ft, Tt)
  
     if validate == False:
         return CLASS
  
-    plt.ion()
+    print("Validating Stain...")
+    plt.ion()    
+    if not maskfile == "":
+        E.close()                                                                  #close the ENVI file
+        mask = scipy.misc.imread(maskfile, flatten=True)
+        print(numpy.count_nonzero(mask))
+        E = envi.envi(envifile, mask=mask)
+    
     Fv = E.loadbatch(batch_size)                                                #load the first batch
     n = 0
     while not Fv == []:                                                         #loop until an empty batch is returned
@@ -41,7 +55,7 @@ def generate_stain(envifile, stainfile, N=5000, batch_size=10000, validate=True)
             Tv = CLASS.predict(Fv.transpose()).transpose()
         else:
             Tv = numpy.append(Tv, CLASS.predict(Fv.transpose()).transpose(), 1)                        #append the predicted labels from this batch to those of previous batches
-        COLORS = spectral.unsift2(Tv, E.batchmask())                                    #convert the matrix of class labels to a 2D array
+        COLORS = hyperspectral.unsift2(Tv, E.batchmask())                                    #convert the matrix of class labels to a 2D array
         RGB = numpy.rollaxis(COLORS, 0, 3).astype(numpy.ubyte)
         plt.imshow(RGB)                                                             #display it
         plt.pause(0.05)
@@ -100,11 +100,11 @@ class envi_header:
             #t = map(str.strip, t)               #strip all of the tokens in the token list
  
             #handle the simple conditions
-            if l[li].startswith("file type"):
-                if not l[li].strip().endswith("ENVI Standard"):
-                    print("ERROR: unsupported ENVI file format: " + l[li].strip())
-                    return
-            elif l[li].startswith("samples"):
+            #if l[li].startswith("file type"):
+            #    if not l[li].strip().endswith("ENVI Standard"):
+            #        print("ERROR: unsupported ENVI file format: " + l[li].strip())
+            #        return
+            if l[li].startswith("samples"):
                 self.samples = int(l[li].split()[-1])
             elif l[li].startswith("lines"):
                 self.lines = int(l[li].split()[-1])
@@ -473,6 +473,146 @@ public:
 		return E[e];									//return the specified edge (casting it to a fiber)
 	}
  
+	/// subdivide current network
+	void subdivision() {
+
+		std::vector<unsigned> ori_index;		// original index
+		std::vector<unsigned> new_index;		// new index
+		std::vector<edge> nE;					// new edge
+		std::vector<vertex> nV;					// new vector
+		unsigned id = 0;
+
+		for (unsigned i = 0; i < num_edge; i++) {
+			if (E[i].size() == 2) {				// if current edge can't be subdivided
+				stim::centerline<T> line(2);
+				for (unsigned k = 0; k < 2; k++)
+					line[k] = E[i][k];
+				line.update();
+
+				edge new_edge(line);
+
+				vertex new_vertex = new_edge[0];
+				id = E[i].v[0];
+				auto position = std::find(ori_index.begin(), ori_index.end(), id);
+				if (position == ori_index.end()) {		 // new vertex
+					ori_index.push_back(id);
+					new_index.push_back(nV.size());
+
+					new_vertex.e[0].push_back(nE.size());
+					new_edge.v[0] = nV.size();
+					nV.push_back(new_vertex);			// push back vertex as a new vertex
+				}
+				else {									// existing vertex
+					int k = std::distance(ori_index.begin(), position);
+					new_edge.v[0] = new_index[k];
+					nV[new_index[k]].e[0].push_back(nE.size());
+				}
+
+				new_vertex = new_edge[1];
+				id = E[i].v[1];
+				position = std::find(ori_index.begin(), ori_index.end(), id);
+				if (position == ori_index.end()) {		 // new vertex
+					ori_index.push_back(id);
+					new_index.push_back(nV.size());
+
+					new_vertex.e[1].push_back(nE.size());
+					new_edge.v[1] = nV.size();
+					nV.push_back(new_vertex);			// push back vertex as a new vertex
+				}
+				else {									// existing vertex
+					int k = std::distance(ori_index.begin(), position);
+					new_edge.v[1] = new_index[k];
+					nV[new_index[k]].e[1].push_back(nE.size());
+				}
+
+				nE.push_back(new_edge);
+
+				nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(0));
+				nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(1));
+			}
+			else {								// subdivide current edge
+				for (unsigned j = 0; j < E[i].size() - 1; j++) {
+					stim::centerline<T> line(2);
+					for (unsigned k = 0; k < 2; k++)
+						line[k] = E[i][j + k];
+					line.update();
+
+					edge new_edge(line);
+
+					if (j == 0) {						// edge contains original starting point
+						vertex new_vertex = new_edge[0];
+						id = E[i].v[0];
+						auto position = std::find(ori_index.begin(), ori_index.end(), id);
+						if (position == ori_index.end()) {		 // new vertex
+							ori_index.push_back(id);
+							new_index.push_back(nV.size());
+
+							new_vertex.e[0].push_back(nE.size());
+							new_edge.v[0] = nV.size();
+							nV.push_back(new_vertex);			// push back vertex as a new vertex
+						}
+						else {									// existing vertex
+							int k = std::distance(ori_index.begin(), position);
+							new_edge.v[0] = new_index[k];
+							nV[new_index[k]].e[0].push_back(nE.size());
+						}
+
+						new_vertex = new_edge[1];
+						new_vertex.e[1].push_back(nE.size());
+						new_edge.v[1] = nV.size();
+						nV.push_back(new_vertex);				// push back internal point as a new vertex
+
+						nE.push_back(new_edge);
+					}
+
+					else if (j == E[i].size() - 2) {	// edge contains original ending point
+
+						vertex new_vertex = new_edge[1];
+						nV[nV.size() - 1].e[0].push_back(nE.size());
+						new_edge.v[0] = nV.size() - 1;
+
+						id = E[i].v[1];
+						auto position = std::find(ori_index.begin(), ori_index.end(), id);
+						if (position == ori_index.end()) {		 // new vertex
+							ori_index.push_back(id);
+							new_index.push_back(nV.size());
+
+							new_vertex.e[1].push_back(nE.size());
+							new_edge.v[1] = nV.size();
+							nV.push_back(new_vertex);			// push back vertex as a new vertex
+						}
+						else {									// existing vertex
+							int k = std::distance(ori_index.begin(), position);
+							new_edge.v[1] = new_index[k];
+							nV[new_index[k]].e[1].push_back(nE.size());
+						}
+
+						nE.push_back(new_edge);
+					}
+
+					else {
+						vertex new_vertex = new_edge[1];
+
+						nV[nV.size() - 1].e[0].push_back(nE.size());
+						new_vertex.e[1].push_back(nE.size());
+						new_edge.v[0] = nV.size() - 1;
+						new_edge.v[1] = nV.size();
+						nV.push_back(new_vertex);
+
+						nE.push_back(new_edge);
+					}
+
+					// get radii
+					nE[nE.size() - 1].cylinder<T>::set_r(0, E[i].cylinder<T>::r(j));
+					nE[nE.size() - 1].cylinder<T>::set_r(1, E[i].cylinder<T>::r(j + 1));
+				}
+			}
+		}
+
+		(*this).E = nE;
+		(*this).V = nV;
+	}
+
 	//load a network from an OBJ file
 	void load_obj(std::string filename){
  
@@ -715,7 +855,7 @@ public:
 			edge new_edge(C3);								// new edge	
  
 			//create an edge from the given centerline
-			unsigned int I = new_edge.size();				//calculate the number of points on the centerline
+			unsigned int I = (unsigned)new_edge.size();				//calculate the number of points on the centerline
  
 			//get the first and last vertex IDs for the line
 			i[0] = S.E[l].front();
@@ -1113,7 +1253,7 @@ public:
 		unsigned int id = 0;									// split value								
 		for(unsigned e = 0; e < E.size(); e++){					// for every edge
 			for(unsigned p = 0; p < E[e].size() - 1; p++){		// for every point in each edge
-				int t = E[e].length() / sigma * 2;
+				int t = (int)(E[e].length() / sigma * 2);
 				if (t <= 20)
 					threshold_fac = E[e].size();
 				else
@@ -206,9 +206,12 @@ public:
 	}
  
  
+	/* This was causing compiler errors. I don't think this function call exists anywhere:
+
 	void read(std::string file, unsigned int X, unsigned int Y, unsigned int Z, unsigned int C = 1, unsigned int header = 0){
 		read(file, stim::vec<unsigned long>(C, X, Y, Z), header);
 	}
+	*/
  
 	T* data(){
 		return ptr;
@@ -653,7 +653,7 @@ public:
  
 	//crop regions given by an array of 1D index values
 	std::vector< image<T> > crop_idx(size_t w, size_t h, std::vector<size_t> idx) {
-		std::vector<image<T>> result(idx.size());										//create an array of image files to return
+		std::vector< image<T> > result(idx.size());										//create an array of image files to return
 		for (size_t i = 0; i < idx.size(); i++) {										//for each specified index point
 			size_t y = idx[i] / X();													//calculate the y coordinate from the 1D index (center of ROI)
 			size_t x = idx[i] - y * X();												//calculate the x coordinate (center of ROI)
@@ -56,13 +56,12 @@ public:
 		init();
 	}
  
-/*
 	cylinder(stim::centerline<T> p, std::vector<T> s)
 	{
-		d = p;
+		//was d = s;
+		p = s;
 		init();
 	}
-*/
  
 	//cylinder(centerline<T>c, T r) : centerline(c) {
 	//	init();
 #ifndef STIM_GL_NETWORK
 #define STIM_GL_NETWORK
  
+#include <GL/glut.h>
 #include <stim/biomodels/network.h>
 #include <stim/visualization/aaboundingbox.h>