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

Pavel Govyadinov
2 parents c4887649 5687cf1b
Showing 12 changed files with 558 additions and 18 deletions Show diff stats
cmake/FindANN.cmake
FindSTIM.cmake → cmake/FindSTIM.cmake
stim/biomodels/network.h
stim/envi/bil.h
stim/envi/binary.h
stim/envi/bip.h
stim/envi/bsq.h
stim/envi/envi.h
stim/image/image.h
stim/math/mathvec.h → stim/math/vector.h
stim/visualization/colormap.h
stim/visualization/obj.h
+# - Try to find ANN
+# Once done this will define
+#  
+#  ANN_FOUND        - system has ANN
+#  ANN_INCLUDE_DIR  - the ANN include directory
+#  ANN_LIBRARY      - Link these to use ANN
+#   
+
+IF (ANN_INCLUDE_DIRS)
+  # Already in cache, be silent
+  SET(ANN_FIND_QUIETLY TRUE)
+ENDIF (ANN_INCLUDE_DIRS)
+
+FIND_PATH( ANN_INCLUDE_DIR ANN/ANN.h
+           PATHS "/usr/include" "C:/libs/ANN/include")
+
+if( WIN32 )
+
+	set(ANN_LIBRARY $ENV{ANN_PATH}\\ANN.lib)
+	set(ANN_INCLUDE_DIR $ENV{ANN_PATH})
+
+               
+	 # Store the library dir. May be used for linking to dll!
+	 # GET_FILENAME_COMPONENT( ANN_LIBRARY_DIR ${ANN_LIBRARY} PATH )
+
+	 find_package_handle_standard_args(ANN DEFAULT_MSG ANN_INCLUDE_DIR)
+ 
+else (WIN32)
+
+ FIND_LIBRARY( ANN_LIBRARY
+               NAMES ann ANN
+               PATHS /lib /usr/lib /usr/lib64 /usr/local/lib )
+
+endif( WIN32)
+
+
+IF (ANN_INCLUDE_DIR AND ANN_LIBRARY)
+  SET(ANN_FOUND TRUE)
+ELSE (ANN_INCLUDE_DIR AND ANN_LIBRARY)
+  SET( ANN_FOUND FALSE )
+ENDIF (ANN_INCLUDE_DIR AND ANN_LIBRARY)
+
+#ifndef STIM_NETWORK_H
+#define STIM_NETWORK_H
+
+#include <stim/math/vector.h>
+#include <stim/visualization/obj.h>
+#include <list>
+#include <ANN/ANN.h>
+
+namespace stim{
+
+/** This class provides an interface for dealing with biological networks.
+ *  It takes the following aspects into account:
+ *  	1) Network geometry and centerlines
+ *  	2) Network connectivity (a graph structure can be extracted)
+ *  	3) Network surface structure (the surface is represented as a triangular mesh and referenced to the centerline)
+ */
+
+template<typename T>
+class network{
+
+	//helper classes
+	/// Stores information about a geometric point on the network centerline (including point position and radius)
+	//template<typename T>
+	class point : public stim::vec<T>{
+
+	public:
+		T r;
+
+		point() : stim::vec<T>(){}
+
+		//casting constructor
+		point(stim::vec<T> rhs) : stim::vec<T>(rhs){}
+	};
+
+	//template<typename T>
+	class t_node;
+	class fiber;
+
+	//create typedefs for the iterators to simplify the network code
+	typedef typename std::list< fiber >::iterator fiber_i;
+	typedef typename std::list< t_node >::iterator t_node_i;
+
+	/// Stores information about a single capillary (a length of vessel between two branch or end points)
+	//template<typename T>
+	class fiber{
+
+	public:
+		std::list< point > P;		//geometric point positions
+
+		typename std::list< t_node >::iterator n[2];				//indices to terminal nodes
+		unsigned int id;
+
+	public:
+
+		/// Calculate the length of the fiber and return it.
+		T length(){
+
+			point p0, p1;
+			T l = 0;				//initialize the length to zero
+
+			//for each point
+			typename std::list< point >::iterator i;	//create a point iterator
+			for(i = P.begin(); i != P.end(); i++){		//for each point in the fiber
+
+				if(i == P.begin())						//if this is the first point, just store it
+					p1 = *i;
+				else{									//if this is any other point
+					p0 = p1;							//shift p1->p0
+					p1 = *i;							//set p1 to the new point
+					l += (p1 - p0).len();				//add the length of p1 - p0 to the running sum
+				}
+			}
+
+			return l;									//return the length
+		}
+
+		T radius(T& length){
+
+			point p0, p1;				//temporary variables to store point positions
+			T r0, r1;					//temporary variables to store radii at points
+			T l, r;						//temporary variable to store the length and average radius of a fiber segment
+			T length_sum = 0;			//initialize the length to zero
+			T radius_sum = 0;			//initialize the radius sum to zero
+
+			//for each point
+			typename std::list< point >::iterator i;	//create a point iterator
+			for(i = P.begin(); i != P.end(); i++){		//for each point in the fiber
+
+				if(i == P.begin()){						//if this is the first point, just store it
+					p1 = *i;
+					r1 = i->r;
+				}
+				else{									//if this is any other point
+					p0 = p1;							//shift p1->p0 and r1->r0
+					r0 = r1;
+					p1 = *i;							//set p1 to the new point
+					r1 = i->r;								//and r1
+
+					l = (p1 - p0).len();				//calculate the length of the p0-p1 segment
+					r = (r0 + r1) / 2;					//calculate the average radius of the segment
+
+					radius_sum += r * l;				//add the radius scaled by the length to a running sum
+					length_sum += l;					//add the length of p1 - p0 to the running sum
+				}
+			}
+
+			length = length_sum;						//store the total length
+			return radius_sum / length;					//return the average radius of the fiber
+		}
+
+		std::string str(){
+			std::stringstream ss;
+
+			//create an iterator for the point list
+			typename std::list<point>::iterator i;
+			for(i = P.begin(); i != P.end(); i++){
+				ss<<i->str()<<"  r = "<<i->r<<std::endl;
+			}
+
+			return ss.str();
+		}
+	};
+
+	/// Terminal node for a capillary. This is analogous to a graph vertex and contains a list of edge indices.
+	//template<typename T>
+	class t_node{
+
+	public:
+
+		unsigned int id;
+
+		//lists of edge indices for capillaries
+			//the "in" and "out" just indicate how the geometry is defined:
+			//		edges in the "in" list are geometrically oriented such that the terminal node is last
+			//		edges in the "out" list are geometrically oriented such that the terminal node is first
+		std::list< fiber_i > in;			//edge indices for incoming capillaries
+		std::list< fiber_i > out;		//edge indices for outgoing capillaries
+
+		std::string str(){
+
+			std::stringstream ss;
+
+			ss<<id<<": ";						//output the node ID
+
+			//output the IDs for both lists
+			typename std::list< fiber_i >::iterator f;
+
+			for(f = in.begin(); f != in.end(); f++){
+
+				if(f != in.begin())
+					ss<<", ";
+				ss<<(*f)->n[0]->id;
+			}
+
+			//if there are nodes in both lists, separate them by a comma
+			if(out.size() > 0 && in.size() > 0)
+				ss<<", ";
+
+			for(f = out.begin(); f != out.end(); f++){
+
+				if(f != out.begin())
+					ss<<", ";
+				ss<<(*f)->n[1]->id;
+			}
+
+
+			return ss.str();
+
+
+
+		}
+	};
+
+
+
+
+protected:
+
+	//list of terminal nodes
+	std::list<t_node> N;
+
+	//list of fibers
+	std::list<fiber> F;
+
+	/// Sets a unique ID for each terminal node and fiber
+	void set_names(){
+
+		unsigned int i;
+
+		i = 0;
+		for(t_node_i ti = N.begin(); ti != N.end(); ti++)
+			ti->id = i++;
+
+		i = 0;
+		for(fiber_i fi = F.begin(); fi != F.end(); fi++)
+			fi->id = i++;
+	}
+
+public:
+
+	std::string str(){
+
+
+		//assign names to elements of the network
+		set_names();
+
+		//create a stringstream for output
+		std::stringstream ss;
+
+		//output the nodes
+		ss<<"Nodes----------------------------"<<std::endl;
+		for(t_node_i i = N.begin(); i != N.end(); i++){
+			ss<<i->str()<<std::endl;
+		}
+
+		//output the fibers
+		ss<<std::endl<<"Fibers---------------------------"<<std::endl;
+
+		T length, radius;
+		//output every fiber
+		for(fiber_i f = F.begin(); f != F.end(); f++){
+
+			//calculate the length and average radius
+			radius = f->radius(length);
+
+			//output the IDs of the terminal nodes
+			ss<<f->n[0]->id<<" -- "<<f->n[1]->id<<": length = "<<length<<",  average radius = "<<radius<<std::endl;
+		}
+
+		return ss.str();
+	}
+
+	/// Load a network from an OBJ object
+	void load( stim::obj<T> object){
+
+		//get the number of vertices in the object
+		unsigned int nV = object.numV();
+
+		//allocate an array of pointers to nodes, which will be used to preserve connectivity
+			//initiate all values to T.end()
+		std::vector< t_node_i > node_hash(nV, N.end());
+
+		unsigned int nL = object.numL();			//get the number of lines in the OBJ
+
+		//for each line in the OBJ structure
+		for(unsigned int li = 0; li < nL; li++){
+
+			F.push_back(fiber());				//push a new fiber onto the fiber list
+
+			fiber_i f = --(F.end());			//get an iterator to the new fiber
+
+			//----------Handle the terminating nodes for the fiber
+
+			//get the indices of the line vertices
+			std::vector< unsigned int > Li = object.getL_Vi(li);
+			unsigned int i0 = Li.front() - 1;
+			unsigned int i1 = Li.back() - 1;
+
+			//deal with the first end point of the capillary
+			if(node_hash[i0] != N.end()){			//if the node has been used before
+				(*f).n[0] = node_hash[i0];			//assign the node to the new capillary
+				(*node_hash[i0]).out.push_back(f);	//add an out pointer to the existing node
+			}
+			else{									//otherwise
+				N.push_back(t_node());				//create a new node and add it to the node list
+				t_node_i t = --(N.end());			//get an iterator to the new node
+				node_hash[i0] = t;					//add a pointer to the new node to the hash list
+				(*f).n[0] = t;						//add a pointer to the new node to the capillary
+				(*t).out.push_back(f);				//add a pointer to the capillary to the new node
+			}
+
+			//deal with the last end point of the capillary
+			if(node_hash[i1] != N.end()){
+				(*f).n[1] = node_hash[i1];
+				(*node_hash[i1]).in.push_back(f);
+			}
+			else{
+				N.push_back(t_node());
+				t_node_i t = --(N.end());
+				node_hash[i1] = t;			//add the new node to the hash list
+				(*f).n[1] = t;
+				(*t).in.push_back(f);
+			}
+
+			//-------------Handle the geometric points for the fiber
+			std::vector< vec<T> > L = object.getL_V(li);
+			std::vector< vec<T> > R = object.getL_VT(li);
+
+			unsigned int nP = L.size();				//get the number of geometric points in the fiber
+			//for each vertex in the fiber
+			for(unsigned int pi = 0; pi < nP; pi++){
+				point p = (point)L[pi];					//move the geometric coordinates into a point structure
+				p.r = R[pi][0];							//store the radius
+				f->P.push_back(p);						//push the point onto the current fiber
+			}
+		}
+
+	}	//end load()
+
+	/// This function returns the information necessary for a simple graph-based physical (ex. fluid) simulation.
+
+	/// @param n0 is a array which will contain the list of source nodes
+	/// @param n1 is a array which will contain the list of destination nodes
+	/// @param length is a array containing the lengths of fibers in the network
+	/// @param radius is a array containing the average radii of fibers in the network
+	void build_simgraph(std::vector<unsigned int>& n0, std::vector<unsigned int>& n1, std::vector<T>& length, std::vector<T>& radius){
+
+		//determine the number of fibers in the network
+		unsigned int nF = F.size();
+
+		//allocate the necessary space to store the fiber information
+		n0.resize(nF);
+		n1.resize(nF);
+		length.resize(nF);
+		radius.resize(nF);
+
+		//assign names (identifiers) to the network components
+		set_names();
+
+		//fill the arrays
+		unsigned int i = 0;
+		T l, r;
+		for(fiber_i f = F.begin(); f != F.end(); f++){
+			n0[i] = f->n[0]->id;	//get the identifiers for the first and second nodes for the current fiber
+			n1[i] = f->n[1]->id;
+
+			r = f->radius(l);		//get the length and radius of the capillary (calculated at the same time)
+
+			radius[i] = r;			//store the radius in the output array
+			length[i] = l;			//store the length in the output array
+
+			i++;					//increment the array index
+		}
+
+
+	}
+
+};
+
+};	//end namespace stim
+
+
+#endif
@@ -826,7 +826,7 @@ public:
 	}
  
 	///Saves to disk only those spectra corresponding to mask values != 0
