segmentation / netmets

Commit d8e684dca35baa092e86aa83ed9ae11d2cda2a20

Authored by David Mayerich
2 parents bf7e73b2 5ac53c9e

Merge branch 'JACK' into 'master' 

Jack

mapping done

See merge request !2
Showing 2 changed files with 540 additions and 33 deletions   Show diff stats
Inline Side-by-side
1   -#include <stdlib.h>
  1 +#include <stdlib.h>
2 2 #include <string>
3 3 #include <fstream>
4 4 #include <algorithm>
... ... @@ -22,6 +22,9 @@
22 22 //ANN includes
23 23 //#include <ANN/ANN.h>
24 24  
  25 +//chrono includes
  26 +//#include <chrono>
  27 +
25 28 //BOOST includes
26 29 #include <boost/tuple/tuple.hpp>
27 30  
... ... @@ -32,16 +35,40 @@ stim::camera cam; //camera object
32 35 unsigned num_nets = 0;
33 36 stim::gl_network<float> GT; //ground truth network
34 37 stim::gl_network<float> T; //test network
  38 +stim::gl_network<float> _GT; //splitted GT
  39 +stim::gl_network<float> _T; //splitted T
  40 +
  41 +unsigned ind = 0; //indicator of mapping
  42 +
  43 +std::vector<unsigned> _gt_t; // store indices of nearest edge points in _T for _GT
  44 +std::vector<unsigned> _t_gt; // store indices of nearest edge points in _GT for _T
35 45  
36 46 //hard-coded parameters
37 47 float resample_rate = 0.5f; //sample rate for the network (fraction of sigma used as the maximum sample rate)
38 48 float camera_factor = 1.2f; //start point of the camera as a function of X and Y size
39 49 float orbit_factor = 0.01f; //degrees per pixel used to orbit the camera
40 50  
  51 +//mouse click
  52 +bool LButtonDown = false; // true when left button down
  53 +bool RButtonDown = false;
  54 +
41 55 //mouse position tracking
42 56 int mouse_x;
43 57 int mouse_y;
44 58  
  59 +bool compareMode = true; // default mode is compare mode
  60 +bool mappingMode = false;
  61 +
  62 +// random color set
  63 +std::vector<float> colormap;
  64 +
  65 +// special key indicator
  66 +int mods;
  67 +
  68 +// create display lists
  69 +GLuint dlist1;
  70 +GLuint dlist2;
  71 +
