segmentation / netmets

Commit db598823d738e9bbd29e5546cf14c89df3869910

Authored by David Mayerich
1 parent 94a07643

fixed GCC errors and warnings associated with NetMets

Showing 1 changed file with 296 additions and 0 deletions   Show diff stats
Inline Side-by-side
main.cu 0 → 100644
  Show/Hide comments   View file @db59882
  1 +#include <stdlib.h>
  2 +#include <string>
  3 +#include <fstream>
  4 +#include <algorithm>
  5 +
  6 +//OpenGL includes
  7 +#include <GL/glut.h>
  8 +#include <GL/freeglut.h>
  9 +
  10 +//STIM includes
  11 +#include <stim/visualization/gl_network.h>
  12 +#include <stim/biomodels/network.h>
  13 +#include <stim/visualization/gl_aaboundingbox.h>
  14 +#include <stim/parser/arguments.h>
  15 +#include <stim/visualization/camera.h>
  16 +
  17 +//ANN includes
  18 +//#include <ANN/ANN.h>
  19 +
  20 +//BOOST includes
  21 +#include <boost/tuple/tuple.hpp>
  22 +
  23 +//visualization objects
  24 +stim::gl_aaboundingbox<float> bb; //axis-aligned bounding box object
  25 +stim::camera cam; //camera object
  26 +
  27 +unsigned num_nets = 0;
  28 +stim::gl_network<float> GT; //ground truth network
  29 +stim::gl_network<float> T; //test network
  30 +
  31 +//hard-coded parameters
  32 +float resample_rate = 0.5; //sample rate for the network (fraction of sigma used as the maximum sample rate)
  33 +float camera_factor = 1.2; //start point of the camera as a function of X and Y size
  34 +float orbit_factor = 0.01; //degrees per pixel used to orbit the camera
  35 +
  36 +//mouse position tracking
  37 +int mouse_x;
  38 +int mouse_y;
  39 +
  40 +//OpenGL objects
  41 +GLuint cmap_tex = 0; //texture name for the color map
  42 +
  43 +//sets an OpenGL viewport taking up the entire window
  44 +void glut_render_single_projection(){
  45 +
  46 + glMatrixMode(GL_PROJECTION); //load the projection matrix for editing
  47 + glLoadIdentity(); //start with the identity matrix
  48 + int X = glutGet(GLUT_WINDOW_WIDTH); //use the whole screen for rendering
  49 + int Y = glutGet(GLUT_WINDOW_HEIGHT);
  50 + glViewport(0, 0, X, Y); //specify a viewport for the entire window
  51 + float aspect = (float)X / (float)Y; //calculate the aspect ratio
  52 + gluPerspective(60, aspect, 0.1, 1000000); //set up a perspective projection
  53 +}
  54 +
  55 +//sets an OpenGL viewport taking up the left half of the window
  56 +void glut_render_left_projection(){
  57 +
  58 + glMatrixMode(GL_PROJECTION); //load the projection matrix for editing
  59 + glLoadIdentity(); //start with the identity matrix
  60 + int X = glutGet(GLUT_WINDOW_WIDTH) / 2; //only use half of the screen for the viewport
  61 + int Y = glutGet(GLUT_WINDOW_HEIGHT);
  62 + glViewport(0, 0, X, Y); //specify the viewport on the left
  63 + float aspect = (float)X / (float)Y; //calculate the aspect ratio
  64 + gluPerspective(60, aspect, 0.1, 1000000); //set up a perspective projection
  65 +}
  66 +
  67 +//sets an OpenGL viewport taking up the right half of the window
  68 +void glut_render_right_projection(){
  69 +
  70 + glMatrixMode(GL_PROJECTION); //load the projection matrix for editing
  71 + glLoadIdentity(); //start with the identity matrix
  72 + int X = glutGet(GLUT_WINDOW_WIDTH) / 2; //only use half of the screen for the viewport
  73 + int Y = glutGet(GLUT_WINDOW_HEIGHT);
  74 + glViewport(X, 0, X, Y); //specify the viewport on the right
  75 + float aspect = (float)X / (float)Y; //calculate the aspect ratio
  76 + gluPerspective(60, aspect, 0.1, 1000000); //set up a perspective projection
  77 +}
  78 +
  79 +void glut_render_modelview(){
  80 +
  81 + glMatrixMode(GL_MODELVIEW); //load the modelview matrix for editing
  82 + glLoadIdentity(); //start with the identity matrix
  83 + stim::vec3<float> eye = cam.getPosition(); //get the camera position (eye point)
  84 + stim::vec3<float> focus = cam.getLookAt(); //get the camera focal point
  85 + stim::vec3<float> up = cam.getUp(); //get the camera "up" orientation
  86 +
  87 + gluLookAt(eye[0], eye[1], eye[2], focus[0], focus[1], focus[2], up[0], up[1], up[2]); //set up the OpenGL camera
  88 +}
  89 +
  90 +//draws the network(s)
  91 +void glut_render(void) {
  92 +
  93 + if(num_nets == 1){ //if a single network is loaded
  94 + glut_render_single_projection(); //fill the entire viewport
  95 + glut_render_modelview(); //set up the modelview matrix with camera details
  96 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  97 + GT.glCenterline(GT.nmags() - 1); //render the GT network (the only one loaded)
  98 + }
  99 +
  100 + if(num_nets == 2){ //if two networks are loaded
  101 +
  102 + glut_render_left_projection(); //set up a projection for the left half of the window
  103 + glut_render_modelview(); //set up the modelview matrix using camera details
  104 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  105 +
  106 + glEnable(GL_TEXTURE_1D); //enable texture mapping
  107 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //texture map will be used as the network color
  108 + glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the Brewer texture map
  109 +
  110 + GT.glCenterline(GT.nmags() - 1); //render the GT network
  111 +
  112 + glut_render_right_projection(); //set up a projection for the right half of the window
  113 + glut_render_modelview(); //set up the modelview matrix using camera details
  114 + T.glCenterline(T.nmags() - 1); //render the T network
  115 +
  116 + }
  117 +
  118 + glutSwapBuffers();
  119 +}
  120 +
  121 +// defines camera motion based on mouse dragging
  122 +void glut_motion(int x, int y){
  123 +
  124 +
  125 + float theta = orbit_factor * (mouse_x - x); //determine the number of degrees along the x-axis to rotate
  126 + float phi = orbit_factor * (y - mouse_y); //number of degrees along the y-axis to rotate
  127 +
  128 + cam.OrbitFocus(theta, phi); //rotate the camera around the focal point
  129 +
  130 + mouse_x = x; //update the mouse position
  131 + mouse_y = y;
  132 +
  133 + glutPostRedisplay(); //re-draw the visualization
  134 +}
  135 +
  136 +// sets the mouse position when clicked
  137 +void glut_mouse(int button, int state, int x, int y){
  138 + mouse_x = x;
  139 + mouse_y = y;
  140 +}
  141 +
  142 +#define BREWER_CTRL_PTS 11 //number of control points in the Brewer map
  143 +void texture_initialize(){
  144 +
  145 + //define the colormap
  146 + static float brewer_map[BREWER_CTRL_PTS][3] = { //generate a Brewer color map (blue to red)
  147 + {0.192157f, 0.211765f, 0.584314f},
  148 + {0.270588f, 0.458824f, 0.705882f},
  149 + {0.454902f, 0.678431f, 0.819608f},
  150 + {0.670588f, 0.85098f, 0.913725f},
  151 + {0.878431f, 0.952941f, 0.972549f},
  152 + {1.0f, 1.0f, 0.74902f},
  153 + {0.996078f, 0.878431f, 0.564706f},
  154 + {0.992157f, 0.682353f, 0.380392f},
  155 + {0.956863f, 0.427451f, 0.262745f},
  156 + {0.843137f, 0.188235f, 0.152941f},
  157 + {0.647059f, 0.0f, 0.14902f}
  158 + };
  159 +
  160 + glGenTextures(1, &cmap_tex); //generate a texture map name
  161 + glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the texture map
  162 +
  163 + glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); //enable linear interpolation
  164 + glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  165 + glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP); //clamp the values at the minimum and maximum
  166 + glTexImage1D(GL_TEXTURE_1D, 0, 3, BREWER_CTRL_PTS, 0, GL_RGB, GL_FLOAT, //upload the texture map to the GPU
  167 + brewer_map);
  168 +}
  169 +
  170 +//Initialize the OpenGL (GLUT) window, including starting resolution, callbacks, texture maps, and camera
  171 +void glut_initialize(){
  172 +
  173 + int myargc = 1; //GLUT requires arguments, so create some bogus ones
  174 + char* myargv[1];
  175 + myargv [0]=strdup ("netmets");
  176 +
  177 + glutInit(&myargc, myargv); //pass bogus arguments to glutInit()
  178 + glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA); //generate a color buffer, depth buffer, and enable double buffering
  179 + glutInitWindowPosition(100,100); //set the initial window position
  180 + glutInitWindowSize(320,320); //set the initial window size
  181 + glutCreateWindow("NetMets - STIM Lab, UH"); //set the dialog box title
  182 +
  183 +
  184 + // register callback functions
  185 + glutDisplayFunc(glut_render); //function executed for rendering - renders networks
  186 + glutMouseFunc(glut_mouse); //executed on a mouse click - sets starting mouse positions for rotations
  187 + glutMotionFunc(glut_motion); //executed when the mouse is moved while a button is pressed
  188 +
  189 + texture_initialize(); //set up texture mapping (create texture maps, enable features)
  190 +
  191 + stim::vec3<float> c = bb.center(); //get the center of the network bounding box
  192 +
  193 + //place the camera along the z-axis at a distance determined by the network size along x and y
  194 + cam.setPosition(c + stim::vec<float>(0, 0, camera_factor * std::max(bb.size()[0], bb.size()[1])));
  195 + cam.LookAt(c[0], c[1], c[2]); //look at the center of the network
  196 +}
  197 +
  198 +//compare both networks and fill the networks with error information
  199 +void compare(float sigma){
  200 +
  201 + GT = GT.compare(T, sigma); //compare the ground truth to the test case - store errors in GT
  202 + T = T.compare(GT, sigma); //compare the test case to the ground truth - store errors in T
  203 +
  204 + //calculate the metrics
  205 + float FPR = GT.average(1); //calculate the metrics
  206 + float FNR = T.average(1);
  207 +
  208 + std::cout << "FNR: " << FPR << std::endl; //print false alarms and misses
  209 + std::cout << "FPR: " << FNR << std::endl;
  210 +}
  211 +
  212 +// 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
  213 +// Pranathi wrote this - saves network features to a CSV file
  214 +void features(std::string filename){
  215 + double avgL_t, avgL_gt, avgT_t, avgT_gt, avgB_t, avgB_gt, avgC_t, avgC_gt, avgFD_t, avgFD_gt;
  216 + unsigned int e_t, e_gt, b_gt, b_t;
  217 + avgL_gt = GT.Lengths();
  218 + avgT_gt = GT.Tortuosities();
  219 + avgL_t = T.Lengths();
  220 + avgT_t = T.Tortuosities();
  221 + avgB_gt = GT.BranchingIndex();
  222 + avgB_t = T.BranchingIndex();
  223 + avgC_gt = GT.Contractions();
  224 + avgFD_gt = GT.FractalDimensions();
  225 + avgC_t = T.Contractions();
  226 + avgFD_t = T.FractalDimensions();
  227 + e_gt = GT.EndP();
  228 + e_t = T.EndP();
  229 + b_gt = GT.BranchP();
  230 + b_t = T.BranchP();
  231 + std::ofstream myfile;
  232 + myfile.open (filename.c_str());
  233 + myfile << "Length, Tortuosity, Contraction, Fractal Dimension, Branch Points, End points, Branching Index, \n";
  234 + myfile << avgL_gt << "," << avgT_gt << "," << avgC_gt << "," << avgFD_gt << "," << b_gt << "," << e_gt << "," << avgB_gt <<std::endl;
  235 + myfile << avgL_t << "," << avgT_t << "," << avgC_t << "," << avgFD_t << "," << b_t << "," << e_t << "," << avgB_t <<std::endl;
  236 + myfile.close();
  237 +}
  238 +
  239 +// Output an advertisement for the lab, authors, and usage information
  240 +void advertise(){
  241 + std::cout<<std::endl<<std::endl;
  242 + std::cout<<"========================================================================="<<std::endl;
  243 + std::cout<<"Thank you for using the NetMets network comparison tool!"<<std::endl;
  244 + std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
  245 + std::cout<<"Developers: Pranathi Vemuri, David Mayerich"<<std::endl;
  246 + std::cout<<"Source: https://git.stim.ee.uh.edu/segmentation/netmets"<<std::endl;
  247 + std::cout<<"========================================================================="<<std::endl<<std::endl;
  248 +
  249 + std::cout<<"usage: netmets file1 file2 --sigma 10"<<std::endl;
  250 + std::cout<<" compare two files with a tolerance of 10 (units defined by the network)"<<std::endl;
  251 + std::cout<<" netmets file1 --gui"<<std::endl<<std::endl;
  252 + std::cout<<" load a file and display it using OpenGL"<<std::endl;
  253 +}
  254 +
  255 +int main(int argc, char* argv[])
  256 +{
  257 + stim::arglist args; //create an instance of arglist
  258 +
  259 + //add arguments
  260 + args.add("help", "prints this help");
  261 + args.add("sigma", "force a sigma value to specify the tolerance of the network comparison", "10");
  262 + args.add("gui", "display the network or network comparison using OpenGL");
  263 + args.add("features", "save features to a CSV file, specify file name");
  264 +
  265 + args.parse(argc, argv); //parse the user arguments
  266 +
  267 + if(args["help"].is_set() || args.nargs() == 0){ //test for help
  268 + advertise(); //output the advertisement
  269 + std::cout<<args.str(); //output arguments
  270 + exit(1); //exit
  271 + }
  272 +
  273 + if(args.nargs() >= 1){ //if at least one network file is specified
  274 + num_nets = 1; //set the number of networks to one
  275 + GT.load_obj(args.arg(0)); //load the specified file as the ground truth
  276 + /*GT.to_txt("Graph.txt");*/
  277 + }
  278 +
  279 + if(args.nargs() == 2){ //if two files are specified, they will be displayed in neighboring viewports and compared
  280 + num_nets = 2; //set the number of networks to two
  281 + float sigma = args["sigma"].as_float(); //get the sigma value from the user
  282 + T.load_obj(args.arg(1)); //load the second (test) network
  283 + if(args["features"].is_set()) //if the user wants to save features
  284 + features(args["features"].as_string());
  285 + GT = GT.resample(resample_rate * sigma); //resample both networks based on the sigma value
  286 + T = T.resample(resample_rate * sigma);
  287 + compare(sigma); //run the comparison algorithm
  288 + }
  289 +
  290 + //if a GUI is requested, display the network using OpenGL
  291 + if(args["gui"].is_set()){
  292 + bb = GT.boundingbox(); //generate a bounding volume
  293 + glut_initialize(); //create the GLUT window and set callback functions
  294 + glutMainLoop(); // enter GLUT event processing cycle
  295 + }
  296 +}
... ...