THESE/codes/graphe/GCmex1.9/GCoptimization.cpp

   1 #include "energy.h"\r
   2 #include "graph.h"\r
   3 #include "GCoptimization.h"\r
   4 #include <stdio.h>\r
   5 #include <time.h>\r
   6 #include <stdlib.h>\r
   7 #define MAX_INTT 1000000000\r
   8 \r
   9 \r
  10 /**************************************************************************************/\r
  11 \r
  12 void GCoptimization::initialize_memory()\r
  13 {\r
  14     int i = 0;\r
  15     \r
  16         m_lookupPixVar = (PixelType *) new PixelType[m_num_pixels];\r
  17         m_labelTable   = (LabelType *) new LabelType[m_num_labels];\r
  18 \r
  19         terminateOnError( !m_lookupPixVar || !m_labelTable,"Not enough memory");\r
  20 \r
  21         for ( i = 0; i < m_num_labels; i++ )\r
  22                 m_labelTable[i] = i;\r
  23 \r
  24         for ( i = 0; i < m_num_pixels; i++ )\r
  25                 m_lookupPixVar[i] = -1;\r
  26 \r
  27 }\r
  28 \r
  29 \r
  30 /**************************************************************************************/\r
  31 \r
  32 void GCoptimization::commonGridInitialization(PixelType width, PixelType height, int nLabels)\r
  33 {\r
  34         terminateOnError( (width < 0) || (height <0) || (nLabels <0 ),"Illegal negative parameters");\r
  35 \r
  36         m_width              = width;\r
  37         m_height             = height;\r
  38         m_num_pixels         = width*height;\r
  39         m_num_labels         = nLabels;\r
  40         m_grid_graph         = 1;\r
  41                 \r
  42         \r
  43         srand(time(NULL));\r
  44 }\r
  45 /**************************************************************************************/\r
  46 \r
  47 void GCoptimization::commonNonGridInitialization(PixelType nupixels, int num_labels)\r
  48 {\r
  49         terminateOnError( (nupixels <0) || (num_labels <0 ),"Illegal negative parameters");\r
  50 \r
  51         m_num_labels             = num_labels;\r
  52         m_num_pixels             = nupixels;\r
  53         m_grid_graph         = 0;\r
  54 \r
  55         m_neighbors = (LinkedBlockList *) new LinkedBlockList[nupixels];\r
  56 \r
  57         terminateOnError(!m_neighbors,"Not enough memory");\r
  58 \r
  59 }\r
  60 \r
  61 /**************************************************************************************/\r
  62 \r
  63 void GCoptimization::commonInitialization(int dataSetup, int smoothSetup)\r
  64 {\r
  65         int i;\r
  66 \r
  67         m_random_label_order = 1;\r
  68         m_dataInput          = dataSetup;\r
  69         m_smoothInput        = smoothSetup;\r
  70 \r
  71 \r
  72         if (m_dataInput == SET_INDIVIDUALLY )\r
  73         {\r
  74                 m_datacost    = (EnergyTermType *) new EnergyTermType[m_num_labels*m_num_pixels];\r
  75                 terminateOnError(!m_datacost,"Not enough memory");\r
  76                 for ( i = 0; i < m_num_labels*m_num_pixels; i++ ) \r
  77                         m_datacost[i] = (EnergyTermType) 0;\r
  78                 m_dataType = ARRAY;\r
  79         }\r
  80         else m_dataType = NONE;\r
  81 \r
  82 \r
  83         if ( m_smoothInput == SET_INDIVIDUALLY )\r
  84         {\r
  85                 m_smoothcost  = (EnergyTermType *) new EnergyTermType[m_num_labels*m_num_labels];\r
  86                 terminateOnError(!m_smoothcost,"Not enough memory");\r
  87                 for ( i = 0; i < m_num_labels*m_num_labels; i++ ) \r
  88                         m_smoothcost[i] = (EnergyTermType) 0;\r
  89 \r
  90                 m_smoothType      = ARRAY;\r
  91                 m_varying_weights = 0;\r
  92         }\r
  93         else m_smoothType  = NONE;\r
  94 \r
  95 \r
  96         initialize_memory();\r
  97         srand(time(NULL));\r
  98     // <!-- bagon\r
  99     // sign class as valid\r
 100     class_sig = VALID_CLASS_SIGNITURE;\r
 101     truncate_flag = false;\r
 102     // bagon -->\r
 103 }\r
 104 \r
 105 /**************************************************************************************/\r
 106 /* Use this constructor only for grid graphs                                          */\r
 107 GCoptimization::GCoptimization(PixelType width,PixelType height,int nLabels,int dataSetup, int smoothSetup )\r
 108 {\r
 109         commonGridInitialization(width,height,nLabels);\r
 110         \r
 111         m_labeling           = (LabelType *) new LabelType[m_num_pixels];\r
 112         terminateOnError(!m_labeling,"out of memory");\r
 113         for ( int i = 0; i < m_num_pixels; i++ ) m_labeling[i] = (LabelType) 0;\r
 114 \r
 115         m_deleteLabeling = 1;\r
 116         commonInitialization(dataSetup,smoothSetup);    \r
 117 }\r
 118 \r
 119 /**************************************************************************************/\r
 120 /* Use this constructor only for grid graphs                                          */\r
 121 GCoptimization::GCoptimization(LabelType *m_answer,PixelType width,PixelType height,int nLabels,\r
 122                                                            int dataSetup, int smoothSetup) \r
 123                                                    \r
 124 {\r
 125         commonGridInitialization(width,height,nLabels);\r
 126 \r
 127         m_labeling = m_answer;\r
 128 \r
 129         \r
 130         for ( int i = 0; i < m_num_pixels; i++ )\r
 131                 terminateOnError(m_labeling[i] < 0 || m_labeling[i] >= nLabels,"Initial labels are out of valid range");\r
 132         \r
 133         m_deleteLabeling = 0;\r
 134                 \r
 135         commonInitialization(dataSetup,smoothSetup);    \r
 136 }\r
 137 \r
 138 /**************************************************************************************/\r
 139 /* Use this constructor for general graphs                                            */\r
 140 GCoptimization::GCoptimization(PixelType nupixels,int nLabels,int dataSetup, int smoothSetup )\r
 141 {\r
 142 \r
 143         commonNonGridInitialization(nupixels, nLabels);\r
 144         \r
 145         m_labeling           = (LabelType *) new LabelType[m_num_pixels];\r
 146         terminateOnError(!m_labeling,"out of memory");\r
 147         for ( int i = 0; i < nupixels; i++ ) m_labeling[i] = (LabelType) 0;\r
 148 \r
 149         m_deleteLabeling = 1;\r
 150 \r
 151         commonInitialization(dataSetup,smoothSetup);    \r
 152 }\r
 153 \r
 154 /**************************************************************************************/\r
 155 /* Use this constructor for general graphs                                            */\r
 156 GCoptimization::GCoptimization(LabelType *m_answer, PixelType nupixels,int nLabels,int dataSetup, int smoothSetup)\r
 157 {\r
 158         commonNonGridInitialization(nupixels, nLabels);\r
 159         \r
 160 \r
 161         m_labeling = m_answer;\r
 162         for ( int i = 0; i < m_num_pixels; i++ )\r
 163                 terminateOnError(m_labeling[i] < 0 || m_labeling[i] >= nLabels,"Initial labels are out of valid range");\r
 164 \r
 165         m_deleteLabeling = 0;\r
 166 \r
 167         commonInitialization(dataSetup,smoothSetup);    \r
 168 }\r
 169 \r
 170 /**************************************************************************************/\r
 171 \r
 172 void GCoptimization::setData(EnergyTermType* dataArray)\r
 173 {\r
 174         terminateOnError(m_dataType != NONE,\r
 175                             "ERROR: you already set the data, or said you'll use member function setDataCost() to set data");   \r
 176 \r
 177     // <!-- bagon\r
 178     /* allocate memory dinamiocally */\r
 179     m_datacost = (EnergyTermType *) new EnergyTermType[m_num_pixels*m_num_labels];\r
 180     for ( int i(0); i < m_num_pixels*m_num_labels; i++ )\r
 181         m_datacost[i] = dataArray[i];\r
 182         // m_datacost  = dataArray;\r
 183     // bagon -->\r
 184         m_dataType  = ARRAY;\r
 185 }\r
 186 \r
 187 /**************************************************************************************/\r
 188 \r
 189 void GCoptimization::setData(dataFnPix dataFn)\r
 190 {\r
 191         terminateOnError(m_dataType != NONE,\r
 192                             "ERROR: you already set the data, or said you'll use member function setDataCost() to set data");   \r
 193         \r
 194         m_dataFnPix = dataFn;\r
 195         m_dataType  = FUNCTION_PIX;\r
 196 }\r
 197 \r
 198 /**************************************************************************************/\r
 199 \r
 200 void GCoptimization::setData(dataFnCoord dataFn)\r
 201 {\r
 202         terminateOnError(m_dataType != NONE,\r
 203                             "ERROR: you already set the data, or said you'll use member function setDataCost() to set data");   \r
 204 \r
 205         terminateOnError( !m_grid_graph,"Cannot use data function based on coordinates for non-grid graph");\r
 206 \r
 207         m_dataFnCoord = dataFn;\r
 208         m_dataType    = FUNCTION_COORD;\r
 209 }\r
 210 \r
 211 /**************************************************************************************/\r
 212 \r
 213 void GCoptimization::setSmoothness(EnergyTermType* V)\r
 214 {\r
 215 \r
 216         terminateOnError(m_smoothType != NONE,\r
 217                             "ERROR: you already set smoothness, or said you'll use member function setSmoothCost() to set Smoothness Costs");   \r
 218 \r
 219 \r
 220         m_smoothType = ARRAY;\r
 221     // <!-- bagon\r
 222     /* allocate memory dinamiocally */\r
 223     m_smoothcost = (EnergyTermType *) new EnergyTermType[m_num_labels*m_num_labels];\r
 224     for ( int i(0); i < m_num_labels*m_num_labels; i++ )\r
 225         m_smoothcost[i] = V[i];\r
 226         // m_smoothcost = V;\r
 227     // bagon -->\r
 228     m_varying_weights = 0;\r
 229 \r
 230 }\r
 231 \r
 232 /**************************************************************************************/\r
 233 void GCoptimization::setSmoothness(EnergyTermType* V,EnergyTermType* hCue, EnergyTermType* vCue)\r
 234 {\r
 235         terminateOnError(m_smoothType != NONE,\r
 236                             "ERROR: you already set smoothness, or said you'll use member function setSmoothCost() to set Smoothness Costs");   \r
 237 \r
 238         terminateOnError(!m_grid_graph,\r
 239                             "ERROR: for a grid graph, you can't use vertical and horizontal cues.  Use setNeighbors() member function to encode spatially varying cues");       \r
 240 \r
 241         m_varying_weights = 1;\r
 242 \r
 243     m_smoothType      = ARRAY;\r
 244     // <!-- bagon\r
 245     /* allocate memory dinamiocally */\r
 246     int i = 0;\r
 247     m_smoothcost = (EnergyTermType *) new EnergyTermType[m_num_labels*m_num_labels];\r
 248     for ( i = 0; i < m_num_labels*m_num_labels; i++ )\r
 249         m_smoothcost[i] = V[i];\r
 250     \r
 251     m_vertWeights = (EnergyTermType *) new EnergyTermType[m_num_pixels];\r
 252     for ( i = 0; i < m_num_pixels; i++ )\r
 253         m_vertWeights[i] = vCue[i];\r
 254     \r
 255     m_horizWeights = (EnergyTermType *) new EnergyTermType[m_num_pixels];\r
 256     for ( i = 0; i < m_num_pixels; i++ )\r
 257         m_horizWeights[i] = hCue[i];\r
 258     \r
 259         //m_vertWeights     = vCue;\r
 260         //m_horizWeights    = hCue;\r
 261         //m_smoothcost      = V;\r
 262     // bagon -->\r
 263 }\r
 264 \r
 265 /**************************************************************************************/\r
 266 \r
 267 void GCoptimization::setSmoothness(smoothFnCoord horz_cost, smoothFnCoord vert_cost)\r
 268 {\r
 269 \r
 270         terminateOnError(m_smoothType != NONE,\r
 271                             "ERROR: you already set smoothness, or said you'll use member function setSmoothCost() to set Smoothness Costs");   \r
 272 \r
 273         terminateOnError( !m_grid_graph,"Cannot use smoothness function based on coordinates for non-grid graph");\r
 274 \r
 275         m_smoothType    = FUNCTION_COORD;\r
 276         m_horz_cost     = horz_cost;\r
 277         m_vert_cost     = vert_cost;\r
 278 }\r
 279 \r
 280 /**************************************************************************************/\r
 281 \r
 282 void GCoptimization::setSmoothness(smoothFnPix cost)\r
 283 {\r
 284         terminateOnError(m_smoothType != NONE,\r
 285                             "ERROR: you already set smoothness, or said you'll use member function setSmoothCost() to set Smoothness Costs");   \r
 286 \r
 287         m_smoothType    = FUNCTION_PIX;\r
 288         m_smoothFnPix   = cost;\r
 289 }\r
 290 \r
 291 /**************************************************************************************/\r
 292 // <!-- bagon\r
 293 void \r
 294 GCoptimization::SetAllLabels(const LabelType* labels)\r
 295 {\r
 296     for( int i(0); i < m_num_pixels; i++ ) {\r
 297         terminateOnError( (labels[i]<0) || (labels[i]>=m_num_labels), "Wrong label value");\r
 298         m_labeling[i] = labels[i];\r
 299     }\r
 300 }\r
 301 /**************************************************************************************/    \r
 302 void \r
 303 GCoptimization::ExportLabels(LabelType* labels)\r
 304 {\r
 305     for( int i(0); i < m_num_pixels; i++)\r
 306         labels[i] = m_labeling[i];\r
 307 }\r
 308 // bagon -->\r
 309 /**************************************************************************************/    \r
 310 GCoptimization::EnergyType GCoptimization::giveDataEnergy()\r
 311 {\r
 312         if ( m_dataType == ARRAY) \r
 313                 return(giveDataEnergyArray());\r
 314         else if ( m_dataType == FUNCTION_PIX )\r
 315                 return(giveDataEnergyFnPix());\r
 316         else if (m_dataType == FUNCTION_COORD ) return(giveDataEnergyFnCoord());\r
 317         else terminateOnError(1,"Did not initialize the data costs yet");\r
 318 \r
 319         return(0);\r
 320 }\r
 321 \r
 322 /**************************************************************************************/\r
 323         \r
 324 GCoptimization::EnergyType GCoptimization::giveDataEnergyArray()\r
 325 {\r
 326         EnergyType eng = (EnergyType) 0;\r
 327 \r
 328 \r
 329         for ( int i = 0; i < m_num_pixels; i++ ) {\r
 330         \r
 331         if (m_labeling[i] != m_num_labels+1) {\r
 332             eng = eng + m_datacost(i,m_labeling[i]);\r
 333             // mexPrintf("node %d, label %d : %f \t (%f)\n", i, m_labeling[i], m_datacost(i,m_labeling[i]), eng);\r
 334         }\r
 335     }\r
 336 \r
 337         return(eng);\r
 338 }\r
 339 \r
 340 /**************************************************************************************/\r
 341         \r
 342 GCoptimization::EnergyType GCoptimization::giveDataEnergyFnPix()\r
 343 {\r
 344 \r
 345         EnergyType eng = (EnergyType) 0;\r
 346 \r
 347         for ( int i = 0; i < m_num_pixels; i++ )\r
 348                 eng = eng + m_dataFnPix(i,m_labeling[i]);\r
 349 \r
 350         return(eng);\r
 351 }\r
 352 /**************************************************************************************/\r
 353 \r
 354 GCoptimization::EnergyType GCoptimization::giveDataEnergyFnCoord()\r
 355 {\r
 356         EnergyType eng = (EnergyType) 0;\r
 357 \r
 358         for ( int y = 0; y < m_height; y++ )\r
 359                 for ( int x = 0; x < m_width; x++ )\r
 360                         eng = eng + m_dataFnCoord(x,y,m_labeling[x+y*m_width]);\r
 361 \r
 362         return(eng);\r
 363 \r
 364 }\r
 365 \r
 366 /**************************************************************************************/\r
 367 \r
 368 GCoptimization::EnergyType GCoptimization::giveSmoothEnergy()\r
 369 {\r
 370 \r
 371         if ( m_grid_graph )\r
 372         {\r
 373                 if ( m_smoothType == ARRAY )\r
 374                 {\r
 375                         if (m_varying_weights) return(giveSmoothEnergy_G_ARRAY_VW());\r
 376                         else return(giveSmoothEnergy_G_ARRAY());\r
 377                 }\r
 378                 else if ( m_smoothType == FUNCTION_PIX ) return(giveSmoothEnergy_G_FnPix());\r
 379                 else if ( m_smoothType == FUNCTION_COORD ) return(giveSmoothEnergy_G_FnCoord());\r
 380                 else terminateOnError(1,"Did not initialize smoothness costs yet, can't compute smooth energy");\r
 381         }\r
 382         else\r
 383         {\r
 384                 if ( m_smoothType == ARRAY ) return(giveSmoothEnergy_NG_ARRAY());\r
 385                 else if ( m_smoothType == FUNCTION_PIX ) return(giveSmoothEnergy_NG_FnPix());\r
 386                 else terminateOnError(1,"Did not initialize smoothness costs yet, can't compute smooth energy");\r
 387         }\r
 388         return(0);\r
 389 \r
 390 }\r
 391 \r
 392 /**************************************************************************************/\r
 393 \r
 394 GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_NG_FnPix()\r
 395 {\r
 396 \r
 397         EnergyType eng = (EnergyType) 0;\r
 398         int i;\r
 399         Neighbor *temp; \r
 400 \r
 401         for ( i = 0; i < m_num_pixels; i++ )\r
 402                 if ( !m_neighbors[i].isEmpty() )\r
 403                 {\r
 404                         m_neighbors[i].setCursorFront();\r
 405                         while ( m_neighbors[i].hasNext() )\r
 406                         {\r
 407                                 temp = (Neighbor *) m_neighbors[i].next();\r
 408                                 if ( i < temp->to_node )\r
 409                                         eng = eng + m_smoothFnPix(i,temp->to_node, m_labeling[i],m_labeling[temp->to_node]);\r
 410                         }\r
 411                 }\r
 412                 \r
 413         return(eng);\r
 414         \r
 415 }\r
 416 \r
 417 /**************************************************************************************/\r
 418 \r
 419 GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_NG_ARRAY()\r
 420 {\r
 421         EnergyType eng = (EnergyType) 0;\r
 422         int i;\r
 423         Neighbor *temp; \r
 424 \r
 425         for ( i = 0; i < m_num_pixels; i++ )\r
 426                 if ( !m_neighbors[i].isEmpty() )\r
 427                 {\r
 428                         m_neighbors[i].setCursorFront();\r
 429                         while ( m_neighbors[i].hasNext() )\r
 430                         {\r
 431                                 temp = (Neighbor *) m_neighbors[i].next();\r
 432 \r
 433                                 if ( i < temp->to_node )\r
 434                                         eng = eng + m_smoothcost(m_labeling[i],m_labeling[temp->to_node])*(temp->weight);\r
 435 \r
 436                         }\r
 437                 }\r
 438 \r
 439         return(eng);\r
 440 }\r
 441 \r
 442 /**************************************************************************************/\r
 443 \r
 444 GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_G_ARRAY_VW()\r
 445 {\r
 446 \r
 447         EnergyType eng = (EnergyType) 0;\r
 448         int x,y,pix;\r
 449 \r
 450         // printf("\nIn right place");\r
 451         \r
 452         for ( y = 0; y < m_height; y++ )\r
 453                 for ( x = 1; x < m_width; x++ ) if ((m_labeling[x+y*m_width-1] != m_num_labels +1) &&(m_labeling[x+y*m_width] != m_num_labels +1) ) \r
 454                 {\r
 455                         pix = x+y*m_width;\r
 456                         eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-1])*m_horizWeights[pix-1];\r
 457                 }\r
 458 \r
 459         for ( y = 1; y < m_height; y++ )\r
 460                 for ( x = 0; x < m_width; x++ )if ((m_labeling[x+y*m_width-m_width] != m_num_labels +1) &&(m_labeling[x+y*m_width] != m_num_labels +1) ) \r
 461                 {\r
 462                         pix = x+y*m_width;\r
 463                         eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-m_width])*m_vertWeights[pix-m_width];\r
 464                 }\r
 465 \r
 466         \r
 467         return(eng);\r
 468         \r
 469 }\r
 470 /**************************************************************************************/\r
 471 \r
 472 GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_G_ARRAY()\r
 473 {\r
 474 \r
 475         EnergyType eng = (EnergyType) 0;\r
 476         int x,y,pix;\r
 477 \r
 478 \r
 479         for ( y = 0; y < m_height; y++ )\r
 480                 for ( x = 1; x < m_width; x++ )\r
 481                 {\r
 482                         pix = x+y*m_width;\r
 483                         eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-1]);\r
 484                 }\r
 485 \r
 486         for ( y = 1; y < m_height; y++ )\r
 487                 for ( x = 0; x < m_width; x++ )\r
 488                 {\r
 489                         pix = x+y*m_width;\r
 490                         eng = eng + m_smoothcost(m_labeling[pix],m_labeling[pix-m_width]);\r
 491                 }\r
 492 \r
 493         return(eng);\r
 494         \r
 495 }\r
 496 /**************************************************************************************/\r
 497 \r
 498 GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_G_FnPix()\r
 499 {\r
 500 \r
 501         EnergyType eng = (EnergyType) 0;\r
 502         int x,y,pix;\r
 503 \r
 504 \r
 505         for ( y = 0; y < m_height; y++ )\r
 506                 for ( x = 1; x < m_width; x++ )\r
 507                 {\r
 508                         pix = x+y*m_width;\r
 509                         eng = eng + m_smoothFnPix(pix,pix-1,m_labeling[pix],m_labeling[pix-1]);\r
 510                 }\r
 511 \r
 512         for ( y = 1; y < m_height; y++ )\r
 513                 for ( x = 0; x < m_width; x++ )\r
 514                 {\r
 515                         pix = x+y*m_width;\r
 516                         eng = eng + m_smoothFnPix(pix,pix-m_width,m_labeling[pix],m_labeling[pix-m_width]);\r
 517                 }\r
 518 \r
 519         return(eng);\r
 520         \r
 521 }\r
 522 /**************************************************************************************/\r
 523 \r
 524 GCoptimization::EnergyType GCoptimization::giveSmoothEnergy_G_FnCoord()\r
 525 {\r
 526 \r
 527         EnergyType eng = (EnergyType) 0;\r
 528         int x,y,pix;\r
 529 \r
 530 \r
 531         for ( y = 0; y < m_height; y++ )\r
 532                 for ( x = 1; x < m_width; x++ )\r
 533                 {\r
 534                         pix = x+y*m_width;\r
 535                         eng = eng + m_horz_cost(x-1,y,m_labeling[pix],m_labeling[pix-1]);\r
 536                 }\r
 537 \r
 538         for ( y = 1; y < m_height; y++ )\r
 539                 for ( x = 0; x < m_width; x++ )\r
 540                 {\r
 541                         pix = x+y*m_width;\r
 542                         eng = eng + m_vert_cost(x,y-1,m_labeling[pix],m_labeling[pix-m_width]);\r
 543                 }\r
 544 \r
 545         return(eng);\r
 546         \r
 547 }\r
 548 \r
 549 /**************************************************************************************/\r
 550  \r
 551 GCoptimization::EnergyType GCoptimization::compute_energy()\r
 552 {\r
 553         return(giveDataEnergy()+giveSmoothEnergy());\r
 554 }\r
 555 \r
 556 /**************************************************************************************/\r
 557 \r
 558 GCoptimization::EnergyType GCoptimization::expansion(int max_num_iterations)\r
 559 {\r
 560         return(start_expansion(max_num_iterations)); \r
 561 }\r
 562 \r
 563 /**************************************************************************************/\r
 564 \r
 565 \r
 566 GCoptimization::EnergyType GCoptimization::expansion()\r
 567 {\r
 568         return(start_expansion(MAX_INTT));\r
 569 }\r
 570 \r
 571 /**************************************************************************************/\r
 572 \r
 573 \r
 574 GCoptimization::EnergyType GCoptimization::start_expansion(int max_num_iterations )\r
 575 {\r
 576         \r
 577         int curr_cycle = 1;\r
 578         EnergyType new_energy,old_energy;\r
 579         \r
 580 \r
 581         new_energy = compute_energy();\r
 582 \r
 583         old_energy = (new_energy+1)*10; // BAGON changed init value to exceed current energy by factor of 10 (thanks to A. Khan)\r
 584 \r
 585         while ( old_energy > new_energy  && curr_cycle <= max_num_iterations)\r
 586         {\r
 587                 old_energy = new_energy;\r
 588                 new_energy = oneExpansionIteration();\r
 589                 \r
 590                 curr_cycle++;   \r
 591         }\r
 592 \r
 593         return(new_energy);\r
 594 }\r
 595 \r
 596 /****************************************************************************/\r
 597 \r
 598 void GCoptimization::scramble_label_table()\r
 599 {\r
 600    LabelType r1,r2,temp;\r
 601    int num_times,cnt;\r
 602 \r
 603 \r
 604    num_times = m_num_labels*2;\r
 605 \r
 606    for ( cnt = 0; cnt < num_times; cnt++ )\r
 607    {\r
 608       r1 = rand()%m_num_labels;  \r
 609       r2 = rand()%m_num_labels;  \r
 610 \r
 611       temp             = m_labelTable[r1];\r
 612       m_labelTable[r1] = m_labelTable[r2];\r
 613       m_labelTable[r2] = temp;\r
 614    }\r
 615 }\r
 616 \r
 617 /**************************************************************************************/\r
 618 \r
 619 GCoptimization::EnergyType GCoptimization::alpha_expansion(LabelType label)\r
 620 {\r
 621         terminateOnError( label < 0 || label >= m_num_labels,"Illegal Label to Expand On");\r
 622 \r
 623         perform_alpha_expansion(label);\r
 624         return(compute_energy());\r
 625 }\r
 626 \r
 627 /**************************************************************************************/\r
 628 \r
 629 GCoptimization::EnergyType GCoptimization::alpha_expansion(LabelType alpha_label, PixelType *pixels, int num )\r
 630 {\r
 631         PixelType i,size = 0; \r
 632         Energy *e = new Energy(mexErrMsgTxt, truncate_flag);\r
 633         \r
 634 \r
 635 \r
 636         for ( i = 0; i<num ; i++ )\r
 637         {\r
 638                 if ( m_labeling[pixels[i]] != alpha_label )\r
 639                 {\r
 640                         m_lookupPixVar[pixels[i]] = i;\r
 641                 }\r
 642         }\r
 643 \r
 644         \r
 645         if ( size > 0 ) \r
 646         {\r
 647                 Energy::Var *variables = (Energy::Var *) new Energy::Var[size];\r
 648 \r
 649                 for ( i = 0; i < size; i++ )\r
 650                         variables[i] = e ->add_variable();\r
 651 \r
 652                 if ( m_dataType == ARRAY ) add_t_links_ARRAY(e,variables,size,alpha_label);\r
 653                 else  if  ( m_dataType == FUNCTION_PIX ) add_t_links_FnPix(e,variables,size,alpha_label);\r
 654                 else  add_t_links_FnCoord(e,variables,size,alpha_label);\r
 655 \r
 656 \r
 657                 if ( m_grid_graph )\r
 658                 {\r
 659                         if ( m_smoothType == ARRAY )\r
 660                         {\r
 661                                 if (m_varying_weights) set_up_expansion_energy_G_ARRAY_VW_pix(size,alpha_label,e,variables,pixels,num);\r
 662                                 else set_up_expansion_energy_G_ARRAY_pix(size,alpha_label,e,variables,pixels,num);\r
 663                         }\r
 664                         else if ( m_smoothType == FUNCTION_PIX ) {\r
 665                 mexErrMsgIdAndTxt("GraphCut:internal", "NOT SUPPORTED YET,exiting!");\r
 666 //                              printf("NOT SUPPORTED YET,exiting!");\r
 667 //                              exit(0);\r
 668                                 //if ( m_smoothType == FUNCTION_PIX ) set_up_expansion_energy_G_FnPix_pix(size,alpha_label,e,variables);\r
 669                         }\r
 670                         else\r
 671                         {\r
 672                 mexErrMsgIdAndTxt("GraphCut:internal", "NOT SUPPORTED YET,exiting!");\r
 673 //                              printf("NOT SUPPORTED YET,exiting!");\r
 674 //                              exit(0);\r
 675                                 //set_up_expansion_energy_G_FnCoord_pix(size,alpha_label,e,variables);\r
 676                         }\r
 677                         \r
 678                 }\r
 679                 else\r
 680                 {\r
 681             mexErrMsgIdAndTxt("GraphCut:internal", "NOT SUPPORTED YET,exiting!");\r
 682 //                      printf("NOT SUPPORTED YET,exiting!");\r
 683 //                      exit(0);\r
 684                         /*if ( m_smoothType == ARRAY ) set_up_expansion_energy_NG_ARRAY_pix(size,alpha_label,e,variables);\r
 685                         else if ( m_smoothType == FUNCTION_PIX ) set_up_expansion_energy_NG_FnPix_pix(size,alpha_label,e,variables);*/\r
 686                 }\r
 687                 \r
 688                 e -> minimize();\r
 689         \r
 690         \r
 691 \r
 692                 for ( i = 0; i < num; i++ )\r
 693                 {\r
 694                         if ( m_labeling[pixels[i]] != alpha_label )\r
 695                         {\r
 696                                 if ( e->get_var(variables[i]) == 0 )\r
 697                                         m_labeling[pixels[i]] = alpha_label;\r
 698                         }\r
 699                         m_lookupPixVar[pixels[i]] = -1;\r
 700                 }\r
 701 \r
 702                 delete [] variables;\r
 703         }\r
 704 \r
 705         delete e;\r
 706 \r
 707         return(compute_energy());\r
 708 }\r
 709 \r
 710 /**************************************************************************************/\r
 711 \r
 712 void GCoptimization::add_t_links_ARRAY(Energy *e,Energy::Var *variables,int size,LabelType alpha_label)\r
 713 {\r
 714         for ( int i = 0; i < size; i++ )\r
 715                 e -> add_term1(variables[i], m_datacost(m_lookupPixVar[i],alpha_label),\r
 716                                              m_datacost(m_lookupPixVar[i],m_labeling[m_lookupPixVar[i]]));\r
 717 \r
 718 }\r
 719 \r
 720 /**************************************************************************************/\r
 721 \r
 722 void GCoptimization::add_t_links_FnPix(Energy *e,Energy::Var *variables,int size,LabelType alpha_label)\r
 723 {\r
 724         for ( int i = 0; i < size; i++ )\r
 725                 e -> add_term1(variables[i], m_dataFnPix(m_lookupPixVar[i],alpha_label),\r
 726                                              m_dataFnPix(m_lookupPixVar[i],m_labeling[m_lookupPixVar[i]]));\r
 727 \r
 728 }\r
 729 /**************************************************************************************/\r
 730 \r
 731 void GCoptimization::add_t_links_FnCoord(Energy *e,Energy::Var *variables,int size,LabelType alpha_label)\r
 732 {\r
 733         int x,y;\r
 734 \r
 735         for ( int i = 0; i < size; i++ )\r
 736         {\r
 737 \r
 738                 y = m_lookupPixVar[i]/m_width;\r
 739                 x = m_lookupPixVar[i] - y*m_width;\r
 740 \r
 741                 e -> add_term1(variables[i], m_dataFnCoord(x,y,alpha_label),\r
 742                                              m_dataFnCoord(x,y,m_labeling[m_lookupPixVar[i]]));\r
 743         }\r
 744 \r
 745 }\r
 746 \r
 747 /**************************************************************************************/\r
 748 \r
 749 void GCoptimization::perform_alpha_expansion(LabelType alpha_label)\r
 750 {\r
 751         PixelType i,size = 0; \r
 752         Energy *e = new Energy(mexErrMsgTxt, truncate_flag);\r
 753         \r
 754 \r
 755         \r
 756         for ( i = 0; i < m_num_pixels; i++ )\r
 757         {\r
 758                 if ( m_labeling[i] != alpha_label )\r
 759                 {\r
 760                         m_lookupPixVar[size] = i;\r
 761                         size++;\r
 762                 }\r
 763         }\r
 764 \r
 765                 \r
 766         if ( size > 0 ) \r
 767         {\r
 768                 Energy::Var *variables = (Energy::Var *) new Energy::Var[size];\r
 769 \r
 770                 for ( i = 0; i < size; i++ )\r
 771                         variables[i] = e ->add_variable();\r
 772 \r
 773                 if ( m_dataType == ARRAY ) add_t_links_ARRAY(e,variables,size,alpha_label);\r
 774                 else  if  ( m_dataType == FUNCTION_PIX ) add_t_links_FnPix(e,variables,size,alpha_label);\r
 775                 else  add_t_links_FnCoord(e,variables,size,alpha_label);\r
 776 \r
 777 \r
 778                 if ( m_grid_graph )\r
 779                 {\r
 780                         if ( m_smoothType == ARRAY )\r
 781                         {\r
 782                                 if (m_varying_weights) set_up_expansion_energy_G_ARRAY_VW(size,alpha_label,e,variables);\r
 783                                 else set_up_expansion_energy_G_ARRAY(size,alpha_label,e,variables);\r
 784                         }\r
 785                         else if ( m_smoothType == FUNCTION_PIX ) set_up_expansion_energy_G_FnPix(size,alpha_label,e,variables);\r
 786                         else  set_up_expansion_energy_G_FnCoord(size,alpha_label,e,variables);\r
 787                         \r
 788                 }\r
 789                 else\r
 790                 {\r
 791                         if ( m_smoothType == ARRAY ) set_up_expansion_energy_NG_ARRAY(size,alpha_label,e,variables);\r
 792                         else if ( m_smoothType == FUNCTION_PIX ) set_up_expansion_energy_NG_FnPix(size,alpha_label,e,variables);\r
 793                 }\r
 794                 \r
 795                 e -> minimize();\r
 796         \r
 797         \r
 798 \r
 799                 for ( i = 0,size = 0; i < m_num_pixels; i++ )\r
 800                 {\r
 801                         if ( m_labeling[i] != alpha_label )\r
 802                         {\r
 803                                 if ( e->get_var(variables[size]) == 0 )\r
 804                                         m_labeling[i] = alpha_label;\r
 805 \r
 806                                 size++;\r
 807                         }\r
 808                 }\r
 809 \r
 810                 delete [] variables;\r
 811         }\r
 812 \r
 813         delete e;\r
 814 }\r
 815 \r
 816 /**********************************************************************************************/\r
 817 /* Performs alpha-expansion for non regular grid graph for case when energy terms are NOT     */\r
 818 /* specified by a function */\r
 819 \r
 820 void GCoptimization::set_up_expansion_energy_NG_ARRAY(int size, LabelType alpha_label,Energy *e,Energy::Var *variables )\r
 821 {\r
 822         EnergyTermType weight;\r
 823         Neighbor *tmp;\r
 824         int i,nPix,pix;;\r
 825 \r
 826 \r
 827 \r
 828         for ( i = size - 1; i >= 0; i-- )\r
 829         {\r
 830                 pix = m_lookupPixVar[i];\r
 831                 m_lookupPixVar[pix] = i;\r
 832 \r
 833                 if ( !m_neighbors[pix].isEmpty() )\r
 834                 {\r
 835                         m_neighbors[pix].setCursorFront();\r
 836                         \r
 837                         while ( m_neighbors[pix].hasNext() )\r
 838                         {\r
 839                                 tmp = (Neighbor *) (m_neighbors[pix].next());\r
 840                                 nPix   = tmp->to_node;\r
 841                                 weight = tmp->weight;\r
 842                                 \r
 843                                 if ( m_labeling[nPix] != alpha_label )\r
 844                                 {\r
 845                                         if ( pix < nPix )\r
 846                                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
 847                                                                   m_smoothcost(alpha_label,alpha_label)*weight,\r
 848                                                                           m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
 849                                                                           m_smoothcost(m_labeling[pix],alpha_label)*weight,\r
 850                                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix])*weight);\r
 851                                 }\r
 852                                 else\r
 853                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,alpha_label)*weight,\r
 854                                                                m_smoothcost(m_labeling[pix],alpha_label)*weight);\r
 855                                 \r
 856                         }\r
 857                 }\r
 858         }\r
 859 \r
 860         \r
 861 }\r
 862 \r
 863 /**********************************************************************************************/\r
 864 /* Performs alpha-expansion for non regular grid graph for case when energy terms are        */\r
 865 /* specified by a function */\r
 866 \r
 867 void GCoptimization::set_up_expansion_energy_NG_FnPix(int size, LabelType alpha_label,Energy *e,Energy::Var *variables )\r
 868 {\r
 869         Neighbor *tmp;\r
 870         int i,nPix,pix;\r
 871         \r
 872 \r
 873 \r
 874         for ( i = size - 1; i >= 0; i-- )\r
 875         {\r
 876                 pix = m_lookupPixVar[i];\r
 877                 m_lookupPixVar[pix] = i;\r
 878 \r
 879 \r
 880                 if ( !m_neighbors[pix].isEmpty() )\r
 881                 {\r
 882                         m_neighbors[pix].setCursorFront();\r
 883                         \r
 884                         while ( m_neighbors[pix].hasNext() )\r
 885                         {\r
 886                                 tmp = (Neighbor *) (m_neighbors[pix].next());\r
 887                                 nPix   = tmp->to_node;\r
 888                                 \r
 889                                 if ( m_labeling[nPix] != alpha_label )\r
 890                                 {\r
 891                                         if ( pix < nPix )\r
 892                                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
 893                                                                   m_smoothFnPix(pix,nPix,alpha_label,alpha_label),\r
 894                                                                           m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
 895                                                                           m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label),\r
 896                                                                           m_smoothFnPix(pix,nPix,m_labeling[pix],m_labeling[nPix]));\r
 897                                 }\r
 898                                 else\r
 899                                         e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
 900                                                                        m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label));\r
 901                                 \r
 902                         }\r
 903                 }\r
 904         }\r
 905 }\r
 906 \r
 907 /**********************************************************************************************/\r
 908 /* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */\r
 909 /* specified by a function */\r
 910 void GCoptimization::set_up_expansion_energy_G_ARRAY_VW(int size, LabelType alpha_label,Energy *e,\r
 911                                                                                                           Energy::Var *variables )\r
 912 {\r
 913         int i,nPix,pix,x,y;\r
 914         EnergyTermType weight;\r
 915 \r
 916 \r
 917         for ( i = size - 1; i >= 0; i-- )\r
 918         {\r
 919                 pix = m_lookupPixVar[i];\r
 920                 y = pix/m_width;\r
 921                 x = pix - y*m_width;\r
 922 \r
 923                 m_lookupPixVar[pix] = i;\r
 924 \r
 925                 if ( x < m_width - 1 )\r
 926                 {\r
 927                         nPix = pix + 1;\r
 928                         weight = m_horizWeights[pix];\r
 929                         if ( m_labeling[nPix] != alpha_label )\r
 930                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
 931                                                   m_smoothcost(alpha_label,alpha_label)*weight,\r
 932                                                           m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
 933                                                           m_smoothcost(m_labeling[pix],alpha_label)*weight,\r
 934                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix])*weight);\r
 935                         else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
 936                                                  m_smoothcost(m_labeling[pix],alpha_label)*weight);\r
 937                 }       \r
 938 \r
 939                 if ( y < m_height - 1 )\r
 940                 {\r
 941                         nPix = pix + m_width;\r
 942                         weight = m_vertWeights[pix];\r
 943                         if ( m_labeling[nPix] != alpha_label )\r
 944                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
 945                                                   m_smoothcost(alpha_label,alpha_label)*weight ,\r
 946                                                           m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
 947                                                           m_smoothcost(m_labeling[pix],alpha_label)*weight ,\r
 948                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix])*weight );\r
 949                         else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
 950                                                  m_smoothcost(m_labeling[pix],alpha_label)*weight);\r
 951                 }       \r
 952                 if ( x > 0 )\r
 953                 {\r
 954                         nPix = pix - 1;\r
 955         \r
 956                         if ( m_labeling[nPix] == alpha_label )\r
 957                            e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*m_horizWeights[nPix],\r
 958                                                  m_smoothcost(m_labeling[pix],alpha_label)*m_horizWeights[nPix]);\r
 959                 }       \r
 960 \r
 961                 if ( y > 0 )\r
 962                 {\r
 963                         nPix = pix - m_width;\r
 964         \r
 965                         if ( m_labeling[nPix] == alpha_label )\r
 966                            e ->add_term1(variables[i],m_smoothcost(alpha_label,alpha_label)*m_vertWeights[nPix],\r
 967                                                  m_smoothcost(m_labeling[pix],alpha_label)*m_vertWeights[nPix]);\r
 968                 }       \r
 969                         \r
 970         }\r
 971         \r
 972 }\r
 973 \r
 974 \r
 975 \r
 976 \r
 977 \r
 978 /**********************************************************************************************/\r
 979 /* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */\r
 980 /* specified by a function */\r
 981 void GCoptimization::set_up_expansion_energy_G_ARRAY(int size, LabelType alpha_label,Energy *e,\r
 982                                                                                                          Energy::Var *variables )\r
 983 {\r
 984         int i,nPix,pix,x,y;\r
 985 \r
 986 \r
 987         for ( i = size - 1; i >= 0; i-- )\r
 988         {\r
 989                 pix = m_lookupPixVar[i];\r
 990                 y = pix/m_width;\r
 991                 x = pix - y*m_width;\r
 992 \r
 993                 m_lookupPixVar[pix] = i;\r
 994 \r
 995                 if ( x < m_width - 1 )\r
 996                 {\r
 997                         nPix = pix + 1;\r
 998 \r
 999                         if ( m_labeling[nPix] != alpha_label )\r
1000                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1001                                                   m_smoothcost(alpha_label,alpha_label),\r
1002                                                           m_smoothcost(alpha_label,m_labeling[nPix]),\r
1003                                                           m_smoothcost(m_labeling[pix],alpha_label),\r
1004                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix]));\r
1005                         else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1006                                                  m_smoothcost(m_labeling[pix],alpha_label));\r
1007                 }       \r
1008 \r
1009                 if ( y < m_height - 1 )\r
1010                 {\r
1011                         nPix = pix + m_width;\r
1012 \r
1013                         if ( m_labeling[nPix] != alpha_label )\r
1014                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1015                                                   m_smoothcost(alpha_label,alpha_label) ,\r
1016                                                           m_smoothcost(alpha_label,m_labeling[nPix]),\r
1017                                                           m_smoothcost(m_labeling[pix],alpha_label) ,\r
1018                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix]) );\r
1019                         else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1020                                                  m_smoothcost(m_labeling[pix],alpha_label));\r
1021                 }       \r
1022                 if ( x > 0 )\r
1023                 {\r
1024                         nPix = pix - 1;\r
1025         \r
1026                         if ( m_labeling[nPix] == alpha_label )\r
1027                            e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1028                                                  m_smoothcost(m_labeling[pix],alpha_label) );\r
1029                 }       \r
1030 \r
1031                 if ( y > 0 )\r
1032                 {\r
1033                         nPix = pix - m_width;\r
1034         \r
1035                         if ( m_labeling[nPix] == alpha_label )\r
1036                            e ->add_term1(variables[i],m_smoothcost(alpha_label,alpha_label),\r
1037                                                  m_smoothcost(m_labeling[pix],alpha_label));\r
1038                 }       \r
1039                         \r
1040         }\r
1041         \r
1042 }\r
1043 \r
1044 /**********************************************************************************************/\r
1045 /* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */\r
1046 /* specified by a function */\r
1047 void GCoptimization::set_up_expansion_energy_G_ARRAY_VW_pix(int size, LabelType alpha_label,Energy *e,\r
1048                                                                                                           Energy::Var *variables, PixelType *pixels, int num )\r
1049 {\r
1050         int i,nPix,pix,x,y;\r
1051         EnergyTermType weight;\r
1052 \r
1053 \r
1054         for ( i = 0; i < num; i++ )\r
1055         {\r
1056                 pix = pixels[i];\r
1057                 if ( m_labeling[pix]!= alpha_label )\r
1058                 {\r
1059                         y = pix/m_width;\r
1060                         x = pix - y*m_width;\r
1061 \r
1062                         if ( x < m_width - 1 )\r
1063                         {\r
1064                                 nPix = pix + 1;\r
1065                                 weight = m_horizWeights[pix];\r
1066                                 if ( m_labeling[nPix] != alpha_label && m_lookupPixVar[nPix] != -1 )\r
1067                                         e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1068                                                   m_smoothcost(alpha_label,alpha_label)*weight,\r
1069                                                           m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
1070                                                           m_smoothcost(m_labeling[pix],alpha_label)*weight,\r
1071                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix])*weight);\r
1072                                 else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
1073                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix])*weight);\r
1074                         }       \r
1075 \r
1076                 if ( y < m_height - 1 )\r
1077                 {\r
1078                         nPix = pix + m_width;\r
1079                         weight = m_vertWeights[pix];\r
1080                         if ( m_labeling[nPix] != alpha_label && m_lookupPixVar[nPix] != -1 )\r
1081                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1082                                                   m_smoothcost(alpha_label,alpha_label)*weight ,\r
1083                                                           m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
1084                                                           m_smoothcost(m_labeling[pix],alpha_label)*weight ,\r
1085                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix])*weight );\r
1086                         else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
1087                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix])*weight);\r
1088                 }       \r
1089                 if ( x > 0 )\r
1090                 {\r
1091                         nPix = pix - 1;\r
1092         \r
1093                         if ( m_labeling[nPix] == alpha_label || m_lookupPixVar[nPix] == -1)\r
1094                            e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*m_horizWeights[nPix],\r
1095                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix])*m_horizWeights[nPix]);\r
1096                 }       \r
1097 \r
1098                 if ( y > 0 )\r
1099                 {\r
1100                         nPix = pix - m_width;\r
1101         \r
1102                         if ( m_labeling[nPix] == alpha_label || m_lookupPixVar[nPix] == -1)\r
1103                            e ->add_term1(variables[i],m_smoothcost(alpha_label,alpha_label)*m_vertWeights[nPix],\r
1104                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix])*m_vertWeights[nPix]);\r
1105                 }       \r
1106                 }       \r
1107         }\r
1108         \r
1109 }\r
1110 \r
1111 /**********************************************************************************************/\r
1112 /* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */\r
1113 /* specified by a function */\r
1114 void GCoptimization::set_up_expansion_energy_G_ARRAY_pix(int size, LabelType alpha_label,Energy *e,\r
1115                                                                                                          Energy::Var *variables, PixelType *pixels,\r
1116                                                                                                          int num)\r
1117 {\r
1118         int i,nPix,pix,x,y;\r
1119 \r
1120 \r
1121         for ( i = 0; i < num; i++ )\r
1122         {\r
1123                 pix = pixels[i];\r
1124                 y = pix/m_width;\r
1125                 x = pix - y*m_width;\r
1126 \r
1127 \r
1128                 if ( m_labeling[pix]!= alpha_label )\r
1129                 {\r
1130                         if ( x < m_width - 1 )\r
1131                         {\r
1132                                 nPix = pix + 1;\r
1133 \r
1134                                 if ( m_labeling[nPix] != alpha_label && m_lookupPixVar[pix] != -1)\r
1135                                         e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1136                                                   m_smoothcost(alpha_label,alpha_label),\r
1137                                                           m_smoothcost(alpha_label,m_labeling[nPix]),\r
1138                                                           m_smoothcost(m_labeling[pix],alpha_label),\r
1139                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix]));\r
1140                                 else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1141                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix]));\r
1142                         }       \r
1143 \r
1144                 if ( y < m_height - 1 )\r
1145                 {\r
1146                         nPix = pix + m_width;\r
1147 \r
1148                         if ( m_labeling[nPix] != alpha_label && m_lookupPixVar[pix] != -1)\r
1149                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1150                                                   m_smoothcost(alpha_label,alpha_label) ,\r
1151                                                           m_smoothcost(alpha_label,m_labeling[nPix]),\r
1152                                                           m_smoothcost(m_labeling[pix],alpha_label) ,\r
1153                                                           m_smoothcost(m_labeling[pix],m_labeling[nPix]) );\r
1154                         else   e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1155                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix]));\r
1156                 }       \r
1157                 if ( x > 0 )\r
1158                 {\r
1159                         nPix = pix - 1;\r
1160         \r
1161                         if ( m_labeling[nPix] == alpha_label || m_lookupPixVar[nPix] == -1)\r
1162                            e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1163                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix]) );\r
1164                 }       \r
1165 \r
1166                 if ( y > 0 )\r
1167                 {\r
1168                         nPix = pix - m_width;\r
1169         \r
1170                         if ( m_labeling[nPix] == alpha_label || m_lookupPixVar[nPix] == -1)\r
1171                            e ->add_term1(variables[i],m_smoothcost(alpha_label,alpha_label),\r
1172                                                  m_smoothcost(m_labeling[pix],m_labeling[nPix]));\r
1173                 }       \r
1174                         \r
1175         }\r
1176         }\r
1177 }\r
1178 \r
1179 \r
1180 /**********************************************************************************************/\r
1181 /* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */\r
1182 /* specified by a function */\r
1183 void GCoptimization::set_up_expansion_energy_G_FnPix(int size, LabelType alpha_label,Energy *e,\r
1184                                                                                                          Energy::Var *variables )\r
1185 {\r
1186         int i,nPix,pix,x,y;\r
1187 \r
1188 \r
1189         for ( i = size - 1; i >= 0; i-- )\r
1190         {\r
1191                 pix = m_lookupPixVar[i];\r
1192                 y = pix/m_width;\r
1193                 x = pix - y*m_width;\r
1194 \r
1195                 m_lookupPixVar[pix] = i;\r
1196 \r
1197                 if ( x < m_width - 1 )\r
1198                 {\r
1199                         nPix = pix + 1;\r
1200 \r
1201                         if ( m_labeling[nPix] != alpha_label )\r
1202                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1203                                                   m_smoothFnPix(pix,nPix,alpha_label,alpha_label),\r
1204                                                           m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1205                                                           m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label),\r
1206                                                           m_smoothFnPix(pix,nPix,m_labeling[pix],m_labeling[nPix]));\r
1207                         else   e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1208                                                  m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label));\r
1209                 }       \r
1210 \r
1211                 if ( y < m_height - 1 )\r
1212                 {\r
1213                         nPix = pix + m_width;\r
1214 \r
1215                         if ( m_labeling[nPix] != alpha_label )\r
1216                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1217                                                   m_smoothFnPix(pix,nPix,alpha_label,alpha_label) ,\r
1218                                                           m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1219                                                           m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label) ,\r
1220                                                           m_smoothFnPix(pix,nPix,m_labeling[pix],m_labeling[nPix]) );\r
1221                         else   e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1222                                                  m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label));\r
1223                 }       \r
1224                 if ( x > 0 )\r
1225                 {\r
1226                         nPix = pix - 1;\r
1227         \r
1228                         if ( m_labeling[nPix] == alpha_label )\r
1229                            e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1230                                                  m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label) );\r
1231                 }       \r
1232 \r
1233                 if ( y > 0 )\r
1234                 {\r
1235                         nPix = pix - m_width;\r
1236         \r
1237                         if ( m_labeling[nPix] == alpha_label )\r
1238                            e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,alpha_label),\r
1239                                                  m_smoothFnPix(pix,nPix,m_labeling[pix],alpha_label));\r
1240                 }       \r
1241                         \r
1242         }\r
1243         \r
1244 }\r
1245 \r
1246 /**********************************************************************************************/\r
1247 /* Performs alpha-expansion for  regular grid graph for case when energy terms are NOT        */\r
1248 /* specified by a function */\r
1249 void GCoptimization::set_up_expansion_energy_G_FnCoord(int size, LabelType alpha_label,Energy *e,\r
1250                                                                                                          Energy::Var *variables )\r
1251 {\r
1252         int i,nPix,pix,x,y;\r
1253 \r
1254 \r
1255         for ( i = size - 1; i >= 0; i-- )\r
1256         {\r
1257                 pix = m_lookupPixVar[i];\r
1258                 y = pix/m_width;\r
1259                 x = pix - y*m_width;\r
1260 \r
1261                 m_lookupPixVar[pix] = i;\r
1262 \r
1263                 if ( x < m_width - 1 )\r
1264                 {\r
1265                         nPix = pix + 1;\r
1266 \r
1267                         if ( m_labeling[nPix] != alpha_label )\r
1268                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1269                                                   m_horz_cost(x,y,alpha_label,alpha_label),\r
1270                                                           m_horz_cost(x,y,alpha_label,m_labeling[nPix]),\r
1271                                                           m_horz_cost(x,y,m_labeling[pix],alpha_label),\r
1272                                                           m_horz_cost(x,y,m_labeling[pix],m_labeling[nPix]));\r
1273                         else   e ->add_term1(variables[i],m_horz_cost(x,y,alpha_label,alpha_label),\r
1274                                                  m_horz_cost(x,y,m_labeling[pix],alpha_label));\r
1275                 }       \r
1276 \r
1277                 if ( y < m_height - 1 )\r
1278                 {\r
1279                         nPix = pix + m_width;\r
1280 \r
1281                         if ( m_labeling[nPix] != alpha_label )\r
1282                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1283                                                   m_vert_cost(x,y,alpha_label,alpha_label) ,\r
1284                                                           m_vert_cost(x,y,alpha_label,m_labeling[nPix]),\r
1285                                                           m_vert_cost(x,y,m_labeling[pix],alpha_label) ,\r
1286                                                           m_vert_cost(x,y,m_labeling[pix],m_labeling[nPix]) );\r
1287                         else   e ->add_term1(variables[i],m_vert_cost(x,y,alpha_label,alpha_label),\r
1288                                                  m_vert_cost(x,y,m_labeling[pix],alpha_label));\r
1289                 }       \r
1290                 if ( x > 0 )\r
1291                 {\r
1292                         nPix = pix - 1;\r
1293         \r
1294                         if ( m_labeling[nPix] == alpha_label )\r
1295                            e ->add_term1(variables[i],m_horz_cost(x-1,y,alpha_label,m_labeling[nPix]),\r
1296                                                  m_horz_cost(x-1,y,m_labeling[pix],alpha_label) );\r
1297                 }       \r
1298 \r
1299                 if ( y > 0 )\r
1300                 {\r
1301                         nPix = pix - m_width;\r
1302         \r
1303                         if ( m_labeling[nPix] == alpha_label )\r
1304                            e ->add_term1(variables[i],m_vert_cost(x,y-1,alpha_label,alpha_label),\r
1305                                                  m_vert_cost(x,y-1,m_labeling[pix],alpha_label));\r
1306                 }       \r
1307         }\r
1308 }\r
1309 \r
1310 \r
1311 \r
1312 /**************************************************************************************/\r
1313 \r
1314 GCoptimization::EnergyType GCoptimization::oneExpansionIteration()\r
1315 {\r
1316         int next;\r
1317    \r
1318         terminateOnError( m_dataType == NONE,"You have to set up the data cost before running optimization");\r
1319         terminateOnError( m_smoothType == NONE,"You have to set up the smoothness cost before running optimization");\r
1320 \r
1321         if (m_random_label_order) scramble_label_table();\r
1322         \r
1323 \r
1324         for (next = 0;  next < m_num_labels;  next++ )\r
1325                 perform_alpha_expansion(m_labelTable[next]);\r
1326         \r
1327         \r
1328         return(compute_energy());\r
1329 }\r
1330 \r
1331 /**************************************************************************************/\r
1332 \r
1333 void GCoptimization::setNeighbors(PixelType pixel1, int pixel2, EnergyTermType weight)\r
1334 {\r
1335 \r
1336         assert(pixel1 < m_num_pixels && pixel1 >= 0 && pixel2 < m_num_pixels && pixel2 >= 0);\r
1337         assert(m_grid_graph == 0);\r
1338 \r
1339         Neighbor *temp1 = (Neighbor *) new Neighbor;\r
1340         Neighbor *temp2 = (Neighbor *) new Neighbor;\r
1341 \r
1342         temp1->weight  = weight;\r
1343         temp1->to_node = pixel2;\r
1344 \r
1345         temp2->weight  = weight;\r
1346         temp2->to_node = pixel1;\r
1347 \r
1348         m_neighbors[pixel1].addFront(temp1);\r
1349         m_neighbors[pixel2].addFront(temp2);\r
1350         \r
1351 }\r
1352 \r
1353 /**************************************************************************************/\r
1354 \r
1355 void GCoptimization::setNeighbors(PixelType pixel1, int pixel2)\r
1356 {\r
1357 \r
1358         assert(pixel1 < m_num_pixels && pixel1 >= 0 && pixel2 < m_num_pixels && pixel2 >= 0);\r
1359         assert(m_grid_graph == 0);\r
1360         \r
1361 \r
1362         Neighbor *temp1 = (Neighbor *) new Neighbor;\r
1363         Neighbor *temp2 = (Neighbor *) new Neighbor;\r
1364 \r
1365         temp1->weight  = (EnergyTermType) 1;\r
1366         temp1->to_node = pixel2;\r
1367 \r
1368         temp2->weight  = (EnergyTermType) 1;\r
1369         temp2->to_node = pixel1;\r
1370 \r
1371         m_neighbors[pixel1].addFront(temp1);\r
1372         m_neighbors[pixel2].addFront(temp2);\r
1373         \r
1374 }\r
1375 \r
1376 /**************************************************************************************/\r
1377 \r
1378 GCoptimization::~GCoptimization()\r
1379 {\r
1380 \r
1381     // <!-- bagon\r
1382     // mexWarnMsgTxt("Calling destructor"); /* bagon added */\r
1383     if ( m_dataType == ARRAY )\r
1384         delete [] m_datacost;\r
1385     if ( m_smoothType == ARRAY )\r
1386         delete [] m_smoothcost;\r
1387     if ( m_varying_weights == 1 ) {\r
1388         delete [] m_vertWeights;\r
1389         delete [] m_horizWeights;\r
1390     }\r
1391     // bagon -->\r
1392     \r
1393         if ( m_deleteLabeling ) \r
1394                 delete [] m_labeling;\r
1395 \r
1396 \r
1397         if ( m_dataInput == SET_INDIVIDUALLY )\r
1398                 delete [] m_datacost;\r
1399 \r
1400         if ( m_smoothInput == SET_INDIVIDUALLY )\r
1401                         delete [] m_smoothcost;\r
1402                 \r
1403         if ( ! m_grid_graph )\r
1404                 delete [] m_neighbors;                  \r
1405 \r
1406 \r
1407         delete [] m_labelTable;\r
1408         delete [] m_lookupPixVar;\r
1409                         \r
1410 }\r
1411 \r
1412 \r
1413 /**************************************************************************************/\r
1414 \r
1415 GCoptimization::EnergyType GCoptimization::swap(int max_num_iterations)\r
1416 {\r
1417         return(start_swap(max_num_iterations)); \r
1418 }\r
1419 \r
1420 /**************************************************************************************/\r
1421 \r
1422 \r
1423 GCoptimization::EnergyType GCoptimization::swap()\r
1424 {\r
1425         return(start_swap(MAX_INTT));\r
1426 }\r
1427 \r
1428 /**************************************************************************************/\r
1429 \r
1430 \r
1431 GCoptimization::EnergyType GCoptimization::start_swap(int max_num_iterations )\r
1432 {\r
1433         \r
1434         int curr_cycle = 1;\r
1435         EnergyType new_energy,old_energy;\r
1436         \r
1437 \r
1438         new_energy = compute_energy();\r
1439 \r
1440         old_energy = (new_energy+1)*10; // BAGON changed init value to exceed current energy by factor of 10 (thanks to A. Khan)    \r
1441 \r
1442         while ( old_energy > new_energy  && curr_cycle <= max_num_iterations)\r
1443         {\r
1444                 old_energy = new_energy;\r
1445                 new_energy = oneSwapIteration();\r
1446                 \r
1447                 curr_cycle++;   \r
1448         }\r
1449 \r
1450         return(new_energy);\r
1451 }\r
1452 \r
1453 /**************************************************************************************/\r
1454 \r
1455 GCoptimization::EnergyType GCoptimization::oneSwapIteration()\r
1456 {\r
1457         int next,next1;\r
1458    \r
1459         if (m_random_label_order) scramble_label_table();\r
1460                 \r
1461 \r
1462         for (next = 0;  next < m_num_labels;  next++ )\r
1463                 for (next1 = m_num_labels - 1;  next1 >= 0;  next1-- )\r
1464                         if ( m_labelTable[next] < m_labelTable[next1] )\r
1465                                 perform_alpha_beta_swap(m_labelTable[next],m_labelTable[next1]);\r
1466 \r
1467         return(compute_energy());\r
1468 }\r
1469 \r
1470 /**************************************************************************************/\r
1471 \r
1472 GCoptimization::EnergyType GCoptimization::alpha_beta_swap(LabelType alpha_label, LabelType beta_label)\r
1473 {\r
1474         terminateOnError( alpha_label < 0 || alpha_label >= m_num_labels || beta_label < 0 || beta_label >= m_num_labels,\r
1475                 "Illegal Label to Expand On");\r
1476         perform_alpha_beta_swap(alpha_label,beta_label);\r
1477         return(compute_energy());\r
1478 }\r
1479 /**************************************************************************************/\r
1480 \r
1481 void GCoptimization::add_t_links_ARRAY_swap(Energy *e,Energy::Var *variables,int size,\r
1482                                                                                         LabelType alpha_label, LabelType beta_label,\r
1483                                                                                         PixelType *pixels)\r
1484 {\r
1485         for ( int i = 0; i < size; i++ )\r
1486                 e -> add_term1(variables[i], m_datacost(pixels[i],alpha_label),\r
1487                                                                          m_datacost(pixels[i],beta_label));\r
1488 \r
1489 }\r
1490         \r
1491 /**************************************************************************************/\r
1492 \r
1493 void GCoptimization::add_t_links_FnPix_swap(Energy *e,Energy::Var *variables,int size,\r
1494                                                                                         LabelType alpha_label, LabelType beta_label,\r
1495                                                                                         PixelType *pixels)\r
1496 {\r
1497         for ( int i = 0; i < size; i++ )\r
1498                 e -> add_term1(variables[i], m_dataFnPix(pixels[i],alpha_label),\r
1499                                                                          m_dataFnPix(pixels[i],beta_label));\r
1500 \r
1501 }\r
1502 /**************************************************************************************/\r
1503 \r
1504 void GCoptimization::add_t_links_FnCoord_swap(Energy *e,Energy::Var *variables,int size,\r
1505                                                                                           LabelType alpha_label, LabelType beta_label,\r
1506                                                                                           PixelType *pixels)\r
1507 {\r
1508         int x,y;\r
1509 \r
1510         for ( int i = 0; i < size; i++ )\r
1511         {\r
1512 \r
1513                 y = pixels[i]/m_width;\r
1514                 x = pixels[i] - y*m_width;\r
1515 \r
1516                 e -> add_term1(variables[i], m_dataFnCoord(x,y,alpha_label),m_dataFnCoord(x,y,beta_label));\r
1517         }\r
1518 }\r
1519 \r
1520 \r
1521 /**************************************************************************************/\r
1522 \r
1523 void  GCoptimization::perform_alpha_beta_swap(LabelType alpha_label, LabelType beta_label)\r
1524 {\r
1525         PixelType i,size = 0;\r
1526         Energy *e = new Energy(mexErrMsgTxt, truncate_flag);\r
1527         PixelType *pixels = new PixelType[m_num_pixels];\r
1528         \r
1529 \r
1530         for ( i = 0; i < m_num_pixels; i++ )\r
1531         {\r
1532                 if ( m_labeling[i] == alpha_label || m_labeling[i] == beta_label)\r
1533                 {\r
1534                         pixels[size]    = i;\r
1535                         m_lookupPixVar[i] = size;\r
1536                         size++;\r
1537                 }\r
1538         }\r
1539 \r
1540         if ( size == 0 )\r
1541         {\r
1542                 delete e;\r
1543                 delete [] pixels;\r
1544                 return;\r
1545         }\r
1546 \r
1547 \r
1548         Energy::Var *variables = (Energy::Var *) new Energy::Var[size];\r
1549 \r
1550 \r
1551         for ( i = 0; i < size; i++ )\r
1552                 variables[i] = e ->add_variable();\r
1553 \r
1554         if ( m_dataType == ARRAY ) add_t_links_ARRAY_swap(e,variables,size,alpha_label,beta_label,pixels);\r
1555         else  if  ( m_dataType == FUNCTION_PIX ) add_t_links_FnPix_swap(e,variables,size,alpha_label,beta_label,pixels);\r
1556         else  add_t_links_FnCoord_swap(e,variables,size,alpha_label,beta_label,pixels);\r
1557 \r
1558 \r
1559 \r
1560         if ( m_grid_graph )\r
1561         {\r
1562                 if ( m_smoothType == ARRAY )\r
1563                 {\r
1564                         if (m_varying_weights) set_up_swap_energy_G_ARRAY_VW(size,alpha_label,beta_label,pixels,e,variables);\r
1565                         else set_up_swap_energy_G_ARRAY(size,alpha_label,beta_label,pixels,e,variables);\r
1566                 }\r
1567                 else if ( m_smoothType == FUNCTION_PIX ) set_up_swap_energy_G_FnPix(size,alpha_label,beta_label,pixels,e,variables);\r
1568                 else  set_up_swap_energy_G_FnCoord(size,alpha_label,beta_label,pixels,e,variables);\r
1569                         \r
1570         }\r
1571         else\r
1572         {\r
1573                 if ( m_smoothType == ARRAY ) set_up_swap_energy_NG_ARRAY(size,alpha_label,beta_label,pixels,e,variables);\r
1574                 else if ( m_smoothType == FUNCTION_PIX ) set_up_swap_energy_NG_FnPix(size,alpha_label,beta_label,pixels,e,variables);\r
1575         }\r
1576                 \r
1577 \r
1578         e -> minimize();\r
1579 \r
1580         for ( i = 0; i < size; i++ )\r
1581                 if ( e->get_var(variables[i]) == 0 )\r
1582                         m_labeling[pixels[i]] = alpha_label;\r
1583                 else m_labeling[pixels[i]] = beta_label;\r
1584 \r
1585 \r
1586         delete [] variables;\r
1587         delete [] pixels;\r
1588         delete e;\r
1589 \r
1590 }\r
1591 \r
1592 /**************************************************************************************/\r
1593 \r
1594 void GCoptimization::set_up_swap_energy_NG_ARRAY(int size,LabelType alpha_label,LabelType beta_label,\r
1595                                                                                                  PixelType *pixels,Energy* e, Energy::Var *variables)\r
1596 {\r
1597         PixelType nPix,pix,i;\r
1598         EnergyTermType weight;\r
1599         Neighbor *tmp;\r
1600         \r
1601 \r
1602 \r
1603         for ( i = 0; i < size; i++ )\r
1604         {\r
1605                 pix = pixels[i];\r
1606                 if ( !m_neighbors[pix].isEmpty() )\r
1607                 {\r
1608                         m_neighbors[pix].setCursorFront();\r
1609                         \r
1610                         while ( m_neighbors[pix].hasNext() )\r
1611                         {\r
1612                                 tmp = (Neighbor *) (m_neighbors[pix].next());\r
1613                                 nPix   = tmp->to_node;\r
1614                                 weight = tmp->weight;\r
1615                                 \r
1616                                 if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1617                                 {\r
1618                                         if ( pix < nPix )\r
1619                                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1620                                                                   m_smoothcost(alpha_label,alpha_label)*weight,\r
1621                                                                           m_smoothcost(alpha_label,beta_label)*weight,\r
1622                                                                           m_smoothcost(beta_label,alpha_label)*weight,\r
1623                                                                           m_smoothcost(beta_label,beta_label)*weight);\r
1624                                 }\r
1625                                 else\r
1626                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
1627                                                                        m_smoothcost(beta_label,m_labeling[nPix])*weight);\r
1628                         }\r
1629                 }\r
1630         }\r
1631 }\r
1632 \r
1633 /**************************************************************************************/\r
1634 \r
1635 void GCoptimization::set_up_swap_energy_NG_FnPix(int size,LabelType alpha_label,LabelType beta_label,\r
1636                                                                                                  PixelType *pixels,Energy* e, Energy::Var *variables)\r
1637 {\r
1638         PixelType nPix,pix,i;\r
1639         Neighbor *tmp;\r
1640         \r
1641 \r
1642         for ( i = 0; i < size; i++ )\r
1643         {\r
1644                 pix = pixels[i];\r
1645                 if ( !m_neighbors[pix].isEmpty() )\r
1646                 {\r
1647                         m_neighbors[pix].setCursorFront();\r
1648                         \r
1649                         while ( m_neighbors[pix].hasNext() )\r
1650                         {\r
1651                                 tmp = (Neighbor *) (m_neighbors[pix].next());\r
1652                                 nPix   = tmp->to_node;\r
1653                                 \r
1654                                 if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1655                                 {\r
1656                                         if ( pix < nPix )\r
1657                                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1658                                                                   m_smoothFnPix(pix,nPix,alpha_label,alpha_label),\r
1659                                                                           m_smoothFnPix(pix,nPix,alpha_label,beta_label),\r
1660                                                                           m_smoothFnPix(pix,nPix,beta_label,alpha_label),\r
1661                                                                           m_smoothFnPix(pix,nPix,beta_label,beta_label) );\r
1662                                 }\r
1663                                 else\r
1664                                         e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1665                                                                        m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));\r
1666                         }\r
1667                 }\r
1668         }\r
1669 }\r
1670 \r
1671 /**************************************************************************************/\r
1672 \r
1673 void GCoptimization::set_up_swap_energy_G_FnPix(int size,LabelType alpha_label,LabelType beta_label,\r
1674                                                                                                 PixelType *pixels,Energy* e, Energy::Var *variables)\r
1675 {\r
1676         PixelType nPix,pix,i,x,y;\r
1677 \r
1678         \r
1679         for ( i = 0; i < size; i++ )\r
1680         {\r
1681                 pix = pixels[i];\r
1682                 y = pix/m_width;\r
1683                 x = pix - y*m_width;\r
1684 \r
1685                 if ( x > 0 )\r
1686                 {\r
1687                         nPix = pix - 1;\r
1688         \r
1689                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1690                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1691                                               m_smoothFnPix(pix,nPix,alpha_label,alpha_label),\r
1692                                                           m_smoothFnPix(pix,nPix,alpha_label,beta_label),\r
1693                                                           m_smoothFnPix(pix,nPix,beta_label,alpha_label),\r
1694                                                           m_smoothFnPix(pix,nPix,beta_label,beta_label) );\r
1695         \r
1696                                 else\r
1697                                         e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1698                                                                m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));\r
1699         \r
1700                 }       \r
1701                 if ( y > 0 )\r
1702                 {\r
1703                         nPix = pix - m_width;\r
1704                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1705                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1706                                               m_smoothFnPix(pix,nPix,alpha_label,alpha_label),\r
1707                                                           m_smoothFnPix(pix,nPix,alpha_label,beta_label),\r
1708                                                           m_smoothFnPix(pix,nPix,beta_label,alpha_label),\r
1709                                                           m_smoothFnPix(pix,nPix,beta_label,beta_label) );\r
1710         \r
1711                                 else\r
1712                                         e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1713                                                                m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));\r
1714                 }       \r
1715 \r
1716                 if ( x < m_width - 1 )\r
1717                 {\r
1718                         nPix = pix + 1;\r
1719         \r
1720                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )\r
1721                                         e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1722                                                                    m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));\r
1723                 }       \r
1724 \r
1725                 if ( y < m_height - 1 )\r
1726                 {\r
1727                         nPix = pix + m_width;\r
1728         \r
1729                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )\r
1730                                 e ->add_term1(variables[i],m_smoothFnPix(pix,nPix,alpha_label,m_labeling[nPix]),\r
1731                                                            m_smoothFnPix(pix,nPix,beta_label,m_labeling[nPix]));\r
1732 \r
1733                 }\r
1734         }\r
1735 }\r
1736 /**************************************************************************************/\r
1737 \r
1738 void GCoptimization::set_up_swap_energy_G_FnCoord(int size,LabelType alpha_label,LabelType beta_label,PixelType *pixels,\r
1739                                                                                          Energy* e, Energy::Var *variables)\r
1740 {  \r
1741         PixelType nPix,pix,i,x,y;\r
1742 \r
1743         \r
1744         for ( i = 0; i < size; i++ )\r
1745         {\r
1746                 pix = pixels[i];\r
1747                 y = pix/m_width;\r
1748                 x = pix - y*m_width;\r
1749 \r
1750                 if ( x > 0 )\r
1751                 {\r
1752                         nPix = pix - 1;\r
1753         \r
1754                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1755                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1756                                               m_horz_cost(x-1,y,alpha_label,alpha_label),\r
1757                                                           m_horz_cost(x-1,y,alpha_label,beta_label),\r
1758                                                           m_horz_cost(x-1,y,beta_label,alpha_label),\r
1759                                                           m_horz_cost(x-1,y,beta_label,beta_label) );\r
1760         \r
1761                                 else\r
1762                                         e ->add_term1(variables[i],m_horz_cost(x-1,y,alpha_label,m_labeling[nPix]),\r
1763                                                                    m_horz_cost(x-1,y,beta_label,m_labeling[nPix]));\r
1764         \r
1765                 }       \r
1766                 if ( y > 0 )\r
1767                 {\r
1768                         nPix = pix - m_width;\r
1769                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1770                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1771                                               m_vert_cost(x,y-1,alpha_label,alpha_label),\r
1772                                                           m_vert_cost(x,y-1,alpha_label,beta_label),\r
1773                                                           m_vert_cost(x,y-1,beta_label,alpha_label),\r
1774                                                           m_vert_cost(x,y-1,beta_label,beta_label) );\r
1775         \r
1776                                 else\r
1777                                         e ->add_term1(variables[i],m_vert_cost(x,y-1,alpha_label,m_labeling[nPix]),\r
1778                                                                m_vert_cost(x,y-1,beta_label,m_labeling[nPix]));\r
1779                 }       \r
1780 \r
1781                 if ( x < m_width - 1 )\r
1782                 {\r
1783                         nPix = pix + 1;\r
1784         \r
1785                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )\r
1786                                         e ->add_term1(variables[i],m_horz_cost(x,y,alpha_label,m_labeling[nPix]),\r
1787                                                                    m_horz_cost(x,y,beta_label,m_labeling[nPix]));\r
1788                 }       \r
1789 \r
1790                 if ( y < m_height - 1 )\r
1791                 {\r
1792                         nPix = pix + m_width;\r
1793         \r
1794                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )\r
1795                                 e ->add_term1(variables[i],m_vert_cost(x,y,alpha_label,m_labeling[nPix]),\r
1796                                                            m_vert_cost(x,y,beta_label,m_labeling[nPix]));\r
1797 \r
1798                 }\r
1799         }\r
1800 }\r
1801 \r
1802 /**************************************************************************************/\r
1803 \r
1804 void GCoptimization::set_up_swap_energy_G_ARRAY_VW(int size,LabelType alpha_label,LabelType beta_label,\r
1805                                                                                                    PixelType *pixels,Energy* e, Energy::Var *variables)\r
1806 {\r
1807         PixelType nPix,pix,i,x,y;\r
1808         EnergyTermType weight;  \r
1809 \r
1810 \r
1811 \r
1812         for ( i = 0; i < size; i++ )\r
1813         {\r
1814                 pix = pixels[i];\r
1815                 y = pix/m_width;\r
1816                 x = pix - y*m_width;\r
1817 \r
1818                 if ( x > 0 )\r
1819                 {\r
1820                         nPix = pix - 1;\r
1821                         weight = m_horizWeights[nPix];\r
1822         \r
1823                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1824                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1825                                               m_smoothcost(alpha_label,alpha_label)*weight,\r
1826                                                           m_smoothcost(alpha_label,beta_label)*weight,\r
1827                                                           m_smoothcost(beta_label,alpha_label)*weight,\r
1828                                                           m_smoothcost(beta_label,beta_label)*weight );\r
1829         \r
1830                                 else\r
1831                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
1832                                                                m_smoothcost(beta_label,m_labeling[nPix])*weight);\r
1833         \r
1834                 }       \r
1835                 if ( y > 0 )\r
1836                 {\r
1837                         nPix = pix - m_width;\r
1838                         weight = m_vertWeights[nPix];\r
1839 \r
1840                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1841                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1842                                               m_smoothcost(alpha_label,alpha_label)*weight,\r
1843                                                           m_smoothcost(alpha_label,beta_label)*weight,\r
1844                                                           m_smoothcost(beta_label,alpha_label)*weight,\r
1845                                                           m_smoothcost(beta_label,beta_label)*weight );\r
1846         \r
1847                                 else\r
1848                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*weight,\r
1849                                                                m_smoothcost(beta_label,m_labeling[nPix])*weight);\r
1850                 }       \r
1851 \r
1852                 if ( x < m_width - 1 )\r
1853                 {\r
1854                         nPix = pix + 1;\r
1855         \r
1856                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )\r
1857                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*m_horizWeights[pix],\r
1858                                                                    m_smoothcost(beta_label,m_labeling[nPix])*m_horizWeights[pix]);\r
1859                 }       \r
1860 \r
1861                 if ( y < m_height - 1 )\r
1862                 {\r
1863                         nPix = pix + m_width;\r
1864         \r
1865                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )\r
1866                                 e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix])*m_vertWeights[pix],\r
1867                                                            m_smoothcost(beta_label,m_labeling[nPix])*m_vertWeights[pix]);\r
1868 \r
1869                 }\r
1870         }\r
1871 }\r
1872 \r
1873 /**************************************************************************************/\r
1874 \r
1875 void GCoptimization::set_up_swap_energy_G_ARRAY(int size,LabelType alpha_label,LabelType beta_label,\r
1876                                                                                            PixelType *pixels,Energy* e, Energy::Var *variables)\r
1877 \r
1878 {\r
1879         PixelType nPix,pix,i,x,y;\r
1880         \r
1881 \r
1882 \r
1883         for ( i = 0; i < size; i++ )\r
1884         {\r
1885                 pix = pixels[i];\r
1886                 y = pix/m_width;\r
1887                 x = pix - y*m_width;\r
1888 \r
1889                 if ( x > 0 )\r
1890                 {\r
1891                         nPix = pix - 1;\r
1892         \r
1893         \r
1894                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1895                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1896                                               m_smoothcost(alpha_label,alpha_label),\r
1897                                                           m_smoothcost(alpha_label,beta_label),\r
1898                                                           m_smoothcost(beta_label,alpha_label),\r
1899                                                           m_smoothcost(beta_label,beta_label) );\r
1900         \r
1901                                 else\r
1902                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1903                                                                m_smoothcost(beta_label,m_labeling[nPix]));\r
1904         \r
1905                 }       \r
1906                 if ( y > 0 )\r
1907                 {\r
1908                         nPix = pix - m_width;\r
1909 \r
1910                         if ( m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label)\r
1911                                 e ->add_term2(variables[i],variables[m_lookupPixVar[nPix]],\r
1912                                               m_smoothcost(alpha_label,alpha_label),\r
1913                                                           m_smoothcost(alpha_label,beta_label),\r
1914                                                           m_smoothcost(beta_label,alpha_label),\r
1915                                                           m_smoothcost(beta_label,beta_label) );\r
1916         \r
1917                                 else\r
1918                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1919                                                                m_smoothcost(beta_label,m_labeling[nPix]));\r
1920                 }       \r
1921 \r
1922                 if ( x < m_width - 1 )\r
1923                 {\r
1924                         nPix = pix + 1;\r
1925         \r
1926                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label) )\r
1927                                         e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1928                                                                    m_smoothcost(beta_label,m_labeling[nPix]));\r
1929                 }       \r
1930 \r
1931                 if ( y < m_height - 1 )\r
1932                 {\r
1933                         nPix = pix + m_width;\r
1934         \r
1935                         if ( !(m_labeling[nPix] == alpha_label || m_labeling[nPix] == beta_label))\r
1936                                 e ->add_term1(variables[i],m_smoothcost(alpha_label,m_labeling[nPix]),\r
1937                                                            m_smoothcost(beta_label,m_labeling[nPix]));\r
1938 \r
1939                 }\r
1940         }\r
1941 }\r
1942 \r
1943 \r
1944 /**************************************************************************************/\r
1945 \r
1946 void GCoptimization::setLabelOrder(bool RANDOM_LABEL_ORDER)\r
1947 {\r
1948         m_random_label_order = RANDOM_LABEL_ORDER;\r
1949 }\r
1950 \r
1951 /****************************************************************************/\r
1952 /* This procedure checks if an error has occured, terminates program if yes */\r
1953 \r
1954 void GCoptimization::terminateOnError(bool error_condition,const char *message)\r
1955 \r
1956 { \r
1957    if  (error_condition) \r
1958    {\r
1959        mexErrMsgIdAndTxt("GraphCut:internal_error","\nGCoptimization error: %s\n", message);\r
1960     }\r
1961 }\r