diff -rupN MRF2.1/BP-S.cpp MRF-v2.1.src-patched/BP-S.cpp
--- MRF2.1/BP-S.cpp	2006-11-15 04:27:02.000000000 +0100
+++ MRF-v2.1.src-patched/BP-S.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -9,6 +9,7 @@
 #include "typeTruncatedQuadratic2D.h"
 #undef private
 
+using namespace mrfLib;
 
 #define m_D(pix,l)  m_D[(pix)*m_nLabels+(l)]
 #define m_V(l1,l2)  m_V[(l1)*m_nLabels+(l2)]
diff -rupN MRF2.1/BP-S.h MRF-v2.1.src-patched/BP-S.h
--- MRF2.1/BP-S.h	2006-11-15 04:27:02.000000000 +0100
+++ MRF-v2.1.src-patched/BP-S.h	2012-08-13 12:25:46.000000000 +0200
@@ -1,5 +1,5 @@
-#ifndef __BPS_H__
-#define __BPS_H__
+#ifndef __MRF_BPS_H__
+#define __MRF_BPS_H__
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -7,6 +7,7 @@
 #include <assert.h>
 #include "mrf.h"
 
+namespace mrfLib {
 
 class BPS : public MRF{
  public:
@@ -44,7 +45,7 @@ class BPS : public MRF{
 	    L2,
 	    FIXED_MATRIX,
 	    GENERAL,
-	    BINARY,
+	    BINARY
 	} m_type;
 
     CostVal m_smoothMax; // used only if
@@ -76,4 +77,6 @@ class BPS : public MRF{
     void optimize_GRID_BINARY(int nIterations);
 };
 
+} // end namespace mrfLib
+
 #endif /*  __BPS_H__ */
diff -rupN MRF2.1/GCoptimization.cpp MRF-v2.1.src-patched/GCoptimization.cpp
--- MRF2.1/GCoptimization.cpp	2007-03-18 23:39:34.000000000 +0100
+++ MRF-v2.1.src-patched/GCoptimization.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -7,7 +7,7 @@
 #include "string.h"
 #define MAX_INTT 1000000000
 
-
+using namespace mrfLib;
 
 /**************************************************************************************/
 
diff -rupN MRF2.1/GCoptimization.h MRF-v2.1.src-patched/GCoptimization.h
--- MRF2.1/GCoptimization.h	2007-03-18 23:53:50.000000000 +0100
+++ MRF-v2.1.src-patched/GCoptimization.h	2012-08-13 12:25:46.000000000 +0200
@@ -25,13 +25,16 @@
 
 
 
-#ifndef __GCOPTIMIZATION_H__
-#define __GCOPTIMIZATION_H__
+#ifndef __MRF_GCOPTIMIZATION_H__
+#define __MRF_GCOPTIMIZATION_H__
 
 #include "LinkedBlockList.h"
 #include <assert.h>
 #include "graph.h"
 #include "energy.h"
+
+namespace mrfLib {
+
 #define m_datacost(pix,lab)     (m_datacost[(pix)*m_nLabels+(lab)] )
 #define m_smoothcost(lab1,lab2) (m_smoothcost[(lab1)+(lab2)*m_nLabels] )
 #define USE_MEMBER_FUNCTION 0
@@ -288,7 +291,7 @@ private:
 };
 
 
-
+} // end namespace mrfLib
 
 #endif
 
diff -rupN MRF2.1/ICM.cpp MRF-v2.1.src-patched/ICM.cpp
--- MRF2.1/ICM.cpp	2006-11-15 04:27:04.000000000 +0100
+++ MRF-v2.1.src-patched/ICM.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -9,7 +9,7 @@
 
 
         
