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

   1 \r
   2 #include "mex.h"\r
   3 #include "GCoptimization.h"\r
   4 #include "GraphCut.h"\r
   5 #include <stdlib.h>\r
   6 \r
   7 /* Declarations */\r
   8 void SetLabels(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);\r
   9 void GetLabels(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);\r
  10 void ABswaps(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);\r
  11 void Expand(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);\r
  12 void Energy(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);\r
  13 void Truncate(GCoptimization *MyGraph, int nout, mxArray *pout[], int nin, const mxArray *pin[]);\r
  14 GCoptimization* GetGCHandle(const mxArray *x);    /* extract ahndle from mxArry */\r
  15 \r
  16 \r
  17 /*\r
  18  * Matlab wrapper for Weksler graph cut implementation\r
  19  * \r
  20  * GCoptimization cde by Olga Veksler.\r
  21  * Wrapper code by Shai Bagon.\r
  22  *\r
  23  *\r
  24  *   Performing Graph Cut operations on a 2D grid.\r
  25  *   \r
  26  *   Usage:\r
  27  *       [gch ...] = GraphCutMex(gch, mode ...);\r
  28  *   \r
  29  *   Notes:\r
  30  *   1. Data types are crucial!\r
  31  *   2. The embedded implementation treat matrices as a row stack, while\r
  32  *       matlab treats them as column stack; thus there is a need to\r
  33  *       transpose the label matrices and the indices passed to expand\r
  34  *       algorithm.\r
  35  *\r
  36  *   Inputs:\r
  37  *   - gch: a valid Graph Cut handle (use GraphCutConstr to create a handle).\r
  38  *   - mode: a char specifying mode of operation. See details below.\r
  39  *\r
  40  *   Output:\r
  41  *   - gch: A handle to the constructed graph. Handle this handle with care\r
  42  *              and don't forget to close it in the end!\r
  43  *\r
  44  *   Possible modes:\r
  45  *\r
  46  *   - 's': Set labels\r
  47  *           [gch] = GraphCutMex(gch, 's', labels)\r
  48  *\r
  49  *       Inputs:\r
  50  *           - labels: a width*height array of type int32, containing a\r
  51  *              label per pixel. note that labels values must be is range\r
  52  *              [0..num_labels-1]. \r
  53  *\r
  54  *   - 'g': Get current labeling\r
  55  *           [gch labels] = GraphCutMex(gch, 'g')\r
  56  *\r
  57  *       Outputs:\r
  58  *           - labels: a width*height array of type int32, containing a\r
  59  *              label per pixel. note that labels values must be is range\r
  60  *              [0..num_labels-1].\r
  61  *\r
  62  *   - 'n': Get current values of energy terms\r
  63  *           [gch se de] = GraphCutMex(gch, 'n')\r
  64  *\r
  65  *       Outputs:\r
  66  *           - se: Smoothness energy term.\r
  67  *           - de: Data energy term.\r
  68  *\r
  69  *   - 'e': Perform labels expansion\r
  70  *           [gch labels] = GraphCutMex(gch, 'e')\r
  71  *           [gch labels] = GraphCutMex(gch, 'e', iter)\r
  72  *           [gch labels] = GraphCutMex(gch, 'e', [], label)\r
  73  *           [gch labels] = GraphCutMex(gch, 'e', [], label, indices)\r
  74  *\r
  75  *       When no inputs are provided, GraphCut performs expansion steps\r
  76  *       until it converges.\r
  77  *\r
  78  *       Inputs:\r
  79  *           - iter: a double scalar, the maximum number of expand\r
  80  *                      iterations to perform.\r
  81  *           - label: int32 scalar denoting the label for which to perfom\r
  82  *                        expand step.\r
  83  *           - indices: int32 array of linear indices of pixels for which\r
  84  *                            expand step is computed. indices _MUST_ be zero offset and not\r
  85  *                            one offset like matlab usuall indexing system, i.e., to include\r
  86  *                            the first pixel in the expand indices array must contain 0 and\r
  87  *                            NOT 1!\r
  88  *\r
  89  *       Outputs:\r
  90  *           - labels: a width*height array of type int32, containing a\r
  91  *              label per pixel. note that labels values must be is range\r
  92  *              [0..num_labels-1].\r
  93  *\r
  94  *   - 'a': Perform alpha - beta swappings\r
  95  *           [gch labels] = GraphCutMex(gch, 'a')\r
  96  *           [gch labels] = GraphCutMex(gch, 'a', iter)\r
  97  *           [gch labels] = GraphCutMex(gch, 'a', label1, label2)\r
  98  *\r
  99  *       When no inputs are provided, GraphCut performs alpha - beta swaps steps\r
 100  *       until it converges.\r
 101  *\r
 102  *       Inputs:\r
 103  *           - iter: a double scalar, the maximum number of swap\r
 104  *                      iterations to perform.\r
 105  *           - label1, label2: int32 scalars denoting two labels for swap\r
 106  *                                       step.\r
 107  *\r
 108  *       Outputs:\r
 109  *           - labels: a width*height array of type int32, containing a\r
 110  *              label per pixel. note that labels values must be is range\r
 111  *              [0..num_labels-1].\r
 112  *\r
 113  *   - 'c': Close the graph and release allocated resources.\r
 114  *       [gch] = GraphCutMex(gch,'c');\r
 115  *\r
 116  *\r
 117  *   See Also:\r
 118  *       GraphCutConstr\r
 119  *\r
 120  *   This wrapper for Matlab was written by Shai Bagon (shai.bagon@weizmann.ac.il).\r
 121  *   Department of Computer Science and Applied Mathmatics\r
 122  *   Wiezmann Institute of Science\r
 123  *   http://www.wisdom.weizmann.ac.il/\r
 124  *\r
 125  *   The core cpp application was written by Veksler Olga:\r
 126  *        \r
 127  *   [1] Efficient Approximate Energy Minimization via Graph Cuts\r
 128  *       Yuri Boykov, Olga Veksler, Ramin Zabih,\r
 129  *       IEEE transactions on PAMI, vol. 20, no. 12, p. 1222-1239, November 2001.\r
 130  * \r
 131  *   [2] What Energy Functions can be Minimized via Graph Cuts?\r
 132  *       Vladimir Kolmogorov and Ramin Zabih.\r
 133  *       To appear in IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI).\r
 134  *       Earlier version appeared in European Conference on Computer Vision (ECCV), May 2002.\r
 135  * \r
 136  *   [3] An Experimental Comparison of Min-Cut/Max-Flow Algorithms\r
 137  *       for Energy Minimization in Vision.\r
 138  *       Yuri Boykov and Vladimir Kolmogorov.\r
 139  *       In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI),\r
 140  *       September 2004\r
 141  *\r
 142  *   [4] Matlab Wrapper for Graph Cut.\r
 143  *       Shai Bagon.\r
 144  *       in www.wisdom.weizmann.ac.il/~bagon, December 2006.\r
 145  *     \r
 146  *      This software can be used only for research purposes, you should cite\r
 147  *      the aforementioned papers in any resulting publication.\r
 148  *      If you wish to use this software (or the algorithms described in the aforementioned paper)\r
 149  *      for commercial purposes, you should be aware that there is a US patent: \r
 150  * \r
 151  *              R. Zabih, Y. Boykov, O. Veksler, \r
 152  *              "System and method for fast approximate energy minimization via graph cuts ", \r
 153  *              United Stated Patent 6,744,923, June 1, 2004\r
 154  *\r
 155  *\r
 156  *   The Software is provided "as is", without warranty of any kind.\r
 157  *\r
 158  *\r
 159  */\r
 160 void mexFunction(\r
 161     int           nlhs,         /* number of expected outputs */\r
 162     mxArray       *plhs[],      /* mxArray output pointer array */\r
 163     int           nrhs,         /* number of inputs */\r
 164     const mxArray         *prhs[]       /* mxArray input pointer array */\r
 165     )\r
 166 {\r
 167     /* building graph is only call without gch */\r
 168     GCoptimization *MyGraph = NULL;\r
 169     char mode ='\0';\r
 170     if ( nrhs < 2 ) {\r
 171         mexErrMsgIdAndTxt("GraphCut","Too few input arguments");\r
 172     }\r
 173 \r
 174     /* get graph handle */\r
 175     MyGraph = GetGCHandle(prhs[0]);\r
 176     \r
 177     /* get mode */\r
 178     GetScalar(prhs[1], mode);\r
 179     \r
 180    \r
 181     switch (mode) {\r
 182         case 'c': /* close */\r
 183             delete MyGraph; /* ->~GCoptimization(); /* explicitly call the destructor */\r
 184             MyGraph = NULL;\r
 185             break;\r
 186         case 's': /* set labels */\r
 187             /* setting the labels: we have an extra parameter - int32 array of size w*l\r
 188              * containing the new labels */\r
 189             SetLabels(MyGraph, nlhs, plhs, nrhs, prhs);\r
 190             break;\r
 191         case 'g': /* get labels */\r
 192             GetLabels(MyGraph, nlhs, plhs, nrhs, prhs);\r
 193             break;\r
 194         case 'a': /* a-b swaps */\r
 195             ABswaps(MyGraph, nlhs, plhs, nrhs, prhs);\r
 196             break;\r
 197         case 'e': /* expand steps */\r
 198             Expand(MyGraph, nlhs, plhs, nrhs, prhs);\r
 199             break;\r
 200         case 'n': /* get the current labeling energy */\r
 201             Energy(MyGraph, nlhs, plhs, nrhs, prhs);\r
 202             break;\r
 203         case 't': /* truncation */\r
 204             Truncate(MyGraph, nlhs, plhs, nrhs, prhs);\r
 205             break;\r
 206         default:\r
 207             mexErrMsgIdAndTxt("GraphCut:mode","unrecognized mode");\r
 208     }\r
 209     /* update the gch output handle */\r
 210     GraphHandle *pgh;\r
 211     plhs[0] = mxCreateNumericMatrix(1, 1, MATLAB_POINTER_TYPE, mxREAL);\r
 212     pgh = (GraphHandle*) mxGetData(plhs[0]);\r
 213     *pgh = (GraphHandle)(MyGraph);\r
 214 }\r
 215 /**************************************************************************************/\r
 216 /* set user defined labels to graph */\r
 217 void SetLabels(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\r
 218 {\r
 219     /* we need exactly 3 input arguments: gch, mode, labels */\r
 220     if (nrhs != 3 ) \r
 221         mexErrMsgIdAndTxt("GraphCut:SetLabels","Wrong number of input arguments");\r
 222     if ( mxGetClassID(prhs[2]) != mxINT32_CLASS ) \r
 223         mexErrMsgIdAndTxt("GraphCut:SetLabels","labels array not int32");\r
 224     if ( mxGetNumberOfElements(prhs[2]) != MyGraph->GetNumPixels() )\r
 225         mexErrMsgIdAndTxt("GraphCut:SetLabels","wrong number of elements in labels array");\r
 226     \r
 227     /* verify only one output parameter */\r
 228     if (nlhs != 1)\r
 229         mexErrMsgIdAndTxt("GraphCut:SetLabels","wrong number of output parameters");\r
 230     \r
 231     MyGraph->SetAllLabels( (GCoptimization::LabelType*) mxGetData(prhs[2]) );\r
 232 }\r
 233 /**************************************************************************************/\r
 234 /* set user defined labels to graph */\r
 235 void GetLabels(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\r
 236 {\r
 237     /* exactly two input arguments */\r
 238     if ( nrhs != 2 ) \r
 239         mexErrMsgIdAndTxt("GraphCut:GetLabels","Wrong number of input arguments");\r
 240         \r
 241     /* we need exactly 2 output arguments: gch, labels */\r
 242     if (nlhs != 2 ) \r
 243         mexErrMsgIdAndTxt("GraphCut:GetLabels","Wrong number of output arguments");\r
 244 \r
 245     /* transposed result */\r
 246     plhs[1] = mxCreateNumericMatrix(MyGraph->GetWidth(), MyGraph->GetNumPixels() / MyGraph->GetWidth(), mxINT32_CLASS, mxREAL);\r
 247     MyGraph->ExportLabels( (GCoptimization::LabelType*) mxGetData(plhs[1]) );\r
 248 }\r
 249 /**************************************************************************************/\r
 250 /* support several kinds of a/b swaps:\r
 251  * 1. defualt - swap till convergence (no extra parameters)\r
 252  * 2. #iterations - performs #iterations:   GraphCut(gch, mode, #iteration)\r
 253  * 3. swap two labels - GraphCut(gch, mode, l1, l2)\r
 254  */\r
 255 void ABswaps(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\r
 256 {\r
 257     int max_iterations(0), li(0);\r
 258     GCoptimization::LabelType label[2] = { 0, 1};\r
 259     GCoptimization::PixelType *indices = NULL;\r
 260     \r
 261     /* check inputs */\r
 262     switch (nrhs) {\r
 263         case 2:\r
 264             /* default - expand till convergence */\r
 265             MyGraph->swap();\r
 266             break;\r
 267         case 3:\r
 268             /* number of iterations */\r
 269             GetScalar(prhs[2], max_iterations);\r
 270             if ( max_iterations == 1 ) \r
 271                 MyGraph->oneSwapIteration();\r
 272             else\r
 273                 MyGraph->swap(max_iterations);\r
 274             break;\r
 275         case 4:\r
 276             for ( li = 2; li<4;li++) {\r
 277                 GetScalar(prhs[li], label[li-2]);\r
 278                 if ( label[li-2] < 0 || label[li-2] >= MyGraph->GetNumLabels() ) \r
 279                     mexErrMsgIdAndTxt("GraphCut:swap","invalid label value");\r
 280             }\r
 281             MyGraph->alpha_beta_swap(label[0], label[1]);\r
 282             break;\r
 283         default:\r
 284             mexErrMsgIdAndTxt("GraphCut:swap","Too many input arguments");\r
 285     }\r
 286         \r
 287     /* output: at least one (gch) can output the labels as well */\r
 288     if ( nlhs > 2 ) \r
 289         mexErrMsgIdAndTxt("GraphCut:swap","too many output arguments");\r
 290 \r
 291     /* output labels */\r
 292     if ( nlhs >= 2 ) {\r
 293         /* set the lables */\r
 294         \r
 295         /* transposed result */\r
 296         plhs[1] = mxCreateNumericMatrix(MyGraph->GetWidth(), MyGraph->GetNumPixels() / MyGraph->GetWidth(), mxINT32_CLASS, mxREAL);\r
 297         MyGraph->ExportLabels( (GCoptimization::LabelType*) mxGetData(plhs[1]) );\r
 298     }\r
 299                     \r
 300 }\r
 301 \r
 302 /**************************************************************************************/\r
 303 /* support several kinds of expansion \r
 304  * 1. default - expand untill convergence no extra parameters\r
 305  * 2. #iterations - performs #iterations:   GraphCut(gch, mode, #iteration)\r
 306  * 3. expand label - expand specific label  GraphCut(gch, mode, [], label)\r
 307  * 4. expand label at specific indices      GraphCut(gch, mode, [], label, indices) indices start with zero and not 1 like usuall matlab inices!!\r
 308  */\r
 309 void Expand(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\r
 310 {   \r
 311     int num(0), max_iterations(0);\r
 312     GCoptimization::LabelType label(0);\r
 313     GCoptimization::PixelType *indices = NULL;\r
 314     \r
 315     /* check inputs */\r
 316     switch (nrhs) {\r
 317         case 2:\r
 318             /* default - expand till convergence */\r
 319             MyGraph->expansion();\r
 320             break;\r
 321         case 3:\r
 322             /* number of iterations */\r
 323             GetScalar(prhs[2], max_iterations);\r
 324             if ( max_iterations == 1 ) \r
 325                 MyGraph->oneExpansionIteration();\r
 326             else\r
 327                 MyGraph->expansion(max_iterations);\r
 328             break;\r
 329         case 5:\r
 330             /* get indices */\r
 331             if (mxGetClassID(prhs[4]) != mxINT32_CLASS)\r
 332                 mexErrMsgIdAndTxt("GraphCut:Expand","indices must be int32");\r
 333             num = mxGetNumberOfElements(prhs[4]);\r
 334             if (num < 0 || num > MyGraph->GetNumPixels())\r
 335                 mexErrMsgIdAndTxt("GraphCut:Expand","too many indices");\r
 336             indices = (GCoptimization::PixelType*)mxGetData(prhs[4]);\r
 337         case 4:\r
 338             /* expand specific label */\r
 339             if (mxGetNumberOfElements(prhs[2]) != 0)\r
 340                 mexErrMsgIdAndTxt("GraphCut:Expand","third argument must empty");\r
 341             GetScalar(prhs[3], label);\r
 342 \r
 343             if ( indices != NULL ) {\r
 344                 /* expand specific label at specific indices */\r
 345                 MyGraph->alpha_expansion(label, indices, num);\r
 346             } else\r
 347                 MyGraph->alpha_expansion(label);\r
 348 \r
 349             break;\r
 350         default:\r
 351             mexErrMsgIdAndTxt("GraphCut:Expand","Too many input arguments");\r
 352     }\r
 353         \r
 354     /* output: at least one (gch) can output the labels as well */\r
 355     if ( nlhs > 2 ) \r
 356         mexErrMsgIdAndTxt("GraphCut:Expand","too many output arguments");\r
 357 \r
 358     /* output labels */\r
 359     if ( nlhs >= 2 ) {\r
 360         /* set the lables */\r
 361         \r
 362         /* transposed result */        \r
 363         plhs[1] = mxCreateNumericMatrix(MyGraph->GetWidth(), MyGraph->GetNumPixels() / MyGraph->GetWidth(), mxINT32_CLASS, mxREAL);\r
 364         MyGraph->ExportLabels( (GCoptimization::LabelType*) mxGetData(plhs[1]) );\r
 365 \r
 366     }\r
 367 }\r
 368 \r
 369 /**************************************************************************************/\r
 370 /* extract energy of labeling  \r
 371  * [gch se de] = GraphCut(gch, 'n');\r
 372  */\r
 373 void Energy(GCoptimization *MyGraph, int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])\r
 374 {\r
 375     if ( nrhs != 2 )\r
 376         mexErrMsgIdAndTxt("GraphCut:Energy","Too many input arguments");\r
 377     if ( nlhs != 3 )\r
 378         mexErrMsgIdAndTxt("GraphCut:Energy","wrong number of output arguments");\r
 379     \r
 380     GCoptimization::EnergyType *e;\r
 381     \r
 382     plhs[1] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);\r
 383     e = (GCoptimization::EnergyType *)mxGetData(plhs[1]);\r
 384     *e = MyGraph->giveSmoothEnergy();\r
 385     \r
 386     plhs[2] = mxCreateNumericMatrix(1, 1, mxSINGLE_CLASS, mxREAL);\r
 387     e = (GCoptimization::EnergyType *)mxGetData(plhs[2]);\r
 388     *e = MyGraph->giveDataEnergy();\r
 389 }\r
 390 /**************************************************************************************/\r
 391 /* get/set truncation flag\r
 392  * [gch tf] = GraphCut(gch, 'truncate', [ntf]) \r
 393  */\r
 394 void Truncate(GCoptimization *MyGraph, int nout, mxArray *pout[], int nin, const mxArray *pin[])\r
 395 {\r
 396     bool tf;\r
 397     if ( nin == 2 ) {\r
 398         // get the state\r
 399         tf = MyGraph->GetTruncateState();\r
 400         \r
 401     } else if ( nin == 3 ) {\r
 402         // set the state        \r
 403         if ( mxIsComplex(pin[2]) || mxIsSparse(pin[2]) || mxGetNumberOfElements(pin[2])!=1 )\r
 404             mexErrMsgIdAndTxt("GraphCut:Truncate","wrong truncate_flag");\r
 405         \r
 406         tf = (mxGetScalar(pin[2]) != 0);\r
 407         \r
 408         mexPrintf("Set truncate %d\n",tf);\r
 409         \r
 410         MyGraph->SetTruncateState(tf);\r
 411     } else {\r
 412         mexErrMsgIdAndTxt("GraphCut:Truncate","wrong number of input arguments");\r
 413     }\r
 414     if ( nout != 2 )\r
 415         mexErrMsgIdAndTxt("GraphCut:Truncate","wrong number of output arguments");\r
 416     \r
 417     pout[1] = mxCreateDoubleScalar( static_cast<double>(tf) );\r
 418 }\r
 419 \r
 420 /**************************************************************************************/\r
 421 GCoptimization* GetGCHandle(const mxArray *x)\r
 422 {\r
 423     GraphHandle gch = 0;\r
 424     GCoptimization* MyGraph = 0;\r
 425     \r
 426     if ( mxGetNumberOfElements(x) != 1 ) {\r
 427         mexErrMsgIdAndTxt("GraphCut:handle",\r
 428         "Too many GC handles");\r
 429     }\r
 430     if ( mxGetClassID(x) != MATLAB_POINTER_TYPE ) {\r
 431         mexErrMsgIdAndTxt("GraphCut:handle",\r
 432         "GC handle argument is not of proper type");\r
 433     }\r
 434     gch = (GraphHandle*)mxGetData(x);\r
 435     MyGraph = (GCoptimization*)(*(POINTER_CAST*)gch);\r
 436     if ( MyGraph==NULL || ! MyGraph->IsClassValid() ) {\r
 437         mexErrMsgIdAndTxt("GraphCut:handle",\r
 438         "GC handle is not valid");\r
 439     }\r
 440 \r
 441     \r
 442     return MyGraph;\r
 443 }\r
 444 /**************************************************************************************/\r
 445 \r