THESE/codes/graphe/GCmex1.9/GraphCutConstrSparse.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 /* Defines */\r
   8 \r
   9 \r
  10 /*\r
  11  * Matlab wrapper for Weksler graph cut implementation\r
  12  *\r
  13  * usage:\r
  14  * [gch] = GraphCutConstrSparse(dc, sc, SparseSc)\r
  15  *\r
  16  * Note that data types are crucials!\r
  17  * \r
  18  * Inputs:\r
  19  *  dc - of type float, array size [#labels*#nodes], the data term for node\r
  20  *             n recieving label l is stroed at [n*#labels + l]\r
  21  *  sc - of type float, array size [#labels.^2] the cost between l1 and l2 is\r
  22  *                   stored at [l1+l2*#labels] = Vpq(lp,lq)\r
  23  *  SparseSc - Sparse matrix of type double defining both the graph structure \r
  24  *              and spatially dependent smoothness term\r
  25  *\r
  26  * Outputs:\r
  27  *  gch - of type int32, graph cut handle - do NOT mess with it!\r
  28  */\r
  29 \r
  30 \r
  31 void mexFunction(\r
  32     int           nlhs,         /* number of expected outputs */\r
  33     mxArray       *plhs[],      /* mxArray output pointer array */\r
  34     int           nrhs,         /* number of inputs */\r
  35     const mxArray         *prhs[]       /* mxArray input pointer array */\r
  36     )\r
  37 {\r
  38     /* check number of arguments */\r
  39     if (nrhs != 3) {\r
  40         mexErrMsgIdAndTxt("GraphCut:NarginError","Wrong number of input argumnets");\r
  41     }\r
  42     if (nlhs != 1) {\r
  43         mexErrMsgIdAndTxt("GraphCut:NargoutError","Wrong number of output argumnets");\r
  44     }\r
  45     /* check inputs */\r
  46     if (mxGetClassID(prhs[0]) != mxSINGLE_CLASS ) {\r
  47         mexErrMsgIdAndTxt("GraphCut:DataCost", "DataCost argument is not of type float");\r
  48     }\r
  49     if (mxGetClassID(prhs[1]) != mxSINGLE_CLASS ) {\r
  50         mexErrMsgIdAndTxt("GraphCut:SmoothnessCost", "SmoothnessCost argument is not of type float");\r
  51     }\r
  52     if (! mxIsSparse(prhs[2]) ) {\r
  53         mexErrMsgIdAndTxt("GraphCut:SmoothnessCost", "Graph Structure Matrix is not sparse");\r
  54     }\r
  55     \r
  56     GCoptimization::PixelType num_pixels;\r
  57     int num_labels;\r
  58     \r
  59     num_pixels = mxGetN(prhs[2]);\r
  60     if (mxGetM(prhs[2]) != num_pixels) {\r
  61         mexErrMsgIdAndTxt("GraphCut:SmoothnessCost", "Graph Structure Matrix is no square");\r
  62     }\r
  63     num_labels = mxGetNumberOfElements(prhs[0])/num_pixels;\r
  64     if (mxGetNumberOfElements(prhs[1]) != num_labels*num_labels) {\r
  65         mexErrMsgIdAndTxt("GraphCut:SmoothnessCost", "Size does not match number of labels");\r
  66     }\r
  67     \r
  68     /* construct the graph */\r
  69     GCoptimization* MyGraph = new GCoptimization(num_pixels, num_labels, SET_ALL_AT_ONCE, SET_ALL_AT_ONCE);\r
  70     \r
  71     /* set the nieghbors and weights according to sparse matrix */\r
  72     \r
  73     double   *pr;\r
  74     mwIndex  *ir, *jc;\r
  75     mwSize   col, total=0;\r
  76     mwIndex  starting_row_index, stopping_row_index, current_row_index;\r
  77 \r
  78   \r
  79     /* Get the starting positions of all four data arrays. */\r
  80     pr = mxGetPr(prhs[2]);\r
  81     ir = mxGetIr(prhs[2]);\r
  82     jc = mxGetJc(prhs[2]);\r
  83     \r
  84     for (col=0; col<num_pixels; col++)  {\r
  85         starting_row_index = jc[col];\r
  86         stopping_row_index = jc[col+1];\r
  87         if (starting_row_index == stopping_row_index) {\r
  88             continue;\r
  89         } else {\r
  90             for (current_row_index = starting_row_index;\r
  91                 current_row_index < stopping_row_index;\r
  92                 current_row_index++)  {\r
  93                     /* use only upper triangle of matrix */\r
  94                     if ( ir[current_row_index] >= col ) {\r
  95                         MyGraph->setNeighbors(ir[current_row_index], col, pr[total++]);\r
  96                     } else {\r
  97                         total++;\r
  98                     }\r
  99                     \r
 100             }\r
 101         }\r
 102     }\r
 103     Graph::captype *DataCost = (Graph::captype*)mxGetData(prhs[0]);\r
 104     Graph::captype *SmoothnessCost = (Graph::captype*)mxGetData(prhs[1]);\r
 105     \r
 106     /* set data term */\r
 107     MyGraph->setData(DataCost);\r
 108     /* set the smoothness term */\r
 109     MyGraph->setSmoothness(SmoothnessCost);\r
 110     \r
 111         \r
 112     /* create a container for the pointer */\r
 113     const mwSize dims[2] = {1,0};\r
 114     plhs[0] = mxCreateNumericArray(1, /*(int*)*/dims, MATLAB_POINTER_TYPE, mxREAL);\r
 115     \r
 116     GraphHandle* gh;\r
 117     gh = (GraphHandle*) mxGetData(plhs[0]);\r
 118     *gh = (GraphHandle)(MyGraph);\r
 119 }\r