-
+using namespace mrfLib;
 
 ICM::ICM(int width, int height, int nLabels,EnergyFunction *eng):MRF(width,height,nLabels,eng)
 {
@@ -33,7 +33,7 @@ ICM::~ICM()
 void ICM::initializeAlg()
 {
     m_answer = (Label *) new Label[m_nPixels];
-    if ( !m_answer ){printf("\nNot enough memory, exiting");exit(0);}
+    if ( !m_answer ) {printf("\nNot enough memory, exiting");exit(0);}
 
     if (!m_grid_graph)
     {
diff -rupN MRF2.1/ICM.h MRF-v2.1.src-patched/ICM.h
--- MRF2.1/ICM.h	2006-11-15 04:27:04.000000000 +0100
+++ MRF-v2.1.src-patched/ICM.h	2012-08-13 12:25:46.000000000 +0200
@@ -1,5 +1,5 @@
-#ifndef __ICM_H__
-#define __ICM_H__
+#ifndef __MRF_ICM_H__
+#define __MRF_ICM_H__
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -8,6 +8,7 @@
 #include "mrf.h"
 #include "LinkedBlockList.h"
 
+namespace mrfLib {
 
 class ICM : public MRF{
 public:
@@ -51,6 +52,7 @@ private:
     LinkedBlockList *m_neighbors;
 };
 
+} // end namespace mrfLib
 
 #endif /*  __ICM_H__ */
 
diff -rupN MRF2.1/LinkedBlockList.cpp MRF-v2.1.src-patched/LinkedBlockList.cpp
--- MRF2.1/LinkedBlockList.cpp	2006-11-15 04:27:04.000000000 +0100
+++ MRF-v2.1.src-patched/LinkedBlockList.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -1,5 +1,7 @@
 #include "LinkedBlockList.h"
 
+using namespace mrfLib;
+
 /*********************************************************************/
 
 void LinkedBlockList::addFront(ListType item) {
diff -rupN MRF2.1/LinkedBlockList.h MRF-v2.1.src-patched/LinkedBlockList.h
--- MRF2.1/LinkedBlockList.h	2006-11-15 04:27:04.000000000 +0100
+++ MRF-v2.1.src-patched/LinkedBlockList.h	2012-08-13 12:25:46.000000000 +0200
@@ -6,8 +6,10 @@
 // For other graphs, it should be set to the average expected number of neighbors
 // Data in linked list for the neighborhood system is allocated in blocks of size GCLL_BLOCK_SIZE 
 
-#ifndef __LINKEDBLOCKLIST_H__
-#define __LINKEDBLOCKLIST_H__
+#ifndef __MRF_LINKEDBLOCKLIST_H__
+#define __MRF_LINKEDBLOCKLIST_H__
+
+namespace mrfLib {
 
 #define GCLL_BLOCK_SIZE 4  
 // GCLL_BLOCKSIZE should "fit" into the type BlockType. That is 
@@ -46,4 +48,6 @@ private:
     LLBlock *m_cursor;
 };
 
+} // end namespace mrfLib
+
 #endif
diff -rupN MRF2.1/Makefile MRF-v2.1.src-patched/Makefile
--- MRF2.1/Makefile	2007-10-22 21:24:18.000000000 +0200
+++ MRF-v2.1.src-patched/Makefile	1970-01-01 01:00:00.000000000 +0100
@@ -1,47 +0,0 @@
-VERSION = MRF2.1
-
-SRC =  mrf.cpp ICM.cpp GCoptimization.cpp graph.cpp maxflow.cpp \
-       MaxProdBP.cpp LinkedBlockList.cpp regions-maxprod.cpp \
-       TRW-S.cpp BP-S.cpp
-
-CC = g++
-
-WARN   = -W -Wall
-OPT ?= -O3
-CPPFLAGS = $(OPT) $(WARN)
-#CPPFLAGS = $(OPT) $(WARN) -DUSE_64_BIT_PTR_CAST
-
-OBJ = $(SRC:.cpp=.o)
-
-all: libMRF.a example
-
-libMRF.a: $(OBJ)
-	rm -f libMRF.a
-	ar ruc libMRF.a $(OBJ)
-	ranlib libMRF.a
-
-example: libMRF.a example.cpp
-	$(CC) -o example example.cpp -L. -lMRF
-
-clean: 
-	rm -f $(OBJ) core *.stackdump *.bak
-
-allclean: clean
-	rm -f libMRF.a example example.exe
-
-depend:
-	@makedepend -Y -- $(CPPFLAGS) -- $(SRC) 2>> /dev/null
-
-# DO NOT DELETE THIS LINE -- make depend depends on it.
-
-mrf.o: mrf.h
-ICM.o: ICM.h mrf.h LinkedBlockList.h
-GCoptimization.o: energy.h graph.h block.h mrf.h GCoptimization.h
-GCoptimization.o: LinkedBlockList.h
-graph.o: graph.h block.h mrf.h
-maxflow.o: graph.h block.h mrf.h
-MaxProdBP.o: MaxProdBP.h mrf.h LinkedBlockList.h regions-new.h
-LinkedBlockList.o: LinkedBlockList.h
-regions-maxprod.o: MaxProdBP.h mrf.h LinkedBlockList.h regions-new.h
-TRW-S.o: TRW-S.h mrf.h typeTruncatedQuadratic2D.h
-BP-S.o: BP-S.h mrf.h typeTruncatedQuadratic2D.h
diff -rupN MRF2.1/MaxProdBP.cpp MRF-v2.1.src-patched/MaxProdBP.cpp
--- MRF2.1/MaxProdBP.cpp	2006-11-27 16:14:00.000000000 +0100
+++ MRF-v2.1.src-patched/MaxProdBP.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -6,6 +6,8 @@
 #include "MaxProdBP.h"
 #include "regions-new.h"
 
+using namespace mrfLib;
+
 #define m_D(pix,l)  m_D[(pix)*m_nLabels+(l)]
 #define m_V(l1,l2)  m_V[(l1)*m_nLabels+(l2)]
 
@@ -39,7 +41,7 @@ void MaxProdBP::initializeAlg()
 void MaxProdBP::BPinitializeAlg()
 {
 	m_answer = (Label *) new Label[m_nPixels];
-	if ( !m_answer ){printf("\nNot enough memory, exiting");exit(0);}
+	if ( !m_answer ) {printf("\nNot enough memory, exiting");exit(0);}
 
 	m_scratchMatrix = new FLOATTYPE[m_nLabels * m_nLabels];
 	// MEMORY LEAK? where does this ever get deleted??
@@ -62,7 +64,7 @@ void MaxProdBP::BPinitializeAlg()
 	  const int clen = 4*m_nPixels * m_nLabels ;
 	  //printf("clen:%d\n",clen/1024/1024);
 	  m_message_chunk = (FloatType *) new FloatType[clen];
-	  if ( !m_message_chunk ){printf("\nNot enough memory for messages, exiting");exit(0);}
+	  if ( !m_message_chunk ) {printf("\nNot enough memory for messages, exiting");exit(0);}
 	  for(int i = 0; i < clen; i++)
 	    m_message_chunk[i]=0;
 	  initOneNodeMsgMem(nodeArray, m_message_chunk, m_nPixels, m_nLabels);
diff -rupN MRF2.1/MaxProdBP.h MRF-v2.1.src-patched/MaxProdBP.h
--- MRF2.1/MaxProdBP.h	2006-11-15 04:27:06.000000000 +0100
+++ MRF-v2.1.src-patched/MaxProdBP.h	2012-08-13 12:25:46.000000000 +0200
@@ -1,5 +1,5 @@
-#ifndef __MAXPRODBP_H__
-#define __MAXPRODBP_H__
+#ifndef __MRF_MAXPRODBP_H__
+#define __MRF_MAXPRODBP_H__
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -9,6 +9,8 @@
 #include "LinkedBlockList.h"
 #include "regions-new.h"
 
+namespace mrfLib {
+
 #define FloatType float
 #define FLOATTYPE float
 class MaxProdBP;
@@ -55,7 +57,7 @@ public:
       L2,
       FIXED_MATRIX,
       GENERAL,
-      BINARY,
+      BINARY
     } m_type;
 
 protected:
@@ -91,6 +93,7 @@ private:
   LinkedBlockList *m_neighbors;
 };
 
+} // end namespace mrfLib
 
 #endif /*  __ICM_H__ */
 
diff -rupN MRF2.1/TRW-S.cpp MRF-v2.1.src-patched/TRW-S.cpp
--- MRF2.1/TRW-S.cpp	2006-11-15 04:27:06.000000000 +0100
+++ MRF-v2.1.src-patched/TRW-S.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -9,6 +9,7 @@
 #include "typeTruncatedQuadratic2D.h"
 #undef private
 
+using namespace mrfLib;
 
 #define m_D(pix,l)  m_D[(pix)*m_nLabels+(l)]
 #define m_V(l1,l2)  m_V[(l1)*m_nLabels+(l2)]
diff -rupN MRF2.1/TRW-S.h MRF-v2.1.src-patched/TRW-S.h
--- MRF2.1/TRW-S.h	2006-11-15 04:27:06.000000000 +0100
+++ MRF-v2.1.src-patched/TRW-S.h	2012-08-13 12:25:46.000000000 +0200
@@ -1,5 +1,5 @@
-#ifndef __TRWS_H__
-#define __TRWS_H__
+#ifndef __MRF_TRWS_H__
+#define __MRF_TRWS_H__
 
 #include <stdio.h>
 #include <stdlib.h>
@@ -7,6 +7,7 @@
 #include <assert.h>
 #include "mrf.h"
 
+namespace mrfLib {
 
 class TRWS : public MRF{
  public:
@@ -49,7 +50,7 @@ class TRWS : public MRF{
 	    L2,
 	    FIXED_MATRIX,
 	    GENERAL,
-	    BINARY,
+	    BINARY
 	} m_type;
 
     CostVal m_smoothMax; // used only if
@@ -83,4 +84,6 @@ class TRWS : public MRF{
     void optimize_GRID_BINARY(int nIterations);
 };
 
+} // end namespace mrfLib
+
 #endif /*  __TRWS_H__ */
diff -rupN MRF2.1/block.h MRF-v2.1.src-patched/block.h
--- MRF2.1/block.h	2006-11-15 04:27:00.000000000 +0100
+++ MRF-v2.1.src-patched/block.h	2012-08-13 12:25:46.000000000 +0200
@@ -106,8 +106,8 @@
     deallocated only when the destructor is called.
 */
 
-#ifndef __BLOCK_H__
-#define __BLOCK_H__
+#ifndef __MRF_BLOCK_H__
+#define __MRF_BLOCK_H__
 
 #include <stdlib.h>
 
@@ -115,6 +115,8 @@
 /***********************************************************************/
 /***********************************************************************/
 
+namespace mrfLib {
+
 template <class Type> class Block
 {
 public:
@@ -122,7 +124,7 @@ public:
        (optionally) the pointer to the function which
        will be called if allocation failed; the message
        passed to this function is "Not enough memory!" */
-    Block(int size, void (*err_function)(char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; }
+    Block(int size, void (*err_function)(const char *) = NULL) { first = last = NULL; block_size = size; error_function = err_function; }
 
     /* Destructor. Deallocates all items added so far */
     ~Block() { while (first) { block *next = first -> next; delete first; first = next; } }
@@ -209,7 +211,7 @@ private:
     block   *scan_current_block;
     Type    *scan_current_data;
 
-    void    (*error_function)(char *);
+    void    (*error_function)(const char *);
 };
 
 /***********************************************************************/
@@ -223,7 +225,7 @@ public:
        (optionally) the pointer to the function which
        will be called if allocation failed; the message
        passed to this function is "Not enough memory!" */
-    DBlock(int size, void (*err_function)(char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; }
+    DBlock(int size, void (*err_function)(const char *) = NULL) { first = NULL; first_free = NULL; block_size = size; error_function = err_function; }
 
     /* Destructor. Deallocates all items added so far */
     ~DBlock() { while (first) { block *next = first -> next; delete first; first = next; } }
@@ -277,9 +279,10 @@ private:
     block       *first;
     block_item  *first_free;
 
-    void    (*error_function)(char *);
+    void    (*error_function)(const char *);
 };
 
+} // end namespace mrfLib
 
 #endif
 
diff -rupN MRF2.1/energy.h MRF-v2.1.src-patched/energy.h
--- MRF2.1/energy.h	2007-09-16 04:16:54.000000000 +0200
+++ MRF-v2.1.src-patched/energy.h	2012-08-13 12:25:46.000000000 +0200
@@ -107,12 +107,14 @@
     ///////////////////////////////////////////////////
 */
 
-#ifndef __ENERGY_H__
-#define __ENERGY_H__
+#ifndef __MRF_ENERGY_H__
+#define __MRF_ENERGY_H__
 
 #include <assert.h>
 #include "graph.h"
 
+namespace mrfLib {
+
 class Energy : Graph
 {
 public:
@@ -132,7 +134,7 @@ public:
        function which will be called if an error occurs;
        an error message is passed to this function. If this
        argument is omitted, exit(1) will be called. */
-    Energy(void (*err_function)(char *) = NULL);
+    Energy(void (*err_function)(const char *) = NULL);
 
     /* Destructor */
     ~Energy();
@@ -196,7 +198,7 @@ private:
     /* internal variables and functions */
 
     TotalValue  Econst;
-    void        (*error_function)(char *);  /* this function is called if a error occurs,
+    void        (*error_function)(const char *);  /* this function is called if a error occurs,
                                             with a corresponding error message
                                             (or exit(1) is called if it's NULL) */
 };
@@ -219,7 +221,7 @@ private:
 /************************  Implementation ******************************/
 /***********************************************************************/
 
-inline Energy::Energy(void (*err_function)(char *)) : Graph(err_function)
+inline Energy::Energy(void (*err_function)(const char *)) : Graph(err_function)
 {
     Econst = 0;
     error_function = err_function;
@@ -368,4 +370,6 @@ inline Energy::TotalValue Energy::minimi
 
 inline int Energy::get_var(Var x) { return (int)what_segment(x); }
 
+} // end namespace mrfLib
+
 #endif
diff -rupN MRF2.1/example.cpp MRF-v2.1.src-patched/example.cpp
--- MRF2.1/example.cpp	2007-10-22 21:45:15.000000000 +0200
+++ MRF-v2.1.src-patched/example.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -1,6 +1,6 @@
 // example.cpp -- illustrates calling the MRF code
 
-static char *usage = "usage: %s [energyType] (a number between 0 and 3)\n";
+static const char *usage = "usage: %s [energyType] (a number between 0 and 3)\n";
 
 // uncomment "#define COUNT_TRUNCATIONS" in energy.h to enable counting of truncations
 
@@ -16,6 +16,8 @@ static char *usage = "usage: %s [energyT
 #include <time.h>
 #include <new>
 
+using namespace mrfLib;
+
 const int sizeX = 50;
 const int sizeY = 50;
 const int numLabels = 20;
diff -rupN MRF2.1/graph.cpp MRF-v2.1.src-patched/graph.cpp
--- MRF2.1/graph.cpp	2006-11-15 04:27:04.000000000 +0100
+++ MRF-v2.1.src-patched/graph.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -21,7 +21,9 @@
 #include <stdio.h>
 #include "graph.h"
 
-Graph::Graph(void (*err_function)(char *))
+using namespace mrfLib;
+
+Graph::Graph(void (*err_function)(const char *))
 {
     error_function = err_function;
     node_block_first = NULL;
diff -rupN MRF2.1/graph.h MRF-v2.1.src-patched/graph.h
--- MRF2.1/graph.h	2006-11-15 04:27:04.000000000 +0100
+++ MRF-v2.1.src-patched/graph.h	2012-08-22 13:39:57.000000000 +0200
@@ -41,8 +41,8 @@
     and memory usage see README.TXT.
 */
 
-#ifndef __GRAPH_H__
-#define __GRAPH_H__
+#ifndef __MRF_GRAPH_H__
+#define __MRF_GRAPH_H__
 
 #ifndef	USE_64_BIT_PTR_CAST
 #define PTR_CAST int
@@ -54,14 +54,16 @@
 #include "block.h"
 #include "mrf.h"
 
+namespace mrfLib {
+
 /*
     Nodes, arcs and pointers to nodes are
     added in blocks for memory and time efficiency.
     Below are numbers of items in blocks
 */
-#define NODE_BLOCK_SIZE 512
-#define ARC_BLOCK_SIZE 1024
-#define NODEPTR_BLOCK_SIZE 128
+const unsigned int NODE_BLOCK_SIZE = 512;
+const unsigned int ARC_BLOCK_SIZE = 1024;
+const unsigned int NODEPTR_BLOCK_SIZE = 128;
 
 class Graph
 {
@@ -85,7 +87,7 @@ public:
        function which will be called if an error occurs;
        an error message is passed to this function. If this
        argument is omitted, exit(1) will be called. */
-    Graph(void (*err_function)(char *) = NULL);
+    Graph(void (*err_function)(const char *) = NULL);
 
     /* Destructor */
     ~Graph();
@@ -233,7 +235,7 @@ private:
     arc_rev_block       *arc_rev_block_first;
     DBlock<nodeptr>     *nodeptr_block;
 
-    void    (*error_function)(char *);  /* this function is called if a error occurs,
+    void    (*error_function)(const char *);  /* this function is called if a error occurs,
                                            with a corresponding error message
                                            (or exit(1) is called if it's NULL) */
 
@@ -258,4 +260,6 @@ private:
     void process_sink_orphan(node *i);
 };
 
+} // end namespace mrfLib
+
 #endif
diff -rupN MRF2.1/maxflow.cpp MRF-v2.1.src-patched/maxflow.cpp
--- MRF2.1/maxflow.cpp	2006-11-15 04:27:04.000000000 +0100
+++ MRF-v2.1.src-patched/maxflow.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -21,13 +21,15 @@
 #include <stdio.h>
 #include "graph.h"
 
+using namespace mrfLib;
+
 /*
     special constants for node->parent
 */
 #define TERMINAL ( (arc_forward *) 1 )      /* to terminal */
 #define ORPHAN   ( (arc_forward *) 2 )      /* orphan */
 
-#define INFINITE_D 1000000000       /* infinite distance to the terminal */
+#define MRF_INFINITE_D 1000000000       /* infinite distance to the terminal */
 
 /***********************************************************************/
 
@@ -289,7 +291,7 @@ void Graph::process_source_orphan(node *
     arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
     arc_forward *a0_min = NULL, *a;
     nodeptr *np;
-    int d, d_min = INFINITE_D;
+    int d, d_min = MRF_INFINITE_D;
 
     /* trying to find a new parent */
     a0_for_first = i -> first_out;
@@ -331,13 +333,13 @@ void Graph::process_source_orphan(node *
                     j -> DIST = 1;
                     break;
                 }
-                if (a==ORPHAN) { d = INFINITE_D; break; }
+                if (a==ORPHAN) { d = MRF_INFINITE_D; break; }
                 if (IS_ODD(a))
                     j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                 else
                     j = NEIGHBOR_NODE(j, a -> shift);
             }
-            if (d<INFINITE_D) /* j originates from the source - done */
+            if (d<MRF_INFINITE_D) /* j originates from the source - done */
             {
                 if (d<d_min)
                 {
@@ -383,13 +385,13 @@ void Graph::process_source_orphan(node *
                         j -> DIST = 1;
                         break;
                     }
-                    if (a==ORPHAN) { d = INFINITE_D; break; }
+                    if (a==ORPHAN) { d = MRF_INFINITE_D; break; }
                     if (IS_ODD(a))
                         j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                     else
                         j = NEIGHBOR_NODE(j, a -> shift);
                 }
-                if (d<INFINITE_D) /* j originates from the source - done */
+                if (d<MRF_INFINITE_D) /* j originates from the source - done */
                 {
                     if (d<d_min)
                     {
@@ -472,7 +474,7 @@ void Graph::process_sink_orphan(node *i)
     arc_reverse *a0_rev, *a0_rev_first, *a0_rev_last;
     arc_forward *a0_min = NULL, *a;
     nodeptr *np;
-    int d, d_min = INFINITE_D;
+    int d, d_min = MRF_INFINITE_D;
 
     /* trying to find a new parent */
     a0_for_first = i -> first_out;
@@ -514,13 +516,13 @@ void Graph::process_sink_orphan(node *i)
                     j -> DIST = 1;
                     break;
                 }
-                if (a==ORPHAN) { d = INFINITE_D; break; }
+                if (a==ORPHAN) { d = MRF_INFINITE_D; break; }
                 if (IS_ODD(a))
                     j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                 else
                     j = NEIGHBOR_NODE(j, a -> shift);
             }
-            if (d<INFINITE_D) /* j originates from the sink - done */
+            if (d<MRF_INFINITE_D) /* j originates from the sink - done */
             {
                 if (d<d_min)
                 {
@@ -566,13 +568,13 @@ void Graph::process_sink_orphan(node *i)
                         j -> DIST = 1;
                         break;
                     }
-                    if (a==ORPHAN) { d = INFINITE_D; break; }
+                    if (a==ORPHAN) { d = MRF_INFINITE_D; break; }
                     if (IS_ODD(a))
                         j = NEIGHBOR_NODE_REV(j, MAKE_EVEN(a) -> shift);
                     else
                         j = NEIGHBOR_NODE(j, a -> shift);
                 }
-                if (d<INFINITE_D) /* j originates from the sink - done */
+                if (d<MRF_INFINITE_D) /* j originates from the sink - done */
                 {
                     if (d<d_min)
                     {
diff -rupN MRF2.1/mrf.cpp MRF-v2.1.src-patched/mrf.cpp
--- MRF2.1/mrf.cpp	2006-11-15 04:27:06.000000000 +0100
+++ MRF-v2.1.src-patched/mrf.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -5,6 +5,7 @@
 #include <assert.h>
 #include "mrf.h"
 
+using namespace mrfLib;
 
 void MRF::initialize()
 {
diff -rupN MRF2.1/mrf.h MRF-v2.1.src-patched/mrf.h
--- MRF2.1/mrf.h	2006-11-27 16:15:42.000000000 +0100
+++ MRF-v2.1.src-patched/mrf.h	2012-08-13 12:25:46.000000000 +0200
@@ -6,6 +6,8 @@
 #define __MRF_H__
 #include <stdio.h>
 
+namespace mrfLib {
+
 class EnergyFunction;
 
 class MRF
@@ -37,9 +39,9 @@ public:
     
 
     // *********** EVALUATING THE ENERGY
-    typedef int Label;
-    typedef int EnergyVal;        /* The total energy of a labeling */
-    typedef int CostVal;          /* costs of individual terms of the energy */
+    typedef MRFLABELVALUE Label;
+    typedef MRFENERGYVALUE EnergyVal;        /* The total energy of a labeling */
+    typedef MRFCOSTVALUE CostVal;          /* costs of individual terms of the energy */
  
     EnergyVal totalEnergy();      /* returns energy of current labeling */
     virtual EnergyVal dataEnergy() = 0;        /* returns the data part of the energy */
@@ -249,6 +251,7 @@ public:
     SmoothnessCost *m_smoothCost;
 };
 
+} // end namespace mrfLib
 
 #endif /*  __MRF_H__ */
 
diff -rupN MRF2.1/mrflib.h MRF-v2.1.src-patched/mrflib.h
--- MRF2.1/mrflib.h	1970-01-01 01:00:00.000000000 +0100
+++ MRF-v2.1.src-patched/mrflib.h	2012-08-13 12:25:46.000000000 +0200
@@ -0,0 +1,11 @@
+#ifndef MRFLIB_H_
+#define MRFLIB_H_
+
+#include "mrf.h"
+#include "ICM.h"
+#include "GCoptimization.h"
+#include "MaxProdBP.h"
+#include "TRW-S.h"
+#include "BP-S.h"
+
+#endif /* MRFLIB_H_ */
diff -rupN MRF2.1/regions-maxprod.cpp MRF-v2.1.src-patched/regions-maxprod.cpp
--- MRF2.1/regions-maxprod.cpp	2006-11-27 16:14:24.000000000 +0100
+++ MRF-v2.1.src-patched/regions-maxprod.cpp	2012-08-13 12:25:46.000000000 +0200
@@ -4,6 +4,9 @@
 #include <stdio.h>
 #include "MaxProdBP.h"
 #include "assert.h"
+
+using namespace mrfLib;
+
 int numIterRun;
 // Some of the GBP code has been disabled here
 
@@ -21,6 +24,8 @@ OneNodeCluster::OneNodeCluster()
 
 int OneNodeCluster::numStates;
 
+namespace mrfLib {
+
 FLOATTYPE vec_min(FLOATTYPE *vec, int length)
 {
 
@@ -292,6 +297,8 @@ inline void l2_dist_trans_comp(FLOATTYPE
 
 }
 
+} // end namespace mrfLib
+
 void OneNodeCluster::ComputeMsgRight(FLOATTYPE *msgDest, int r, int c, MaxProdBP *mrf)
 {
 
@@ -331,7 +338,7 @@ void OneNodeCluster::ComputeMsgRight(FLO
   {
     FLOATTYPE *psiMat, var_weight;
   
-    getPsiMat(*this,psiMat,r,c,mrf,RIGHT, var_weight);
+    mrfLib::getPsiMat(*this,psiMat,r,c,mrf,RIGHT, var_weight);
     FLOATTYPE *cmessage = msgDest;    
     for(int rightNodeInd = 0; rightNodeInd < numStates; rightNodeInd++)
     {
@@ -406,7 +413,7 @@ void OneNodeCluster::ComputeMsgLeft(FLOA
   {
     FLOATTYPE *psiMat, var_weight;
     
-    getPsiMat(*this,psiMat,r,c,mrf,LEFT, var_weight);
+    mrfLib::getPsiMat(*this,psiMat,r,c,mrf,LEFT, var_weight);
   
     FLOATTYPE *cmessage = msgDest;
     
@@ -632,6 +639,7 @@ int OneNodeCluster::getBeliefMaxInd()
 
 }
 
+namespace mrfLib {
 
 void computeMessagesLeftRight(OneNodeCluster *nodeArray, const int numCols, const int /*numRows*/, const int currRow, const FLOATTYPE alpha, MaxProdBP *mrf)
 {
@@ -691,4 +699,5 @@ void computeMessagesUpDown(OneNodeCluste
 
 }
 
+} // end namespace mrfLib
 
diff -rupN MRF2.1/regions-new.h MRF-v2.1.src-patched/regions-new.h
--- MRF2.1/regions-new.h	2006-11-15 04:27:06.000000000 +0100
+++ MRF-v2.1.src-patched/regions-new.h	2012-08-13 12:25:46.000000000 +0200
@@ -1,7 +1,7 @@
 // (C) 2002 Marshall Tappen, MIT AI Lab
 
-#ifndef _reg_h
-#define _reg_h
+#ifndef _MRF_reg_h
+#define _MRF_reg_h
 
 #define FLOATTYPE float
 #define UP 0
@@ -10,6 +10,8 @@
 #define RIGHT 3
 #include "MaxProdBP.h"
 
+namespace mrfLib {
+
 class MaxProdBP;
 class OneNodeCluster
 {
@@ -53,4 +55,6 @@ void computeOneNodeMessagesPeriodic(OneN
                                     const int numCols, const FLOATTYPE alpha);
 
 
+} // end namespace mrfLib
+
 #endif
diff -rupN MRF2.1/typeTruncatedQuadratic2D.h MRF-v2.1.src-patched/typeTruncatedQuadratic2D.h
--- MRF2.1/typeTruncatedQuadratic2D.h	2007-03-20 04:22:00.000000000 +0100
+++ MRF-v2.1.src-patched/typeTruncatedQuadratic2D.h	2006-11-18 15:45:12.000000000 +0100
@@ -1,32 +1,574 @@
-#ifndef __AJALSOQJAJSDFASD_H__
-#define __AJALSOQJAJSDFASD_H__
-
-#include <stdio.h>
-
-// This file is a stub, to make TRWS compile.
-
-
-struct TypeTruncatedQuadratic2D
-{
-    typedef double REAL;
-
-    struct Edge
-    {
-	void DistanceTransformL2(int /*K*/, int /*stride*/, REAL /*alpha*/, REAL* /*source*/, REAL* /*dest*/,
-				 int* /*parabolas*/, int* /*intersections*/)
-	{
-	    printf("\n\
-+-------------------------------------------------------------+\n\
-|   In order to run TRW-S with truncted L2 terms,             |\n\
-|   you need to download the implementation from              |\n\
-|      http://www.adastral.ucl.ac.uk/~vnk/papers/TRW-S.html   |\n\
-|   and copy file  typeTruncatedQuadratic2D.h                 |\n\
-|   to the main directory (i.e. replace the existing file)    |\n\
-+-------------------------------------------------------------+\n\
-			");
-	    exit(1);
-	}
-    };
-};
-
-#endif
+/******************************************************************
+typeTruncatedQuadratic2D.h
+
+Energy function with 2-dimensional truncated quadratic interactions:
+   E(x)   =   \sum_i D_i(x_i)   +   \sum_ij V_ij(x_i,x_j)
+   where x_i \in {0, 1, ..., KX-1} x {0, 1, ..., KY-1}
+   (i.e. x_i = (x_i(1), x_i(2))).
+   V_ij(ki, kj) = min { alpha1_ij*(ki(1)-kj(2))^2 + alpha2_ij*(ki(2)-kj(2))^2, lambda_ij }.
+   alpha1_ij, alpha2_ij and lambda_ij must be non-negative.
+
+Example usage:
+
+Minimize function E(x,y) = Dx(x) + Dy(y) + min { alpha1*(x(1) - y(1))^2 + alpha2*(x(2) - y(2))^2, lambda } 
+where 
+  x,y \in {0,1} x {0,1,2}
+  Dx(0,0) = 0, Dx(0,1) = 1, Dx(0,2) = 2,
+  Dx(1,0) = 3, Dx(1,1) = 4, Dx(1,2) = 5,
+  Dy(y) = 0 for all y,
+  alpha1 = 6,
+  alpha2 = 7,
+  lambda = 8
+
+
+
+
+#include <stdio.h>
+#include "MRFEnergy.h"
+
+void testTruncatedQuadratic2D()
+{
+	MRFEnergy<TypeTruncatedQuadratic2D>* mrf;
+	MRFEnergy<TypeTruncatedQuadratic2D>::NodeId* nodes;
+	MRFEnergy<TypeTruncatedQuadratic2D>::Options options;
+	TypeTruncatedQuadratic2D::REAL energy, lowerBound;
+
+	const int nodeNum = 2; // number of nodes
+	const int KX = 2; //   label
+	const int KY = 3; // dimensions
+	const int K = KX*KY; // number of labels
+	TypeTruncatedQuadratic2D::REAL D[K];
+	TypeTruncatedQuadratic2D::Label x, y;
+
+	mrf = new MRFEnergy<TypeTruncatedQuadratic2D>(TypeTruncatedQuadratic2D::GlobalSize(KX, KY));
+	nodes = new MRFEnergy<TypeTruncatedQuadratic2D>::NodeId[nodeNum];
+
+	// construct energy
+	D[0] = 0; D[1] = 1; D[2] = 2;
+	D[3] = 3; D[4] = 4; D[5] = 5;
+	nodes[0] = mrf->AddNode(TypeTruncatedQuadratic2D::LocalSize(), TypeTruncatedQuadratic2D::NodeData(D));
+	D[0] = 0; D[1] = 0; D[2] = 0;
+	D[3] = 0; D[4] = 0; D[5] = 0;
+	nodes[1] = mrf->AddNode(TypeTruncatedQuadratic2D::LocalSize(), TypeTruncatedQuadratic2D::NodeData(D));
+	mrf->AddEdge(nodes[0], nodes[1], TypeTruncatedQuadratic2D::EdgeData(6, 7, 8));
+
+	// Function below is optional - it may help if, for example, nodes are added in a random order
+	// mrf->SetAutomaticOrdering();
+
+	/////////////////////// TRW-S algorithm //////////////////////
+	options.m_iterMax = 30; // maximum number of iterations
+	mrf->Minimize_TRW_S(options, lowerBound, energy);
+
+	// read solution
+	x = mrf->GetSolution(nodes[0]);
+	y = mrf->GetSolution(nodes[1]);
+
+	printf("Solution: %d %d\n", x, y);
+
+	//////////////////////// BP algorithm ////////////////////////
+	mrf->ZeroMessages(); // in general not necessary - it may be faster to start 
+	                     // with messages computed in previous iterations
+
+	options.m_iterMax = 30; // maximum number of iterations
+	mrf->Minimize_BP(options, energy);
+
+	// read solution
+	x = mrf->GetSolution(nodes[0]);
+	y = mrf->GetSolution(nodes[1]);
+
+	printf("Solution: (%d, %d) (%d, %d)\n", x.m_kx, x.m_ky, y.m_kx, y.m_ky);
+
+	// done
+	delete nodes;
+	delete mrf;
+}
+
+*******************************************************************/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#ifndef __TYPETRUNCATEDQUADRATIC2D_H__
+#define __TYPETRUNCATEDQUADRATIC2D_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeTruncatedQuadratic2D
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	// types declarations
+	struct Label
+	{
+		int		m_kx, m_ky;
+	};
+	typedef double REAL;
+	struct GlobalSize; // global information about number of labels
+	struct LocalSize; // local information about number of labels (stored at each node)
+	struct NodeData; // argument to MRFEnergy::AddNode()
+	struct EdgeData; // argument to MRFEnergy::AddEdge()
+
+
+	struct GlobalSize
+	{
+		GlobalSize(int KX, int KY);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		int		m_KX, m_KY; // label dimensions
+		int		m_K; // number of labels
+	};
+
+	struct LocalSize // number of labels is stored at MRFEnergy::m_Kglobal
+	{
+	};
+
+	struct NodeData
+	{
+		NodeData(REAL* data); // data = pointer to array of size MRFEnergy::m_Kglobal
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL*		m_data;
+	};
+
+	struct EdgeData
+	{
+		EdgeData(REAL alphaX, REAL alphaY, REAL lambda);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_alphaX, m_alphaY;
+		REAL		m_lambda;
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypeTruncatedQuadratic2D>;
+
+	struct Vector
+	{
+		static int GetSizeInBytes(GlobalSize Kglobal, LocalSize K); // returns -1 if invalid K's
+		void Initialize(GlobalSize Kglobal, LocalSize K, NodeData data);  // called once when user adds a node
+		void Add(GlobalSize Kglobal, LocalSize K, NodeData data); // called once when user calls MRFEnergy::AddNodeData()
+
+		void SetZero(GlobalSize Kglobal, LocalSize K);                            // set this[k] = 0
+		void Copy(GlobalSize Kglobal, LocalSize K, Vector* V);                    // set this[k] = V[k]
+		void Add(GlobalSize Kglobal, LocalSize K, Vector* V);                     // set this[k] = this[k] + V[k]
+		REAL GetValue(GlobalSize Kglobal, LocalSize K, Label k);                  // return this[k]
+		REAL ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin);            // return vMin = min_k { this[k] }, set kMin
+		REAL ComputeAndSubtractMin(GlobalSize Kglobal, LocalSize K);              // same as previous, but additionally set this[k] -= vMin (and kMin is not returned)
+
+		static int GetArraySize(GlobalSize Kglobal, LocalSize K);
+		REAL GetArrayValue(GlobalSize Kglobal, LocalSize K, int k); // note: k is an integer in [0..GetArraySize()-1].
+		                                                            // For Potts functions GetArrayValue() and GetValue() are the same,
+		                                                            // but they are different for, say, 2-dimensional labels.
+		void SetArrayValue(GlobalSize Kglobal, LocalSize K, int k, REAL x);
+
+	private:
+	friend struct Edge;
+		REAL		m_data[1]; // actual size is MRFEnergy::m_Kglobal
+	};
+
+	struct Edge
+	{
+		static int GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data); // returns -1 if invalid data
+		static int GetBufSizeInBytes(int vectorMaxSizeInBytes); // returns size of buffer need for UpdateMessage()
+		void Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj); // called once when user adds an edge
+		Vector* GetMessagePtr();
+		void Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj); // if the client calls this function, then the meaning of 'dir'
+								                                               // in distance transform functions is swapped
+
+		// When UpdateMessage() is called, edge contains message from dest to source.
+		// The function must replace it with the message from source to dest.
+		// The update rule is given below assuming that source corresponds to tail (i) and dest corresponds
+		// to head (j) (which is the case if dir==0).
+		//
+		// 1. Compute Di[ki] = gamma*source[ki] - message[ki].  (Note: message = message from j to i).
+		// 2. Compute distance transform: set
+		//       message[kj] = min_{ki} (Di[ki] + V(ki,kj)). (Note: message = message from i to j).
+		// 3. Compute vMin = min_{kj} m_message[kj].
+		// 4. Set m_message[kj] -= vMin.
+		// 5. Return vMin.
+		//
+		// If dir==1 then source corresponds to j, sink corresponds to i. Then the update rule is
+		//
+		// 1. Compute Dj[kj] = gamma*source[kj] - message[kj].  (Note: message = message from i to j).
+		// 2. Compute distance transform: set
+		//       message[ki] = min_{kj} (Dj[kj] + V(ki,kj)). (Note: message = message from j to i).
+		// 3. Compute vMin = min_{ki} m_message[ki].
+		// 4. Set m_message[ki] -= vMin.
+		// 5. Return vMin.
+		//
+		// If Edge::Swap has been called odd number of times, then the meaning of dir is swapped.
+		//
+		// Vector 'source' must not be modified. Function may use 'buf' as a temporary storage.
+		REAL UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* buf);
+
+		// If dir==0, then sets dest[kj] += V(ksource,kj).
+		// If dir==1, then sets dest[ki] += V(ki,ksource).
+		// If Swap() has been called odd number of times, then the meaning of dir is swapped.
+		void AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir);
+
+	private:
+		// parabolas must be array of size K
+		// intersections must be array of size K+1
+		void DistanceTransformL2(int K, int stride, REAL alpha, REAL* source, REAL* dest, int* parabolas, int* intersections);
+
+		// edge information
+		REAL		m_alphaX, m_alphaY;
+		REAL		m_lambda;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+inline TypeTruncatedQuadratic2D::GlobalSize::GlobalSize(int KX, int KY)
+{
+	m_KX = KX;
+	m_KY = KY;
+	m_K = KX*KY;
+}
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeTruncatedQuadratic2D::NodeData::NodeData(REAL* data)
+{
+	m_data = data;
+}
+
+inline TypeTruncatedQuadratic2D::EdgeData::EdgeData(REAL alphaX, REAL alphaY, REAL lambda)
+{
+	m_alphaX = alphaX;
+	m_alphaY = alphaY;
+	m_lambda = lambda;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeTruncatedQuadratic2D::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	if (Kglobal.m_KX < 1 || Kglobal.m_KY < 1)
+	{
+		return -1;
+	}
+	return Kglobal.m_K*sizeof(REAL);
+}
+inline void TypeTruncatedQuadratic2D::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	memcpy(m_data, data.m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::Add(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	for (int k=0; k<Kglobal.m_K; k++)
+	{
+		m_data[k] += data.m_data[k];
+	}
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	memset(m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	memcpy(m_data, V->m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::Add(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	for (int k=0; k<Kglobal.m_K; k++)
+	{
+		m_data[k] += V->m_data[k];
+	}
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k.m_kx>=0 && k.m_kx<Kglobal.m_KX && k.m_ky>=0 && k.m_ky<Kglobal.m_KY);
+	return m_data[k.m_kx + k.m_ky*Kglobal.m_KX];
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& _kMin)
+{
+	REAL vMin = m_data[0];
+	int kMin = 0;
+	for (int k=1; k<Kglobal.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+			kMin = k;
+		}
+	}
+
+	_kMin.m_ky = kMin / Kglobal.m_KX;
+	_kMin.m_kx = kMin - _kMin.m_ky * Kglobal.m_KX;
+
+	return vMin;
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::ComputeAndSubtractMin(GlobalSize Kglobal, LocalSize K)
+{
+	REAL vMin = m_data[0];
+	for (int k=1; k<Kglobal.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+		}
+	}
+	for (int k=0; k<Kglobal.m_K; k++)
+	{
+		m_data[k] -= vMin;
+	}
+
+	return vMin;
+}
+
+inline int TypeTruncatedQuadratic2D::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return Kglobal.m_K;
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::SetArrayValue(GlobalSize Kglobal, LocalSize K, int k, REAL x)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	m_data[k] = x;
+}
+
+///////////////////// EdgeDataAndMessage implementation /////////////////////////
+
+inline int TypeTruncatedQuadratic2D::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	if (data.m_alphaX < 0 || data.m_alphaY < 0 || data.m_lambda < 0)
+	{
+		return -1;
+	}
+	return sizeof(Edge) - sizeof(Vector) + Kglobal.m_K*sizeof(REAL);
+}
+
+inline int TypeTruncatedQuadratic2D::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	int K = vectorMaxSizeInBytes / sizeof(REAL);
+	return K*sizeof(REAL) + (2*K+1)*sizeof(int);
+}
+
+inline void TypeTruncatedQuadratic2D::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	m_alphaX = data.m_alphaX;
+	m_alphaY = data.m_alphaY;
+	m_lambda = data.m_lambda;
+	memset(m_message.m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline TypeTruncatedQuadratic2D::Vector* TypeTruncatedQuadratic2D::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypeTruncatedQuadratic2D::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* _buf)
+{
+	REAL* buf = (REAL*) _buf;
+	int* parabolas = (int*) ((char*)_buf + Kglobal.m_K*sizeof(REAL));
+	int* intersections = (int*) ((char*)_buf + Kglobal.m_K*(sizeof(REAL) + sizeof(int)));
+	Label k;
+	REAL vMin;
+	REAL* sourcePtr;
+	REAL* destPtr;
+
+	m_message.m_data[0] = gamma*source->m_data[0] - m_message.m_data[0];
+	vMin = m_message.m_data[0];
+
+	sourcePtr = source->m_data;
+	destPtr = m_message.m_data;
+	for (sourcePtr++, destPtr++ ; sourcePtr<&source->m_data[Kglobal.m_K]; sourcePtr++, destPtr++)
+	{
+		destPtr[0] = gamma*sourcePtr[0] - destPtr[0];
+		if (vMin > destPtr[0])
+		{
+			vMin = destPtr[0];
+		}
+	}
+
+	// distance transform in 'X' direction
+	for (k.m_ky=0; k.m_ky<Kglobal.m_KY; k.m_ky++)
+	{
+		DistanceTransformL2(Kglobal.m_KX, 1, m_alphaX, &m_message.m_data[k.m_ky*Kglobal.m_KX], buf+k.m_ky*Kglobal.m_KX, parabolas, intersections);
+	}
+
+	// distance transform in 'Y' direction
+	for (k.m_kx=0; k.m_kx<Kglobal.m_KX; k.m_kx++)
+	{
+		DistanceTransformL2(Kglobal.m_KY, Kglobal.m_KX, m_alphaY, buf+k.m_kx, &m_message.m_data[k.m_kx], parabolas, intersections);
+	}
+
+	for (destPtr=m_message.m_data; destPtr<&m_message.m_data[Kglobal.m_K]; destPtr++)
+	{
+		destPtr[0] -= vMin;
+		if (destPtr[0] > m_lambda)
+		{
+			destPtr[0] = m_lambda;
+		}
+	}
+
+	return vMin;
+}
+
+inline void TypeTruncatedQuadratic2D::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	assert(ksource.m_kx>=0 && ksource.m_kx<Kglobal.m_KX && ksource.m_ky>=0 && ksource.m_ky<Kglobal.m_KY);
+
+	REAL* destPtr;
+	Label k;
+
+	destPtr = dest->m_data;
+	for (k.m_ky=0; k.m_ky<Kglobal.m_KY; k.m_ky++)
+	for (k.m_kx=0; k.m_kx<Kglobal.m_KX; k.m_kx++, destPtr++)
+	{
+		REAL cost = m_alphaX*(k.m_kx - ksource.m_kx)*(k.m_kx - ksource.m_kx)
+		          + m_alphaY*(k.m_ky - ksource.m_ky)*(k.m_ky - ksource.m_ky);
+		destPtr[0] += (cost < m_lambda) ? cost : m_lambda;
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+
+
+inline void TypeTruncatedQuadratic2D::Edge::DistanceTransformL2(int K, int stride, REAL alpha, REAL* source, REAL* dest, int* parabolas, int* intersections)
+{
+	assert(alpha >= 0);
+
+	if (alpha == 0)
+	{
+		REAL* ptr;
+		REAL vMin = source[0];
+
+		for (ptr=source+stride; ptr<source+K*stride; ptr+=stride)
+		{
+			if (vMin > ptr[0])
+			{
+				vMin = ptr[0];
+			}
+		}
+		for (ptr=dest; ptr<dest+K*stride; ptr+=stride)
+		{
+			ptr[0] = vMin;
+		}
+		return;
+	}
+
+	int i, j, k, p;
+	int r = 0; // = number of parabolas minus 1
+	// parabolas[p] will be base of parabola p (0<=p<=r)
+	// intersections[p] will be intersection between parabolas p-1 and p (1<=p<=r)
+	// intersections[0] will be always 0
+	parabolas[0] = 0;
+	intersections[0] = 0;
+
+	for (i=1; i<K; i++)
+	{
+		while ( 1 )
+		{
+			j = parabolas[r]; // base of previous rightmost visible parabola
+			// k is intersection of parabolas i and j
+			k = (int) ( 1 + ( (i + j)  +  (source[stride*j] - source[stride*i]) / (alpha * (j - i)) ) / 2);
+			if (k >= K)
+			{
+				// i is not visible
+				break;
+			}
+			if (k < 0)
+			{
+				k = 0;
+			}
+
+			if (k > intersections[r])
+			{
+				// intersection is rightmost, add it to end
+				r ++;
+				parabolas[r] = i;
+				intersections[r] = k;
+				break;
+			}
+			// j is not visible
+			if (r == 0)
+			{
+				parabolas[0] = i;
+				break;
+			}
+			r --;
+		}
+	}
+
+	intersections[r + 1] = K;
+
+	i = 0;
+	for (p=0; p<=r; p++)
+	{
+		j = parabolas[p];
+		// i values in [intersections[p], intersections[p+1]) are assigned to j
+		for ( ; i<intersections[p+1]; i++)
+		{
+			dest[stride*i] = source[stride*j] + alpha*(i-j)*(i-j);
+		}
+	}
+}
+
+#endif
diff -rupN MRF2.1/typeTruncatedQuadratic2D_old.h MRF-v2.1.src-patched/typeTruncatedQuadratic2D_old.h
--- MRF2.1/typeTruncatedQuadratic2D_old.h	1970-01-01 01:00:00.000000000 +0100
+++ MRF-v2.1.src-patched/typeTruncatedQuadratic2D_old.h	2012-08-13 12:25:46.000000000 +0200
@@ -0,0 +1,574 @@
+/******************************************************************
+typeTruncatedQuadratic2D.h
+
+Energy function with 2-dimensional truncated quadratic interactions:
+   E(x)   =   \sum_i D_i(x_i)   +   \sum_ij V_ij(x_i,x_j)
+   where x_i \in {0, 1, ..., KX-1} x {0, 1, ..., KY-1}
+   (i.e. x_i = (x_i(1), x_i(2))).
+   V_ij(ki, kj) = min { alpha1_ij*(ki(1)-kj(2))^2 + alpha2_ij*(ki(2)-kj(2))^2, lambda_ij }.
+   alpha1_ij, alpha2_ij and lambda_ij must be non-negative.
+
+Example usage:
+
+Minimize function E(x,y) = Dx(x) + Dy(y) + min { alpha1*(x(1) - y(1))^2 + alpha2*(x(2) - y(2))^2, lambda } 
+where 
+  x,y \in {0,1} x {0,1,2}
+  Dx(0,0) = 0, Dx(0,1) = 1, Dx(0,2) = 2,
+  Dx(1,0) = 3, Dx(1,1) = 4, Dx(1,2) = 5,
+  Dy(y) = 0 for all y,
+  alpha1 = 6,
+  alpha2 = 7,
+  lambda = 8
+
+
+
+
+#include <stdio.h>
+#include "MRFEnergy.h"
+
+void testTruncatedQuadratic2D()
+{
+	MRFEnergy<TypeTruncatedQuadratic2D>* mrf;
+	MRFEnergy<TypeTruncatedQuadratic2D>::NodeId* nodes;
+	MRFEnergy<TypeTruncatedQuadratic2D>::Options options;
+	TypeTruncatedQuadratic2D::REAL energy, lowerBound;
+
+	const int nodeNum = 2; // number of nodes
+	const int KX = 2; //   label
+	const int KY = 3; // dimensions
+	const int K = KX*KY; // number of labels
+	TypeTruncatedQuadratic2D::REAL D[K];
+	TypeTruncatedQuadratic2D::Label x, y;
+
+	mrf = new MRFEnergy<TypeTruncatedQuadratic2D>(TypeTruncatedQuadratic2D::GlobalSize(KX, KY));
+	nodes = new MRFEnergy<TypeTruncatedQuadratic2D>::NodeId[nodeNum];
+
+	// construct energy
+	D[0] = 0; D[1] = 1; D[2] = 2;
+	D[3] = 3; D[4] = 4; D[5] = 5;
+	nodes[0] = mrf->AddNode(TypeTruncatedQuadratic2D::LocalSize(), TypeTruncatedQuadratic2D::NodeData(D));
+	D[0] = 0; D[1] = 0; D[2] = 0;
+	D[3] = 0; D[4] = 0; D[5] = 0;
+	nodes[1] = mrf->AddNode(TypeTruncatedQuadratic2D::LocalSize(), TypeTruncatedQuadratic2D::NodeData(D));
+	mrf->AddEdge(nodes[0], nodes[1], TypeTruncatedQuadratic2D::EdgeData(6, 7, 8));
+
+	// Function below is optional - it may help if, for example, nodes are added in a random order
+	// mrf->SetAutomaticOrdering();
+
+	/////////////////////// TRW-S algorithm //////////////////////
+	options.m_iterMax = 30; // maximum number of iterations
+	mrf->Minimize_TRW_S(options, lowerBound, energy);
+
+	// read solution
+	x = mrf->GetSolution(nodes[0]);
+	y = mrf->GetSolution(nodes[1]);
+
+	printf("Solution: %d %d\n", x, y);
+
+	//////////////////////// BP algorithm ////////////////////////
+	mrf->ZeroMessages(); // in general not necessary - it may be faster to start 
+	                     // with messages computed in previous iterations
+
+	options.m_iterMax = 30; // maximum number of iterations
+	mrf->Minimize_BP(options, energy);
+
+	// read solution
+	x = mrf->GetSolution(nodes[0]);
+	y = mrf->GetSolution(nodes[1]);
+
+	printf("Solution: (%d, %d) (%d, %d)\n", x.m_kx, x.m_ky, y.m_kx, y.m_ky);
+
+	// done
+	delete nodes;
+	delete mrf;
+}
+
+*******************************************************************/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#ifndef __TYPETRUNCATEDQUADRATIC2D_H__
+#define __TYPETRUNCATEDQUADRATIC2D_H__
+
+#include <string.h>
+#include <assert.h>
+
+
+template <class T> class MRFEnergy;
+
+
+class TypeTruncatedQuadratic2D
+{
+private:
+	struct Vector; // node parameters and messages
+	struct Edge; // stores edge information and either forward or backward message
+
+public:
+	// types declarations
+	struct Label
+	{
+		int		m_kx, m_ky;
+	};
+	typedef double REAL;
+	struct GlobalSize; // global information about number of labels
+	struct LocalSize; // local information about number of labels (stored at each node)
+	struct NodeData; // argument to MRFEnergy::AddNode()
+	struct EdgeData; // argument to MRFEnergy::AddEdge()
+
+
+	struct GlobalSize
+	{
+		GlobalSize(int KX, int KY);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		int		m_KX, m_KY; // label dimensions
+		int		m_K; // number of labels
+	};
+
+	struct LocalSize // number of labels is stored at MRFEnergy::m_Kglobal
+	{
+	};
+
+	struct NodeData
+	{
+		NodeData(REAL* data); // data = pointer to array of size MRFEnergy::m_Kglobal
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL*		m_data;
+	};
+
+	struct EdgeData
+	{
+		EdgeData(REAL alphaX, REAL alphaY, REAL lambda);
+
+	private:
+	friend struct Vector;
+	friend struct Edge;
+		REAL		m_alphaX, m_alphaY;
+		REAL		m_lambda;
+	};
+
+
+
+
+
+
+
+	//////////////////////////////////////////////////////////////////////////////////
+	////////////////////////// Visible only to MRFEnergy /////////////////////////////
+	//////////////////////////////////////////////////////////////////////////////////
+
+private:
+friend class MRFEnergy<TypeTruncatedQuadratic2D>;
+
+	struct Vector
+	{
+		static int GetSizeInBytes(GlobalSize Kglobal, LocalSize K); // returns -1 if invalid K's
+		void Initialize(GlobalSize Kglobal, LocalSize K, NodeData data);  // called once when user adds a node
+		void Add(GlobalSize Kglobal, LocalSize K, NodeData data); // called once when user calls MRFEnergy::AddNodeData()
+
+		void SetZero(GlobalSize Kglobal, LocalSize K);                            // set this[k] = 0
+		void Copy(GlobalSize Kglobal, LocalSize K, Vector* V);                    // set this[k] = V[k]
+		void Add(GlobalSize Kglobal, LocalSize K, Vector* V);                     // set this[k] = this[k] + V[k]
+		REAL GetValue(GlobalSize Kglobal, LocalSize K, Label k);                  // return this[k]
+		REAL ComputeMin(GlobalSize Kglobal, LocalSize K, Label& kMin);            // return vMin = min_k { this[k] }, set kMin
+		REAL ComputeAndSubtractMin(GlobalSize Kglobal, LocalSize K);              // same as previous, but additionally set this[k] -= vMin (and kMin is not returned)
+
+		static int GetArraySize(GlobalSize Kglobal, LocalSize K);
+		REAL GetArrayValue(GlobalSize Kglobal, LocalSize K, int k); // note: k is an integer in [0..GetArraySize()-1].
+		                                                            // For Potts functions GetArrayValue() and GetValue() are the same,
+		                                                            // but they are different for, say, 2-dimensional labels.
+		void SetArrayValue(GlobalSize Kglobal, LocalSize K, int k, REAL x);
+
+	private:
+	friend struct Edge;
+		REAL		m_data[1]; // actual size is MRFEnergy::m_Kglobal
+	};
+
+	struct Edge
+	{
+		static int GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data); // returns -1 if invalid data
+		static int GetBufSizeInBytes(int vectorMaxSizeInBytes); // returns size of buffer need for UpdateMessage()
+		void Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj); // called once when user adds an edge
+		Vector* GetMessagePtr();
+		void Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj); // if the client calls this function, then the meaning of 'dir'
+								                                               // in distance transform functions is swapped
+
+		// When UpdateMessage() is called, edge contains message from dest to source.
+		// The function must replace it with the message from source to dest.
+		// The update rule is given below assuming that source corresponds to tail (i) and dest corresponds
+		// to head (j) (which is the case if dir==0).
+		//
+		// 1. Compute Di[ki] = gamma*source[ki] - message[ki].  (Note: message = message from j to i).
+		// 2. Compute distance transform: set
+		//       message[kj] = min_{ki} (Di[ki] + V(ki,kj)). (Note: message = message from i to j).
+		// 3. Compute vMin = min_{kj} m_message[kj].
+		// 4. Set m_message[kj] -= vMin.
+		// 5. Return vMin.
+		//
+		// If dir==1 then source corresponds to j, sink corresponds to i. Then the update rule is
+		//
+		// 1. Compute Dj[kj] = gamma*source[kj] - message[kj].  (Note: message = message from i to j).
+		// 2. Compute distance transform: set
+		//       message[ki] = min_{kj} (Dj[kj] + V(ki,kj)). (Note: message = message from j to i).
+		// 3. Compute vMin = min_{ki} m_message[ki].
+		// 4. Set m_message[ki] -= vMin.
+		// 5. Return vMin.
+		//
+		// If Edge::Swap has been called odd number of times, then the meaning of dir is swapped.
+		//
+		// Vector 'source' must not be modified. Function may use 'buf' as a temporary storage.
+		REAL UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* buf);
+
+		// If dir==0, then sets dest[kj] += V(ksource,kj).
+		// If dir==1, then sets dest[ki] += V(ki,ksource).
+		// If Swap() has been called odd number of times, then the meaning of dir is swapped.
+		void AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir);
+
+	private:
+		// parabolas must be array of size K
+		// intersections must be array of size K+1
+		void DistanceTransformL2(int K, int stride, REAL alpha, REAL* source, REAL* dest, int* parabolas, int* intersections);
+
+		// edge information
+		REAL		m_alphaX, m_alphaY;
+		REAL		m_lambda;
+
+		// message
+		Vector		m_message;
+	};
+};
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////////
+/////////////////////////////// Implementation ///////////////////////////////////
+//////////////////////////////////////////////////////////////////////////////////
+
+
+inline TypeTruncatedQuadratic2D::GlobalSize::GlobalSize(int KX, int KY)
+{
+	m_KX = KX;
+	m_KY = KY;
+	m_K = KX*KY;
+}
+
+///////////////////// NodeData and EdgeData ///////////////////////
+
+inline TypeTruncatedQuadratic2D::NodeData::NodeData(REAL* data)
+{
+	m_data = data;
+}
+
+inline TypeTruncatedQuadratic2D::EdgeData::EdgeData(REAL alphaX, REAL alphaY, REAL lambda)
+{
+	m_alphaX = alphaX;
+	m_alphaY = alphaY;
+	m_lambda = lambda;
+}
+
+///////////////////// Vector ///////////////////////
+
+inline int TypeTruncatedQuadratic2D::Vector::GetSizeInBytes(GlobalSize Kglobal, LocalSize K)
+{
+	if (Kglobal.m_KX < 1 || Kglobal.m_KY < 1)
+	{
+		return -1;
+	}
+	return Kglobal.m_K*sizeof(REAL);
+}
+inline void TypeTruncatedQuadratic2D::Vector::Initialize(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	memcpy(m_data, data.m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::Add(GlobalSize Kglobal, LocalSize K, NodeData data)
+{
+	for (int k=0; k<Kglobal.m_K; k++)
+	{
+		m_data[k] += data.m_data[k];
+	}
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::SetZero(GlobalSize Kglobal, LocalSize K)
+{
+	memset(m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::Copy(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	memcpy(m_data, V->m_data, Kglobal.m_K*sizeof(REAL));
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::Add(GlobalSize Kglobal, LocalSize K, Vector* V)
+{
+	for (int k=0; k<Kglobal.m_K; k++)
+	{
+		m_data[k] += V->m_data[k];
+	}
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::GetValue(GlobalSize Kglobal, LocalSize K, Label k)
+{
+	assert(k.m_kx>=0 && k.m_kx<Kglobal.m_KX && k.m_ky>=0 && k.m_ky<Kglobal.m_KY);
+	return m_data[k.m_kx + k.m_ky*Kglobal.m_KX];
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::ComputeMin(GlobalSize Kglobal, LocalSize K, Label& _kMin)
+{
+	REAL vMin = m_data[0];
+	int kMin = 0;
+	for (int k=1; k<Kglobal.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+			kMin = k;
+		}
+	}
+
+	_kMin.m_ky = kMin / Kglobal.m_KX;
+	_kMin.m_kx = kMin - _kMin.m_ky * Kglobal.m_KX;
+
+	return vMin;
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::ComputeAndSubtractMin(GlobalSize Kglobal, LocalSize K)
+{
+	REAL vMin = m_data[0];
+	for (int k=1; k<Kglobal.m_K; k++)
+	{
+		if (vMin > m_data[k])
+		{
+			vMin = m_data[k];
+		}
+	}
+	for (int k=0; k<Kglobal.m_K; k++)
+	{
+		m_data[k] -= vMin;
+	}
+
+	return vMin;
+}
+
+inline int TypeTruncatedQuadratic2D::Vector::GetArraySize(GlobalSize Kglobal, LocalSize K)
+{
+	return Kglobal.m_K;
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Vector::GetArrayValue(GlobalSize Kglobal, LocalSize K, int k)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	return m_data[k];
+}
+
+inline void TypeTruncatedQuadratic2D::Vector::SetArrayValue(GlobalSize Kglobal, LocalSize K, int k, REAL x)
+{
+	assert(k>=0 && k<Kglobal.m_K);
+	m_data[k] = x;
+}
+
+///////////////////// EdgeDataAndMessage implementation /////////////////////////
+
+inline int TypeTruncatedQuadratic2D::Edge::GetSizeInBytes(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data)
+{
+	if (data.m_alphaX < 0 || data.m_alphaY < 0 || data.m_lambda < 0)
+	{
+		return -1;
+	}
+	return sizeof(Edge) - sizeof(Vector) + Kglobal.m_K*sizeof(REAL);
+}
+
+inline int TypeTruncatedQuadratic2D::Edge::GetBufSizeInBytes(int vectorMaxSizeInBytes)
+{
+	int K = vectorMaxSizeInBytes / sizeof(REAL);
+	return K*sizeof(REAL) + (2*K+1)*sizeof(int);
+}
+
+inline void TypeTruncatedQuadratic2D::Edge::Initialize(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj, EdgeData data, Vector* Di, Vector* Dj)
+{
+	m_alphaX = data.m_alphaX;
+	m_alphaY = data.m_alphaY;
+	m_lambda = data.m_lambda;
+	memset(m_message.m_data, 0, Kglobal.m_K*sizeof(REAL));
+}
+
+inline TypeTruncatedQuadratic2D::Vector* TypeTruncatedQuadratic2D::Edge::GetMessagePtr()
+{
+	return &m_message;
+}
+
+inline void TypeTruncatedQuadratic2D::Edge::Swap(GlobalSize Kglobal, LocalSize Ki, LocalSize Kj)
+{
+}
+
+inline TypeTruncatedQuadratic2D::REAL TypeTruncatedQuadratic2D::Edge::UpdateMessage(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Vector* source, REAL gamma, int dir, void* _buf)
+{
+	REAL* buf = (REAL*) _buf;
+	int* parabolas = (int*) ((char*)_buf + Kglobal.m_K*sizeof(REAL));
+	int* intersections = (int*) ((char*)_buf + Kglobal.m_K*(sizeof(REAL) + sizeof(int)));
+	Label k;
+	REAL vMin;
+	REAL* sourcePtr;
+	REAL* destPtr;
+
+	m_message.m_data[0] = gamma*source->m_data[0] - m_message.m_data[0];
+	vMin = m_message.m_data[0];
+
+	sourcePtr = source->m_data;
+	destPtr = m_message.m_data;
+	for (sourcePtr++, destPtr++ ; sourcePtr<&source->m_data[Kglobal.m_K]; sourcePtr++, destPtr++)
+	{
+		destPtr[0] = gamma*sourcePtr[0] - destPtr[0];
+		if (vMin > destPtr[0])
+		{
+			vMin = destPtr[0];
+		}
+	}
+
+	// distance transform in 'X' direction
+	for (k.m_ky=0; k.m_ky<Kglobal.m_KY; k.m_ky++)
+	{
+		DistanceTransformL2(Kglobal.m_KX, 1, m_alphaX, &m_message.m_data[k.m_ky*Kglobal.m_KX], buf+k.m_ky*Kglobal.m_KX, parabolas, intersections);
+	}
+
+	// distance transform in 'Y' direction
+	for (k.m_kx=0; k.m_kx<Kglobal.m_KX; k.m_kx++)
+	{
+		DistanceTransformL2(Kglobal.m_KY, Kglobal.m_KX, m_alphaY, buf+k.m_kx, &m_message.m_data[k.m_kx], parabolas, intersections);
+	}
+
+	for (destPtr=m_message.m_data; destPtr<&m_message.m_data[Kglobal.m_K]; destPtr++)
+	{
+		destPtr[0] -= vMin;
+		if (destPtr[0] > m_lambda)
+		{
+			destPtr[0] = m_lambda;
+		}
+	}
+
+	return vMin;
+}
+
+inline void TypeTruncatedQuadratic2D::Edge::AddColumn(GlobalSize Kglobal, LocalSize Ksource, LocalSize Kdest, Label ksource, Vector* dest, int dir)
+{
+	assert(ksource.m_kx>=0 && ksource.m_kx<Kglobal.m_KX && ksource.m_ky>=0 && ksource.m_ky<Kglobal.m_KY);
+
+	REAL* destPtr;
+	Label k;
+
+	destPtr = dest->m_data;
+	for (k.m_ky=0; k.m_ky<Kglobal.m_KY; k.m_ky++)
+	for (k.m_kx=0; k.m_kx<Kglobal.m_KX; k.m_kx++, destPtr++)
+	{
+		REAL cost = m_alphaX*(k.m_kx - ksource.m_kx)*(k.m_kx - ksource.m_kx)
+		          + m_alphaY*(k.m_ky - ksource.m_ky)*(k.m_ky - ksource.m_ky);
+		destPtr[0] += (cost < m_lambda) ? cost : m_lambda;
+	}
+}
+
+//////////////////////////////////////////////////////////////////////////////////
+
+
+
+inline void TypeTruncatedQuadratic2D::Edge::DistanceTransformL2(int K, int stride, REAL alpha, REAL* source, REAL* dest, int* parabolas, int* intersections)
+{
+	assert(alpha >= 0);
+
+	if (alpha == 0)
+	{
+		REAL* ptr;
+		REAL vMin = source[0];
+
+		for (ptr=source+stride; ptr<source+K*stride; ptr+=stride)
+		{
+			if (vMin > ptr[0])
+			{
+				vMin = ptr[0];
+			}
+		}
+		for (ptr=dest; ptr<dest+K*stride; ptr+=stride)
+		{
+			ptr[0] = vMin;
+		}
+		return;
+	}
+
+	int i, j, k, p;
+	int r = 0; // = number of parabolas minus 1
+	// parabolas[p] will be base of parabola p (0<=p<=r)
+	// intersections[p] will be intersection between parabolas p-1 and p (1<=p<=r)
+	// intersections[0] will be always 0
+	parabolas[0] = 0;
+	intersections[0] = 0;
+
+	for (i=1; i<K; i++)
+	{
+		while ( 1 )
+		{
+			j = parabolas[r]; // base of previous rightmost visible parabola
+			// k is intersection of parabolas i and j
+			k = (int) ( 1 + ( (i + j)  +  (source[stride*j] - source[stride*i]) / (alpha * (j - i)) ) / 2);
+			if (k >= K)
+			{
+				// i is not visible
+				break;
+			}
+			if (k < 0)
+			{
+				k = 0;
+			}
+
+			if (k > intersections[r])
+			{
+				// intersection is rightmost, add it to end
+				r ++;
+				parabolas[r] = i;
+				intersections[r] = k;
+				break;
+			}
+			// j is not visible
+			if (r == 0)
+			{
+				parabolas[0] = i;
+				break;
+			}
+			r --;
+		}
+	}
+
+	intersections[r + 1] = K;
+
+	i = 0;
+	for (p=0; p<=r; p++)
+	{
+		j = parabolas[p];
+		// i values in [intersections[p], intersections[p+1]) are assigned to j
+		for ( ; i<intersections[p+1]; i++)
+		{
+			dest[stride*i] = source[stride*j] + alpha*(i-j)*(i-j);
+		}
+	}
+}
+
+#endif
