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

   1 \r
   2 #include "mex.h"\r
   3 #include "GCoptimization.h"\r
   4 #include "GraphCut.h"\r
   5 #include <stdlib.h>\r
   6 #include <math.h>\r
   7 \r
   8 /* Defines */\r
   9 // define A64BITS for compiling on a 64 bits machine...\r
  10 \r
  11 /*\r
  12  * Matlab wrapper for Weksler graph cut implementation\r
  13  *\r
  14  * usage:\r
  15  * [gch] = GraphCut3dConstr(R, C, Z, num_labels, DataCost, SmoothnessCost, [Contrast])\r
  16  *\r
  17  * Note that data types are crucials!\r
  18  * \r
  19  * Inputs:\r
  20  *  R, C, Z - size of 3D grid. \r
  21  *  num_labels - number of labels.  \r
  22  *  DataCost - of type float, array size [width*height*depth*#labels], the data term for pixel\r
  23  *             x,y,z recieving label l is stroed at [(x+y*width+z*width*depth)*#labels + l]\r
  24  *  SmoothnessCost - of type float, array size [#labels.^2] the cost between l1 and l2 is\r
  25  *                   stored at [l1+l2*#labels] = Vpq(lp,lq)\r
  26  *  Contrast - of type float, array size[width*height*depth], the weight Wpq will be determined by the contrast:\r
  27  *                  Wpq = exp(-|C(p)-C(q)|)\r
  28  *\r
  29  * Outputs:\r
  30  *  gch - of type int32, graph cut handle - do NOT mess with it!\r
  31  */\r
  32 void mexFunction(\r
  33     int           nlhs,         /* number of expected outputs */\r
  34     mxArray       *plhs[],      /* mxArray output pointer array */\r
  35     int           nrhs,         /* number of inputs */\r
  36     const mxArray         *prhs[]       /* mxArray input pointer array */\r
  37     )\r
  38 {\r
  39     \r
  40     int num_labels;\r
  41     Graph::captype *DataCost;\r
  42     Graph::captype *SmoothnessCost;\r
  43     Graph::captype *Contrast = NULL;\r
  44     GCoptimization::LabelType *Labels;\r
  45     GCoptimization::PixelType R, C, Z; \r
  46     \r
  47     GCoptimization *MyGraph = NULL;\r
  48         \r
  49     /* check number of inout arguments - must be 6 */\r
  50     if ((nrhs != 6)&&(nrhs!=7)) {\r
  51         mexErrMsgIdAndTxt("GraphCut:NarginError","Wrong number of input argumnets");\r
  52     }\r
  53     GetScalar(prhs[0], R);\r
  54     GetScalar(prhs[1], C);\r
  55     GetScalar(prhs[2], Z);\r
  56     GetScalar(prhs[3], num_labels);\r
  57     \r
  58   \r
  59     /* check second input DataCost: must have #labels*height*width elements of type float */\r
  60     if ( mxGetNumberOfElements(prhs[4]) != num_labels*R*C*Z ) {\r
  61         mexErrMsgIdAndTxt("GraphCut:DataCost",\r
  62         "DataCost argument does not contains the right number of elements");\r
  63     }\r
  64     if (mxGetClassID(prhs[4]) != mxSINGLE_CLASS ) {\r
  65         mexErrMsgIdAndTxt("GraphCut:DataCost",\r
  66         "DataCost argument is not of type float");\r
  67     }\r
  68     DataCost = (Graph::captype*)mxGetData(prhs[4]);\r
  69 \r
  70     /* check fifth input SmoothnessCost: must have #labels.^2 elements of type float */\r
  71     if ( mxGetNumberOfElements(prhs[5]) != num_labels*num_labels ) {\r
  72         mexErrMsgIdAndTxt("GraphCut:SmoothnessCost",\r
  73         "SmoothnessCost argument does not contains the right number of elements");\r
  74     }\r
  75     if (mxGetClassID(prhs[5]) != mxSINGLE_CLASS ) {\r
  76         mexErrMsgIdAndTxt("GraphCut:SmoothnessCost",\r
  77         "SmoothnessCost argument is not of type float");\r
  78     }\r
  79     SmoothnessCost = (Graph::captype*)mxGetData(prhs[5]);\r
  80 \r
  81     if ( nrhs == 7 ) { \r
  82         /* add Contrast cue */\r
  83         if ( mxGetNumberOfElements(prhs[6]) != R*C*Z ) {\r
  84             mexErrMsgIdAndTxt("GraphCut:SmoothnessCost",\r
  85             "Contrast argument does not contains the right number of elements");\r
  86         }\r
  87         if (mxGetClassID(prhs[6]) != mxSINGLE_CLASS ) {\r
  88             mexErrMsgIdAndTxt("GraphCut:SmoothnessCost",\r
  89             "Contrast argument is not of type float");\r
  90         }\r
  91         Contrast = (Graph::captype*)mxGetData(prhs[6]);\r
  92         \r
  93     }\r
  94     /* prepare the output argument */\r
  95     if ( nlhs != 1 ) {\r
  96         mexErrMsgIdAndTxt("GraphCut:OutputArg","Wrong number of output arguments");\r
  97     }\r
  98     \r
  99     MyGraph = new GCoptimization(R*C*Z, num_labels, SET_ALL_AT_ONCE, SET_ALL_AT_ONCE);\r
 100     /* neighborhod setup */\r
 101     GCoptimization::PixelType c(0), r(0), z(0), p(0), q(0);\r
 102     if (Contrast) {\r
 103         for ( r = 0 ; r <= R - 2;  r++ ) {\r
 104             for ( c = 0 ; c <= C - 2; c++ ) {\r
 105                 for ( z = 0 ; z <= Z - 2; z++ ) {\r
 106                     /* all nodes with 3 nieghbors */\r
 107                     p = r+c*R+z*R*C;\r
 108                     q = r+1+c*R+z*R*C;\r
 109                     MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 110                     q = r+(c+1)*R+z*R*C;\r
 111                     MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 112                     q = r+c*R+(z+1)*R*C;\r
 113                     MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 114                 }\r
 115                /* top of Z has 2 n in c and r */\r
 116                 p = r+c*R+z*R*C;\r
 117                 q = r+1+c*R+z*R*C;\r
 118                 MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 119                 q = r+(c+1)*R+z*R*C;\r
 120                 MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 121             }\r
 122             /* end of c has 2 n in z and r */\r
 123             for ( z = 0 ; z <= Z -2 ; z++ ) {\r
 124                 p = r+c*R+z*R*C;\r
 125                 q = r+1+c*R+z*R*C;\r
 126                 MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 127                 q = r+c*R+(z+1)*R*C;\r
 128                 MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 129             }\r
 130             /* end of c abd z has n in r */\r
 131             p = r+c*R+z*R*C;\r
 132             q = r+1+c*R+z*R*C;\r
 133             MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 134         }\r
 135         /* end of r has n in z and c */\r
 136         for ( c = 0 ; c <= C - 2; c++ ) {\r
 137             for ( z = 0 ; z <= Z - 2; z++ ) {\r
 138                 p = r+c*R+z*R*C;\r
 139                 q = r+c*R+(z+1)*R*C;\r
 140                 MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 141                 q = r+(c+1)*R+z*R*C;\r
 142                 MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 143             }\r
 144             p = r+c*R+z*R*C;\r
 145             q = r+(c+1)*R+z*R*C;\r
 146             MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 147         }\r
 148         for ( z = 0 ; z <= Z - 2; z++ ) {\r
 149             p = r+c*R+z*R*C;\r
 150             q = r+c*R+(z+1)*R*C;\r
 151             MyGraph->setNeighbors(p, q, exp(-fabs(Contrast[p]-Contrast[q])));\r
 152         }\r
 153         /* end of graph construction with contrast */\r
 154     } else {\r
 155         for ( r = 0 ; r <= R - 2;  r++ ) {\r
 156             for ( c = 0 ; c <= C - 2; c++ ) {\r
 157                 for ( z = 0 ; z <= Z - 2; z++ ) {\r
 158                 /* all nodes with 3 nieghbors */\r
 159                     MyGraph->setNeighbors(r+c*R+z*R*C,r+1+c*R+z*R*C);\r
 160                     MyGraph->setNeighbors(r+c*R+z*R*C,r+(c+1)*R+z*R*C);\r
 161                     MyGraph->setNeighbors(r+c*R+z*R*C,r+c*R+(z+1)*R*C);\r
 162                 }\r
 163             /* top of Z has 2 n in c and r */\r
 164                 MyGraph->setNeighbors(r+c*R+z*R*C,r+1+c*R+z*R*C);\r
 165                 MyGraph->setNeighbors(r+c*R+z*R*C,r+(c+1)*R+z*R*C);\r
 166             }\r
 167         /* end of c has 2 n in z and r */\r
 168             for ( z = 0 ; z <= Z -2 ; z++ ) {\r
 169                 MyGraph->setNeighbors(r+c*R+z*R*C,r+1+c*R+z*R*C);\r
 170                 MyGraph->setNeighbors(r+c*R+z*R*C,r+c*R+(z+1)*R*C);\r
 171             }\r
 172         /* end of c abd z has n in r */\r
 173             MyGraph->setNeighbors(r+c*R+z*R*C,r+1+c*R+z*R*C);\r
 174         }\r
 175     /* end of r has n in z and c */\r
 176         for ( c = 0 ; c <= C - 2; c++ ) {\r
 177             for ( z = 0 ; z <= Z - 2; z++ ) {\r
 178                 MyGraph->setNeighbors(r+c*R+z*R*C,r+c*R+(z+1)*R*C);\r
 179                 MyGraph->setNeighbors(r+c*R+z*R*C,r+(c+1)*R+z*R*C);\r
 180             }\r
 181             MyGraph->setNeighbors(r+c*R+z*R*C,r+(c+1)*R+z*R*C);\r
 182         }\r
 183         for ( z = 0 ; z <= Z - 2; z++ ) {\r
 184             MyGraph->setNeighbors(r+c*R+z*R*C,r+c*R+(z+1)*R*C);\r
 185         }\r
 186     }\r
 187     MyGraph->setData(DataCost);\r
 188 /*    if ( vCue != NULL && vCue != NULL ) \r
 189         MyGraph->setSmoothness(SmoothnessCost, hCue, vCue);\r
 190     else*/\r
 191     MyGraph->setSmoothness(SmoothnessCost);\r
 192     \r
 193         \r
 194     /* create a container for the pointer */\r
 195     plhs[0] = mxCreateNumericMatrix(1, 1, MATLAB_POINTER_TYPE, mxREAL);\r
 196     \r
 197     GraphHandle* gh;\r
 198     gh = (GraphHandle*) mxGetData(plhs[0]);\r
 199     *gh = (GraphHandle)(MyGraph);\r
 200 }\r
 201     \r
 202 \r