segmentation / netmets

Commit 8d519daa169a5f29537f18faff72190beb28d73f

Authored by David Mayerich
2 parents 484e01e8 f001495e

Merge branch 'JACK' into 'master' 

fix minor errors for loading and rendering networks from swc files

See merge request !7
Showing 1 changed file with 182 additions and 56 deletions   Show diff stats
Inline Side-by-side
1   -#include <stdlib.h>
  1 +#ifndef GL_MULTISAMPLE
  2 +#define GL_MULTISAMPLE 0x809D
  3 +#endif
  4 +
  5 +#include <stdlib.h>
2 6 #include <string>
3 7 #include <fstream>
4   -#include <algorithm>
  8 +#include <algorithm>
5 9  
6 10 //OpenGL includes
7 11 #include <GL/glut.h>
... ... @@ -19,12 +23,6 @@
19 23 #include <cuda.h>
20 24 #endif
21 25  
22   -//ANN includes
23   -//#include <ANN/ANN.h>
24   -
25   -//chrono includes
26   -//#include <chrono>
27   -
28 26 //BOOST includes
29 27 #include <boost/tuple/tuple.hpp>
30 28  
... ... @@ -32,15 +30,20 @@
32 30 stim::gl_aaboundingbox<float> bb; //axis-aligned bounding box object
33 31 stim::camera cam; //camera object
34 32  
  33 +// number of networks
35 34 unsigned num_nets = 0;
  35 +
  36 +// networks
36 37 stim::gl_network<float> GT; //ground truth network
37 38 stim::gl_network<float> T; //test network
38 39 stim::gl_network<float> _GT; //splitted GT
39 40 stim::gl_network<float> _T; //splitted T
40 41  
  42 +// indicator
41 43 unsigned ind = 0; //indicator of mapping
42 44 unsigned swc_ind = 0; //indicator of rendering swc file as networks
43 45  
  46 +// relationships
44 47 std::vector<unsigned> _gt_t; // store indices of nearest edge points in _T for _GT
45 48 std::vector<unsigned> _t_gt; // store indices of nearest edge points in _GT for _T
46 49  
... ... @@ -57,6 +60,10 @@ bool RButtonDown = false;
57 60 int mouse_x;
58 61 int mouse_y;
59 62  
  63 +// mouse wheel move
  64 +float des = 0.0f;
  65 +
  66 +// render modes
