THESE/codes/graphe/GCmex1.9/GraphCut.m

   1 function [gch, varargout] = GraphCut(mode, varargin)\r
   2 %\r
   3 %   Performing Graph Cut energy minimization operations on a 2D grid.\r
   4 %   \r
   5 %   Usage:\r
   6 %       [gch ...] = GraphCut(mode, ...);\r
   7 %   \r
   8 %\r
   9 %   Inputs:\r
  10 %   - mode: a string specifying mode of operation. See details below.\r
  11 %\r
  12 %   Output:\r
  13 %   - gch: A handle to the constructed graph. Handle this handle with care\r
  14 %          and don't forget to close it in the end!\r
  15 %\r
  16 %   Possible modes:\r
  17 %   - 'open': Create a new graph object\r
  18 %           [gch] = GraphCut('open', DataCost, SmoothnessCost);\r
  19 %           [gch] = GraphCut('open', DataCost, SmoothnessCost, vC, hC);\r
  20 %           [gch] = GraphCut('open', DataCost, SmoothnessCost, SparseSmoothness);\r
  21 %\r
  22 %       Inputs:\r
  23 %           - DataCost a height by width by num_labels matrix where\r
  24 %             Dc(r,c,l) equals the cost for assigning label l to pixel at (r,c)\r
  25 %             Note that the graph dimensions, and the number of labels are deduced \r
  26 %             form the size of the DataCost matrix.\r
  27 %             When using SparseSmoothness Dc is of (L)x(P) where L is the\r
  28 %             number of labels and P is the number of nodes/pixels in the\r
  29 %             graph.\r
  30 %           - SmoothnessCost a #labels by #labels matrix where Sc(l1, l2)\r
  31 %             is the cost of assigning neighboring pixels with label1 and\r
  32 %             label2. This cost is spatialy invariant\r
  33 %           - vC, hC:optional arrays defining spatialy varying smoothness cost. \r
  34 %                       Single precission arrays of size width*height.\r
  35 %                       The smoothness cost is computed using:\r
  36 %                       V_pq(l1, l2) = V(l1, l2) * w_pq\r
  37 %                       where V is the SmoothnessCost matrix\r
  38 %                       w_pq is spatialy varying parameter:\r
  39 %                       if p=(r,c) and q=(r+1,c) then w_pq = vCue(r,c)\r
  40 %                       if p=(r,c) and q=(r,c+1) then w_pq = hCue(r,c)\r
  41 %                       (therefore in practice the last column of vC and\r
  42 %                       the last row of vC are not used).\r
  43 %           - SparseSmoothness: a sparse matrix defining both the graph\r
  44 %               structure (might be other than grid) and the spatialy varying\r
  45 %               smoothness term. Must be real positive sparse matrix of size\r
  46 %               num_pixels by num_pixels, each non zero entry (i,j) defines a link\r
  47 %               between pixels i and j with w_pq = SparseSmoothness(i,j).\r
  48 %\r
  49 %   - 'set': Set labels\r
  50 %           [gch] = GraphCut('set', gch, labels)\r
  51 %\r
  52 %       Inputs:\r
  53 %           - labels: a width by height array containing a label per pixel.\r
  54 %             Array should be the same size of the grid with values\r
  55 %             [0..num_labels].\r
  56 %\r
  57 %\r
  58 %   - 'get': Get current labeling\r
  59 %           [gch labels] = GraphCut('get', gch)\r
  60 %\r
  61 %       Outputs:\r
  62 %           - labels: a height by width array, containing a label per pixel. \r
  63 %             note that labels values are in range [0..num_labels-1].\r
  64 %\r
  65 %\r
  66 %   - 'energy': Get current values of energy terms\r
  67 %           [gch se de] = GraphCut('energy', gch)\r
  68 %           [gch e] = GraphCut('energy', gch)\r
  69 %\r
  70 %       Outputs:\r
  71 %           - se: Smoothness energy term.\r
  72 %           - de: Data energy term.\r
  73 %           - e = se + de\r
  74 %\r
  75 %\r
  76 %   - 'expand': Perform labels expansion\r
  77 %           [gch labels] = GraphCut('expand', gch)\r
  78 %           [gch labels] = GraphCut('expand', gch, iter)\r
  79 %           [gch labels] = GraphCut('expand', gch, [], label)\r
  80 %           [gch labels] = GraphCut('expand', gch, [], label, indices)\r
  81 %\r
  82 %       When no inputs are provided, GraphCut performs expansion steps\r
  83 %       until it converges.\r
  84 %\r
  85 %       Inputs:\r
  86 %           - iter: a double scalar, the maximum number of expand\r
  87 %                    iterations to perform.\r
  88 %           - label: scalar denoting the label for which to perfom\r
  89 %                     expand step (labels are [0..num_labels-1]).\r
  90 %           - indices: array of linear indices of pixels for which\r
  91 %                        expand step is computed. \r
  92 %\r
  93 %       Outputs:\r
  94 %           - labels: a width*height array of type int32, containing a\r
  95 %              label per pixel. note that labels values must be is range\r
  96 %              [0..num_labels-1].\r
  97 %\r
  98 %\r
  99 %   - 'swap': Perform alpha - beta swappings\r
 100 %           [gch labels] = GraphCut('swap', gch)\r
 101 %           [gch labels] = GraphCut('swap', gch, iter)\r
 102 %           [gch labels] = GraphCut('swap', gch, label1, label2)\r
 103 %\r
 104 %       When no inputs are provided, GraphCut performs alpha - beta swaps steps\r
 105 %       until it converges.\r
 106 %\r
 107 %       Inputs:\r
 108 %           - iter: a double scalar, the maximum number of swap\r
 109 %                      iterations to perform.\r
 110 %           - label1, label2: int32 scalars denoting two labels for swap\r
 111 %                                       step.\r
 112 %\r
 113 %       Outputs:\r
 114 %           - labels: a width*height array of type int32, containing a\r
 115 %              label per pixel. note that labels values must be is range\r
 116 %              [0..num_labels-1].\r
 117 %\r
 118 %\r
 119 %   - 'truncate': truncating (or not) violating expansion terms \r
 120 %                 (see Rother etal. Digital Tapestry, CVPR2005)\r
 121 %       [gch truncate_flag] = GraphCut('truncate', gch, trancate_flag);\r
 122 %   \r
 123 %       When no truncate_flag is provided the function returns the current\r
 124 %       state of truncation\r
 125 %       \r
 126 %       Inputs:\r
 127 %           - trancate_flag: set truncate_flag to this state\r
 128 %\r
 129 %       Outputs:\r
 130 %           - trancate_flag: current state (after modification if\r
 131 %                            applicable)\r
 132 %\r
 133 %   - 'close': Close the graph and release allocated resources.\r
 134 %       [gch] = GraphCut('close', gch);\r
 135 %\r
 136 %\r
 137 %\r
 138 %   This wrapper for Matlab was written by Shai Bagon (shai.bagon@weizmann.ac.il).\r
 139 %   Department of Computer Science and Applied Mathmatics\r
 140 %   Wiezmann Institute of Science\r
 141 %   http://www.wisdom.weizmann.ac.il/\r
 142 %\r
 143 %       The core cpp application was written by Olga Veksler\r
 144 %       (available from http://www.csd.uwo.ca/faculty/olga/code.html):\r
 145 %\r
 146 %   [1] Efficient Approximate Energy Minimization via Graph Cuts\r
 147 %        Yuri Boykov, Olga Veksler, Ramin Zabih,\r
 148 %        IEEE transactions on PAMI, vol. 20, no. 12, p. 1222-1239, November\r
 149 %        2001.\r
 150 % \r
 151 %   [2] What Energy Functions can be Minimized via Graph Cuts?\r
 152 %        Vladimir Kolmogorov and Ramin Zabih.\r
 153 %        IEEE Transactions on Pattern Analysis and Machine Intelligence\r
 154 %        (PAMI), vol. 26, no. 2,\r
 155 %        February 2004, pp. 147-159.\r
 156 % \r
 157 %   [3] An Experimental Comparison of Min-Cut/Max-Flow Algorithms\r
 158 %        for Energy Minimization in Vision.\r
 159 %        Yuri Boykov and Vladimir Kolmogorov.\r
 160 %        In IEEE Transactions on Pattern Analysis and Machine Intelligence\r
 161 %        (PAMI),\r
 162 %        vol. 26, no. 9, September 2004, pp. 1124-1137.\r
 163 % \r
 164 %   [4] Matlab Wrapper for Graph Cut.\r
 165 %        Shai Bagon.\r
 166 %        in www.wisdom.weizmann.ac.il/~bagon, December 2006.\r
 167 % \r
 168 %   This software can be used only for research purposes, you should  cite ALL of\r
 169 %   the aforementioned papers in any resulting publication.\r
 170 %   If you wish to use this software (or the algorithms described in the\r
 171 %   aforementioned paper)\r
 172 %   for commercial purposes, you should be aware that there is a US patent:\r
 173 %\r
 174 %       R. Zabih, Y. Boykov, O. Veksler,\r
 175 %       "System and method for fast approximate energy minimization via\r
 176 %       graph cuts ",\r
 177 %       United Stated Patent 6,744,923, June 1, 2004\r
 178 %\r
 179 %\r
 180 %   The Software is provided "as is", without warranty of any kind.\r
 181 %\r
 182 %\r
 183 \r
 184 switch lower(mode)\r
 185     case {'o', 'open'}\r
 186         % open a new graph cut\r
 187         if nargout ~= 1\r
 188             error('GraphCut:Open: wrong number of output arguments');\r
 189         end\r
 190 \r
 191         gch = OpenGraph(varargin{:});\r
 192 \r
 193     case {'c', 'close'}\r
 194         % close the GraphCut handle - free memory.\r
 195         if nargin ~= 2\r
 196             error('GraphCut:Close: Too many inputs');\r
 197         end\r
 198         gch = varargin{1};\r
 199         [gch] = GraphCutMex(gch, 'c');\r
 200         \r
 201     case {'g', 'get'}\r
 202         % get current labeling\r
 203         \r
 204         if nargout ~= 2\r
 205             error('GraphCut:GetLabels: wrong number of outputs');\r
 206         end\r
 207         [gch labels] = GetLabels(varargin{:});\r
 208 \r
 209         varargout{1} = labels;\r
 210         \r
 211     case {'s', 'set'}\r
 212         % set user defined labeling\r
 213         if nargout ~= 1\r
 214             error('GraphCut:SetLabels: Too many outputs');\r
 215         end\r
 216         \r
 217         [gch] = SetLabels(varargin{:});\r
 218         \r
 219     case {'en', 'n', 'energy'}\r
 220         % get current energy values\r
 221         if nargin ~= 2\r
 222             error('GraphCut:GetEnergy: too many input arguments');\r
 223         end\r
 224         gch = varargin{1};\r
 225         [gch se de] = GraphCutMex(gch, 'n');\r
 226         switch nargout\r
 227             case 2\r
 228                 varargout{1} = se+de;\r
 229             case 3\r
 230                 varargout{1} = se;\r
 231                 varargout{2} = de;\r
 232             case 1\r
 233             otherwise\r
 234                 error('GraphCut:GetEnergy: wrong number of output arguments')\r
 235         end\r
 236    \r
 237     \r
 238     case {'e', 'ex', 'expand'}\r
 239         % use expand steps to minimize energy\r
 240         if nargout > 2\r
 241             error('GraphCut:Expand: too many output arguments');\r
 242         end\r
 243         [gch labels] = Expand(varargin{:});\r
 244         if nargout == 2\r
 245             varargout{1} = labels;\r
 246         end\r
 247         \r
 248     case {'sw', 'a', 'ab', 'swap'}\r
 249         % use alpha beta swapping to minimize energy\r
 250         if nargout > 2\r
 251             error('GraphCut:Expand: too many output arguments');\r
 252         end\r
 253         [gch labels] = Swap(varargin{:});\r
 254         if nargout == 2\r
 255             varargout{1} = labels;\r
 256         end\r
 257         \r
 258     case {'truncate'}\r
 259         \r
 260         if numel(varargin) == 2\r
 261             gch = varargin{1};\r
 262             [gch tf] = GraphCutMex(gch, 't', varargin{2});\r
 263         elseif numel(varargin) == 1\r
 264             gch = varargin{1};\r
 265             [gch tf] = GraphCutMex(gch, 't');\r
 266         else\r
 267             error('GraphCut:Truncate: wrong number of input arguments');\r
 268         end\r
 269         if nargout > 2\r
 270             error('GraphCut:Truncate: too many output arguments');\r
 271         end\r
 272         if nargout == 2\r
 273             varargout{1} = tf;\r
 274         end\r
 275             \r
 276     otherwise\r
 277         error('GraphCut: Unrecognized mode %s', mode);\r
 278 end\r
 279 \r
 280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
 281 %\r
 282 %   Aux functions\r
 283 %\r
 284 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
 285 function gch = OpenGraph(varargin)\r
 286 % Usage [gch] = OpenGraph(Dc, Sc, [vC, hC]) - 2D grid\r
 287 % or    [gch] = OpenGraph(Dc, Sc, [Contrast]) -3D grid\r
 288 % or    [gch] = GraphCut(DataCost, SmoothnessCost, SparseSmoothness) - any graph\r
 289 nin = numel(varargin);\r
 290 if (nin~=2)  && (nin ~= 3) && (nin~=4) \r
 291     error('GraphCut:Open: wrong number of inputs');\r
 292 end\r
 293 \r
 294 % Data cost\r
 295 Dc = varargin{1};\r
 296 if ndims(Dc) == 4\r
 297     %%% 3D graph\r
 298     [R C Z L] = size(Dc);\r
 299     if ~strcmp(class(Dc),'single')\r
 300         Dc = single(Dc);\r
 301     end\r
 302     Dc = permute(Dc,[4 1 2 3]);\r
 303     Dc = Dc(:)';\r
 304     \r
 305     % smoothness cost\r
 306     Sc = varargin{2};\r
 307     if any( size(Sc) ~= [L L] )\r
 308         error('GraphCut:Open: smoothness cost has incorrect size');\r
 309     end\r
 310     if ~strcmp(class(Sc),'single')\r
 311         Sc = single(Sc);\r
 312     end\r
 313     Sc = Sc(:)';\r
 314     if nin == 3\r
 315         Contrast = varargin{3};\r
 316         if any( size(Contrast) ~= [R C Z] )\r
 317             error('GraphCut:Open: Contrast term is of wrong size');\r
 318         end\r
 319         if ~strcmp(class(Contrast),'single')\r
 320             Contrast = single(Contrast);\r
 321         end\r
 322         Contrast = Contrast(:);\r
 323         \r
 324         gch = GraphCut3dConstr(R, C, Z, L, Dc, Sc, Contrast);\r
 325     elseif nin == 2\r
 326         gch = GraphCut3dConstr(R, C, Z, L, Dc, Sc);\r
 327     else\r
 328         error('GraphCut:Open: wrong number of inputs for 3D graph');\r
 329     end\r
 330 elseif ndims(Dc) == 3\r
 331     %%% 2D graph\r
 332     [h w l] = size(Dc);\r
 333     if ~strcmp(class(Dc),'single')\r
 334         Dc = single(Dc);\r
 335     end\r
 336     Dc = permute(Dc,[3 2 1]);\r
 337     Dc = Dc(:)';\r
 338 \r
 339     % smoothness cost\r
 340     Sc = varargin{2};\r
 341     if any( size(Sc) ~= [l l] )\r
 342         error('GraphCut:Open: smoothness cost has incorrect size');\r
 343     end\r
 344     if ~strcmp(class(Sc),'single')\r
 345         Sc = single(Sc);\r
 346     end\r
 347     Sc = Sc(:)';\r
 348 \r
 349     if nin==4\r
 350         vC = varargin{3};\r
 351         if any( size(vC) ~= [h w] )\r
 352             error('GraphCut:Open: vertical cue size incorrect');\r
 353         end\r
 354         if ~strcmp(class(vC),'single')\r
 355             vC = single(vC);\r
 356         end\r
 357         vC = vC';\r
 358 \r
 359         hC = varargin{4};\r
 360         if any( size(hC) ~= [h w] )\r
 361             error('GraphCut:Open: horizontal cue size incorrect');\r
 362         end\r
 363         if ~strcmp(class(hC),'single')\r
 364             hC = single(hC);\r
 365         end\r
 366         hC = hC';\r
 367         gch = GraphCutConstr(w, h, l, Dc, Sc, vC(:), hC(:));\r
 368     elseif nin == 2\r
 369         gch = GraphCutConstr(w, h, l, Dc, Sc);\r
 370     else\r
 371         error('GraphCut:Open: wrong number of input for 2D grid');\r
 372     end\r
 373 elseif ndims(Dc) == 2\r
 374     %%% arbitrary graph\r
 375     if nin ~= 3\r
 376         error('GraphCut:Open', 'incorect number of inputs');\r
 377     end\r
 378     \r
 379     [nl np] = size(Dc);\r
 380     Sc = varargin{2};\r
 381     if any(size(Sc) ~= [nl nl])\r
 382         error('GraphCut:Open', 'Wrong size of Dc or Sc');\r
 383     end\r
 384     \r
 385     SparseSc = varargin{3};\r
 386     if any(size(SparseSc) ~= [np np])\r
 387         error('GraphCut:Open', 'Wrong size of SparseSc');\r
 388     end\r
 389         \r
 390     gch = GraphCutConstrSparse(single(Dc(:)), single(Sc(:)), SparseSc);\r
 391         \r
 392 else\r
 393     %%% Unknown dimensionality...\r
 394     error('GraphCut:Open: data cost has incorect dimensionality');\r
 395 end\r
 396 \r
 397 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
 398 function [gch] = SetLabels(varargin)\r
 399 % usgae [gch] = SetLabels(gch, labels)\r
 400 \r
 401 if nargin ~= 2\r
 402     error('GraphCut:SetLabels: wrong number of inputs');\r
 403 end\r
 404 gch = varargin{1};\r
 405 labels = varargin{2};\r
 406 \r
 407 if ~strcmp(class(labels), 'int32')\r
 408     labels = int32(labels);\r
 409 end\r
 410 labels = labels';\r
 411 [gch] = GraphCutMex(gch, 's', labels(:));\r
 412 \r
 413 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
 414 function [gch labels] = GetLabels(varargin)\r
 415 \r
 416 if nargin ~= 1\r
 417     error('GraphCut:GetLabels: wrong number of inputs');\r
 418 end\r
 419 gch = varargin{1};\r
 420 [gch labels] = GraphCutMex(gch, 'g');\r
 421 labels = labels';\r
 422 \r
 423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
 424 function [gch labels] = Expand(varargin)\r
 425 gch = varargin{1};\r
 426 switch nargin\r
 427     case 1\r
 428         [gch labels] = GraphCutMex(gch, 'e');\r
 429     case 2\r
 430         [gch labels] = GraphCutMex(gch, 'e', varargin{2});\r
 431     case 3\r
 432         [gch labels] = GraphCutMex(gch, 'e', varargin{2}, varargin{3});\r
 433     case 4\r
 434         ind = varargin{4};\r
 435         ind = int32(ind(:)-1)'; % convert to int32\r
 436         [gch labels] = GraphCutMex(gch, 'e', varargin{2}, varargin{3}, ind);\r
 437     otherwise\r
 438         error('GraphCut:Expand: wrong number of inputs');\r
 439 end\r
 440 labels = labels';\r
 441 \r
 442 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
 443 function [gch labels] = Swap(varargin)\r
 444 gch = varargin{1};\r
 445 switch nargin\r
 446     case 1\r
 447         [gch labels] = GraphCutMex(gch, 'a');\r
 448     case 2\r
 449         [gch labels] = GraphCutMex(gch, 'a', varargin{2});\r
 450     case 3\r
 451         [gch labels] = GraphCutMex(gch, 'a', varargin{2}, varargin{3});\r
 452     otherwise\r
 453         error('GraphCut:Swap: wrong number of inputarguments');\r
 454 end\r
 455 labels = labels';\r
 456 \r
 457 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r