diff --git a/main.cpp b/main.cpp
index 1e8293a..7c7526d 100644
--- a/main.cpp
+++ b/main.cpp
@@ -8,6 +8,7 @@
 
 //STIM includes
 #include <stim/visualization/gl_network.h>
+#include <stim/biomodels/network.h>
 #include <stim/visualization/gl_aaboundingbox.h>
 #include <stim/parser/arguments.h>
 #include <stim/visualization/camera.h>
@@ -207,6 +208,46 @@ void compare(float sigma){
 	std::cout << "FPR: " << FNR << std::endl;
 }
 
+// writes features of the networks i.e average segment length, tortuosity, branching index, contraction, fractal dimension, number of end and branch points to a csv file
+void features(){
+			double avgL_t, avgL_gt, avgT_t, avgT_gt, avgB_t, avgB_gt, avgC_t, avgC_gt, avgFD_t, avgFD_gt;
+		unsigned int e_t, e_gt, b_gt, b_t;
+		avgL_gt = GT.Lengths();
+		avgT_gt = GT.Tortuosities();
+		//std::cout<<avgL_gt<<"---average segment length in the ground truth network"<<std::endl;
+		//std::cout<<avgT_gt<<"---average segment tortuosity in the ground truth network"<<std::endl;
+		avgL_t = T.Lengths();
+		avgT_t = T.Tortuosities();
+		//std::cout<<avgL_t<<"---average segment length in the truth network"<<std::endl;
+		//std::cout<<avgT_t<<"---average segment tortuosity in the truth network"<<std::endl;
+		avgB_gt = GT.BranchingIndex();
+		//std::cout<<avgB_gt<<"---average branching index in the ground truth"<<std::endl;
+		avgB_t = T.BranchingIndex();
+		//std::cout<<avgB_t<<"---average branching index in the truth case"<<std::endl;
+		avgC_gt = GT.Contractions();
+		avgFD_gt = GT.FractalDimensions();
+		//std::cout<<avgC_gt<<"---average segment contraction in the ground truth network"<<std::endl;
+		//std::cout<<avgFD_gt<<"---average segment fractal dimension in the ground truth network"<<std::endl;
+		avgC_t = T.Contractions();
+		avgFD_t = T.FractalDimensions();
+		//std::cout<<avgC_t<<"---average segment contraction in the truth network"<<std::endl;
+		//std::cout<<avgFD_t<<"---average segment fractal dimension in the truth network"<<std::endl;
+		e_gt = GT.EndP();
+		//std::cout<<e_gt<<"---Number of tips in the ground truth"<<std::endl;
+		e_t = T.EndP();
+		//std::cout<<e_t<<"---Number of tips in the truth case"<<std::endl;
+		b_gt = GT.BranchP();
+		//std::cout<<b_gt<<"---Number of branch points in the ground truth"<<std::endl;
+		b_t = T.BranchP();
+		//std::cout<<b_t<<"---Number of branch points in the truth case"<<std::endl;
+		std::ofstream myfile;
+		myfile.open ("features.csv");
+		myfile << "Length, Tortuosity, Contraction, Fractal Dimension, Branch Points, End points, Branching Index, \n";
+		myfile << avgL_gt << "," << avgT_gt << "," << avgC_gt << "," << avgFD_gt << "," << b_gt << "," << e_gt << "," << avgB_gt <<std::endl;
+		myfile << avgL_t << "," << avgT_t << "," << avgC_t << "," << avgFD_t << "," << b_t << "," << e_t << "," << avgB_t <<std::endl;
+		myfile.close();
+}
+
 // Output an advertisement for the lab, authors, and usage information
 void advertise(){
 	std::cout<<std::endl<<std::endl;
@@ -243,24 +284,16 @@ int main(int argc, char* argv[])
 	if(args.nargs() >= 1){					//if at least one network file is specified
 		num_nets = 1;						//set the number of networks to one
 		GT.load_obj(args.arg(0));			//load the specified file as the ground truth
+		/*GT.to_txt("Graph.txt");*/
 	}
 	
 	if(args.nargs() == 2){			//if two files are specified, they will be displayed in neighboring viewports and compared
 		num_nets = 2;										//set the number of networks to two
 		float sigma = args["sigma"].as_float();				//get the sigma value from the user
-		T.load_obj(args.arg(1));							//load the second (test) network
-		double avgL, avgT;
-		avgL = GT.Lengths();
-		avgT = GT.Tortuosities();
-		std::cout<<avgL<<"---average segment length in the ground truth network"<<std::endl;
-		std::cout<<avgT<<"---average segment tortuosity in the ground truth network"<<std::endl;
-		avgL = T.Lengths();
-		avgT = T.Tortuosities();
-		std::cout<<avgL<<"---average segment length in the truth network"<<std::endl;
-		std::cout<<avgT<<"---average segment tortuosity in the truth network"<<std::endl;
+		T.load_obj(args.arg(1));                           //load the second (test) network
+		features();
 		GT = GT.resample(resample_rate * sigma);			//resample both networks based on the sigma value
 		T = T.resample(resample_rate * sigma);
-
 		compare(sigma);										//run the comparison algorithm
 	}
 
--
libgit2 0.21.4