45 72 //OpenGL objects
46 73 GLuint cmap_tex = 0; //texture name for the color map
47 74  
... ... @@ -95,37 +122,102 @@ void glut_render_modelview(){
95 122 //draws the network(s)
96 123 void glut_render(void) {
97 124  
98   - if(num_nets == 1){ //if a single network is loaded
99   - glut_render_single_projection(); //fill the entire viewport
100   - glut_render_modelview(); //set up the modelview matrix with camera details
101   - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
102   - GT.glCenterline(GT.nmags() - 1); //render the GT network (the only one loaded)
103   - }
104   -
105   - if(num_nets == 2){ //if two networks are loaded
  125 + if(ind == 0){
  126 + if(num_nets == 1){ //if a single network is loaded
  127 + glut_render_single_projection(); //fill the entire viewport
  128 + glut_render_modelview(); //set up the modelview matrix with camera details
  129 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  130 + GT.glCenterline0(); //render the GT network (the only one loaded)
  131 + }
106 132  
107   - glut_render_left_projection(); //set up a projection for the left half of the window
108   - glut_render_modelview(); //set up the modelview matrix using camera details
109   - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  133 + if(num_nets == 2){ //if two networks are loaded
110 134  
111   - glEnable(GL_TEXTURE_1D); //enable texture mapping
112   - glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //texture map will be used as the network color
113   - glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the Brewer texture map
  135 + glut_render_left_projection(); //set up a projection for the left half of the window
  136 + glut_render_modelview(); //set up the modelview matrix using camera details
  137 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
114 138  
115   - GT.glCenterline(GT.nmags() - 1); //render the GT network
  139 + glEnable(GL_TEXTURE_1D); //enable texture mapping
  140 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //texture map will be used as the network color
  141 + glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the Brewer texture map
116 142  
117   - glut_render_right_projection(); //set up a projection for the right half of the window
118   - glut_render_modelview(); //set up the modelview matrix using camera details
119   - T.glCenterline(T.nmags() - 1); //render the T network
  143 + GT.glCenterline(GT.nmags() - 1); //render the GT network
120 144  
  145 + glut_render_right_projection(); //set up a projection for the right half of the window
  146 + glut_render_modelview(); //set up the modelview matrix using camera details
  147 + T.glCenterline(T.nmags() - 1); //render the T network
  148 + }
  149 + }
  150 + else{
  151 + if(num_nets == 1){ //if a single network is loaded
  152 + std::cout << "You should have at least two networks to do mapping." << std::endl; //exit program because there isn't enough network
  153 + exit(1);
  154 + }
  155 + if(num_nets == 2){ //if two networks are loaded
  156 + if(compareMode){
  157 + glut_render_left_projection(); //set up a projection for the left half of the window
  158 + glut_render_modelview(); //set up the modelview matrix using camera details
  159 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  160 +
  161 + glEnable(GL_TEXTURE_1D); //enable texture mapping
  162 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //texture map will be used as the network color
  163 + glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the Brewer texture map
  164 +
  165 + _GT.glCenterline(_GT.nmags() - 1); //render the GT network
  166 +
  167 + glut_render_right_projection(); //set up a projection for the right half of the window
  168 + glut_render_modelview(); //set up the modelview matrix using camera details
  169 + _T.glCenterline(_T.nmags() - 1); //render the T network
  170 +
  171 + }
  172 + else{
  173 + glut_render_left_projection(); //set up a projection for the left half of the window
  174 + glut_render_modelview(); //set up the modelview matrix using camera details
  175 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  176 +
  177 + _GT.glRandColorCenterlineGT(dlist1, _gt_t, colormap);
  178 +
  179 + glut_render_right_projection(); //set up a projection for the right half of the window
  180 + glut_render_modelview(); //set up the modelview matrix using camera details
  181 + _T.glRandColorCenterlineT(dlist2, _t_gt, colormap);
  182 + }
  183 + }
121 184 }
122 185  
  186 + if(num_nets == 2){
  187 + std::ostringstream ss;
  188 + if (mappingMode) // if it is in mapping mode
  189 + ss << "Mapping Mode";
  190 + else
  191 + ss << "Compare Mode"; // default mode is compare mode
  192 +
  193 + glDisable(GL_TEXTURE_1D);
  194 + glMatrixMode(GL_PROJECTION); //Set up the 2d viewport for mode text printing
  195 + glPushMatrix();
  196 + glLoadIdentity();
  197 + int X = glutGet(GLUT_WINDOW_WIDTH);
  198 + int Y = glutGet(GLUT_WINDOW_HEIGHT);
  199 + glViewport(0, 0, X / 2, Y); // locate to left bottom corner
  200 + gluOrtho2D(0, X, 0, Y); // define othogonal aspect
  201 + glColor3f(0.0, 1.0, 0.0); // using green to show mode
  202 +
  203 + glMatrixMode(GL_MODELVIEW);
  204 + glPushMatrix();
  205 + glLoadIdentity();
  206 +
  207 + glRasterPos2f(0, 5); //print text in the bottom left corner
  208 + glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_24, (const unsigned char*)(ss.str().c_str()));
  209 +
  210 + glPopMatrix();
  211 + glMatrixMode(GL_PROJECTION);
  212 + glPopMatrix();
  213 + }
123 214 glutSwapBuffers();
124 215 }
125 216  
126 217 // defines camera motion based on mouse dragging
127 218 void glut_motion(int x, int y){
128 219  
  220 + if(LButtonDown == true && RButtonDown == false && mods != GLUT_ACTIVE_CTRL){
129 221  
130 222 float theta = orbit_factor * (mouse_x - x); //determine the number of degrees along the x-axis to rotate
131 223 float phi = orbit_factor * (y - mouse_y); //number of degrees along the y-axis to rotate
... ... @@ -136,12 +228,71 @@ void glut_motion(int x, int y){
136 228 mouse_y = y;
137 229  
138 230 glutPostRedisplay(); //re-draw the visualization
  231 + }
139 232 }
140 233  
141 234 // sets the mouse position when clicked
142 235 void glut_mouse(int button, int state, int x, int y){
143   - mouse_x = x;
144   - mouse_y = y;
  236 +
  237 + if(button == GLUT_LEFT_BUTTON && state == GLUT_DOWN){
  238 + mouse_x = x;
  239 + mouse_y = y;
  240 + LButtonDown = true;
  241 + }
  242 + else if(button == GLUT_RIGHT_BUTTON && state == GLUT_DOWN){
  243 + mouse_x = x;
  244 + mouse_y = y;
  245 + RButtonDown = true;
  246 + }
  247 + else if(button == GLUT_LEFT_BUTTON && state == GLUT_UP){
  248 + mouse_x = x;
  249 + mouse_y = y;
  250 + LButtonDown = false;
  251 + }
  252 + else if(button == GLUT_RIGHT_BUTTON && state == GLUT_UP){
  253 + mouse_x = x;
  254 + mouse_y = y;
  255 + RButtonDown = false;
  256 + }
  257 +
  258 + /// implementation of mouse click mapping feedback
  259 + mods = glutGetModifiers(); // get modifier keys
  260 + if (mods == GLUT_ACTIVE_CTRL) // if the CTRL key is pressed
  261 + if (button == GLUT_LEFT_BUTTON && state == GLUT_DOWN) {
  262 + std::cout << "( " << x << ", " << y << " )" << std::endl; // if the CTRL key is pressed and LEFT BUTTON is DOWN, print the window coordinates
  263 +
  264 + GLint viewport[4];
  265 + GLdouble modelview[16];
  266 + GLdouble projection[16];
  267 + GLdouble winX, winY, winZ;
  268 + GLdouble posX, posY, posZ;
  269 +
  270 + glGetIntegerv(GL_VIEWPORT, viewport);
  271 + glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
  272 + glGetDoublev(GL_PROJECTION_MATRIX, projection);
  273 +
  274 + winX = (GLdouble)x;
  275 + winY = viewport[3] - (GLdouble)y;
  276 + glReadPixels((GLint)winX, (GLint)winY, (GLsizei)1, (GLsizei)1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ); // need frame buffer FBO
  277 + gluUnProject(winX, winY, winZ, modelview, projection, viewport, &posX, &posY, &posZ); // not sure why it should add 1 to the winZ
  278 +
  279 + std::cout << "( " << posX << ", " << posY << ", "<< posZ <<" )" << std::endl;
  280 + }
  281 +}
  282 +
  283 +void glut_keyboard(unsigned char key, int x, int y){
  284 + if(key == 'm') // press m to change mode
  285 + {
  286 + if(compareMode && !mappingMode){ // if it is in compare mode
  287 + compareMode = false;
  288 + mappingMode = true;
  289 + }
  290 + else{ // if it is in mapping mode
  291 + compareMode = true;
  292 + mappingMode = false;
  293 + }
  294 + }
  295 + glutPostRedisplay();
145 296 }
146 297  
147 298 #define BREWER_CTRL_PTS 11 //number of control points in the Brewer map
... ... @@ -182,7 +333,7 @@ void glut_initialize(){
182 333 glutInit(&myargc, myargv); //pass bogus arguments to glutInit()
183 334 glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA); //generate a color buffer, depth buffer, and enable double buffering
184 335 glutInitWindowPosition(100,100); //set the initial window position
185   - glutInitWindowSize(320,320); //set the initial window size
  336 + glutInitWindowSize(320, 320); //set the initial window size
186 337 glutCreateWindow("NetMets - STIM Lab, UH"); //set the dialog box title
187 338  
188 339  
... ... @@ -190,16 +341,16 @@ void glut_initialize(){
190 341 glutDisplayFunc(glut_render); //function executed for rendering - renders networks
191 342 glutMouseFunc(glut_mouse); //executed on a mouse click - sets starting mouse positions for rotations
192 343 glutMotionFunc(glut_motion); //executed when the mouse is moved while a button is pressed
  344 + if(ind == 1) //only in mapping mode, keyboard will be used
  345 + glutKeyboardFunc(glut_keyboard);
193 346  
194   - texture_initialize(); //set up texture mapping (create texture maps, enable features)
  347 + texture_initialize(); //set up texture mapping (create texture maps, enable features)
195 348  
196 349 stim::vec3<float> c = bb.center(); //get the center of the network bounding box
197 350  
198 351 //place the camera along the z-axis at a distance determined by the network size along x and y
199 352 cam.setPosition(c + stim::vec<float>(0, 0, camera_factor * std::max(bb.size()[0], bb.size()[1])));
200 353 cam.LookAt(c[0], c[1], c[2]); //look at the center of the network
201   -
202   - glClearColor(1, 1, 1, 1);
203 354 }
204 355  
205 356 #ifdef __CUDACC__
... ... @@ -231,6 +382,27 @@ void compare(float sigma, int device){
231 382 std::cout << "FPR: " << FNR << std::endl;
232 383 }
233 384  
  385 +void map(float sigma, int device){
  386 +
  387 + _GT.split(GT, T, sigma, device);
  388 + _T.split(T, GT, sigma, device);
  389 +
  390 + _GT.mapping(_T, _gt_t, device);
  391 + _T.mapping(_GT, _t_gt, device);
  392 +
  393 + size_t num = _gt_t.size(); // also create random color for unmapping edge, but won't be used though
  394 + colormap.resize(3 * num); // 3 portions compound RGB
  395 + for(int i = 0; i < 3 * num; i++)
  396 + colormap[i] = rand()/(float)RAND_MAX;
  397 +
  398 + //calculate the metrics
  399 + float FPR = _GT.average(0); //calculate the metrics
  400 + float FNR = _T.average(0);
  401 +
  402 + std::cout << "FNR: " << FPR << std::endl; //print false alarms and misses
  403 + std::cout << "FPR: " << FNR << std::endl;
  404 +}
  405 +
234 406 // 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
235 407 // Pranathi wrote this - saves network features to a CSV file
236 408 void features(std::string filename){
... ... @@ -265,7 +437,7 @@ void advertise(){
265 437 std::cout<<"Thank you for using the NetMets network comparison tool!"<<std::endl;
266 438 std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
267 439 std::cout<<"Developers: Pranathi Vemuri, David Mayerich"<<std::endl;
268   - std::cout<<"Source: https://git.stim.ee.uh.edu/segmentation/netmets"<<std::endl;
  440 + std::cout<<"Source: https://git.stim.ee.uh.edu/segmentation/netmets" <<std::endl;
269 441 std::cout<<"========================================================================="<<std::endl<<std::endl;
270 442  
271 443 std::cout<<"usage: netmets file1 file2 --sigma 3"<<std::endl;
... ... @@ -286,6 +458,7 @@ int main(int argc, char* argv[])
286 458 args.add("gui", "display the network or network comparison using OpenGL");
287 459 args.add("device", "choose specific device to run", "0");
288 460 args.add("features", "save features to a CSV file, specify file name");
  461 + args.add("mapping", "mapping input according to similarity");
289 462  
290 463 args.parse(argc, argv); //parse the user arguments
291 464  
... ... @@ -310,14 +483,25 @@ int main(int argc, char* argv[])
310 483 features(args["features"].as_string());
311 484 GT = GT.resample(resample_rate * sigma); //resample both networks based on the sigma value
312 485 T = T.resample(resample_rate * sigma);
313   - setdevice(device);
314   - compare(sigma, device); //run the comparison algorithm
  486 + if(args["mapping"].is_set()){
  487 + map(sigma, device);
  488 + }
  489 + else
  490 + compare(sigma, device); //run the comparison algorithm
315 491 }
316 492  
317 493 //if a GUI is requested, display the network using OpenGL
318   - if(args["gui"].is_set()){
319   - bb = GT.boundingbox(); //generate a bounding volume
320   - glut_initialize(); //create the GLUT window and set callback functions
321   - glutMainLoop(); // enter GLUT event processing cycle
322   - }
  494 + if(args["gui"].is_set()){
  495 + if(args["mapping"].is_set()){
  496 + ind = 1;
  497 + bb = _GT.boundingbox(); //generate a bounding volume
  498 + glut_initialize(); //create the GLUT window and set callback functions
  499 + glutMainLoop(); // enter GLUT event processing cycle
  500 + }
  501 + else{
  502 + bb = GT.boundingbox(); //generate a bounding volume
  503 + glut_initialize(); //create the GLUT window and set callback functions
  504 + glutMainLoop(); // enter GLUT event processing cycle
  505 + }
  506 + }
323 507 }
... ...
main_dep.cu 0 → 100644
  Show/Hide comments   View file @d8e684d
  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 +#ifdef __CUDACC__
  18 +//CUDA includes
  19 +#include <cuda.h>
  20 +#endif
  21 +
  22 +//ANN includes
  23 +//#include <ANN/ANN.h>
  24 +
  25 +//BOOST includes
  26 +#include <boost/tuple/tuple.hpp>
  27 +
  28 +//visualization objects
  29 +stim::gl_aaboundingbox<float> bb; //axis-aligned bounding box object
  30 +stim::camera cam; //camera object
  31 +
  32 +unsigned num_nets = 0;
  33 +stim::gl_network<float> GT; //ground truth network
  34 +stim::gl_network<float> T; //test network
  35 +
  36 +//hard-coded parameters
  37 +float resample_rate = 0.5f; //sample rate for the network (fraction of sigma used as the maximum sample rate)
  38 +float camera_factor = 1.2f; //start point of the camera as a function of X and Y size
  39 +float orbit_factor = 0.01f; //degrees per pixel used to orbit the camera
  40 +
  41 +//mouse position tracking
  42 +int mouse_x;
  43 +int mouse_y;
  44 +
  45 +//OpenGL objects
  46 +GLuint cmap_tex = 0; //texture name for the color map
  47 +
  48 +//sets an OpenGL viewport taking up the entire window
  49 +void glut_render_single_projection(){
  50 +
  51 + glMatrixMode(GL_PROJECTION); //load the projection matrix for editing
  52 + glLoadIdentity(); //start with the identity matrix
  53 + int X = glutGet(GLUT_WINDOW_WIDTH); //use the whole screen for rendering
  54 + int Y = glutGet(GLUT_WINDOW_HEIGHT);
  55 + glViewport(0, 0, X, Y); //specify a viewport for the entire window
  56 + float aspect = (float)X / (float)Y; //calculate the aspect ratio
  57 + gluPerspective(60, aspect, 0.1, 1000000); //set up a perspective projection
  58 +}
  59 +
  60 +//sets an OpenGL viewport taking up the left half of the window
  61 +void glut_render_left_projection(){
  62 +
  63 + glMatrixMode(GL_PROJECTION); //load the projection matrix for editing
  64 + glLoadIdentity(); //start with the identity matrix
  65 + int X = glutGet(GLUT_WINDOW_WIDTH) / 2; //only use half of the screen for the viewport
  66 + int Y = glutGet(GLUT_WINDOW_HEIGHT);
  67 + glViewport(0, 0, X, Y); //specify the viewport on the left
  68 + float aspect = (float)X / (float)Y; //calculate the aspect ratio
  69 + gluPerspective(60, aspect, 0.1, 1000000); //set up a perspective projection
  70 +}
  71 +
  72 +//sets an OpenGL viewport taking up the right half of the window
  73 +void glut_render_right_projection(){
  74 +
  75 + glMatrixMode(GL_PROJECTION); //load the projection matrix for editing
  76 + glLoadIdentity(); //start with the identity matrix
  77 + int X = glutGet(GLUT_WINDOW_WIDTH) / 2; //only use half of the screen for the viewport
  78 + int Y = glutGet(GLUT_WINDOW_HEIGHT);
  79 + glViewport(X, 0, X, Y); //specify the viewport on the right
  80 + float aspect = (float)X / (float)Y; //calculate the aspect ratio
  81 + gluPerspective(60, aspect, 0.1, 1000000); //set up a perspective projection
  82 +}
  83 +
  84 +void glut_render_modelview(){
  85 +
  86 + glMatrixMode(GL_MODELVIEW); //load the modelview matrix for editing
  87 + glLoadIdentity(); //start with the identity matrix
  88 + stim::vec3<float> eye = cam.getPosition(); //get the camera position (eye point)
  89 + stim::vec3<float> focus = cam.getLookAt(); //get the camera focal point
  90 + stim::vec3<float> up = cam.getUp(); //get the camera "up" orientation
  91 +
  92 + gluLookAt(eye[0], eye[1], eye[2], focus[0], focus[1], focus[2], up[0], up[1], up[2]); //set up the OpenGL camera
  93 +}
  94 +
  95 +//draws the network(s)
  96 +void glut_render(void) {
  97 +
  98 + if(num_nets == 1){ //if a single network is loaded
  99 + glut_render_single_projection(); //fill the entire viewport
  100 + glut_render_modelview(); //set up the modelview matrix with camera details
  101 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  102 + GT.glCenterline(GT.nmags() - 1); //render the GT network (the only one loaded)
  103 + }
  104 +
  105 + if(num_nets == 2){ //if two networks are loaded
  106 +
  107 + glut_render_left_projection(); //set up a projection for the left half of the window
  108 + glut_render_modelview(); //set up the modelview matrix using camera details
  109 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  110 +
  111 + glEnable(GL_TEXTURE_1D); //enable texture mapping
  112 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //texture map will be used as the network color
  113 + glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the Brewer texture map
  114 +
  115 + GT.glCenterline(GT.nmags() - 1); //render the GT network
  116 +
  117 + glut_render_right_projection(); //set up a projection for the right half of the window
  118 + glut_render_modelview(); //set up the modelview matrix using camera details
  119 + T.glCenterline(T.nmags() - 1); //render the T network
  120 +
  121 + }
  122 +
  123 + glutSwapBuffers();
  124 +}
  125 +
  126 +// defines camera motion based on mouse dragging
  127 +void glut_motion(int x, int y){
  128 +
  129 +
  130 + float theta = orbit_factor * (mouse_x - x); //determine the number of degrees along the x-axis to rotate
  131 + float phi = orbit_factor * (y - mouse_y); //number of degrees along the y-axis to rotate
  132 +
  133 + cam.OrbitFocus(theta, phi); //rotate the camera around the focal point
  134 +
  135 + mouse_x = x; //update the mouse position
  136 + mouse_y = y;
  137 +
  138 + glutPostRedisplay(); //re-draw the visualization
  139 +}
  140 +
  141 +// sets the mouse position when clicked
  142 +void glut_mouse(int button, int state, int x, int y){
  143 + mouse_x = x;
  144 + mouse_y = y;
  145 +}
  146 +
  147 +#define BREWER_CTRL_PTS 11 //number of control points in the Brewer map
  148 +void texture_initialize(){
  149 +
  150 + //define the colormap
  151 + static float brewer_map[BREWER_CTRL_PTS][3] = { //generate a Brewer color map (blue to red)
  152 + {0.192157f, 0.211765f, 0.584314f},
  153 + {0.270588f, 0.458824f, 0.705882f},
  154 + {0.454902f, 0.678431f, 0.819608f},
  155 + {0.670588f, 0.85098f, 0.913725f},
  156 + {0.878431f, 0.952941f, 0.972549f},
  157 + {1.0f, 1.0f, 0.74902f},
  158 + {0.996078f, 0.878431f, 0.564706f},
  159 + {0.992157f, 0.682353f, 0.380392f},
  160 + {0.956863f, 0.427451f, 0.262745f},
  161 + {0.843137f, 0.188235f, 0.152941f},
  162 + {0.647059f, 0.0f, 0.14902f}
  163 + };
  164 +
  165 + glGenTextures(1, &cmap_tex); //generate a texture map name
  166 + glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the texture map
  167 +
  168 + glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_LINEAR); //enable linear interpolation
  169 + glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
  170 + glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP); //clamp the values at the minimum and maximum
  171 + glTexImage1D(GL_TEXTURE_1D, 0, 3, BREWER_CTRL_PTS, 0, GL_RGB, GL_FLOAT, //upload the texture map to the GPU
  172 + brewer_map);
  173 +}
  174 +
  175 +//Initialize the OpenGL (GLUT) window, including starting resolution, callbacks, texture maps, and camera
  176 +void glut_initialize(){
  177 +
  178 + int myargc = 1; //GLUT requires arguments, so create some bogus ones
  179 + char* myargv[1];
  180 + myargv [0]=strdup ("netmets");
  181 +
  182 + glutInit(&myargc, myargv); //pass bogus arguments to glutInit()
  183 + glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA); //generate a color buffer, depth buffer, and enable double buffering
  184 + glutInitWindowPosition(100,100); //set the initial window position
  185 + glutInitWindowSize(320,320); //set the initial window size
  186 + glutCreateWindow("NetMets - STIM Lab, UH"); //set the dialog box title
  187 +
  188 +
  189 + // register callback functions
  190 + glutDisplayFunc(glut_render); //function executed for rendering - renders networks
  191 + glutMouseFunc(glut_mouse); //executed on a mouse click - sets starting mouse positions for rotations
  192 + glutMotionFunc(glut_motion); //executed when the mouse is moved while a button is pressed
  193 +
  194 + texture_initialize(); //set up texture mapping (create texture maps, enable features)
  195 +
  196 + stim::vec3<float> c = bb.center(); //get the center of the network bounding box
  197 +
  198 + //place the camera along the z-axis at a distance determined by the network size along x and y
  199 + cam.setPosition(c + stim::vec<float>(0, 0, camera_factor * std::max(bb.size()[0], bb.size()[1])));
  200 + cam.LookAt(c[0], c[1], c[2]); //look at the center of the network
  201 +
  202 + glClearColor(1, 1, 1, 1);
  203 +}
  204 +
  205 +#ifdef __CUDACC__
  206 +void setdevice(int &device){
  207 + int count;
  208 + cudaGetDeviceCount(&count); // numbers of device that are available
  209 + if(count < device + 1){
  210 + std::cout<<"No such device available, please set another device"<<std::endl;
  211 + exit(1);
  212 + }
  213 +}
  214 +#else
  215 +void setdevice(int &device){
  216 + device = -1;
  217 +}
  218 +#endif
  219 +
  220 +//compare both networks and fill the networks with error information
  221 +void compare(float sigma, int device){
  222 +
  223 + GT = GT.compare(T, sigma, device); //compare the ground truth to the test case - store errors in GT
  224 + T = T.compare(GT, sigma, device); //compare the test case to the ground truth - store errors in T
  225 +
  226 + //calculate the metrics
  227 + float FPR = GT.average(0); //calculate the metrics
  228 + float FNR = T.average(0);
  229 +
  230 + std::cout << "FNR: " << FPR << std::endl; //print false alarms and misses
  231 + std::cout << "FPR: " << FNR << std::endl;
  232 +}
  233 +
  234 +// 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
  235 +// Pranathi wrote this - saves network features to a CSV file
  236 +void features(std::string filename){
  237 + double avgL_t, avgL_gt, avgT_t, avgT_gt, avgB_t, avgB_gt, avgC_t, avgC_gt, avgFD_t, avgFD_gt;
  238 + unsigned int e_t, e_gt, b_gt, b_t;
  239 + avgL_gt = GT.Lengths();
  240 + avgT_gt = GT.Tortuosities();
  241 + avgL_t = T.Lengths();
  242 + avgT_t = T.Tortuosities();
  243 + avgB_gt = GT.BranchingIndex();
  244 + avgB_t = T.BranchingIndex();
  245 + avgC_gt = GT.Contractions();
  246 + avgFD_gt = GT.FractalDimensions();
  247 + avgC_t = T.Contractions();
  248 + avgFD_t = T.FractalDimensions();
  249 + e_gt = GT.EndP();
  250 + e_t = T.EndP();
  251 + b_gt = GT.BranchP();
  252 + b_t = T.BranchP();
  253 + std::ofstream myfile;
  254 + myfile.open (filename.c_str());
  255 + myfile << "Length, Tortuosity, Contraction, Fractal Dimension, Branch Points, End points, Branching Index, \n";
  256 + myfile << avgL_gt << "," << avgT_gt << "," << avgC_gt << "," << avgFD_gt << "," << b_gt << "," << e_gt << "," << avgB_gt <<std::endl;
  257 + myfile << avgL_t << "," << avgT_t << "," << avgC_t << "," << avgFD_t << "," << b_t << "," << e_t << "," << avgB_t <<std::endl;
  258 + myfile.close();
  259 +}
  260 +
  261 +// Output an advertisement for the lab, authors, and usage information
  262 +void advertise(){
  263 + std::cout<<std::endl<<std::endl;
  264 + std::cout<<"========================================================================="<<std::endl;
  265 + std::cout<<"Thank you for using the NetMets network comparison tool!"<<std::endl;
  266 + std::cout<<"Scalable Tissue Imaging and Modeling (STIM) Lab, University of Houston"<<std::endl;
  267 + std::cout<<"Developers: Pranathi Vemuri, David Mayerich"<<std::endl;
  268 + std::cout<<"Source: https://git.stim.ee.uh.edu/segmentation/netmets"<<std::endl;
  269 + std::cout<<"========================================================================="<<std::endl<<std::endl;
  270 +
  271 + std::cout<<"usage: netmets file1 file2 --sigma 3"<<std::endl;
  272 + std::cout<<" compare two files with a tolerance of 3 (units defined by the network)"<<std::endl<<std::endl;
  273 + std::cout<<" netmets file1 --gui"<<std::endl;
  274 + std::cout<<" load a file and display it using OpenGL"<<std::endl<<std::endl;
  275 + std::cout<<" netmets file1 file2 --device 0"<<std::endl;
  276 + std::cout<<" compare two files using device 0 (if there isn't a gpu, use cpu)"<<std::endl<<std::endl;
  277 +}
  278 +
  279 +int main(int argc, char* argv[])
  280 +{
  281 + stim::arglist args; //create an instance of arglist
  282 +
  283 + //add arguments
  284 + args.add("help", "prints this help");
  285 + args.add("sigma", "force a sigma value to specify the tolerance of the network comparison", "3");
  286 + args.add("gui", "display the network or network comparison using OpenGL");
  287 + args.add("device", "choose specific device to run", "0");
  288 + args.add("features", "save features to a CSV file, specify file name");
  289 +
  290 + args.parse(argc, argv); //parse the user arguments
  291 +
  292 + if(args["help"].is_set() || args.nargs() == 0){ //test for help
  293 + advertise(); //output the advertisement
  294 + std::cout<<args.str(); //output arguments
  295 + exit(1); //exit
  296 + }
  297 +
  298 + if(args.nargs() >= 1){ //if at least one network file is specified
  299 + num_nets = 1; //set the number of networks to one
  300 + GT.load_obj(args.arg(0)); //load the specified file as the ground truth
  301 + /*GT.to_txt("Graph.txt");*/
  302 + }
  303 +
  304 + if(args.nargs() == 2){ //if two files are specified, they will be displayed in neighboring viewports and compared
  305 + int device = args["device"].as_int(); //get the device value from the user
  306 + num_nets = 2; //set the number of networks to two
  307 + float sigma = args["sigma"].as_float(); //get the sigma value from the user
  308 + T.load_obj(args.arg(1)); //load the second (test) network
  309 + if(args["features"].is_set()) //if the user wants to save features
  310 + features(args["features"].as_string());
  311 + GT = GT.resample(resample_rate * sigma); //resample both networks based on the sigma value
  312 + T = T.resample(resample_rate * sigma);
  313 + setdevice(device);
  314 + compare(sigma, device); //run the comparison algorithm
  315 + }
  316 +
  317 + //if a GUI is requested, display the network using OpenGL
  318 + if(args["gui"].is_set()){
  319 + bb = GT.boundingbox(); //generate a bounding volume
  320 + glut_initialize(); //create the GLUT window and set callback functions
  321 + glutMainLoop(); // enter GLUT event processing cycle
  322 + }
  323 +}
... ...