-	bool sift_mask(std::string outfile, unsigned char* p){
+	bool sift(std::string outfile, unsigned char* p){
 		// Assume X() = X, Y() = Y, Z() = Z.
 		std::ofstream target(outfile.c_str(), std::ios::binary);
  
@@ -100,7 +100,7 @@ public:
 	/// @param filename is the name of the binary file
 	/// @param r is a STIM vector specifying the size of the binary file along each dimension
 	/// @param h is the length (in bytes) of any header file (default zero)
-	bool open(std::string filename, vec<unsigned int, D> r, unsigned int h = 0){
+	bool open(std::string filename, vec<unsigned int> r, unsigned int h = 0){
  
 		for(unsigned int i = 0; i < D; i++)		//set the dimensions of the binary file object
 			R[i] = r[i];
@@ -117,7 +117,7 @@ public:
 	/// @param filename is the name of the binary file to be created
 	/// @param r is a STIM vector specifying the size of the file along each dimension
 	/// @offset specifies how many bytes to offset the file (used to leave room for a header)
-	bool create(std::string filename, vec<unsigned int, D> r, unsigned int offset = 0){
+	bool create(std::string filename, vec<unsigned int> r, unsigned int offset = 0){
  
 		std::ofstream target(filename.c_str(), std::ios::binary);
  
@@ -836,8 +836,8 @@ public:
 	}
  
  
-	///Saves to disk only those spectra corresponding to mask values != 0
-	bool sift_mask(std::string outfile, unsigned char* p){
+	/// Saves to disk only those spectra corresponding to mask values != 0
+	bool sift(std::string outfile, unsigned char* p){
 		// Assume X() = X, Y() = Y, Z() = Z.
 		std::ofstream target(outfile.c_str(), std::ios::binary);
  
@@ -759,8 +759,8 @@ public:
 		return true;
 	}
  
-	///Saves to disk only those spectra corresponding to mask values != 0
-	bool sift_mask(std::string outfile, unsigned char* p){
+	/// Saves to disk only those spectra corresponding to mask values != 0
+	bool sift(std::string outfile, unsigned char* p){
 		std::ofstream target(outfile.c_str(), std::ios::binary);
 		// open a band (XY plane)
 		unsigned XY = X() * Y(); //Number of XY pixels
@@ -788,6 +788,57 @@ public:
 		return true;
 	}
  
+	/// Generates a spectral image from a matrix of spectral values in lexicographic order and a mask
+	bool unsift(std::string outfile, unsigned char* p, unsigned int samples, unsigned int lines){
+
+		//create a binary output stream
+		std::ofstream target(outfile.c_str(), std::ios::binary);
+
+		//make sure that there's only one line
+		if(Y() != 1){
+			std::cout<<"ERROR in stim::bsq::sift() - number of lines does not equal 1"<<std::endl;
+			return false;
+		}
+
+		std::cout<<"started sifting"<<std::endl;
+
+		//get the number of pixels and bands in the input image
+		unsigned int P = X(); 					//Number of pixels
+		unsigned int B = Z();					//number of bands
+		unsigned int XY = samples * lines;		//total number of pixels in an unsifted image
+
+		// allocate memory for a sifted band
+		T * sifted = (T*)malloc(P * sizeof(T)); //allocate memory for a band
+
+		//allocate memory for an unsifted band image
+		T* unsifted = (T*) malloc(XY * sizeof(T));
+
+		//for each band
+		for(unsigned int b = 0; b < B; b++){
+
+			//set the unsifted index value to zero
+			unsigned int i = 0;
+
+			//retrieve the sifted band (masked pixels only)
+			band_index(sifted, b);
+
+			//for each pixel in the final image (treat it as a 1D image)
+			for(unsigned int xi = 0; xi < XY; xi++){
+				if( p[xi] == 0 )
+					unsifted[xi] = 0;
+				else{
+					unsifted[xi] = sifted[i];
+					i++;
+				}
+			}
+
+			//write the band image to disk
+			target.write(reinterpret_cast<const char*>(unsifted), sizeof(T) * XY);
+		}
+
+		return true;
+	}
+
  
 	/// Calculate the mean band value (average along B) at each pixel location.
  
@@ -373,32 +373,48 @@ public:
 	}
  
 	/// sift-mask saves in an array only those spectra corresponding to nonzero values of the mask. 
-	bool sift_mask(std::string outfile, unsigned char* p)
+	bool sift(std::string outfile, unsigned char* p)
 	{
  
+		//calculate the number of non-zero values in the mask
+		unsigned int nnz = 0;
+		unsigned int npixels = header.lines * header.samples;
+		for(unsigned int i = 0; i < npixels; i++)
+			if( p[i] > 0 ) nnz++;
+
+		//create a new header
+		envi_header new_header = header;
+
+		//set the number of lines to 1 (this is a matrix with 1 line and N samples)
+		new_header.lines = 1;
+		new_header.samples = nnz;
+		new_header.save(outfile + ".hdr");
+
+		std::cout<<"Saved: "<<outfile+".hdr"<<std::endl;
+
 		if (header.interleave == envi_header::BSQ){		//if the infile is bsq file
 			if (header.data_type == envi_header::float32)
-				return ((bsq<float>*)file)->sift_mask(outfile, p);
+				return ((bsq<float>*)file)->sift(outfile, p);
 			else if (header.data_type == envi_header::float64)
-				return ((bsq<double>*)file)->sift_mask(outfile, p);
+				return ((bsq<double>*)file)->sift(outfile, p);
 			else
 				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
  
 		else if (header.interleave == envi_header::BIL){		//if the infile is bil file
 			if (header.data_type == envi_header::float32)
-				return ((bil<float>*)file)->sift_mask(outfile, p);
+				return ((bil<float>*)file)->sift(outfile, p);
 			else if (header.data_type == envi_header::float64)
-				return ((bil<double>*)file)->sift_mask(outfile, p);
+				return ((bil<double>*)file)->sift(outfile, p);
 			else
 				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
  
 		else if (header.interleave == envi_header::BIP){		//if the infile is bip file
 			if (header.data_type == envi_header::float32)
-				return ((bip<float>*)file)->sift_mask(outfile, p);
+				return ((bip<float>*)file)->sift(outfile, p);
 			else if (header.data_type == envi_header::float64)
-				return ((bip<double>*)file)->sift_mask(outfile, p);
+				return ((bip<double>*)file)->sift(outfile, p);
 			else
 				std::cout << "ERROR: unidentified data type" << std::endl;
 		}
@@ -407,9 +423,48 @@ public:
 			std::cout << "ERROR: unidentified file type" << std::endl;
 			exit(1);
 		}
+
+		
 		return false;
 	}
  
+	bool unsift(std::string outfile, unsigned char* mask, unsigned int samples, unsigned int lines){
+
+		//create a new header
+		envi_header new_header = header;
+
+		//set the number of lines and samples in the output file (that's all that changes)
+		new_header.lines = lines;
+		new_header.samples = samples;
+		new_header.save(outfile + ".hdr");
+
+
+		if (header.interleave == envi_header::BSQ){		//if the infile is bsq file
+			if (header.data_type == envi_header::float32)
+				return ((bsq<float>*)file)->unsift(outfile, mask, samples, lines);
+			else if (header.data_type == envi_header::float64)
+				return ((bsq<double>*)file)->unsift(outfile, mask, samples, lines);
+			else
+				std::cout << "ERROR: unidentified data type" << std::endl;
+		}
+
+		else if (header.interleave == envi_header::BIL){		//if the infile is bil file
+			
+				std::cout << "ERROR in stim::envi::unsift - BIL files aren't supported yet" << std::endl;
+		}
+
+		else if (header.interleave == envi_header::BIP){		//if the infile is bip file
+			
+				std::cout << "ERROR in stim::envi::unsift - BIP files aren't supported yet" << std::endl;
+		}
+
+		else{
+			std::cout << "ERROR: unidentified file type" << std::endl;
+			exit(1);
+		}
+
+	}
+
 	/// Compute the ratio of two baseline-corrected peaks. The result is stored in a pre-allocated array.
  
 	/// @param lb1 is the label value for the left baseline point for the first peak (numerator)
@@ -28,6 +28,11 @@ public:
 		img.load(filename.c_str());
 	}
  
+	/// Constructor initializes an image to a given size
+	image(unsigned int x, unsigned int y = 1, unsigned int z = 1){
+		img = cimg_library::CImg<T>(x, y, z);
+	}
+
 	//Load an image from a file
 	void load(std::string filename){
 		img.load(filename.c_str());
@@ -72,6 +77,17 @@ public:
 				data[x * C + c] = ptr[c * X + x];
 	}
  
+	image<T> channel(unsigned int c){
+
+		//create a new image
+		image<T> single;
+
+		single.img = img.channel(c);
+
+		return single;
+
+	}
+
 	unsigned int channels(){
 		return (unsigned int)img.spectrum();
 	}
@@ -287,7 +287,7 @@ static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, T 
 }
  
 template<class T>
-static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale)//, bool positive = false)
+static void cpu2cpu(T* cpuSource, unsigned char* cpuDest, unsigned int nVals, colormapType cm = cmGrayscale)
 {
     //computes the max and min range automatically
  
@@ -329,13 +329,13 @@ static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, u
 }
  
 template<typename T>
-static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, colormapType cm = cmGrayscale, bool positive = false)
+static void cpu2image(T* cpuSource, std::string fileDest, unsigned int x_size, unsigned int y_size, colormapType cm = cmGrayscale)
 {
     //allocate a color buffer
 	unsigned char* cpuBuffer = (unsigned char*) malloc(sizeof(unsigned char) * 3 * x_size * y_size);
  
 	//do the mapping
-	cpu2cpu<T>(cpuSource, cpuBuffer, x_size * y_size, cm, positive);
+	cpu2cpu<T>(cpuSource, cpuBuffer, x_size * y_size, cm);
  
 	//copy the buffer to an image
 	buffer2image(cpuBuffer, fileDest, x_size, y_size);
@@ -6,7 +6,7 @@
 #include <fstream>
 #include <stdlib.h>
 #include <stim/parser/parser.h>
-#include <stim/math/mathvec.h>
+#include <stim/math/vector.h>
  
 namespace stim{
  
@@ -600,9 +600,42 @@ public:
 		}
  
 		return l;
+	}
+
+	/// Returns a vector containing the list of texture coordinates associated with each point of a line.
+
+	/// @param i is the index of the desired line
+	std::vector< stim::vec<T> > getL_VT(unsigned int i){
+
+		//get the number of points in the specified line
+		unsigned int nP = L[i].size();
+
+		//create a vector of points
+		std::vector< stim::vec<T> > l;
+
+		//set the size of the vector
+		l.resize(nP);
+
+		//copy the points from the point list to the stim vector
+		unsigned int pi;
+		for(unsigned int p = 0; p < nP; p++){
+
+			//get the index of the geometry point
+			pi = L[i][p][1] - 1;
+
+			//get the coordinates of the current point
+			stim::vec<T> newP = VT[pi];
+
+			//copy the point into the vector
+			l[p] = newP;
+		}
+
+		return l;
  
 	}
  
+
+
 };