60 67 bool compareMode = true; // default mode is compare mode
61 68 bool mappingMode = false;
62 69  
... ... @@ -117,6 +124,7 @@ void glut_render_modelview(){
117 124 stim::vec3<float> focus = cam.getLookAt(); //get the camera focal point
118 125 stim::vec3<float> up = cam.getUp(); //get the camera "up" orientation
119 126  
  127 + eye[2] += des; //get camera closer to target by factor des
120 128 gluLookAt(eye[0], eye[1], eye[2], focus[0], focus[1], focus[2], up[0], up[1], up[2]); //set up the OpenGL camera
121 129 }
122 130  
... ... @@ -125,9 +133,10 @@ void glut_render(void) {
125 133  
126 134 //no mapping, just comparing
127 135 if (ind == 0) {
128   - //working with obj files
  136 + //-----working with obj files-----
129 137 if (swc_ind == 0) {
130 138 if (num_nets == 1) { //if a single network is loaded
  139 + glEnable(GL_DEPTH_TEST); //enable depth
131 140 glut_render_single_projection(); //fill the entire viewport
132 141 glut_render_modelview(); //set up the modelview matrix with camera details
133 142 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
... ... @@ -135,7 +144,7 @@ void glut_render(void) {
135 144 }
136 145  
137 146 if (num_nets == 2) { //if two networks are loaded
138   -
  147 + glEnable(GL_DEPTH_TEST); //enable depth
139 148 glut_render_left_projection(); //set up a projection for the left half of the window
140 149 glut_render_modelview(); //set up the modelview matrix using camera details
141 150 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
... ... @@ -144,16 +153,17 @@ void glut_render(void) {
144 153 glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //texture map will be used as the network color
145 154 glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the Brewer texture map
146 155  
147   - GT.glCenterline(); //render the GT network
  156 + GT.glCylinder(); //render the GT network
148 157  
149 158 glut_render_right_projection(); //set up a projection for the right half of the window
150 159 glut_render_modelview(); //set up the modelview matrix using camera details
151   - T.glCenterline(); //render the T network
  160 + T.glCylinder(); //render the T network
152 161 }
153 162 }
154   - //working with swc files
  163 + //-----working with swc files-----
155 164 else {
156 165 if (num_nets == 1) { //if a single network is loaded
  166 + glEnable(GL_DEPTH_TEST); //enable depth
157 167 glut_render_single_projection(); //fill the entire viewport
158 168 glut_render_modelview(); //set up the modelview matrix with camera details
159 169 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
... ... @@ -161,7 +171,7 @@ void glut_render(void) {
161 171 }
162 172  
163 173 if (num_nets == 2) { //if two networks are loaded
164   -
  174 + glEnable(GL_DEPTH_TEST); //enable depth
165 175 glut_render_left_projection(); //set up a projection for the left half of the window
166 176 glut_render_modelview(); //set up the modelview matrix using camera details
167 177 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
... ... @@ -179,9 +189,9 @@ void glut_render(void) {
179 189 }
180 190 }
181 191  
182   - //do mapping
  192 + //do comparing and mapping
183 193 else {
184   - //working with obj files
  194 + //-----working with obj files-----
185 195 if (swc_ind == 0) {
186 196 if (num_nets == 1) { //if a single network is loaded
187 197 std::cout << "You should have at least two networks to do mapping." << std::endl; //exit program because there isn't enough network
... ... @@ -189,6 +199,48 @@ void glut_render(void) {
189 199 }
190 200 if (num_nets == 2) { //if two networks are loaded
191 201 if (compareMode) {
  202 + glEnable(GL_DEPTH_TEST); //enable depth
  203 + glut_render_left_projection(); //set up a projection for the left half of the window
  204 + glut_render_modelview(); //set up the modelview matrix using camera details
  205 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  206 +
  207 + glEnable(GL_TEXTURE_1D); //enable texture mapping
  208 + glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_REPLACE); //texture map will be used as the network color
  209 + glBindTexture(GL_TEXTURE_1D, cmap_tex); //bind the Brewer texture map
  210 +
  211 + _GT.glCylinder(); //render the GT network
  212 +
  213 + glut_render_right_projection(); //set up a projection for the right half of the window
  214 + glut_render_modelview(); //set up the modelview matrix using camera details
  215 + _T.glCylinder(); //render the T network
  216 +
  217 + }
  218 + else {
  219 + glEnable(GL_DEPTH_TEST); //enable depth
  220 + glut_render_left_projection(); //set up a projection for the left half of the window
  221 + glut_render_modelview(); //set up the modelview matrix using camera details
  222 + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
  223 +
  224 + //_GT.glRandColorCenterlineGT(dlist1, _gt_t, colormap);
  225 + _GT.glRandColorCylinder1(dlist1, _gt_t, colormap);
  226 +
  227 + glut_render_right_projection(); //set up a projection for the right half of the window
  228 + glut_render_modelview(); //set up the modelview matrix using camera details
  229 + //_T.glRandColorCenterlineT(dlist2, _t_gt, colormap);
  230 + _T.glRandColorCylinder2(dlist2, _t_gt, colormap);
  231 + }
  232 + }
  233 + }
  234 + //-----working with swc files-----
  235 + else {
  236 + if (num_nets == 1) { //if a single network is loaded
  237 + std::cout << "You should have at least two networks to do mapping." << std::endl; //exit program because there isn't enough network
  238 + exit(1);
  239 + }
  240 + if (num_nets == 2) { //if two networks are loaded
  241 + //in compare mode
  242 + if (compareMode) {
  243 + glEnable(GL_DEPTH_TEST); //enable depth
192 244 glut_render_left_projection(); //set up a projection for the left half of the window
193 245 glut_render_modelview(); //set up the modelview matrix using camera details
194 246 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
... ... @@ -204,37 +256,38 @@ void glut_render(void) {
204 256 _T.glCenterline(); //render the T network
205 257  
206 258 }
  259 + //in mapping mode
207 260 else {
208   - glEnable(GL_DEPTH_TEST);
  261 + glEnable(GL_DEPTH_TEST); //enable depth
209 262 glut_render_left_projection(); //set up a projection for the left half of the window
210 263 glut_render_modelview(); //set up the modelview matrix using camera details
211 264 glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); //clear the screen
212 265  
213 266 //_GT.glRandColorCenterlineGT(dlist1, _gt_t, colormap);
214   - _GT.glCylinderGT(dlist1, _gt_t, colormap);
  267 + _GT.glRandColorCylinder1_swc(dlist1, _gt_t, colormap);
215 268  
216 269 glut_render_right_projection(); //set up a projection for the right half of the window
217 270 glut_render_modelview(); //set up the modelview matrix using camera details
218 271 //_T.glRandColorCenterlineT(dlist2, _t_gt, colormap);
219   - _T.glCylinderT(dlist2, _t_gt, colormap);
  272 + _T.glRandColorCylinder2_swc(dlist2, _t_gt, colormap);
220 273 }
221 274 }
222 275 }
223 276 }
224 277  
225   - if (num_nets == 2) {
  278 + if (num_nets == 2) { // works only with two networks
226 279 std::ostringstream ss;
227 280 if (mappingMode) // if it is in mapping mode
228 281 ss << "Mapping Mode";
229 282 else
230   - ss << "Compare Mode"; // default mode is compare mode
  283 + ss << "Compare Mode"; // default mode is compare mode
231 284  
232 285 glDisable(GL_TEXTURE_1D);
233   - glMatrixMode(GL_PROJECTION); //Set up the 2d viewport for mode text printing
  286 + glMatrixMode(GL_PROJECTION); // set up the 2d viewport for mode text printing
234 287 glPushMatrix();
235 288 glLoadIdentity();
236   - int X = glutGet(GLUT_WINDOW_WIDTH);
237   - int Y = glutGet(GLUT_WINDOW_HEIGHT);
  289 + int X = glutGet(GLUT_WINDOW_WIDTH); // get the current window width
  290 + int Y = glutGet(GLUT_WINDOW_HEIGHT); // get the current window height
238 291 glViewport(0, 0, X / 2, Y); // locate to left bottom corner
239 292 gluOrtho2D(0, X, 0, Y); // define othogonal aspect
240 293 glColor3f(0.0, 1.0, 0.0); // using green to show mode
... ... @@ -243,7 +296,7 @@ void glut_render(void) {
243 296 glPushMatrix();
244 297 glLoadIdentity();
245 298  
246   - glRasterPos2f(0, 5); //print text in the bottom left corner
  299 + glRasterPos2f(0, 5); //print text in the left bottom corner
247 300 glutBitmapString(GLUT_BITMAP_TIMES_ROMAN_24, (const unsigned char*)(ss.str().c_str()));
248 301  
249 302 glPopMatrix();
... ... @@ -270,6 +323,22 @@ void glut_motion(int x, int y){
270 323 }
271 324 }
272 325  
  326 +// sets the menu options
  327 +void glut_menu(int value) {
  328 + if (value == 1) { // menu 1 represents comparing mode
  329 + compareMode = true;
  330 + mappingMode = false;
  331 + }
  332 + if (value == 2) { // menu 2 represents mapping mode
  333 + compareMode = false;
  334 + mappingMode = true;
  335 + }
  336 + if (value == 3) {
  337 + exit(0);
  338 + }
  339 + glutPostRedisplay();
  340 +}
  341 +
273 342 // sets the mouse position when clicked
274 343 void glut_mouse(int button, int state, int x, int y){
275 344  
... ... @@ -319,8 +388,23 @@ void glut_mouse(int button, int state, int x, int y){
319 388 }
320 389 }
321 390  
  391 +// define camera move based on mouse wheel move(actually we can combine this with glut_mouse)
  392 +void glut_wheel(int wheel, int direction, int x, int y) {
  393 + float cam_move_fac; // camera move unit length
  394 + stim::vec3<float> eye = cam.getPosition(); // get the camera position (eye point)
  395 + stim::vec3<float> focus = cam.getLookAt(); // get the camera focal point
  396 + cam_move_fac = fabs(focus[2] - eye[2]) / 50; // divided by 50
  397 +
  398 + if (direction > 0) // if it is button 3(up)
  399 + des -= cam_move_fac;
  400 + else // if it is button 4(down)
  401 + des += cam_move_fac;
  402 + glutPostRedisplay();
  403 +}
  404 +
  405 +// define keyboard inputs
322 406 void glut_keyboard(unsigned char key, int x, int y){
323   - if(key == 'm') // press m to change mode
  407 + if(key == 'm') // if keyboard 'm' is pressed
324 408 {
325 409 if(compareMode && !mappingMode){ // if it is in compare mode
326 410 compareMode = false;
... ... @@ -331,6 +415,9 @@ void glut_keyboard(unsigned char key, int x, int y){
331 415 mappingMode = false;
332 416 }
333 417 }
  418 + if (key == 27) { // if keyboard "ESC" is pressed
  419 + exit(0); // exit
  420 + }
334 421 glutPostRedisplay();
335 422 }
336 423  
... ... @@ -370,18 +457,27 @@ void glut_initialize(){
370 457 myargv [0]=strdup ("netmets");
371 458  
372 459 glutInit(&myargc, myargv); //pass bogus arguments to glutInit()
373   - glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA); //generate a color buffer, depth buffer, and enable double buffering
  460 + glutSetOption(GLUT_MULTISAMPLE, 8);
  461 + glutInitDisplayMode(GLUT_DEPTH | GLUT_DOUBLE | GLUT_RGBA | GLUT_MULTISAMPLE); //generate a color buffer, depth buffer, and enable double buffering
374 462 glutInitWindowPosition(100,100); //set the initial window position
375 463 glutInitWindowSize(320, 320); //set the initial window size
376 464 glutCreateWindow("NetMets - STIM Lab, UH"); //set the dialog box title
377   -
378 465  
  466 + glEnable(GL_MULTISAMPLE);
379 467 // register callback functions
380 468 glutDisplayFunc(glut_render); //function executed for rendering - renders networks
381 469 glutMouseFunc(glut_mouse); //executed on a mouse click - sets starting mouse positions for rotations
382 470 glutMotionFunc(glut_motion); //executed when the mouse is moved while a button is pressed
383   - if(ind == 1) //only in mapping mode, keyboard will be used
384   - glutKeyboardFunc(glut_keyboard);
  471 + if (ind == 1) { //only in mapping mode, keyboard will be used
  472 + glutKeyboardFunc(glut_keyboard); //register keyboard callback
  473 + glutCreateMenu(glut_menu); //register menu option callback
  474 + glutAddMenuEntry("Comparing Mode", 1);//register menu 1 as comparing mode
  475 + glutAddMenuEntry("Mapping Mode", 2); //register menu 2 as mapping mode
  476 + glutAddMenuEntry("Exit", 3); //register menu 3 as exiting
  477 + glutAttachMenu(GLUT_RIGHT_BUTTON); //register right mouse to open menu option
  478 + }
  479 + if (swc_ind == 1) //only in rendering swc files, mouse wheel will be used
  480 + glutMouseWheelFunc(glut_wheel);
385 481  
386 482 texture_initialize(); //set up texture mapping (create texture maps, enable features)
387 483  
... ... @@ -393,6 +489,7 @@ void glut_initialize(){
393 489 }
394 490  
395 491 #ifdef __CUDACC__
  492 +// set specific device to work on
396 493 void setdevice(int &device){
397 494 int count;
398 495 cudaGetDeviceCount(&count); // numbers of device that are available
... ... @@ -403,7 +500,7 @@ void setdevice(int &amp;device){
403 500 }
404 501 #else
405 502 void setdevice(int &device){
406   - device = -1;
  503 + device = -1; // set to default -1
407 504 }
408 505 #endif
409 506  
... ... @@ -421,18 +518,22 @@ void compare(float sigma, int device){
421 518 std::cout << "FPR: " << FNR << std::endl;
422 519 }
423 520  
  521 +//split and map two networks and fill the networks' R with metric information
424 522 void map(float sigma, int device, float threshold){
425 523  
  524 + // compare and split two networks
426 525 _GT.split(GT, T, sigma, device, threshold);
427 526 _T.split(T, GT, sigma, device, threshold);
428 527  
  528 + // mapping two new splitted networks and get their edge relation
429 529 _GT.mapping(_T, _gt_t, device, threshold);
430 530 _T.mapping(_GT, _t_gt, device, threshold);
431 531  
  532 + // generate random color set based on the number of edges in GT
432 533 size_t num = _gt_t.size(); // also create random color for unmapping edge, but won't be used though
433 534 colormap.resize(3 * num); // 3 portions compound RGB
434 535 for(int i = 0; i < 3 * num; i++)
435   - colormap[i] = rand()/(float)RAND_MAX;
  536 + colormap[i] = rand()/(float)RAND_MAX; // set to [0, 1]
436 537  
437 538 //calculate the metrics
438 539 float FPR = _GT.average(0); //calculate the metrics
... ... @@ -442,6 +543,30 @@ void map(float sigma, int device, float threshold){
442 543 std::cout << "FPR: " << FNR << std::endl;
443 544 }
444 545  
  546 +void map_swc(float sigma, int device, float threshold) {
  547 +
  548 + // compare two networks
  549 + _GT = GT.compare(T, sigma, device); //compare the ground truth to the test case - store errors in _GT
  550 + _T = T.compare(_GT, sigma, device); //compare the test case to the ground truth - store errors in _T
  551 +
  552 + // mapping two networks and get their edge relation
  553 + _GT.mapping(_T, _gt_t, device, threshold);
  554 + _T.mapping(_GT, _t_gt, device, threshold);
  555 +
  556 + // generate random color set based on the number of edges in GT
  557 + size_t num = _gt_t.size(); // also create random color for unmapping edge, but won't be used though
  558 + colormap.resize(3 * num); // 3 portions compound RGB
  559 + for (int i = 0; i < 3 * num; i++)
  560 + colormap[i] = rand() / (float)RAND_MAX; // set to [0, 1]
  561 +
  562 + //calculate the metrics
  563 + float FPR = _GT.average(0); //calculate the metrics
  564 + float FNR = _T.average(0);
  565 +
  566 + std::cout << "FNR: " << FPR << std::endl; //print false alarms and misses
  567 + std::cout << "FPR: " << FNR << std::endl;
  568 +}
  569 +
445 570 // 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
446 571 // Pranathi wrote this - saves network features to a CSV file
447 572 void features(std::string filename){
... ... @@ -480,7 +605,7 @@ void advertise(){
480 605 std::cout<<"========================================================================="<<std::endl<<std::endl;
481 606  
482 607 std::cout<<"usage: netmets file1 file2 --sigma 3"<<std::endl;
483   - std::cout<<" compare two files with a tolerance of 3 (units defined by the network)"<<std::endl<<std::endl;
  608 + std::cout<<" compare two .obj files with a tolerance of 3 (units defined by the network)"<<std::endl<<std::endl;
484 609 std::cout<<" netmets file1 --gui"<<std::endl;
485 610 std::cout<<" load a file and display it using OpenGL"<<std::endl<<std::endl;
486 611 std::cout<<" netmets file1 file2 --device 0"<<std::endl;
... ... @@ -500,7 +625,6 @@ int main(int argc, char* argv[])
500 625 args.add("device", "choose specific device to run", "0");
501 626 args.add("features", "save features to a CSV file, specify file name");
502 627 args.add("mapping", "mapping input according to similarity");
503   - args.add("swc", "load swc file instead of obj file");
504 628  
505 629 args.parse(argc, argv); //parse the user arguments
506 630  
... ... @@ -510,15 +634,26 @@ int main(int argc, char* argv[])
510 634 exit(1); //exit
511 635 }
512 636  
513   - if (args["swc"].is_set()) { // if it is to load a swc file, right now we only load two files(.swc and .obj)
514   - if (args.nargs() >= 1) { // if at least one network file is specified
515   - num_nets = 1; // set the number of networks to one
516   - GT.load_swc(args.arg(0)); // load the specified file as the ground truth
517   - }
518   - if (args.nargs() == 2) { //if two files are specified, they will be displayed in neighboring viewports and compared
  637 + if (args.nargs() >= 1) { // if at least one network file is specified
  638 + num_nets = 1; // set the number of networks to one
  639 + std::vector<std::string> tmp = stim::parser::split(args.arg(0), '.'); // split the filename at '.'
  640 + if ("swc" == tmp[1]) { // loading swc file
  641 + GT.load_swc(args.arg(0)); // load the specified file as the ground truth
  642 + swc_ind = 1; // set the indicator of swc file to 1
  643 + }
  644 + else if ("obj" == tmp[1]) // loading obj file
  645 + GT.load_obj(args.arg(0)); // load the specified file as the ground truth
  646 + else {
  647 + std::cout << "Invalid loading file" << std::endl;
  648 + exit(1);
  649 + }
  650 + }
  651 +
  652 + if (args.nargs() == 2) { //if two files are specified, they will be displayed in neighboring viewports and compared
  653 + if (1 == swc_ind) { //loading swc files
519 654 int device = args["device"].as_int(); //get the device value from the user
520 655 num_nets = 2; //set the number of networks to two
521   - //does it need to be resampled??
  656 + //does it need to be resampled??
522 657 float sigma = args["sigma"].as_float(); //get the sigma value from the user
523 658 T.load_swc(args.arg(1)); //load the second (test) network
524 659 if (args["features"].is_set()) //if the user wants to save features
... ... @@ -526,21 +661,15 @@ int main(int argc, char* argv[])
526 661 GT = GT.resample(resample_rate * sigma); //resample both networks based on the sigma value
527 662 T = T.resample(resample_rate * sigma);
528 663 if (args["mapping"].is_set()) {
529   - std::cout << "right now networks that are loaded from swc files do not need to be mapped with each other!" << std::endl;
530   - exit(1);
  664 + float threshold = args["mapping"].as_float();
  665 + map_swc(sigma, device, threshold);
  666 + //std::cout << "right now networks that are loaded from swc files do not need to be mapped with each other!" << std::endl;
  667 + //exit(1);
531 668 }
532 669 else
533 670 compare(sigma, device); //run the comparison algorithm
534 671 }
535   - }
536   -
537   - else { //if it is to load a obj file
538   - if (args.nargs() >= 1) { //if at least one network file is specified
539   - num_nets = 1; //set the number of networks to one
540   - GT.load_obj(args.arg(0)); //load the specified file as the ground truth
541   - /*GT.to_txt("Graph.txt");*/
542   - }
543   - if (args.nargs() == 2) { //if two files are specified, they will be displayed in neighboring viewports and compared
  672 + else {
544 673 int device = args["device"].as_int(); //get the device value from the user
545 674 num_nets = 2; //set the number of networks to two
546 675 float sigma = args["sigma"].as_float(); //get the sigma value from the user
... ... @@ -554,15 +683,12 @@ int main(int argc, char* argv[])
554 683 map(sigma, device, threshold);
555 684 }
556 685 else
557   - compare(sigma, device); //run the comparison algorithm
  686 + compare(sigma, device); //run the comparison algorithm
558 687 }
559 688 }
560 689  
561 690 //if a GUI is requested, display the network using OpenGL
562 691 if(args["gui"].is_set()){
563   - if (args["swc"].is_set()) {
564   - swc_ind = 1;
565   - }
566 692 if (args["mapping"].is_set()) {
567 693 ind = 1; //set indicator of mapping to 1(true)
568 694 bb = _GT.boundingbox(); //generate a bounding volume
... ...