************************************************** * MRF energy minimization software * * Version 1.6 * * May 5, 2006 * ************************************************** This directory contains the MRF energy minimization software accompanying the paper [1] A Comparative Study of Energy Minimization Methods for Markov Random Fields. R. Szeliski, R. Zabih, D. Scharstein, O. Veksler, V. Kolmogorov, A. Agarwala, M. Tappen, and C. Rother. In Ninth European Conference on Computer Vision (ECCV 2006), volume 2, pages 16-29, Graz, Austria, May 2006. Optimization methods included: 1) ICM 2) Graph Cuts 3) Max-Product Belief Propagation The TRW / TRW-S method is currently not included, but hopefully will be in a future version. Instructions for compiling and using the software are included below. CREDITS: * MRF code, graph cut interface, and example code by Olga Veksler * Graph cut library by Yuri Boykov and Vladimir Kolmogorov * Belief propagation code by Marshall Tappen ================================ INSTRUCTIONS FOR CITATIONS ======================== If you use this software, you should cite the paper [1]. In addition, since this software builds on the algorithms and libraries developed by many different people, you should also follow the instructions about proper citing below. (a) If you use the GraphCuts optimization (written by Olga Veksler, using the libraries provided by Yuri Boykov and Vladimir Kolmogorov), you should cite the following papers: [2] Fast Approximate Energy Minimization via Graph Cuts. Y. Boykov, O. Veksler, and R. Zabih. In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), vol. 23, no. 11, pages 1222-1239, November 2001. [3] What Energy Functions can be Minimized via Graph Cuts? V. Kolmogorov and R. Zabih. In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), vol. 26, no. 2, pages 147-159, February 2004. An earlier version appeared in European Conference on Computer Vision (ECCV), May 2002. [4] An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision. Y. Boykov and Vladimir Kolmogorov. In IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), vol. 26, no. 9, pages 1124-1137, September 2004. The graph-cuts software can be used only for research purposes. If you wish to use the graph-cuts software for commercial purposes, you should be aware that there is a US patent: R. Zabih, Y. Boykov, O. Veksler "System and method for fast approximate energy minimization via graph cuts", United Stated Patent 6,744,923, June 1, 2004 Citation [2] develops the graph-cuts algorithms (alpha-expansion and the swap algorithms), citation [3] gives a simpler graph construction for the methods described in [2], and citation [3] gives an efficient min-cut/max-flow algorithm for computing the minimum graph cut. The file energy.h included here is a slight modification of the file energy.h in energy-v1.1.src by Vladimir Kolmogorov, available at http://www.adastral.ucl.ac.uk/~vladkolm/software.html The files block.h, graph.cpp, graph.h, and maxflow.cpp included here are slight modifications of the files in maxflow-v2.2.src/forward_star, by Vladimir Kolmogorov and Yuri Boykov, available at http://www.adastral.ucl.ac.uk/~vladkolm/software.html (b) If you are using the Belief Propagation software (provided by Marshall Tappen), you should cite [5] M. F. Tappen and W. T. Freeman. Comparison of Graph Cuts with Belief Propagation for Stereo, using Identical MRF Parameters. In Proceedings of the Ninth IEEE International Conference on Computer Vision (ICCV), pages 900-907, 2003 ============================= INSTRUCTIONS FOR USAGE ============================================= A simple Makefile is included - typing "make" will compile the library and a sample driver file "example.cpp" into an executable "example". Here is a brief description of how to use the code; see also example.cpp and mrf.h. The default type for energy is float. If you want to change this type to other numerical type, such as int, double, etc., in mrf.h, change the following 2 lines to the desired type: typedef int EnergyVal; /* The total energy of a labeling */ typedef int CostVal; /* costs of individual terms of the energy */ Note that CostVal is the type of individual energy terms (such as the data cost for pixel p and label l), and EnergyVal is the type for the total energy value. Therefore, the sum of all individual energy terms should not "overflow" the EnergyVal type. Step 1: Set up an energy function An energy function is specified by setting data costs and smoothness costs. Data costs and smoothness costs can be specified by an array or function. In addition, smoothness cost can be specified by parameters lambda, k, V_max, as explained in section 2 of paper [1], in which case the smoothness cost V(|l1-l2|) = lambda*min(V_max,|l1-l2|^k). Spatially varying per-pairing weights w_pq can be also specified. (a) Setting up Data Terms There are 2 ways to set up data terms: (i) you can allocate an array dCost of size (numberofPixels)*(numberofLabels) of the CostVal type, where dCost[pix*numberofLabels + l] is the cost of assigning label l to pixel pix (ii) or you can set up a function MRF::CostVal dCost(int pix, int i), which takes pixel pix and label i and returns the data cost After you initialized data cost array, or setup a function data cost, to initialize data cost for optimization, call: DataCost *data = new DataCost(dCost); (b) Setting up Smoothness Terms There are 3 ways to set up smoothness terms: (i) Set up a symmetric array smooth of size numberofLabels*numberofLabels, where V(l1,l2) = smooth[l1+numberofLabels*l2] = smooth[l2+numberofLabels*l1] (ii) Specify V_max and k (iii) The most general way is to set up afunction MRF::CostVal fnCost(int pix1, int pix2, int i, int j) which takes pixels pix1 and pix2, labels i and j and returns the smoothness penalty V for assigning labels i to pix1 and j to pix2. In addition, if you used (i) or (ii) to specify the smoothness penalty V, and if the neighborhood system is the standard 4 connected grid, you can specify the spacial varying weights w_pq in arrays vCue, hCue, each of which has size width*height, where width is the width of the grid and height is the height of the set of the grid. hCue[x+y*width] holds the variable weight for edge between pixels (x+1,y) and (x,y) and vCue[x+y*width] holds the variable weight for edge between pixels (x,y+1) and (x,y) After you've done steps (i), or (ii), or (iii), to initialize smoothness cost for optimization, call: In case of (i) and NO spacially varying weights w_pq: SmoothnessCost *smooth = new SmoothnessCost(smooth); In case of (i) and spacially varying weights w_pq: SmoothnessCost *smooth = new SmoothnessCost(smooth,hCue,vCue); Note that in this case, V(pix1,pix2,l1,l2) = smooth[l1+l2*numofLabels]*hCue[min(pix1,pix2)] if pix1 and pix2 are horizontal neighbors, and V(pix1,pix2,l1,l2) = smooth[l1+l2*numofLabels]*vCue[min(pix1,pix2)] if pix1 and pix2 are vertical neighbors In case of (ii) and NO spacially varying weights w_pq: SmoothnessCost *smooth = new SmoothnessCost(k,V_max,lambda); In case of (ii) and spacially varying weights w_pq: SmoothnessCost *smooth = new SmoothnessCost(k,V_max,lambda,hCue,vCue); Note that in this case, V(pix1,pix2,l1,l2) = lambda*min(V_max,|l1-l2]^k)*hCue[min(pix1,pix2)] if pix1 and pix2 are horizontal neighbors, and V(pix1,pix2,l1,l2) = lambda*min(V_max,|l1-l2]^k)*vCue[min(pix1,pix2)] if pix1 and pix2 are vertical neighbors In case of (iii): SmoothnessCost *smooth = new SmoothnessCost(fnCost); After the data and smoothness terms are set-up, the energy can be specified by calling, EnergyFunction *eng = new EnergyFunction(data,smooth); Note that we directly use all the arrays that you make, that is we do not copy them. So don't free any memory in the arrays that you allocate until you are done with optimization. Step 2: Invoke an optimization algorithm, assuming the grid has dimensions width by height, and number of labels is numberOfLabels, and eng is set up as above: For ICM: float t; MRF* mrf = new ICM(width,height,numberOfLabels,eng); mrf->initialize(); mrf->clearAnswer(); mrf->optimize(5,t); // run for 5 iterations, store time t it took MRF::EnergyVal E_smooth = mrf->smoothnessEnergy(); MRF::EnergyVal E_data = mrf->dataEnergy(); printf("Total Energy = %d (Smoothness energy %d, Data Energy %d)\n", E_smooth+E_data,E_smooth,E_data); for (int pix =0; pix < width*height; pix++ ) printf("Label of pixel %d is %d",pix, mrf->getLabel(pix)); delete mrf; For Graph-cuts with expansion algorithm: float t; MRF* mrf = new Expansion(width,height,numberOfLabels,eng); mrf->initialize(); mrf->clearAnswer(); mrf->optimize(5,t); // run for 5 iterations, store time t it took MRF::EnergyVal E_smooth = mrf->smoothnessEnergy(); MRF::EnergyVal E_data = mrf->dataEnergy(); printf("Total Energy = %d (Smoothness energy %d, Data Energy %d)\n", E_smooth+E_data,E_smooth,E_data); for (int pix =0; pix < width*height; pix++ ) printf("Label of pixel %d is %d",pix, mrf->getLabel(pix)); delete mrf; For Graph-cuts with swap algorithm: float t; MRF* mrf = new Swap(width,height,numberOfLabels,eng); mrf->initialize(); mrf->clearAnswer(); mrf->optimize(5,t); // run for 5 iterations, store time t it took MRF::EnergyVal E_smooth = mrf->smoothnessEnergy(); MRF::EnergyVal E_data = mrf->dataEnergy(); printf("Total Energy = %d (Smoothness energy %d, Data Energy %d)\n", E_smooth+E_data,E_smooth,E_data); for (int pix =0; pix < width*height; pix++ ) printf("Label of pixel %d is %d",pix, mrf->getLabel(pix)); delete mrf; For Max-product Belief Propagation: float t; MRF* mrf = new MaxProdBP(width,height,numberOfLabels,eng); mrf->initialize(); mrf->clearAnswer(); mrf->optimize(5,t); // run for 5 iterations, store time t it took MRF::EnergyVal E_smooth = mrf->smoothnessEnergy(); MRF::EnergyVal E_data = mrf->dataEnergy(); printf("Total Energy = %d (Smoothness energy %d, Data Energy %d)\n", E_smooth+E_data,E_smooth,E_data); for (int pix =0; pix < width*height; pix++ ) printf("Label of pixel %d is %d",pix, mrf->getLabel(pix)); delete mrf; The above code works for grid graphs. Currently, the graph cuts algorithms (swap and expansion) are also implemented for general neighborhood systems, in which case after specifying energy you also have to specify the neighborhood system using function: mrf->setNeighbors(int pix1, int pix2, CostVal weight); which makes pix1 and pix2 neighbors with spacially varying weight. Note that in this case, you can't use hCue and vCue, when setting up the smoothness terms, since they only work for 4-connected grid. Also, if you use general neighborhood system, you have to use a different constructor: For Graph-cuts with expansion algorithm: float t; MRF* mrf = new Expansion(numberofPixels,numberOfLabels,eng); for (int i = 0; i setNeighbors(i,j, (i-j)*(i-j)); mrf->initialize(); mrf->clearAnswer(); mrf->optimize(5,t); // run for 5 iterations, store time t it took MRF::EnergyVal E_smooth = mrf->smoothnessEnergy(); MRF::EnergyVal E_data = mrf->dataEnergy(); printf("Total Energy = %d (Smoothness energy %d, Data Energy %d)\n", E_smooth+E_data,E_smooth,E_data); for (int pix =0; pix < width*height; pix++ ) printf("Label of pixel %d is %d",pix, mrf->getLabel(pix)); delete mrf; For Graph-cuts with swap algorithm: float t; MRF* mrf = new Swap(numberofPixels,numberOfLabels,eng); for (int i = 0; i setNeighbors(i,j, (i-j)*(i-j)); mrf->initialize(); mrf->clearAnswer(); mrf->optimize(5,t); // run for 5 iterations, store time t it took MRF::EnergyVal E_smooth = mrf->smoothnessEnergy(); MRF::EnergyVal E_data = mrf->dataEnergy(); printf("Total Energy = %d (Smoothness energy %d, Data Energy %d)\n", E_smooth+E_data,E_smooth,E_data); for (int pix =0; pix < width*height; pix++ ) printf("Label of pixel %d is %d",pix, mrf->getLabel(pix)); delete mrf;