THESE/codes/graphe/GCmex1.9/energy.h

   1 /* energy.h */\r
   2 /* Vladimir Kolmogorov (vnk@cs.cornell.edu), 2003. */\r
   3 \r
   4 /*\r
   5         This software implements an energy minimization technique described in\r
   6 \r
   7         What Energy Functions can be Minimized via Graph Cuts?\r
   8         Vladimir Kolmogorov and Ramin Zabih. \r
   9         To appear in IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI). \r
  10         Earlier version appeared in European Conference on Computer Vision (ECCV), May 2002. \r
  11 \r
  12         More specifically, it computes the global minimum of a function E of binary\r
  13         variables x_1, ..., x_n which can be written as a sum of terms involving\r
  14         at most three variables at a time:\r
  15 \r
  16                 E(x_1, ..., x_n) = \sum_{i}     E^{i}    (x_i)\r
  17                                  + \sum_{i,j}   E^{i,j}  (x_i, x_j)\r
  18                                  + \sum_{i,j,k} E^{i,j,k}(x_i, x_j, x_k)\r
  19 \r
  20         The method works only if each term is "regular". Definitions of regularity\r
  21         for terms E^{i}, E^{i,j}, E^{i,j,k} are given below as comments to functions\r
  22         add_term1(), add_term2(), add_term3(). \r
  23 \r
  24         This software can be used only for research purposes. IF YOU USE THIS SOFTWARE,\r
  25         YOU SHOULD CITE THE AFOREMENTIONED PAPER IN ANY RESULTING PUBLICATION.\r
  26 \r
  27         In order to use it, you will also need a MAXFLOW software which can be\r
  28         obtained from http://www.cs.cornell.edu/People/vnk/software.html\r
  29 \r
  30 \r
  31         Example usage\r
  32         (Minimizes the following function of 3 binary variables:\r
  33         E(x, y, z) = x - 2*y + 3*(1-z) - 4*x*y + 5*|y-z|):\r
  34 \r
  35         ///////////////////////////////////////////////////\r
  36 \r
  37         #include <stdio.h>\r
  38         #include "energy.h"\r
  39 \r
  40         void main()\r
  41         {\r
  42                 // Minimize the following function of 3 binary variables:\r
  43                 // E(x, y, z) = x - 2*y + 3*(1-z) - 4*x*y + 5*|y-z|\r
  44                    \r
  45                 Energy::Var varx, vary, varz;\r
  46                 Energy *e = new Energy();\r
  47 \r
  48                 varx = e -> add_variable();\r
  49                 vary = e -> add_variable();\r
  50                 varz = e -> add_variable();\r
  51 \r
  52                 e -> add_term1(varx, 0, 1);  // add term x \r
  53                 e -> add_term1(vary, 0, -2); // add term -2*y\r
  54                 e -> add_term1(varz, 3, 0);  // add term 3*(1-z)\r
  55 \r
  56                 e -> add_term2(x, y, 0, 0, 0, -4); // add term -4*x*y\r
  57                 e -> add_term2(y, z, 0, 5, 5, 0); // add term 5*|y-z|\r
  58 \r
  59                 Energy::TotalValue Emin = e -> minimize();\r
  60                 \r
  61                 printf("Minimum = %d\n", Emin);\r
  62                 printf("Optimal solution:\n");\r
  63                 printf("x = %d\n", e->get_var(varx));\r
  64                 printf("y = %d\n", e->get_var(vary));\r
  65                 printf("z = %d\n", e->get_var(varz));\r
  66 \r
  67                 delete e;\r
  68         }\r
  69 \r
  70         ///////////////////////////////////////////////////\r
  71 */\r
  72 \r
  73 #ifndef __ENERGY_H__\r
  74 #define __ENERGY_H__\r
  75 \r
  76 #include "bagon_assert.h" // <assert.h>\r
  77 #include "graph.h"\r
  78 \r
  79 \r
  80 \r
  81 class Energy : Graph\r
  82 {\r
  83 public:\r
  84         typedef node_id Var;\r
  85     \r
  86     \r
  87         /* Types of energy values.\r
  88            Value is a type of a value in a single term\r
  89            TotalValue is a type of a value of the total energy.\r
  90            By default Value = short, TotalValue = int.\r
  91            To change it, change the corresponding types in graph.h */\r
  92         typedef captype Value;\r
  93         typedef flowtype TotalValue;\r
  94 \r
  95         /* interface functions */\r
  96 \r
  97         /* Constructor. Optional argument is the pointer to the\r
  98            function which will be called if an error occurs;\r
  99            an error message is passed to this function. If this\r
 100            argument is omitted, exit(1) will be called. */\r
 101         Energy(void (*err_function)(const char *) = NULL, bool tf = false);\r
 102     \r
 103         \r
 104         /* Destructor */\r
 105         ~Energy();\r
 106 \r
 107         /* Adds a new binary variable */\r
 108         Var add_variable();\r
 109 \r
 110         /* Adds a constant E to the energy function */\r
 111         void add_constant(Value E);\r
 112 \r
 113         /* Adds a new term E(x) of one binary variable\r
 114            to the energy function, where\r
 115                E(0) = E0, E(1) = E1\r
 116            E0 and E1 can be arbitrary */\r
 117         void add_term1(Var x,\r
 118                        Value E0, Value E1);\r
 119 \r
 120         /* Adds a new term E(x,y) of two binary variables\r
 121            to the energy function, where\r
 122                E(0,0) = E00, E(0,1) = E01\r
 123                E(1,0) = E10, E(1,1) = E11\r
 124            The term must be regular, i.e. E00 + E11 <= E01 + E10 */\r
 125         void add_term2(Var x, Var y,\r
 126                        Value E00, Value E01,\r
 127                        Value E10, Value E11);\r
 128 \r
 129         /* Adds a new term E(x,y,z) of three binary variables\r
 130            to the energy function, where\r
 131                E(0,0,0) = E000, E(0,0,1) = E001\r
 132                E(0,1,0) = E010, E(0,1,1) = E011\r
 133                E(1,0,0) = E100, E(1,0,1) = E101\r
 134                E(1,1,0) = E110, E(1,1,1) = E111\r
 135            The term must be regular. It means that if one\r
 136            of the variables is fixed (for example, y=1), then\r
 137            the resulting function of two variables must be regular.\r
 138            Since there are 6 ways to fix one variable\r
 139            (3 variables times 2 binary values - 0 and 1),\r
 140            this is equivalent to 6 inequalities */\r
 141         void add_term3(Var x, Var y, Var z,\r
 142                        Value E000, Value E001,\r
 143                        Value E010, Value E011,\r
 144                        Value E100, Value E101,\r
 145                        Value E110, Value E111);\r
 146 \r
 147         /* After the energy function has been constructed,\r
 148            call this function to minimize it.\r
 149            Returns the minimum of the function */\r
 150         TotalValue minimize();\r
 151 \r
 152         /* After 'minimize' has been called, this function\r
 153            can be used to determine the value of variable 'x'\r
 154            in the optimal solution.\r
 155            Returns either 0 or 1 */\r
 156         int get_var(Var x);\r
 157 \r
 158 /***********************************************************************/\r
 159 /***********************************************************************/\r
 160 /***********************************************************************/\r
 161 \r
 162 private:\r
 163         /* internal variables and functions */\r
 164 \r
 165         TotalValue      Econst;\r
 166         void            (*error_function)(const char *);        /* this function is called if a error occurs,\r
 167                                                                                         with a corresponding error message\r
 168                                                                                         (or exit(1) is called if it's NULL) */\r
 169     // BAGON\r
 170     bool truncate_flag;\r
 171     // \r
 172 };\r
 173 \r
 174 \r
 175 \r
 176 \r
 177 \r
 178 \r
 179 \r
 180 \r
 181 \r
 182 \r
 183 \r
 184 \r
 185 \r
 186 \r
 187 \r
 188 /***********************************************************************/\r
 189 /************************  Implementation ******************************/\r
 190 /***********************************************************************/\r
 191 \r
 192 \r
 193 inline Energy::Energy(void (*err_function)(const char *), bool tf) : Graph(err_function), truncate_flag(tf)\r
 194 {\r
 195         Econst = 0;\r
 196         error_function = err_function;\r
 197     \r
 198 }\r
 199 \r
 200 \r
 201 inline Energy::~Energy() {}\r
 202 \r
 203 \r
 204 inline Energy::Var Energy::add_variable() {     return add_node(); }\r
 205 \r
 206 \r
 207 inline void Energy::add_constant(Value A) { Econst += A; }\r
 208 \r
 209 \r
 210 inline void Energy::add_term1(Var x,\r
 211                               Value A, Value B)\r
 212 {\r
 213         add_tweights(x, B, A);\r
 214 }\r
 215 \r
 216 \r
 217 inline void Energy::add_term2(Var x, Var y,\r
 218                               Value A, Value B,\r
 219                               Value C, Value D)\r
 220 {\r
 221     \r
 222     \r
 223     \r
 224     /* truncation */\r
 225     if ( truncate_flag && ( A + D > C + B ) ) {\r
 226         \r
 227         Value delta = .5*(D-B-C);\r
 228         B += delta;\r
 229         C += delta;\r
 230     }\r
 231         /* \r
 232            E = A A  +  0   B-A\r
 233                D D     C-D 0\r
 234            Add edges for the first term\r
 235         */\r
 236         add_tweights(x, D, A);\r
 237         B -= A; C -= D;\r
 238 \r
 239         /* now need to represent\r
 240            0 B\r
 241            C 0\r
 242         */\r
 243 \r
 244         assert(B + C >= 0); /* check regularity (i.e., triangle inequality) */\r
 245         if (B < 0)\r
 246         {\r
 247                 /* Write it as\r
 248                    B B  +  -B 0  +  0   0\r
 249                    0 0     -B 0     B+C 0\r
 250                 */\r
 251                 add_tweights(x, 0, B); /* first term */\r
 252                 add_tweights(y, 0, -B); /* second term */\r
 253                 add_edge(x, y, 0, B+C); /* third term */\r
 254         }\r
 255         else if (C < 0)\r
 256         {\r
 257                 /* Write it as\r
 258                    -C -C  +  C 0  +  0 B+C\r
 259                     0  0     C 0     0 0\r
 260                 */\r
 261                 add_tweights(x, 0, -C); /* first term */\r
 262                 add_tweights(y, 0, C); /* second term */\r
 263                 add_edge(x, y, B+C, 0); /* third term */\r
 264         }\r
 265         else /* B >= 0, C >= 0 */\r
 266         {\r
 267                 add_edge(x, y, B, C);\r
 268         }\r
 269 }\r
 270 \r
 271 inline void Energy::add_term3(Var x, Var y, Var z,\r
 272                               Value E000, Value E001,\r
 273                               Value E010, Value E011,\r
 274                               Value E100, Value E101,\r
 275                               Value E110, Value E111)\r
 276 {\r
 277         register Value pi = (E000 + E011 + E101 + E110) - (E100 + E010 + E001 + E111);\r
 278         register Value delta;\r
 279         register Var u;\r
 280 \r
 281         if (pi >= 0)\r
 282         {\r
 283                 Econst += E111 - (E011 + E101 + E110);\r
 284 \r
 285                 add_tweights(x, E101, E001);\r
 286                 add_tweights(y, E110, E100);\r
 287                 add_tweights(z, E011, E010);\r
 288 \r
 289                 delta = (E010 + E001) - (E000 + E011); /* -pi(E[x=0]) */\r
 290                 assert(delta >= 0); /* check regularity */\r
 291                 add_edge(y, z, delta, 0);\r
 292 \r
 293                 delta = (E100 + E001) - (E000 + E101); /* -pi(E[y=0]) */\r
 294                 assert(delta >= 0); /* check regularity */\r
 295                 add_edge(z, x, delta, 0);\r
 296 \r
 297                 delta = (E100 + E010) - (E000 + E110); /* -pi(E[z=0]) */\r
 298                 assert(delta >= 0); /* check regularity */\r
 299                 add_edge(x, y, delta, 0);\r
 300 \r
 301                 if (pi > 0)\r
 302                 {\r
 303                         u = add_variable();\r
 304                         add_edge(x, u, pi, 0);\r
 305                         add_edge(y, u, pi, 0);\r
 306                         add_edge(z, u, pi, 0);\r
 307                         add_tweights(u, 0, pi);\r
 308                 }\r
 309         }\r
 310         else\r
 311         {\r
 312                 Econst += E000 - (E100 + E010 + E001);\r
 313 \r
 314                 add_tweights(x, E110, E010);\r
 315                 add_tweights(y, E011, E001);\r
 316                 add_tweights(z, E101, E100);\r
 317 \r
 318                 delta = (E110 + E101) - (E100 + E111); /* -pi(E[x=1]) */\r
 319                 assert(delta >= 0); /* check regularity */\r
 320                 add_edge(z, y, delta, 0);\r
 321 \r
 322                 delta = (E110 + E011) - (E010 + E111); /* -pi(E[y=1]) */\r
 323                 assert(delta >= 0); /* check regularity */\r
 324                 add_edge(x, z, delta, 0);\r
 325 \r
 326                 delta = (E101 + E011) - (E001 + E111); /* -pi(E[z=1]) */\r
 327                 assert(delta >= 0); /* check regularity */\r
 328                 add_edge(y, x, delta, 0);\r
 329 \r
 330                 u = add_variable();\r
 331                 add_edge(u, x, -pi, 0);\r
 332                 add_edge(u, y, -pi, 0);\r
 333                 add_edge(u, z, -pi, 0);\r
 334                 add_tweights(u, -pi, 0);\r
 335         }\r
 336 }\r
 337 \r
 338 inline Energy::TotalValue Energy::minimize() { return Econst + maxflow(); }\r
 339 \r
 340 inline int Energy::get_var(Var x) { return (int) what_segment(x); }\r
 341 \r
 342 #endif\r